analysis_notebook.Rmd

---
title: Near-term forecasts of NEON lakes reveal gradients of environmental predictability
  across the U.S.
author:
  - R. Quinn Thomas
  - Ryan P. McClure
  - Tadhg N. Moore
  - Whitney M. Woelmer
  - Carl Boettiger
  - Renato J. Figueiredo
  - Robert T. Hensley
  - Cayelan C. Carey
date: "`r Sys.Date()`"
output:
  word_document: default
  html_notebook: default
  pdf_document: default
---

# Set up analysis

The analysis needs the following libraries, environment variables, and sourced R functions.

```{r}
library(tidyverse)
library(lubridate)
library(patchwork)
library(egg)
library(broom)
library(arrow)
library(dplyr)
Sys.setenv("AWS_EC2_METADATA_DISABLED"="TRUE")
Sys.unsetenv("AWS_ACCESS_KEY_ID")
Sys.unsetenv("AWS_SECRET_ACCESS_KEY")
Sys.unsetenv("AWS_DEFAULT_REGION")
Sys.unsetenv("AWS_S3_ENDPOINT")
lake_directory <- here::here()
use_archive <- TRUE
source(file.path(lake_directory,"workflows","neon_lakes_ms","scoring.R"))
dir.create("notebook", showWarnings = FALSE)
```

# Generation of forecasts

The analysis below assumes that forecasts have already been generated using `workflows/neon_lakes_ms/01_combined_paper_workflow.R`.  Forecasts can be re-generated by running the script and setting `use_archive <- FALSE` on line 16.  Forecasts can be also be download from Zenado at the following DOI (Pending)

# Reading in scores

The scores (i.e., metrics of how well a forecast predicts an observation) are already calculated from the forecasts using the the `workflows/neon_lakes_ms/02_score_forecasts.R` script. The following code reads the individual forecast files into memory and stacks the the scored forecasts together. The scores are located in a remote s3 bucket and are read into local memory without needing to download the files. The filtering to select only the "teams" of interest (i.e., models run) occurs remotely, thus reducing the size of the data that are transferred.

```{r}
s <- score_schema()
if(!use_archive){
  s3 <- arrow::s3_bucket(bucket = "scores",
                         endpoint_override = "s3.flare-forecast.org",
                         anonymous=TRUE)
}else{
  dir.create("scores")
  obj <- contentid::resolve("hash://md5/786ccc31d5721656e938dba409c656b1", store=TRUE)
  unzip(obj,exdir = "scores")
  s3 <- arrow::SubTreeFileSystem$create("scores")
}

ds <- arrow::open_dataset(s3, schema=s, format = "csv", skip_rows = 1)
combined <- ds %>% filter((team == "ms2_glm_flare" | team == "ms2_doymean") & 
                          time < as_date("2021-10-24")) %>% collect()
```

# Clean-up up table of scores

The combined table is cleaned up by renaming of the columns, defining the NEON ecoclimatic domains, expanding the model names, and defining the start of "fall" that is used in later analyses. The squared error is calculated for each forecast-observation pair.

```{r}
combined1 <- combined %>%
  rename("siteID" = theme) %>%
  mutate(resid = observed - mean,
         sq_error = (resid)^2) %>%
  select(siteID, team, time, forecast_start_time, issue_date, horizon, depth, mean, observed,sq_error, crps) %>%
  rename(Model = team) %>%
  mutate(Model = ifelse(Model == "ms2_doymean", "Day-of-year mean", Model),
         Model = ifelse(Model == "ms2_glm_flare","FLARE-GLM", Model)) %>%
  mutate(Model = factor(Model, levels = c("FLARE-GLM","Day-of-year mean")))
```

