Tornado Paths in Oklahoma

R packages used to prepare data and to create these maps include sf, rgdal, ggplot2, dplyr, scales, RColorBrewer, units, and cowplot.

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

# Load files------
okcounty <- st_read("Chapter5/ok_counties.shp", quiet = TRUE)
tpoint <- st_read("Chapter5/ok_tornado_point.shp", quiet = TRUE)
tpath <- st_read("Chapter5/ok_tornado_path.shp", quiet = TRUE)


# As sp
okcounty_sp <- as(okcounty, 'Spatial')

okcounty_sf <- st_as_sf(okcounty_sp)

# Create subset of tornado points between 2016 and 2021
tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

# Create subset of tornado paths between 2016 and 2021
tpath_16_21 <- tpath %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

# Join points to OK counties
countypnt <- st_join(tpoint_16_21, okcounty)

# Drop geometry and summarize count by county
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>% 
  group_by(GEOID) %>%
  summarize(tcnt = n())


# Join county sum
countymap <- okcounty %>%
  left_join(countysum, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()


okcntrd <- st_centroid(countymap)

Tornado Paths and Points

Tornadoes in Oklahoma are common, and there were 434 tornadoes in Oklahoma between 2016 and 2021. The median number of tornadoes in Oklahoma during this time was 60 tornadoes per year.

# create a map of tornado paths
tpathcolored <- ggplot() +
  geom_sf(data = tpath_16_21, 
          aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) +
    scale_color_manual(name = "Year", 
                      values = c("#c7e9b4","#41b6c4","#1d91c0","#225ea8","#253494", "#081d58")) +
  labs(title = "a. Tornado Paths") +
  coord_sf(datum = NA) +
  theme_void()

# create a map of tornado points
tpointcolored <- ggplot() +
  geom_sf(data = tpoint_16_21, 
          aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_manual(name = "Year", 
                      values = c("#c7e9b4","#41b6c4","#1d91c0","#225ea8","#253494", "#081d58")) +
  labs(title = "b. Tornado Points") +
  coord_sf(datum = NA) +
  theme_void()

# Use plot_grid to combine paths and points

plot_grid(tpathcolored, tpointcolored,
          ncol = 1,
          hjust = 0,
          label_x = 0,
          align = "hv")

Figure 1. Tornadoes in Oklahoma between 2016 and 2021 are displayed as paths (left) and points (right). Each year is distinguished using an ordinal scale.

Density of Tornadoes by County and Year

The density of tornadoes varies by county within Oklahoma. Some counties are more likely to be hit by tornadoes, such as Tulsa county. Others had few or tornadoes during this time period, such as Grant and Alfalfa. Tornadoes also vary by year; using the data from this dataset, there were 148 tornadoes in 2019, but only 38 tornadoes in 2020.

countysum2 <- countypnt %>%
  group_by(GEOID, yr) %>%
  summarize(tcnt = n())

# Calculate density and add area column
countymap2 <- okcounty %>%
  left_join(countysum2, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(.),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

countymap2 <- countymap2 %>% filter(yr != 0)


# Create a map of tornadoes by year
ggplot(data = countymap2) +
  geom_sf(data = okcounty, fill = NA) +
  geom_sf(aes(fill = tdens)) +
  scale_fill_distiller(name = expression("Tornadoes/1000 km"^2),
                       palette = "YlOrRd",
                       breaks = pretty_breaks(),
                       direction = 1) +
  labs(title = "Tornado Density by Oklahoma County",
       subtitle = "Between 2016 and 2021\n") +
  theme_void() +
  theme(legend.position = "bottom") +
  facet_wrap(~yr)

Figure 2. Tornado density by county between 2016 and 2021 in Oklahoma. Tornadoes were most frequent in 2019 and least frequent in 2020.

Tornado Density - Different Classes

First, I set the number of classes and breaks for each map. Next, I create a new density variable using the breaks for each class.

# Question 3 -----

# Set number of classes an breaks for each map
numclas_3 <- 3
qbrks_3 <- seq(0, 1, length.out = numclas_3 + 1)

numclas_4 <- 4
qbrks_4 <- seq(0, 1, length.out = numclas_4 + 1)

numclas_5 <- 5
qbrks_5 <- seq(0, 1, length.out = numclas_5 + 1)

numclas_6 <- 6
qbrks_6 <- seq(0, 1, length.out = numclas_6 + 1)

# Create new density variable using the breaks
countymap_quant <- countymap %>%
  mutate(tdens_c3 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks_3,
                                          include.lowest = T)),
         tdens_c4 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks_4,
                                          include.lowest = T)),
         tdens_c5 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks_5,
                                          include.lowest = T)),
         tdens_c6 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks_6,
                                          include.lowest = T)))

After the breaks and classes have been set, I create four maps within the same figure. These maps depict the difference between density set at four different breakpoints. The maps contain different numbers of classes - each map is different, but depicts three, four, five, and six classes.

# A  map featuring 3 classes       
q_3_map <- ggplot(data = countymap_quant) +
  geom_sf(aes(fill = tdens_c3)) +
  scale_fill_brewer(name = "", palette = "YlOrRd") +
  labs(title = "a. Three Classes",
       subtitle = "Tornadoes/1000 km_2") +
  theme_void() +
  theme(legend.position = "bottom",
                 text = element_text(color = "#525252"))

# A map featuring 4 classes  
q_4_map <- ggplot(data = countymap_quant) +
  geom_sf(aes(fill = tdens_c4)) +
  scale_fill_brewer(name = "", palette = "YlOrRd") +
  labs(title = "b. Four Classes",
       subtitle = "Tornadoes/1000 km_2") +
  theme_void() +
  theme(legend.position = "bottom",
        text = element_text(color = "#525252"))

# A map featuring 5 classes 
q_5_map <- ggplot(data = countymap_quant) +
  geom_sf(aes(fill = tdens_c5)) +
  scale_fill_brewer(name = "", palette = "YlOrRd") +
  labs(title = "c. Five Classes",
       subtitle = "Tornadoes/1000 km_2") +
  theme_void() +
  theme(legend.position = "bottom",
        text = element_text(color = "#525252"))

# A map featuring 6 classes 
q_6_map <- ggplot(data = countymap_quant) +
  geom_sf(aes(fill = tdens_c6)) +
  scale_fill_brewer(name = "", palette = "YlOrRd") +
  labs(title = "d. Six Classes",
       subtitle = "Tornadoes/1000 km_2") +
  theme_void() +
  theme(legend.position = "bottom",
        text = element_text(color = "#525252"))

plot_grid(q_3_map, q_4_map, q_5_map, q_6_map,
          labels = "\nTornado Density by County",
          vjust = -.5,
          # labels = c("a. Three Classes",
          #            "b. Four Classes",
          #            "c. Five Classes",
          #            "d. Six Classes",
          #            label_size = 10),
          ncol = 2,
          align = "hv")

Figure 3. Tornado density by county between 2016 and 2021 using different breaks. Each map depicts a different number of classes: map A depicts 3 breaks, map B depicts 4 breaks, map C depicts 5 breaks, and map D depicts 6 breaks. An increase in breaks better reveals the highest tornado counties in Oklahoma.