Lecture 6

Map visualization with geom_sf(): wildfire smoke

Byeong-Hak Choe

bchoe@geneseo.edu

SUNY Geneseo

September 30, 2026

🗺️ Draw maps

📁 Files for this lecture

Use read_sf() for GeoJSON. Keep county codes, site IDs, and acquisition times as text when reading the CSVs.

# Boundaries and shapes (GeoJSON)
counties <- read_sf("https://bcdanl.github.io/data/wildfire_us_counties_2020.geojson")
states <- read_sf("https://bcdanl.github.io/data/wildfire_us_states_2020.geojson")
provinces <- read_sf("https://bcdanl.github.io/data/wildfire_canada_provinces_2021.geojson")
plumes <- read_sf("https://bcdanl.github.io/data/wildfire_hms_smoke_2023_06.geojson")
nbac <- read_sf("https://bcdanl.github.io/data/wildfire_nbac_burned_area_2023.geojson")

# Data tables (CSV)
smoke_day <- read_csv(
  "https://bcdanl.github.io/data/wildfire_smoke_county_day_2023_06.csv",
  col_types = cols(GEOID = col_character())
)
smoke_month <- read_csv(
  "https://bcdanl.github.io/data/wildfire_smoke_county_month_2006_2025.csv",
  col_types = cols(GEOID = col_character())
)
ember <- read_csv(
  "https://bcdanl.github.io/data/wildfire_ember_county_day_2023_06.csv",
  col_types = cols(GEOID = col_character())
)
aqs <- read_csv(
  "https://bcdanl.github.io/data/wildfire_aqs_pm25_2023_06.csv",
  col_types = cols(GEOID = col_character(), site_id = col_character())
)
fires <- read_csv(
  "https://bcdanl.github.io/data/wildfire_firms_canada_2023.csv",
  col_types = cols(acq_time = col_character())
)

These ten files are the complete input set for today’s lecture. The GeoJSON files contain geometry; the CSV files provide smoke, monitor, and fire observations. Character column types preserve leading zeros in identifiers and acquisition times.

🗺️ Maps with geospatial data

• A map shows where smoke, fires, and people are located.
• A familiar boundary can make a regional pattern easy to recognize.
• A chart often makes rankings and time trends easier to compare.

Lecture 8 began with election results in a dot plot before drawing a map. We take the same route with wildfire smoke. A map is useful when the locations and neighboring areas matter, but it is not automatically the clearest way to compare numbers.

📋 State smoke averages from county data

state_smoke_june7 <- smoke_day |>
  filter(date == as.Date("2023-06-07")) |>
  left_join(
    counties |>
      st_drop_geometry() |>
      select(GEOID, STUSPS, state_name, region, pop_2020),
    by = "GEOID"
  ) |>
  summarize(
    smoke_pm25 = weighted.mean(smoke_pm25, pop_2020),
    people = sum(pop_2020),
    .by = c(STUSPS, state_name, region)
  )

• Each state value weights its county estimates by population.

The statewide average describes a population-weighted county estimate, not what every resident inhaled.

📋 Select states and a national reference

top_state_smoke <- state_smoke_june7 |>
  filter(STUSPS != "DC") |>
  group_by(region) |>
  slice_max(smoke_pm25, n = 5) |>
  ungroup()
national_smoke <- weighted.mean(
  state_smoke_june7$smoke_pm25,
  state_smoke_june7$people
)

• The chart uses the five states with the highest averages in each Census region.
• The reference line uses the population-weighted national average.

We keep all four regions in the chart, rather than selecting only the highest states nationally.

🔵 Dot plot: set up the axes

p_state0 <- ggplot(
  top_state_smoke,
  aes(
    x = smoke_pm25,
    y = fct_reorder(state_name, smoke_pm25),
    color = region
  )
)

• fct_reorder() sorts states by their smoke average.
• The axes compare values without asking readers to compare areas on a map.

➕ Dot plot: add points and a reference line

p_state1 <- p_state0 +
  geom_vline(xintercept = national_smoke,
             color = "grey45",
             linewidth = 0.5) +
  geom_point(size = 2.4) +
  scale_color_brewer(palette = "Dark2") +
  theme_minimal(base_size = 13) +
  labs(x = "Smoke PM2.5 (µg/m³)", y = NULL,
       color = "Census region")

p_state1

• The line marks the population-weighted national average: 28.9 µg/m³.
• Dark2 gives the four regions distinct colors before we fine-tune them.

🎨 Dot plot: apply regional colors

region_colors <- c(
  Midwest = "#d79a2b", Northeast = "#b2182b",
  South = "#3f78a8", West = "#578b67"
)

p_state2 <- p_state1 +
  scale_color_manual(values = region_colors)
p_state2

• A named color vector replaces the preset palette and keeps each region tied to the same color.

🗂️ Dot plot: facet by Census region

p_state2 +
  facet_wrap(~ region, ncol = 2,
             scales = "free_y") +
  guides(color = "none")

• The shared x-axis keeps the regional values comparable.
• A map will show whether the high-smoke states form a geographic cluster.

🛰️ June 7, 2023, seen from above

Where did this smoke come from? What would you check before trusting this map?

Let the map sit for a moment. Students usually notice the dark counties in Pennsylvania and New Jersey, the red fire dots in Quebec, and the dashed plume outlines. Ask what each layer is (points, polygons), where each number comes from, and what could mislead: a plume outline is drawn from satellite images and may be high in the air; county fills are averages. By the end of class every student should be able to build a map like this one. The composite slide in section “Points, plumes, and layers” builds a simpler version from the same layers (fire detections, Heavy plumes, cities); this hook adds the county smoke fill and draws the plumes as dashed outlines.

🗺️ A map is geometry plus a projection

• Points: a fire detection, an air-quality monitor, a city
• Lines: a road, a river, a flight path
• Polygons: a county, a smoke plume, a burned area
• A projection flattens the round Earth onto a flat page, and every projection distorts something

Healy: a map is just shapes drawn in the right order on an x–y plane, plus a projection. The rest is ordinary ggplot.

Healy, Data Visualization, ch. 7: maps are polygons made of latitude–longitude points, drawn with a projection. Albers is the usual choice for the lower 48 states. The same logic applies to all our layers today.

📋 Plan for today

• Simple features and projections: read, draw, and flatten boundaries
• Choropleths and small multiples: join smoke data to county shapes and compare dates
• Hexagons, points, plumes, and spatial joins
• When a map is not the best chart

Everything today uses the tidyverse, sf, scales, and ggrepel. The classwork uses a different smoke episode (the Midwest, June 27–29, 2023) with the same tools.

📐 Draw maps with sf

📂 The course files

All files are at https://bcdanl.github.io/data/.

File	One row is
wildfire_us_counties_2020.geojson	a lower-48 county (3,108, incl. DC)
wildfire_us_states_2020.geojson	a state or DC (49)
wildfire_canada_provinces_2021.geojson	a province or territory (13)
wildfire_smoke_county_day_2023_06.csv	county × day, June 2023
wildfire_smoke_county_month_2006_2025.csv	county × month, 2006–2025
wildfire_aqs_pm25_2023_06.csv	air monitor × day, June 2023
wildfire_ember_county_day_2023_06.csv	county × day, modeled fire PM2.5
wildfire_hms_smoke_2023_06.geojson	smoke plume drawn by an analyst, June 2023
wildfire_firms_canada_2023.csv	satellite fire detection in Canada, 2023
wildfire_nbac_burned_area_2023.geojson	2023 Canadian fire, outer outline

