usgs_metab_long-term_analysis_v4.Rmd

---
output:
  pdf_document: default
  html_document: default
---

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

# Data for: Is the short-term phenology and long-term trajectory of river productivity changing over time?

## Nick Marzolf, 2024-01-29

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

# Introduction

RMarkdown file to accompany: Is the short-term phenology and long-term trajectory of river productivity changing over time?, a manuscript submitted to Limnology and Oceanography: Letters Special Issue: "Changing Phenology in Aquatic Ecosystems".

Nicholas S. Marzolf^1,\*^, Michael J. Vlah^1^, Heili E. Lowman^1^, Weston M. Slaughter^1^, Emily S. Bernhardt^1^

Affiliations

^1^ Department of Biology, Duke University, Durham, NC, USA

\*corresponding author: [nick.marzolf\@jonesctr.org](mailto:nick.marzolf@jonesctr.org){.email} (ORCiD: 0000-0001-9146-1643)

To recreate this HTML, load the data files into a folder named 'data' at the same path as the .rmd file, open the .rmd file, and click 'Knit'.

# Load packages

```{r load-packages, echo = TRUE, message = FALSE, warning = FALSE}
# data manipulation
library(tidyverse)
library(dplyr)
library(lubridate)
library(purrr)
library(tibble)
library(broom)

# plotting
library(ggplot2)
library(ggpubr)
library(ggExtra)
library(GGally)
library(gghighlight)
source("code/theme_nick.R")
ggplot2::theme_set(theme_nick())

# spatial
library(sf)
library(spData)

# statistics
library(lme4)
library(nlme)
library(performance)
library(EnvStats)
library(EflowStats)
library(car)
library(trend)
library(merTools)

# reproducability
library(pander)
```

```{r define-constants, echo=TRUE, message=FALSE, warning=FALSE}
# define constants and axis labels
metab_units_area <- (expression(paste('GPP (g C ', m^-2, ' ',y^-1,')', sep = ' ')))

low_p <- '#103801'
high_p <- '#43fa00'

# quotient to convert GPP from O2 to C
# g C = g O2 * (1 mol O2/32 g O2) * (2 mol O/1 mol O2) * (1 mol C/1 mol O) * (12 g C/1 mol C)
rq <- (1/32)*(2/1)*(1/1)*(12/1) 

# where to save figures
save_figs <- TRUE

dpi <- 400

version <- 'v4'

fig_dir <- glue::glue('manuscript/long-term_GPP/figures/{version}/')

if(!dir.exists(fig_dir)){
  dir.create(fig_dir)
}
```



# Load GPP data

```{r load-data, warning = FALSE, echo = TRUE}

# modified function from 'code/metabolism_synthesis/Part II- Filtering and gap filling.Rmd
## this function source calls code/functions/fillMiss3.R, which is the 'workhorse' of this function
source('code/functions/fill_and_normalize_metab.R')

table_s1 <- readr::read_csv('data/data_citation/1_site_info.csv')
river_metab <- readr::read_csv('data/data_citation/5_usgs_metabolism.csv') 

# gap fill
river_metab_filled <- fill_and_normalize_metab(river_metab) 

head(river_metab_filled)

n_obs <- length(river_metab_filled$GPP_filled)

n_sites <- length(unique(river_metab_filled %>% 
                           dplyr::filter(!is.na(site)) %>% 
                           dplyr::pull(site)))

```

# 1) Calculate cumulative daily GPP and total annual GPP

