Data Science II

Session 3: Choropleths, Classification, and Interactive Maps

Bogdan G. Popescu
bogdan.popescu@johncabot.edu

John Cabot University

What You’ll Learn Today

  1. Building a choropleth of Mexico from a table and a set of boundaries
  2. Choosing colour — and why viridis rather than a rainbow
  3. Classification: the same data, four schemes, four different stories
  4. The count-versus-rate trap, and the small-numbers trap behind it
  5. Interactive maps with leaflet, embedded in your Quarto report

Our variable today: intentional homicides per 100,000 inhabitants, by state and by municipio.

Before We Start

Code reads ../data/..., exactly as in Sessions 1 and 2:

#Step1: Libraries for today
library(sf)          # spatial objects
library(dplyr)       # data wrangling
library(ggplot2)     # static maps
library(classInt)    # classification schemes
library(viridis)     # colour palettes
library(leaflet)     # interactive maps
library(patchwork)   # side-by-side plots

sf_use_s2(FALSE)

The complete, ordered code for this session is in session3_code.R.

Today’s Data

../data/mexico_homicides.gpkg — a GeoPackage with two layers:

Layer Rows Contents
states 32 boundaries, population, homicides, rate
municipios 2,469 the same, at municipio level
  • Homicides — SESNSP homicidio doloso, annual mean over 2021–2023
  • Population — INEGI Censo 2020, via the mxmaps package
  • Boundaries — INEGI geometry with CVEGEO codes

A three-year mean, not a single year, because one bad year in a small municipio is noise, not signal.

Reading It In

#Step1: A GeoPackage can hold several layers -- ask what is inside
st_layers("../data/mexico_homicides.gpkg")$name
[1] "municipios" "states"    
#Step2: Read the state layer
sta <- st_read("../data/mexico_homicides.gpkg", layer = "states", quiet = TRUE)

#Step3: Look at it
head(st_drop_geometry(sta), 3)
  cve_ent          state_name     pop  homicides      rate
1      01      Aguascalientes 1425607   76.66667  5.377826
2      02     Baja California 3769020 2389.66667 63.402865
3      03 Baja California Sur  798447   37.00000  4.633996

Sanity-Check Before You Map

#Step1: Do the national totals look right?
nrow(sta)
[1] 32
format(sum(sta$pop), big.mark = ",")
[1] "126,014,024"
round(sum(sta$homicides))
[1] 26612

126 million people, about 26,600 homicides a year. Both match the published national figures.

Important

Always do this. A choropleth will happily draw nonsense. Checking that your totals reproduce a number you can look up elsewhere is the cheapest error-catching you will ever do.

Part A: From a Table to a Choropleth

Choropleth Maps

A choropleth is a thematic map where regions are coloured according to a statistical variable.

Three decisions, and every one of them changes the story:

  1. What variable — a count, or a rate?
  2. What colours — which palette, in which direction?
  3. What breaks — how do continuous values become discrete classes?

Most bad maps are bad because of decision 1 or decision 3.

The Join Key

Homicides come from SESNSP. Population comes from the census. Boundaries come from INEGI. Three sources, three files.

They are joined on CVEGEO — INEGI’s geostatistical code. 31050 is Mérida, Yucatán: 31 for the state, 050 for the municipio.

Warning

Never join Mexican data on place names. GADM calls the capital Distrito Federal; INEGI calls it Ciudad de México. There are two municipios called Guadalupe, and several spellings of Tláhuac. Codes do not have this problem.

Your First Choropleth

ggplot(sta) +
  geom_sf(aes(fill = rate)) +
  coord_sf(crs = 6372) +
  theme_void()

It works. The north-west and the Pacific centre are visibly worse off than the south-east.

Now we make it readable.

Part B: Colour

Three Kinds of Palette

Type Use when Example
Sequential low → high homicide rate
Diverging departure from a midpoint change since 2015
Qualitative unordered categories dominant party, land use

Using a qualitative palette for ordered data — or a rainbow for anything — invents structure that is not in the data.

The viridis Family

ggplot(sta) +
  geom_sf(aes(fill = rate)) +
  scale_fill_viridis_c(option = "magma") +
  coord_sf(crs = 6372) +
  theme_void()

The viridis Family

ggplot(sta) +
  geom_sf(aes(fill = rate)) +
  scale_fill_viridis_c(option = "inferno") +
  coord_sf(crs = 6372) +
  theme_void()

The viridis Family

ggplot(sta) +
  geom_sf(aes(fill = rate)) +
  scale_fill_viridis_c(option = "plasma") +
  coord_sf(crs = 6372) +
  theme_void()

The viridis Family

ggplot(sta) +
  geom_sf(aes(fill = rate)) +
  scale_fill_viridis_c(option = "viridis") +
  coord_sf(crs = 6372) +
  theme_void()

