Mapping Tornado Activity in Oklahoma: Spatial Patterns and County-Level Density (2016–2021)

Part A: Static Mapping in R

Introduction

This analysis explores the spatial distribution of tornadoes across Oklahoma from 2016 to 2021 using historical tornado data from NOAA’s Storm Prediction Center and county boundary data. Using the sf and ggplot2 packages in R, three mapping exercises examine tornado touchdown locations and paths, annual variations in county-level tornado density, and the effects of different quantile classification schemes on map interpretation. Together, these exercises demonstrate how spatial data processing and cartographic design choices influence the presentation of geographic patterns.

Data and Packages

The analysis uses Oklahoma county boundaries and historical tornado touchdown and path data. Several R packages are used to import and manipulate spatial data, calculate tornado densities, and create maps. The following code loads the necessary packages and prepares the datasets for analysis.

# Load required packages
library(here)
library(sf)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(RColorBrewer)
library(units)
library(cowplot)

# Import spatial datasets
okcounty <- st_read(here("data", "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(here("data", "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(here("data", "ok_tornado_path.shp"), quiet = TRUE)

# Filter tornado points and paths 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)

Exercise 1: Tornado Touchdown Points and Paths by Year

The first exercise compares the spatial distribution of tornado touchdown points and paths across Oklahoma from 2016 to 2021. Each year is represented by a different color, with consistent colors used across both maps. The two maps are then combined using plot_grid() to compare tornado locations with the geographic extent of their paths.

# Map tornado touchdown points
point_map <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "gray60") +
  geom_sf(data = tpoint_16_21,
          aes(color = as.factor(yr)),
          size = 1) +
  scale_color_brewer(palette = "Dark2", name = "Year") +
  coord_sf(datum = NA) +
  labs(title = "   Tornado Touchdown Points") +
  theme_void() +
  theme(legend.position = "bottom")

# Map tornado paths
path_map <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "gray60") +
  geom_sf(data = tpath_16_21,
          aes(color = as.factor(yr)),
          linewidth = 0.8) +
  scale_color_brewer(palette = "Dark2", name = "Year") +
  coord_sf(datum = NA) +
  labs(title = "   Tornado Paths") +
  theme_void() +
  theme(legend.position = "bottom")

# Combine the two maps
exercise1_map <- plot_grid(
  point_map,
  path_map,
  labels = c("A", "B"),
  ncol = 1,
  align = "v"
)

exercise1_map

Figure 1. Comparison of tornado touchdown points (A) and tornado paths (B) in Oklahoma from 2016 to 2021. Colors distinguish the year of each event, allowing the spatial distribution of tornado origins and paths to be compared.

The two maps reveal similar geographic patterns, with tornado activity appearing particularly concentrated in eastern Oklahoma. The 2019 tornadoes are especially prominent compared with other years. While the touchdown point map identifies where tornadoes began, the path map provides additional information about the distance and direction of tornado movement. Viewing these maps together demonstrates how different geometry types can provide complementary information about the same events.

Exercise 2: Annual Tornado Density by County

The second exercise examines how tornado density varies across Oklahoma counties from 2016 to 2021. Tornado touchdown points are spatially joined to county polygons and counted separately for each year. These counts are then divided by county area to calculate tornadoes per 1,000 square kilometers. The resulting maps use a consistent color scale to allow comparisons across years.

# Spatially join tornado points to counties
county_year_points <- st_join(tpoint_16_21, okcounty) %>%
  st_drop_geometry()

# Count tornadoes by county and year
county_year_counts <- county_year_points %>%
  filter(!is.na(GEOID)) %>%
  group_by(GEOID, yr) %>%
  summarize(tcnt = n(), .groups = "drop")

# Create all county-year combinations while preserving geometry
county_year_map <- okcounty %>%
  mutate(area = as.numeric(st_area(geometry))) %>%
  slice(rep(seq_len(n()), each = 6)) %>%
  mutate(yr = rep(2016:2021, times = nrow(okcounty))) %>%
  left_join(county_year_counts, by = c("GEOID", "yr")) %>%
  mutate(
    tcnt = replace_na(tcnt, 0),
    tdens = 10^9 * tcnt / area
  )