```{r GPP-cumsum, warning = FALSE, echo = TRUE}

# calculate site-year cumulative GPP
river_metab_cumsum <- river_metab_filled %>% 
  dplyr::mutate(year = lubridate::year(date)) %>% 
  dplyr::group_by(site, year) %>% 
  dplyr::mutate(jday = lubridate::yday(date),
                GPP_C = GPP_filled*rq,
                GPP_C_clean = replace(GPP_C, GPP_C < 0, 0), # replace negative values with 0 in cumulative counts
                GPP_C_cumsum = cumsum(GPP_C_clean)) 

# calculate annual GPP for each site-year
metab_all_annual <- river_metab_cumsum %>% 
  #dplyr::mutate(year = lubridate::year(date)) %>% 
  dplyr::group_by(site, year) %>% 
  dplyr::summarise(GPP_daily_mean = mean(GPP_C_clean, na.rm = TRUE),
                   GPP_ann_sum = sum(GPP_C_clean, na.rm = TRUE),
                   GPP_daily_cv = EnvStats::cv(GPP_C_clean, 
                                               na.rm = TRUE, 
                                               method = 'l.moments')*100,
                   Q_daily_mean  = mean(discharge, na.rm = TRUE),
                   Q_daily_cv = EnvStats::cv(discharge,
                                             na.rm = TRUE,
                                             method = 'l.moments')*100,
                   temp_daily_mean = mean(temp.water, na.rm = TRUE),
                   temp_daily_cv = EnvStats::cv(temp.water,
                                                na.rm = TRUE,
                                                method = 'l.moments')*100) %>% 
  dplyr::mutate(clean_site_num = as.numeric(substr(site, 6,30))) %>% 
  dplyr::left_join(table_s1 %>% 
                     rename(clean_site_num = `Site Number`,
                            clean_name = `Clean Name`) , 
                   'clean_site_num')


# calculate median annual and CV for each site time series
site_mean_ann_gpp <- metab_all_annual %>% 
  dplyr::group_by(site, clean_name) %>% 
  dplyr::summarise(median_ann_GPP = median(GPP_ann_sum, na.rm = TRUE),
                   sd_ann_GPP = sd(GPP_ann_sum, na.rm = TRUE),
                   min_ann_GPP = min(GPP_ann_sum, na.rm = TRUE),
                   max_ann_GPP = max(GPP_ann_sum, na.rm = TRUE),
                   GPP_percent_diff = ((max_ann_GPP - min_ann_GPP)/min_ann_GPP)*100,
                   GPP_ann_cv = EnvStats::cv(GPP_ann_sum, na.rm = TRUE, method = 'l.moments')*100,
                   mean_ann_temp = mean(temp_daily_mean, na.rm = TRUE),
                   mean_sd_temp = sd(temp_daily_mean, na.rm = TRUE),
                   cv_mean_temp = EnvStats::cv(temp_daily_mean, method = 'l.moments')*100,
                   median_ann_discharge = median(Q_daily_mean, na.rm = TRUE),
                   median_sd_discharge = sd(Q_daily_mean, na.rm = TRUE),
                   cv_median_discharge = EnvStats::cv(Q_daily_mean, na.rm = TRUE, method = 'l.moments')*100) %>% 
  dplyr::arrange(desc(max_ann_GPP)) %>% 
  dplyr::mutate(clean_site_num = as.numeric(substr(site, 6,30))) %>% 
  dplyr::left_join(table_s1 %>% 
                     rename(clean_site_num = `Site Number`) , 
                   'clean_site_num')

```

```{r fig-1, fig.width=15, fig.height=15}
# cumulative site-year GPP
river_metab_cumsum_names <- river_metab_cumsum %>% 
  dplyr::mutate(clean_site_num = as.numeric(substr(site, 6,30))) %>% 
  dplyr::left_join(table_s1 %>% 
                     rename(clean_site_num = `Site Number`) , 
                   'clean_site_num') %>% 
  dplyr::rename(clean_name = `Clean Name`)


river_metab_cumsum_names$Name <- suppressWarnings(forcats::fct_relevel(river_metab_cumsum_names$Name,
                                                                       site_mean_ann_gpp %>% 
                                                                         dplyr::pull(clean_name)))

y_axis_title <- bquote(Cumulative ~ .(metab_units_area))


# subset of 9 sites to show as figure 1 (for size/space concerns)
highest_4 <- site_mean_ann_gpp %>% 
  arrange(desc(median_ann_GPP)) %>% 
  head(., n = 4) %>% 
  dplyr::pull(site)

lowest_4 <- site_mean_ann_gpp %>% 
  arrange(median_ann_GPP) %>% 
  head(., n = 4) %>% 
  dplyr::pull(site)

extremes <- c(lowest_4, highest_4)


fig_1 <- ggplot(data = river_metab_cumsum_names %>% 
                  dplyr::filter(site %in% extremes),
                aes(x = jday,
                    y = GPP_C_cumsum,
                    group = interaction(year, site),
                    color = GPP_C_cumsum))+
  geom_line(show.legend = FALSE)+
  facet_wrap(forcats::fct_relevel(clean_name,
                                  site_mean_ann_gpp %>% 
                                    dplyr::filter(site %in% extremes) %>% 
                                    dplyr::pull(clean_name)) ~ .,
             scales = 'free_y',
             ncol = 4)+
  scale_color_gradient(low = low_p,
                       high = high_p)+
  scale_y_continuous(breaks = function(y) {
    max_value <- max(y, na.rm = TRUE)*0.9
    tick_marks <- seq(0, round(max_value, digits = -1), length.out = 3)
    return(tick_marks)
  }
  )+
  labs(x = 'Day of Year',
       y = expression(paste('Cumulative GPP (g C  ', m^-2, '',y^-1,')')))+
  theme(strip.text = element_text(size = 7),
        strip.background = element_rect(fill = 'white'))
fig_1

if(save_figs){
  ggsave(plot = fig_1,
         glue::glue(fig_dir, 'fig_1.png'),
         dpi = 800, 
         width = 12, height = 4)
}

```

# 2) Annual productivity range

