Introduction

This R markdown document contains several map plots to visualize the incidence of tornados in Oklahoma by location and across time (years 2016-2021). Practice Exercise 1 shows the tornado data by both points and paths. Then Practice Exercises 2 and 3 use the point data to calculate densities in each county and display them continuously by year (Exercise 2) or by different quantiles during all six years (Exercise 3).

Data Preparation

This chunk imports the three shapefiles: 1) Oklahoma counties, 2) tornado points, and 3) tornado paths. Then it filters (2) and (3) so that we use tornados that occurred in years 2016 through 2021 in the analysis below.

# Packages needed for this report
library(tidyverse)
library(sf)
library(here)
library(cowplot)
library(units)

# Import three datasets
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 the data to only years 2016-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)

Practice Exercise 1

This chunk produces a composite figure of tornado points and paths, symbolized by years 2016-2021.

# Generate a map of tornado paths where the paths from each year are displayed as a different color
ok_t_paths_16_21 <- ggplot() +
  geom_sf(data = tpath_16_21, 
          aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void() +
  labs(title = "Oklahoma Tornado Paths by Year, 2016-2021")

#ggsave(here("outputs","ok_t_paths_16_21.png"))

# Generate a map of tornado points where each year has as a different color
ok_t_points_16_21 <- 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()+
  labs(title = "Oklahoma Tornado Points by Year, 2016-2021")

#ggsave(here("outputs","ok_t_points_16_21.png"))

# Create a composite figure containing the map of tornado paths and the map of tornado points using plot_grid()
plot_grid(ok_t_points_16_21, ok_t_paths_16_21, 
          ncol = 1, 
          hjust = 0, 
          label_x = 0, 
          align = "hv")

Practice Exercise 2

This chunk produces a cloropleth map of tornado point density by county, faceted by years 2016-2021. Of the six years, 2019 clearly had more tornado activity, which was concentrated in the NE corner of the state. Odd numbered years had more tornado activity in the western half of the state.

# Summarize the number of tornado points in each county
countyptsum <- st_join(tpoint_16_21, okcounty) %>% #spatial join of tornado points to the county polygons
  st_drop_geometry() %>% #removes geometry column for tabular manipulation
  group_by(GEOID, yr) %>%
  summarize(t_count = n())

# Join the resulting summary counts to the county polygon layer 
countyptsum_map <- okcounty %>%
  left_join(countyptsum, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(.),
         tdens = 10^6 * 10^3 * t_count / area) %>%
  drop_units() %>%
  filter(yr > 0) #filter out 2 records with yr==0 (dk where they came from)

# Map the tornado point density in each county by year - faceted plot
ggplot(data = countyptsum_map) +
  geom_sf(data = okcounty, fill = NA) + #first plot blank counties
  geom_sf(aes(fill = tdens)) + #then layer on the densities
  facet_wrap(vars(yr), ncol = 2) +
  theme_void() +
  labs(title = "Oklahoma Tornado Point Density by County, for Years 2016-2021")

Practice Exercise 3

This chunk uses a for loop to create four choropleth maps of tornado density with four different quantile breaks: 3, 4, 5, and 6. Then the plot_grid() function creates a composite figure. The quintile (breaks=5) map appears to have more counties with higher tornado density; however, this perception is also affected by which colors are used for which density value ranges in each map.

# Use for loop to generate the 4 quantile dataframes needed to create the maps
for(i in (3:6)) {
  numclas <- i
  qbrks <- seq(0, 1, length.out = numclas + 1)
  
  # Calculate quantile values
  countyptsum_temp <- countyptsum_map %>%
    mutate(tdens_c1 = cut(tdens,
                          breaks = quantile(tdens, probs = qbrks),
                          include.lowest = T))
  
  # Plot the quantile values on map
  countydensity_temp <- ggplot() +
    geom_sf(data = okcounty, fill = NA) +
    geom_sf(data = countyptsum_temp, 
            aes(fill = tdens_c1)) +
    scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                      palette = "YlOrRd") +
    coord_sf(datum = NA) +
    theme_void() +
    theme(legend.position = "bottom")

  # Save the quantile values to a dataframe for potential future use
  assign(paste0('countyptsum_map', i), countyptsum_temp)
  
  # Save the map plot to an object for plot_grid to use
  assign(paste0('countydensity_map', i), countydensity_temp)
  
  # Cleanup temp objects
  remove(countyptsum_temp)
  remove(countydensity_temp)
}

# Use plot_grid() function to create composite figure 
plot4 <- plot_grid(countydensity_map3, 
          countydensity_map4, 
          countydensity_map5, 
          countydensity_map6,
          labels = c("3 Quantiles", 
                     "4 Quantiles", 
                     "5 Quantiles", 
                     "6 Quantiles"
                     ),
          ncol = 1, 
          hjust = 0, 
          label_x = 0, 
          align = "hv")

# Add title to composite figure
plot4title <- ggdraw() + 
  draw_label(
    "Different Quantiles of OK Tornado Pt Density by County, during 2016-2021",
    fontface = 'bold',
    x = 0,
    hjust = 0
  ) +
  theme(
    # add margin on the left of the drawing canvas,
    # so title is aligned with left edge of first plot
    plot.margin = margin(0, 0, 0, 7)
  )

plot_grid(
  plot4title, plot4,
  ncol = 1,
  # rel_heights values control vertical title margins
  rel_heights = c(0.1, 1)
)