LabA

Introduction

This report documents information pertaining to Tornadoes which occurred within Oklahoma during the 2016 to 2021 time frame. Information regarding the locations of tornadoes includes the starting points and paths for each year as well as the density for each county.

Prepare The Date

Load Packages

The first step was to load packages.

library(sf) #to work with spatial data
library(rgdal) #to work with spatial data
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales) #pretty_breaks functions
library(RColorBrewer)
library(units) # with sf
library(cowplot) #plot_grid function
library(here)

Read In Oklahoma Tornadoes Geospatial Data

Next the data was imported into the R platform. The data derived from national-level datsets generated by NOAA National Weather Service Storm Prediction Center.

okcounty <- st_read("ok_counties.shp", quiet = TRUE)
tpoint <- st_read("ok_tornado_point.shp", quiet = TRUE)
tpath <- st_read("ok_tornado_path.shp", quiet = TRUE)
class(okcounty)

## [1] "sf"         "data.frame"

Prepare the Data

In order to review only the recent data, the dataset was modify to include tornado point and path data from the 2016 to 2021 time frame

tpoint_16_21 <- tpoint %>% ##Filter the tornado point data to only include data from 2016-2021
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

tpath_16_21 <- tpath %>% ##Filter the tornado path data to only include data from 2016-2021
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

Using spatial join tools, the number of tornado points within each county was calculated. Before grouping and summarizing could be preformed, countypmt needed to be converted from a sf object to a normal data fram using st_drop_geometry(). The the okcounty and countysum were joined to link the polygon to the tornado summary.

countypnt <- st_join(tpoint_16_21, okcounty)  ##Link row
countypnt <- st_drop_geometry(countypnt) ##Compute the total number of tornadoes per county
countysum <- countypnt %>% ## Calculate number of tornadoes in each county using summarize
  group_by(GEOID) %>%
  summarize(tcnt = n())

countymap <- okcounty %>%
  left_join(countysum, by = "GEOID") %>% ##Use left_join to ensure all county polygons are retained in the join.
  replace(is.na(.), 0) %>%  ##since some counties had no tornadoes, replace NS with zero.
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%  ##Compute the area of each county and calculate the density of tornadoes per county.
  drop_units() ##Remove the explicit measurements (like [m^2]) returned from st_area()

Results of Practice Exercises

A) Tornado Points & Paths Color Coded by Year

Generate a 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().

pathcolor <- ggplot() + ##Path by year in color
  geom_sf(data = tpath_16_21,  #map the tornadoes paths with color by year
          aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) + #map the county boundaries
  scale_color_discrete(name = "Year") + #set legend name 
  theme_bw()+ #display the map over a white background while retaining the graticules
  theme_void() # removes the plot axes and labels

pointcolor <- ggplot() +  ##Point by year in color
  geom_sf(data = tpoint_16_21, #map the tornadoes points with color by year
          aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) + #map the county boundaries
  scale_color_discrete(name = "Year") + #set legend name 
  coord_sf(datum = NA) +
  theme_void() #removes the plot axes and labels

##Combine Point & Path Maps
plot_grid(pointcolor, pathcolor, #Combine point and path colored maps
          labels = c("A) Tornado Starting Point by Year", #add labels
                     "B) Tornado Paths by Year",
                     label_size = 12), #specify label size
          ncol = 1, #single column
          hjust = 0, #justify labels
          label_x = 0, #move label to top of each map
          align = "hv") # align maps horizontally and vertically

## B) Summary of Density for Tornado Points by County & Year 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.

##Count by Year
countypnt <- st_join(tpoint_16_21, okcounty) #point data joined with county shapefile
countypnt_yr <- st_drop_geometry(countypnt) #geometry is dropped so operation can be preformed by date
countysum_yr <- countypnt_yr %>% #group by GEIOD and year
  group_by(GEOID, yr) %>%
  summarize(tcnt=n()) %>% #count each point and summarized for each county and year
  aggregate(tcnt ~ yr + GEOID, FUN=sum, na.rm=TRUE)

##Summarized Data Joined with County Map
countymap_sumyr <- okcounty %>%
  left_join(countysum_yr, by = "GEOID") %>% #summarized data is re-joined with county shapefile
  replace(is.na(.), 0) %>%  ##since some counties had no tornadoes, replace NS with zero.
  mutate(area = st_area(.),
         tdens2 = 10^6 * 10^3 * tcnt / area) %>%  ##Compute the area of each county and calculate the density of tornadoes per county.
  drop_units() ##Remove the explicit measurements (like [m^2]) returned from st_area()


