Introduction

This document shows maps of tornado paths and densities in Oklahoma.

1. Generate a map of tornado paths where paths from each year are displayed as a different color.

Read in data and packages

knitr::opts_chunk$set(message = FALSE, warning = FALSE)
# import packages
packages <- c("tidyverse", "here", "ggplot2", "dplyr", "tidyr", "scales", "RColorBrewer", "units", "cowplot", "sf")
invisible(lapply(packages, library, character.only=TRUE)) # dont show package import
# read in tornado data
okcounty <- st_read(here("Data/Chapter5", "ok_counties.shp"))
## Reading layer `ok_counties' from data source 
##   `C:\Users\kmayo\Documents\588\data\Chapter5\ok_counties.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 77 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -103.0025 ymin: 33.62184 xmax: -94.43151 ymax: 37.00163
## Geodetic CRS:  NAD83
tpoint <- st_read(here("Data/Chapter5", "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(here("Data/Chapter5","ok_tornado_path.shp"), quiet = TRUE)

Create map of tornado paths by year

# Create map with tornado paths distinguished by year (different color)
ggplot() +
  geom_sf(data=tpath, # use path data
          aes(color=as.factor(yr))) + # use yr as color determinant
  geom_sf(data=okcounty, fill=NA) + # show ok county map on map too
  coord_sf(datum=NA)+
  theme_void() # remove plot axes

### Create plot for tornado points & paths, and join them into one with plot_grid

# Create map for tornado points
pointmap<-ggplot() + # name new plot,
  geom_sf(data=okcounty, fill=NA) + # using the ok county map outlines,
  geom_sf(data=tpoint, # show tornado point data
          color="red", # in red
          size=.5)+ # size 0.5
  theme_bw() # use black and white backdrop theme
# Create map for tornado paths
pathmap<-ggplot() +
  geom_sf(data=okcounty, fill=NA) + # using the ok county map, unfilled as background,
  geom_sf(data=tpath,# show tornado path data
          color="blue", # in blue
          size=.75) + # size 0.75
  theme_bw() # use black and white backdrop theme.

# Plot both maps into one output
plot_grid(pointmap,pathmap, # plot the two maps we just made
          labels=c("A) Tornado Point Map", # name map 1
                   "B) Tornado Path Map"), # name map 2
                   label_size=12, # font size
          ncol=1, # number of columns
          hjust=0, # put label on left
          align="hv")

2. Summarise tornado point density by county and year & create composite figure

Summarise tornado density from 2016-2021

# Re read tpoint, coming in as data frame?
class(tpoint)
## [1] "sf"         "data.frame"
tpoint <- st_read(here("Data/Chapter5", "ok_tornado_point.shp"), quiet = TRUE)
# Filter only years 2016-2021
tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date) # om = OID
# Spatial join so points are by county
countypnt <- st_join(tpoint_16_21, okcounty)
# Group geom and count points per county & year
countysum_yr<-countypnt%>%
  st_drop_geometry() %>%
  group_by(GEOID,yr) %>% # for each county how many tornadoes in each year
  summarize(tcnt = n())
# Drop geom, summarise by county so each polygon has above sum data
countymap_yr<-okcounty%>%
  left_join(countysum_yr, by = "GEOID") %>%
  replace(is.na(.), 0) %>% # change missing na values to 0
  mutate(area = st_area(geometry), # add a new column for county area
         tdens = 10^6 * 10^3 * tcnt / area) %>% # calc tornado density per 1000km, fill column
  drop_units()
# Generate a faceted plot to display county-level tornado density from 2016-2021. 1 map for each year
ggplot(data=countymap_yr)+
  geom_sf(aes(fill=tdens))+ # county fill color based on tornado density stat
  scale_fill_distiller(name=expression("Tornadoes/1000km"^2), # Add scale name
                       palette="YlOrRd",
                       breaks=pretty_breaks(), # auto break scale steps
                       direction=1)+
  facet_wrap(vars(yr),ncol=2)+ # create map for each year
  coord_sf(datum=NA)+
  theme_void()+
  theme(legend.position="bottom")

3. Create 4 choropleth maps of tornado density based on quantile breaks.

make_map_quantile <- function(numclas) {
  # quantile 0-1 w 4 breaks
  qbrks <- seq(0, 1, length.out = numclas + 1)

  # separate tdens into numclas
  choropleth <- countymap_yr %>%
    mutate(
      tdens_c = cut(tdens,
                    breaks = quantile(tdens, probs = qbrks, na.rm = TRUE),
                    include.lowest = TRUE))
  # make map with tdense_c quantile
  ggplot(choropleth) +
    geom_sf(aes(fill = tdens_c)) +
    scale_fill_brewer(
      name = expression("Tornadoes/1000 km"^2),
      palette = "YlOrRd") +
    theme(
      legend.key.size=unit(0.3,"cm"),
      legend.text = element_text(size=4)) +
    theme(legend.position = "bottom")}

# Run through function for each quantile
p3 <- make_map_quantile(3)
p4 <- make_map_quantile(4)
p5 <- make_map_quantile(5)
p6 <- make_map_quantile(6)

# Composite maps
cowplot::plot_grid(
  p3, p4, p5, p6,
  labels = c("A) 3", "B) 4", "C) 5", "D) 6"),
  label_size = 8,
  ncol = 2, hjust = 0, label_x = 0, align = "hv")

Discussion

Fewer quantile classes make differences look stronger but lose detail. More classes give more information and detail but end up being more difficult to interpret with more classes.