```{r viz-ann-prod, warning = FALSE, echo = TRUE, fig.width=11, fig.height=5}

cv <- site_mean_ann_gpp %>% 
  dplyr::select(site, GPP_ann_cv) %>% 
  left_join(.,metab_all_annual, 'site')

cv_sum <- cv %>% 
  group_by(site) %>% 
  summarise(mean_ann_GPP_C = mean(GPP_ann_sum, na.rm = TRUE), 
            sd_ann_GPP_C = sd(GPP_ann_sum, na.rm = TRUE),
            median_daily_CV = median(GPP_daily_cv, na.rm = TRUE),
            min_daily_CV = min(GPP_daily_cv), 
            max_daily_CV = max(GPP_daily_cv),
            GPP_ann_cv = first(GPP_ann_cv),
            GPP_range = max(GPP_ann_sum) - min(GPP_ann_sum),
            CV_ratio = median_daily_CV/GPP_ann_cv)

# boxplot of annual GPP for each site
fig_2_a <- ggplot(metab_all_annual,
                  aes(x = forcats::fct_relevel(clean_name,
                                               site_mean_ann_gpp %>% 
                                                 dplyr::arrange(desc(median_ann_GPP)) %>% 
                                                 dplyr::pull(clean_name)), 
                      y = GPP_ann_sum))+
  geom_boxplot(width = 0.75,
               outlier.shape = NA,
               show.legend = FALSE)+
  geom_point(aes(color = GPP_ann_sum),
             size = 1)+
  coord_flip()+
  scale_y_continuous(name = metab_units_area,
                     limits = c(-1, 4000))+
  scale_color_gradient(name = metab_units_area,
                       low = low_p,
                       high = high_p)+
  theme(axis.text.y = element_text(angle = 0, size = 8,
                                   hjust = 0.95,vjust = 0.2),
        axis.title.y = element_blank(),
        legend.position = c(0.9,0.9),
        legend.justification = c(0.9,0.9),
        legend.background = element_rect(fill = NA, 
                                         color = 'black'))
fig_2_a

# Are more variable sites more productive?
summary(lm(data = site_mean_ann_gpp,
           median_ann_GPP ~ GPP_ann_cv))
# No, p = 0.50

# what if we remove the highly variable outliers?
summary(lm(data = site_mean_ann_gpp %>% 
             dplyr::filter(sd_ann_GPP < 800),
           median_ann_GPP ~ GPP_ann_cv))
# Yes, p = 0.06


fig_2_b <- ggplot2::ggplot(data = site_mean_ann_gpp,
                           aes(y = median_ann_GPP, 
                               x = GPP_ann_cv,
                               group = site,
                               color = median_ann_GPP))+
  ggplot2::geom_point(size = 2)+
  geom_errorbar(aes(ymin = median_ann_GPP - sd_ann_GPP,
                    ymax = median_ann_GPP + sd_ann_GPP))+
  ggplot2::labs(x = expression(paste(CV[GPP-Inter],' (%)')),
                y = expression(paste('Median Annual GPP (g C  ',m^-2, ' ',y^-1,')')))+
  lims(y = c(0, 3500),
       x = c(0, 40))+
  scale_color_gradient(low = low_p,
                       high = high_p)+
  theme(legend.position = 'none')
fig_2_b


fig_2_c <- ggplot(cv,
                  aes(x = GPP_ann_cv,
                      y = GPP_daily_cv,
                      color = GPP_ann_sum,
                      group = site))+
  geom_point()+
  geom_line()+
  geom_smooth(method = 'lm')+
  geom_abline(intercept = 0, slope = 1, linetype = 'dashed')+
  scale_color_gradient(name = metab_units_area,
                       low = low_p,
                       high = high_p,)+
  lims(x = c(5,40),
       y = c(0, 100))+
  labs(x = expression(paste(CV[GPP-Inter], ' (%)')),
       y = expression(paste(CV[GPP-Intra], ' (%)')))+
  theme(legend.position = 'none')



fig_2 <- ggpubr::ggarrange(fig_2_a,NA,
                           ggpubr::ggarrange(fig_2_b,
                                             fig_2_c, 
                                             nrow = 2,
                                             labels = c('b', 'c'),
                                             align = 'v'),
                           widths = c(1.5,0,1),
                           labels = 'a',
                           # hjust = c(-8.5,-8),
                           # vjust = c(1.7, 1.7),
                           ncol = 3,
                           align = 'v')
fig_2

if(save_figs){
  ggsave(plot = fig_2,
         glue::glue(fig_dir, 'fig_2.png'),
         dpi = 650,
         width = 11, height = 8)
}
```

# 3) When are rivers productive

