Classwork 5

Map Canadian Wildfire Smoke with geom_sf()

Read boundary and smoke data from course URLs, join them, and build projected maps with geom_sf(): choropleths, points, plumes, small multiples, and spatial joins.
Published

September 30, 2026

ImportantGoal

Map the Canadian wildfire smoke episode that reached the U.S. Midwest on June 27–29, 2023, and connect it to the 2023 Canadian fire season. Every map should state its date, unit, geography, and data source.

Work with a partner. One person drives the code while the other checks join keys, projections, color scales, and labels. Switch roles after Part 2.

NoteCore path for one 75-minute meeting

This classwork is longer than one class meeting. In class, complete the core path: Parts 0–4 and Part 7A. Parts 5, 6, 7B, 8, 9, and 10 are extensions to finish after class or to use as practice for your project.

By the end, you should be able to:

  1. read boundary files with read_sf() and data tables with read_csv() directly from their URLs, keeping GEOID as text,
  2. draw and layer maps with geom_sf() and choose a projection with coord_sf(),
  3. join a data table to county polygons by GEOID, build continuous and binned choropleth maps, and keep missing values missing,
  4. turn longitude–latitude tables into points with st_as_sf() and compare days with facet_wrap() small multiples,
  5. combine layers with spatial joins (st_filter(), st_join()) and combine counties with group_by() and summarize(), and
  6. explain in plain words what each map can and cannot show.

🧰 Part 0: Read the course files

The course publishes every file used in this classwork at https://bcdanl.github.io/data/.

File One row is
wildfire_us_counties_2020.geojson a county of the contiguous U.S. (incl. DC), with 2020 population
wildfire_us_states_2020.geojson a state (48 states + DC)
wildfire_canada_provinces_2021.geojson a Canadian province or territory
wildfire_hms_smoke_2023_06.geojson a smoke plume drawn by NOAA analysts from satellite images, June 2023
wildfire_nbac_burned_area_2023.geojson the burned-area outline of one 2023 Canadian fire
wildfire_smoke_county_day_2023_06.csv a county on one day of June 2023 (smoke and total PM2.5)
wildfire_smoke_county_month_2006_2025.csv a county in one month, 2006–2025
wildfire_ember_county_day_2023_06.csv a county on one day of June 2023 (EPA model estimates)
wildfire_aqs_pm25_2023_06.csv an air-quality monitor on one day of June 2023
wildfire_firms_canada_2023.csv a satellite fire detection in Canada, 2023

PM2.5 is fine particulate matter, measured in micrograms per cubic meter (µg/m³). The files describe fires, smoke, air quality, boundaries, and population.

In the smoke tables, total_pm25 is PM2.5 from all sources. smoke_pm25 is the part of PM2.5 above a place’s usual no-smoke level, counted only on days when satellite images show a smoke plume overhead; on other days it is set to 0. Both are estimated on a 15 km grid that fills in the gaps between air monitors (researchers call this method kriging, hence “kriged” in the map captions), then averaged over each county’s residents (weighted by population). Because smoke_pm25 is a difference, it can be slightly negative.

0A. Load packages

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

0B. Read the files from their URLs

Use read_sf() for the five GeoJSON boundary and shape files and read_csv() for the five tables. The county code GEOID is a five-character FIPS code with leading zeros ("01001" is Autauga County, Alabama), so tell read_csv() to keep it as text with col_character(). Do the same for the monitor ID site_id and the fire detection time acq_time ("0450" means 04:50 UTC).

# 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())
)

0C. Check rows, keys, and the coordinate reference system

c(
  counties = nrow(___),
  states = nrow(___),
  provinces = nrow(___),
  smoke_counties = n_distinct(smoke_day$___)
)

st_crs(counties)$epsg
class(counties)

as.numeric("01001")
  1. Do the row counts match 3,108 counties, 49 states (48 + DC), and 13 provinces and territories? Does the daily smoke table cover every county?
  2. Which coordinate reference system (CRS) do the files use, and in what units are the coordinates?
  3. What does as.numeric("01001") return? Why would a numeric GEOID break a join with counties?

0D. Set a map theme once

theme_void() removes axes and grid lines that a map does not need. Setting it once with theme_set() applies it to every later plot.

