---
title: "`r paste(params$hurdat_id, params$storm_name, params$storm_year)` STORM REPORT"
subtitle: "Generated on: `r Sys.Date()`"
format:
  pdf:
    theme: cosmo
    keep-tex: false
    dev: png
    geometry:
      - top=0.75in
      - bottom=0.75in
      - left=0.75in
      - right=0.75in
      - footskip=0.3in
    embed-resourced: true
    fig-format: png
execute:
  warning: false
  message: false
  dpi: 100
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)
library(ggrepel)
library(ggtext)

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()$latest_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.text = element_text(size = 10)
    )

  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"
    ) +
    plot_theme

  print(p)
```

## Track Observations

```{r track_table, echo=FALSE}
landfall_rows <- which(storm_track$record_identifier == "L")
landfall_colors <- paste0(storm_track$line_color[landfall_rows], "A0")

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 = 11,
    position = "left"
  ) %>%
  column_spec(
    2,
    color = "white",
    background = paste0(storm_track$line_color, "A0"),
    bold = T
  ) %>%
  row_spec(
    landfall_rows,
    color = "white",
    background = landfall_colors
  ) %>% 
  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)

# year, ccn, gdp, pop, housing
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)

metrics_table <- data.frame(
  Metric = c("CCNSFACDG", "GDP Deflator", "US Population", "US Housing"),
  Base_Year = c(
    base_yearly_metrics$ccnsfacdg,
    base_yearly_metrics$gdp_deflator,
    base_yearly_metrics$population,
    base_yearly_metrics$housing
  ),
  Latest_Year = c(
    ref_yearly_metrics$ccnsfacdg,
    ref_yearly_metrics$gdp_deflator,
    ref_yearly_metrics$population,
    ref_yearly_metrics$housing
  )
) %>% 
  mutate(
    Ratio = Latest_Year / Base_Year
  )

metrics_table %>%
  kable(
    col.names = c(
      "Metric",
      paste("Base Year (", storm$storm_year, ")", sep = ""),
      paste("Ref Year (", latest_normalization_year, ")", sep = ""),
      "Ratio"
    ),
    align = c("l", "r", "r", "r"),
    booktabs = TRUE,
    digits = 2
  ) %>%
  kable_styling(
    font_size = 11,
    position = "left"
  )

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
    )

  last_year <- filtered_costs %>% filter(normalization_year == max(normalization_year))

  cost_plot <- ggplot(filtered_costs, aes(x = normalization_year)) +
    geom_line(aes(y = mmh_loss), color = "blue") +
    geom_line(aes(y = mmp_loss), color = "red") +
    geom_text(
      data = last_year,
      aes(x = normalization_year + 1, y = mmh_loss, label = "MMH"),
      hjust = 0,
      size = 3,
      color = "blue"
    ) +
    geom_text_repel(
      data = last_year,
      aes(x = normalization_year + 1, y = mmp_loss, label = "MMP"),
      hjust = 0,
      size = 3,
      color = "red"
    ) +
    scale_y_continuous(labels = scales::label_dollar(scale_cut = cut_short_scale())) +
    coord_cartesian(clip = "off") +
    labs(
      title = paste(lf, "MMH and MMP Loss"),
      x = "Year",
      y = "Aggregate Loss (USD)"
    ) +
    theme_minimal() +
    plot_theme +
    theme(plot.margin = margin(5.5, 40, 5.5, 5.5))
  
  #mmh_plot <- ggplot(filtered_costs, aes(x = normalization_year, y = mmh_loss)) +
  #  geom_line(color = "blue") +
  #  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_line(color = "red") +
  #  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(cost_plot)
}
```

## Population and Housing Growth

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

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"])

  lon_width <- diff(lon_range)
  lat_height <- diff(lat_range)

  target_ratio <- 16/9
  current_ratio <- lon_width / lat_height

  if (current_ratio < target_ratio) {
    needed_width <- lat_height * target_ratio
    extra_width <- (needed_width - lon_width) / 2
    lon_range <- c(lon_range[1] - extra_width, lon_range[2] + extra_width)
  }

  if (current_ratio > target_ratio) {
    needed_height <- lon_width / target_ratio
    extra_height <- (needed_height - lat_height) / 2
    lat_range <- c(lat_range[1] - extra_height, lat_range[2] + extra_height)
  }

  padding_percent <- 0.1

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

  base_growth_map <- 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_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() +
    plot_theme

  h_growth_map <- base_growth_map +
    geom_sf(data = growth_sf, fill = "blue", color = "#e5e5e5", size = 0.5, alpha = growth_sf$housing_opacity) +
    geom_sf_label(data = growth_sf, aes(label = growth_sf$name)) +
    labs(
      title = paste(lf, "Housing Growth Map By County"),
    ) +
    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 
    ) 

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

  h_plot <- ggplot(growth_colored, aes(x = year, y = normalized_housing,
                                          color = name,
                                          group = county_fips)) +
    geom_line() +
    scale_color_manual(values = housing_color_values) +
    geom_text_repel(
      data = growth_colored %>% filter(year == max(year)),
      aes(label = name, x = year + 1),
      direction = "y",
      hjust = 0,
      size = 3
    ) +
    coord_cartesian(clip = "off") +
    guides(color = "none") +
    labs(
      title = paste(lf, "Housing Growth"),
      x = "Year",
      y = "Normalized Growth"
    ) +
    theme_minimal() +
    plot_theme +
    theme(plot.margin = margin(5.5, 40, 5.5, 5.5))

  print(h_growth_map)
  print(h_plot)

  p_growth_map <- base_growth_map +
    geom_sf(data = growth_sf, fill = "red", color = "#e5e5e5", size = 0.5, alpha = growth_sf$population_opacity) +
    geom_sf_label(data = growth_sf, aes(label = growth_sf$name)) +
    labs(
      title = paste(lf, "Population Growth Map By County"),
    ) + 
    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 
    ) 

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

  p_plot <- ggplot(growth_colored, aes(x = year, y = normalized_population,
                                          color = name,
                                          group = county_fips)) +
    geom_line() +
    scale_color_manual(values = population_color_values) +
    geom_text_repel(
      data = growth_colored %>% filter(year == max(year)),
      aes(label = name, x = year + 1),
      direction = "y",
      hjust = 0,
      size = 3
    ) +
    coord_cartesian(clip = "off") +
    guides(color = "none") +
    labs(
      title = paste(lf, "Population Growth"),
      x = "Year",
      y = "Normalized Growth"
    ) +
    theme_minimal() +
    plot_theme +
    theme(plot.margin = margin(5.5, 40, 5.5, 5.5))


  print(p_growth_map)
  print(p_plot)
}
```