The viridis Family

ggplot(sta) +
  geom_sf(aes(fill = rate)) +
  scale_fill_viridis_c(option = "cividis") +
  coord_sf(crs = 6372) +
  theme_void()

The viridis Family

ggplot(sta) +
  geom_sf(aes(fill = rate)) +
  scale_fill_viridis_c(option = "turbo") +
  coord_sf(crs = 6372) +
  theme_void()

Why viridis?

  • Perceptually uniform — equal steps in data look like equal steps in colour
  • Colourblind-safe — about 8% of men have some form of colour vision deficiency
  • Prints in greyscale — the order survives a black-and-white printer

Warning

turbo is the exception — it is a rainbow, and it fails all three tests. Look back at the turbo map: the eye jumps to the green states, which are in the middle of the distribution. Do not use it for a policy map.

Direction Matters

ggplot(sta) +
  geom_sf(aes(fill = rate)) +
  scale_fill_viridis_c(direction = -1) +
  coord_sf(crs = 6372) +
  theme_void()

direction = -1 makes dark mean high. For a variable like homicide, dark-is-bad matches what readers expect — and it survives photocopying better.

_c, _b, or _d?

Function Gives you Use for
scale_fill_viridis_c() continuous gradient smooth variables
scale_fill_viridis_b() binned steps choropleths
scale_fill_viridis_d() discrete categories factors

For choropleths you almost always want binned — the eye cannot read a continuous ramp accurately. Which raises the real question: where do the bins go?

Part C: Classification

Continuous to Discrete

A binned choropleth reduces a continuous variable to a handful of colours.

That reduction is a choice, and it is the single biggest lever on what the reader sees.

Never pick breaks before you have looked at the distribution.

The Distribution

ggplot(sta, aes(rate)) +
  geom_histogram(bins = 20,
                 fill = "#5a6b7d") +
  labs(x = "Homicides per 100,000",
       y = NULL) +
  theme_minimal()

Right-skewed, but not degenerate. Mean 23.7, standard deviation 21.7. Eleven states are below 10; four are above 50.

classInt: Four Ways to Cut

#Step1: Same data, same number of classes, very different cut points
round(classIntervals(sta$rate, n = 5, style = "equal")$brks, 1)
[1]  1.7 19.6 37.5 55.5 73.4 91.3
round(classIntervals(sta$rate, n = 5, style = "quantile")$brks, 1)
[1]  1.7  7.9 12.3 19.2 42.7 91.3
round(classIntervals(sta$rate, n = 5, style = "jenks")$brks, 1)
[1]  1.7  9.6 21.1 49.0 63.4 91.3

Now let us see what each one does to the map.

Equal Interval

#Step1: Cut the RANGE into five equal-width slices
brks <- classIntervals(sta$rate, n = 5, style = "equal")$brks

#Step2: Use those breaks as the colour steps
ggplot(sta) +
  geom_sf(aes(fill = rate)) +
  scale_fill_viridis_b(breaks = brks, direction = -1) +
  coord_sf(crs = 6372) + theme_void()

Equal Interval: What It Does

brks <- classIntervals(sta$rate, n = 5, style = "equal")$brks
table(cut(sta$rate, breaks = brks, include.lowest = TRUE))

[1.65,19.6] (19.6,37.5] (37.5,55.5] (55.5,73.4] (73.4,91.3] 
         21           3           4           3           1

Twenty-one of thirty-two states in the first class. Equal interval divides the range, and the range is stretched by Colima at 91.

It is honest about magnitude — the classes really are equally wide — but it wastes four of your five colours on eight states.

Use it when data are roughly uniform, or when round-number breaks matter more than resolution.

Quantile

#Step1: Cut so each class holds the same NUMBER of states
brks <- classIntervals(sta$rate, n = 5, style = "quantile")$brks
brks <- classIntervals(sta$rate, n = 5, style = "quantile")$brks
table(cut(sta$rate, breaks = brks, include.lowest = TRUE))

[1.65,7.87] (7.87,12.3] (12.3,19.2] (19.2,42.7] (42.7,91.3] 
          7           6           6           6           7

Quantile: The Trade-Off

Six or seven states per class, so the whole map varies and you can read the middle of the distribution.

But look at the colour bar: the swatches are wildly unequal in width. The top class runs from 43 to 91; the bottom from 1.7 to 8.

Warning

Quantile maps make rank visible and magnitude invisible. A state at 44 and a state at 91 get the same colour, though one is twice as violent.

Standard Deviation

#Step1: Breaks measured in standard deviations from the mean
round(classIntervals(sta$rate, n = 5, style = "sd")$brks, 1)
[1] -19.8   2.0  23.7  45.4  67.2  88.9 110.6