theme_set(
  theme_void(base_size = 12) +
    theme(
      plot.title = element_text(face = "bold"),
      plot.caption = element_text(hjust = 0, color = "grey40"),
      plot.margin = margin(8, 8, 8, 8)
    )
)

🗺️ Part 1: First maps and layers

1A. Counties with state borders on top

Draw the counties first, then the states with no fill so that only their borders show. coord_sf(crs = "ESRI:102003") draws the map in an Albers equal-area projection for the contiguous U.S.; Part 2 compares projections.

ggplot() +
  geom_sf(data = ___, fill = "grey92", color = "white", linewidth = 0.05) +
  geom_sf(data = ___, fill = ___, color = "grey25", linewidth = 0.3) +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "Counties and states of the contiguous United States",
    subtitle = "3,108 counties (including DC) under 48 state borders and DC",
    caption = "Source: U.S. Census Bureau 2020 cartographic boundaries. Albers equal-area projection."
  )

What happens if you swap the order of the two geom_sf() layers? Why do the states need fill = NA?

1B. Dissolve counties into Census divisions

group_by() followed by summarize() on an sf object also merges (unions) the polygons in each group. Combine the counties into the nine Census divisions, adding up the population and counting the counties.

divisions <- counties |>
  group_by(region, ___) |>
  summarize(
    pop_2020 = sum(___),
    n_counties = ___,
    .groups = "drop"
  )

divisions |>
  st_drop_geometry() |>
  arrange(desc(pop_2020))

st_drop_geometry() removes the geometry column, so the result prints as an ordinary table; use it whenever you want numbers rather than a map.

Now map the divisions and label each one with its name and population. The course county file was simplified one shared border at a time, so two neighboring counties still share exactly the same border line, and the dissolve leaves no gaps inside a division. A thin white outline (color = "white") separates neighboring divisions cleanly.

ggplot(divisions) +
  geom_sf(aes(fill = ___), color = "white", linewidth = 0.3, show.legend = FALSE) +
  geom_sf_label(
    aes(label = paste0(str_wrap(division, 12), "\n", round(pop_2020 / 1e6, 1), "M")),
    size = 2.6,
    lineheight = 0.9
  ) +
  scale_fill_brewer(palette = "Set3") +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "The nine Census divisions, dissolved from 3,108 counties",
    subtitle = "Labels show each division's 2020 population in millions",
    caption = "Source: U.S. Census Bureau 2020 cartographic boundaries and population. Albers equal-area projection."
  )

Which division has the most people, and how many counties does it contain? Which division has the most counties? Do the dissolved polygons still show county borders? Which lines does the white outline draw?

🧭 Part 2: Projections

2A. The same map in two projections

Save the base map of Canadian provinces and U.S. states as an object, then add two different coord_sf() projections: 4326 (unprojected longitude/latitude) and "ESRI:102009" (North America Lambert Conformal Conic).

north_america <- ggplot() +
  geom_sf(data = provinces, fill = "grey85", color = "white", linewidth = 0.2) +
  geom_sf(data = states, fill = "grey70", color = "white", linewidth = 0.2) +
  labs(
    subtitle = "13 Canadian provinces and territories; 48 U.S. states and DC",
    caption = "Sources: Statistics Canada 2021 and U.S. Census Bureau 2020 boundaries."
  )

north_america +
  coord_sf(crs = ___) +
  labs(title = "Unprojected longitude and latitude (EPSG:4326)")

north_america +
  coord_sf(crs = ___) +
  labs(title = "North America Lambert Conformal Conic (ESRI:102009)")

Compare the Arctic islands and the U.S.–Canada border along the 49th parallel in the two maps. Which projection makes northern Canada look too large?

2B. coord_sf() versus st_transform()

coord_sf(crs = ...) changes only how a plot is drawn. st_transform() changes the coordinates stored in the data. Transform the states to ESRI:102009 and compare the bounding boxes.

states_lcc <- st_transform(states, crs = ___)

st_bbox(states)
st_bbox(states_lcc)
st_crs(states_lcc)$units

In what units are the two bounding boxes? When would you want to transform the data rather than only the plot?

2C. Zoom in on the Great Lakes