```{r}
combined2 <- combined1 %>% 
  select(siteID, Model, time, forecast_start_time, issue_date, horizon, depth, sq_error) %>%
  pivot_wider(names_from = Model, values_from = sq_error) %>%
  na.omit() %>%
  select(siteID, time, horizon, depth, `Day-of-year mean`, `FLARE-GLM`) %>% 
  pivot_longer(cols = -c("siteID","time", "horizon", "depth"), names_to = "Model",values_to = "sq_error") %>%
  pivot_wider(names_from = horizon, values_from = sq_error) %>% 
  na.omit() %>% 
  pivot_longer(cols = -c("siteID","time", "Model", "depth"), names_to = "horizon", values_to = "sq_error")

combined_top <- combined2 %>% 
  filter(depth == 0.0) %>% 
  group_by(Model, siteID, horizon) %>%  # average over siteID
  summarise(rmse = mean(sq_error, na.rm =TRUE),
            .groups = "drop") %>%
  mutate(rmse = sqrt(rmse)) %>% 
  mutate(domain = NA,
         domain = ifelse(siteID %in% c("BARC","SUGG"), "Southeast", domain),
         domain = ifelse(siteID %in% c("CRAM","LIRO"), "Great Lakes", domain),
         domain = ifelse(siteID %in% c("PRLA","PRPO"), "Northern Plains", domain),
         domain = factor(domain, levels = c("Northern Plains","Great Lakes","Southeast")),
         siteID = factor(siteID, levels = c("PRLA","PRPO","CRAM","LIRO","BARC","SUGG")))

focal_depth_top <- 0.0
focal_depth_bottom <- 1.0
combined_1m <- combined2 %>% 
  filter(depth >= focal_depth_top,
         depth <= focal_depth_bottom) %>% 
  group_by(Model, siteID, horizon) %>%  # average over siteID
  summarise(rmse = mean(sq_error, na.rm =TRUE),
            .groups = "drop") %>%
  mutate(rmse = sqrt(rmse)) %>% 
  mutate(domain = NA,
         domain = ifelse(siteID %in% c("BARC","SUGG"), "Southeast", domain),
         domain = ifelse(siteID %in% c("CRAM","LIRO"), "Great Lakes", domain),
         domain = ifelse(siteID %in% c("PRLA","PRPO"), "Northern Plains", domain),
         domain = factor(domain, levels = c("Northern Plains","Great Lakes","Southeast")),
         siteID = factor(siteID, levels = c("PRLA","PRPO","CRAM","LIRO","BARC","SUGG")))

combined_1m_all_sites <- combined2 %>% 
  filter(depth >= focal_depth_top,
         depth <= focal_depth_bottom) %>% 
  group_by(Model, horizon) %>%  # average over siteID
  summarise(rmse = mean(sq_error, na.rm =TRUE),
            .groups = "drop") %>%
  mutate(rmse = sqrt(rmse))

combined_all_depths <- combined2 %>% 
  group_by(Model, siteID, horizon) %>%  # average over siteID
  summarise(rmse = mean(sq_error, na.rm =TRUE),
            .groups = "drop") %>%
  mutate(rmse = sqrt(rmse)) %>% 
  mutate(domain = NA,
         domain = ifelse(siteID %in% c("BARC","SUGG"), "Southeast", domain),
         domain = ifelse(siteID %in% c("CRAM","LIRO"), "Great Lakes", domain),
         domain = ifelse(siteID %in% c("PRLA","PRPO"), "Northern Plains", domain),
         domain = factor(domain, levels = c("Northern Plains","Great Lakes","Southeast")),
         siteID = factor(siteID, levels = c("PRLA","PRPO","CRAM","LIRO","BARC","SUGG")))

combined_1m_crps <- combined1 %>% 
  select(siteID, Model, time, forecast_start_time, issue_date, horizon, depth, crps) %>%
  pivot_wider(names_from = Model, values_from = crps) %>%
  na.omit() %>%
  select(siteID, time, horizon, depth, `Day-of-year mean`, `FLARE-GLM`) %>% 
  pivot_longer(cols = -c("siteID","time", "horizon", "depth"), names_to = "Model",values_to = "crps") %>%
  pivot_wider(names_from = horizon, values_from = crps) %>% 
  na.omit() %>% 
  pivot_longer(cols = -c("siteID","time", "Model", "depth"), names_to = "horizon", values_to = "crps") %>% 
  filter(depth >= focal_depth_top,
         depth <= focal_depth_bottom) %>% 
  group_by(Model, siteID, horizon) %>%  # average over siteID
  summarise(crps = mean(crps, na.rm =TRUE),
            .groups = "drop") %>%
  mutate(domain = NA,
         domain = ifelse(siteID %in% c("BARC","SUGG"), "Southeast", domain),
         domain = ifelse(siteID %in% c("CRAM","LIRO"), "Great Lakes", domain),
         domain = ifelse(siteID %in% c("PRLA","PRPO"), "Northern Plains", domain),
         domain = factor(domain, levels = c("Northern Plains","Great Lakes","Southeast")),
         siteID = factor(siteID, levels = c("PRLA","PRPO","CRAM","LIRO","BARC","SUGG")))
```

