state_vars.Rmd

---
title: "Fourmile state variables"
author: "Megan Sears"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: yes
    toc_float: yes
    theme: journal
editor_options:
  chunk_output_type: console
---

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

library(tidyverse)
library(ggplot2); theme_set(theme_bw(base_size = 20))
library(plotly)
library(knitr)
library(kableExtra)
library(cowplot)

```

```{r, include = F}

statn <- read_csv('./38HRU_dec2024_updates/catchment_output/Catchment_statnall.csv') %>%
  mutate(date = mdy(time),
         source = 'statn') %>%
  dplyr::select(-time)

grid <- read_csv('./38HRU_dec2024_updates/catchment_output/Catchment_gridall.csv') %>%
  mutate(date = mdy(time),
         source = 'grid') %>%
  dplyr::select(-time)

both <- bind_rows(statn, grid)  
  
```

# Catchment

## Inputs

```{r}

temp <- ggplot(both, 
               aes(date, 
                   tmean,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(temp)

tmean_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_tmean = mean(tmean))
  
ggplot(tmean_avg, aes(x=month, y=mean_tmean,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean tmean (C)')


precip <- ggplot(both, 
               aes(date, 
                   precip,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(precip)

precip_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month, year) %>%
  summarize(sum_mo = sum(precip)) %>%
  ungroup() %>% 
  group_by(source, month) %>%
  summarize(mean_p = mean(sum_mo))
  
ggplot(precip_avg, aes(x=month, y=mean_p,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean precip (mm)')

solrad_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_solrad = mean(solRad))
  
ggplot(solrad_avg, aes(x=month, y=mean_solrad,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean solrad')


```

## Interception

- Monthly average: Statn > grid 1:4 and 10:12, grid > statn 5:9
- During the summer rain interception is greater from grid
- In the winter snow interception is greater from statn

```{r}
int_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_interception = mean(interception))
  
ggplot(int_avg, aes(x=month, y=mean_interception,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean interception')


interception <- ggplot(both, 
               aes(date, 
                   interception,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(interception)

```

## ET 

- Average ET: grid > statn
- Monthly average ET: grid > statn 1:5 and 9:12, statn > grid 6:8

```{r}

et_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_et = mean(actET...21))

ggplot(et_avg, aes(x=month, y=mean_et,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean ET (mm)')

et <- ggplot(both, 
               aes(date, 
                   actET...21,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='ET (mm)')

ggplotly(et)

```

## SWE 

- Despite grid being higher all winter, it melts out faster than statn
- grid generally has higher SWE, but statn starts accumulating earlier & melts slower


```{r}

swe_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_swe = mean(snowTotSWE))

ggplot(swe_avg, aes(x=month, y=mean_swe,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean SWE (l/day)')

snowTotSWE <- ggplot(both, 
                     aes(date, 
                         snowTotSWE,
                         color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='SWE (l/day)')

ggplotly(snowTotSWE)

```

## Soil sat (whole profile)

- 4,5,6,10 statn > SM than grid
- SM higher for statn

```{r}

soil_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_soil = mean(soilSat))

ggplot(soil_avg, aes(x=month, y=mean_soil,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean soil sat')

soil <- ggplot(both, 
               aes(date, 
                   soilSat,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(soil)

```

## Surface runoff

- RG1: very little other than Sept, where grid > statn


```{r}

rd1_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rd1 = mean(outRD1))

ggplot(rd1_avg, aes(x=month, y=mean_rd1,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean surface runoff (l)')

rd1 <- ggplot(both, 
               aes(date, 
                   outRD1,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='Surface runoff (l)')

ggplotly(rd1)

```

## Interflow

- Statn > grid
- very little grid

```{r}

rd2_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rd2 = mean(outRD2))

ggplot(rd2_avg, aes(x=month, y=mean_rd2,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean interflow (l)')

rd2 <- ggplot(both, 
               aes(date, 
                   outRD2,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='Interflow (l)')

ggplotly(rd2)

```

## Groundwater

- RG1 = shallow groundwater
- RG1: grid > statn

- RG2 = groundwater baseflow
- RG2: Very small relative to RG1
- RG2: Grid > statn

```{r}
gw1_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rg1 = mean(outRG1))

ggplot(gw1_avg, aes(x=month, y=mean_rg1,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean RG1 (unit?)')

gw1 <- ggplot(both, 
               aes(date, 
                   outRG1,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(gw1)

################################################################
gw2_avg <- both %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rg2 = mean(outRG2))

ggplot(gw2_avg, aes(x=month, y=mean_rg2,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean RG2 (unit?)')

gw2 <- ggplot(both, 
               aes(date, 
                   outRG2,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(gw2)

```

NOTHING BELOW HAS BEEN UPDATED WITH DEC2024 UPDATES

# High elevation HRU - old

HRU18

```{r}

statn_hi <- read_csv('./38HRU_dec2024_updates/HRU/statn_HRU_18.csv') %>% 
  mutate(date = mdy(time),
         source = 'statn') %>%
  dplyr::select(-time)

grid_hi <- read_csv('./38HRU_dec2024_updates/HRU/grid_HRU_18.csv') %>%
  mutate(date = mdy(time),
         source = 'grid') %>%
  dplyr::select(-time)

both_hi <- bind_rows(statn_hi, grid_hi)  

```

## Inputs

```{r}

tmean_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_tmean = mean(tmean))
  
ggplot(tmean_avg, aes(x=month, y=mean_tmean,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean tmean (C)')

temp <- ggplot(both_hi, 
               aes(date, 
                   tmean,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(temp)
############
precip_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month, year) %>%
  summarize(sum_mo = sum(precip)) %>%
  ungroup() %>% 
  group_by(source, month) %>%
  summarize(mean_p = mean(sum_mo))
  
ggplot(precip_avg, aes(x=month, y=mean_p,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean precip (mm)')

precip <- ggplot(both_hi, 
               aes(date, 
                   precip,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(precip)
#################

solrad_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_solrad = mean(solRad))
  
ggplot(solrad_avg, aes(x=month, y=mean_solrad,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean solrad (MJ/m2/day)')

```

## Interception

- Station > grid except for Jul

```{r}
int_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_interception = mean(interception))
  
ggplot(int_avg, aes(x=month, y=mean_interception,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean interception')


interception <- ggplot(both_hi, 
               aes(date, 
                   interception,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(interception)

```

## ET 

- Grid > statn except for Mar

```{r}

et_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_et = mean(actET))

ggplot(et_avg, aes(x=month, y=mean_et,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean ET (mm)')

et <- ggplot(both_hi, 
               aes(date, 
                   actET,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='ET (mm)')

ggplotly(et)

```

## SWE 

- Grid > statn

```{r}

swe_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_swe = mean(snowTotSWE))

ggplot(swe_avg, aes(x=month, y=mean_swe,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean SWE (l/day)')

snowTotSWE <- ggplot(both_hi, 
                     aes(date, 
                         snowTotSWE,
                         color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='SWE (l/day)')

ggplotly(snowTotSWE)

```

## Soil sat (whole profile)

- Statn > grid

```{r}

soil_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_soil = mean(soilSat))

ggplot(soil_avg, aes(x=month, y=mean_soil,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean soil sat')

soil <- ggplot(both_hi, 
               aes(date, 
                   soilSat,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(soil)

```

## Surface runoff

- No surface runoff

```{r}

rd1_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rd1 = mean(outRD1))

ggplot(rd1_avg, aes(x=month, y=mean_rd1,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean surface runoff (l)')

rd1 <- ggplot(both_hi, 
               aes(date, 
                   outRD1,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='Surface runoff (l)')

ggplotly(rd1)

```

## Interflow

- Statn > grid Jan to Apr
- Grid > statn May to Jun
- Other months are similar

```{r}

rd2_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rd2 = mean(outRD2))

ggplot(rd2_avg, aes(x=month, y=mean_rd2,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean interflow (l)')

rd2 <- ggplot(both_hi, 
               aes(date, 
                   outRD2,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='Interflow (l)')

ggplotly(rd2)

```

## Groundwater

- RG1 = shallow groundwater
- RG1: Jan to May statn > grid
- RG1: grid > statn Jun, Jul, Oct
- RG1: similar rest of year 

- RG2 = groundwater baseflow
- RG2: Very small relative to RG1
- RG2: Grid > statn

```{r}
gw1_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rg1 = mean(outRG1))

ggplot(gw1_avg, aes(x=month, y=mean_rg1,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean RG1 (unit?)')

gw1 <- ggplot(both_hi, 
               aes(date, 
                   outRG1,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(gw1)

################################################################
gw2_avg <- both_hi %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rg2 = mean(outRG2))

ggplot(gw2_avg, aes(x=month, y=mean_rg2,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean RG2 (unit?)')

gw2 <- ggplot(both_hi, 
               aes(date, 
                   outRG2,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(gw2)

```

# Low elevation HRU - old

HRU19

```{r}

statn_lo <- read_csv('./38HRU_sept2013_updates/HRU_19_statn.csv') %>% 
  mutate(date = mdy(time),
         source = 'statn') %>%
  select(-time)

grid_lo <- read_csv('./38HRU_sept2013_updates/HRU_19_grid.csv') %>%
  mutate(date = mdy(time),
         source = 'grid') %>%
  select(-time)

both_lo <- bind_rows(statn_lo, grid_lo)  

```

## Inputs

```{r}

tmean_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_tmean = mean(tmean))
  
ggplot(tmean_avg, aes(x=month, y=mean_tmean,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean tmean (C)')

temp <- ggplot(both_lo, 
               aes(date, 
                   tmean,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(temp)
############
precip_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month, year) %>%
  summarize(sum_mo = sum(precip)) %>%
  ungroup() %>% 
  group_by(source, month) %>%
  summarize(mean_p = mean(sum_mo))
  
ggplot(precip_avg, aes(x=month, y=mean_p,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean precip (mm)')

precip <- ggplot(both_lo, 
               aes(date, 
                   precip,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(precip)
#################

solrad_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_solrad = mean(solRad))
  
ggplot(solrad_avg, aes(x=month, y=mean_solrad,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean solrad (MJ/m2/day)')

```

## Interception

- Statn > grid Oct to Apr
- Grid > statn May to Sep

```{r}
int_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_interception = mean(interception))
  
ggplot(int_avg, aes(x=month, y=mean_interception,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean interception')


interception <- ggplot(both_lo, 
               aes(date, 
                   interception,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(interception)

```

## ET 

- Grid > statn Sep to May

```{r}

et_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_et = mean(actET))

ggplot(et_avg, aes(x=month, y=mean_et,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean ET (mm)')

et <- ggplot(both_lo, 
               aes(date, 
                   actET,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='ET (mm)')

ggplotly(et)

```

## SWE 

- Grid > statn Oct to Feb
- Statn > grid Mar to Apr and Sep

```{r}

swe_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_swe = mean(snowTotSWE))

ggplot(swe_avg, aes(x=month, y=mean_swe,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean SWE (l/day)')

snowTotSWE <- ggplot(both_lo, 
                     aes(date, 
                         snowTotSWE,
                         color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='SWE (l/day)')

ggplotly(snowTotSWE)

```

## Soil sat (whole profile)

- Statn > grid

```{r}

soil_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_soil = mean(soilSat))

ggplot(soil_avg, aes(x=month, y=mean_soil,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean soil sat')

soil <- ggplot(both_lo, 
               aes(date, 
                   soilSat,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(soil)

```

## Surface runoff

- Very little except for Sep 2013 (grid > statn)

```{r}

rd1_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rd1 = mean(outRD1))

ggplot(rd1_avg, aes(x=month, y=mean_rd1,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean surface runoff (l)')

rd1 <- ggplot(both_lo, 
               aes(date, 
                   outRD1,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='Surface runoff (l)')

ggplotly(rd1)

```

## Interflow

- Statn > grid Jan, May, June, July, Sep

```{r}

rd2_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rd2 = mean(outRD2))

ggplot(rd2_avg, aes(x=month, y=mean_rd2,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean interflow (l)')

rd2 <- ggplot(both_lo, 
               aes(date, 
                   outRD2,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y='Interflow (l)')

ggplotly(rd2)

```

## Groundwater

- RG1 = shallow groundwater
- RG1: Jan, Feb, May, Aug, Sep, Dec statn > grid

- RG2 = groundwater baseflow
- RG2: Very small relative to RG1
- RG2: Grid > statn

```{r}
gw1_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rg1 = mean(outRG1))

ggplot(gw1_avg, aes(x=month, y=mean_rg1,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean RG1 (unit?)')

gw1 <- ggplot(both_lo, 
               aes(date, 
                   outRG1,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(gw1)

################################################################
gw2_avg <- both_lo %>%
  mutate(month = month(date),
         year = year(date)) %>% 
  group_by(source, month) %>%
  summarize(mean_rg2 = mean(outRG2))

ggplot(gw2_avg, aes(x=month, y=mean_rg2,
                 color=source)) +
  geom_point(size=3) +
  scale_color_manual(values = c("#0047AB", "orange")) +
  labs(y= 'Monthly mean RG2 (unit?)')

gw2 <- ggplot(both_lo, 
               aes(date, 
                   outRG2,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(gw2)

```


# 2016 for low and high HRU

```{r}

statn_hi <- read_csv('./38HRU_dec2024_updates/HRU/statn_HRU_18.csv') %>% 
  mutate(date = mdy(time),
         source = 'statn') %>%
  dplyr::select(-time)

grid_hi <- read_csv('./38HRU_dec2024_updates/HRU/grid_HRU_18.csv') %>%
  mutate(date = mdy(time),
         source = 'grid') %>%
  dplyr::select(-time)

both_hi <- bind_rows(statn_hi, grid_hi) %>%
  mutate(year = year(date),
         month = month(date)) %>%
  filter(year == 2016) #
         #month %in% c(4,5,6))


temp <- ggplot(both_hi, 
               aes(date, 
                   tmean,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(temp)

swe <- ggplot(both_hi, 
               aes(date, 
                   snowTotSWE,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(swe)

```

```{r}

statn_low <- read_csv('./38HRU_dec2024_updates/HRU/statn_HRU_19.csv') %>% 
  mutate(date = mdy(time),
         source = 'statn') %>%
  dplyr::select(-time)

grid_low <- read_csv('./38HRU_dec2024_updates/HRU/grid_HRU_19.csv') %>%
  mutate(date = mdy(time),
         source = 'grid') %>%
  dplyr::select(-time)

both_low <- bind_rows(statn_low, grid_low) %>%
  mutate(year = year(date),
         month = month(date)) %>%
  filter(year == 2016,
         month %in% c(4,5,6))


temp <- ggplot(both_low, 
               aes(date, 
                   tmean,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(temp)

swe <- ggplot(both_low, 
               aes(date, 
                   snowTotSWE,
                   color = source)) +
  geom_line() +
  scale_color_manual(values = c("#0047AB", "orange"))

ggplotly(swe)

## combine low and high together
hi <- both_hi %>%
  dplyr::select(tmean, snowTotSWE, precip, source, date) %>%
  mutate(elevation = 'High')

low <- both_low %>%
  dplyr::select(tmean, snowTotSWE, precip, source, date) %>%
  mutate(elevation = 'Low')

both <- bind_rows(hi, low) %>%
  mutate(elev_source =
           case_when(source == 'statn' & elevation == 'High' ~ 'Station - High',
                     source == 'statn' & elevation == 'Low' ~ 'Station - Low',
                     source == 'grid' & elevation == 'High' ~ 'Gridded - High',
                     source == 'grid' & elevation == 'Low' ~ 'Gridded - Low'))

# figure
tmean <- ggplot(both,
       aes(x=date, y = tmean, 
           color = elev_source,
           linetype=elev_source)) +
   geom_line(size=1, alpha=0.8) +  
    scale_color_manual(values = c("Gridded - High" = "#0047AB", 
                                  "Gridded - Low" = "#0047AB", 
                                  "Station - High" = "orange",
                                  "Station - Low" = "orange")) +
  scale_linetype_manual(values = c("Gridded - High" = "dashed", 
                                  "Gridded - Low" = "solid", 
                                  "Station - High" = "dashed",
                                  "Station - Low" = "solid")) + 
  scale_x_date(date_labels = "%b") +
  labs(x = element_blank(),
       y = 'Mean temperature (°C)') +
  theme_bw(base_size = 20) +
  theme(legend.position = 'bottom',
    #legend.position = c(0.24, 0.852), 
    legend.text = element_text(size = 10),
    legend.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    strip.background = element_blank(),
    text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.key.width = unit(1, "cm")
    #plot.margin = margin(t = 10, r = 20, b = 10, l = 10)
  ) 

tmean

swe <- ggplot(both,
       aes(x=date, y = snowTotSWE, 
           color = elev_source,
           linetype=elev_source)) +
   geom_line(size=1, alpha=0.8) +  
    scale_color_manual(values = c("Gridded - High" = "#0047AB", 
                                  "Gridded - Low" = "#0047AB", 
                                  "Station - High" = "orange",
                                  "Station - Low" = "orange")) +
  scale_linetype_manual(values = c("Gridded - High" = "dashed", 
                                  "Gridded - Low" = "solid", 
                                  "Station - High" = "dashed",
                                  "Station - Low" = "solid")) + 
  scale_x_date(date_labels = "%b") +
  #scale_y_log10() +
  labs(x = element_blank(),
       y = 'SWE (mm)') +
  theme_bw(base_size = 20) +
  theme(
    #legend.position = c(0.89, 0.87),
    legend.position = 'bottom',
    legend.text = element_text(size = 10),
    legend.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    strip.background = element_blank(),
    text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.key.width = unit(1, "cm")
    #plot.margin = margin(t = 10, r = 20, b = 10, l = 10)
  ) 

swe

p <- both %>%
  group_by(elev_source) %>%
  mutate(cumP = cumsum(precip)) %>% 
  ggplot(.,
       aes(x=date, y = cumP, 
           color = elev_source,
           linetype=elev_source)) +
   geom_line(size=1, alpha=0.8) +  
    scale_color_manual(values = c("Gridded - High" = "#0047AB", 
                                  "Gridded - Low" = "#0047AB", 
                                  "Station - High" = "orange",
                                  "Station - Low" = "orange")) +
  scale_linetype_manual(values = c("Gridded - High" = "dashed", 
                                  "Gridded - Low" = "solid", 
                                  "Station - High" = "dashed",
                                  "Station - Low" = "solid")) + 
  scale_x_date(date_labels = "%b") +
  #scale_y_log10() +
  labs(x = element_blank(),
       y = 'Precip (mm)') +
  theme_bw(base_size = 20) +
  theme(
    legend.position = 'bottom',
    #legend.text = element_text(size = 12),
    #legend.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    strip.background = element_blank(),
    text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.key.width = unit(1, "cm")
    #plot.margin = margin(t = 10, r = 20, b = 10, l = 10)
  ) 

p

# plot_grid(tmean, swe)
# 
# ggsave('./figures/tmean_swe_2016.png',
#        width = 10,
#        height = 5,
#        dpi=600)

plot_grid(tmean, p, swe, nrow=1)

ggsave('./figures/tmean_swe_p_2016.png',
       width = 10,
       height = 5,
       dpi=600)

library(cowplot)


legend <- cowplot::get_plot_component(tmean, "guide-box", return_all = TRUE)[[3]]

tmean_noleg <- tmean + theme(legend.position = "none")
p_noleg <- p + theme(legend.position = "none")
swe_noleg <- swe + theme(legend.position = "none")

combined <- plot_grid(tmean_noleg, p_noleg, swe_noleg, nrow = 1)
final_plot <- plot_grid(combined, legend, ncol = 1, rel_heights = c(1, 0.1))
final_plot

ggsave('./figures/tmean_swe_p_2016.png',
       width = 10,
       height = 5,
       dpi=600)

# stats for precip, tmean, and swe
both_avg <- both %>%
  group_by(elev_source) %>%
  summarize(mean_temp = mean(tmean),
            max_swe = max(snowTotSWE),
            mean_p = mean(precip))

```

# Burn vs Unburn HRU

```{r}

unburn <- read_csv('./38HRU_dec2024_updates/HRU/statn_HRU_2_unburn.csv') %>%
  mutate(date = mdy(time),
         year = year(date),
         type = 'unburn') %>%
  filter(year == 2016)

burn <- read_csv('./38HRU_dec2024_updates/HRU/statn_HRU_5_burn.csv') %>%
  mutate(date = mdy(time),
         year = year(date),
         type = 'burn') %>%
  filter(year == 2016)

both <- bind_rows(unburn, burn)

# interception
inter <- both %>%
  dplyr::select(date, type, intercStorage)

int_plot <- ggplot(inter,
       aes(x=date,,
           y=intercStorage,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)

ggplotly(int_plot)

# swe
swe <- both %>%
  dplyr::select(date, type, snowTotSWE)

ggplot(swe,
       aes(x=date,,
           y=snowTotSWE,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)

# swe
p <- both %>%
  dplyr::select(date, type, precip) %>%
  group_by(type) %>%
  mutate(p_cumu = cumsum(precip))

ggplot(p,
       aes(x=date,
           y=p_cumu,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)

# soil sat
soilsat <- both %>%
  dplyr::select(date, type, soilSat)

ggplot(soilsat,
       aes(x=date,
           y=soilSat,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)


# surface runoff
runoff <- both %>%
  dplyr::select(date, type, outRD1)

ggplot(runoff,
       aes(x=date,,
           y=outRD1,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)

# interflow
interflow <- both %>%
  dplyr::select(date, type, outRD2) %>%
  mutate(interflow_cum = cumsum(outRD2))

ggplot(interflow,
       aes(x=date,,
           y=interflow_cum,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)


gw <- both %>%
  dplyr::select(date, type, outRG1) %>%
  mutate(gw_cumu = cumsum(outRG1))

ggplot(gw,
       aes(x=date,,
           y=gw_cumu,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)

# dps
dps <- both %>%
  dplyr::select(date, type, actDPS)

ggplot(dps,
       aes(x=date,,
           y=actDPS,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)

# mps
mps <- both %>%
  dplyr::select(date, type, soilSatMPS)

ggplot(mps,
       aes(x=date,,
           y=soilSatMPS,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)

# lps
lps <- both %>%
  dplyr::select(date, type, soilSatLPS)

ggplot(lps,
       aes(x=date,,
           y=soilSatLPS,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)

# soil sat
soilsat <- both %>%
  dplyr::select(date, type, soilSat)

ggplot(soilsat,
       aes(x=date,,
           y=soilSat,
           color=type)) +
  geom_line() +
  theme_bw(base_size=20)

## make the actual plot
p_plot <- ggplot(p, aes(x = date,
                        y = p_cumu, 
                        color = type)) +
  geom_line(size=1, alpha=0.8) +
    scale_color_manual(values = c("burn" = "red", 
                                  "unburn" = "black")) +
  scale_x_date(date_labels = "%b") +
  labs(x = element_blank(),
       y = 'Cumulative precipitation (mm)') +
  theme_bw(base_size = 18) +
  theme(legend.position = 'bottom',
    #legend.position = c(0.24, 0.852), 
    legend.text = element_text(size = 10),
    legend.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    strip.background = element_blank(),
    text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.key.width = unit(1, "cm")
  ) 

p_plot

inter <- inter %>%
  mutate(type = factor(type, levels = c("unburn", "burn")))

int_plot <- ggplot(inter, aes(x = date,
                        y = intercStorage, 
                        color = type)) +
  geom_line(size=1, alpha=0.7) +
    scale_color_manual(values = c("burn" = "red", 
                                  "unburn" = "black")) +
  scale_x_date(date_labels = "%b") +
  labs(x = element_blank(),
       y = 'Interception (mm)') +
  theme_bw(base_size = 20) +
  theme(legend.position = 'bottom',
    #legend.position = c(0.24, 0.852), 
    legend.text = element_text(size = 10),
    legend.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    strip.background = element_blank(),
    text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.key.width = unit(1, "cm")
  ) 

int_plot

soilsat_plot <- ggplot(soilsat, aes(x = date,
                        y = soilSat, 
                        color = type)) +
  geom_line(size=1, alpha=0.7) +
    scale_color_manual(values = c("burn" = "red", 
                                  "unburn" = "black")) +
  scale_x_date(date_labels = "%b") +
  labs(x = element_blank(),
       y = 'Soil saturation') +
  theme_bw(base_size = 20) +
  theme(legend.position = 'bottom',
    #legend.position = c(0.24, 0.852), 
    legend.text = element_text(size = 10),
    legend.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    strip.background = element_blank(),
    text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.key.width = unit(1, "cm")
  ) 

soilsat_plot

interflow_plot <- ggplot(interflow, aes(x = date,
                        y = interflow_cum, 
                        color = type)) +
  geom_line(size=1, alpha=0.7) +
    scale_color_manual(values = c("burn" = "red", 
                                  "unburn" = "black")) +
  scale_x_date(date_labels = "%b") +
  labs(x = element_blank(),
       y = 'Interflow (L)') +
  theme_bw(base_size = 20) +
  theme(legend.position = 'bottom',
    #legend.position = c(0.24, 0.852), 
    legend.text = element_text(size = 10),
    legend.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    strip.background = element_blank(),
    text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.key.width = unit(1, "cm")
  ) 

interflow_plot

gw_plot <- ggplot(gw, aes(x = date,
                        y = gw_cumu, 
                        color = type)) +
  geom_line(size=1, alpha=0.7) +
    scale_color_manual(values = c("burn" = "red", 
                                  "unburn" = "black")) +
  scale_x_date(date_labels = "%b") +
  labs(x = element_blank(),
       y = 'Shallow groundwater (L)') +
  theme_bw(base_size = 20) +
  theme(legend.position = 'bottom',
    #legend.position = c(0.24, 0.852), 
    legend.text = element_text(size = 10),
    legend.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    strip.background = element_blank(),
    text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.key.width = unit(1, "cm")
  ) 

gw_plot

swe_plot <- ggplot(swe, aes(x = date,
                        y = snowTotSWE, 
                        color = type)) +
  geom_line(size=1, alpha=0.7) +
    scale_color_manual(values = c("burn" = "red", 
                                  "unburn" = "black")) +
  scale_x_date(date_labels = "%b") +
  labs(x = element_blank(),
       y = 'SWE (mm)') +
  theme_bw(base_size = 20) +
  theme(legend.position = 'bottom',
    #legend.position = c(0.24, 0.852), 
    legend.text = element_text(size = 10),
    legend.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    strip.background = element_blank(),
    text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.key.width = unit(1, "cm")
  ) 

swe_plot

legend <- cowplot::get_plot_component(p_plot, "guide-box", return_all = TRUE)[[3]]

p_plot <- p_plot + theme(legend.position = "none")
int_plot <- int_plot + theme(legend.position = "none")
soilsat_plot <- soilsat_plot + theme(legend.position = "none")
interflow_plot <- interflow_plot + theme(legend.position = "none")
gw_plot <- gw_plot + theme(legend.position = "none")
swe_plot <- swe_plot + theme(legend.position = "none")


plots <- plot_grid(p_plot,
          int_plot,
          swe_plot,
          soilsat_plot,
          interflow_plot,
          gw_plot,
          ncol = 3)

final_plot <- plot_grid(
  plots,
  legend,
  ncol = 1,
  rel_heights = c(1, 0.07)  # Adjust height of legend
)

final_plot

ggsave('./figures/burn_unburn_eg.png')

soilsat_mean <- soilsat %>% 
  group_by(type) %>% 
  summarize(mean_ss = mean(soilSat))

soilsat_mean

intercep_mean <- inter %>% 
  group_by(type) %>% 
  summarize(mean_inter = mean(intercStorage))

intercep_mean

```