Zoom in to longitudes −97 to −74 and latitudes 36 to 50. Because these limits are in degrees while the map is projected, tell coord_sf() that the limits use EPSG:4326 with default_crs = st_crs(4326).

north_america +
  coord_sf(
    crs = "ESRI:102009",
    xlim = c(___, ___),
    ylim = c(___, ___),
    default_crs = ___
  ) +
  labs(title = "The Great Lakes and the Midwest")

What happens if you delete the default_crs line? Why?

🔗 Part 3: Join smoke data to counties — June 28, 2023

3A. Filter one day and check the keys

Before any join, check that the keys match in both directions with anti_join(), which keeps the rows that have no match.

smoke_0628 <- smoke_day |>
  filter(date == as.Date(___))

nrow(smoke_0628)

smoke_0628 |> anti_join(counties, by = ___) |> nrow()
counties |> anti_join(smoke_0628, by = ___) |> nrow()

How many rows does smoke_0628 have? How many county codes fail to match in each direction?

3B. Join the table to the polygons

smoke_0628_map <- ___ |>
  left_join(___, by = "GEOID")

class(smoke_0628_map)
nrow(smoke_0628_map)

Why should the sf object counties be on the left side of left_join()?

3C. A continuous color scale

Map the smoke PM2.5 column with a continuous color scale.

ggplot(smoke_0628_map) +
  geom_sf(aes(fill = ___), color = NA) +
  geom_sf(data = states, fill = NA, color = "grey30", linewidth = 0.2) +
  scale_fill_viridis_c(option = "inferno", direction = -1) +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "Wildfire smoke PM2.5 on June 28, 2023",
    subtitle = "Daily smoke PM2.5, county population-weighted average, contiguous U.S.",
    fill = "µg/m³",
    caption = "Source: CSU kriged smoke PM2.5 (O'Dell et al.). Albers equal-area projection."
  )

Where is the smoke heaviest? Can you tell a county at 10 µg/m³ from one at 30 µg/m³?

3D. A binned color scale with meaningful breaks

Use breaks at 5, 15, 35, 75, and 150 µg/m³; 35 µg/m³ is the EPA 24-hour PM2.5 standard. Small negative smoke values (measurement noise from subtracting the usual level; see Part 0) fall into the lowest bin.

ggplot(smoke_0628_map) +
  geom_sf(aes(fill = smoke_pm25), color = NA) +
  geom_sf(data = states, fill = NA, color = "grey30", linewidth = 0.2) +
  scale_fill_fermenter(
    palette = "YlOrRd",
    direction = 1,
    breaks = ___
  ) +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "Wildfire smoke PM2.5 on June 28, 2023",
    subtitle = "Daily smoke PM2.5, county population-weighted average; 35 = EPA 24-hour standard",
    fill = "µg/m³",
    caption = "Source: CSU kriged smoke PM2.5 (O'Dell et al.). Albers equal-area projection."
  )

smoke_0628_map |>
  st_drop_geometry() |>
  summarize(
    counties_over_35 = sum(smoke_pm25 > ___),
    people_over_35_millions = sum(pop_2020[smoke_pm25 > ___]) / 1e6
  )

Which map shows more clearly where smoke exceeded 35 µg/m³? How many counties, and how many people, had smoke above 35 µg/m³ on June 28?

Optional: replace the breaks with quantile(smoke_0628_map$smoke_pm25, c(0.2, 0.4, 0.6, 0.8)), which puts one-fifth of the counties in each color. Which counties change color? Which set of breaks would you show a reader, and why?

3E. The five smokiest counties

smoke_0628_map |>
  st_drop_geometry() |>
  slice_max(___, n = ___) |>
  select(county_name, STUSPS, smoke_pm25, total_pm25, pop_2020)

Which state holds the five smokiest counties? Are they large or small counties by population?

3F. Missing is not zero

Air-quality monitors exist in only some counties. When you join monitor readings to all counties with left_join(), a county without a monitor gets NA. na.value sets the color for NA values; a binned legend does not show it, so say in the subtitle what grey means.

monitor_0628 <- aqs |>
  filter(date == as.Date("2023-06-28")) |>
  group_by(___) |>
  summarize(pm25 = mean(pm25))

