library(tidyverse)
library(tigris)
library(scales)
library(sf)

### Plots the NWS Categorical Risk map for Middle TN for the given day
### Data from:  https://www.spc.noaa.gov/gis/

# Which day X categorical map do you want?  Valid input is "1", "2", or "3"
# Note:  This yields the weather map for Thursday April 18, 2024
day <- "2"

### Limiting the map to just neighboring states to speed up rendering
tn_neighboring_states <-
  c("01", "05", "13", "17", "18", "21", "28", "29", "37", "39", "45", "47", "51")

### Cache TIGRIS data to speed up future map renders
options(tigris_use_cache = TRUE)

### Bounding box for Mid-TN
map_x_min <- -81.5
map_x_max <- -90.5
map_y_min <-  34.5
map_y_max <-  37.0

### Load shapefile with the given day's outlook from NOAA
map_url <- paste("https://www.spc.noaa.gov/products/outlook/day", day, "otlk_cat.kmz", sep = "")

### Create a local map using a bounding box, and change colors to their standard representation
map_weather <- map_url %>% read_sf() %>% st_transform(crs = "NAD83")

# Pull the start/end valid dates for the map and convert from UTC to CST
weather_start <- map_weather %>% st_drop_geometry() %>% select(begin) %>% unique() %>% pull(begin)
weather_end   <- map_weather %>% st_drop_geometry() %>% select(end)   %>% unique() %>% pull(end)

### Load state/county outline maps from TIGRIS
map_county <-
  tigris::counties() %>%
  filter(STATEFP %in% tn_neighboring_states)

map_state <-
  tigris::states()

weather_map <-
  ggplot() +
  theme_bw() +
  theme(
    plot.title      = element_text(hjust = 0.5, size = 24),
    plot.subtitle   = element_text(hjust = 0.5, size = 16),
    legend.title    = element_text(size = 16),
    legend.text     = element_text(size = 15),
    legend.position = "bottom",
    legend.margin   = margin(-10, 0, 0, 0)
  ) +
  labs(x = "", y = "",
       title = paste("NWS Day ", day, " Categorical Risk Map", sep = ""),
       subtitle = paste("(valid from ", weather_start, " until ", weather_end, ")", sep = "")) +
  geom_sf(data = map_weather, aes(fill = Name), color = "black", size = 0.2) +
  geom_sf(data = map_county, fill = NA, size = 0.15) +
  geom_sf(data = map_state, fill = NA, color = "black", size = 3) +
  scale_fill_manual(
    name = "Categorical Risk", drop = FALSE,
    values = c(
      "General Thunder" = "#c0e8c0",
      "Marginal Risk"   = "#7fc57f",
      "Slight Risk"     = "#f6f67f",
      "Enhanced Risk"   = "#e6c27f",
      "Moderate Risk"   = "#e67f7f",
      "High Risk"       = "#ff7fff"
    ),
    limits = c(
      "General Thunder", "Marginal Risk", "Slight Risk", "Enhanced Risk", "Moderate Risk", "High Risk"
    ),
    labels = c(
      "General Thunder", "Marginal Risk", "Slight Risk", "Enhanced Risk", "Moderate Risk", "High Risk"
    )
  ) +
  coord_sf(xlim = c(map_x_min, map_x_max), ylim = c(map_y_min, map_y_max), expand = FALSE)

print(weather_map)