NPS_WBM_Comparison_Report.Rmd

---
output: pdf_document
---

```{r include=FALSE}
knitr::opts_chunk$set(fig.height = 10, fig.width = 8, echo = FALSE, message = FALSE, warning = FALSE)
```

::: {.center data-latex=""}
\begin{center}
\fontsize{20pt}{28pt}\selectfont NPS Water Balance Model Comparison Report
\end{center}
:::

<br>

------------------------------------------------------------------------

# *Comparison to other CONUS-wide hydrologic component estimates (Saxe et al. 2021)*

\vspace{0.5in}

Caitlin Mothes^1^, Kristen Cognac^2^, Katie Willi^2^

\vspace{0.75in}

^1^ Geospatial Centroid, CSU Libraries, Colorado State University, Fort Collins, Colorado

*Corresponding Author: ccmothes\@colostate.edu*



^2^ Radical Open Science Syndicate (ROSS) Lab, Department of Ecosystem Science and Sustainability, Colorado State University, Fort Collins, Colorado

\vspace{3.5in}

December 2024

\newpage

## Background

Saxe et al. 2021 compiled estimates of precipitation, actual evapotranspiration, runoff, snow water equivalent and soil moisture from published datasets representing 47 different hydrologic models, reanalysis data sets and remote sensing products across CONUS. This study produced an open-source dataset of monthly values for the above mentioned variables for each of the 47 models averaged across each U.S. Level 1 ecoregions (Omernik and Griffith 2014) from 1900-2018 (though data availability varies by model).

## Methods

We conducted an analysis of NPS WBM products for Runoff, AET and SWE following the aggregated area-weighted mean methods in Saxe et al. 2021 to calculate monthly values for each ecoregion. For annual average comparisons, values for Runoff and AET were summed, while SWE was averaged. Below we report the comparisons between the NPS WBM values and other models visually, while also highlighting other hydrologic models of interest and ecoregions of interest based on the target park unit for this report.The other models of interest include **NLDAS2-Noah, NLDAS2-VIC, and Livneh-VIC.**



## Ecoregions

![Map of Level 1 ecoregions (*from Saxe et al. 2021*).](www/ecoregions_saxe_etal_2021.png)


<br>


### Data Visualization Tool
You can also interactively view the data presented in this report here: https://geocentroid.shinyapps.io/nps_wbm_comparison/



\newpage

# **Findings**

```{r setup, include=FALSE}
source("setup.R")

################################## INTERNAL NOTE ##################################
# All DATA is in the data/Saxe_etal folder. CODE used to generate NPS WBM matching data is in the saxe_comparison_job.R script

# READ IN DATASETS

## NPS WBM generated datasets
wbm_eco_means <- purrr::map_df(list.files("data/Saxe_etal", pattern = "eco_mean", full.names = TRUE), read_csv)

## Saxe et al. dataset
saxe_data <- read_csv("data/Saxe_etal/Ecoregion Aggregates.csv")

## set regions of interest for specific park/NPS unit
roi <- c("NORTH AMERICAN DESERTS", "NORTHWESTERN FORESTED MOUNTAINS")


# JOIN DATASETS

annual_combined <- saxe_data %>% 
  #filter data to match wbm
  filter(ecoregion %in% wbm_eco_means$ecoregion,
         date %in% wbm_eco_means$date) %>% 
  bind_rows(wbm_eco_means) %>% 
  # summarize annual data for comparison
  mutate(year = year(date)) %>% 
  pivot_wider(names_from = component, values_from = value) %>% 
  group_by(model, ecoregion, year) %>% 
  dplyr::summarize(across(AET:R, \(x) sum(x, na.rm = TRUE), .names = "{col}_annual"),
                   SWE_annual = mean(SWE, na.rm = TRUE)) %>% 
  mutate_all(function(x) ifelse(is.nan(x), NA, x)) %>% 
  # drop NAs across all vars
  #drop_na(AET_annual, R_annual, P_annual, SWE_annual) %>% 
  # add color column for plots
  mutate(color = case_when(model == "NPS_WBM"~ "Mod1", 
                           model %in% c("NLDAS2-Noah", "NLDAS2-VIC", "Livneh-VIC") ~ "Mod2",
                           .default = "Mod3"))


monthly_combined <- saxe_data %>% 
  #filter data to match wbm
  filter(ecoregion %in% wbm_eco_means$ecoregion,
         date %in% wbm_eco_means$date) %>% 
  bind_rows(wbm_eco_means) %>% 
  # add color column for plots
  mutate(color = case_when(model == "NPS_WBM"~ "Mod1", 
                           model %in% c("NLDAS2-Noah", "NLDAS2-VIC", "Livneh-VIC") ~ "Mod2",
                           .default = "Mod3"))


```