monitor_map <- counties |>
  left_join(monitor_0628, by = ___)

ggplot(monitor_map) +
  geom_sf(aes(fill = pm25), color = NA) +
  geom_sf(data = states, fill = NA, color = "grey30", linewidth = 0.2) +
  scale_fill_fermenter(
    palette = "YlOrRd", direction = 1,
    breaks = c(5, 15, 35, 75, 150),
    na.value = ___
  ) +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "Only some counties have an air-quality monitor, June 28, 2023",
    subtitle = "Daily total PM2.5, mean of the EPA monitors in each county (µg/m³); grey = no monitor reported",
    fill = "µg/m³",
    caption = "Source: EPA Air Quality System. Albers equal-area projection."
  )

monitor_map |>
  st_drop_geometry() |>
  summarize(
    counties = n(),
    with_monitor = sum(!is.na(pm25)),
    share_of_people = sum(pop_2020[!is.na(pm25)]) / sum(pop_2020),
    mean_na_left_out = mean(pm25, na.rm = ___),
    mean_na_as_zero = mean(replace_na(pm25, ___))
  )

How many counties have a monitor, and what share of people live in them? Why is a grey county not a county with clean air? What happens to the average if you code grey counties as 0, and why is that wrong?

⚖️ Part 4: Counts versus rates

The monthly table has days_gt15, the number of days in a month with smoke PM2.5 above 15 µg/m³. Multiplying by population gives person-days: the total number of smoky days lived by all residents of a county. One is a rate per resident; the other is a total count.

june_2023 <- smoke_month |>
  filter(year == ___, month == ___)

june_map <- counties |>
  left_join(june_2023, by = "GEOID") |>
  mutate(person_days = ___ * ___)

Map the rate (days per resident) and the count (person-days, in millions).

ggplot(june_map) +
  geom_sf(aes(fill = ___), color = NA) +
  geom_sf(data = states, fill = NA, color = "grey30", linewidth = 0.2) +
  scale_fill_fermenter(palette = "YlOrRd", direction = 1, breaks = c(1, 4, 8, 12, 16)) +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "Rate: smoky days in June 2023",
    subtitle = "Days with smoke PM2.5 above 15 µg/m³, by county, contiguous U.S.",
    fill = "Days",
    caption = "Source: CSU kriged smoke PM2.5 (O'Dell et al.). Albers equal-area projection."
  )

ggplot(june_map) +
  geom_sf(aes(fill = ___ / 1e6), color = NA) +
  geom_sf(data = states, fill = NA, color = "grey30", linewidth = 0.2) +
  scale_fill_fermenter(palette = "YlOrRd", direction = 1, breaks = c(0.1, 0.5, 2, 5, 20)) +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "Count: smoky person-days in June 2023",
    subtitle = "Days with smoke PM2.5 above 15 µg/m³ × 2020 population, by county",
    fill = "Million\nperson-days",
    caption = "Sources: CSU kriged smoke PM2.5 (O'Dell et al.); U.S. Census Bureau 2020 population."
  )

june_map |>
  st_drop_geometry() |>
  slice_max(___, n = 5) |>
  select(county_name, STUSPS, days_gt15, pop_2020, person_days)

june_map |>
  st_drop_geometry() |>
  slice_max(days_gt15, n = 3) |>
  select(county_name, STUSPS, days_gt15, pop_2020, person_days)

Which map looks more like a map of where people live? Which county had the most person-days, and did it also have the most smoky days? Which map would you use to compare how smoky June was for a typical resident of each county?

📍 Part 5: Air-quality monitors as points

5A. Monitors on June 28 over the county map

A monitor table has longitude and latitude columns. st_as_sf() turns them into point geometry; list longitude (x) before latitude (y), and set the CRS to 4326.

aqs_0628 <- aqs |>
  filter(date == as.Date("2023-06-28")) |>
  st_as_sf(coords = c(___, ___), crs = ___)

nrow(aqs_0628)

Monitors measure PM2.5 from all sources, so compare them with the county total_pm25, not smoke_pm25. Both layers map a variable to fill (the points use shape = 21, a circle with a fill), so one binned scale colors the counties and the monitors the same way.