```{r gpp-quantiles, warning = FALSE, echo = TRUE, fig.width=6, fig.height=10}
df_out <- data.frame(
  site = rep(rep(unique(river_metab_cumsum_names$site), each = 14),times = 4),
  year = rep(rep(unique(river_metab_cumsum_names$year), times = 59),each = 4),
  quantile = rep(c(0.25, 0.5, 0.75, 0.95), times = 14*59)
)

for(i in 1:nrow(df_out)){
  
  deets <- df_out[i,] 
  site <- deets[,'site']
  year <- deets[,'year']
  quantile <- deets[,'quantile']
  
  df_use <- river_metab_cumsum_names %>% 
    dplyr::filter(site %in% !!site,
                  year %in% !!year)
  
  if(nrow(df_use) == 0)
    next
  
  cumsum <- max(df_use$GPP_C_cumsum, na.rm = TRUE)
  df_out[i, 'GPP_C_cumsum'] <- cumsum
  
  df_math <- df_use %>% 
    dplyr::mutate(GPP_C_cumsum_cdf = GPP_C_cumsum/cumsum) 
  
  jday <- which.min(abs(df_math$GPP_C_cumsum_cdf - quantile))
  df_out[i,'jday'] <- jday
  
  
  row <- data.frame(df_math %>% 
                      dplyr::filter(jday == !!jday))
  df_out[i, 'clean_name'] <- row[,'clean_name']
  df_out[i, 'date'] <- row[,'date']
  df_out[i, 'GPP_C_cumsum'] <- row[,'GPP_C_cumsum']
  df_out[i, 'cdf'] <- row[,'GPP_C_cumsum_cdf']
  
} # end for loop

df_quants <- df_out %>% 
  dplyr::filter(!is.na(jday)) %>% 
  dplyr::group_by(site) %>% 
  tidyr::fill(clean_name, .direction = 'down') 

df_quants %>% 
  dplyr::group_by(quantile) %>% 
  dplyr::summarise(jday_median = median(jday, na.rm = TRUE),
                   jday_sd = sd(jday, na.rm = TRUE),
                   jday_min = min(jday, na.rm = TRUE),
                   jday_max = max(jday, na.rm = TRUE))

df_quants %>% 
  dplyr::group_by(quantile, clean_name) %>% 
  dplyr::summarise(jday_median = median(jday, na.rm = TRUE),
                   jday_sd = sd(jday, na.rm = TRUE),
                   jday_min = min(jday, na.rm = TRUE),
                   jday_max = max(jday, na.rm = TRUE)) 


# by quantile plots 

quant_25_rank <- df_quants %>%
  dplyr::filter(!is.na(jday),
                quantile == 0.25) %>% 
  dplyr::group_by(site, clean_name) %>% 
  dplyr::summarise(early = min(jday)) %>% 
  dplyr::arrange(early) %>% 
  dplyr::pull(clean_name)

quant_50_rank <- df_quants %>%
  dplyr::filter(!is.na(jday),
                quantile == 0.5) %>% 
  dplyr::group_by(site) %>% 
  dplyr::summarise(iqr = IQR(jday)) %>% 
  dplyr::arrange(desc(iqr)) %>% 
  dplyr::pull(site)

quant_75_rank <- df_quants %>%
  dplyr::filter(!is.na(jday),
                quantile == 0.75) %>% 
  dplyr::group_by(site) %>% 
  dplyr::summarise(iqr = IQR(jday)) %>% 
  dplyr::arrange(desc(iqr)) %>% 
  dplyr::pull(site)


quant_95_rank <- df_quants %>%
  dplyr::filter(!is.na(jday),
                quantile == 0.95) %>% 
  dplyr::group_by(site) %>% 
  dplyr::summarise(iqr = IQR(jday)) %>% 
  dplyr::arrange(desc(iqr)) %>% 
  dplyr::pull(site)



# temp changers
# temp_changers <- river_metab_cumsum_names %>% 
#   dplyr::filter(!is.na(temp.water)) %>% 
#   dplyr::group_by(clean_name) %>% 
#   tidyr::nest() %>% 
#   dplyr::mutate(sens = purrr::map(data, ~trend::sens.slope(x = .$temp.water)),
#                 mk_test = purrr::map(data, ~trend::mk.test(x = .$temp.water))) %>% 
#   dplyr::mutate(sens_p = purrr::map(sens, ~.$p.value), 
#                 sens_s = purrr::map_dbl(sens, ~.$estimates),
#                 mk_p = purrr::map_dbl(mk_test, ~.$p.value),
#                 mk_t = purrr::map_dbl(mk_test, ~.$estimates['tau']),
#                 sens_sig = ifelse(sens_p <= 0.05, 'significant', 'non-significant'),
#                 sens_slope = ifelse(sens_s > 0, 'warming', 'cooling')
#                 # mk_sig = ifelse(mk_p <= 0.05, 'significant', 'non-significant'),
#                 # mk_slope = ifelse(mk_t > 0, 'warming', 'cooling')
#                 ) 
# 
# temp_changers %>% 
#   filter(mk_p < 0.05) %>% 
#   group_by(sens_slope) %>% 
#   summarise(n = n())
# 
# warmers <- temp_changers %>% 
#   filter(mk_p < 0.05,
#          sens_slope == 'warming')
# 
# coolers <- temp_changers %>% 
#   filter(mk_p < 0.05,
#          sens_slope == 'cooling')

```