# Building derived tables

The analysis uses a set of metric that summarize each 35-day forecast at a specific depth (the top 1m). The first is the RMSE at day 1. We choice the 1-day horizon because weather forecasts are best one-day out so it isolates the model error. The second is different between the max and min RMSE. This represents the magnitude that the forecast degrades.

Read in lake characteristic data

```{r}
lake_info <- read_csv(file.path(lake_directory,"workflows","neon_lakes_ms", "LakeTable_15Nov2021.csv"))

variables <- tibble(variable = c("max_depth_m", "fetch_m", "volume_m3", "surface_area_km2", "mean_secchi_depth_m",
         "mean_annual_temperature_degreeC", "mean_annual_precipitation_mm", "air_temp_10day_sd_degreeC",
         "mean_hydraulic_residence_time_yrs", "catchment_size_km2", "latitude"),
         full_name = c("Maximum\ndepth", "Fetch", "Volume", "Surface\narea","Water clarity\n(Secchi depth)","Mean annual air\ntemperature", "Mean annual\nprecipitation","10-day variance in\n air temperature","Mean hydrologic\nresidence time","Catchment\nsize", "Latitude"),
         full_name_1line = c("Maximum depth", "Fetch", "Volume", "Surface area","Water clarity (Secchi depth)","Mean annual air temperature", "Mean annual precipitation","10-day variance in air temperature","Mean hydrologic residence time","Catchment size", "Latitude") )
```

## RMSE at 1-day ahead

```{r}
day1_max_min <- combined_1m %>% 
  filter(Model == "FLARE-GLM",
         horizon == 1) %>% 
  mutate(horizon = as.character(horizon)) %>% 
  select(Model, siteID, horizon, rmse)
```

## Max - Min RMSE over all horizons

```{r}
max_min35day <- combined_1m %>% 
  filter(Model == "FLARE-GLM") %>% 
  group_by(siteID, Model) %>% 
  summarize(rmse = max(rmse) - min(rmse), .groups = "drop") %>% 
  mutate(horizon = "max_min") %>% 
  select(Model,siteID, horizon, rmse)
```

## Combine tables together

```{r}

day1_max_min <- bind_rows(day1_max_min, max_min35day) %>% 
left_join(lake_info, by = "siteID") %>% 
   mutate(domain = NA,
         domain = ifelse(siteID %in% c("BARC","SUGG"), "Southeast", domain),
         domain = ifelse(siteID %in% c("CRAM","LIRO"), "Great Lakes", domain),
         domain = ifelse(siteID %in% c("PRLA","PRPO"), "Northern Plains", domain),
         domain = factor(domain, levels = c("Northern Plains","Great Lakes","Southeast")),
         siteID = factor(siteID, levels = c("PRLA","PRPO","CRAM","LIRO","BARC","SUGG")))
```

## Calculate correlations with lake characteristics

Correlations with 1-day ahead RMSE and max - min RMSE. 1-day ahead RMSE is multiplied by -1 so that smaller values equal less accuracy. This allows a positive correlation to mean that accuracy increases with the lake characteristic of interest

Calculate Spearmean's Rho