ggplot() +
  geom_sf(data = smoke_0628_map, aes(fill = ___), color = NA) +
  geom_sf(data = states, fill = NA, color = "grey30", linewidth = 0.2) +
  geom_sf(
    data = aqs_0628, aes(fill = ___),
    shape = 21, color = "black", size = 1.8, stroke = 0.3
  ) +
  scale_fill_fermenter(palette = "YlOrRd", direction = 1, breaks = c(5, 15, 35, 75, 150)) +
  coord_sf(
    crs = "ESRI:102003",
    xlim = c(-97, -70), ylim = c(35, 48),
    default_crs = st_crs(4326)
  ) +
  labs(
    title = "County estimates and EPA monitors, June 28, 2023",
    subtitle = "Daily total PM2.5: counties (population-weighted) and EPA monitors (circles)",
    fill = "µg/m³",
    caption = "Sources: CSU kriged PM2.5 (O'Dell et al.); EPA Air Quality System. Albers equal-area projection."
  )

How many monitors reported on June 28? Do most circles match the color of the county around them?

5B. Monitors versus county values

Join each monitor to its county’s total_pm25 by GEOID. Put the county value on the x-axis and the monitor value on the y-axis. The dashed 45° line marks perfect agreement. cor() computes the correlation, a number from −1 to 1 that measures how closely the points follow one straight line.

aqs_compare <- aqs |>
  filter(date == as.Date("2023-06-28")) |>
  left_join(
    smoke_0628 |> select(GEOID, total_pm25),
    by = ___
  )

