Introduction

Oklahoma sits in the middle of Tornado Alley, and between 2016 and 2021 the NOAA Storm Prediction Center recorded 434 tornadoes touching down in the state. We can do a lot more with this data than just count the tornadoes, though: a count alone doesn’t show where in the state tornadoes are concentrated, whether that pattern is consistent from one year to the next, or which paths they take. This report works through those three themes with the practice exercises from Chapter 5 of Wimberly’s (2023) Geographic Data Science with R.

The first exercise maps each tornado’s path and touchdown point, colored by year, to show where tornadoes occurred. The second converts touchdown points into density for every county and year to show how the pattern changes over time. The third holds the data fixed and changes only the number of classes on a choropleth map, to demonstrate the effect of symbology on a map’s perception.

Data

All three layers come with the Chapter 5 materials and are stored as ESRI shapefiles in NAD83 geographic coordinates (longitude and latitude):

  • Counties: the 77 Oklahoma county boundaries from the U.S. Census Bureau’s cartographic boundary files.
  • Tornado touchdown points: one point per tornado where it first touched down, from the NOAA Storm Prediction Center (SPC) Severe Weather Database from 1950-2021.
  • Tornado paths: one line per tornado from its touchdown point to where it lifted off, also from the SPC database. Each path shares an ID (om) with its touchdown point.

Each row in the point and path layers represents one tornado. Both layers also carry the SPC attributes, like the date, EF-scale magnitude, injuries, and the recorded path length in miles (len).

Data preparation

Load packages and data

Besides the packages from earlier labs, sf reads, joins, and measures spatial data, and cowplot provides plot_grid() for combining several maps into one figure. Every map uses theme_void(), which removes the axes and gray background, so in order to avoid redundancy I set it once for the whole report with theme_set() instead of adding it to each map individually.

library(tidyverse) # data wrangling and ggplot2
library(sf) # simple features: spatial vector data
library(cowplot) # plot_grid() to combine maps
library(here) # builds file paths from project root

# set one theme for every map in the report: no axes or background, centered titles
theme_set(
  theme_void() +
    theme(
      plot.title  = element_text(hjust = 0.5),
      plot.margin = margin(10, 10, 10, 10)
    )
)
counties   <- st_read(here("GEOG588", "lab4", "data", "Chapter5", "ok_counties.shp"), quiet = TRUE)
points_all <- st_read(here("GEOG588", "lab4", "data", "Chapter5", "ok_tornado_point.shp"), quiet = TRUE)
paths_all  <- st_read(here("GEOG588", "lab4", "data", "Chapter5", "ok_tornado_path.shp"), quiet = TRUE)

Keep 2016-2021

Like the chapter, I kept the six most recent years and only the columns the maps need. I found the chapter’s variable names difficult to read, and select() can rename columns as it keeps them, so the SPC’s short names become tornado_id and year. For the paths, I also kept the SPC’s recorded path length, which is needed in the next step.

tornado_points <- points_all |>
  filter(yr >= 2016 & yr <= 2021) |>
  select(tornado_id = om, year = yr, date)

tornado_paths <- paths_all |>
  filter(yr >= 2016 & yr <= 2021) |>
  select(tornado_id = om, year = yr, date, recorded_miles = len)

That leaves 434 tornadoes in each layer.

Check the path geometries

When I coded along with the chapter, the map of tornado paths strangely looked more like a map of points than lines, although sf does report every path as a line:

st_geometry_type(tornado_paths, by_geometry = FALSE)
## [1] LINESTRING
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING MULTIPOLYGON GEOMETRYCOLLECTION ... TRIANGLE

The geometry type alone doesn’t tell the whole story, however. If I use st_length() to measure each line, it shows that some paths have no length at all:

tornado_paths <- tornado_paths |>
  mutate(path_km = as.numeric(st_length(geometry)) / 1000) # st_length() returns meters

# paths whose start and end points are the same
zero_length <- tornado_paths |>
  filter(path_km == 0)

zero_length |>
  st_drop_geometry() |>
  count(year, name = "zero_length_paths") |>
  knitr::kable(col.names = c("Year", "Paths with zero length"))
Year Paths with zero length
2016 4
2017 21
2018 8
2019 47
2020 7
2021 13

