Introduction: This project focuses on visualizing and analyzing tornado activity in Oklahoma from 2016 to 2021 using spatial data and various visualization techniques. The primary objectives of this project are to map tornado paths, points, and densities across Oklahoma counties to identify areas most affected by tornadoes.

Data

Part 1 Generate a map of tornado paths where paths from each year are displayed in different colors & Create the map of tornado points for reference

knitr::opts_chunk$set(echo = TRUE)
library(sf)
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
library(RColorBrewer)
library(units)
library(cowplot)

Part 1 Tornado Paths and Points

folder_path <- "C:\\PENNSTATE\\GEOG588_Analytical Approaches_\\R_Projects\\Lab4\\Data"

# Read the shapefiles using full file paths
okcounty <- st_read(file.path(folder_path, "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(file.path(folder_path, "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(file.path(folder_path, "ok_tornado_path.shp"), quiet = TRUE)

# Filter and create the necessary tpoint_16_21 dataset
tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  mutate(year = as.factor(yr))

# Filter and create the necessary tpath_16_21 dataset
tpath_16_21 <- tpath %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  mutate(year = as.factor(yr))

# Generate a map of tornado paths where paths from each year are displayed in different colors
tornado_paths <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "gray") +
  geom_sf(data = tpath_16_21, aes(color = year)) +
  scale_color_brewer(palette = "Set3", name = "Year") +
  theme_void() +
  ggtitle("A) Oklahoma Tornado Paths ") +
  theme(plot.title = element_text(hjust = 0.5, margin = margin(b = 10)))

# Create the map of tornado points 
tornado_points <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "gray") +
  geom_sf(data = tpoint_16_21, aes(color = year)) +
  scale_color_brewer(palette = "Set3", name = "Year") +
  theme_void() +
  ggtitle("B) Tornado Points In Oklahoma") +
  theme(plot.title = element_text(hjust = 0.5, margin = margin(b = 10)))

# Create a composite figure using plot_grid()
composite_figure <- plot_grid(tornado_paths, tornado_points, 
                              labels = NULL, ncol = 1)

# Display the composite figure
print(composite_figure)

Figure 1: Tornado activity in Oklahoma from 2016 to 2021. Panel A shows the tornado paths, color-coded by year, while Panel B displays the tornado points within the same timeframe.

Part 2 density of tornado points by both county and year

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

# Specify the folder path
folder_path <- "C:\\PENNSTATE\\GEOG588_Analytical Approaches_\\R_Projects\\Lab4\\Data"

