Introduction

The story these maps tell is that tornadoes in Oklahoma between 2016 and 2021 were not spread evenly across space or across time. Touchdowns and paths are concentrated in the central and eastern parts of the state, the number of tornadoes varies a lot from year to year (2019 was far more active than any other year), and the picture of “where tornadoes are most common” depends on how the data are summarized and classified. I build this story in three steps, following the practice exercises in Wimberly’s (2023) chapter on vector geospatial data: (1) maps of tornado paths and touchdown points by year, (2) county-level tornado density for each year, and (3) a comparison of density maps that use 3 to 6 quantile classes.

Setup

The sf package stores spatial data as data frames with a geometry column, so the usual dplyr and tidyr functions work on spatial data. ggplot2 draws the maps, units handles measurement units (areas), and cowplot combines several maps into one figure.

# Spatial data
library(sf)

# Data wrangling and plotting
library(dplyr)
library(tidyr)
library(ggplot2)

# Area units and multi-panel figures
library(units)
library(cowplot)

I read the three shapefiles with st_read(): the Oklahoma county boundaries (polygons), tornado touchdown points, and tornado paths (lines). The tornado files come from the NOAA Storm Prediction Center.

# Folder that holds the unzipped Chapter 5 shapefiles
data_dir <- "C:/Users/Wiest/Documents/chapter5_data/Chapter5/"

