vulnerability_index_calculation.Rmd

---
title: "NPS Water Supply Vulnerability Index Calculation"
author: "Caitlin Mothes"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    theme: flatly
    code_folding: show
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

source('setup.R')
```

# Overview

This workflow calculates a **Total Vulnerability Index** for each NPS water supply system following the hierarchical Euclidean distance framework from [Michalak et al. (2021)](https://irma.nps.gov/DataStore/Reference/Profile/2287636).

**Calculation steps at every level:**

1.  Min-max normalize each indicator to 0–1 (higher = more vulnerable)
2.  Aggregate indicators within a group via Euclidean distance: $\sqrt{\sum x_i^2}$
3.  Re-normalize the aggregated score to 0–1
4.  Pass result to the next hierarchical level

**Hierarchy:** Indicators → Vulnerability Factors → Components (Exposure / Sensitivity) → Total Vulnerability Index (Exposure + Sensitivity)

## Indicator Framework Table:

|  |  |  |  |  |  |
|------------|------------|------------|------------|------------|------------|
| COMPONENT | FACTOR | INDICATOR | description | data_source | function_name |
| Exposure | Runoff | Change in runoff | \% change in mean annual runoff |  | calc_runoff_exposure.R |
| Runoff model agreement | \% gcms predicting decrease in runoff |  | calc_runoff_exposure.R |  |  |
| Precipitation | Change in precipitation | \% change in mean annual precipitation |  | calc_precip_exposure.R |  |
| Precipitation model agreement | \% models predicting decrease in precipitation |  | calc_precip_exposure.R |  |  |
| Flood | Change in flood risk |  |  | TBD |  |
| Sea Level Rise | Inundation from sea level rise | \% area inundated |  | calc_inundation_exposure.R |  |
| Inundation from storm surge |  |  | TBD |  |  |
| Projected saltwater intrusion |  |  | TBD |  |  |
| Wildfire | Change in fire probability | \% change in probability of wildfire |  | calc_fire_exposure.R |  |
| Sensitivity | Demand | Historic visitation trend |  |  | calc_visitation_sensitivity.R |
| Competition | trend in water use for water supply county |  | calc_nearby_use_sensitivity.R |  |  |
| Water Supply | Water supply type |  |  | WSD |  |
| Treatment | Water supply treatment type |  |  | WSD |  |
| Wildfire | Current wildfire risk |  |  | calc_fire_sensitivity.R |  |
| Flood | Current flood risk | \% area in high risk flood zone | FEMA flood maps | calc_flood_sensitivity.R |  |
| Sea Level Rise | Distance to coast | water supply euclidean distance to nearest coastline |  | TBD |  |
| Current inundation | \% area currently inundated |  | calc_inundation_exposure.R |  |  |
| Runoff | Historic runoff | historic 30 year trend in runoff |  | calc_runoff_sensitivity.R |  |
| Precipitation | Historic precipitation | historic 40 year trend in precipitation |  | calc_precip_sensitivity.R |  |

# 1. Helper Functions

```{r helper-functions}
#' Min-max normalize a vector to [0, 1]
#' @param x numeric vector
#' @param na.rm logical; remove NAs before computing min/max
minmax_norm <- function(x, na.rm = TRUE) {
  rng <- range(x, na.rm = na.rm)
  if (rng[2] == rng[1]) return(rep(0, length(x)))
  (x - rng[1]) / (rng[2] - rng[1])
}


#' Euclidean distance aggregation (row-wise)
#' @param df data.frame of normalized indicator columns (one row per water supply)
#' @return numeric vector of Euclidean distances
euclid_agg <- function(df) {
  sqrt(rowSums(df^2, na.rm = TRUE))
}


#' Full aggregation pipeline: Euclidean distance + re-normalize
#' @param df data.frame with only the indicator columns to aggregate
#' @return numeric vector, rescaled 0–1
aggregate_indicators <- function(df) {
  raw <- euclid_agg(df)
  minmax_norm(raw)
}
```

# 2. Assemble Indicators

Pull the most recent output file for each indicator/component combination listed in the indicators table above. Files with extra leading or trailing text in the name (e.g., `kw_`, `moderate`) are excluded by filtering to exact stem matches.

```{r assemble-indicators}
# ── Valid stems derived from the indicators table (excludes TBD / WSD rows) ──
valid_stems <- c(
  # Exposure
  "fire_exposure",
  "inundation_exposure",
  "swe_exposure",
  "precipitation_exposure",
  "runoff_exposure",
  # Sensitivity
  #"competition_sensitivity",
  "fire_sensitivity",
  "flood_sensitivity",
  "nearby_use_sensitivity",
  "precip_sensitivity",
  "runoff_sensitivity",
  "visitation_sensitivity",
  "swe_sensitivity",
  "inundation_sensitivity"
)

# ── Find the most recent file for each valid stem ────────────────────────────
latest_files <- list.files("data/rapid/outputs", pattern = "\\.csv$", full.names = TRUE) |>
  tibble(file = _) |>
  mutate(
    stem = str_remove(basename(file), "_\\d{4}-\\d{2}-\\d{2}\\.csv$"),
    date = as.Date(str_extract(basename(file), "\\d{4}-\\d{2}-\\d{2}"))
  ) |>
  filter(stem %in% valid_stems) |>
  slice_max(date, by = stem)