```{r figure-functions}
# annual averages boxplot
boxplot <- function(data, variable, variable_name) {
  data %>%
    filter({{variable}} != 0) %>%
    drop_na({{variable}}) %>% 
    ggplot() +
    geom_boxplot(aes(
      x = {{variable}},
      y = reorder(model, {{variable}}, fun = median),
      fill = color
    )) +
    scale_fill_manual(
      values = c("red", "royalblue1", "white"),
      labels = c("NPS WBM", "Models of Interest", "Other Models")
    ) +
    facet_wrap( ~ ecoregion, scales = "free_x", ncol = 3) +
    labs(
      y = "",
      title = variable_name,
      x = paste0("Annual ", variable_name, " (mm)")
    ) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      legend.position = "bottom",
      axis.text.y = element_text(size = 8),
      title = element_text(size = 14),
      text = element_text(family = "serif")
    )
}


# time series function
timeseries <- function(data,
                       regions = unique(saxe_data$ecoregion),
                       ncol,
                       variable,
                       variable_name) {
  data %>%
    filter(component == variable, ecoregion %in% regions) %>%
    ggplot() +
    geom_line(aes(
      x = date,
      y = value,
      color = factor(color, levels = c("Mod3", "Mod2", "Mod1"))
    ), alpha = 0.5) +
    scale_color_manual(
      values = c("gray", "blue", "red"),
      labels = c("Other Models", "Models of Interest", "NPS WBM")
    ) +
    facet_wrap( ~ ecoregion, 
                scales = "free_y", 
                ncol = ncol) +
    labs(
      title = paste0("Monthly ", variable_name),
      y =  paste0("Total ", variable_name, " (mm)")
    ) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      legend.position = "bottom",
      text = element_text(family = "serif")
    )
  
}

```

## Annual Averages

```{r}
# Runoff
boxplot(annual_combined, R_annual, "Runoff")
```

```{r}
# AET
boxplot(annual_combined, AET_annual, "AET")
```

```{r}
# SWE
boxplot(annual_combined, SWE_annual, "SWE")
```

<br>

## Monthly Time-Series

```{r}
# Runoff all regions
timeseries(monthly_combined, ncol = 2, variable = "R", variable_name = "Runoff")

```


```{r eval = FALSE}
# Runoff regions of interest
timeseries(
  monthly_combined,
  ncol = 1,
  regions = roi,
  variable = "R",
  variable_name = "Runoff"
)
```




```{r}
# AET all regions
timeseries(monthly_combined, ncol = 2, variable = "AET", variable_name = "AET")

```

```{r eval=FALSE}
# AET regions of interest
timeseries(
  monthly_combined,
  ncol = 1,
  regions = roi,
  variable = "AET",
  variable_name = "AET"
)
```





```{r}
# SWE all regions
timeseries(monthly_combined, ncol = 2, variable = "SWE", variable_name = "SWE")

```

```{r eval = FALSE}
# AET regions of interest
timeseries(
  monthly_combined,
  ncol = 1,
  regions = roi,
  variable = "SWE",
  variable_name = "SWE"
)
```

<br>

\newpage

# Literature Cited

Omernik, J. M. and Griffith, G. E.: Ecoregions of the conterminous United States: Evolution of a hierarchical spatial framework, Environ. Manage., 54, 1249–1266, <https://doi.org/10.1007/s00267-014-0364-1>, 2014.

Saxe, S., Farmer, W., Driscoll, J., and Hogue, T. S.: Implications of model selection: a comparison of publicly available, conterminous US-extent hydrologic component estimates, Hydrol. Earth Syst. Sci., 25, 1529–1568, <https://doi.org/10.5194/hess-25-1529-2021>, 2021.