Every file is small (the largest is about 15 MB). GeoJSON files store geometry; CSV files store numbers that we join to geometry by a key (GEOID for counties). The files cover fires, smoke, air quality, boundaries, and population only.

📁 Read boundaries with read_sf()

library(tidyverse)
library(sf)
library(scales)
library(ggrepel)

course_url <- "https://bcdanl.github.io/data/"

counties  <- read_sf(paste0(course_url, "wildfire_us_counties_2020.geojson"))
states    <- read_sf(paste0(course_url, "wildfire_us_states_2020.geojson"))
provinces <- read_sf(paste0(course_url, "wildfire_canada_provinces_2021.geojson"))

• read_sf() reads GeoJSON, shapefiles (.shp), and GeoPackages into a tibble.
• Census and Statistics Canada publish boundary files like these for free.

read_sf() is the tidyverse-friendly version of st_read(): quiet, returns a tibble, keeps strings as characters. In DANL 310 we read Census shapefiles; GeoJSON is one text file instead of four .shp/.dbf/.shx/.prj files.

🔍 An sf object is a tibble with a geometry column

counties |>
  select(GEOID, county_name, pop_2020) |>
  print(n = 3)

Simple feature collection with 3108 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -124.7627 ymin: 24.5213 xmax: -66.9499 ymax: 49.3844
Geodetic CRS:  WGS 84
# A tibble: 3,108 × 4
  GEOID county_name    pop_2020                     geometry
  <chr> <chr>             <dbl>           <MULTIPOLYGON [°]>
