---
title: "`r params$hurdat_id` `r params$storm_name` `r params$storm_year` STORM REPORT"
subtitle: "Generated on: `r Sys.Date()`"
format:
  pdf:
    theme: cosmo
    geometry:
      - top=0.75in
      - bottom=0.75in
      - left=0.75in
      - right=0.75in
      - footskip=0.3in
execute:
  warning: false
  message: false
params:
    storm_basin: AL
    storm_year: 1900
    storm_name: GALVESTON
    hurdat_id: AL011900
---

```{r setup, echo=FALSE}
library(ggplot2)
library(sf)
sf_use_s2(FALSE)
library(rnaturalearth)
library(dplyr)
library(scales)
library(knitr)
library(kableExtra)
library(patchwork)
library(lubridate)

source(file = "queries.R")

storm <- list(
    storm_basin = params$storm_basin,
    storm_year = params$storm_year,
    storm_name = params$storm_name
)

latest_normalization_year <- get_latest_normalization_year()

unique_lfs <- get_unique_lf_ids(storm)
hurdat_id <- get_hurdat_id(storm)

plot_theme <- theme(
  plot.title = element_text(hjust = 0.5, size = 10),
  axis.title = element_text(size = 10),
  axis.text = element_text(size = 8),
  legend.title = element_text(face = "bold"),
  legend.text = element_text(face = "bold")
)
```

```{r storm_overview, echo=FALSE}
# storm stats
# mmh, mmp, landfalls, hurdat summary (max cat, min press, max wind, hurricane points, total points)
# fatality data?

hurdat_summary <- get_best_track_summary(storm)

losses <- get_latest_aggregate_loss(storm)

storm_track <- get_hurdat_track(storm)

```

# Storm Overview

- **HURDAT2 ID** `r hurdat_id$hurdatid`
- **HURDAT2 Revision Date** `r hurdat_id$hurdat_rev`

## Normalized Cost

- **MMH24** `r dollar(losses$mmh)`
- **MMP24** `r dollar(losses$mmp)`

## HURDAT2 Summary

- **Date Range** `r min(storm_track$formatted_datetime)` to `r max(storm_track$formatted_datetime)`
- **Total Observations** `r hurdat_summary$total_observations`
- **Hurricane Observations** `r hurdat_summary$hurricane_observations`
- **Landfall Observations** `r hurdat_summary$hurdat_landfalls`
- **Maximum Saffir-Simpson Category** `r hurdat_summary$max_category`
- **Minimum Pressure** `r hurdat_summary$min_pressure` mb
- **Maximum Windspeed** `r hurdat_summary$max_windspeed` kts

## Storm Track