# Create faceted tornado density maps
exercise2_map <- ggplot(county_year_map) +
  geom_sf(aes(fill = tdens), color = "gray70", linewidth = 0.15) +
  scale_fill_distiller(
    name = expression("Tornadoes/1000 km"^2),
    palette = "YlOrRd",
    direction = 1
  ) +
  facet_wrap(vars(yr), ncol = 2) +
  coord_sf(datum = NA) +
  labs(title = "Tornado Density by County, 2016–2021") +
  theme_void() +
  theme(
    legend.position = "bottom",
    strip.text = element_text(size = 11, face = "bold"),
    plot.title = element_text(hjust = 0.5)
  )

exercise2_map

Figure 2. Annual tornado touchdown density across Oklahoma counties from 2016 to 2021, expressed as the number of tornadoes per 1,000 square kilometers. A consistent color scale allows comparisons between years.

The faceted maps show considerable variation in tornado density between years. Tornado activity appears particularly widespread in 2019, with several counties in central and northeastern Oklahoma exhibiting relatively high densities. In contrast, 2020 shows lower tornado densities across much of the state. Some counties have no recorded tornado touchdown points in particular years, illustrating the uneven spatial and temporal distribution of tornado activity. Using a consistent color scale across all six maps makes these annual differences easier to compare.

Exercise 3: Comparing Quantile Classification Schemes

The final exercise examines how the number of quantile classes affects the visualization of tornado density. County-level tornado densities are calculated for the full 2016–2021 period and classified into three, four, five, and six quantile classes. Four choropleth maps are then combined using plot_grid() to compare how different classification schemes influence the appearance and interpretation of spatial patterns.

# Count tornado touchdown points by county
county_points <- st_join(tpoint_16_21, okcounty) %>%
  st_drop_geometry()

countysum <- county_points %>%
  filter(!is.na(GEOID)) %>%
  group_by(GEOID) %>%
  summarize(tcnt = n(), .groups = "drop")

# Calculate tornado density per 1,000 square kilometers
countymap <- okcounty %>%
  left_join(countysum, by = "GEOID") %>%
  mutate(
    tcnt = replace_na(tcnt, 0),
    area = st_area(geometry),
    tdens = 10^9 * tcnt / area
  ) %>%
  drop_units()

# Function to create quantile classes
make_quantiles <- function(x, n) {
  breaks <- quantile(
    x,
    probs = seq(0, 1, length.out = n + 1),
    na.rm = TRUE
  )
  cut(x, breaks = breaks, include.lowest = TRUE)
}

# Classify tornado density into 3, 4, 5, and 6 classes
countymap <- countymap %>%
  mutate(
    q3 = make_quantiles(tdens, 3),
    q4 = make_quantiles(tdens, 4),
    q5 = make_quantiles(tdens, 5),
    q6 = make_quantiles(tdens, 6)
  )

# Function to create choropleth maps

make_choropleth <- function(data, column, title) {
  ggplot(data) +
    geom_sf(
      aes(fill = .data[[column]]),
      color = "gray70",
      linewidth = 0.15
    ) +
    scale_fill_brewer(
      palette = "YlOrRd",
      name = "Tornado Density",
      guide = guide_legend(ncol = 1, reverse = TRUE)
    ) +
    coord_sf(datum = NA) +
    labs(title = title) +
    theme_void() +
    theme(
      legend.position = "bottom",
      legend.text = element_text(size = 7),
      legend.title = element_text(size = 8),
      plot.title = element_text(
        hjust = 0.5,
        face = "bold",
        size = 11
      ),
      plot.background = element_rect(fill = "white", color = NA)
    )
}


# Create four maps
map_q3 <- make_choropleth(countymap, "q3", "3 Quantile Classes")
map_q4 <- make_choropleth(countymap, "q4", "4 Quantile Classes")
map_q5 <- make_choropleth(countymap, "q5", "5 Quantile Classes")
map_q6 <- make_choropleth(countymap, "q6", "6 Quantile Classes")

# Combine maps into a 2 x 2 grid
exercise3_map <- plot_grid(
  map_q3, map_q4, map_q5, map_q6,
  ncol = 2,
  labels = c("A", "B", "C", "D"),
  label_size = 12
)

exercise3_map

Figure 3. Comparison of county-level tornado density in Oklahoma for 2016–2021 using three (A), four (B), five (C), and six (D) quantile classes. Although the underlying data remain the same, increasing the number of classes changes how differences in tornado density are represented.