```{r}
spear_corr_table <- day1_max_min %>% 
  select(rmse, horizon, max_depth_m, fetch_m, volume_m3, surface_area_km2, mean_secchi_depth_m,
         mean_annual_temperature_degreeC, mean_annual_precipitation_mm, air_temp_10day_sd_degreeC,
         mean_hydraulic_residence_time_yrs, catchment_size_km2, latitude) %>% 
  pivot_longer(cols = -c("rmse","horizon"), names_to = "variable", values_to = "value") %>% 
  filter(horizon %in% c("1", "max_min")) %>% 
  group_by(horizon, variable) %>% 
  summarize(rho = cor(rmse, value, method = "spearman"), 
            .groups = "drop") %>% 
    mutate(Metric = ifelse(horizon == "1",  "accuracy", "accuracy\ndegradation"),
           rho = ifelse(Metric == "accuracy", -rho, rho)) %>%
   arrange(desc(abs(rho))) %>% 
    left_join(variables, by = "variable")

corr_names <- rev(as.character(unique(spear_corr_table$full_name)))

spear_corr_table <- spear_corr_table %>%
  mutate(full_name = factor(full_name, levels = corr_names))
```

# Generate Figures

## Figure 1

Figure 1 is the map.

## Figure 2

```{r}
y_label <- expression(paste('RMSE (',degree,'C)', sep = ""))
cbPalette <- c("#a50026", "#f46d43", "#fdae61", "#fee090","#74add1", "#313695")
  
Fig2a <- combined_1m %>% 
  filter(Model == "FLARE-GLM") %>% 
  mutate(horizon = as.numeric(horizon)) %>% 
  ggplot(aes(x = horizon, y = rmse, col = siteID, linetype = domain)) +
  geom_line(size=1) +
  ylim(0,2.75) +
  labs(x = "Horizon (days)", y = y_label, tag = "(a)", color = "NEON Site ID", linetype = "NEON Domain") +
  scale_linetype_manual(values=c("solid", "dashed", "dotted")) +
  scale_color_manual(values = cbPalette) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 8),
        legend.key.width = unit(1, "cm"),
        legend.key.height = unit(0.3, "cm"),
        legend.position = "bottom",
        legend.margin=margin(0,0,0,0),
        plot.tag.position = c(0.25,0.88))

Fig2b <- combined_1m %>% 
  mutate(horizon = as.numeric(horizon)) %>% 
  mutate(domain = factor(domain, levels = c("Northern Plains","Great Lakes","Southeast"))) %>% 
  pivot_wider(names_from = Model, values_from = rmse) %>% 
  mutate(rmse_skill =  (`Day-of-year mean` - `FLARE-GLM`)) %>% 
  ggplot(aes(horizon, rmse_skill, color = siteID, linetype = domain)) +
  geom_line(size=1) +
  geom_hline(aes(yintercept = 0)) +
  ylim(-1.5,3.4) +
  labs(x = "Horizon (days)", y = "RMSE Skill (Null - FLARE)", tag = "(b)", title = " ", color = "NEON Site ID", linetype = "NEON Domain") +
  annotate("text", x=20, y=3.15, label = "FLARE performs better", size = 3) +
  annotate("text", x=20, y=-1.45, label = "Null performs better", size = 3) +
  scale_color_manual(values = cbPalette) +
  theme_bw() +
  scale_linetype_manual(values=c("solid", "dashed", "dotted")) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        legend.key.height = unit(0.3, "cm"),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 8),
        legend.margin=margin(0,0,0,0),
        plot.tag.position = c(0.25,0.88))

Fig2 <- Fig2a + Fig2b + plot_layout(guides = 'collect') & theme(legend.position = "bottom") &   guides(color=guide_legend(nrow=2, title.position = "top"),linetype = guide_legend(nrow=3,byrow=TRUE, title.position = "top"))

Fig2

ggsave(filename = file.path(lake_directory, "notebook","Figure2.jpg"), device = "jpeg", plot = Fig2, width = 4.5, height = 3.5, units = "in", dpi = 350)
```

## Figure 3