```{r quant-trends, echo=TRUE, warning = FALSE, fig.height=7, fig.width=7}
# trend analysis for quantiles
quantile_trends <- df_quants %>% 
  dplyr::arrange(site, year) %>% 
  dplyr::group_by(site, quantile) %>% 
  tidyr::nest() %>% 
  dplyr::mutate(sens = purrr::map(data, ~trend::sens.slope(x = .$jday)),
                mk_test = purrr::map(data, ~trend::mk.test(x = .$jday))) %>% 
  dplyr::mutate(sens_s = purrr::map_dbl(sens, ~.$estimates),
                sens_p = purrr::map(sens, ~.$p.value), 
                mk_p = purrr::map_dbl(mk_test, ~.$p.value),
                # mk_t = purrr::map_dbl(mk_test, ~.$estimates['tau']),
                sens_sig = ifelse(sens_p <= 0.05, 'significant', 'non-significant'),
                sens_slope = ifelse(sens_s > 0, 'increasing', 'decreasing'),
                mk_sig = ifelse(mk_p <= 0.05, 'significant', 'non-significant')
                #mk_slope = ifelse(mk_t > 0, 'increasing', 'decreasing')
  ) 



quantile_sum_sigs <- quantile_trends %>% 
  dplyr::ungroup() %>% 
  dplyr::select(quantile, site, 
                mk_p, #mk_t, 
                mk_sig, #mk_slope,
                sens_s, sens_sig, sens_slope)

quants_n_sig <- quantile_sum_sigs %>% 
  dplyr::filter(mk_p <= 0.05) %>% 
  dplyr::group_by(quantile, sens_slope) %>% 
  summarise(n = n()) %>% 
  data.frame()

quantile_trends_sum <- df_quants %>% 
  dplyr::left_join(quantile_sum_sigs) %>% 
  dplyr::mutate(sens_meaning = dplyr::case_when(sens_slope == 'increasing' ~ 'Later',
                                                sens_slope == 'decreasing' ~ 'Earlier'),
                # mk_meaning = case_when(mk_slope == 'increasing' ~ 'Later',
                #                        mk_slope == 'decreasing' ~ 'Earlier'),
  )

quantile_trenders <- quantile_trends_sum %>% 
  dplyr::filter(mk_p <= 0.05) %>% 
  dplyr::group_by(quantile, clean_name) %>%
  dplyr::filter(row_number() == 1)


fig_3 <- ggplot(quantile_trends_sum,
                aes(x = year,
                    y = jday,
                    group = site,
                    color = sens_meaning))+
  geom_line(linewidth = 0.8)+
  gghighlight::gghighlight(mk_sig == 'significant',
                           use_direct_label = FALSE,
                           keep_scales = TRUE,
                           unhighlighted_params = list(alpha("grey10", 0.4),
                                                       linewidth = 0.4),
                           calculate_per_facet = TRUE)+
  facet_wrap(. ~ quantile,
             scales = 'free_y', 
             nrow = 2, ncol = 2)+
  #scale_y_continuous(limits = c(0, 366))+
  scale_color_manual(name = element_blank(),
                     values = c(high_p, low_p))+
  labs(x = 'Year',
       y = 'Day of Year')+
  theme(legend.title = element_blank())
fig_3

if(save_figs){
  ggsave(plot = fig_3,
         glue::glue(fig_dir, 'fig_3.png'),
         dpi = 800,
         width = 7, height = 5)
}

```

# 4) Annual productivity trends