okcounty <- st_read(paste0(data_dir, "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(paste0(data_dir, "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(paste0(data_dir, "ok_tornado_path.shp"), quiet = TRUE)

The tornado files cover 1950-2021. The study period here is 2016-2021, so I filter both layers by year and keep only the columns needed for mapping.

# Touchdown points, 2016-2021
tpoint_16_21 <- tpoint |>
  filter(yr >= 2016, yr <= 2021) |>
  select(om, yr, date)

# Tornado paths, 2016-2021
tpath_16_21 <- tpath |>
  filter(yr >= 2016, yr <= 2021) |>
  select(om, yr, date)

# Number of tornadoes per year
table(tpoint_16_21$yr)
## 
## 2016 2017 2018 2019 2020 2021 
##   57   85   43  148   38   63

Exercise 1: Tornado paths and points by year

To show the year of each tornado, I map yr to the color aesthetic. Wrapping it in as.factor() makes ggplot use a discrete color scale (one color per year) instead of a continuous gradient. County boundaries are added last so they are drawn on top, and theme_void() removes axes so only the map is shown.

# Map of touchdown points, colored by year
point_map <- ggplot() +
  geom_sf(data = tpoint_16_21, aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

# Map of tornado paths, colored by year
path_map <- ggplot() +
  geom_sf(
    data = tpath_16_21,
    aes(color = as.factor(yr)),
    linewidth = 0.8
  ) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

plot_grid() from cowplot stacks the two maps into a single figure. The hjust and label_x arguments left-align the panel labels so they do not cover the map, and align = "hv" makes the panels the same size.

plot_grid(
  path_map, point_map,
  labels = c("A) Tornado Paths", "B) Tornado Points"),
  label_size = 12,
  ncol = 1,
  hjust = 0,
  label_x = 0,
  align = "hv"
)

Exercise 2: County-level tornado density by year

This exercise needs one row for every county and year. First, a spatial join (st_join()) attaches the county of each tornado touchdown point. Then I count tornadoes for each county-year combination.

# Attach county attributes to each tornado point (join by location)
countypnt <- st_join(tpoint_16_21, okcounty)

# Count tornadoes per county and year (drop geometry for a plain summary)
countysum_year <- countypnt |>
  st_drop_geometry() |>
  group_by(GEOID, yr) |>
  summarize(tcnt = n(), .groups = "drop")

A county-year with no tornadoes does not appear in that summary. To avoid missing rows, expand_grid() builds every county-year combination, the counts are joined on, and missing counts are set to zero.

county_year <- expand_grid(GEOID = okcounty$GEOID, yr = 2016:2021) |>
  # Keep all combinations; counties without tornadoes get NA counts
  left_join(countysum_year, by = c("GEOID", "yr")) |>
  # NA here means "no tornadoes", so replace with zero
  replace_na(list(tcnt = 0))

Next I join the counts to the county polygons and compute density. st_area() returns square meters, so the area is converted to square kilometers and density is expressed as tornadoes per 1000 km².

county_density <- okcounty |>
  left_join(county_year, by = "GEOID") |>
  mutate(
    # Area of each county polygon in km^2 (plain number, no units attached)
    area_km2 = as.numeric(set_units(st_area(geometry), km^2)),
    # Tornadoes per 1000 km^2
    tdens = tcnt / area_km2 * 1000
  )

I map density with a sequential palette (YlOrRd, darker means higher density) and one facet per year. All facets share one color scale, so colors can be compared directly across years.

ggplot(county_density) +
  geom_sf(aes(fill = tdens)) +
  # One small map per year
  facet_wrap(vars(yr), ncol = 2) +
  scale_fill_distiller(
    name = "Tornadoes/1000 km\u00b2",
    palette = "YlOrRd",
    direction = 1
  ) +
  coord_sf(datum = NA) +
  theme_void()

Exercise 3: Number of classes in quantile choropleth maps

For this exercise I need density for the whole 2016-2021 period, with one value per county.

countymap <- okcounty |>
  left_join(
    # Total tornadoes per county across all years
    countypnt |> st_drop_geometry() |> count(GEOID, name = "tcnt"),
    by = "GEOID"
  ) |>
  # Counties with no tornadoes are missing from the counts
  replace_na(list(tcnt = 0)) |>
  mutate(
    area_km2 = as.numeric(set_units(st_area(geometry), km^2)),
    tdens = tcnt / area_km2 * 1000
  )

Two small helper functions avoid repeating code for the four maps. make_classes() splits a variable into n classes using quantile breaks, so each class holds about the same number of counties. class_map() draws the choropleth for a chosen number of classes.

# Cut a numeric variable into n quantile classes
make_classes <- function(x, n) {
  breaks <- quantile(x, probs = seq(0, 1, length.out = n + 1), na.rm = TRUE)
  # unique() guards against duplicate breaks (many counties may share a value)
  cut(x, breaks = unique(breaks), include.lowest = TRUE)
}

# Draw a choropleth of tornado density with n quantile classes
class_map <- function(n) {
  countymap |>
    mutate(density_class = make_classes(tdens, n)) |>
    ggplot() +
    geom_sf(aes(fill = density_class)) +
    scale_fill_brewer(palette = "YlOrRd") +
    labs(fill = "Density") +
    theme_void()
}
plot_grid(
  class_map(3), class_map(4), class_map(5), class_map(6),
  labels = c("A) 3 Classes", "B) 4 Classes", "C) 5 Classes", "D) 6 Classes"),
  label_size = 12,
  ncol = 2,
  align = "hv"
)

Discussion of results

Exercise 1: where and when. Both the point and the path maps show tornadoes concentrated in central and eastern Oklahoma, with far fewer in the western Panhandle. The yearly colors show that no single year dominates the whole state. Paths are mostly short, so at this scale many appear as small dashes, and the point map is easier to read for overall location. Using the same year colors in both maps lets readers compare the two views of the same events.

Exercise 2: density by year. The yearly pattern varies a lot. There are 57, 85, 43, 148, 38 and 63 tornadoes in 2016-2021, so 2019 stands out as the most active year and has the darkest map; the highest county-year density is Tulsa County in 2019. In contrast, 2018 and 2020 are quiet years in which more than half of Oklahoma’s 77 counties had no tornado at all (56 and 52 counties). Because each county-year density comes from a small count divided by the county’s area, a single tornado can change a county’s color noticeably, so I would not draw conclusions about one county in one year.

Exercise 3: number of classes. All four maps show the same broad pattern: higher density across central and eastern Oklahoma and in parts of the southwest, and lower density in the Panhandle and parts of the southeast. What changes is how much detail the map shows and how the classes are read:

  • With 3 classes, the map is the easiest to read, but each class covers a wide range of density values (the upper class starts at about 3.2 and goes to 10.4 tornadoes per 1000 km²), so counties with quite different densities share a color.
  • With 4 and 5 classes, the middle of the distribution is separated into more groups and the top class gets narrower (starting at about 3.7 and 3.8).
  • With 6 classes, the map is the most detailed (the top class starts at about 4.0), but neighboring colors are harder to tell apart, and small differences between counties look more important than they may be.

Because the breaks are quantiles, each class holds about the same number of counties, so the top class always contains the top 1/n of counties no matter how large the gaps between values are. A class boundary can also fall between counties with almost identical values. More classes are therefore not automatically better: the choice is a tradeoff between detail and readability, and it changes what a reader takes away from the same data.

Caveats. Each tornado is assigned to the county where it touched down, even though its path may cross county lines. Density is a count divided by area, so small counties such as Tulsa can rank high with only a modest number of tornadoes. The study period is six years, which is short for describing a long-term climatology.