Try to fill different counties with different colors for okcounty data.

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

Plot tornado from different years in different colors in the county map.

tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  dplyr::select(om, yr, date)
tpath_16_21 <- tpath %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  dplyr::select(om, yr, date)
ggplot() +
  geom_sf(data = okcounty, fill = NA) +
  geom_sf(data = tpoint_16_21, aes(color = as.factor(yr))) +
  theme_bw()

countypnt <- st_join(tpoint_16_21, okcounty)

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

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(countymap) +
  geom_sf(aes( fill = tdens))+
  theme_bw()

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

tpoint_ok <- tpoint %>%
  filter(yr >= 2000 & yr <= 2021)

countypnt <- st_join(okcounty, st_as_sf(tpoint_ok, 
                     coords = c("slon", "slat"), 
                     crs = 4326)) %>%
  dplyr::select(GEOID, NAME, yr, om, inj, fat)

injurysum <- st_drop_geometry(countypnt) %>%
  group_by(GEOID) %>%
  summarize(
    total_injuries = sum(inj, na.rm = TRUE),
    years          = n_distinct(yr),
    avg_inj_year   = total_injuries / years
  )

injurymap <- okcounty %>%
  left_join(injurysum, by = "GEOID") %>%
  mutate(avg_inj_year = replace_na(avg_inj_year, 0))

ggplot(injurymap) +
  geom_sf(aes(fill = avg_inj_year), color = "white", linewidth = 0.3) +
  scale_fill_distiller(
    palette  = "YlOrRd",
    direction = 1,
    name     = "Avg Injuries\nper Year",
    labels   = scales::comma
    
  )

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.

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

ggplot(data = countymap) +
  geom_sf(aes(fill = tdens_c1), color = "white", linewidth = 0.3) +
  scale_fill_brewer(
    name    = expression("Tornadoes / 1000 km"^2),
    palette = "YlOrRd",
    drop    = FALSE       
  ) 

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

map_paths <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "grey60", linewidth = 0.3) +
  geom_sf(data = tpath_16_21, aes(color = as.factor(yr)), linewidth = 0.8) +
  scale_color_brewer(
    name    = "Year",
    palette = "Set1"
  ) +
  labs(title = "Tornado Paths (2016–2021)") +
  theme_bw() +
  theme(
    legend.position  = "bottom",
    legend.direction = "horizontal",
    plot.title       = element_text(face = "bold", size = 11, hjust = 0.5)
  )


map_points <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "grey60", linewidth = 0.3) +
  geom_sf(data = tpoint_16_21, aes(color = as.factor(yr)), size = 1.5, alpha = 0.7) +
  scale_color_brewer(
    name    = "Year",
    palette = "Set1"
  ) +
  labs(title = "Tornado Points (2016–2021)") +
  theme_bw() +
  theme(
    legend.position  = "bottom",
    legend.direction = "horizontal",
    plot.title       = element_text(face = "bold", size = 11, hjust = 0.5)
  )
shared_legend <- get_legend(
  map_paths +
    theme(legend.position = "bottom",
          legend.direction = "horizontal")
)

plot_grid(
  plot_grid(
    map_paths  + theme(legend.position = "none"),
    map_points + theme(legend.position = "none"),
    ncol       = 2,
    labels     = c("A", "B"),
    label_size = 12
  ),
  shared_legend,
  ncol         = 1,
  rel_heights  = c(1, 0.1))