```{r track_map, echo=FALSE}
#| fig-width: 10
#| fig-height: 8
#| out-width: "100%"
  storm_track <- storm_track %>%
    mutate(
      hurricane_category = case_when(
        storm_status == "HU" & windspeed >= 64 & windspeed <= 82 ~ 1,
        storm_status == "HU" & windspeed >= 83 & windspeed <= 95 ~ 2,
        storm_status == "HU" & windspeed >= 96 & windspeed <= 112 ~ 3,
        storm_status == "HU" & windspeed >= 113 & windspeed <= 136 ~ 4,
        storm_status == "HU" & windspeed >= 137 ~ 5,
      ),
      status_label = case_when(
        storm_status == "TD" ~ "Tropical Depression (TD)",
        storm_status == "TS" ~ "Tropical Storm (TS)",
        hurricane_category == 1 ~ "Hurricane Category 1",
        hurricane_category == 2 ~ "Hurricane Category 2",
        hurricane_category == 3 ~ "Hurricane Category 3",
        hurricane_category == 4 ~ "Hurricane Category 4",
        hurricane_category == 5 ~ "Hurricane Category 5",
        storm_status == "EX" ~ "Extratropical Cyclone (EX)",
        storm_status == "SD" ~ "Subtropical Depression (SD)",
        storm_status == "SS" ~ "Subtropical Storm (SS)",
        storm_status %in% c("LO", "WV", "DB") ~ "Low/Wave/Disturbance",
        TRUE ~ "Missing Data"
      ),
      line_color = case_when(
        storm_status == "TD" ~ "#2AFF00",
        storm_status == "TS" ~ "#FFD020",
        hurricane_category == 1 ~ "#FF4343",
        hurricane_category == 2 ~ "#FF6FFF",
        hurricane_category == 3 ~ "#FF23D3",
        hurricane_category == 4 ~ "#C916FF",
        hurricane_category == 5 ~ "#FFFFFF",
        storm_status == "EX" ~ "#202020",
        storm_status == "SD" ~ "#0055FF",
        storm_status == "SS" ~ "#6CE2FF",
        storm_status %in% c("LO", "WV", "DB") ~ "#A1A1A1",
        TRUE ~ "#FF5C00"
      )
    ) %>%
    arrange(datetime)

  storm_track$status_label <- factor(
    storm_track$status_label,
    levels = c(
      "Tropical Depression (TD)",
      "Tropical Storm (TS)",
      "Hurricane Category 1",
      "Hurricane Category 2",
      "Hurricane Category 3",
      "Hurricane Category 4",
      "Hurricane Category 5",
      "Extratropical Cyclone (EX)",
      "Subtropical Depression (SD)",
      "Subtropical Storm (SS)",
      "Low/Wave/Disturbance",
      "Missing Data"
    )
  )

  world <- ne_countries(scale = "medium", returnclass = "sf")
  states <- ne_states(returnclass = "sf")

  min_lon_range <- range(-100, -40)
  min_lat_range <- range(10, 50)

  #lon_range <- range(storm_track$lon)
  #lat_range <- range(storm_track$lat)
  #lon_padding <- diff(lon_range) * 0.1
  #lat_padding <- diff(lat_range) * 0.1

  p <- ggplot() +
    geom_sf(data = world, fill = "#E5E5E5", color = "#999999", size = 0.3) +
    geom_sf(data = states, fill = NA, color = "#CCCCCC", size = 0.2) +
    #coord_sf(
      #xlim = c(lon_range[1] - lon_padding, lon_range[2] + lon_padding),
      #ylim = c(lat_range[1] - lat_padding, lat_range[2] + lat_padding),
      #expand = FALSE
    #) +
    coord_sf(
      xlim = c(min_lon_range[1], min_lon_range[2]),
      ylim = c(min_lat_range[1], min_lat_range[2]),
      expand = FALSE
    ) +
    theme_minimal() +
    theme(
      panel.background = element_rect(fill = "#D4E6F1"),
      panel.grid.major = element_line(color = "#BBBBBB", size = 0.2),
      legend.position = "right",
      legend.key.size = unit(0.4, "cm"),
      legend.text = element_text(size = 8)
    ) +
    plot_theme

  if (nrow(storm_track) >= 2) {
    for (i in 1:(nrow(storm_track) - 1)) {
      segment_data <- storm_track[i:(i + 1), ]
      p <- p +
        geom_path(
          data = segment_data,
          aes(x = lon, y = lat),
          color = storm_track$line_color[i],
          size = 1
        )
    }
  }

  p <- p +
    geom_point(
      data = storm_track,
      aes(x = lon, y = lat, color = status_label),
      size = 1.5
    )

  landfall_data <- storm_track %>%
    filter(record_identifier == "L")

  if (nrow(landfall_data) > 0) {
    landfall_data <- landfall_data %>%
      mutate(rmw_degrees = rmw_meters / 111320)

    p <- p +
      geom_point(
        data = landfall_data,
        aes(x = lon, y = lat),
        color = landfall_data$line_color,
        size = 4,
        shape = 19
      ) +
      geom_point(
        data = landfall_data,
        aes(x = lon, y = lat),
        color = landfall_data$line_color,
        size = landfall_data$rmw_degrees * 100,
        shape = 1,
        stroke = 1.5,
        alpha = 0.5
      )
  }

  color_values <- c(
    "Tropical Depression (TD)" = "#2AFF00",
    "Tropical Storm (TS)" = "#FFD020",
    "Hurricane Category 1" = "#FF4343",
    "Hurricane Category 2" = "#FF6FFF",
    "Hurricane Category 3" = "#FF23D3",
    "Hurricane Category 4" = "#C916FF",
    "Hurricane Category 5" = "#FFFFFF",
    "Extratropical Cyclone (EX)" = "#202020",
    "Subtropical Depression (SD)" = "#0055FF",
    "Subtropical Storm (SS)" = "#6CE2FF",
    "Low/Wave/Disturbance" = "#A1A1A1",
    "Missing Data" = "#FF5C00"
  )

  p <- p +
    scale_color_manual(
      name = "Storm Status",
      values = color_values,
      breaks = names(color_values),
      drop = FALSE
    ) +
    labs(
      title = paste(params$storm_name, params$storm_year),
      x = "Longitude",
      y = "Latitude"
    )

  print(p)
```

## Track Observations

