Introduction

This analysis uses historical tornado data to examine tornado activity across Oklahoma counties from 2016 to 2021. Tornado points and paths are mapped by year to show where tornadoes occurred and how their locations varied over time. County-level tornado density is then calculated and mapped to compare tornado activity between counties. Finally, four choropleth maps using different numbers of quantile classes are created to explore how classification choices affect the way tornado density patterns are displayed and interpreted.

Load packages

library(sf)
library(ggplot2)
library(dplyr)        
library(tidyr)        
library(scales)       
library(RColorBrewer) 
library(units)
library(cowplot)

Import datasets

# Import Oklahoma county boundaries 
okcounty <- st_read("ok_counties.shp", quiet = TRUE)
# Import historical tornado datasets 
tpoint <- st_read("ok_tornado_point.shp", quiet = TRUE)
tpath <- st_read("ok_tornado_path.shp", quiet = TRUE)

Modify the point and path datasets to filter to a smaller year range

# Filter tornado datasets to 2016-2021
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)

2.A. Map tornado paths and points by year

Step 1. Generate a map of tornado paths where the paths from each year are displayed as a different color

# Map tornado paths 
pathmap <- ggplot() +
  geom_sf(data = tpath_16_21, 
          aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  ggtitle("OK Tornado Paths 2016-2021") +
  theme(plot.title = element_text(
    hjust = 0.5, 
    size = 12, 
    face = "bold"
  ))

Step 2. Generate a map of tornado points where the points from each year are displayed as a different color

# Map tornado points 
pointmap <- 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) +
  ggtitle("OK Tornado Points 2016-2021") + 
  theme(plot.title = element_text(
    hjust = 0.5, 
    size = 12, 
    face = "bold"
  ))

Step 3. Create a composite figure containing the map of path and and the map of points

# Composite figure of tornado points and paths 
plot_grid(
  pathmap, 
  pointmap,
  ncol = 1, 
  hjust = 0, 
  label_x = 0, 
  align = "hv")

Figure 1

The above maps show the spatial distribution of tornado paths and tornado points across Oklahoma from 2016 to 2021. The path map illustrates the routes tornadoes traveled, while the point map shows their recorded locations. Comparing the maps also highlights how tornado paths can extend across multiple counties, while tornado points represent individual recorded locations.

2.B Calculate county-level tornado density by year

# Calculate county area 
okcounty <- okcounty %>%
  mutate(area = st_area(okcounty))

# Spatial join tornado points to counties 
countypnt <- st_join(tpoint_16_21, okcounty)

# Remove geometry for summarizing
countypnt <- st_drop_geometry(countypnt)

# Count tornado points by county and year
countysum_yr <- countypnt %>%
  group_by(GEOID, yr) %>%
  summarize(tcnt = n(), .groups = "drop")

# Create complete county-year combinations, including zero counts 
county_yr <- countysum_yr %>%
  complete(
    GEOID = okcounty$GEOID,
    yr = 2016:2021,
    fill = list(tcnt = 0)
  )

# Join tornado summaries to county polygons
countymap_yr <- okcounty %>%
  left_join(county_yr, by = "GEOID") %>%
  mutate(tcnt = replace_na(tcnt, 0)) %>% 
  mutate(tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units() %>% 
  filter(!is.na(yr))

# Create a faceted map to show all years 
ggplot(data = countymap_yr) +
  geom_sf(aes(fill = tdens), color = "black", linewidth = 0.2) +
  facet_wrap(vars(yr), ncol = 2) +
  scale_fill_distiller(name = expression("Tornado Density (3per 1000 km^2)"),
                       palette = "YlOrRd",
                       breaks = pretty_breaks(),
                       direction = 1) +
  coord_sf(datum = NA) +
  labs(
    title = "Tornado Density by County 2016-2021") +
  theme_void()

Figure 2

The maps above show tornado density by county from 2016 to 2021. Tornado density varies between counties and years, with higher concentrations appearing in central and eastern Oklahoma. 2019 shows particularly noticeable areas of higher density compared to other years.

2.C. Choropleth maps of tornado density based on quantile breaks

Step 1. Calculate the total tornado density per county

# Count tornado points by county across all years
countysum <- countypnt %>%
  group_by(GEOID) %>%
  summarize(tcnt = n(), .groups = "drop")

# Join overall tornado counts to county polygons
countymap <- okcounty %>%
  left_join(countysum, by = "GEOID") %>%
  mutate(tcnt = replace_na(tcnt, 0)) %>%
  mutate(tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

Step 2. Create a function for quantile classification

classify_density <- function(data, numclas) {
  qbrks <- seq(0, 1, length.out = numclas + 1)
  
  data %>%
    mutate(class = cut(
      tdens,
      breaks = unique(quantile(tdens, probs = qbrks, na.rm = TRUE)),
      include.lowest = TRUE
    ))
}

Step 3. Generate maps with 3,4,5, and 6 quantile classes

# 3 classes
countymap_3 <- classify_density(countymap, 3)

# 4 classes
countymap_4 <- classify_density(countymap, 4)

# 5 classes
countymap_5 <- classify_density(countymap, 5)

# 6 classes
countymap_6 <- classify_density(countymap, 6)

Map 1: 3 classes

p1 <- ggplot(countymap_3) +
  geom_sf(aes(fill = class), color = "black") +
  scale_fill_brewer(
    name = expression("Tornadoes/1000 km"^2),
    palette = "YlGnBu"
  ) +
  ggtitle("3 Classes") +
  theme_void() +
  theme(legend.position = "right")

Map 2: 4 classes

p2 <- ggplot(countymap_4) +
  geom_sf(aes(fill = class), color = "black") +
  scale_fill_brewer(
    name = expression("Tornadoes/1000 km"^2),
    palette = "YlGnBu"
  ) +
  ggtitle("4 Classes") +
  theme_void() +
  theme(legend.position = "right")

Map 3: 5 classes

p3 <- ggplot(countymap_5) +
  geom_sf(aes(fill = class), color = "black") +
  scale_fill_brewer(
    name = expression("Tornadoes/1000 km"^2),
    palette = "YlGnBu"
  ) +
  ggtitle("5 Classes") +
  theme_void() +
  theme(legend.position = "right")

Map 4: 6 classes

p4 <- ggplot(countymap_6) +
  geom_sf(aes(fill = class), color = "black") +
  scale_fill_brewer(
    name = expression("Tornadoes/1000 km"^2),
    palette = "YlGnBu"
  ) +
  ggtitle("6 Classes") +
  theme_void() +
  theme(legend.position = "right")

Composite figure containing four classification maps

plot_grid(
  p1, p2, p3, p4,
  ncol = 2,
  align = "hv"
)

Figure 3

Increasing the number of classes creates more differentiation between counties with similar tornado densities. The 3-class map provides a simpler overview of spatial patterns, while the 6-class map reveals more detailed differences but can make the map harder to interpret.