```{r} 
Fig3 <- spear_corr_table %>% 
  mutate(rho = ifelse(variable == "mean_secchi_depth_m" & horizon == "1", rho + 0.006, rho)) %>% 
  filter(horizon %in% c("1","max_min")) %>% 
  ggplot(aes(x = rho, y = full_name, col = Metric)) +
  geom_point() +
  geom_vline(xintercept = 0.5) +
  geom_vline(xintercept = -0.5) +
  geom_hline(yintercept = 5.5, linetype = "dashed") +
  xlim(-1,1) +
  labs(x = "Spearman's correlation coefficient", y = "variable") +
  scale_color_manual(values=c("red", "blue")) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position= "bottom",
        axis.text.x = element_text(color = "black"),
        axis.text.y = element_text(color = "black"))
Fig3

ggsave(filename = file.path(lake_directory, "notebook","Figure3.jpg"), device = "jpeg", plot = Fig3, width = 4.5, height = 4.5, units = "in", dpi = 350)

```

# Values in text

Number of days that had forecasts produced

```{r}
length(unique(combined1$forecast_start_time))
```

What is the max RMSE for all horizon and lakes

```{r}
rmse_all_mean <- combined_1m %>% 
  filter(Model == "FLARE-GLM") %>% 
  arrange(desc(rmse))

round(max(rmse_all_mean$rmse),4)
```

Horizon where climatology starts to perform better than FLARE.  If a site is not included in the table below that means that FLARE is better across all horizon.

```{r}
clim_vs_flare_horizon <- combined_1m %>% 
  select(-domain) %>% 
  pivot_wider(names_from = Model, values_from = rmse) %>%
  mutate(diff = `Day-of-year mean` - `FLARE-GLM`) %>%
  filter(diff < 0) %>%
  group_by(siteID) %>%
  summarize(min = min(as.numeric(horizon), na.rm = TRUE))
clim_vs_flare_horizon
```

What are the RMSEs at different horizons for the different lakes? The 1 and 35 day means are in the abstract

```{r}
rmse_site_mean <- combined_1m %>% 
  filter(Model == "FLARE-GLM") %>%
  select(-domain) %>% 
  filter(horizon %in% c(1, 7, 21, 35)) %>% 
  mutate(horizon = as.numeric(horizon)) %>% 
  group_by(horizon) %>% 
  summarize(mean = round(mean(rmse), 2),
            min = round(min(rmse), 2),
            max = round(max(rmse), 2))
rmse_site_mean %>% 
  arrange(horizon)

rmse_site_mean
```

How much more do forests in the northern sites degrade than the southern sites?

```{r}
d <- day1_max_min %>% 
  filter(horizon == "max_min") %>% 
  mutate(domain = ifelse(domain == "Southeast", "south", "north")) %>% 
  group_by(domain) %>% 
  summarize(mean = mean(rmse)) %>% 
  pivot_wider(names_from = domain, values_from = mean)

round(d$north/d$south, 2)

```


```{r}
TableS1 <- spear_corr_table %>% 
  filter(abs(rho) > 0.5) %>%   
  mutate(analysis = "spearmen")%>% 
  select(variable, horizon, rho) %>% 
  rename(metric = horizon) %>% 
  mutate(metric = ifelse(metric == 1, "accuracy","degradation")) %>% 
  arrange(metric, variable) %>% 
  left_join(variables) %>% 
  select(full_name_1line, metric, rho) %>% 
  rename(variable = full_name_1line) %>% 
  mutate(rho = round(rho, 2))

knitr::kable(TableS1)

```

#Supplemental Information Figures

## Figure S1 

