Lesson 4 Part A

The first part of this lab examines point and line features representing tornado touchdown locations and travel paths respectively. The simple features and ggplot libraries were used to map and visualize the tornado dataset. Figure 1 plots a composite of touchdown points and paths symbolized with varying color for each year from 1950-2021. Figure 2 maps tornado touchdown density by county for each year from 2016-2021. Figure 3 shows an attempt at setting custom symbology breaks. This attempt was unsuccessful.

#Load package libraries

library(sf)
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
library(RColorBrewer)
library(units)
library(cowplot)
library(here)
## here() starts at C:/PSU/GEOG588/WS_Lesson4

#Read in 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)

#Map tornado paths & points then save as objects

pathmap <- ggplot() +
  geom_sf(data = tpath, 
          aes(color = as.factor(yr)), show.legend = FALSE) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

pointmap <- ggplot() +
  geom_sf(data = tpoint, 
          aes(color = as.factor(yr)), show.legend = FALSE) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

#Create composite plot

plot_grid(pathmap, pointmap, 
          labels = c("A) Tornado paths 1950-2021", 
                     "B) Tornado points 1950-2021",
                     label_size = 12),
          ncol = 1, 
          hjust = 0, 
          label_x = 0, 
          align = "hv")



Figure 1.) A composite plot of tornado touchdown points and paths with year values from 1950-2021 symbolized with unique color grouping.

#Iterate over annual data, creating a density layer for each year from 2016 to 2021

# filter year 2016
tpoint_16 <- tpoint %>%
  filter(yr == 2016) %>%
  select(om, yr, date)

# summarize county w/ points for 2016 (will have no data values)
countypnt16 <- st_join(tpoint_16, okcounty)
countypnt16 <- st_drop_geometry(countypnt16)
countysum16 <- countypnt16 %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())

# build county layer summarized for 2016 (fill no data values from previous step)
countymap16 <- okcounty %>%
  left_join(countysum16, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

# append year column lost during summary with appropriate value
countymap16['yr'] <- '2016'

# filter year 2017
tpoint_17 <- tpoint %>%
  filter(yr == 2017) %>%
  select(om, yr, date)

# summarize county w/ points for 2017 (will have no data values)
countypnt17 <- st_join(tpoint_17, okcounty)
countypnt17 <- st_drop_geometry(countypnt17)
countysum17 <- countypnt17 %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())

# build county layer summarized for 2017 (fill no data values from previous step)
countymap17 <- okcounty %>%
  left_join(countysum17, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

# append year column lost during summary with appropriate value
countymap17['yr'] <- '2017'

# filter year 2018
tpoint_18 <- tpoint %>%
  filter(yr == 2018) %>%
  select(om, yr, date)

# summarize county w/ points for 2018 (will have no data values)
countypnt18 <- st_join(tpoint_18, okcounty)
countypnt18 <- st_drop_geometry(countypnt18)
countysum18 <- countypnt18 %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())

# build county layer summarized for 2018 (fill no data values from previous step)
countymap18 <- okcounty %>%
  left_join(countysum18, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

# append year column lost during summary with appropriate value
countymap18['yr'] <- '2018'

# filter year 2019
tpoint_19 <- tpoint %>%
  filter(yr == 2019) %>%
  select(om, yr, date)

# summarize county w/ points for 2019 (will have no data values)
countypnt19 <- st_join(tpoint_19, okcounty)
countypnt19 <- st_drop_geometry(countypnt19)
countysum19 <- countypnt19 %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())

# build county layer summarized for 2019 (fill no data values from previous step)
countymap19 <- okcounty %>%
  left_join(countysum19, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

# append year column lost during summary with appropriate value
countymap19['yr'] <- '2019'

# filter year 2020
tpoint_20 <- tpoint %>%
  filter(yr == 2020) %>%
  select(om, yr, date)

# summarize county w/ points for 2020 (will have no data values)
countypnt20 <- st_join(tpoint_20, okcounty)
countypnt20 <- st_drop_geometry(countypnt20)
countysum20 <- countypnt20 %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())

# build county layer summarized for 2020 (fill no data values from previous step)
countymap20 <- okcounty %>%
  left_join(countysum20, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

# append year column lost during summary with appropriate value
countymap20['yr'] <- '2020'

# filter year 2021
tpoint_21 <- tpoint %>%
  filter(yr == 2021) %>%
  select(om, yr, date)

# summarize county w/ points for 2021 (will have no data values)
countypnt21 <- st_join(tpoint_21, okcounty)
countypnt21 <- st_drop_geometry(countypnt21)
countysum21 <- countypnt21 %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())

# build county layer summarized for 2021 (fill no data values from previous step)
countymap21 <- okcounty %>%
  left_join(countysum21, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

# append year column lost during summary with appropriate value
countymap21['yr'] <- '2021'

# compile county layers using rbind()
countymap <- rbind(countymap16, countymap17, countymap18, countymap19, countymap20, countymap21)

# sum dens by county & year
ggplot(data = countymap) +
  geom_sf(aes(fill = tdens)) +
  facet_wrap(vars(yr), ncol = 2) +
  coord_sf(datum = NA) +
  theme_void()

Figure 2.) Here tornado touchdown points are joined with county boundaries to create a summarized density by county and year.

#Create 4 choropleth maps
#First break manipulation

# filter year range ----
tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

countypnt <- st_join(tpoint_16_21, okcounty)

countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())

countymap <- okcounty %>%
  left_join(countysum, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

numclas <- 4
qbrks <- seq(0, 1, length.out = numclas + 1)
qbrks
## [1] 0.00 0.25 0.50 0.75 1.00
countymap <- countymap %>%
  mutate(tdens_c1 = cut(tdens,
                        breaks = quantile(tdens, breaks = qbrks),
                        include.lowest = T))

densTry1 <- ggplot(data = countymap) +
  geom_sf(aes(fill = tdens_c1)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom")

#Second manipulation

countypnt <- st_join(tpoint_16_21, okcounty)

countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())

countymap <- okcounty %>%
  left_join(countysum, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

numclas <- 6
qbrks <- seq(0, .5, length.out = numclas + 1)
qbrks
## [1] 0.00000000 0.08333333 0.16666667 0.25000000 0.33333333 0.41666667 0.50000000
countymap <- countymap %>%
  mutate(tdens_c1 = cut(tdens,
                        breaks = quantile(tdens, breaks = qbrks),
                        include.lowest = T))

densTry2 <- ggplot(data = countymap) +
  geom_sf(aes(fill = tdens_c1)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom")

#Plot density attempts

plot_grid(densTry1, densTry2, 
          labels = c("Density Attempt 1", 
                     "Density Attempt 2",
                     label_size = 12),
          ncol = 1, 
          hjust = 0, 
          label_x = 0, 
          align = "hv")


Figure 3.) This plot attempts to show variation among breaks. In each attempt, the number class and quantile break values were manipulated but I did not observe variation within the plots.