The four maps demonstrate how the number of quantile classes influences the interpretation of geographic patterns. With three classes, broader regional differences in tornado density are easier to identify, but variation between neighboring counties is less apparent. Increasing the number of classes reveals more local differences, particularly in central and northeastern Oklahoma. However, the six-class map is also more visually complex and requires greater attention to the legend. Although the underlying density values remain unchanged, the classification scheme affects how strongly certain counties stand out.

Conclusion

These exercises demonstrate how different spatial representations and cartographic choices influence the interpretation of tornado activity in Oklahoma. Point and line geometries provide complementary information about tornado locations and paths, while county-level density maps reveal differences across space and time. Comparing quantile classifications further illustrates how map design can emphasize or obscure geographic patterns. Together, these approaches highlight the importance of selecting appropriate spatial analysis and visualization methods for communicating geographic information.

Part B: Interactive Mapping with Leaflet

Introduction

For Part B, I used the leaflet package to create an interactive map of tornado touchdown locations in Oklahoma from 2016 to 2021. Building on the static maps from Part A, this map allows users to zoom into different areas and explore individual tornado events. Each touchdown point is colored according to its year, and interactive labels and popups provide additional information about each tornado.

# Load required packages
library(leaflet)

Map Creation

To create the interactive map, I used the tornado touchdown point dataset from Part A, filtered to events between 2016 and 2021. I first transformed the data to the WGS 84 coordinate system (EPSG:4326), which is compatible with Leaflet. I then created a categorical color palette using colorFactor() to assign a different color to each year.

Using the leaflet package, I added an Esri WorldTopoMap and plotted the tornado touchdown locations as colored circle markers. I also included a legend identifying the years, hover labels displaying tornado dates, and clickable popups containing the year, date, and tornado ID. These features allow users to explore individual events while maintaining a broader view of tornado activity across Oklahoma. Additionally, I overlaid Oklahoma county boundaries using the county shapefile from Part A. The boundaries are displayed as thin gray lines, with county names appearing when users hover over them, providing additional geographic context for the tornado touchdown locations.

# Prepare county boundaries for Leaflet
counties_leaflet <- okcounty %>%
  st_transform(4326)

# Prepare tornado points for Leaflet
tornado_leaflet <- tpoint_16_21 %>%
  st_transform(4326) %>%
  mutate(
    year_label = as.character(yr),
    popup_info = paste0(
      "<b>Tornado Touchdown</b>",
      "<br>Year: ", yr,
      "<br>Date: ", date,
      "<br>Tornado ID: ", om
    ),
    hover_label = paste("Tornado:", date)
  )

# Create a categorical color palette for each year
year_pal <- colorFactor(
  palette = "Dark2",
  domain = as.character(2016:2021)
)

# Create interactive map
tornado_map <- leaflet(tornado_leaflet,
                       width = "100%",
                       height = 550) %>%
  addProviderTiles(providers$Esri.WorldTopoMap) %>%
  addPolygons(
  data = counties_leaflet,
  fillColor = "transparent",
  fillOpacity = 0,
  color = "grey",
  weight = 0.5,
  opacity = 1,
  label = ~NAME
) %>%
  setView(lng = -97.3, lat = 35.5, zoom = 7) %>%
  addCircleMarkers(
    radius = 4,
    color = ~year_pal(year_label),
    fillColor = ~year_pal(year_label),
    fillOpacity = 0.7,
    weight = 1,
    popup = ~popup_info,
    label = ~hover_label
  ) %>%
  addLegend(
    position = "bottomright",
    pal = year_pal,
    values = ~year_label,
    title = "Tornado Year",
    opacity = 1
  )

tornado_map

Figure 4. Interactive Map of Tornado Touchdown Locations in Oklahoma (2016–2021). Tornado touchdown points are colored by year and displayed over Oklahoma county boundaries. Users can zoom and pan to explore spatial patterns, hover over points to view tornado dates, and click individual markers for additional event information, including year, date, and tornado ID

Results and Discussion

The interactive map shows tornado touchdown locations distributed across Oklahoma between 2016 and 2021, with multiple years represented throughout the state. The year-based color scheme makes it possible to identify areas where tornadoes occurred in different years and examine concentrations of events more closely.

Compared with the static maps in Part A, the interactive map provides greater flexibility for exploring individual tornadoes. Users can zoom into specific counties, examine the locations of nearby events, and access additional information without overcrowding the map with labels. While the static maps are better suited for comparing overall spatial patterns and tornado density, the interactive map is useful for investigating individual events and local geographic patterns.