```{r}

y_label <- expression(paste('Water temprature (',degree,'C)', sep = ""))
cbPalette <- c("#CC79A7", "#009E73", "#E69F00")

FigS1 <- combined %>%
  rename("siteID" = theme,
         Model = team) %>%
  select(siteID, Model, time, forecast_start_time, depth, mean, observed, quantile10, quantile90) %>%
  mutate(domain = NA,
         domain = ifelse(siteID %in% c("BARC","SUGG"), "Southeast", domain),
         domain = ifelse(siteID %in% c("CRAM","LIRO"), "Great Lakes", domain),
         domain = ifelse(siteID %in% c("PRLA","PRPO"), "Northern Plains", domain),
         domain = factor(domain, levels = c("Northern Plains","Great Lakes","Southeast")),
         siteID = factor(siteID, levels = c("PRLA","PRPO","CRAM","LIRO","BARC","SUGG"))) %>% 
  filter(as_date(forecast_start_time) %in% as_date(c("2021-06-15", "2021-07-01")),
         depth == 0.0,
         Model == "ms2_glm_flare") %>%
  mutate(forecast_start_time = factor(forecast_start_time)) %>% 
  ggplot(aes(x = time)) +
  geom_ribbon(aes(ymin = quantile10, ymax = quantile90, color = forecast_start_time, fill = forecast_start_time), alpha = 0.3) +
  geom_line(aes(y = mean, color = forecast_start_time)) +
  geom_point(aes(y = observed), color = "black", size = 0.5) +
  facet_wrap(~siteID,nrow = 3) +
  labs(x = "Date", y = y_label, color = "Forecast start date", fill = "Forecast start date") +
  scale_color_manual(values = cbPalette) +
  scale_fill_manual(values = cbPalette) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 8),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        strip.text = element_text(size = 10),
        strip.background = element_blank(),
        strip.text.x = element_text(margin = margin(0.0,0,0.1,0, "cm")),
        legend.position = "bottom")

ggsave(filename = file.path(lake_directory, "notebook","FigureS1.jpg"), device = "jpeg", plot = FigS1, width = 4.5, height = 6.5, units = "in", dpi = 350)
```


## Figure S2

Averaging all scores for all depths within each lake (Panel a)

```{r}
y_label <- expression(paste('RMSE (',degree,'C)', sep = ""))
cbPalette <- c("#a50026", "#f46d43", "#fdae61", "#fee090","#74add1", "#313695")

FigS2 <- combined_all_depths %>% 
  filter(Model == "FLARE-GLM") %>% 
  mutate(horizon = as.numeric(horizon)) %>% 
  ggplot(aes(x = horizon, y = rmse, col = siteID, linetype = domain)) +
  geom_line(size = 1) +
  ylim(0,3) +
  labs(x = "Horizon (days)", y = y_label, color = "NEON Site ID", linetype = "NEON Domain") +
  scale_color_manual(values = cbPalette) +
  scale_linetype_manual(values=c("solid", "dashed", "dotted")) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        legend.key.height = unit(0.3, "cm"),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 8),
        legend.margin=margin(0,0,0,0),
        plot.tag.position = c(0.30,0.95))

FigS2 <- FigS2 + plot_layout(guides = 'collect') & theme(legend.position = "bottom") &   guides(color=guide_legend(nrow=2, title.position = "top"),linetype = guide_legend(nrow=3,byrow=TRUE, title.position = "top"))


ggsave(filename = file.path(lake_directory, "notebook","FigureS2.jpg"), device = "jpeg", plot = FigS2, width = 4.5, height = 4.1, units = "in", dpi = 350)
```

## Figure S3

Relationship between the 1-day RMSE and the lake characteristic with colors representing the sites and shapes representing the 3 domains

```{r}
y_label <- expression(paste('Accurancy (-RMSE) (',degree,'C)', sep = ""))
cbPalette <- c("#a50026", "#f46d43", "#fdae61", "#fee090","#74add1", "#313695")
FigS3 <- day1_max_min %>% 
  select(horizon, siteID, rmse, max_depth_m, fetch_m, volume_m3, surface_area_km2, mean_secchi_depth_m,
         mean_annual_temperature_degreeC, mean_annual_precipitation_mm, air_temp_10day_sd_degreeC,
         mean_hydraulic_residence_time_yrs, catchment_size_km2, latitude) %>% 
  pivot_longer(cols = -c("rmse","siteID", "horizon"), names_to = "variable", values_to = "value") %>% 
  left_join(variables, by = "variable") %>% 
  filter(horizon == "1") %>% 
    mutate(domain = NA,
         domain = ifelse(siteID %in% c("BARC","SUGG"), "Southeast", domain),
         domain = ifelse(siteID %in% c("CRAM","LIRO"), "Great Lakes", domain),
         domain = ifelse(siteID %in% c("PRLA","PRPO"), "Northern Plains", domain),
         domain = factor(domain, levels = c("Northern Plains","Great Lakes","Southeast")),
         siteID = factor(siteID, levels = c("PRLA","PRPO","CRAM","LIRO","BARC","SUGG"))) %>% 
  ggplot(aes(x = value, y = -rmse, col = siteID, shape = domain)) +
  geom_point() +
  scale_color_manual(values = cbPalette) +
  facet_wrap(~full_name_1line, scale = "free", ncol = 2) +
  labs(y = "accuracy (-RMSE)", color = "NEON Site ID", shape = "NEON Domain") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 7),
        legend.spacing.x = unit(0.01,"cm"),
        legend.key.width = unit(1, "cm"),
        legend.key.height = unit(0.2, "cm"),
        legend.text = element_text(size = 7),
        legend.title = element_text(size = 7),
        strip.text = element_text(size = 7),
        strip.background = element_blank(),
        strip.text.x = element_text(margin = margin(0.0,0,0.1,0, "cm")),
        legend.position = "bottom",
        legend.margin=margin(0,0,0,0),
        plot.tag.position = c(0.3,0.93)) +
  guides(color=guide_legend(nrow=2,byrow=TRUE, 
                            title.position = "top"),
         shape = guide_legend(nrow=2,byrow=TRUE, 
                                 title.position = "top")) 
ggsave(filename = file.path(lake_directory, "notebook","FigureF3.jpg"), device = "jpeg", plot = FigS3, width = 4.5, height = 6.5, units = "in", dpi = 350)
```