##Map Tornado Density per county by year
ggplot() +
  geom_sf(data = okcounty, #layer for the county borders
          fill = NA,
          color = "gray")+
  geom_sf(data = countymap_sumyr, aes(fill=tdens2))+  #layer for the density per county
  facet_wrap(vars(yr), ncol = 2)+ #organize by year
  coord_sf(datum = NA)+
  scale_fill_distiller(name=expression("Tornadoes/1000 km"^2), #specify color ramp from the RColorBrewer and set the name of of the argument
                       palette = "YlOrRd", #Yellow to red
                       breaks=pretty_breaks(), #automatically set breaks for the legend
                       direction = 1)+
  theme_void() #removes the plot axes and labels

## C) Tornado Density Based on Quantile Breaks 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.

#3b - Three discrete classes
numclas_3b <- 3
qbrks_3b <- seq(0, 1, length.out = numclas_3b + 1)
qbrks_3b

## [1] 0.0000000 0.3333333 0.6666667 1.0000000

#Results: [1] 0.0000000 0.3333333 0.6666667 1.0000000

countymap_3b <- countymap %>%
  mutate(tdens_c1 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks_3b),
                        include.lowest = T))

#3b - Three Discrete Map - saved as an R object
q_3_map <- ggplot(data = countymap_3b) +
  geom_sf(aes(fill = tdens_c1)) +
  scale_fill_brewer(name = "", palette = "YlOrRd") +
  labs(title = "Tornadoes/1000 km^2")+
  theme_void() +
  theme(legend.position = "bottom",
        text = element_text(colour = "#525252"))

#4b - Four discrete classes
numclas_4b <- 4
qbrks_4b <- seq(0, 1, length.out = numclas_4b + 1)
qbrks_4b

## [1] 0.00 0.25 0.50 0.75 1.00

#Results: [1] 0.00 0.25 0.50 0.75 1.00

countymap_4b <- countymap %>%
  mutate(tdens_c1 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks_4b),
                        include.lowest = T))

#4b - Four Discrete Map - saved as an R object
q_4_map <- ggplot(data = countymap_4b) +
  geom_sf(aes(fill = tdens_c1)) +
  scale_fill_brewer(name = "", palette = "YlOrRd") +
  labs(title = "Tornadoes/1000 km^2")+
  theme_void() +
  theme(legend.position = "bottom",
        text = element_text(colour = "#525252"))

#5b - Five discrete classes
numclas_5b <- 5
qbrks_5b <- seq(0, 1, length.out = numclas_5b + 1)
qbrks_5b

## [1] 0.0 0.2 0.4 0.6 0.8 1.0

#Results: [1] 0.0 0.2 0.4 0.6 0.8 1.0

countymap_5b <- countymap %>%
  mutate(tdens_c1 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks_5b),
                        include.lowest = T))

#5b - Five Discrete Map - saved as an R object
q_5_map <- ggplot(data = countymap_5b) +
  geom_sf(aes(fill = tdens_c1)) +
  scale_fill_brewer(name = "", palette = "YlOrRd") +
  labs(title = "Tornadoes/1000 km^2")+
  theme_void() +
  theme(legend.position = "bottom",
        text = element_text(colour = "#525252"))

#6b - Six discrete classes
numclas_6b <- 6
qbrks_6b <- seq(0, 1, length.out = numclas_6b + 1)
qbrks_6b

## [1] 0.0000000 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333 1.0000000

#Results: [1] 0.0000000 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333 1.0000000

countymap_6b <- countymap %>%
  mutate(tdens_c1 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks_6b),
                        include.lowest = T))

#6b - Six Discrete Map - saved as an R object
q_6_map <- ggplot(data = countymap_6b) +
  geom_sf(aes(fill = tdens_c1)) +
  scale_fill_brewer(name = "", palette = "YlOrRd") +
  labs(title = "Tornadoes/1000 km^2")+
  theme_void() +
  theme(legend.position = "bottom",
        text = element_text(colour = "#525252"))

#Combine Discrete Map 4 through 6 

plot_grid(q_3_map, q_4_map, q_5_map, q_6_map,
          labels = c("Three Classes",
                     "Four Classes", 
                     "Five Classes",
                     "Six Classes",
                     label_size = 10),
          ncol = 2,
          hjust = 0,
          label_x = 0, 
          align = "hv")

## Unknown Issues

The tiles for the figures appeared ontop of data within the figres. I attempted to change this within the code but could not figure it out. Based on what I observed, my code was similar (or the exact same) as depicted in the lesson.