Tornadoes in Oklahoma
GEOG588 Lab 4 Part A
Anna Betros
March 26, 2026
Introduction
Oklahoma averages more than 50 tornadoes a year and is one of the most tornado-prone areas in the world, according to the Oklahoma Emergency Management Department (2026). This report visualizes the initial touch down points of tornadoes and the paths that tornadoes follow on the ground in Oklahoma from 2016 to 2021. This report also visualizes the density of tornadoes by county and by year. The tornado data are derived from larger, national-level datasets generated by the National Oceanographic and Atmospheric Administration (NOAA) National Weather Service Storm Prediction Center (Wimberly, 2023; NOAA, 2025). Visualizing this data answers the three questions posed in Chapter 5 Vector Geospatial Data (Wimberly, 2023).
Data Preparation
Prepare the work space by loading numerous packages.
library(sf)
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
library(RColorBrewer)
library(units)
library(cowplot)Data preparation for this report is minimal. The tornado paths and initial points data are already limited to Oklahoma; the only preparation remaining is to filter the data from 2016 to 2021.
# Read in shapefiles and suppress information messages
okcounty <- st_read("ok_counties.shp", quiet = TRUE)
tpoint <- st_read("ok_tornado_point.shp", quiet = TRUE)
tpath <- st_read("ok_tornado_path.shp", quiet = TRUE)
# Filter tornado initial point data from 2016 to 2021 and assign to a new variable
tpoint_16_21 <- tpoint %>% filter(yr >= 2016 & yr <= 2021)
# Filter tornado path data from 2016 to 2021 and assign to a new variable
tpath_16_21 <- tpath %>% filter(yr >= 2016 & yr <= 2021) Question 1
Generate a map of tornado paths where the paths from each year are displayed as a different color, similar to the map of tornado points in Figure 5.4. Create a composite figure containing the map of tornado paths and the map of tornado points using plot_grid() (Wimberly, 2023).
Method: the two maps are created separately and are plotted together using plot_grid().
# Map the tornado paths by year overlaid on OK counties
torpath <- ggplot() +
# Map OK counties
geom_sf(data = okcounty, fill = "whitesmoke") +
# Map tornado paths and distinguish tornadoes by year with color
geom_sf(data = tpath_16_21,
aes(color = as.factor(yr)),
size = 1) +
# Apply symbology to the paths and title the legend
scale_color_discrete(name = "Year", palette = "YlOrRd") +
# Do not display graticules
coord_sf(datum = NA) +
# Remove extraneous map details
theme_void() +
# Provide a title and caption
labs(title = "Path of Tornadoes in Oklahoma, 2016 - 2021",
caption = "Figure 1: Tornado paths. This map illustrates the path of each tornado that
touched down in Oklahoma between 2016 and 2021. Data is obtained
from NOAA's National Weather Service Storm Prediction Center.") +
# change text size
theme(plot.caption = element_text(size = 8))
# Map the tornado start points by year overlaid on OK counties
torstart <- ggplot() +
# Map OK counties
geom_sf(data = okcounty, fill = "whitesmoke") +
# Map tornado start points and distinguish the year with color
geom_sf(data = tpoint_16_21,
aes(color = as.factor(yr))) +
# Apply symbology to the points and title the legend
scale_color_discrete(name = "Year", palette = "YlOrRd") +
# Remove extraneous information
coord_sf(datum = NA) +
theme_void() +
# Add map title and caption
labs(title = "Start Point of Tornadoes in Oklahoma, 2016 - 2021",
caption = "Figure 2: Tornado start points. This map illustrates the inital locations
of tornadoes that touched down in Oklahoma between 2016 and 2021. Data
is obtained from NOAA's National Weather Service Storm Prediction Center.") +
# change text size
theme(plot.caption = element_text(size = 8))
# Show two maps in one plot window
plot_grid(torpath, torstart,
ncol = 1, # specify number of columns
hjust = 0, # horizontal justify labels
label_x = 0,
align = "hv") # align maps horizontally and vertically
Results: Both figures 1 and 2 show a southwest-to-northeast pattern of tornadoes. Figure 1 shows that tornadoes appear to have longer paths in northeast Oklahoma than elsewhere in the state.
Question 2
Summarize the density of tornado points by both county and year and generate a faceted plot that displays maps of county-level tornado density from 2016-2021 (Wimberly, 2023).
Method: a for loop calculates the density of tornadoes per county by year and is captured in a new column (tdens). Each iteration of the loop runs on one year of data. The result of the first iteration is appended to an empty dataframe. The results of the following iterations are appended to the same dataframe. Six plots are shown together with a facet_wrap().
# Create empty dataframe
df <- data.frame()
# Create FOR loop to calculate tornado density by year and by county
# Iterate every year from 2016 to 2021
for (year in 2016:2021){
# filter point data by year (2016 - 2021) and assign to new variable
tpointyear <- tpoint_16_21 %>% filter(yr == year)
# Spatial-join OK county data and tornado point data, assign to new variable
countypnt <- st_join(tpointyear, okcounty)
# Drop geometry
countypnt <- st_drop_geometry(countypnt)
# Summarize the tornado count by county and save to variable (tcnt)
countysum <- countypnt %>%
group_by(GEOID) %>%
summarize(tcnt = n())
# Join datasets and calculate tornado density (tdens) by county
countymap <- okcounty %>%
left_join(countysum, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area,
yr = year) %>%
drop_units()
# Append yearly results to empty or prior dataframes
df <- rbind(df,countymap)
} # end loop
# Map county-density data
ggplot(data = df) +
# Identify tornado density value (tdens) to symbolize
geom_sf(aes(fill = tdens)) +
# Symbolize values and name the legend
scale_fill_distiller(name = expression("Tornadoes/1,000 km"^2),
palette = "YlOrRd",
trans = "reverse",
breaks = pretty_breaks()) +
# Create six maps each showing one year, and organize in two columns
facet_wrap(vars(yr), ncol = 2) +
# Relocate legend and change its size
theme_void() +
theme(legend.position = "top",
legend.title = element_text(size = 11),
legend.text = element_text(size = 9),
legend.key.size = unit(0.5, "cm"),
legend.direction = "horizontal",
plot.title = element_text(face = "bold")) +
# Add map title and caption
labs(title = "Yearly Tornado Density by OK County, 2016 to 2021",
subtitle = "",
caption = "Figure 3: Tornado densities by year. This map series illustrates the
densities of tornadoes in Oklahoma counties between 2016 and 2021. Data is
obtained from NOAA's National Weather Service Storm Prediction Center.")
Results: the yearly maps in figure 3 also show a southwest-to-northeast trend in tornadoes, specifically in years 2017 and 2021. The year 2020 shows a southeastern shift in the tornado trendline. The year 2019 has the high density of tornadoes, located in the northeastern counties, compared to the other years.
Question 3
Generate four choropleth maps of tornado density based on quantile breaks with numbers of classes ranging from 3 to 6. Create a composite figure containing the four maps using plot_grid() and examine how the number of classes changes the interpretation of the map (Wimberly, 2023).
Method: four maps showing the same data but using different class numbers are plotted together with plot_grid().
# Calculate tornado density by county for all years
# Spatial-join OK county data and tornado point data, assign to new variable
countypnt <- st_join(tpoint_16_21, okcounty)
# Drop geometry so that the data will be a normal dataframe
countypnt <- st_drop_geometry(countypnt)
# Group tornadoes by county (GEOID) and count the tornadoes and assign to new variable (tcnt)
countysum <- countypnt %>%
group_by(GEOID) %>%
summarize(tcnt = n())
# Join county tornado sum to OK county
countymap <- okcounty %>%
# join field is county ID (GEOID)
left_join(countysum, by = "GEOID") %>%
# replace NA values with 0 because some counties did not have tornadoes
replace(is.na(.), 0) %>%
# Calculate density by...
# Calculating county area
mutate(area = st_area(okcounty),
# tornado count / area -> convert to km2 -> assign to new variable
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
# Create four choropleth maps by county and tornado density
# Create four variables with different class numbers
numclas3 <- 3
numclas4 <- 4
numclas5 <- 5
numclas6 <- 6
# Define quantiles using classes defined in (numclas#)
qbrks3 <- seq(0, 1, length.out = numclas3 + 1)
qbrks4 <- seq(0, 1, length.out = numclas4 + 1)
qbrks5 <- seq(0, 1, length.out = numclas5 + 1)
qbrks6 <- seq(0, 1, length.out = numclas6 + 1)
# Create four new columns (tdens_c#) for each quantile (qbrks#) by cutting
# the tornado density column (tdens).
countymap <- countymap %>%
mutate(tdens_c3 = cut(tdens, breaks = quantile(tdens, qbrks3), include.lowest = T)) %>%
mutate(tdens_c4 = cut(tdens, breaks = quantile(tdens, qbrks4), include.lowest = T)) %>%
mutate(tdens_c5 = cut(tdens, breaks = quantile(tdens, qbrks5), include.lowest = T)) %>%
mutate(tdens_c6 = cut(tdens, breaks = quantile(tdens, qbrks6), include.lowest = T))
# Create four choropleth maps with each map plotting a column with different class
# numbers (tdens_c#)
map3 <- ggplot(data = countymap) +
geom_sf(aes(fill = tdens_c3)) +
scale_fill_brewer(name = expression("Tornadoes/1,000 km"^2),
palette = "YlOrRd") +
theme_void() +
theme(legend.position = "right",
legend.title = element_text(size =9),
legend.text = element_text(size =7),
legend.key.size = unit(0.4, "cm"),
legend.direction = "vertical",
plot.title = element_text(face = "bold")) +
labs(title = "Tornado Densities by County, Oklahoma 2016 - 2021",
subtitle = "Three classes")
map4 <- ggplot(data = countymap) +
geom_sf(aes(fill = tdens_c4)) +
scale_fill_brewer(name = expression("Tornadoes/1,000 km"^2),
palette = "YlOrRd") +
theme_void() +
theme(legend.position = "right",
legend.title = element_text(size =9),
legend.text = element_text(size =7),
legend.key.size = unit(0.4, "cm"),
legend.direction = "vertical") +
labs(subtitle = "Four classes")
map5 <- ggplot(data = countymap) +
geom_sf(aes(fill = tdens_c5)) +
scale_fill_brewer(name = expression("Tornadoes/1,000 km"^2),
palette = "YlOrRd") +
theme_void() +
theme(legend.position = "right",
legend.title = element_text(size =9),
legend.text = element_text(size =7),
legend.key.size = unit(0.4, "cm"),
legend.direction = "vertical") +
labs(subtitle = "Five classes")
map6 <- ggplot(data = countymap) +
geom_sf(aes(fill = tdens_c6)) +
scale_fill_brewer(name = expression("Tornadoes/1,000 km"^2),
palette = "YlOrRd") +
theme_void() +
theme(legend.position = "right",
legend.title = element_text(size =9),
legend.text = element_text(size =7),
legend.key.size = unit(0.4, "cm"),
legend.direction = "vertical",
plot.caption = element_text(size = 7)) +
labs(subtitle = "Six classes",
caption = "Figure 4: Tornado densities. This map series illustrates the densities of tornadoes in Oklahoma
counties between 2016 and 2021 as choropleth maps with different class numbers. These maps
illustrate how using a different number of classes can affect the interpretation of the counties'
tornado densities. Data is obtained from NOAA's National Weather Service Storm Prediction Center.")
# Show the four maps in one plot window
plot_grid(map3, map4, map5, map6,
ncol = 2, # specify number of columns
align = "hv") # align maps horizontally and vertically
Results: figure 4 also shows the southwest-to-northeastern trend of tornadoes despite the four different visualization methods. The six-classed map (lower right) highlights the counties with the lowest and highest values whereas the three-classed map (upper left) shows averaged values that tend to hide the extremes.
Sources
National Weather Service Storm Prediction Center. (2025). Tornado Paths and Initial Points. [Data]. National Oceanographic and Atmospheric Administration (NOAA). Retrieved March 25, 2026 from https://www.spc.noaa.gov/gis/svrgis/.
Oklahoma Emergency Management Department. (2026). Tornado Awareness. Retrieved March 26, 2026 from https://oklahoma.gov/oem/readyok/be-aware/tornadoes.html.
Wimberly, M. 2023. Geographic Data Science. [Online edition]. Retrieved March 25, 2026 from https://bookdown.org/mcwimberly/gdswr-book/vector-geospatial-data.html.