The first break is negative — but a homicide rate cannot be below zero.

Standard-deviation breaks assume a roughly symmetric distribution. Ours is right-skewed, so the scheme spends a class on values that cannot exist.

Use it for genuinely bell-shaped data, or for a diverging map around a meaningful centre — change since last year, or deviation from the national average.

Jenks Natural Breaks

#Step1: Let the algorithm find the gaps in the distribution
brks <- classIntervals(sta$rate, n = 5, style = "jenks")$brks

Jenks minimises variance within classes and maximises it between them — it puts the cuts where the data already has gaps. Here: 11, 11, 6, 3, 1.

Which Should You Use?

Scheme Good when Fails when
Equal interval data are uniform; breaks must be round numbers outliers stretch the range
Quantile you care about rank; want a map that varies magnitude matters
Std. deviation data are symmetric; diverging maps data are skewed
Jenks you want the data’s own structure maps must be comparable

Important

If you make several maps that must be compared — 2015 next to 2023, say — fix the breaks by hand and reuse them. Jenks chooses different cuts for every map, and the maps stop being comparable.

Part D: Counts, Rates, and Small Numbers

The Most Common Mistake

#Step1: The five states with the most homicides
sta %>% st_drop_geometry() %>% arrange(desc(homicides)) %>%
  transmute(state_name, homicides = round(homicides), rate = round(rate, 1)) %>% head(5)
       state_name homicides rate
1      Guanajuato      2679 43.4
2 Baja California      2390 63.4
3          México      2286 13.5
4       Michoacán      1880 39.6
5       Chihuahua      1835 49.0
#Step2: The five states with the highest rate
sta %>% st_drop_geometry() %>% arrange(desc(rate)) %>%
  transmute(state_name, homicides = round(homicides), rate = round(rate, 1)) %>% head(5)
       state_name homicides rate
1          Colima       668 91.3
2 Baja California      2390 63.4
3       Zacatecas       941 58.0
4         Morelos      1126 57.1
5       Chihuahua      1835 49.0

Two Different Countries

By count, the State of México is third — 2,286 homicides a year.

By rate, it is 13.5 per 100,000, well below the national average. It has 17 million people.

Colima has 668 homicides — barely a quarter of México’s — and the worst rate in the country at 91.3.

Important

A count map of almost any human event is, mostly, a map of where the people are. Almost always, you want a rate.

Count and Rate Compared

#Step1: Jenks breaks for each variable separately
b1 <- classIntervals(sta$homicides, n = 5, style = "jenks")$brks
b2 <- classIntervals(sta$rate,      n = 5, style = "jenks")$brks

#Step2: One map each, side by side with patchwork
g1 <- ggplot(sta) + geom_sf(aes(fill = homicides)) +
  scale_fill_viridis_b(breaks = unique(b1), direction = -1) + ggtitle("Count")
g2 <- ggplot(sta) + geom_sf(aes(fill = rate)) +
  scale_fill_viridis_b(breaks = unique(b2), direction = -1) + ggtitle("Per 100,000")

g1 | g2

Count and Rate Compared

But Rates Have Their Own Trap

#Step1: Read the municipio layer -- 2,469 units instead of 32
mun <- st_read("../data/mexico_homicides.gpkg", layer = "municipios", quiet = TRUE)

#Step2: The five worst rates in the country
mun %>% st_drop_geometry() %>% arrange(desc(rate)) %>%
  transmute(municipio_name, state_name, pop,
            homicides = round(homicides, 1), rate = round(rate)) %>% head(5)
           municipio_name state_name  pop homicides rate
1             Doctor Coss Nuevo León 1360      10.0  735
2                 Oquitoa     Sonora  496       1.7  336
3     Sitio de Xitlapehua     Oaxaca  713       1.7  234
4 Santa María Mixtequilla     Oaxaca 4690       9.3  199
5 General Enrique Estrada  Zacatecas 6644      13.0  196

The Small-Numbers Problem

Doctor Coss, Nuevo León has the worst homicide rate in Mexico: 735 per 100,000 — eight times Colima.

It has 1,360 inhabitants and averaged 10 homicides a year. One more incident moves the rate by 74 points.

Warning

A rate with a tiny denominator is unstable, not extreme. Ranking municipios by raw rate puts the smallest ones on top almost every time — and a map of it is a map of small populations.

What To Do About It

  • Pool years — we already used a three-year mean rather than one year
  • Suppress small denominators — grey out units below a population threshold and say so in the caption
  • Aggregate — report at state level where municipio counts are thin
  • Smooth — empirical Bayes shrinks unreliable rates toward the regional mean

The map on the next slide does the simplest of these: it shows municipios, but greys out those under 10,000 people.

