Wimberly Tornado Practice Exercises

load in packages

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

the example code was applied but hidden for this report

Practice Exercise 1

generate map of tornado paths where the paths from each year are displayed as a different color. Create a composite figure containing the map of tornado paths and the map of tornado points using plot_grid.
#plot tornado paths by year
paths <- ggplot() +
  geom_sf(data = tpath_16_21, 
          aes(color = as.factor(yr)), 
          size = 2) +
  #creates county layer
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

#plot tornado points by year
points <- ggplot() +
  geom_sf(data = tpoint_16_21, 
          aes(color = as.factor(yr))) +
  #creates county layer
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

#save the paths plot
ggsave("paths.tiff", 
       plot = paths,
       width = 6, 
       height = 4, 
       dpi = 300, 
       compression = "lzw")

#save the points plot
ggsave("points.tiff",
       plot = points,
       width = 6, 
       height = 4, 
       dpi = 300, 
       compression = "lzw")

#creates grid of paths and points plot
plot_grid(labels = c("A) Tornado Paths Map", 
                     "B) Tornado Points Map",
                     label_size = 12),
          paths, points, 
          ncol = 1, 
          hjust = 0, 
          label_x = 0, 
          align = "hv")
Fig 1. Tornado paths and points in Oklahoma from 2016 to 2021.

Fig 1. Tornado paths and points in Oklahoma from 2016 to 2021.

Practice Exercise 2

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.
#join tornado points to counties
countypnt <- st_join(tpoint_16_21, okcounty)
glimpse(countypnt)
## Rows: 434
## Columns: 11
## $ om       <dbl> 613662, 613675, 613676, 613677, 613678, 613727, 613728, 61372…
## $ yr       <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2…
## $ date     <chr> "2016-03-23", "2016-03-30", "2016-03-30", "2016-03-30", "2016…
## $ STATEFP  <chr> "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "…
## $ COUNTYFP <chr> "001", "113", "105", "131", "035", "139", "139", "139", "139"…
## $ COUNTYNS <chr> "01101788", "01101844", "01101840", "01101853", "01101805", "…
## $ AFFGEOID <chr> "0500000US40001", "0500000US40113", "0500000US40105", "050000…
## $ GEOID    <chr> "40001", "40113", "40105", "40131", "40035", "40139", "40139"…
## $ NAME     <chr> "Adair", "Osage", "Nowata", "Rogers", "Craig", "Texas", "Texa…
## $ LSAD     <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "…
## $ geometry <POINT [°]> POINT (-94.5042 35.6916), POINT (-96.0151 36.2151), POI…
#drop geometry from countypnt file
countypnt <- st_drop_geometry(countypnt)
#group by county GEOID and year
countysum2 <- countypnt %>%
  group_by(GEOID, yr) %>%
  #summarize by tornado count
  summarize(tcnt = n())
glimpse(countysum2)
## Rows: 217
## Columns: 3
## Groups: GEOID [75]
## $ GEOID <chr> "40001", "40001", "40001", "40001", "40001", "40005", "40005", "…
## $ yr    <dbl> 2016, 2017, 2019, 2020, 2021, 2016, 2019, 2017, 2018, 2019, 2017…
## $ tcnt  <int> 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 3, 1, 3, 1, 1, 1, 1, 2, 1, 5, 4, 1…
#left join county sum by GEOID and year to county shapefile
okcountyyr <- okcounty %>%
  left_join(countysum2, by = "GEOID")