```{r track_table, echo=FALSE}
landfall_rows <- which(storm_track$record_identifier == "L")

storm_track %>%
  select(
    formatted_datetime,
    storm_status,
    record_identifier,
    lon,
    lat,
    rmw,
    pressure,
    windspeed
  ) %>%
  kable(
    col.names = c(
      "Date",
      "Status",
      "Record ID",
      "Lon",
      "Lat",
      "RMW",
      "Pressure",
      "Windspeed"
    ),
    align = c("l", "c", "c", "r", "r", "r", "r", "r"),
    booktabs = TRUE
  ) %>%
  kable_styling(
    #latex_options = c("striped"),
    font_size = 9,
    position = "center"
  ) %>%
  column_spec(
    2,
    color = "white",
    background = paste0(storm_track$line_color, "A0"),
    bold = T
  ) %>%
  row_spec(
    landfall_rows,
    color = "white",
    background = paste0(storm_track$line_color[landfall_rows], "A0")
  ) %>% 
  column_spec(
    6,
    color = "white",
    background = spec_color(storm_track$rmw, end = 0.7),
    bold = T
  ) %>% 
  column_spec(
    7,
    color = "white",
    background = spec_color(storm_track$pressure, end = 0.7),
    bold = T
  ) %>% 
  column_spec(
    8,
    color = "white",
    background = spec_color(storm_track$windspeed, end = 0.7),
    bold = T
  )

```

## Cost Normalization

```{r normalization, echo=FALSE}
#| out-width: "100%"

base_yearly_econ <- get_yearly_economics(storm$storm_year) 
base_yearly_usa <- get_yearly_usa_pop_hu(storm$storm_year)

base_yearly_metrics <- base_yearly_usa %>% 
  left_join(base_yearly_econ, by = "year") %>% 
  select(-aggregate_storm_loss)

ref_yearly_econ <- get_yearly_economics(latest_normalization_year)
ref_yearly_usa <- get_yearly_usa_pop_hu(latest_normalization_year)

ref_yearly_metrics <- ref_yearly_usa %>% 
  left_join(ref_yearly_econ, by = "year") %>% 
  select(-aggregate_storm_loss)

# year, ccn, gdp, ag loss
# base...
# ref...

# table cols metric, base year, ref year, ratio
# ccn
# gdp
# us pop
# us hous
# rwpc/hu

costs <- get_all_normalized_cost_index(storm)

# for each lf, create charts of mmh/mmp

unique_lf_ids <- unique_lfs$full_lf_id

for(lf in unique_lf_ids) {
  filtered_costs <- costs %>% 
    filter(
      full_lf_id == lf
    )
  
  mmh_plot <- ggplot(filtered_costs, aes(x = normalization_year, y = mmh_loss)) +
    geom_area(fill = "blue", alpha = 0.7) +
    scale_y_continuous(labels = scales::label_dollar(scale_cut = cut_short_scale())) +
    labs(
      title = paste(lf, "MMH Loss"),
      x = "Year",
      y = "Aggregate Loss (USD)"
    ) +
    theme_minimal() +
    plot_theme

  mmp_plot <- ggplot(filtered_costs, aes(x = normalization_year, y = mmp_loss)) +
    geom_area(fill = "red", alpha = 0.7) +
    scale_y_continuous(labels = scales::label_dollar(scale_cut = cut_short_scale())) +
    labs(
      title = paste(lf, "MMP Loss"),
      x = "Year",
      y = "Aggregate Loss (USD)"
    ) +
    theme_minimal() +
    plot_theme

  print(mmh_plot)
  print(mmp_plot)
}
```

## Population and Housing Growth