ggplot(aqs_compare, aes(x = ___, y = ___)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey40") +
  geom_point(alpha = 0.5) +
  labs(
    title = "Monitors versus their county's estimate, June 28, 2023",
    subtitle = "Daily total PM2.5 (µg/m³); dashed line = perfect agreement",
    x = "County estimate (population-weighted, µg/m³)",
    y = "Monitor reading (µg/m³)",
    caption = "Sources: EPA Air Quality System; CSU kriged PM2.5 (O'Dell et al.)."
  ) +
  theme_minimal(base_size = 12)

aqs_compare |>
  summarize(
    monitors = n(),
    correlation = cor(pm25, total_pm25),
    share_within_10 = mean(abs(pm25 - total_pm25) <= 10)
  )

aqs_compare |>
  mutate(gap = pm25 - total_pm25) |>
  slice_max(abs(gap), n = 3)

How closely do monitors and county values agree? Look at the largest gaps. Give two reasons why a single monitor can differ from its county’s value.

🔥 Part 6: Fires and burned areas in Canada

6A. Fire detections by month

Most 2023 detections occurred from May through September (months 5 to 9). local_date is the local solar date of each detection (the UTC time shifted by longitude), so a nighttime detection counts on the local day it happened. Keep those months, label them, and turn the table into points. The fire table is sorted by date, so geom_sf() would draw September last, on top of earlier months in the same places. Shuffle the rows with slice_sample(prop = 1) so that no month is always drawn on top.

set.seed(399)
fires_sf <- fires |>
  filter(month(local_date) %in% ___) |>
  mutate(month = month(local_date, label = TRUE)) |>
  slice_sample(prop = 1) |>  # shuffle the drawing order
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

nrow(fires_sf)

ggplot() +
  geom_sf(data = provinces, fill = "grey95", color = "grey60", linewidth = 0.2) +
  geom_sf(data = fires_sf, aes(color = ___), size = 0.05, alpha = 0.4) +
  scale_color_viridis_d(option = "plasma", end = 0.9) +
  guides(color = guide_legend(override.aes = list(size = 3, alpha = 1))) +
  coord_sf(crs = "ESRI:102009") +
  labs(
    title = "Where and when Canada burned in 2023",
    subtitle = "MODIS active-fire detections, May–September 2023, colored by local month",
    color = "Month",
    caption = "Sources: NASA FIRMS MODIS Collection 6.1; Statistics Canada 2021 boundaries. Lambert conformal conic projection."
  )

Count the detections by province and month. pivot_wider() spreads the months into columns, so each province becomes one row.

fires_sf |>
  st_drop_geometry() |>
  count(prov_abbr, ___) |>
  pivot_wider(names_from = month, values_from = n, values_fill = 0) |>
  arrange(desc(Jun))

How many detections are on the map? Using the map and the table, which parts of Canada burned early in the season (May–June), and which burned later (July–September)? Why can a single map with many overlapping points make one month look bigger than it was?

6B. Check a point-in-polygon join

The fire table already has a province label, prov_abbr. Use st_join() to find which province polygon each point falls in, and compare the two labels.

fires_joined <- fires_sf |>
  ___(provinces |> select(prov_join = prov_abbr))

nrow(fires_sf)
nrow(fires_joined)

fires_joined |>
  st_drop_geometry() |>
  count(match = prov_abbr == ___)

fires_joined |>
  st_drop_geometry() |>
  filter(is.na(prov_join) | prov_abbr != prov_join) |>
  count(prov_abbr, prov_join, sort = TRUE)

Does fires_joined have the same number of rows as fires_sf? (st_join() returns one row for every polygon a point falls in, so polygons that overlap would add rows.) How many detections got a different province, or no province at all (NA)? Where are they, and why can a spatial join disagree with the label?

6C. Burned-area outlines by cause

Each outline is the fire’s outer edge with the unburned patches inside it filled in, so it can look larger than the area that actually burned. area_ha is the official burned area in hectares (100 ha = 1 km²).

nbac |>
  st_drop_geometry() |>
  group_by(___) |>
  summarize(fires = n(), burned_million_ha = sum(area_ha) / 1e6)

ggplot() +
  geom_sf(data = provinces, fill = "grey95", color = "grey60", linewidth = 0.2) +
  geom_sf(data = nbac, aes(fill = ___), color = NA) +
  scale_fill_manual(
    values = c(Natural = "#d95f02", Human = "#1b9e77", Undetermined = "grey45")
  ) +
  coord_sf(crs = "ESRI:102009") +
  labs(
    title = "Burned areas of 2023 Canadian fires, by cause",
    subtitle = "Outer outlines of 968 fires with parts of at least 1 km²",
    fill = "Cause",
    caption = "Source: CWFIS National Burned Area Composite (Natural Resources Canada). Lambert conformal conic projection."
  )

Which cause accounts for most of the burned area? Which cause produced the most fires in this file?

🎞️ Part 7: Small multiples

7A. Daily smoke, June 26–30, 2023

Keep five days, join them to the counties, and draw one small map per day with facet_wrap(). Zoom to the Midwest and East.

episode <- smoke_day |>
  filter(date >= as.Date("2023-06-26"), date <= as.Date(___))

episode_map <- counties |>
  left_join(episode, by = "GEOID")

nrow(episode_map)

ggplot(episode_map) +
  geom_sf(aes(fill = smoke_pm25), color = NA) +
  geom_sf(data = states, fill = NA, color = "grey30", linewidth = 0.15) +
  facet_wrap(~ ___, nrow = 2) +
  scale_fill_fermenter(palette = "YlOrRd", direction = 1, breaks = c(5, 15, 35, 75, 150)) +
  coord_sf(
    crs = "ESRI:102003",
    xlim = c(-100, -67), ylim = c(33, 49),
    default_crs = st_crs(4326)
  ) +
  labs(
    title = "Wildfire smoke PM2.5, June 26–30, 2023",
    subtitle = "Daily smoke PM2.5, county population-weighted average, Midwest and East",
    fill = "µg/m³",
    caption = "Source: CSU kriged smoke PM2.5 (O'Dell et al.). Albers equal-area projection."
  )

How many rows does episode_map have, and why? Where was the heaviest smoke on each day, and which way did it move from June 26 to June 30?

7B. Annual smoky days, 2019–2025

Add up the monthly days_gt15 into a yearly total for each county from 2019 onward, then draw one map per year. Finally, compute the share of counties in each Census region with 10 or more smoky days in each year.

annual_days <- smoke_month |>
  filter(year >= ___) |>
  group_by(GEOID, ___) |>
  summarize(days_gt15 = sum(days_gt15), .groups = "drop")

annual_map <- counties |>
  left_join(annual_days, by = "GEOID")

ggplot(annual_map) +
  geom_sf(aes(fill = days_gt15), color = NA) +
  facet_wrap(~ ___, ncol = 4) +
  scale_fill_fermenter(palette = "YlOrRd", direction = 1, breaks = c(1, 5, 10, 20, 40)) +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "Smoky days per year, 2019–2025",
    subtitle = "Days with smoke PM2.5 above 15 µg/m³, by county, contiguous U.S.",
    fill = "Days",
    caption = "Source: CSU kriged smoke PM2.5 (O'Dell et al.). Albers equal-area projection."
  )
