Part A: Mapping tornadoes in Oklahoma with sf and ggplot2

In this part I use the Oklahoma county boundaries together with NOAA Storm Prediction Center tornado data (touchdown points and tornado paths) from Wimberly’s (2023) chapter on vector geospatial data. I restrict the analysis to tornadoes from 2016 to 2021.

Setup and data

I load the packages and read the three shapefiles with st_read(). Each shapefile becomes an sf object, which is a data frame with a special geometry column. That means the usual dplyr verbs work on spatial data.

library(sf)
library(ggplot2)
library(dplyr)
library(tidyr)
library(units)
library(cowplot)
library(leaflet)
# 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 full tornado data cover 1950-2021, so I filter both the point and path layers to 2016-2021 and keep only the columns I need.

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

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

This leaves 434 tornado touchdown points. Most years have between about 40 and 85 tornadoes, and 2019 stands out with 148.

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 paths by year, I map the yr column to color. I convert it to a factor with as.factor() so that ggplot uses a discrete color scale for the year instead of a continuous one. I use the same approach for the points. County boundaries are drawn last so they sit on top of the tornado layers, and theme_void() removes the axes so that only the map is shown.

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

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 combines the two maps into one figure with labels. The hjust and label_x arguments left-align the labels so they do not cover the map.

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"
)

Both maps show tornadoes concentrated in central and eastern Oklahoma, with far fewer in the western Panhandle. Most paths are short, so at this scale many look like small dashes. Paths and points use the same year colors, which makes it easy to compare the two views of the same events.

Exercise 2: County-level tornado density by year

This exercise needs a table with one row per county and year. I first spatially join the tornado points to the counties with st_join(), then count tornadoes for each county-year combination.

countypnt <- st_join(tpoint_16_21, okcounty)

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, I use expand_grid() to build every possible county-year combination, join the counts, and set the missing counts to zero.

county_year <- expand_grid(GEOID = okcounty$GEOID, yr = 2016:2021) |>
  left_join(countysum_year, by = c("GEOID", "yr")) |>
  replace_na(list(tcnt = 0))

Next I join this table to the county polygons and compute density. st_area() returns square meters, so I convert to square kilometers and express density as tornadoes per 1000 km².

county_density <- okcounty |>
  left_join(county_year, by = "GEOID") |>
  mutate(area_km2 = as.numeric(set_units(st_area(geometry), km^2)),
         tdens = tcnt / area_km2 * 1000)

I map density with a sequential YlOrRd palette (darker means higher density), with one facet for each year. All facets share one color scale so they can be compared directly.

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

The pattern changes a lot from year to year. 2019 has the most tornadoes (148) and the darkest maps, with the highest county-year density in Tulsa County. Years like 2018 and 2020 have many counties at or near zero. Because a county’s density depends on a small count, divided by its area, single events can move a county’s color a lot, so I would not read much into one county in one year.

Exercise 3: Number of classes in quantile choropleth maps

For this exercise I first compute county density for the full 2016-2021 period.

countymap <- okcounty |>
  left_join(
    countypnt |> st_drop_geometry() |> count(GEOID, name = "tcnt"),
    by = "GEOID"
  ) |>
  replace_na(list(tcnt = 0)) |>
  mutate(area_km2 = as.numeric(set_units(st_area(geometry), km^2)),
         tdens = tcnt / area_km2 * 1000)

I then write two small helper functions. make_classes() cuts 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.

make_classes <- function(x, n) {
  breaks <- quantile(x, probs = seq(0, 1, length.out = n + 1), na.rm = TRUE)
  cut(x, breaks = unique(breaks), include.lowest = TRUE)
}

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"
)

All four maps show the same broad pattern: higher tornado 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.

  • With 3 classes, the map is easy to read, but each class covers a wide range of density values, so counties with fairly different densities share a color.
  • With 4 and 5 classes, the middle of the distribution separates, and the highest-density counties form a tighter group.
  • With 6 classes, the map is the most detailed, but neighboring colors become hard to tell apart, and small differences between counties look more important than they may be.

Because the breaks are quantiles, a class boundary can fall between counties with almost identical values, and the top class always spans a very large range (the highest density is about 10.4 tornadoes per 1000 km²). A map with more classes is not automatically better: it is a tradeoff between detail and readability, and the choice of classes can change what a reader takes away.

Part B: An interactive map of Oklahoma tornadoes

