Introduction

Oklahoma gets a lot of tornadoes, but they don’t happen evenly across the state or evenly from year to year. In this report, I mapped Oklahoma tornadoes from 2016 to 2021 to show where tornadoes happened, how much activity changed between years, and how the way data is grouped into classes can change what a map tells the reader.

This report follows the three practice exercises at the end of Wimberly’s (2023) chapter on Vector Geospatial Data.

Data

This report uses three datasets, all in shapefile format:

The tornado data come from the NOAA National Weather Service Storm Prediction Center. All three datasets were provided with the chapter data from Wimberly (2023).

Setup

Load packages

These are the R packages used in this report. Packages add extra tools to R.

library(sf)       # read and work with spatial data
library(ggplot2)  # make maps and plots
library(dplyr)    # filter, select, and summarize data
library(tidyr)    # replace missing values
library(units)    # remove units from area values
library(cowplot)  # combine multiple maps into one figure

Import the data

Read in the county boundaries, tornado points, and tornado paths.

# Read in the three shapefiles from the data folder
okcounty <- st_read(here::here("data", "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(here::here("data", "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(here::here("data", "ok_tornado_path.shp"), quiet = TRUE)

Filter to 2016-2021

The tornado data covers many years, but this report only looks at 2016 to 2021. I kept three columns: the tornado ID (om), the year (yr), and the date. The location of each tornado was kept automatically.

# Keep tornado points from 2016-2021
tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

# Keep tornado paths from 2016-2021
tpath_16_21 <- tpath %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

Practice 1: Tornado Paths and Points by Year

In the first exercise, I compared two ways of showing the same tornadoes: as paths and as touchdown points. Each year was given its own color so the years could be told apart.

Make the maps

Make one map of tornado paths and one map of tornado points. The county boundaries are drawn as outlines on top.

# Map tornado paths, colored by year
pathmap <- ggplot() +
  geom_sf(data = tpath_16_21, aes(color = as.factor(yr))) +  # paths
  geom_sf(data = okcounty, fill = NA, color = "gray60") +    # county outlines
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +                                     # remove grid lines
  theme_void() +                                             # remove axes
  theme(text = element_text(color = "white")) +
  ggtitle("Tornado Paths")

# Map tornado points, colored by year
pointmap <- ggplot() +
  geom_sf(data = tpoint_16_21, aes(color = as.factor(yr))) + # points
  geom_sf(data = okcounty, fill = NA, color = "gray60") +    # county outlines
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void() +
  theme(text = element_text(color = "white")) +
  ggtitle("Tornado Touchdown Points")

Combine the maps

Put the two maps side by side with plot_grid() and add an overall title on top.

# Put the path and point maps side by side
maps_p1 <- plot_grid(
  pathmap,
  pointmap,
  labels = c("A", "B"),
  label_colour = "white",
  ncol = 2
)

# Make an overall title for the figure
title_p1 <- ggdraw() +
  draw_label(
    "Oklahoma Tornadoes by Year, 2016-2021",
    colour = "white",
    fontface = "bold",
    size = 16
  )

# Stack the title on top of the maps
plot_grid(title_p1, maps_p1, ncol = 1, rel_heights = c(0.1, 1))
Figure 1. Oklahoma tornadoes from 2016 to 2021, colored by year. (A) Tornado paths, showing the track of each tornado after touchdown. (B) Tornado touchdown points. County boundaries are shown in gray. Data: NOAA Storm Prediction Center, via Wimberly (2023).

Figure 1. Oklahoma tornadoes from 2016 to 2021, colored by year. (A) Tornado paths, showing the track of each tornado after touchdown. (B) Tornado touchdown points. County boundaries are shown in gray. Data: NOAA Storm Prediction Center, via Wimberly (2023).

Results

Map A shows the tornado paths and Map B shows the tornado points. Most tornadoes were in central and eastern Oklahoma, with fewer in the panhandle. Most paths were short, and the longer ones mostly went from southwest to northeast. The points are easier to see, but the paths show how far and in what direction each tornado moved. 2019 had a lot of tornadoes compared to the other years.

Practice 2: Tornado Density by County and Year

In the second exercise, I looked at how tornado activity changed from year to year. Instead of mapping each tornado, I calculated tornado density for each county. Density is the number of tornadoes divided by the county’s area, so large and small counties can be compared fairly.

Count tornadoes by county and year

A spatial join matches each tornado point to the county it falls inside. After the join, I counted the tornadoes in each county for each year.

# Match each tornado point to its county
countypnt <- st_join(tpoint_16_21, okcounty)

# Remove the point locations since only the counts are needed now
countypnt <- st_drop_geometry(countypnt)

# Count tornadoes in each county for each year
countysum <- countypnt %>%
  group_by(GEOID, yr) %>%
  summarize(tcnt = n()) %>%
  ungroup()

Add years with no tornadoes

Some counties had no tornadoes in some years. To keep those counties on the map, I made a copy of every county for each year.

# Repeat the counties 6 times, once for each year
countymap <- okcounty[rep(seq_len(nrow(okcounty)), times = 6), ]

# Add the year to each copy
countymap$yr <- rep(2016:2021, each = nrow(okcounty))

Calculate density

Join the counts to the counties, set counties with no tornadoes to 0, and calculate density as tornadoes per 1,000 km².

countymap <- countymap %>%
  left_join(countysum, by = c("GEOID", "yr")) %>%  # add counts to counties
  mutate(tcnt = replace_na(tcnt, 0)) %>%           # no tornadoes = 0
  mutate(
    area = st_area(geometry),                      # county area in m²
    tdens = 10^6 * 10^3 * tcnt / area              # tornadoes per 1,000 km²
  ) %>%
  drop_units()                                     # remove units for mapping

Map density by year

Make one map for each year using facet_wrap(). All six maps use the same color scale so the years can be compared.

# Map tornado density with one panel per year
ggplot(data = countymap) +
  geom_sf(aes(fill = tdens)) +
  facet_wrap(vars(yr), ncol = 3) +
  scale_fill_distiller(
    name = expression("Tornadoes/1000 km" ^ 2),
    palette = "YlOrRd",
    direction = 1
  ) +
  coord_sf(datum = NA) +
  theme_void() +
  theme(
    text = element_text(color = "white"),
    plot.title = element_text(           # center and bold the title
      hjust = 0.5,
      face = "bold",
      size = 16,
      margin = margin(b = 15)            # add space below the title
    )
  ) +
  ggtitle("Oklahoma Tornado Density by County and Year, 2016-2021")
Figure 2. Tornado density by county in Oklahoma for each year from 2016 to 2021, shown as tornadoes per 1,000 km². Darker red means more tornadoes for the county's size. All years use the same color scale. Data: NOAA Storm Prediction Center, via Wimberly (2023).

Figure 2. Tornado density by county in Oklahoma for each year from 2016 to 2021, shown as tornadoes per 1,000 km². Darker red means more tornadoes for the county’s size. All years use the same color scale. Data: NOAA Storm Prediction Center, via Wimberly (2023).

Results

2019 had the highest tornado density, mostly in central and northeastern Oklahoma, with one county over 5 tornadoes per 1,000 km². 2016 and 2020 had the lowest, with most counties having few or no tornadoes. The other years fell in between. This shows that tornado activity can change a lot from one year to the next.

Practice 3: Changing the Number of Classes

In the third exercise, I looked at how the number of classes changes a map. I calculated tornado density for the whole 2016-2021 period and grouped the counties using quantile breaks. Quantile breaks split the counties into groups with about the same number of counties in each group. I made four maps using 3, 4, 5, and 6 classes.

Calculate density for 2016-2021

Count tornadoes in each county for all six years together and calculate density.

# Count tornadoes in each county for all years combined
countysum <- countypnt %>%
  group_by(GEOID) %>%
  summarize(tcnt = n()) %>%
  ungroup()

# Join counts to counties and calculate density per 1,000 km²
countymap <- okcounty %>%
  left_join(countysum, by = "GEOID") %>%
  mutate(tcnt = replace_na(tcnt, 0)) %>%
  mutate(
    area = st_area(geometry),
    tdens = 10^6 * 10^3 * tcnt / area
  ) %>%
  drop_units()

Create the classes

Divide density into 3, 4, 5, and 6 quantile classes.

# quantile() finds the break points and cut() puts each county in a class
countymap <- countymap %>%
  mutate(
    tdens_c3 = cut(
      tdens,
      breaks = quantile(tdens, probs = seq(0, 1, length.out = 4)),
      include.lowest = TRUE
    ),
    tdens_c4 = cut(
      tdens,
      breaks = quantile(tdens, probs = seq(0, 1, length.out = 5)),
      include.lowest = TRUE
    ),
    tdens_c5 = cut(
      tdens,
      breaks = quantile(tdens, probs = seq(0, 1, length.out = 6)),
      include.lowest = TRUE
    ),
    tdens_c6 = cut(
      tdens,
      breaks = quantile(tdens, probs = seq(0, 1, length.out = 7)),
      include.lowest = TRUE
    )
  )

Make the maps

Make one map for each number of classes.

# Map with 3 classes
map3 <- ggplot(data = countymap) +
  geom_sf(aes(fill = tdens_c3)) +
  scale_fill_brewer(name = "Tornadoes/1000 km²", palette = "YlOrRd") +
  coord_sf(datum = NA) +
  theme_void() +
  theme(text = element_text(color = "white")) +
  ggtitle("3 Classes")

# Map with 4 classes
map4 <- ggplot(data = countymap) +
  geom_sf(aes(fill = tdens_c4)) +
  scale_fill_brewer(name = "Tornadoes/1000 km²", palette = "YlOrRd") +
  coord_sf(datum = NA) +
  theme_void() +
  theme(text = element_text(color = "white")) +
  ggtitle("4 Classes")

# Map with 5 classes
map5 <- ggplot(data = countymap) +
  geom_sf(aes(fill = tdens_c5)) +
  scale_fill_brewer(name = "Tornadoes/1000 km²", palette = "YlOrRd") +
  coord_sf(datum = NA) +
  theme_void() +
  theme(text = element_text(color = "white")) +
  ggtitle("5 Classes")

# Map with 6 classes
map6 <- ggplot(data = countymap) +
  geom_sf(aes(fill = tdens_c6)) +
  scale_fill_brewer(name = "Tornadoes/1000 km²", palette = "YlOrRd") +
  coord_sf(datum = NA) +
  theme_void() +
  theme(text = element_text(color = "white")) +
  ggtitle("6 Classes")

Combine the maps

Put the four maps together with plot_grid() and add an overall title on top.

# Put the four maps in a 2 x 2 grid
maps_p3 <- plot_grid(
  map3,
  map4,
  map5,
  map6,
  labels = c("A", "B", "C", "D"),
  label_colour = "white",
  ncol = 2
)

# Make an overall title for the figure
title_p3 <- ggdraw() +
  draw_label(
    "Oklahoma Tornado Density by Number of Quantile Classes, 2016-2021",
    colour = "white",
    fontface = "bold",
    size = 16
  )

# Stack the title on top of the maps
plot_grid(title_p3, maps_p3, ncol = 1, rel_heights = c(0.06, 1))
Figure 3. Tornado density by county in Oklahoma from 2016 to 2021 (tornadoes per 1,000 km²), grouped using quantile breaks with (A) 3, (B) 4, (C) 5, and (D) 6 classes. The data is the same in all four maps; only the number of classes changes. Data: NOAA Storm Prediction Center, via Wimberly (2023).

Figure 3. Tornado density by county in Oklahoma from 2016 to 2021 (tornadoes per 1,000 km²), grouped using quantile breaks with (A) 3, (B) 4, (C) 5, and (D) 6 classes. The data is the same in all four maps; only the number of classes changes. Data: NOAA Storm Prediction Center, via Wimberly (2023).

Results

With 3 classes, the map shows a simple high, medium, and low pattern, but the top class covers a big range (3.21 to 10.4), so very different counties end up the same color. With more classes, the map shows more detail, and the darkest color only shows the highest counties in central, northeastern, and southwestern Oklahoma. With 6 classes, the colors start to get hard to tell apart. 4 or 5 classes seemed like the best balance between detail and readability.

Discussion

These maps show that tornadoes in Oklahoma from 2016 to 2021 were most common in central and northeastern parts of the state, and least common in the panhandle. Tornado activity also changed a lot from year to year, with 2019 being much more active than the other years. The last exercise shows that map design choices matter. The same data can look different depending on how many classes are used, so it’s important to choose classes that show the pattern clearly without hiding or exaggerating differences.

References

Wimberly, M. C. (2023). Geographic data science with R: Visualizing and analyzing environmental change. Chapman & Hall/CRC. https://bookdown.org/mcwimberly/gdswr-book/

NOAA National Weather Service Storm Prediction Center. Severe weather GIS (SVRGIS) data. https://www.spc.noaa.gov/gis/svrgis/