annual_map |>
  st_drop_geometry() |>
  group_by(region, year) |>
  summarize(share_10_plus = mean(days_gt15 >= ___), .groups = "drop") |>
  pivot_wider(names_from = region, values_from = share_10_plus)

Which year stands out in the Midwest and East? Which year stands out in the West?

☁️ Part 8: Smoke overhead versus smoke at breathing height

NOAA analysts draw smoke plumes from satellite images. A plume shows smoke somewhere in the air above a place, often high up; it does not show whether that smoke reached the ground. The plume date is the UTC date of the plume’s start time.

8A. Counties under Medium or Heavy plumes on June 28

st_filter() keeps the rows of one sf object that overlap (intersect) any shape in another: a county is kept if any part of it lies inside a plume, whether it sits fully inside or crosses the plume’s edge.

plumes_0628 <- plumes |>
  filter(
    date == as.Date("2023-06-28"),
    density %in% c(___, ___)
  ) |>
  mutate(density = factor(density, levels = c("Medium", "Heavy")))

counties_under <- counties |>
  ___(plumes_0628)

nrow(counties_under)

ggplot() +
  geom_sf(data = counties, fill = "grey95", color = NA) +
  geom_sf(data = counties_under, fill = "#fdae6b", color = NA) +
  geom_sf(data = states, fill = NA, color = "grey40", linewidth = 0.15) +
  geom_sf(data = plumes_0628, aes(color = density), fill = NA, linewidth = 0.4) +
  scale_color_manual(values = c(Medium = "#3182bd", Heavy = "black")) +
  coord_sf(
    crs = "ESRI:102003",
    xlim = c(-125, -66.5), ylim = c(24, 50),
    default_crs = st_crs(4326)
  ) +
  labs(
    title = "Counties under a Medium or Heavy smoke plume, June 28, 2023 (UTC)",
    subtitle = "Orange = county overlaps a plume (inside it or crossing its edge); lines = plume outlines",
    color = "Plume density",
    caption = "Sources: NOAA Hazard Mapping System smoke polygons; U.S. Census Bureau 2020 boundaries."
  ) +
  theme(plot.subtitle = element_text(margin = margin(b = 8)))

About how much of the country sits under the plumes on June 28? How many counties overlap a plume?

8B. Population under the plumes versus population with smoky air

plume_check <- smoke_0628_map |>
  st_drop_geometry() |>
  mutate(
    under_plume = GEOID %in% counties_under$___,
    over_35 = smoke_pm25 > ___
  )

plume_check |>
  summarize(
    people_under_plume_millions = sum(pop_2020[under_plume]) / 1e6,
    people_over_35_millions = sum(pop_2020[over_35]) / 1e6
  )

plume_check |>
  group_by(under_plume, over_35) |>
  summarize(
    counties = n(),
    people_millions = sum(___) / 1e6,
    .groups = "drop"
  )

How many people lived under a Medium or Heavy plume, and how many lived in counties with smoke PM2.5 above 35 µg/m³? Why are the two numbers so different?

🍁 Part 9: How much fire smoke came from Canada?

The EPA EMBER project ran an air-quality computer model three ways: with all emissions, with no fires, and with no Canadian fires. The differences give modeled PM2.5 from all fires (pm25_all_fires) and from Canadian fires (pm25_canada_fires). These are model estimates, not measurements.

9A. County shares for June 2023

Average each county’s daily values over June, then divide Canadian-fire PM2.5 by all-fire PM2.5.

ember_county <- ember |>
  group_by(___) |>
  summarize(
    all_fires = mean(pm25_all_fires),
    canada_fires = mean(pm25_canada_fires)
  ) |>
  mutate(canada_share = ___ / ___)

summary(ember_county$canada_share)