In total, 100 of the 434 paths (23%) are lines whose end point is the same as their start point. The SPC still recorded a length for these tornadoes (a median of 0.2 miles), so they did travel, but no end location was stored, so the line ends where it starts. Because a line with no length can’t show a path, I left these 100 tornadoes off the paths map, but they still appear on the touchdown points map and in every density calculation.

Of the rest of the paths some are short. Half of them are under 4.6 km long, and the longest is 75 km, compared to a state that is roughly 750 km wide. At this scale, most paths look like short dashes, which is why the chapter’s path map can still look like there are points.

paths_with_length <- tornado_paths |>
  filter(path_km > 0)

County area and the county of each tornado

Exercises 2 and 3 both need each county’s area and the county each tornado touched down in, so I computed them ahead of time here. st_area() returns square meters, so dividing by 1,000,000 gives square kilometers. The data are in longitude and latitude, and sf measures area on the curved surface of the earth in that case, so the result is accurate without reprojecting first.

counties <- counties |>
  mutate(area_km2 = as.numeric(st_area(geometry)) / 1e6)

st_join() is a spatial join: each touchdown point gets the columns of the county polygon it falls inside. After the join I dropped the point geometry, since the counting that follows only requires the table.

tornado_counties <- tornado_points |>
  st_join(counties) |>
  st_drop_geometry()

All 434 touchdown points fall inside a county, so every tornado is counted.

Exercise 1: Where did tornadoes touch down and travel?

The two maps are built the same way as Figure 5.4 in the chapter, with the county outlines drawn in gray so the tornadoes stand out. Since year is an ordered variable, I used the viridis palette instead of the chapter’s default colors so that early and late years read as a sequence instead of six unrelated categories.

# map A: tornado paths, colored by year
paths_map <- ggplot() +
  geom_sf(data = counties, fill = NA, color = "gray70") +
  geom_sf(data = paths_with_length, aes(color = as.factor(year)), linewidth = 0.8) +
  scale_color_viridis_d(name = "Year") +
  coord_sf(datum = NA) + # no graticule lines
  labs(title = "Tornado paths")

# map B: touchdown points, colored by year
points_map <- ggplot() +
  geom_sf(data = counties, fill = NA, color = "gray70") +
  geom_sf(data = tornado_points, aes(color = as.factor(year)), size = 1) +
  scale_color_viridis_d(name = "Year") +
  coord_sf(datum = NA) +
  labs(title = "Tornado touchdown points")
plot_grid(
  paths_map, points_map,
  ncol       = 1,
  labels     = c("A", "B"),
  label_size = 12
)
Figure 1. Tornado paths (A) and touchdown points (B) in Oklahoma, 2016-2021, colored by year. Each line runs from touchdown to lift-off. 100 tornadoes with no recorded end point appear in B only. Data: NOAA Storm Prediction Center; county boundaries from the U.S. Census Bureau.

Figure 1. Tornado paths (A) and touchdown points (B) in Oklahoma, 2016-2021, colored by year. Each line runs from touchdown to lift-off. 100 tornadoes with no recorded end point appear in B only. Data: NOAA Storm Prediction Center; county boundaries from the U.S. Census Bureau.


The two maps show the same tornadoes in different ways: the paths map shows how tornadoes moved, and that most paths run from southwest to northeast, the direction storms in the southern Plains usually travel. The points map is easier to read for where tornadoes happened, since every tornado is the same size regardless of how far it traveled. Overall, tornadoes touched down across the whole state, and 2019 has the most points of any year, at 148 of the 434.

Exercise 2: How did tornado density change from year to year?

Count tornadoes by county and year

The counting works the same way as in the chapter, but grouped by year as well as county. Since grouping only returns combinations that actually occurred in the data, a county with no tornadoes in 2018 has no 2018 row at all. complete() fixes that by adding every county and year combination and fills in the missing counts with 0.

county_year_counts <- tornado_counties |>
  group_by(GEOID, year) |>
  summarize(tornado_count = n(), .groups = "drop") |>
  # add the county-years with no tornadoes, using the full list of counties
  complete(GEOID = counties$GEOID, year = 2016:2021, fill = list(tornado_count = 0))

This gives 462 rows, one for each of the 77 counties in each of the six years, and 53% of these county-years had no tornadoes.

Join the counts to the counties and compute density

