Lab 4 - Tornadoes in Oklahoma (2016-2021)
Matthew Pawlenko
2023-06-08
Loading Packages and Reading in Data
# Load Packages ----
library(sf)
library(ggplot2)
library(dplyr)
library(tidyr)
library(RColorBrewer)
library(cowplot)
library(units)
library(here)
# Load 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 - 2016-2021
After the data is loaded, the tornado path and points data are filtered to onyl include 2016 through 2021.
# Filter by year
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)Plot Tornado Paths and Points
Plot of Just Paths
First, we plot just the tornado paths. The line width is increased to make the paths more visible and the pahts are colored by year. The plot is then saved to an object.
# Plotting Tornado Paths and Points ----
# Plot the path data, paths colored by year
paths <- ggplot() +
geom_sf(data = okcounty, fill = NA) +
geom_sf(data = tpath_16_21,
linewidth = .75, # increase width of tracks
aes(color = as.factor(yr))) +
scale_color_discrete(name = "Year") +
theme_void()Plot of Just Points
Next, we plot the tornado points. The size of each point is increased and the points are again colored by year. The plot is also saved to an object.
# Plot Point data, points colored by year
points <- ggplot() +
geom_sf(data = tpoint_16_21,
aes(color = as.factor(yr)),
size = 1) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
coord_sf(datum = NA) +
theme_void()Plot Both Paths and Points Maps Together
The plots of both tornado paths and points are plotted together using “plot_grid” from the cowplot package. The label_x and label_y are both adjusted to move the plot labels into more readable location.
# Plot both charts together
plot_grid(paths, points,
labels = c("A) Tornado Paths by Year",
"B) Tornado Locations by Year",
label_size = 12),
ncol = 1,
hjust = 0,
label_x = .15,
label_y = 1.03,
align = "hv")
Plotting Tornado Density by Year
Prepping the Data
Prior to plotting tornado density per year, per county must be calculated. To do this, the tornado point data is joined with the OK county data and the geometry is removed. The data is then grouped first by year, then by county ID (GEOID) in order to account for tornados by year.
# Plotting Tornado Density ----
# Calculate density per county with an st_join
# then summarizing by county
countypnt <- st_join(tpoint_16_21, okcounty)
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
group_by(yr, GEOID) %>%
summarize(tcnt = n())Calculating Tornado Density
The density is then calculated for each county based on the county area. To do this, the county area is first calculated and then the county data is joined to the summary of tornadoes per county. If the area is calculated after the join, there would be an error do to a mismatch between the amount of records in okcounty and those in countysum. Because two counties had “0” tornadoes during this time, both year and count were “NA.” Sicne setting all “NA” to zero would lead to a zero for a year, the NA values were dropped.
# Join tornado count per county and calculate tornado density
countymap <- okcounty %>%
mutate(area = st_area(okcounty)) %>% # Calculate area first
left_join(countysum, by = "GEOID") %>%
na.omit() %>% # Dropping NA values - counties with no tornadoes 2016-2021
mutate(tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()Plotting Tornado Density by Year
Finally, the density is plotted for each county and facet wrapped based on year.
# Plot Density per county by year
ggplot() +
geom_sf(data = okcounty) +
geom_sf(data = countymap, aes(fill = tdens)) +
scale_fill_distiller(name = expression("Tornadoes/1000 km"^2),
palette = "YlOrRd",
direction = 1) +
facet_wrap(~ yr) +
theme_void() +
ggtitle("Tornado Density of OK Counties from 2016-2021") +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5, vjust = 5))
Tornado Density Shown in Choropleth Maps
Prepping the Data
First, the tornado density needs to be recalculated for the entire 5 year span vice broken up by year.
# Create 4 choropleth maps ----
# Calculate density per county with an st_join
# then summarizing by county
countypnt <- st_join(tpoint_16_21, okcounty)
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
group_by(GEOID) %>%
summarize(tcnt = n())
# Join tornado count per county and calculate tornado density
countymap <- okcounty %>%
mutate(area = st_area(okcounty)) %>% # Calculate area first
left_join(countysum, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()Next, the tornado density is then broken into different quantiles - 3, 4, 5, and 6. This is done by “cutting” the tornado density vasriable based on user defined quantiles. The breaks for the quantiles are set using the “probs” attribute in quantile.
# Quantile Breaks = 3
numclas <- c(3, 4, 5, 6)
# Creating Density Variables based on different quantiles
countymap_quant <- countymap %>%
mutate(tdens_c3 = cut(tdens,
breaks = quantile(tdens,
probs = seq(0, 1, length.out = 4)),
include.lowest = T)) %>%
mutate(tdens_c4 = cut(tdens,
breaks = quantile(tdens,
probs = seq(0, 1, length.out = 5)),
include.lowest = T)) %>%
mutate(tdens_c5 = cut(tdens,
breaks = quantile(tdens,
probs = seq(0, 1, length.out = 6)),
include.lowest = T)) %>%
mutate(tdens_c6 = cut(tdens,
breaks = quantile(tdens,
probs = seq(0, 1, length.out = 7)),
include.lowest = T))Plotting the Choropleth Maps
Each choropleth map is then created and saved as an object. Once all 4 are created, they are plotted together. Unfortunately, due to layout, the titles and legends overlap and are unreadable.
# Plotting Choropleth map with 3 quantiles using "probs"
choropleth3 <- ggplot(data = countymap_quant) +
geom_sf(aes(fill = tdens_c3)) +
scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),
palette = "PuBuGn") +
theme_void() +
theme(legend.position = "bottom")
# Plotting Choropleth map with 4 quantiles using "probs"
choropleth4 <- ggplot(data = countymap_quant) +
geom_sf(aes(fill = tdens_c4)) +
scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),
palette = "PuBuGn") +
theme_void() +
theme(legend.position = "bottom")
# Plotting Choropleth map with 5 quantiles using "probs"
choropleth5 <- ggplot(data = countymap_quant) +
geom_sf(aes(fill = tdens_c5)) +
scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),
palette = "PuBuGn") +
theme_void() +
theme(legend.position = "bottom")
# Plotting Choropleth map with 6 quantiles using "probs"
choropleth6 <- ggplot(data = countymap_quant) +
geom_sf(aes(fill = tdens_c6)) +
scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),
palette = "PuBuGn") +
theme_void() +
theme(legend.position = "bottom")
# Plotting all 4 choropleth maps on one chart
plot_grid(choropleth3, choropleth4, choropleth5, choropleth6,
labels = c("A) 3 Quantiles (2016-2021)",
"B) 4 Quantiles (2016-2021)",
"C) 5 Quantiles (2016-2021)",
"D) 6 Quantiles (2016-2021)",
label_size = 1),
ncol = 2,
hjust = 0,
label_x = 0.05,
label_y = 1.025,
align = "hv")