# Read the shapefiles using full file paths
okcounty <- st_read(file.path(folder_path, "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(file.path(folder_path, "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(file.path(folder_path, "ok_tornado_path.shp"), quiet = TRUE)

# Filter and create the necessary tpoint_16_21 dataset
tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  mutate(year = as.factor(yr))

# Filter and create the necessary tpath_16_21 dataset
tpath_16_21 <- tpath %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  mutate(year = as.factor(yr))

# Summarize the density of tornado points by county and year
countypnt_year <- st_join(tpoint_16_21, okcounty) %>%
  st_drop_geometry() %>%
  group_by(GEOID, year) %>%
  summarize(tcnt = n())

# Merge the summarized data with the county shapefile
countymap <- okcounty %>%
  left_join(countypnt_year, by = "GEOID")

# Convert area to numeric values and remove NA values
countymap <- countymap %>%
  mutate(area_km2 = as.numeric(st_area(geometry) / 1e6),
         tdens = tcnt / area_km2 * 1000) %>%
  filter(!is.na(tdens))

# Generate a faceted plot displaying county-level tornado density for each year
tornado_density <- ggplot(data = countymap) +
  geom_sf(aes(fill = tdens)) +
  scale_fill_distiller(name = expression("Tornadoes/1000 km"^2), 
                       palette = "YlOrRd", 
                       breaks = pretty_breaks(), 
                       direction = 1) +
  facet_wrap(vars(year), ncol = 2) +
  ggtitle("County Level Tornado Density in Oklahoma 2016-2021") +
  theme_void() +
  theme(legend.position = "bottom",
        plot.title = element_text(hjust = 0.5, size = 14, margin = margin(b = 10)))

print(tornado_density)

Figure 2: County-Level Tornado Density in Oklahoma (2016-2021): This map visualizes the density of tornado occurrences per square kilometer across different counties in Oklahoma over the years 2016 to 2021.

Part 3 Four choropleth maps of tornado density based on quantile breaks with numbers of classes ranging from 3 to 6

knitr::opts_chunk$set(echo = TRUE)
library(sf)
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
library(RColorBrewer)
library(units)
library(cowplot)

folder_path <- "C:\\PENNSTATE\\GEOG588_Analytical Approaches_\\R_Projects\\Lab4\\Data"
okcounty <- st_read(file.path(folder_path, "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(file.path(folder_path, "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(file.path(folder_path, "ok_tornado_path.shp"), quiet = TRUE)

tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  mutate(year = as.factor(yr))

tpath_16_21 <- tpath %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  mutate(year = as.factor(yr))

countypnt_year <- st_join(tpoint_16_21, okcounty) %>%
  st_drop_geometry() %>%
  group_by(GEOID, year) %>%
  summarize(tcnt = n())

countymap <- okcounty %>%
  left_join(countypnt_year, by = "GEOID")

countymap <- countymap %>%
  mutate(area_km2 = as.numeric(st_area(geometry) / 1e6),
         tdens = tcnt / area_km2 * 1000) %>%
  filter(!is.na(tdens))

# Create choropleth maps with different numbers of quantile breaks
create_choropleth <- function(count, data) {
  breaks <- quantile(data$tdens, probs = seq(0, 1, length.out = count + 1), include.lowest = TRUE)
  data <- data %>%
    mutate(tdens_class = cut(tdens, breaks = breaks, include.lowest = TRUE))

  plot <- ggplot(data = data) +
    geom_sf(aes(fill = tdens_class)) +
    scale_fill_brewer(palette = "YlOrRd", name = expression("Tornadoes/1000 km"^2), na.value = "transparent") +
    theme_void() +
    theme(legend.position = "bottom") +
    labs(title = NULL)  # Remove subtitle
  
  list(plot = plot, breaks = breaks)
}

# Create the four choropleth maps and get the breaks
choropleth_3 <- create_choropleth(3, countymap)
choropleth_4 <- create_choropleth(4, countymap)
choropleth_5 <- create_choropleth(5, countymap)
choropleth_6 <- create_choropleth(6, countymap)

# Function to create a text annotation for the breaks
create_text_annotation <- function(breaks) {
  annotation <- paste0("Classes: ", paste0(round(breaks, 2), collapse = ", "))
  ggplot() + annotate("text", x = 0.5, y = 0.5, label = annotation, size = 4, hjust = 0.5) + theme_void()
}

# Create text annotations for each set of breaks
text_3 <- create_text_annotation(choropleth_3$breaks)
text_4 <- create_text_annotation(choropleth_4$breaks)
text_5 <- create_text_annotation(choropleth_5$breaks)
text_6 <- create_text_annotation(choropleth_6$breaks)

# Extract the legend from one of the plots
combined_legend <- cowplot::get_legend(choropleth_3$plot)

# Create the composite figure without legends and with text annotations
plot_3 <- plot_grid(choropleth_3$plot + theme(legend.position = "none"), text_3, ncol = 1, rel_heights = c(0.9, 0.1))
plot_4 <- plot_grid(choropleth_4$plot + theme(legend.position = "none"), text_4, ncol = 1, rel_heights = c(0.9, 0.1))
plot_5 <- plot_grid(choropleth_5$plot + theme(legend.position = "none"), text_5, ncol = 1, rel_heights = c(0.9, 0.1))
plot_6 <- plot_grid(choropleth_6$plot + theme(legend.position = "none"), text_6, ncol = 1, rel_heights = c(0.9, 0.1))

composite_figure <- plot_grid(
  plot_3, plot_4, plot_5, plot_6,
  labels = c("A) 3 Quantile Classes", "B) 4 Quantile Classes", "C) 5 Quantile Classes", "D) 6 Quantile Classes"),
  label_size = 12, ncol = 2
)

# Combine the composite figure with the combined legend
final_figure <- plot_grid(composite_figure, combined_legend, ncol = 1, rel_heights = c(1, 0.1))

# Display the final figure
print(final_figure)

Figure 3: Distribution of Tornado Density in Oklahoma showing the distribution oftornado density.