Municipios, With Small Ones Suppressed

#Step1: Flag the units we do not trust
mun <- mun %>% mutate(rate_shown = ifelse(pop < 10000, NA, rate))

#Step2: na.value gives the suppressed units their own colour
ggplot(mun) +
  geom_sf(aes(fill = rate_shown), colour = NA) +
  scale_fill_viridis_b(breaks = c(0, 10, 25, 50, 100), direction = -1,
                       na.value = "grey85", name = "per 100,000") +
  coord_sf(crs = 6372) + theme_void()

Municipios, With Small Ones Suppressed

Part E: Making It Publishable

Everything a Reader Needs

brks <- classIntervals(sta$rate, n = 5,
                       style = "jenks")$brks

ggplot(sta) +
  geom_sf(aes(fill = rate),
          colour = "white",
          linewidth = 0.15) +
  scale_fill_viridis_b(
    breaks = unique(brks),
    direction = -1,
    name = "per 100,000") +
  coord_sf(crs = 6372) +
  labs(
    title = "Homicide is concentrated in the north-west and the Pacific centre",
    subtitle = "Intentional homicides per 100,000, annual mean 2021-2023",
    caption = "Sources: SESNSP; INEGI Censo 2020") +
  theme_void(base_size = 9)

The Checklist

  • Title states the finding, not the variable name — “Homicide is concentrated in…”, not “Homicide rate by state”
  • Legend has units — per 100,000, not rate
  • Sources in the caption, always, and the period in the subtitle
  • theme_void() — latitude and longitude gridlines are chart junk on a thematic map
  • coord_sf(crs = 6372) — never publish a Mexican map in Web Mercator

Part F: Interactive Maps with leaflet

Why Interactive?

A static map of 32 states cannot label them all, and a map of 2,469 municipios certainly cannot.

An interactive one does not have to — the reader hovers over their own state and reads the number.

Warning

leaflet needs unprojected coordinates. It expects EPSG:4326 and reprojects to Web Mercator itself. Hand it EPSG:6372 and you get an empty map.

A Minimal leaflet Map

#Step1: Make sure we are in WGS84 -- leaflet requires it
sta_ll <- st_transform(sta, 4326)

#Step2: A basemap and the polygons
leaflet(sta_ll) %>%
  addTiles() %>%
  addPolygons()

Three lines and you have a pannable, zoomable map. Now we make it say something.

Adding Colour, Labels, and a Legend

#Step1: A palette function built from our jenks breaks
pal <- colorBin("viridis", domain = sta_ll$rate,
                bins = unique(brks), reverse = TRUE)

#Step2: An HTML label per state
labels <- sprintf("<strong>%s</strong><br/>%.1f per 100,000<br/>%s homicides/year",
                  sta_ll$state_name, sta_ll$rate,
                  format(round(sta_ll$homicides), big.mark = ",")) %>%
  lapply(htmltools::HTML)

#Step3: Put it together. The ~ means "look this column up in the data".
leaflet(sta_ll) %>%
  addTiles() %>%
  addPolygons(fillColor = ~pal(rate), fillOpacity = 0.8,
              color = "white", weight = 1,
              highlightOptions = highlightOptions(weight = 3, color = "#333"),
              label = labels) %>%
  addLegend(pal = pal, values = ~rate,
            title = "per 100,000", position = "bottomright")

The Finished Map

Hover over a state. That is the whole point.

Do Not Forget This

A leaflet map is JavaScript. In a Quarto HTML document it depends on external files — unless you say otherwise:

---
format:
  html:
    embed-resources: true
---

Important

This is the trap from Session 1, and this is the week it bites. Without embed-resources: true, your socio formador opens the file and sees a blank white rectangle.

ggplot or leaflet?

ggplot

  • Goes in a PDF or a printed report
  • Full control of design
  • Works in a slide deck
  • Small file

leaflet

  • Reader explores it themselves
  • Basemap gives instant context
  • HTML only
  • Adds ~1 MB to the file

Tip

Put both in your reto report: a static map for the printed summary, an interactive one for the people who want to find their own municipio.

What Have We Learned?

  • A choropleth is three decisions: variable, colour, and breaks
  • Join Mexican data on CVEGEO, never on place names
  • Check your national totals before you draw anything
  • viridis is perceptually uniform, colourblind-safe, and prints in grey — turbo is not
  • Classification changes the story: equal interval put 21 of 32 states in one class; quantile spread them 7-6-6-6-7
  • A count map is a population map in disguise — México is 3rd by count, 13.5 by rate; Colima is 1st by rate with a quarter the homicides
  • But rates on small denominators are unstable: Doctor Coss, population 1,360, “leads” the country at 735

Next session: spatial joins — doing this by location rather than by code — plus buffers and intersections.