lab_hw14

library(cowplot)

## Warning: 程序包'cowplot'是用R版本4.5.3 来建造的

library(sf)           

## Warning: 程序包'sf'是用R版本4.5.3 来建造的

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(ggplot2)      
library(dplyr)        

## 
## 载入程序包：'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)        
library(scales)       
library(RColorBrewer) 
library(units)

## Warning: 程序包'units'是用R版本4.5.3 来建造的

## udunits database from C:/Users/User/AppData/Local/R/win-library/4.5/units/share/udunits/udunits2.xml

library(cowplot)

path <- "D:/lilith/Chapter5"
okcounty <- st_read(paste0(path,"/ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(paste0(path,"/ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(paste0(path,"/ok_tornado_path.shp"), quiet = TRUE)

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)

countypnt <- st_join(tpoint_16_21, okcounty)

ggplot(data = okcounty) +
  geom_sf(aes(fill = NAME))

ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "gray") +
  
  geom_sf(data = tpoint_16_21,
          aes(color = factor(yr)),size = 1) +
  
  theme_void()

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

## Rows: 75
## Columns: 2
## $ GEOID <chr> "40001", "40005", "40007", "40009", "40011", "40013", "40015", "~
## $ tcnt  <int> 6, 3, 4, 8, 1, 4, 10, 5, 7, 5, 3, 12, 10, 5, 5, 1, 7, 9, 7, 8, 2~

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()

ggplot(data = countymap) +
  
  geom_sf(aes(fill = tdens)) +
  
  scale_fill_distiller(
    palette = "YlOrRd",
    direction = 1,
    name = expression("Tornadoes / 1000 km"^2)
  ) +
  
  theme_void() +
  
  theme(legend.position = "bottom")

。

countypnt <- st_join(tpoint, okcounty)

countypnt <- st_drop_geometry(countypnt)

countysum <- countypnt %>%
  group_by(GEOID) %>%
  summarize(avg_inj = sum(inj) / 6)

countymap <- okcounty %>%
  left_join(countysum, by = "GEOID") %>%
  replace(is.na(.), 0)

ggplot(data = countymap) +
  geom_sf(aes(fill = avg_inj)) +
  
  scale_fill_distiller(
    palette = "YlOrRd",
    direction = 1,
    name = "Average Injuries per Year"
  ) +
  
  theme_void() +
  theme(legend.position = "bottom")

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

## Rows: 78
## Columns: 2
## $ GEOID <chr> "40001", "40003", "40005", "40007", "40009", "40011", "40013", "~
## $ tcnt  <int> 28, 39, 42, 67, 69, 39, 33, 117, 95, 56, 37, 26, 49, 66, 26, 60,~

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()

countymap <- countymap %>%
  mutate(tdens_c1 = cut(
    tdens,
    breaks = c(0, 2, 4, 6, 8, 10, Inf),
    include.lowest = TRUE
  ))

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

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)

path_map <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "gray") +
  
  geom_sf(data = tpath_16_21,
          aes(color = factor(yr)),
          size = 0.6) +
  
  theme_void() +
  labs(color = "Year")

point_map <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "gray") +
  
  geom_sf(data = tpoint_16_21,
          aes(color = factor(yr)),
          size = 1) +
  
  theme_void() +
  labs(color = "Year")

plot_grid(
  path_map,
  point_map,
  labels = c("A) Tornado Paths", "B) Tornado Points"),
  ncol = 1,
  align = "v"
)