```{r prod-trends, warning = FALSE, echo = TRUE, fig.width=5, fig.height=5}

annual_C_sum_trends <- metab_all_annual %>% 
  # dplyr::select(site, year, GPP_ann_C) %>% 
  dplyr::group_by(site, clean_name) %>% 
  tidyr::nest() %>% 
  dplyr::mutate(sens = purrr::map(data, ~trend::sens.slope(x = .$GPP_ann_sum)),
                mk_test = purrr::map(data, ~trend::mk.test(x = .$GPP_ann_sum))) %>% 
  dplyr::mutate(sens_p = purrr::map(sens, ~.$p.value), 
                sens_s = purrr::map_dbl(sens, ~.$estimates),
                mk_p = purrr::map_dbl(mk_test, ~.$p.value),
                mk_t = purrr::map_dbl(mk_test, ~.$estimates['tau']),
                mk_sig = ifelse(mk_p <= 0.05, 'significant', 'non-significant'),
                mk_slope = ifelse(mk_t > 0, 'increasing', 'decreasing'),
                sens_sig = ifelse(sens_p <= 0.05, 'significant', 'non-significant'),
                sens_slope = ifelse(sens_s > 0, 'increasing', 'decreasing')) 


annual_C_sum_sigs <- annual_C_sum_trends %>% 
  dplyr::ungroup() %>% 
  dplyr::select(site,clean_name, 
                mk_p, mk_t, mk_sig, mk_slope,
                sens_s, sens_sig, sens_slope)


n_sig <- annual_C_sum_sigs %>% 
  dplyr::filter(mk_p <= 0.05) %>% 
  dplyr::group_by(mk_slope) %>% 
  summarise(n = n()) %>% 
  data.frame()

lt_changers <- annual_C_sum_sigs %>% 
  dplyr::filter(mk_p <= 0.05) %>% 
  dplyr::pull(site)

river_trends <- metab_all_annual %>% 
  dplyr::left_join(annual_C_sum_sigs) 


river_trends %>% 
  dplyr::filter(mk_p <= 0.05) %>% 
  dplyr::group_by(clean_name) %>%
  dplyr::filter(row_number() == 1)

fig_4_a <- ggplot2::ggplot(river_trends,
                           ggplot2::aes(x = year, 
                                        y = GPP_ann_sum,
                                        color = mk_slope))+
  ggplot2::geom_line(linewidth = 0.8,
                     aes(group = clean_name))+
  #geom_label_repel()+
  scale_y_continuous(limits = c(10, 4000))+
  scale_x_continuous(n.breaks = 10)+
  scale_color_manual(name = element_blank(),
                     values = c(low_p, high_p),
                     labels = c(glue::glue('Decreasing (n = ', n_sig[1,2],')'),
                                glue::glue('Increasing (n = ', n_sig[2,2],')')))+
  gghighlight::gghighlight(mk_sig == 'significant',
                           use_direct_label = FALSE,
                           keep_scales = TRUE,
                           unhighlighted_params = list(alpha("grey40", 0.4),
                                                       linewidth = 0.4))+
  ggplot2::labs(y = metab_units_area)+
  theme(axis.title.x = element_blank(),
        legend.title = element_blank(),
        legend.key.size = unit(0.3,'cm'),
        legend.justification = c(0.01, 0.99),
        legend.position = c(0.01, 0.99),
        legend.text = element_text(size = 8),
        legend.background = element_rect(color = 'black',
                                         linewidth = 0.25)
  )


fig_4_b <- river_trends %>% 
  filter(site %in% lt_changers) %>% 
  ggplot(aes(x = year,
             y = GPP_daily_cv,
             group = clean_name,
             color = mk_slope))+
  geom_line(linewidth = 0.8,
            show.legend = FALSE)+
  # geom_label(data = river_trends[!duplicated(label_indices), ],
  #            aes(x = year,
  #                y = GPP_daily_cv,
  #                label = clean_name),
  #            show.legend = FALSE)+
  # scale_color_gradient(name = element_blank(),
  #                      low = low_p, high = high_p)+
  scale_color_manual(name = element_blank(),
                     values = c(low_p, high_p),
                     # labels = c(glue::glue('Decreasing (n = ', n_sig[1,2],')'),
                     #            glue::glue('Increasing (n = ', n_sig[2,2],')'))
  )+
  labs(y = expression(paste(CV[GPP-Intra])))+
  lims(y = c(0, 100))+
  theme(axis.title.x = element_blank())


temp_stats <- river_trends %>% 
  group_by(site) %>% 
  summarise(mk_slope = dplyr::first(mk_slope),
            temp_out = list(tidy(lm(log10(GPP_ann_sum) ~ temp_daily_mean)))) %>% 
  unnest(temp_out) %>% 
  filter(term == 'temp_daily_mean')


metab_all_temp_stats <- dplyr::left_join(metab_all_annual, temp_stats, 'site') %>% 
  mutate(is_sig = ifelse(p.value <= 0.05,
                         'significant', 
                         'non-significant'), 
         direction = ifelse(estimate >= 0, 
                            'positive', 
                            'negative'))

metab_all_temp_stats %>% 
  dplyr::filter(site %in% lt_changers,
                is_sig == 'significant')

fig_4_c <- metab_all_temp_stats %>% 
  dplyr::filter(site %in% lt_changers) %>% 
  ggplot(.,
         aes(x = temp_daily_mean,
             y = GPP_ann_sum,
             color = mk_slope))+
  geom_point()+
  geom_smooth(aes(group = site),
              method = 'lm', 
              se = FALSE, 
              color = 'grey40')+
  # directlabels::geom_dl(label = 'text', method = "maxvar.points", inherit.aes = TRUE)+
  # geom_abline(slope = temp_all$coefficients[2],
  #             intercept = temp_all$coefficients[1],
  #             linetype = 'dashed',
  #             color = 'black')+
  gghighlight::gghighlight(is_sig == 'significant',
                           use_direct_label = FALSE,
                           unhighlighted_params = list(alpha("grey10", 0.4),
                                                       linewidth = 0.2))+
  scale_y_log10()+
  scale_color_manual(name = 'Long Term GPP Trend',
                     values = c(low_p, high_p),
                     labels = c('Decreasing','Increasing'))+
  labs(x = 'Mean Daily Temperature (°C)',
       y = metab_units_area)+
  theme(legend.position = 'none')



cv_Q_stats <- river_trends %>% 
  group_by(site) %>% 
  summarise(mk_slope = first(mk_slope),
            out = list(tidy(lm(log10(GPP_ann_sum) ~ Q_daily_cv)))) %>% 
  unnest(out) %>% 
  filter(term == 'Q_daily_cv')


metab_all_cvQ_stats <- left_join(metab_all_annual, cv_Q_stats, 'site') %>% 
  mutate(is_sig = ifelse(p.value <= 0.05,
                         'significant', 
                         'non-significant'), 
         direction = ifelse(estimate >= 0, 
                            'positive', 
                            'negative'))

metab_all_cvQ_stats %>% 
  filter(site %in% lt_changers,
         is_sig == 'significant') 

fig_4_d <- metab_all_cvQ_stats %>% 
  filter(site %in% lt_changers) %>% 
  ggplot(.,
         aes(x = Q_daily_cv,
             y = GPP_ann_sum,
             color = mk_slope))+
  geom_point()+
  geom_smooth(aes(group = site),
              method = 'lm', 
              se = FALSE, 
              color = 'grey40')+
  # geom_abline(slope = temp_all$coefficients[2],
  #             intercept = temp_all$coefficients[1],
  #             linetype = 'dashed',
  #             color = 'black')+
  gghighlight::gghighlight(is_sig == 'significant',
                           use_direct_label = FALSE,
                           unhighlighted_params = list(alpha("grey10", 0.4),
                                                       linewidth = 0.2))+
  # directlabels::geom_dl(label = as.factor(metab_all_cvQ_stats %>%
  #                                           filter(site %in% lt_changers) %>%
  #                                           pull(clean_name)),
  #                       method = "maxvar.points", inherit.aes = TRUE)+
  # facet_grid(. ~ is_sig)+
  scale_y_log10()+
  scale_color_manual(name = 'Long Term GPP Trend',
                     values = c(low_p, high_p),
                     labels = c('Decreasing','Increasing'))+
  labs(x = expression(paste(CV[Q-Intra], ' (%)')),
       y = metab_units_area)+
  theme(legend.position = c(0.01, 0.01),
        legend.justification = c(0.01, 0.01),
        legend.key.size = unit(0.3,'cm'),
        legend.text = element_text(size = 8),
        legend.background = element_rect(color = 'black',
                                         linewidth = 0.25),
        legend.title = element_text(size = 8))

fig_4 <- ggpubr::ggarrange(fig_4_a,
                           fig_4_b,
                           fig_4_c,
                           fig_4_d,
                           nrow = 2,
                           ncol = 2,
                           align = 'hv',
                           labels = 'auto')

if(save_figs){
  ggsave(plot = fig_4,
         glue::glue(fig_dir, 'fig_4.png'),
         dpi = 700,
         width = 10, height = 8)
}

```