# ── Read and full-join all indicators by wsd_source_id ──────────────────────
indicators_assembled <- latest_files |>
  select(stem, file) |>
  deframe() |>
  imap(\(f, nm) {
    df <- read_csv(f, show_col_types = FALSE) %>% 
      select(wsd_source_id, starts_with("raw_"), starts_with("norm_"))
    # this was pre-filtered in rapid_workflow
    #if (nm == "inundation_exposure") df <- filter(df, buffer_km == 4)
    # if (nm == "runoff_exposure") df <- rename(df, raw_runoff_models_showing_decrease_pct = raw_models_showing_decrease_pct)
    # if(nm == "precipitation_exposure") df <- rename(df, raw_precipitation_models_showing_decrease_pct = raw_models_showing_decrease_pct,                                        norm_precipitation_models_showing_decrease_pct = norm_models_showing_decrease_pct)
    df
  }) |>
  reduce(full_join, by = "wsd_source_id")


```

```{r}
# Save the file if updated
write_csv(indicators_assembled, "data/rapid/indicators_assembled.csv")
```

# 3. Prep Indicators

The expected input is a single data frame with **one row per water supply system** and columns for each indicator.

```{r}
# ── Select & rename to standardized indicator columns 
# All norm_ columns are already min-max normalized 0–1.
# For multi-metric indicators, choose the appropriate percentile:
#   - Runoff/precip/swe exposure: p10 (largest decrease = most vulnerable)
#   - Other exposure:         p90 (largest increase = most vulnerable)

# read in assembled indicators
data <- read_csv("data/rapid/indicators_assembled.csv")

# clean
indicators_clean <- data %>%
  transmute(
    # ── ID columns ──
    wsd_source_id,
    
    # EXPOSURE INDICATORS
    
    # Runoff
    #exp_norm_runoff_change = norm_runoff_change_p10,
    exp_runoff_change = raw_runoff_change_p10,
    #exp_norm_runoff_model_agree = norm_runoff_models_showing_decrease,
    exp_runoff_model_agree = raw_runoff_models_showing_decrease,
    
    # Precipitation
    #exp_precip_change = norm_precipitation_change_p10,
    exp_precip_change = raw_precipitation_change_p10,
    #exp_precip_model_agree = norm_precip_models_showing_decrease,
    exp_precip_model_agree = raw_precip_models_showing_decrease,
    
    # SWE
    #exp_swe_change = norm_swe_change_p10,
    exp_swe_change = raw_swe_change_p10,
    #exp_swe_model_agree = norm_swe_models_showing_decrease,
    exp_swe_model_agree = raw_swe_models_showing_decrease,
    
    # Sea Level Rise
    ## Inundation
    #exp_inundation_slr = norm_pct_point_change,
    exp_inundation_slr = raw_pct_point_change,
    
    # Fire
    #exp_fire_prob_change = norm_fire_prob_change_p90,
    exp_fire_prob_change = raw_fire_prob_change_p90,
    
    # SENSITIVITY INDICATORS
    
    # Visitation
    #sen_visitation_trend = norm_scaled_trend_visitation,
    sen_visitation_trend = raw_scaled_trend_visitation,
    
    # Competition
    #sen_competition = norm_slope_nearby_use_per_km2,
    sen_competition = raw_slope_nearby_use_per_km2,
    
    # Placeholders (no raw to add yet)
    sen_water_supply_type = NA_real_,
    sen_treatment_type  = NA_real_,
    
    # Wildfire
    #sen_wildfire_hazard = norm_wildfire_hazard_mean,
    sen_wildfire_hazard = raw_wildfire_hazard_mean,
    
    # Flood
    #sen_flood_risk = norm_flood_risk_percent,
    sen_flood_risk = raw_flood_risk_percent,
    
    # Current inundation
    #sen_inundation_current = norm_coverage_pct_ref,
    sen_inundation_current = raw_coverage_pct_ref,
    
    # Historic runoff
    #sen_runoff_trend = norm_slope_runoff,
    sen_runoff_trend = raw_slope_runoff,
    
    # Historic precipitation
    #sen_precip_trend  = norm_slope_precip,
    sen_precip_trend  = raw_slope_precip,
    
    # Historic SWE
    sen_swe_trend = raw_slope_swe
  ) %>% 
  distinct(wsd_source_id, .keep_all = TRUE)

# ── Summary check ───────────────────────────────────────────────────────
indicators_clean %>%
  distinct(wsd_source_id, .keep_all = TRUE) %>% 
  summarise(across(everything(), ~sum(!is.na(.)))) %>%
  tidyr::pivot_longer(everything(), names_to = "indicator", values_to = "n_valid") %>%
  print(n = 30)


```

# SAVE for app

Save `indicators_clean`, this is the input of raw vars that is filtered and fed to calc_vulnerability_index in the rapid app

```{r}
write_csv(indicators_clean, "rapid_app/app_data/final_indicators.csv")
```

# 4. Calculate Vulnerability Index

Test calculating the vulnerability index

```{r}
# Run the vulnerability index calculation function 
results <- calc_vulnerability_index(
  indicators_clean,
  id_cols = c("wsd_source_id", "park_unit", "park_name")
)

#write_csv(results, "data/rapid/final_index.csv")
#write_csv(results, "rapid_app/app_data/final_index.csv")
```

# Reference

**Framework:** Michalak, J.L., Lawler, J.J., Gross, J.E. & Littlefield, C.E. (2021). *A Strategic Analysis of Climate Vulnerability of National Park Resources and Values.* NPS Natural Resource Report NPS/NRSS/CCRP/NRR—2021/2293.

**Key methodological notes:**

-   All indicators standardized via min-max normalization (0 = least vulnerable, 1 = most vulnerable)
-   Aggregation at every level uses Euclidean distance: $\sqrt{\sum x_i^2}$
-   Aggregated scores are re-normalized to 0–1 before passing to next level
-   No weighting applied (equal weight to all factors); Michalak et al. found weighting had minimal effect on total scores
-   Final total vulnerability score is **not** re-normalized (can exceed 1) per Michalak et al., but we rescale for interpretability since our framework has only 2 components