Joining the counts onto the county polygons makes each county appear six times, once per year, which is what facet_wrap() needs as an input to draw one map per year. Density is in tornadoes per 1,000 km².

county_year_density <- counties |>
  left_join(county_year_counts, by = "GEOID") |>
  mutate(tornado_density = 1000 * tornado_count / area_km2)

# highest county density in each year, used in the text below
yearly_max <- county_year_density |>
  st_drop_geometry() |>
  group_by(year) |>
  summarize(max_density = max(tornado_density))
ggplot(county_year_density) +
  geom_sf(aes(fill = tornado_density), color = "gray60", linewidth = 0.1) +
  scale_fill_distiller(
    name      = "Tornadoes per 1,000 km²",
    palette   = "YlOrRd",
    direction = 1 # light = low, dark = high
  ) +
  facet_wrap(vars(year), ncol = 3) +
  coord_sf(datum = NA) +
  labs(title = "Tornado density by county and year") +
  theme(legend.position = "bottom", legend.key.width = unit(1.5, "cm"))
Figure 2. Tornado touchdowns per 1,000 km² in each Oklahoma county by year, 2016-2021. All six maps share one color scale, so the same color means the same density in every year. Data: NOAA Storm Prediction Center; county boundaries from the U.S. Census Bureau.

Figure 2. Tornado touchdowns per 1,000 km² in each Oklahoma county by year, 2016-2021. All six maps share one color scale, so the same color means the same density in every year. Data: NOAA Storm Prediction Center; county boundaries from the U.S. Census Bureau.


2019 stands out from every other year. It had 148 tornadoes, more than any other year, and the highest county density of any year, 5.9 tornadoes per 1,000 km². Much of that came from one stretch of storms: 77 tornadoes touched down between May 20th and 25th, 2019, spread across central and northeastern Oklahoma, with the most (7) in Tulsa County. In the other five years, no county went above 3 tornadoes per 1,000 km². Most of each map is pale because most counties have zero or one tornado in a given year. At this level of detail, one tornado can change a county’s color, so a single year’s map says more about that year’s storms than about where tornadoes tend to happen across longer periods of time.

Exercise 3: How does the number of classes change the map?

Density for the full period

This is the same county-level density as section 5.3 of the chapter. Two counties had no tornadoes in any of the six years, so they have no row after the join. replace_na() sets their count to 0.

county_counts <- tornado_counties |>
  group_by(GEOID) |>
  summarize(tornado_count = n())

county_density <- counties |>
  left_join(county_counts, by = "GEOID") |>
  mutate(
    tornado_count   = replace_na(tornado_count, 0),
    tornado_density = 1000 * tornado_count / area_km2
  )

One function for four maps

The four maps are identical except for the number of classes, so instead of copying the same code four times, I wrote it once as a function that takes the number of classes as its input.

The chapter passes the breaks to quantile() with breaks =, but quantile() doesn’t have an argument by that name. As a result R ignores it and returns quartiles, so the chapter’s code gives four classes no matter what number is asked for. The argument that worked is probs =.

# choropleth of tornado density with n_classes quantile classes
quantile_map <- function(n_classes) {
  # e.g. 4 classes -> probabilities 0, 0.25, 0.5, 0.75, 1
  quantile_breaks <- quantile(
    county_density$tornado_density,
    probs = seq(0, 1, length.out = n_classes + 1)
  )

  classed_counties <- county_density |>
    mutate(density_class = cut(
      tornado_density,
      breaks         = quantile_breaks,
      include.lowest = TRUE # keep the lowest value (0) in the first class
    ))

  ggplot(classed_counties) +
    geom_sf(aes(fill = density_class), color = "gray60", linewidth = 0.1) +
    scale_fill_brewer(name = "Tornadoes per\n1,000 km²", palette = "YlOrRd") +
    coord_sf(datum = NA) +
    labs(title = paste(n_classes, "classes"))
}
plot_grid(
  quantile_map(3), quantile_map(4),
  quantile_map(5), quantile_map(6),
  ncol       = 2,
  labels     = "AUTO", # letters the panels A-D
  label_size = 12
)
Figure 3. Tornado touchdowns per 1,000 km² by Oklahoma county, 2016-2021, classified by quantiles into 3 (A), 4 (B), 5 (C), and 6 (D) classes. Each class holds about the same number of counties. Data: NOAA Storm Prediction Center; county boundaries from the U.S. Census Bureau.