A share runs from 0 (no fire smoke from Canada) to 1 (all of it from Canada), so use a sequential color scale with fixed limits. Save diverging scales for differences, such as June 2023 minus a usual June, where 0 is the natural middle.

counties |>
  left_join(ember_county, by = "GEOID") |>
  ggplot() +
  geom_sf(aes(fill = canada_share), color = NA) +
  geom_sf(data = states, fill = NA, color = "grey30", linewidth = 0.15) +
  scale_fill_gradient(
    low = "grey97", high = "#b2182b",
    limits = c(___, ___), labels = label_percent()
  ) +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "Share of fire PM2.5 from Canadian fires, June 2023",
    subtitle = "Modeled June mean, by county; darker red = larger share from Canada",
    fill = "Share",
    caption = "Source: EPA EMBER CMAQ model runs (modeled, not measured). Albers equal-area projection."
  )

What is the median county share? Where did Canadian fires supply less than half of the fire smoke?

9B. Aggregate counties to states with population weights

A state value should describe the typical resident, so weight each county by its population with weighted.mean(). Then join the state table to the state polygons.

ember_state <- ember_county |>
  left_join(st_drop_geometry(counties), by = "GEOID") |>
  group_by(STUSPS) |>
  summarize(
    all_fires = weighted.mean(all_fires, w = ___),
    canada_fires = weighted.mean(canada_fires, w = ___)
  ) |>
  mutate(canada_share = canada_fires / all_fires)

state_map <- states |>
  left_join(ember_state, by = ___)

ggplot(state_map) +
  geom_sf(aes(fill = canada_share), color = "white", linewidth = 0.3) +
  geom_sf_text(
    data = filter(state_map, STUSPS != "DC"), aes(label = STUSPS),
    size = 2.3, color = "grey20", check_overlap = TRUE
  ) +
  scale_fill_gradient(
    low = "grey97", high = "#b2182b",
    limits = c(0, 1), labels = label_percent()
  ) +
  coord_sf(crs = "ESRI:102003") +
  labs(
    title = "Share of fire PM2.5 from Canadian fires, by state, June 2023",
    subtitle = "Population-weighted state means of modeled county values",
    fill = "Share",
    caption = "Source: EPA EMBER CMAQ model runs (modeled, not measured); U.S. Census Bureau 2020 population."
  )

ember_state |>
  filter(canada_share < 0.5) |>
  arrange(canada_share)

ember_state |> filter(STUSPS %in% c("NM", "IL"))

Which states received less than half of their fire smoke from Canadian fires? Why weight by population instead of taking a simple average over counties? New Mexico has the lowest share. Compare its all_fires with Illinois’s: why does a low share in New Mexico say little about how smoky its air was?

📈 Part 10: Is the data really spatial?

A map shows one moment at a time. To compare when smoke peaked in different regions, a line chart can work better. Compute population-weighted daily smoke for each Census region.

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

ggplot(regional_daily, aes(x = date, y = smoke_pm25, color = region)) +
  geom_line(linewidth = 1) +
  scale_color_viridis_d(option = "D", end = 0.85) +
  scale_x_date(date_breaks = "1 week", date_labels = "%b %d") +
  labs(
    title = "Two smoke episodes in June 2023, one region at a time",
    subtitle = "Daily smoke PM2.5, population-weighted mean by Census region",
    x = NULL,
    y = "Smoke PM2.5 (µg/m³)",
    color = "Region",
    caption = "Sources: CSU kriged smoke PM2.5 (O'Dell et al.); U.S. Census Bureau 2020 population."
  ) +
  theme_minimal(base_size = 12)

regional_daily |>
  group_by(region) |>
  slice_max(smoke_pm25, n = 1)

On which day did each region peak? What does this chart show more clearly than the maps in Part 7?

✅ Exit ticket

Submit:

  1. Projection. Which CRS would you use for one map that shows Canadian fire detections and U.S. county smoke together, and why?
  2. Rate or count. To show where June 2023 was smokiest for a typical resident, would you map days_gt15 or person_days? Explain in one sentence.
  3. Limits of a plume map. In one sentence, state what an HMS plume map cannot show.
  4. One revised map from this classwork whose title, subtitle, and caption state the date, unit, geography, and source.

Sources

Back to top