Analysis of tornaodes in Oklahoma, 2016 - 2021

Part A: 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().

load necessary R libraries

library(sf)
library(sp)
library(here)
library(tidyverse)
library(units)
library(scales)
library(RColorBrewer)
library(cowplot)

load in tornado and county data

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 data to only include storms between 2016 and 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)

create maps showing tornado touchdown locations and paths, colored by year

tor_paths <- ggplot() +
  geom_sf(data = tpath_16_21,
          aes(color = as.factor(yr))) + # set year as factor instead of int
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

tor_points <- 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) +
  theme_void()

use plot_grid to plot both points and paths in a similar format

plot_grid(tor_paths, tor_points,
          labels = c("A) Tornado paths",
                     "B) Tornado touchdowns",
                     label_size = 12),
          ncol = 1,
          hjust = 0,
          label_x = 0,
          align = "hv")

Part B: 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.

perform a spatial join between tornado data and county boundaries

tpoint_16_21 <- tpoint_16_21 %>%
  st_join(okcounty) %>%
  st_drop_geometry()

calculate area of each county in square kilometers

okcounty$AREA_KM2 <- st_area(okcounty) / 10^6

aggregate tornado data by county and year

tdens <- tpoint_16_21 %>%
  group_by(yr, GEOID) %>%
  summarise(tcnt = n()) %>%
  ungroup()

merge torndado data with county boundaries

okcounty <- left_join(okcounty, tdens, by = "GEOID")

calculate tornado density per 1,000 square kilometers

okcounty$tdens <- okcounty$tcnt / okcounty$AREA_KM2 * 1000

map the data by year using facet_wrap (I can’t get the density to plot, only the count)

ggplot() +
  geom_sf(data = na.omit(okcounty), aes(fill = tcnt)) +
  scale_fill_gradient(name = "Tornado Count", low = "yellow", high = "red") +
  facet_wrap(~yr, nrow = 2) +
  theme_void() +
  labs(title = "Tornado Count by County for Each Year (2016-2021)")

Part C: 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.

generate map for each of the four different break structures

class_3 <- ggplot() +
  geom_sf(data = na.omit(okcounty), aes(fill = tcnt)) +
  scale_fill_gradient(name = "Legend", breaks = quantile(okcounty$tcnt, probs = c(0, 0.33, 0.66, 1), na.rm = TRUE), labels = c("0%", "33%", "66%", "100%")) +
                        theme_void()

class_4 <- ggplot() +
  geom_sf(data = na.omit(okcounty), aes(fill = tcnt)) +
  scale_fill_gradient(name = "Legend", breaks = quantile(okcounty$tcnt, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE), labels = c("0%", "25%", "50%", "75%","100%")) +
                        theme_void()

class_5 <- ggplot() +
  geom_sf(data = na.omit(okcounty), aes(fill = tcnt)) +
  scale_fill_gradient(name = "Legend", breaks = quantile(okcounty$tcnt, probs = c(0, 0.2, 0.4, 0.6, 0.8, 1), na.rm = TRUE), labels = c("0%", "20%", "40%", "60%","80%", "100%")) +
                        theme_void()

class_6 <- ggplot() +
  geom_sf(data = na.omit(okcounty), aes(fill = tcnt)) +
  scale_fill_gradient(name = "Legend", breaks = quantile(okcounty$tcnt, probs = c(0, 0.16, 0.33, 0.50, 0.66, 0.83, 1), na.rm = TRUE), labels = c("0%", "16%", "33%", "50%","66%", "83%", "100%")) +
                        theme_void()

use plot grid to plot all four maps in a 2x2 grid

plot_grid(class_3, class_4, class_5, class_6,
          labels = c("A) 3 classes",
                     "B) 4 classes",
                     "C) 5 classes",
                     "D) 6 classes",
                     label_size = 12),
          ncol = 2,
          hjust = 0,
          label_x = 0,
          align = "hv")