For Part B I chose the same tornado point data, this time as an interactive leaflet map. This data set suits the lab because each tornado is a point with a magnitude on the Enhanced Fujita (EF) scale, which is an ordered variable.

Preparing the data

I keep the columns I need: date, magnitude (mag), injuries (inj), and fatalities (fat). In the NOAA data a magnitude of -9 means the magnitude is unknown, so I recode it to NA. If I left it as -9, it would be treated as the lowest value and would distort the color scale. About a quarter of the tornadoes (105 of 434) are in this unknown group, so I keep them on the map and show them in gray instead of dropping them. Leaflet needs WGS84 longitude and latitude, so I also transform the points with st_transform(4326).

tornado_map <- tpoint |>
  filter(yr >= 2016, yr <= 2021) |>
  select(om, yr, date, mag, inj, fat) |>
  mutate(mag = na_if(mag, -9)) |>
  st_transform(4326)

table(tornado_map$mag, useNA = "ifany")
## 
##    0    1    2    3    4 <NA> 
##  118  178   27    5    1  105

Palette

Magnitude is ordered, so I use a sequential palette. colorFactor() gives each of the five magnitudes (0-4) its own color from the YlOrRd ramp, so darker and redder means a stronger tornado. Unknown magnitude is shown in gray, which is neutral and not part of the ramp.

pal <- colorFactor(palette = "YlOrRd", domain = 0:4, na.color = "gray")

County layer

For context, I add county outlines. The county layer reuses the tornado counts per county from Part A, and I transform it to WGS84 to match.

okcounty_ll <- st_transform(countymap, 4326)

The map

The map is built in layers:

  1. A light gray basemap (Esri.WorldGrayCanvas) keeps the colors of the points easy to see.
  2. County polygons with a transparent fill are drawn first, so they sit underneath the points. Hovering over a county shows its name and its tornado count.
  3. Tornado points are drawn as circle markers. The inside is colored by magnitude, with a thin black outline so that pale points stay visible. Hovering shows the magnitude, and clicking opens a popup with the date, magnitude, injuries, and fatalities.
  4. A legend, a title, and a starting view that fits the whole state complete the map.
  5. A layers control in the top right lets readers switch the county outlines and the tornado points on and off. Each layer is given a group name, and addLayersControl() lists those groups as checkboxes.
leaflet(tornado_map) |>
  addProviderTiles(providers$Esri.WorldGrayCanvas) |>
  addPolygons(data = okcounty_ll,
              group = "County outlines",
              fillOpacity = 0,
              color = "gray20",
              weight = 1.5,
              opacity = 0.8,
              label = ~paste0(NAME, " County: ", tcnt,
                              ifelse(tcnt == 1, " tornado", " tornadoes")),
              highlightOptions = highlightOptions(weight = 3, color = "black")) |>
  addCircleMarkers(
    group = "Tornadoes",
    fillColor = ~pal(mag),
    color = "black",
    weight = 0.5,
    radius = 4,
    fillOpacity = 0.9,
    label = ~paste0("Magnitude: ", ifelse(is.na(mag), "Unknown", mag)),
    popup = ~paste0("Date: ", date, "<br>",
                    "Magnitude: ", ifelse(is.na(mag), "Unknown", mag), "<br>",
                    "Injuries: ", inj, "<br>",
                    "Fatalities: ", fat)
  ) |>
  addLegend(position = "bottomleft",
            pal = pal, values = tornado_map$mag,
            title = "Tornado magnitude (EF)",
            na.label = "Unknown") |>
  addLayersControl(overlayGroups = c("County outlines", "Tornadoes"),
                   options = layersControlOptions(collapsed = FALSE)) |>
  addControl("Oklahoma tornadoes, 2016-2021", position = "bottomright") |>
  fitBounds(lng1 = -103.1, lat1 = 33.6, lng2 = -94.4, lat2 = 37.1)

Results

The map shows the same geographic pattern as the static maps in Part A: tornadoes are concentrated in central and eastern Oklahoma, with fewer in the western Panhandle. Most tornadoes with a known magnitude are weak (EF0 and EF1), shown in pale yellow and light orange. The strongest tornadoes (EF3 and EF4) are rare, with only six in this period, and stand out as the darker red points. The gray points have unknown magnitude and should not be read as weak. Overlap among points in the densest areas can hide some of them, so zooming in helps. Hovering over a county shows its total count for 2016-2021, including tornadoes of unknown magnitude.