Figure 3. Tornado touchdowns per 1,000 km² by Oklahoma county, 2016-2021, classified by quantiles into 3 (A), 4 (B), 5 (C), and 6 (D) classes. Each class holds about the same number of counties. Data: NOAA Storm Prediction Center; county boundaries from the U.S. Census Bureau.


The counties with the highest density over the full period are:

county_density |>
  st_drop_geometry() |>
  slice_max(tornado_density, n = 5) |>
  select(NAME, tornado_count, area_km2, tornado_density) |>
  knitr::kable(digits = 1, col.names = c("County", "Tornadoes", "Area (km²)", "Tornadoes per 1,000 km²"))
County Tornadoes Area (km²) Tornadoes per 1,000 km²
Tulsa 16 1535.6 10.4
Cleveland 10 1424.5 7.0
Oklahoma 13 1860.5 7.0
Tillman 16 2310.7 6.9
Mayes 12 1766.7 6.8

How the number of classes changes the map

Quantile classes put about the same number of counties in each class, which spreads the colors evenly across the map, but this means the classes cover very different ranges of values. The top class shows this most clearly. With 3 classes it starts at 3.21, and with 6 classes it starts at 3.97, but it always runs up to 10.42, Tulsa County’s density. On every map, Tulsa gets the same dark red as counties with less than half its density, so the outlier in the data is hidden inside a wide top class.

Adding classes adds detail mostly at the low end of the range, where most of the counties are. With 3 classes, everything from 0 to 1.68 is one color, and with 6 classes that same range is split into two. Going from 3 to 6 classes makes the map look more varied, but that doesn’t necessarily make it easier to read, since six shades of one color are hard to tell apart and match to the legend. For this data, 4 or 5 classes seems like the best balance between detail and readability.

Discussion

Going back to the three questions from the introduction:

  • Where: tornadoes touched down in 75 of Oklahoma’s 77 counties between 2016 and 2021, so almost no part of the state went untouched. Over the full period, density was highest in Tulsa, Cleveland, and Oklahoma County.
  • Paths: most paths run from southwest to northeast (276 of the 334 paths with a length), the same direction the storms that produce them usually travel.
  • Year to year: the pattern was not consistent from one year to the next. Much of the six-year total came from one year, 2019, and just over half of that year’s tornadoes came from one week of storms in late May. Since most counties had zero or one tornado in a given year, each year’s map reflects that year’s storms more than a lasting pattern.

A few limitations should still be kept in mind when reading these maps:

  • Tulsa, Cleveland, and Oklahoma County contain the Tulsa and Oklahoma City metro areas, and they are also smaller than most counties in the state, so a handful of tornadoes raises their density a lot. Tornado reports are also known to be biased toward populated areas, since more people means more tornadoes are seen and reported (Elsner et al. 2013), so part of their high density may come from reporting rather than from more tornadoes. Reporting can’t be the whole story, however: rural Tillman and Mayes counties rank fourth and fifth, at nearly the same density as Oklahoma County.
  • Density is based on where each tornado first touched down, so a tornado that crossed into a second county only counts in the first.
  • 100 tornadoes had no recorded end point, so the paths map is missing about a quarter of the tornadoes, mostly short ones. This is why it’s worth checking the geometries themselves (and how they visually appear) and not only their type, as every path was a valid line, but not every line had length.

Finally, Exercise 3 shows that the number of classes is not a neutral choice. The same data can look evenly spread or concentrated depending on how it is classified, and with quantile breaks in particular, the most extreme values can end up grouped with much lower ones.

References

  • NOAA Storm Prediction Center. Severe Weather Database: Tornadoes, 1950-2021 (touchdown points and paths). SPC GIS data. Provided with the Chapter 5 materials.
  • U.S. Census Bureau. Cartographic Boundary Files: Counties. Census cartographic boundary files. Provided with the Chapter 5 materials.
  • Wimberly, M. C. 2023. Geographic Data Science with R: Visualizing and Analyzing Environmental Change. CRC Press. Online edition.
  • Elsner, J. B., L. E. Michaels, K. N. Scheitlin, and I. J. Elsner. 2013. The decreasing population bias in tornado reports across the central Plains. Weather, Climate, and Society 5(2):221-232.
  • R packages: tidyverse, sf, cowplot, here, and knitr.