Plot the okcounty map with each county color-filled by the density of tornadoes between 2000 and 2021:

path <- "/Users/HoangDucVinh/Downloads/gdswr_data/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_00_21 <- tpoint %>%
  filter(yr >= 2000 & yr <= 2021)

countypnt_00_21 <- st_join(tpoint_00_21, okcounty) %>%
  st_drop_geometry() %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())

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

ggplot(data = countymap_00_21) +
  geom_sf(aes(fill = tdens)) +
  scale_fill_distiller(name = "Tornadoes/1000 km^2", 
                       palette = "YlOrRd", 
                       direction = 1)


Choose a proper palette map to visualize the average injuries per year by tornadoes in each county between 2000 and 2021:

injury_sum <- tpoint %>%
  filter(yr >= 2000 & yr <= 2021) %>%
  st_join(okcounty) %>%
  st_drop_geometry() %>%
  group_by(GEOID) %>%
  summarize(total_inj = sum(inj, na.rm = TRUE)) %>%
  mutate(avg_inj_yr = total_inj / 22)

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

ggplot(data = injury_map) +
  geom_sf(aes(fill = avg_inj_yr)) +
  scale_fill_distiller(palette = "Reds", direction = 1) +
  labs(
    title = "Average Annual Tornado Injuries by County",
    fill = "Injuries/Year"
  )


Reproduce the map in the previous page with breaking density with 6 groups - 0 to 2, 2 to 4, 4 to 6, 6 to 8, 8 to 10 and above 10:

brks <- c(0, 2, 4, 6, 8, 10, Inf)
labs <- c("0-2", "2-4", "4-6", "6-8", "8-10", ">10")

countymap_00_21 <- countymap_00_21 %>%
  mutate(tdens_discrete = cut(tdens, breaks = brks, labels = labs, include.lowest = TRUE))

ggplot(data = countymap_00_21) +
  geom_sf(aes(fill = tdens_discrete)) +
  scale_fill_brewer(palette = "YlOrRd") +
  theme(legend.position = "bottom") +
  labs(
    title = "Tornado Density (Categorized)",
    fill = "Tornadoes/1000 km^2"
  )


Generate a map of tornado paths where the paths from each year (2016-2021) 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():

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_plot <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "black") +
  geom_sf(data = tpath_16_21, aes(color = as.factor(yr)), size = 0.8) +
  labs(
    title = "Tornado Paths (2016-2021)",
    color = "Year"
  )

point_plot <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "black") +
  geom_sf(data = tpoint_16_21, aes(color = as.factor(yr)), size = 1) +
  labs(
    title = "Tornado Touchdown Points (2016-2021)",
    color = "Year"
  )

plot_grid(path_plot, point_plot, 
          ncol = 1)