## Figure S4

Relationship between the 1-day RMSE and the lake characteristic with colors representing the sites and shapes representing the 3 domains

```{r}
y_label <- expression(paste('Degradation (Max - Min RMSE: ',degree,'C)', sep = ""))
cbPalette <- c("#a50026", "#f46d43", "#fdae61", "#fee090","#74add1", "#313695")
FigS4 <- day1_max_min %>% 
  select(horizon, siteID, rmse, max_depth_m, fetch_m, volume_m3, surface_area_km2, mean_secchi_depth_m,
         mean_annual_temperature_degreeC, mean_annual_precipitation_mm, air_temp_10day_sd_degreeC,
         mean_hydraulic_residence_time_yrs, catchment_size_km2, latitude) %>% 
  pivot_longer(cols = -c("rmse","siteID", "horizon"), names_to = "variable", values_to = "value") %>% 
  left_join(variables, by = "variable") %>% 
  filter(horizon == "max_min") %>% 
    mutate(domain = NA,
         domain = ifelse(siteID %in% c("BARC","SUGG"), "Southeast", domain),
         domain = ifelse(siteID %in% c("CRAM","LIRO"), "Great Lakes", domain),
         domain = ifelse(siteID %in% c("PRLA","PRPO"), "Northern Plains", domain),
         domain = factor(domain, levels = c("Northern Plains","Great Lakes","Southeast")),
         siteID = factor(siteID, levels = c("PRLA","PRPO","CRAM","LIRO","BARC","SUGG"))) %>% 
  ggplot(aes(x = value, y = -rmse, col = siteID, shape = domain)) +
  geom_point() +
  scale_color_manual(values = cbPalette) +
  facet_wrap(~full_name_1line, scale = "free", ncol = 2) +
  labs(y = y_label, color = "NEON Site ID", shape = "NEON Domain") +
  theme_bw() +
 theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 7),
        legend.spacing.x = unit(0.01,"cm"),
        legend.key.width = unit(1, "cm"),
        legend.key.height = unit(0.2, "cm"),
        legend.text = element_text(size = 7),
        legend.title = element_text(size = 7),
        strip.text = element_text(size = 7),
        strip.background = element_blank(),
        strip.text.x = element_text(margin = margin(0.0,0,0.1,0, "cm")),
        legend.position = "bottom",
        legend.margin=margin(0,0,0,0),
        plot.tag.position = c(0.3,0.93)) +
  guides(color=guide_legend(nrow=2,byrow=TRUE, 
                            title.position = "top"),
         shape = guide_legend(nrow=2,byrow=TRUE, 
                                 title.position = "top")) 
ggsave(filename = file.path(lake_directory, "notebook","FigureS4.jpg"), device = "jpeg", plot = FigS4, width = 4.5, height = 6.5, units = "in", dpi = 350)
```

```