1 01001 Autauga County    58805 (((-86.9176 32.6642, -86.71…
2 01003 Baldwin County   231767 (((-87.9729 31.1627, -87.94…
3 01005 Barbour County    25223 (((-85.4164 31.6195, -85.74…
# ℹ 3,105 more rows

• sf = simple features: shapes stored as rows of a table.
• geometry is a list column: one shape per county.
• Geometry is sticky: select() keeps it.
• CRS (coordinate reference system): WGS 84 longitude and latitude.

Simple features is the common standard GIS software uses to store points, lines, and polygons; sf is R’s version of it. Point to three things in the header: 3,108 features (rows), the geometry type (MULTIPOLYGON), and the CRS (WGS 84, i.e., longitude and latitude in degrees). The geometry column came along even though we did not select it. Healy contrasts the older map_data() approach (thousands of rows per state, one per polygon corner) with sf, where each unit is one row, so joins keep one row per county.

🖊️ Your first geom_sf()

ggplot(counties) +
  geom_sf(fill = "#efe7d8",
          color = "#8b7767",
          linewidth = 0.1)

• No aes(x, y): geom_sf() finds the geometry column itself.
• The axes are longitude and latitude in degrees.

The land and county-boundary colors are chosen explicitly. The grey panel, graticule lines, and degree labels still show the default plot theme. Note that the U.S.–Canada border along the 49th parallel is a straight line. That is a sign of unprojected longitude–latitude coordinates. We come back to this in the projection section.

🎨 Set a map theme once

theme_set(
  theme_void(base_size = 15) +
    theme(
      legend.position = "bottom",
      legend.title.position = "top",
      legend.key.width = unit(2.4, "lines"),
      plot.title = element_text(face = "bold"),
      strip.text = element_text(size = rel(1)),
      plot.caption = element_text(color = "grey40", hjust = 0),
      plot.caption.position = "plot"
    )
)

• theme_void() drops axes, grid lines, and the grey panel.
• theme_set() applies it to every plot after this line.

Healy uses a custom theme_map(); theme_void() is the built-in equivalent. Longitude and latitude tick labels rarely help a reader of a thematic map. Line charts later in the deck add theme_minimal() back.

➕ Style and layer: states over counties

state_lines <- geom_sf(
  data = states, fill = NA,
  color = "grey25", linewidth = 0.25
)

ggplot(counties) +
  geom_sf(fill = "grey92",
          color = "white",
          linewidth = 0.05) +
  state_lines

• Layers are drawn in order; fill = NA makes states see-through.
• A saved layer (state_lines) can be added to any later map.

linewidth replaces the old size argument for polygon borders. Very thin white county borders keep 3,108 counties from turning into a grey smudge. We reuse state_lines on almost every U.S. map today.

🔗 Dissolve counties into Census regions

regions <- counties |>
  group_by(region) |>
  summarize(pop = sum(pop_2020))

ggplot(regions) +
  geom_sf(aes(fill = region),
          color = "white", linewidth = 0.3) +
  scale_fill_manual(values = region_colors) +
  geom_sf_label(
    aes(label = paste0(region, "\n",
      number(pop / 1e6, accuracy = 0.1),
      "M")),
    size = 3.2,
    color = "#243b53", fill = "#fffaf0"
  ) +
  guides(fill = "none")

• summarize() on an sf object unions the polygons in each group.

In DANL 310 we wrote summarize(geometry = st_union(geometry)) to merge New York counties into a state outline. With group_by(), the same step makes one polygon per region and sums population. Census 2020 populations: South 126.3 M, West 76.4 M, Midwest 69.0 M, Northeast 57.6 M. The thin white outline shows only the four region borders: no county lines, slivers, or hairline gaps are left inside a region. That works because the course county file was simplified so that neighboring counties still share each border exactly. If each county had been simplified on its own, neighbors would no longer meet perfectly, the union would leave many tiny gaps, and any outline color would trace them.

🌎 North America: Different Projections

🧭 Which coordinate reference system?

st_crs(counties)$Name

[1] "WGS 84"

st_crs(counties)$epsg

[1] 4326

Code	Name	Use today
EPSG:4326	WGS 84 longitude/latitude	how the files are stored
EPSG:4269	NAD83 longitude/latitude	many U.S. federal files
ESRI:102003	USA Contiguous Albers Equal Area	lower-48 maps
EPSG:5070	NAD83 / Conus Albers (meters)	lower-48 data stored in meters
ESRI:102009	North America Lambert Conformal Conic	Canada + U.S. maps

A CRS code has two parts: the catalog (EPSG or ESRI) and the number. epsg.io lists thousands. Longitude–latitude is great for storing data and poor for drawing large areas, because a degree of longitude shrinks as you move north.

🌐 Longitude–latitude vs. Albers

us_map <- ggplot(states) +
  geom_sf(fill = "grey80", color = "white")

us_map + coord_sf(crs = 4326)

us_map +
  coord_sf(crs = "ESRI:102003")

• Albers (right) is equal-area: no region gets extra ink.

Left: in longitude–latitude the 49th parallel is a straight line and the country looks stretched. Healy recommends Albers for U.S. maps. Equal-area matters for choropleths: color covers area, so an inflated area gets extra visual weight.

📐 coord_sf() vs. st_transform()

us_albers <- coord_sf(crs = "ESRI:102003")

counties_5070 <- st_transform(counties, crs = 5070)
st_crs(counties_5070)$Name

[1] "NAD83 / Conus Albers"

st_bbox(counties_5070) |> round()

    xmin     ymin     xmax     ymax 
-2356011   269579  2258155  3172569 

• coord_sf(): draw in a projection; the data stay in degrees.
• st_transform(): changes the stored coordinates (here to meters); layers must share a CRS before a spatial join.

geom_sf() reprojects every layer to the coordinate system of the plot, so layers in different CRSs still line up. Spatial joins such as st_join() require both inputs to share a CRS; st_transform() is how you make that true. We save us_albers as an object and add it to U.S. maps. You do not need st_transform() just to measure: sf measures areas and distances on the sphere for longitude–latitude data, so st_area() and st_distance() return meters without a transform (the “containers” slide uses st_area() this way). All of today’s files are stored in EPSG:4326, so our joins need no transform either. The bounding box above is in meters, which shows the data themselves changed.

🌎 North America in one projection

na_lcc <- coord_sf(crs = "ESRI:102009")

ggplot() +
  geom_sf(data = provinces,
          fill = "#e7d9c4",
          color = "white") +
  geom_sf(data = states,
          fill = "grey80",
          color = "white") +
  na_lcc

• One coord_sf() reprojects both layers.

ESRI:102009 is a Lambert conformal conic projection centered on North America. It keeps shapes locally correct, which suits a map of smoke moving from Canada into the United States.

🔎 Zoom with longitude–latitude limits

east_zoom <- coord_sf(
  crs = "ESRI:102009",
  xlim = c(-95, -60), ylim = c(36, 56),
  default_crs = st_crs(4326)
)

ggplot() +
  geom_sf(data = provinces,
          fill = "#e7d9c4", color = "white") +
  geom_sf(data = states,
          fill = "grey80", color = "white") +
  east_zoom

• default_crs = st_crs(4326): “my limits are in degrees.”

Without default_crs, xlim and ylim are read in the plot’s CRS (meters for ESRI:102009). With it, we can type familiar longitudes and latitudes. The corners are projected, so the frame is slightly larger than the box you typed.

🎨 Choropleth Maps

📂 Read the smoke tables with GEOID as text

smoke_day <- read_csv(
  paste0(course_url, "wildfire_smoke_county_day_2023_06.csv"),
  col_types = cols(GEOID = col_character()))
smoke_month <- read_csv(
  paste0(course_url, "wildfire_smoke_county_month_2006_2025.csv"),
  col_types = cols(GEOID = col_character()))
aqs <- read_csv(
  paste0(course_url, "wildfire_aqs_pm25_2023_06.csv"),
  col_types = cols(GEOID = col_character(), site_id = col_character()))
ember <- read_csv(
  paste0(course_url, "wildfire_ember_county_day_2023_06.csv"),
  col_types = cols(GEOID = col_character()))

• GEOID is a FIPS code: 2-digit state + 3-digit county. "01001" is Autauga County, AL.

PM2.5 = tiny particles (2.5 micrometers or smaller) in the air, measured in micrograms per cubic meter (µg/m³). total_pm25 counts particles from every source; smoke_pm25 is the part above a normal background on days a satellite saw smoke.

readr happens to guess character for GEOID here because the first rows have a leading zero, but read.csv(), Excel, and a file sorted differently would all turn “01001” into the number 1001. Say the type explicitly.

⚠️ What a numeric GEOID does to a join

bad_ids <- as.character(as.numeric(smoke_day$GEOID))
head(unique(bad_ids), 3)

[1] "1001" "1003" "1005"

sum(!counties$GEOID %in% bad_ids)

[1] 287

• 287 counties in six states (AL, AZ, AR, CA, CO, CT) would silently drop out of the join.
• Join on codes, not names, and keep codes as text (Healy).

The six states have FIPS codes 01 to 09, so their county codes start with a zero. Joining a numeric GEOID to a character GEOID fails with a type error; converting it back to text without the zero fails silently, which is worse. Names are not keys either: 30 counties are called “Washington County”. If a file has already lost its zeros, mutate(GEOID = str_pad(GEOID, width = 5, pad = "0")) puts them back, a key repair like the France/Norway fix in DANL 310.

🔍 Check the keys before you join

n_distinct(smoke_day$GEOID)

[1] 3108

anti_join(counties, smoke_day, by = "GEOID") |> nrow()

[1] 0

setdiff(smoke_day$GEOID, counties$GEOID)

character(0)

• Every county has smoke data, and every smoke row has a county.
• anti_join() returns the rows with no match; we want zero.

Boundaries change. Census data from 2022 on use 9 Connecticut planning regions (09110–09190) in place of the 8 old counties (09001–09015) in our 2020 shapes. Joining newer Connecticut data would leave all 8 counties unmatched, and this anti_join() would show them. Match the boundary year to the data year (Healy).

In DANL 310 the World Bank join lost France and Norway because their ISO codes were missing in Natural Earth. The anti_join() check would have caught it before the map did. Healy warns that county boundaries and codes change over time, and Connecticut is the standard example: the state ended county government in 1960, and the Census Bureau began using its nine planning regions as county equivalents in 2022, with new FIPS codes. Older county data cannot be joined directly to newer shapes (and vice versa). Our shapes and smoke tables both use the 2020 counties, so the check returns zero here.

🔗 Join the table to the geometry

june7 <- smoke_day |>
  filter(date == as.Date("2023-06-07"))

smoke_june7 <- counties |>
  left_join(june7, by = "GEOID")

class(smoke_june7)

[1] "sf"         "tbl_df"     "tbl"        "data.frame"

nrow(smoke_june7)

[1] 3108

• Put the sf object on the left; the result stays sf, one row per county.

If the tibble is on the left, left_join() returns a tibble with a geometry column but without the sf class; st_as_sf() turns it back into an sf object.

🔵 Fill counties by smoke category

smoke_category <- smoke_june7 |>
  mutate(smoke_band = case_when(
    smoke_pm25 < 15 ~ "<15",
    smoke_pm25 <= 35 ~ "15–35",
    TRUE ~ ">35"
  ) |> factor(levels = c(
    "<15", "15–35", ">35"
  )))

p_band0 <- ggplot(
  smoke_category,
  aes(fill = smoke_band)
) +
  geom_sf(color = NA) + state_lines +
  scale_fill_brewer(
    palette = "YlOrRd",
    name = "Smoke PM2.5 (µg/m³)"
  ) +
  us_albers +
  theme(legend.position = "right")
p_band0

• fill maps an ordered smoke category to color.
• The category cutoffs are in µg/m³ of the smoke component.
• YlOrRd orders the colors from pale yellow to dark red.

🎨 Apply colors with scale_fill_manual()

smoke_band_colors <- c(
  "<15" = "#fff7bc",
  "15–35" = "#fdae61",
  ">35" = "#b2182b"
)

p_band1 <- p_band0 +
  scale_fill_manual(
    values = smoke_band_colors,
    name = "Smoke PM2.5 (µg/m³)"
  )
p_band1

• The named vector replaces the preset palette and assigns colors by category label.

🔖 Move the legend below the map

p_band2 <- p_band1 +
  guides(fill = guide_legend(nrow = 1)) +
  theme(legend.position = "bottom")
p_band2

• guide_legend(nrow = 1) keeps the three categories in one row.
• The legend states both the variable and its units.

🔄 Reversing colors changes the message

reversed_colors <- setNames(
  rev(unname(smoke_band_colors)),
  names(smoke_band_colors)
)

p_band0 +
  scale_fill_manual(
    values = reversed_colors,
    name = "Smoke PM2.5 (µg/m³)"
  )

• The data stayed the same. The dark color now means low smoke.
• Use a color order that matches what the legend asks readers to notice.

📈 A first choropleth: June 7, 2023

p_smoke0 <- ggplot(smoke_june7,
                   aes(fill = smoke_pm25)) +
  geom_sf(color = NA) + state_lines +
  scale_fill_viridis_c(
    option = "inferno", direction = -1,
    name = "Smoke PM2.5 (µg/m³)"
  ) +
  us_albers +
  labs(caption = "June 7, 2023 · CSU")
p_smoke0

Most of the country looks the same pale color. Why?

Smoke PM2.5 is the part of fine-particle pollution attributed to smoke (total minus a normal background, on days with a satellite smoke plume). The reversed Inferno palette makes high smoke dark, but its scale still stretches from the minimum to the maximum (229), so most counties get nearly the same pale color.

🎨 Set a sequential color gradient

p_smoke0 +
  scale_fill_gradient(
    low = "#fff7bc", high = "#b2182b",
    name = "Smoke PM2.5 (µg/m³)"
  )

• scale_fill_gradient() controls the low and high colors.
• The few extreme counties still compress most of the scale.

🚫 An extreme value stretches the gradient

p_smoke0 +
  scale_fill_gradient(
    low = "#fff7bc", high = "#b2182b",
    limits = c(0, 75), oob = squish,
    name = "Smoke PM2.5 (µg/m³)"
  ) +
  labs(caption = paste(
    "June 7, 2023 · values above 75",
    "share one color"
  ))

• Capping the scale reveals mid-range variation but hides differences above 75.
• The next map uses labeled bins so readers can see the cutoffs.

📊 Most counties sit at the bottom of the scale

ggplot(june7, aes(x = smoke_pm25)) +
  geom_histogram(binwidth = 5,
                 fill = "grey40") +
  theme_minimal(base_size = 13) +
  labs(x = "Smoke PM2.5 (µg/m³), June 7",
       y = "Counties")

• Median 9.7, maximum 229 µg/m³: a long right tail.
• A few small negatives come from subtracting the background from a noisy total: measurement noise, not “negative smoke.”

The top counties on June 7 were all in eastern Pennsylvania (Lehigh 229, Carbon 214, Northampton 212). Negative values come from subtracting a background from a noisy total; they are why a log or square-root scale would fail here.

🎨 Bin the scale at meaningful breaks

pm25_bins <- scale_fill_fermenter(
  palette = "YlOrRd", direction = 1,
  breaks = c(5, 15, 35, 75, 150),
  name = "PM2.5 (µg/m³)"
)

ggplot(smoke_june7) +
  geom_sf(aes(fill = smoke_pm25),
          color = NA) +
  state_lines + pm25_bins + us_albers +
  labs(caption =
    "Smoke PM2.5, June 7, 2023 · CSU")

• 35 µg/m³ is the level of EPA’s 24-hour standard for total PM2.5; here it is only a reference for the smoke part. Negatives fall in the lowest bin.
• We save the scale as pm25_bins and reuse it for every PM2.5 map.

About 74.8 million people lived in counties with smoke above 35 on June 7, and 41.6 million in counties above 75. Breaks tied to a standard are easier to explain than breaks chosen by software.

🔄 Different bins, different story

q_breaks <- quantile(june7$smoke_pm25,
                     c(0.2, 0.4, 0.6, 0.8))

ggplot(smoke_june7) +
  geom_sf(aes(fill = smoke_pm25),
          color = NA) +
  state_lines + us_albers +
  scale_fill_fermenter(
    palette = "YlOrRd", direction = 1,
    breaks = q_breaks,
    labels = label_number(0.1),
    name = "Smoke PM2.5 (quintiles)"
  )

Quintiles always put 40% of counties in the top two bins. Was 40% of the country smoky?

Healy shows that the choice of bins changes the picture. Quantile bins spread colors evenly across counties, which hides how extreme the Northeast was: the top quintile starts near 22 µg/m³, well below the 35 standard.

↔︎️ A diverging scale needs a meaningful midpoint

june_diff <- smoke_month |>
  filter(month == 6) |>
  group_by(GEOID) |>
  summarize(diff = smoke_pm25[year == 2023] -
    mean(smoke_pm25[year <= 2022]))

counties |>
  left_join(june_diff, by = "GEOID") |>
  ggplot(aes(fill = diff)) +
  geom_sf(color = NA) + state_lines +
  scale_fill_gradient2(midpoint = 0,
    low = "#2166ac", high = "#b2182b",
    name = "Change vs. usual June (µg/m³)") +
  us_albers

• White (0) means “a typical June”: the West was near normal; the East far above.

“Usual June” is the 2006–2022 average of each county’s June smoke. Use a diverging scale only when the middle means something (zero difference, a standard, a national average). Here the values run from about −1.6 to 26.7 µg/m³. Because scale_fill_gradient2() is symmetric around 0, the few small negatives come out almost white, and the map is nearly all white-to-red because June 2023 was above normal almost everywhere. That is honest; stretching the blue side to its own minimum would exaggerate the negatives. 8% of counties were below their usual June.

⚖️ Counts vs. rates: smoky days

june23 <- smoke_month |>
  filter(year == 2023, month == 6)

days_map <- counties |>
  left_join(june23, by = "GEOID") |>
  mutate(person_days = days_gt15 * pop_2020)

ggplot(days_map, aes(fill = days_gt15)) +
  geom_sf(color = NA) + state_lines +
  scale_fill_viridis_c(option = "magma",
    direction = -1,
    name = "Days with smoke > 15 µg/m³") +
  us_albers

• Rate: days_gt15 = smoky days for a typical resident, so small and large counties compare fairly. Up to 17 days (Michigan, Wisconsin).
• Count: person_days = days × population, a total (next slide).

days_gt15 counts the days in June 2023 with smoke PM2.5 above 15 µg/m³. It describes the typical person in each county, so it is comparable across small and large counties. person_days multiplies it by population. Here the rate is days_gt15 (smoky days for a typical resident); the count is person_days, a county total.

🏙️ Inside the smoke, totals follow population

ggplot(days_map,
       aes(fill = person_days / 1e6)) +
  geom_sf(color = NA) + state_lines +
  scale_fill_viridis_c(option = "magma",
    direction = -1,
    name = "Person-days > 15 (millions)") +
  us_albers

• Big cities light up wherever smoke reached: Healy’s ur-choropleth.
• 9 of the 10 darkest counties are among the 50 most populous; smoke-free Los Angeles stays pale.

The darkest county is Cook County, IL (9 days × 5.3 million people), then Wayne, MI (Detroit) and Queens, NY. The one top-10 exception is Milwaukee County, WI (57th by population, 15 smoky days). Person-days are useful for totals (“how many person-days of exposure?”), but inside the smoky Midwest and East a map of them mostly shows big cities. Counties with zero smoky days have zero person-days however many people live there: 1,051 counties, including Los Angeles, Harris (TX), Maricopa, San Diego, and Orange, 5 of the 6 largest. So totals follow population only where the smoke went, and the smoky-days map on the previous slide does not follow population at all. If students ask for a number: among the 2,057 counties with at least one smoky day, the correlation with population (1 = move together, 0 = unrelated) is 0.91 for person-days and 0.06 for days. Healy’s ur-choropleth point is broader than totals: many county maps, even some rates, end up showing where people live; totals such as person-days are the most extreme case.

❔ Missing is not zero

monitor_june7 <- aqs |>
  filter(date == as.Date("2023-06-07")) |>
  group_by(GEOID) |>
  summarize(pm25 = mean(pm25))

counties |>
  left_join(monitor_june7, by = "GEOID") |>
  ggplot(aes(fill = pm25)) +
  geom_sf(color = NA) + state_lines +
  scale_fill_fermenter(palette = "YlOrRd",
    direction = 1, na.value = "grey85",
    breaks = c(5, 15, 35, 75, 150),
    name = "Monitor PM2.5 (µg/m³)") +
  us_albers

• Only 506 of 3,108 counties had a monitor reporting. Grey means no data, not clean air.

Never replace NA with 0 to “fill in” a map. The same rule applies to other public data that suppress small counts for privacy: a suppressed count is unknown, not zero (Healy discusses this problem in county maps).

📅 Small-Multiple Maps

📅 Small multiples by day

days_5_10 <- smoke_day |>
  filter(between(date,
    as.Date("2023-06-05"),
    as.Date("2023-06-10")))

counties |>
  left_join(days_5_10,
            by = "GEOID") |>
  ggplot(aes(fill = smoke_pm25)) +
  geom_sf(color = NA) +
  state_lines +
  pm25_bins +
  facet_wrap(~ date, ncol = 3) +
  labs(caption =
    "Smoke PM2.5 · CSU") +
  coord_sf(crs = "ESRI:102003",
    xlim = c(-92, -67),
    ylim = c(35, 48),
    default_crs = st_crs(4326))

The joined object has 3,108 × 6 = 18,648 rows: one county shape per day. Facets share one scale by default, which is what we want: the same color means the same smoke level in every panel. The worst smoke moves from the New York–Pennsylvania border (around Binghamton: Tioga, Broome, Chemung, and Tompkins in NY; Wyoming, Lackawanna, Susquehanna, and Bradford in PA, all about 87–89 µg/m³) on June 6 to eastern Pennsylvania and New Jersey on June 7, then south toward Maryland, Delaware, and DC on June 8.

🗂️ Small multiples by year

yearly_days <- smoke_month |>
  filter(year >= 2017) |>
  summarize(days = sum(days_gt15),
            .by = c(GEOID, year))

counties |>
  left_join(yearly_days,
            by = "GEOID") |>
  ggplot(aes(fill = days)) +
  geom_sf(color = NA) +
  scale_fill_fermenter(
    palette = "Reds",
    direction = 1,
    breaks = c(1, 5, 10, 20, 40),
    name = paste(
      "Days with smoke > 15",
      "µg/m³"
    )) +
  facet_wrap(~ year, ncol = 3) +
  labs(caption =
    "Smoke PM2.5 · CSU") +
  us_albers

2020 and 2021 stand out in the West (California, Oregon, and the northern Rockies); 2023 stands out in the East and Midwest. Weighted by population, the average person had 9.2 days above 15 in 2023, about twice the next highest year (2021, 4.4 days). 2025 also stands out in the upper Midwest and Northeast (4.2 days per person nationally, third-highest since 2006), and 2017–2018 were bad years in the West. summarize(.by = c(GEOID, year)) groups for this one step only, so there is no “regrouping” message.

⚠️ Cautions for small-multiple maps

• Keep one shared scale across panels (the facet_wrap() default).
• Each panel is small, so fine detail disappears.
• Big, sparsely populated western counties dominate the eye.
• To compare timing in a few places, a line chart often works better.

Statebins draw every state as a same-sized tile; cartograms can resize states by population. Both make small, crowded states more visible while changing their geographic shapes or positions (Healy).

Healy makes the same points with small-multiple state maps by year: they show broad change, but the unequal sizes of the units make the picture hard to read. The statebins package (Bob Rudis; used in Healy’s chapter) draws these tile maps; it is not installed for this course and is optional further reading.

🔷 Hexbin Maps

🔷 A regular grid changes the map’s units

• County boundaries vary greatly in size. A hex grid uses cells of equal area.
• The map still preserves location and neighboring cells.
• We assign each county to the cell containing a point inside its boundary.

Lecture 8 uses statebins as an alternative to state choropleths. Here we make an sf hex grid so every map layer can still be drawn with geom_sf(). A hex is a geographic area, not an equal-population tile. This view reduces the visual dominance of huge county outlines but does not make each person equally visible.

📐 Build an sf hex grid

states_m <- st_transform(states, 5070)
hex_geom <- st_make_grid(
  states_m, cellsize = 180000,
  square = FALSE
)
hex_grid <- st_sf(
  hex_id = seq_along(hex_geom),
  geometry = hex_geom
)

• EPSG:5070 lets us set the grid spacing to 180 km.
• The grid consists of equal-area cells in this projection.

📍 Assign counties to hexagons

county_points <- smoke_june7 |>
  st_transform(5070) |>
  st_point_on_surface()

hex_counties <- county_points |>
  select(GEOID, smoke_pm25, pop_2020) |>
  st_join(hex_grid |> select(hex_id),
          join = st_within)

• st_point_on_surface() puts one representative point inside each county.
• st_join() assigns that point to one hexagon.

🔥 Summarize smoke within hexagons

hex_summary <- hex_counties |>
  st_drop_geometry() |>
  filter(!is.na(hex_id)) |>
  summarize(
    smoke_pm25 = weighted.mean(
      smoke_pm25, pop_2020
    ),
    counties = n(), .by = hex_id
  )

hex_smoke <- hex_grid |>
  left_join(hex_summary, by = "hex_id") |>
  filter(!is.na(smoke_pm25))

• Each cell summarizes its assigned counties using population weights.
• A county is represented in one cell, even when its area spans several cells.

🔥 Fill hexagons with smoke estimates

ggplot(hex_smoke) +
  geom_sf(aes(fill = smoke_pm25),
          color = "white",
          linewidth = 0.1) +
  geom_sf(data = states_m, fill = NA,
          color = "grey35",
          linewidth = 0.2) +
  pm25_bins + us_albers +
  labs(caption = paste(
    "June 7, 2023 · county estimates",
    "in hexagons"
  ))

• The same county estimates now appear in regular spatial units.
• State outlines give viewers a geographic reference.

⚖️ Hexagons and counties answer different questions

• County maps match the units in the source table and show local boundaries.
• Hexagons make area sizes consistent but discard county shapes.
• Neither map identifies the concentration that each person actually inhaled.

The hex grid is an aggregation choice, not a new measurement. Large counties can cover multiple hexagons, but the representative point assigns each county to just one. The same source smoke and population values underlie both maps. The grid can be useful for regional pattern recognition; use the county map when county identity matters.

🔁 Plotting Multiple Layers

📂 Read the fire and plume layers

fires <- read_csv(
  paste0(course_url, "wildfire_firms_canada_2023.csv"),
  col_types = cols(acq_time = col_character()))
nbac <- read_sf(paste0(course_url, "wildfire_nbac_burned_area_2023.geojson"))
plumes <- read_sf(paste0(course_url, "wildfire_hms_smoke_2023_06.geojson"))

• fires: one row per satellite fire detection (latitude, longitude, fire power).
• nbac: one outer outline per 2023 fire; plumes: one smoke plume per row.
• The outlines fill in unburned islands, so for burned area use area_ha, not st_area(): 14.8 vs. 16.9 million ha.

FIRMS: NASA MODIS detections, a roughly 1 km pixel that looked hot at the moment a satellite passed over. A detection is any hot pixel, so some may be farm or slash burning rather than wildfire. NBAC: the Canadian Forest Service’s final mapped fire perimeters. HMS: NOAA analysts draw smoke outlines by hand on satellite images.

📍 Fire detections become points

firms_sf <- fires |>
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

firms_sf |>
  select(local_date, frp, prov_abbr) |>
  print(n = 3)

Simple feature collection with 280457 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -140.7173 ymin: 42.102 xmax: -55.498 ymax: 81.8968
Geodetic CRS:  WGS 84
# A tibble: 280,457 × 4
  local_date   frp prov_abbr            geometry
  <date>     <dbl> <chr>             <POINT [°]>
1 2023-01-02     0 AB        (-117.9582 53.4929)
2 2023-01-02     0 AB        (-117.9781 53.4971)
3 2023-01-02     0 AB        (-117.9604 53.4998)
# ℹ 280,454 more rows

• coords: x (longitude) first, then y (latitude). crs = 4326: the numbers are degrees.

Swapping longitude and latitude is the most common mistake: the points land in Antarctica or off the map. frp is fire radiative power in megawatts, a rough measure of how intensely the fire burned. The first rows printed happen to be January detections in Alberta with frp = 0; only about 1% of detections have frp = 0, and the typical summer detection is around 40–60 MW.

🔥 The 2023 fire season

fires_season <- firms_sf |>
  filter(month(local_date) %in% 5:9) |>
  mutate(month = month(local_date,
                       label = TRUE))

ggplot() +
  geom_sf(data = provinces, fill = "grey92") +
  geom_sf(data = fires_season, size = 0.05,
          aes(color = month), alpha = 0.3) +
  scale_color_viridis_d(
    option = "turbo", end = 0.9
  ) +
  guides(color = guide_legend(override.aes =
    list(size = 3, alpha = 1))) +
  na_lcc + labs(color = NULL)

275,742 of the 280,457 detections fall in May–September. Alberta and British Columbia burned from May through September, the Northwest Territories peaked in July and August, and Quebec burned mostly in June. Small, transparent points keep 275,000 dots from becoming one blob; override.aes makes the legend keys readable.

📅 One week: June 1–7

fires_week <- firms_sf |>
  filter(between(local_date,
    as.Date("2023-06-01"),
    as.Date("2023-06-07")))

ggplot() +
  geom_sf(data = provinces,
          fill = "#e7d9c4", color = "white") +
  geom_sf(data = states,
          fill = "grey85", color = "white") +
  geom_sf(data = fires_week, size = 0.3,
          color = "#b2182b", alpha = 0.4) +
  east_zoom

• Quebec had 5,619 detections in these seven days, more than any other province.

local_date is the local solar date (UTC plus longitude/15 hours). Canada is roughly 4–9 hours behind UTC, so an evening satellite pass (about 8–10 p.m. local solar time) happens after midnight UTC and its acq_date is already the next day. About 84,000 detections (30%) are like this; early-morning and daytime passes keep the same date. In Quebec, 8,678 night detections have a UTC date one day later than local_date. Filtering on local_date keeps each detection on the local day it was seen: Quebec has 5,619 detections for June 1–7 by local_date and 5,565 by acq_date, and 3,345 of the 14,379 detections in this week have an acq_date one day later than their local_date. Many Quebec fires started around June 1 (NBAC lists most as natural causes).

🟫 Burned-area perimeters are polygons

ggplot() +
  geom_sf(data = provinces,
          fill = "grey92", color = "white") +
  geom_sf(data = nbac, aes(fill = cause),
          color = NA) +
  scale_fill_manual(values = c(
    Human = "#3f78a8",
    Natural = "#b2182b",
    Undetermined = "grey45")) +
  na_lcc + labs(fill = NULL)

• 755 of 968 mapped fires had a natural cause (mostly lightning).

The largest 2023 perimeter is a Quebec fire (fire_id 346, 969,354 ha, started June 1). The full 2023 NBAC record has about 14.8 million ha burned across 2,244 fires; our extract keeps the 968 fires with parts of at least 1 km². Detections (points) show where a satellite saw heat at one moment; perimeters show the outer outline of the final burned area. Unburned islands inside a fire are filled in, parts under 1 km² are dropped, and edges are simplified (500 m), so the shapes total about 16.9 million ha with st_area(), about 14% more than the official 14.8 million ha in area_ha. Use area_ha for burned area.

🌫️ Smoke plumes are overlapping polygons

plumes_june7 <- plumes |>
  filter(date == as.Date("2023-06-07")) |>
  mutate(density = factor(density, levels =
    c("Light", "Medium", "Heavy")))
plumes_by_density <- plumes_june7 |>
  group_by(density) |> summarize()

ggplot() +
  geom_sf(data = provinces, fill = "grey96") +
  geom_sf(data = plumes_by_density, color = NA,
          aes(fill = density), alpha = 0.7) +
  state_lines + scale_fill_manual(values = c(
    "#fdd49e", "#fc8d59", "#7f0000")) +
  east_zoom + labs(fill = NULL)

• group_by(density) |> summarize() merges overlapping plumes, so each color is one density level.

There are 88 plumes on June 7 (UTC date): 37 Heavy, 18 Medium, 33 Light. The 88 outlines overlap many times, and see-through layers stacked on each other would make the color show how many outlines overlap, not how dense the smoke is. summarize() merges the plumes of each density into one shape. Heavy outlines usually sit inside Medium and Light ones, and summarize() returns Light, Medium, Heavy in factor order, so Heavy is drawn last and ends up on top. state_lines goes last so borders stay visible through the smoke. Plumes are drawn from satellite images and are often high in the air, so a plume overhead is not the same as smoke at breathing height.

✏️ Make your own points: cities

cities <- tribble(
  ~city,          ~lat,    ~lon,
  "New York",     40.713, -74.006,
  "Philadelphia", 39.953, -75.165,
  "Montréal",     45.502, -73.567,
  "Toronto",      43.653, -79.383,
  "Québec City",  46.814, -71.208
) |>
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326)
ggplot(cities) +
  geom_sf(color = "#1c4982", size = 2) +
  geom_sf_label(aes(label = city),
    color = "#1c4982", fill = "#f7f2e9",
    hjust = 0, nudge_x = 0.4) +
  coord_sf(xlim = c(-80, -64))

Same idea as the New York cities example in DANL 310: type coordinates into a tribble(), then st_as_sf(). hjust = 0 starts each label at its point, and nudge_x = 0.4 moves it 0.4 degrees to the right (this map is not projected, so x is longitude), so no label hides its dot. coord_sf(xlim = c(-80, -64)) leaves room on the right so Québec City is not cut off; geom_sf_label() does not widen the panel to fit its labels. When students run this, geom_sf_label() warns that st_point_on_surface may not give correct results for longitude/latitude data; for points (and for the regions labels earlier) the warning is harmless at this scale. Without other layers the map is meaningless, which is the point of the next slide.

🔁 Fire → wind → smoke: put the layers together

heavy_june7 <- plumes_june7 |>
  filter(density == "Heavy")

ggplot() +
  geom_sf(data = provinces, fill = "grey95") +
  geom_sf(data = states, fill = "grey95") +
  geom_sf(data = heavy_june7, color = NA,
          fill = "grey30", alpha = 0.2) +
  geom_sf(data = fires_week, size = 0.3,
          color = "#b2182b", alpha = 0.4) +
  geom_sf(data = cities, size = 1.5,
          color = "#1c4982") +
  geom_label_repel(
    data = cities, size = 3, seed = 42,
    aes(label = city, geometry = geometry),
    stat = "sf_coordinates",
    color = "#1c4982", fill = "#f7f2e9"
  ) + east_zoom +
  labs(caption = paste0(
    "Red: fires, June 1–7 (FIRMS)\n",
    "Grey: Heavy plumes, June 7 UTC (HMS)"
  ))

• geom_label_repel(stat = "sf_coordinates") works like geom_sf_label() but pushes labels apart.

Read it as a story: fires in Quebec (red), Heavy smoke plumes drifting south (grey), and the cities under them. geom_label_repel() from ggrepel pushes labels apart so New York and Philadelphia do not collide; stat = “sf_coordinates” tells it to take x and y from the geometry column.

📍 Monitors and counties share one scale

monitors_june7 <- aqs |>
  filter(date == as.Date("2023-06-07")) |>
  st_as_sf(coords = c("longitude",
    "latitude"), crs = 4326)

ggplot() +
  geom_sf(data = smoke_june7, color = NA,
          aes(fill = total_pm25)) +
  state_lines +
  geom_sf(data = monitors_june7, shape = 21,
          aes(fill = pm25), size = 2) +
  pm25_bins + coord_sf(crs = "ESRI:102003",
    xlim = c(-81, -71), ylim = c(38.5, 44),
    default_crs = st_crs(4326))

• Fill is total PM2.5 for both: county estimates and monitor readings.
• Both layers map fill, so one binned scale colors counties and monitors alike.

County fill here is total PM2.5 (all sources), which is what a monitor measures. 752 monitors reported on June 7; the highest daily mean was 258 µg/m³ in Northampton County, PA. If the points used the colour aesthetic instead, scale_fill_fermenter(…, aesthetics = c(“fill”, “colour”)) shares one scale across both aesthetics.

🔗 Spatial Joins

🔎 Which counties touched a Heavy plume?

counties_heavy <- counties |>
  st_filter(heavy_june7)

ggplot() +
  geom_sf(data = counties_heavy,
          fill = "#b6533c", color = NA) +
  state_lines + us_albers

• st_filter() keeps counties whose shapes intersect a Heavy plume (June 7, UTC).
• nrow(counties_heavy): 721 counties touched, with 113.3 million residents in those counties.

st_filter() is filter() for geometry: the default test is st_intersects(), so even a small overlap marks the whole county. The population total counts all residents of those counties; it is not a count of people directly under the plume. HMS outlines smoke from every fire, not only Canadian ones: the red county in northern Arizona (Coconino) touches a small, separate Heavy plume near Flagstaff (about 220 km²) drawn over a local fire, not Canadian smoke. It shows up again as “Plume only” on the next slide.

🌫️ Plume outlines and surface estimates differ

exposure <- smoke_june7 |>
  mutate(
    touched = GEOID %in% counties_heavy$GEOID,
    group = case_when(
      touched & smoke_pm25 > 35 ~ "Both",
      touched ~ "Plume only",
      smoke_pm25 > 35 ~ "Surface > 35 only",
      TRUE ~ "Neither"))

ggplot(exposure, aes(fill = group)) +
  geom_sf(color = NA) + state_lines +
  scale_fill_manual(values = c(
    Both = "#7f2704", Neither = "grey90",
    "Plume only" = "#fdae6b",
    "Surface > 35 only" = "#3f78a8")) +
  guides(fill = guide_legend(nrow = 2,
    title = NULL)) + us_albers +
  labs(caption = paste0(
    "June 7, 2023. Plume: Heavy HMS plume\n",
    "(UTC date). Surface: smoke > 35 µg/m³"))

“Both” means a county touched a Heavy plume and its estimated county-level surface smoke was above 35. Named colors tie each color to its group, so renaming or reordering groups cannot shift the colors. The small orange patch in northern Arizona is a local Heavy plume near Flagstaff, not Canadian smoke. Orange counties touched a Heavy plume but had a county-level surface estimate at or below 35. Possible reasons include smoke aloft, a small overlap, different time windows, or uncertainty in either data source; the map alone cannot tell which.

👥 Count residents by county group

exposure |>
  st_drop_geometry() |>
  group_by(group) |>
  summarize(counties = n(),
            people_millions = sum(pop_2020) / 1e6)

# A tibble: 4 × 3
  group             counties people_millions
  <chr>                <int>           <dbl>
1 Both                   296           71.9 
2 Neither               2351          213.  
3 Plume only             425           41.4 
4 Surface > 35 only       36            2.89

• st_drop_geometry(): a plain table, so summarize() won’t merge shapes.
• Residents of plume-touched counties: 113.3 million. Residents of counties with estimated surface smoke > 35: 74.8 million.
• Both are county-level indicators, not counts of what each resident inhaled.

HMS plumes come from satellite images of the smoke column; they do not directly measure ground-level exposure. The CSU smoke estimate interpolates monitor readings and also uses HMS outlines, so the two indicators are not independent. Because their time windows and spatial definitions differ, this comparison cannot establish how many people inhaled heavy smoke. Without st_drop_geometry(), summarize() on an sf object would also union the county shapes in each group (as on the regions slide), which is slow and not needed for a table. The 36 “surface only” counties (in Virginia and North Carolina) sit just south of a Heavy plume outline; outlines are hand-drawn and approximate.

📍 Points in polygons with st_join()

fires_prov <- firms_sf |>
  mutate(id = row_number()) |>
  st_join(provinces |> select(prov_join = prov_abbr))

c(detections = nrow(firms_sf), rows = nrow(fires_prov))

detections       rows 
    280457     280457 

• st_join() returns one row for each point–polygon match. It adds the polygon’s columns to the point.
• Same number of rows here: no detection matched two provinces. A point inside two overlapping polygons would get two rows, so compare nrow() before and after every join.

mutate(id = row_number()) gives each detection an id so we can trace it after the join. A point-in-polygon join is a left join by location: points outside every polygon keep NA, and a point inside two overlapping polygons gets two rows. The course province file was simplified so that neighboring provinces still share each border exactly, so no detection matches two provinces: 280,457 rows for 280,457 detections. Layers that do overlap, such as the HMS smoke plumes, would add extra rows, so the nrow() check is a habit worth keeping.

🔍 Check the join before you count

fires_prov |>
  st_drop_geometry() |>
  summarize(
    matches = sum(!is.na(prov_join)),
    .by = id) |>
  count(matches)

# A tibble: 2 × 2
  matches      n
    <int>  <int>
1       0     22
2       1 280435

join_check <- fires_prov |>
  st_drop_geometry() |>
  count(check = case_when(
    is.na(prov_join) ~ "outside",
    prov_join == prov_abbr ~ "same",
    TRUE ~ "different"))
join_check

# A tibble: 3 × 2
  check          n
  <chr>      <int>
1 different    236
2 outside       22
3 same      280199

• matches: 0 = outside every polygon, 1 = one province. No 2s, so each row is one detection.
• check: does the polygon agree with the label prov_abbr? The few misses sit on borders and coasts.

Each detection has exactly one row, so counting rows counts detections. 280,199 detections (99.9% of those inside a polygon) agree with the provided label. Of the 236 that land in a different province, 228 are Northwest Territories detections at most 1.6 km north of the 60th parallel that match British Columbia (150), Alberta (67), or Saskatchewan (11); the other 8 sit on the Alberta–Saskatchewan, Ontario–Quebec, and Northwest Territories–Nunavut borders. Why just north of the line? sf connects the corner points of a longitude–latitude outline along the shortest path on the globe, not along a line of latitude. The course file stores the 60th-parallel border with a corner only every few degrees, so between corners the border bows slightly north of 60°N. The same bow along the U.S. border at the 49th parallel leaves 17 British Columbia and Manitoba detections, all within 0.9 km north of 49°N, outside every polygon (this layer has nothing south of the border). The other 5 of the 22 outside detections are on coasts or small islands (the course file drops islands under 250 km²). Points this close to a border can land on either side, so always check nrow() and matches per point before trusting counts from a spatial join.

📊 Count points per polygon, then map the counts

fire_counts <- fires_prov |>
  st_drop_geometry() |>
  count(prov_join, name = "detections")

provinces |>
  left_join(fire_counts,
    by = c("prov_abbr" = "prov_join")) |>
  mutate(detections =
    replace_na(detections, 0)) |>
  ggplot(aes(fill = detections)) +
  geom_sf(color = "white") +
  scale_fill_viridis_c(
    option = "inferno", labels = comma,
    name = "Fire detections, 2023"
  ) +
  na_lcc

• Here NA means no detections: a true zero (Prince Edward Island).

Nunavut is three times Alberta’s size but had 50 detections. Counts or detections per km²?

count() on the joined label, then left_join() the counts back onto the polygons. The one NA row in fire_counts holds the 22 detections that fell outside every province polygon; it matches no province, so it drops out. Prince Edward Island gets NA after the join because no detection fell inside it. That is a real zero, unlike the monitor map, where NA meant no measurement. So replace_na(0) is correct here and wrong there. Counts: Alberta 74,229, then British Columbia, the Northwest Territories, and Quebec. Per 10,000 km² (st_area() of these polygons) the order changes: Alberta about 1,124, British Columbia 761, Northwest Territories 491, Saskatchewan 307, Quebec 214. Nunavut is about 0.2. Bigger provinces have more room for fires, which is the same count-versus-rate question as person-days.

⚖️ From counties to states with population weights

state_june7 <- smoke_june7 |>
  st_drop_geometry() |>
  group_by(STUSPS) |>
  summarize(smoke = weighted.mean(
    smoke_pm25, pop_2020))

states |>
  left_join(state_june7, by = "STUSPS") |>
  ggplot(aes(fill = smoke)) +
  geom_sf(color = "white",
          linewidth = 0.2) +
  pm25_bins + us_albers +
  labs(caption = "Smoke PM2.5, June 7, 2023")

• Weighting by population gives the smoke level of the average resident.

Top states on June 7: New Jersey 143, then Pennsylvania, Delaware, and New York. A state average hides the variation among its counties, which is why the county map is more informative for this episode.

🇨🇦 Modeled share of fire smoke from Canada

ember_june <- ember |>
  group_by(GEOID) |>
  summarize(all = mean(pm25_all_fires),
    canada = mean(pm25_canada_fires)) |>
  mutate(share = canada / all)

counties |>
  left_join(ember_june, by = "GEOID") |>
  ggplot(aes(fill = share)) +
  geom_sf(color = NA) + state_lines +
  scale_fill_viridis_c(
    option = "inferno", limits = c(0, 1),
    labels = percent,
    name = "Canadian share of fire PM2.5"
  ) +
  us_albers + labs(caption = "EPA CMAQ model")

• Median county share: 91%. This is a model, not a measurement.

EPA’s EMBER project ran an air-quality model three ways: with all emissions, with no fires, and with no Canadian fires. The differences give modeled PM2.5 from all fires and from Canadian fires. The share is above 75% across most of the Midwest and Northeast and lower across the West (Pacific and Mountain counties about one-third) and in Texas and Louisiana, where other fires contributed more.

📊 Is the data really spatial?

📈 Timing is easier to read in a line chart

regional_daily <- smoke_day |>
  left_join(st_drop_geometry(counties),
            by = "GEOID") |>
  group_by(region, date) |>
  summarize(smoke = weighted.mean(
    smoke_pm25, pop_2020), .groups = "drop")

ggplot(regional_daily,
       aes(date, smoke, color = region)) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = region_colors) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top") +
  labs(x = NULL, color = NULL,
       y = "Smoke PM2.5 (µg/m³)")

• Northeast peak: June 7 (104 µg/m³). Midwest peak: June 28 (88).

Healy’s “Is your data really spatial?” section: when the question is about change over time, a line per region can say more than 30 maps. The West stayed below 1 µg/m³ all month.

🗺️ Counties are containers, not the smoke itself

june7_area <- smoke_june7 |>
  # st_area() gives m² on the sphere; / 1e6 gives km²
  mutate(area_km2 = as.numeric(st_area(geometry)) / 1e6) |>
  st_drop_geometry()

june7_area |>
  summarize(
    people_weighted = weighted.mean(smoke_pm25, pop_2020),
    area_weighted   = weighted.mean(smoke_pm25, area_km2),
    simple_mean     = mean(smoke_pm25))

# A tibble: 1 × 3
  people_weighted area_weighted simple_mean
            <dbl>         <dbl>       <dbl>
1            28.9          10.3        18.0

• Smoke is a continuous field; a county value averages a 15 km grid, weighted by population.
• Same day, three “national averages”: say which weight you used.

Healy calls it the fallacy of misplaced concreteness: treating the unit where a number is stored (a county) as if it were the thing itself (a smoke field). The people-weighted average is high because smoke sat over the crowded Northeast; the area-weighted average is low because big, empty western counties had little smoke. st_area() works on longitude–latitude data because sf measures on the sphere; the result is in m², so dividing by 1e6 gives km². Another fix for unequal unit sizes is a statebins or hex cartogram, which gives every state the same size (see the small-multiples cautions).

✅ Wrap-up

✅ Checklist for an honest map: before you draw

• Projection fits the area: Albers for the lower 48, Lambert for North America.
• Keys checked: GEOID as text, boundary year matches data year, anti_join() returns zero rows.
• Rate or count? County totals mostly map population.
• Missing stays missing (grey); never code it as zero.

✅ Checklist for an honest map: on the page

• Scale and bins are stated, with reasons (35 µg/m³ = EPA daily standard for total PM2.5).
• Time convention is stated: HMS date is the UTC date; FIRMS local_date is the local solar date; smoke, EMBER, and AQS are calendar days.
• Units, date, geography, and source are in the caption.
• Models and plumes are labelled as what they are: modeled, or overhead.

Every classwork figure should carry labs(title, subtitle, caption) that states the date, unit, geography, and source.

📚 Data sources and licences

Data	Source	Licence
County and state boundaries, population	U.S. Census Bureau, 2020 cartographic boundaries and centers of population	Public domain
Province boundaries	Statistics Canada, 2021 boundary files	Statistics Canada Open Licence
Smoke PM2.5 (daily, monthly)	CSU kriged smoke PM2.5, Pierce et al., Dryad DOI 10.5061/dryad.k0p2ngfhv	CC0
Modeled fire PM2.5	EPA EMBER CMAQ runs, Zenodo 13737754	CC BY 4.0
Monitor PM2.5	EPA Air Quality System, daily PM2.5	Public domain
Smoke plumes	NOAA Hazard Mapping System	Public domain
Fire detections	NASA FIRMS MODIS C6.1, doi:10.5067/FIRMS/MODIS/MCD14ML	Open; cite NASA FIRMS
Burned-area perimeters	Canadian Forest Service, CWFIS National Burned Area Composite	Open Government Licence – Canada

Methods: Healy, Data Visualization, ch. 7 “Draw Maps” (socviz.co); DANL 310 Lecture 8. Build script: wildfire-maps-lecture/R/build_wildfire_map_data.R.

📌 Main takeaways

• sf keeps geometry in one column, so dplyr verbs and joins just work.
• Draw in a projection with coord_sf(); st_transform() changes the data, and layers must share a CRS before a join.
• A choropleth’s message depends on keys, rates, bins, and missing values.
• Layer points, plumes, and polygons, then ask whether a map is the best chart.

If you could add one more layer to the June 7 map to show who breathed the smoke, what would it be, and how could it still mislead?

Possible answers: monitors (sparse), population density (ur-choropleth risk), hourly smoke (timing), indoor/outdoor time (not in our data). Next: Classwork 5 maps the Midwest episode of June 27–29, 2023 with the same tools.