countymap2 <- okcountyyr %>%
  #replace 0 with no values
  replace(is.na(.), 0) %>%
  #calculate tornado density
  mutate(area = st_area(okcountyyr),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()
glimpse(countymap2) 
## Rows: 219
## Columns: 12
## $ STATEFP  <chr> "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "…
## $ COUNTYFP <chr> "077", "025", "025", "025", "025", "011", "107", "107", "107"…
## $ COUNTYNS <chr> "01101826", "01101800", "01101800", "01101800", "01101800", "…
## $ AFFGEOID <chr> "0500000US40077", "0500000US40025", "0500000US40025", "050000…
## $ GEOID    <chr> "40077", "40025", "40025", "40025", "40025", "40011", "40107"…
## $ NAME     <chr> "Latimer", "Cimarron", "Cimarron", "Cimarron", "Cimarron", "B…
## $ LSAD     <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "…
## $ yr       <dbl> 2020, 2016, 2017, 2019, 2021, 2017, 2017, 2019, 2020, 2021, 2…
## $ tcnt     <dbl> 1, 2, 2, 1, 7, 1, 2, 6, 1, 1, 1, 1, 4, 1, 1, 1, 1, 1, 1, 2, 0…
## $ geometry <POLYGON [°]> POLYGON ((-95.50766 35.0292..., POLYGON ((-103.0025 3…
## $ area     <dbl> 1890663261, 4766283042, 4766283042, 4766283042, 4766283042, 2…
## $ tdens    <dbl> 0.5289149, 0.4196142, 0.4196142, 0.2098071, 1.4686497, 0.4120…
#filter out years with no value
countymap3 <- filter(countymap2, yr %in% c(2016, 2017, 2018, 2019, 2020, 2021))

#facet plot for tornado density by year
ggplot() +
  ggtitle("Tornado Density by Year and County") +
  geom_sf(data = okcounty,
          fill = NA,
          color = "gray") +
  geom_sf(data = countymap3,
          aes(fill = tdens)) +
  facet_wrap(vars(yr), ncol = 2) +
  scale_fill_distiller(name = expression("Tornadoes/100 km"^2),
                       palette = "YlOrRd",
                       breaks = pretty_breaks(),
                       direction = 1) +
  coord_sf(datum = NA) +
  theme_void() + 
  theme(legend.position = "bottom")
Fig 2. Tornado density in Oklahoma from 2016 to 2021.

Fig 2. Tornado density in Oklahoma from 2016 to 2021.

Practice Exercise 3

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.

Practice Exercise 3: 3 class break plot

#calculate 3 class breaks
numclas3 <- 3
#calculates quantiles directly
qbrks3 <- quantile(countymap3$tdens, probs = seq(0, 1, length.out = numclas3 + 1))
qbrks3
##        0% 33.33333% 66.66667%      100% 
## 0.1679519 0.5038558 0.9933283 5.8610691
#make column for 3 class breaks
countymap3 <- countymap3 %>%
  mutate(tdens_c3 = cut(tdens,
                        breaks = qbrks3,
                        include.lowest = T))

#plot of 3 class breaks for tornado density
Breaks3 <- ggplot() +
  geom_sf(data = okcounty,
          fill = NA,
          color = "gray") +
  geom_sf(data = countymap3,
          aes(fill = tdens_c3)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom")

Practice Exercise 3: 5 class break plot

#calculate 5 class breaks
numclas5 <- 5
#calculates quantiles directly
qbrks5 <- quantile(countymap3$tdens, probs = seq(0, 1, length.out = numclas5 + 1))
qbrks5
##        0%       20%       40%       60%       80%      100% 
## 0.1679519 0.4250639 0.5445558 0.8186329 1.3892310 5.8610691
#make column for 5 class breaks
countymap3 <- countymap3 %>%
  mutate(tdens_c5 = cut(tdens,
                        breaks = qbrks5,
                        include.lowest = T))

#plot of 5 class breaks for tornado density
Breaks5 <- ggplot() +
  geom_sf(data = okcounty,
          fill = NA,
          color = "gray") +
  geom_sf(data = countymap3,
          aes(fill = tdens_c5)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom")

Practice Exercise 3: 6 class break plot

#calculate 6 class breaks
numclas6 <- 6
#calculates quantiles directly
qbrks6 <- quantile(countymap3$tdens, probs = seq(0, 1, length.out = numclas6 + 1))
qbrks6
##        0% 16.66667% 33.33333%       50% 66.66667% 83.33333%      100% 
## 0.1679519 0.4120108 0.5038558 0.6590064 0.9933283 1.5385952 5.8610691
#make column for 6 class breaks
countymap3 <- countymap3 %>%
  mutate(tdens_c6 = cut(tdens,
                        breaks = qbrks6,
                        include.lowest = T))

#plot of 6 class breaks for tornado density
Breaks6 <- ggplot() +
  geom_sf(data = okcounty,
          fill = NA,
          color = "gray") +
  geom_sf(data = countymap3,
          aes(fill = tdens_c6)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom")

Practice Exercise 3: build composite plot

#build composite
ggsave("Breaks3.tiff",
       plot = points,
       width = 10, 
       height = 8, 
       dpi = 300, 
       compression = "lzw")

ggsave("Breaks5.tiff",
       plot = points,
       width = 6, 
       height = 4, 
       dpi = 300, 
       compression = "lzw")

ggsave("Breaks6.tiff",
       plot = points,
       width = 6, 
       height = 4, 
       dpi = 300, 
       compression = "lzw")

plot_grid(Breaks3, Breaks5, Breaks6,
          labels = c("A) Three Cleass Breaks", 
                     "B) Five Class Breaks",
                     "C) Six Class Breaks",
                     label_size = 12),
          ncol = 1, 
          hjust = 0, 
          label_x = 0, 
          align = "hv")
Fig 3. Tornado density (2016 to 2021) plotted with three different class breaks. First, a class break of 3. Second, a class break of 5. Third, a class break of 6.

Fig 3. Tornado density (2016 to 2021) plotted with three different class breaks. First, a class break of 3. Second, a class break of 5. Third, a class break of 6.

Please note that chaning the number of class breaks will change how each counties tornado density is displayed, and interpreted. This can lead to different interpretations of the tornado density of each county.