```{r growth, echo=FALSE}
#| out-width: "100%"

for(lf in unique_lf_ids) {
  growth <- get_normalized_metric_growth(storm, lf) %>%
    mutate(
      population_opacity = rescale(
        normalized_population,
        to = c(0.2, 0.8),
        from = range(normalized_population, na.rm = T)
      ),
      housing_opacity = rescale(
        normalized_housing,
        to = c(0.2, 0.8),
        from = range(normalized_housing, na.rm = T)
      )
    ) 

  growth <- get_county_and_state(growth)

  growth_sf <- growth %>%
    st_as_sf(wkt = "geom_wkt", crs = 4326) %>% 
    filter(year == max(year))

  combined_geom <- st_union(growth_sf)

  bbox <- st_bbox(combined_geom)

  lon_range <- c(bbox["xmin"], bbox["xmax"])
  lat_range <- c(bbox["ymin"], bbox["ymax"])

  padding_percent <- 0.1

  lon_padding <- diff(lon_range) * padding_percent
  lat_padding <- diff(lat_range) * padding_percent

  h_growth <- ggplot() +
    geom_sf(data = world, fill = "#E5E5E5", color = "#999999", size = 0.3) +
    geom_sf(data = states, fill = NA, color = "#CCCCCC", size = 0.2) +
    geom_sf(data = growth_sf, fill = "blue", color = "blue", alpha = growth_sf$housing_opacity) +
    geom_sf_label(data = growth_sf, aes(label = growth_sf$name)) +
    coord_sf(
      xlim = c(lon_range[1] - lon_padding, lon_range[2] + lon_padding),
      ylim = c(lat_range[1] - lat_padding, lat_range[2] + lat_padding),
      expand = FALSE
    ) +
    labs(
      x = "",
      y = ""
    )
    theme_minimal() +
    theme(
      panel.background = element_rect(fill = "#D4E6F1"),
      panel.grid.major = element_line(color = "#BBBBBB", size = 0.2),
      legend.position = "right",
      legend.key.size = unit(0.4, "cm"),
      legend.text = element_text(size = 8)
    ) +
    plot_theme

  latest_opacity <- growth %>%
    filter(year == max(year)) %>%
    select(county_fips, name, housing_opacity, population_opacity)

  growth_colored <- growth %>%
    left_join(latest_opacity, by = c("county_fips", "name"), suffix = c("", "_latest"))

  housing_colors <- latest_opacity %>%
    distinct(name, housing_opacity) %>%
    arrange(desc(housing_opacity)) %>%
    mutate(
      alpha_hex = sprintf("%02X", round(housing_opacity * 255)),
      color_with_alpha = paste0("#0000FF", alpha_hex)
    )

  housing_color_values <- setNames(housing_colors$color_with_alpha, housing_colors$name)

  growth_colored <- growth_colored %>%
    mutate(name = factor(name, levels = housing_colors$name))

  mmh_plot <- ggplot(growth_colored, aes(x = year, y = normalized_housing,
                                          color = name,
                                          group = county_fips)) +
    geom_line() +
    scale_color_manual(values = housing_color_values) +
    guides(color = guide_legend(override.aes = list(linewidth = 3))) +
    labs(
      title = paste(lf, "Housing Growth"),
      x = "Year",
      y = "Normalized Growth",
      color = "County"
    ) +
    theme_minimal() +
    plot_theme

  print(h_growth)
  print(mmh_plot)

  p_growth <- ggplot() +
    geom_sf(data = world, fill = "#E5E5E5", color = "#999999", size = 0.3) +
    geom_sf(data = states, fill = NA, color = "#CCCCCC", size = 0.2) +
    geom_sf(data = growth_sf, fill = "red", color = "red", alpha = growth_sf$population_opacity) +
    geom_sf_label(data = growth_sf, aes(label = growth_sf$name)) +
    coord_sf(
      xlim = c(lon_range[1] - lon_padding, lon_range[2] + lon_padding),
      ylim = c(lat_range[1] - lat_padding, lat_range[2] + lat_padding),
      expand = FALSE
    ) +
    labs(
      x = "",
      y = ""
    )
    theme_minimal() +
    theme(
      panel.background = element_rect(fill = "#D4E6F1"),
      panel.grid.major = element_line(color = "#BBBBBB", size = 0.2),
      legend.position = "right",
      legend.key.size = unit(0.4, "cm"),
      legend.text = element_text(size = 8)
    ) +
    plot_theme

  population_colors <- latest_opacity %>%
    distinct(name, population_opacity) %>%
    arrange(desc(population_opacity)) %>%
    mutate(
      alpha_hex = sprintf("%02X", round(population_opacity * 255)),
      color_with_alpha = paste0("#FF0000", alpha_hex)
    )

  population_color_values <- setNames(population_colors$color_with_alpha, population_colors$name)

  growth_colored <- growth_colored %>%
    mutate(name = factor(name, levels = population_colors$name))

  mmp_plot <- ggplot(growth_colored, aes(x = year, y = normalized_population,
                                          color = name,
                                          group = county_fips)) +
    geom_line() +
    scale_color_manual(values = population_color_values) +
    guides(color = guide_legend(override.aes = list(linewidth = 3))) +
    labs(
      title = paste(lf, "Population Growth"),
      x = "Year",
      y = "Normalized Growth",
      color = "County"
    ) +
    theme_minimal() +
    plot_theme

  print(p_growth)
  print(mmp_plot)
}
```