# 5) Supplemental tables and figures
```{r supp-figs, warning = FALSE, echo = TRUE, fig.width=5, fig.height=5}

# Table S1
table_s1 <- readr::read_csv('data/data_citation/1_site_info.csv')
table_s1

# readr::write_csv(table_s1,
#                  'data/data_citation/ms_usgs_annual_metab/table_s1_site_info.csv')

# Table S2
sum_NA <- readr::read_csv('data/data_citation/ms_usgs_annual_metab/summary_site-year_NAs.csv') 

sum_NA_wide <- sum_NA %>% 
  dplyr::select(-start_date, -complete_year) %>% 
  tidyr::pivot_wider(names_from = year,
                     values_from = per_NA) %>% 
  dplyr::mutate(`Site Number` = as.numeric(substr(site, 6, 30))) %>% 
  left_join(.,table_s1, by = 'Site Number') %>% 
  dplyr::select(Name, '2008', '2009', '2010','2011','2012','2013','2014','2015','2016','2017','2018','2019','2020','2021')

readr::write_csv(sum_NA_wide,
                 'data/data_citation/ms_usgs_annual_metab/table_s2_sum_NA_wide.csv')

# Table S3
table_s3 <- quantile_trenders %>% 
  dplyr::select(quantile, clean_name, mk_p, sens_s, sens_meaning) %>% 
  dplyr::arrange(quantile)

readr::write_csv(table_s3,
                 'data/data_citation/ms_usgs_annual_metab/table_s3_quantiles.csv')


# table_s4 <- quantile_trends_sum %>% 
#   dplyr::filter(mk_p < 0.05) %>% 
#   dplyr::select(quantile, site, mk_t, mk_p,mk_meaning) %>% 
#   dplyr::mutate(across(where(is.numeric), round, 3)) %>% 
#   dplyr::group_by(quantile, site) %>% 
#   slice_head()
# 
# readr::write_csv(table_s4,
#                  'data/data_citation/ms_usgs_annual_metab/table_s4_quantile_trends.csv')

# Table S4
table_s4 <- annual_C_sum_sigs %>% 
  dplyr::select(clean_name, mk_p,sens_s) %>% 
  dplyr::mutate(across(where(is.numeric), round, 3))
colnames(table_s4) <- c('Site',
                        "Mann-Kendall p-value",
                        "Sen's Slope")

readr::write_csv(table_s4,
                 'data/data_citation/ms_usgs_annual_metab/table_s4_trend_stats.csv')

# Supp figures

# Figure S1
table_s1_plot <- table_s1 %>%
  sf::st_as_sf(.,
               coords = c('Longitude', 'Latitude'),
               crs = 4326) %>%
  rename(clean_name = `Clean Name`) %>%
  left_join(site_mean_ann_gpp, 'clean_name')

fig_s1_map <- ggplot()+
  geom_sf(data = world %>% 
            dplyr::filter(name_long == 'United States') %>% 
            sf::st_set_crs(4326),
          fill = 'white')+
  geom_sf(data = table_s1_plot,
          aes(color = median_ann_GPP,
              size = median_ann_GPP))+
  scale_color_gradient(name = metab_units_area,
                       low = low_p,
                       high = high_p)+
  scale_size_continuous(guide = 'none')+
  coord_sf(xlim = c(-130, -65), ylim = c(23, 55), expand = FALSE)
fig_s1_map


fig_s1_hist <- ggplot(table_s1,
                      aes(x = `Drainage Area (km2)`*2.58999))+ # convert to km2
  geom_histogram()+
  labs(x = expression(paste('Watershed area (', km^2,')')),
       y = 'Count')+
  scale_x_log10()
fig_s1_hist


fig_s1 <- ggpubr::ggarrange(fig_s1_map,NULL,
                            fig_s1_hist,
                            ncol = 3,
                            widths = c(3,-0.5,1.5),
                            align = 'v'
)
fig_s1

if(save_figs) {
  ggsave(plot = fig_s1,
         glue::glue(fig_dir,'fig_s1.png'),
         dpi = 1200,
         width = 13, height = 3)
}



# Figure S2: all sites as supplemental figure
fig_s2 <- ggplot(river_metab_cumsum_names, 
                 aes(x = jday,
                     y = GPP_C_cumsum,
                     group = interaction(year, site),
                     color = GPP_C_cumsum))+
  geom_line(show.legend = FALSE)+
  facet_wrap(forcats::fct_relevel(clean_name,
                                  site_mean_ann_gpp %>% 
                                    dplyr::pull(clean_name)) ~ .,
             scales = 'free_y',
             ncol = 5)+
  scale_color_gradient(low = low_p,
                       high = high_p)+
  scale_y_continuous(breaks = function(y) {
    max_value <- max(y, na.rm = TRUE)
    tick_marks <- seq(0, round(max_value, digits = -2), length.out = 3)
    return(tick_marks)
  }
  )+
  labs(x = 'Day of Year',
       y = expression(paste('Cumulative GPP (g C  ', m^-2, '',y^-1,')')))+
  theme(strip.text = element_text(size = 7),
        strip.background = element_rect(fill = 'white'))
fig_s2

if(save_figs){
  ggsave(plot = fig_s2,
         glue::glue(fig_dir, 'fig_s2.png'),
         dpi = 400, 
         width = 15, height = 13)
}




# Figure S3
temp_all <- lm(data = metab_all_annual %>% 
                 filter(temp_daily_mean > 5),
               log10(GPP_ann_sum) ~ temp_daily_mean)
summary(temp_all)

temp_by_site <- lm(data = metab_all_annual,
                   log10(GPP_ann_sum) ~ temp_daily_mean + site)
summary(temp_by_site)

temp_by_site_ranef <- lme4::lmer(data = metab_all_annual,
                                 log10(GPP_ann_sum) ~ temp_daily_mean + (1|site))
summary(temp_by_site_ranef)
car::Anova(temp_by_site_ranef)

fig_s3 <- metab_all_temp_stats %>% 
  filter(temp_daily_mean > 5) %>% 
  ggplot(.,
         aes(x = temp_daily_mean,
             y = GPP_ann_sum,
             color = mk_slope))+
  geom_point()+
  geom_smooth(aes(group = site),
              method = 'lm', 
              se = FALSE, 
              color = 'grey40')+
  geom_abline(slope = temp_all$coefficients[2],
              intercept = temp_all$coefficients[1],
              linetype = 'dashed',
              color = 'black')+
  gghighlight::gghighlight(is_sig == 'significant',
                           use_direct_label = FALSE,
                           unhighlighted_params = list(alpha("grey10", 0.4),
                                                       linewidth = 0.2))+
  scale_y_log10()+
  scale_color_manual(name = element_blank(),
                     values = c(low_p, high_p),
                     labels = c('Increasing','Decreasing'))+
  labs(x = 'Mean Daily Temperature (°C)',
       y = metab_units_area)+
  theme(legend.position = 'none')

if(save_figs){
  ggsave(plot = fig_s3,
         glue::glue(fig_dir,'fig_s3.png'),
         dpi = dpi, 
         width = 6, height = 6)
}



# Figure S4
cvQ_all <- lm(data = metab_all_annual, 
              log10(GPP_ann_sum) ~ Q_daily_cv)
summary(cvQ_all)

cvQ_by_site <- lm(data = metab_all_annual,
                  log10(GPP_ann_sum) ~ Q_daily_cv + site)
summary(cvQ_by_site)

cvQ_by_site_ranef <- lme4::lmer(data = metab_all_annual,
                                log10(GPP_ann_sum) ~ Q_daily_cv + (1|site))
summary(cvQ_by_site_ranef)
car::Anova(cvQ_by_site_ranef)

fig_s4 <- metab_all_temp_stats %>% 
  ggplot(.,
         aes(x = Q_daily_cv,
             y = GPP_ann_sum,
             color = mk_slope))+
  geom_point()+
  geom_smooth(aes(group = site),
              method = 'lm', 
              se = FALSE, 
              color = 'grey40')+
  geom_abline(slope = cvQ_all$coefficients[2],
              intercept = cvQ_all$coefficients[1],
              linetype = 'dashed',
              color = 'black')+
  gghighlight::gghighlight(is_sig == 'significant',
                           use_direct_label = FALSE,
                           unhighlighted_params = list(alpha("grey10", 0.4),
                                                       linewidth = 0.2))+
  scale_y_log10()+
  scale_color_manual(name = element_blank(),
                     values = c(low_p, high_p),
                     labels = c('Increasing','Decreasing'))+
  labs(x = expression(paste(CV[Q-Intra], ' (%)')),
       y = metab_units_area)+
  theme(legend.position = 'none')

if(save_figs){
  ggsave(plot = fig_s4,
         glue::glue(fig_dir,'fig_s4.png'),
         dpi = dpi, 
         width = 6, height = 6)
}
```



# Session Info

```{r session-info}
pander::pander(sessionInfo())
```