workflow.Rmd

---
title: "High Discharge Event Workflow"
author: "Taryn Waite"
date: "8/9/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(magrittr)
library(tidyverse)
library(lubridate)
library(dataRetrieval)
library(ggsci)
library(ggthemes)

source("event_functions.R")
source("C_Q_plotting_functions.R")
source("event_stats_functions.R")
```

This document contains the entire high discharge event (HDE) workflow, including reading in the discharge, water quality, and precipitation data, performing data wrangling, identifying HDEs, calculating slope and hysteresis index for each constituent, and statistically analyzing the resulting data. 

# Data retrieval and wrangling

Reading in the data
```{r data retrieval}
# main channel WQ data
mc <- read_csv("data/MC_WQ.csv")
# backwater WQ data
si <- read_csv("data/SI_WQ.csv")

# discharge data
unitDischarge <- readNWISuv(siteNumbers = "05378500", parameterCd="00060", startDate = "2015-04-01", endDate = "2018-10-31") %>% 
  mutate(discharge = X_00060_00000)

# precip data from winona, mn 
precip2 <- read_csv("data/daily_precip_winona.csv") %>% 
  rename(precip = PRCP, snow = SNOW) %>% 
  mutate(precip = as.numeric(precip), snow = as.numeric(snow),
         Date = mdy(Date))

```

WQ data wrangling: Converting the date column into Date objects and also creating a date-time column
```{r WQ date and time}
mc <- mc %>% 
  mutate(date = mdy(`YYYY-MM-DD`)) %>% 
  mutate(dateTime = mdy_hms(paste(`YYYY-MM-DD`, as.character(`hh:mm`))))

si <- si %>% 
  mutate(date = ymd(`YYYY-MM-DD`)) %>% 
  mutate(dateTime = ymd_hms(paste(`YYYY-MM-DD`, as.character(`hh:mm`))))
```


Combining the SI and MC data into one dataframe with a column called "site" taking the value "SI" or "MC".
```{r combining si and mc}
# first we need a new column specifying if si or mc
si_with_ID <- si %>% mutate(site = "SI")
mc_with_ID <- mc %>% mutate(site = "MC")

# get rid of original date and time columns
# (don't need these anymore with new date and dateTime columns)
si_with_ID <- si_with_ID %>% select(-c(1,2, 23))
mc_with_ID <- mc_with_ID %>% select(-c(1,2, 23))

# next we can merge the two data frames by row
all_WQ <- bind_rows(si_with_ID, mc_with_ID)
```

The WQ data is at one- or two- hour frequency (depending on the year), so we need to summarize the 15-minute discharge data by hour.
```{r discharge data wrangling}
# add a date column and an hour column to the discharge dataframe
unitDischarge2 <- unitDischarge %>% 
  mutate(date = date(dateTime), hour = hour(dateTime))

# create a mean hourly discharge dataframe
discharge_mean_hourly <- unitDischarge2 %>% 
  group_by(date, hour) %>% 
  summarize(hourlyDischarge = mean(X_00060_00000))
```

Finally, we can merge the discharge data with the water quality data, joining the dataframes by hour.
```{r merging discharge and WQ}
# add an hour column to the WQ dataframe
all_WQ_hour <- all_WQ %>% 
  mutate(hour = hour(dateTime))


# merge discharge and WQ for all years
WQ_hourly_discharge <- full_join(all_WQ_hour, discharge_mean_hourly, by = c("date", "hour"))

```

```{r constituent summaries}
##### MC 2015 #####
study_2015_MC <- WQ_hourly_discharge %>% 
  filter(site == "MC",
         date >= mdy("06-02-2015"),
         date <= mdy("10-21-2015"))

mean(study_2015_MC$Turb, na.rm = T)
sd(study_2015_MC$Turb, na.rm = T)

mean(study_2015_MC$CHLugL, na.rm = T)
sd(study_2015_MC$CHLugL, na.rm = T)

mean(study_2015_MC$BGAugL, na.rm = T)
sd(study_2015_MC$BGAugL, na.rm = T)

mean(study_2015_MC$FDOMqsu, na.rm = T)
sd(study_2015_MC$FDOMqsu, na.rm = T)

mean(study_2015_MC$NO3_mgL, na.rm = T)
sd(study_2015_MC$NO3_mgL, na.rm = T)

##### MC 2016 #####
study_2016_MC <- WQ_hourly_discharge %>% 
  filter(site == "MC",
         date >= mdy("04-13-2016"),
         date <= mdy("10-31-2016"))

mean(study_2016_MC$Turb, na.rm = T)
sd(study_2016_MC$Turb, na.rm = T)

mean(study_2016_MC$CHLugL, na.rm = T)
sd(study_2016_MC$CHLugL, na.rm = T)

mean(study_2016_MC$BGAugL, na.rm = T)
sd(study_2016_MC$BGAugL, na.rm = T)

mean(study_2016_MC$FDOMqsu, na.rm = T)
sd(study_2016_MC$FDOMqsu, na.rm = T)

mean(study_2016_MC$NO3_mgL, na.rm = T)
sd(study_2016_MC$NO3_mgL, na.rm = T)

##### MC 2017 #####
study_2017_MC <- WQ_hourly_discharge %>% 
  filter(site == "MC",
         date >= mdy("06-09-2017"),
         date <= mdy("10-26-2017"))

mean(study_2017_MC$Turb, na.rm = T)
sd(study_2017_MC$Turb, na.rm = T)

mean(study_2017_MC$CHLugL, na.rm = T)
sd(study_2017_MC$CHLugL, na.rm = T)

mean(study_2017_MC$BGAugL, na.rm = T)
sd(study_2017_MC$BGAugL, na.rm = T)

mean(study_2017_MC$FDOMqsu, na.rm = T)
sd(study_2017_MC$FDOMqsu, na.rm = T)

mean(study_2017_MC$NO3_mgL, na.rm = T)
sd(study_2017_MC$NO3_mgL, na.rm = T)

##### MC 2018 #####
study_2018_MC <- WQ_hourly_discharge %>% 
  filter(site == "MC",
         date >= mdy("05-08-2018"),
         date <= mdy("10-12-2018"))

mean(study_2018_MC$Turb, na.rm = T)
sd(study_2018_MC$Turb, na.rm = T)

mean(study_2018_MC$CHLugL, na.rm = T)
sd(study_2018_MC$CHLugL, na.rm = T)

mean(study_2018_MC$BGAugL, na.rm = T)
sd(study_2018_MC$BGAugL, na.rm = T)

mean(study_2018_MC$FDOMqsu, na.rm = T)
sd(study_2018_MC$FDOMqsu, na.rm = T)

mean(study_2018_MC$NO3_mgL, na.rm = T)
sd(study_2018_MC$NO3_mgL, na.rm = T)


##### SI 2015 #####
study_2015_SI <- WQ_hourly_discharge %>% 
  filter(site == "SI",
         date >= mdy("06-02-2015"),
         date <= mdy("10-21-2015"))

mean(study_2015_SI$Turb, na.rm = T)
sd(study_2015_SI$Turb, na.rm = T)

mean(study_2015_SI$CHLugL, na.rm = T)
sd(study_2015_SI$CHLugL, na.rm = T)

mean(study_2015_SI$BGAugL, na.rm = T)
sd(study_2015_SI$BGAugL, na.rm = T)

mean(study_2015_SI$FDOMqsu, na.rm = T)
sd(study_2015_SI$FDOMqsu, na.rm = T)

mean(study_2015_SI$NO3_mgL, na.rm = T)
sd(study_2015_SI$NO3_mgL, na.rm = T)

##### SI 2016 #####
study_2016_SI <- WQ_hourly_discharge %>% 
  filter(site == "SI",
         date >= mdy("04-13-2016"),
         date <= mdy("10-31-2016"))

mean(study_2016_SI$Turb, na.rm = T)
sd(study_2016_SI$Turb, na.rm = T)

mean(study_2016_SI$CHLugL, na.rm = T)
sd(study_2016_SI$CHLugL, na.rm = T)

mean(study_2016_SI$BGAugL, na.rm = T)
sd(study_2016_SI$BGAugL, na.rm = T)

mean(study_2016_SI$FDOMqsu, na.rm = T)
sd(study_2016_SI$FDOMqsu, na.rm = T)

mean(study_2016_SI$NO3_mgL, na.rm = T)
sd(study_2016_SI$NO3_mgL, na.rm = T)

##### SI 2017 #####
study_2017_SI <- WQ_hourly_discharge %>% 
  filter(site == "SI",
         date >= mdy("06-09-2017"),
         date <= mdy("10-26-2017"))

mean(study_2017_SI$Turb, na.rm = T)
sd(study_2017_SI$Turb, na.rm = T)

mean(study_2017_SI$CHLugL, na.rm = T)
sd(study_2017_SI$CHLugL, na.rm = T)

mean(study_2017_SI$BGAugL, na.rm = T)
sd(study_2017_SI$BGAugL, na.rm = T)

mean(study_2017_SI$FDOMqsu, na.rm = T)
sd(study_2017_SI$FDOMqsu, na.rm = T)

mean(study_2017_SI$NO3_mgL, na.rm = T)
sd(study_2017_SI$NO3_mgL, na.rm = T)

##### MC 2018 #####
study_2018_SI <- WQ_hourly_discharge %>% 
  filter(site == "SI",
         date >= mdy("05-08-2018"),
         date <= mdy("10-12-2018"))

mean(study_2018_SI$Turb, na.rm = T)
sd(study_2018_SI$Turb, na.rm = T)

mean(study_2018_SI$CHLugL, na.rm = T)
sd(study_2018_SI$CHLugL, na.rm = T)

mean(study_2018_SI$BGAugL, na.rm = T)
sd(study_2018_SI$BGAugL, na.rm = T)

mean(study_2018_SI$FDOMqsu, na.rm = T)
sd(study_2018_SI$FDOMqsu, na.rm = T)

mean(study_2018_SI$NO3_mgL, na.rm = T)
sd(study_2018_SI$NO3_mgL, na.rm = T)


##### comparing SI and MC #####

# nitrate sig higher in MC than SI all 4 years
t.test(study_2015_MC$NO3_mgL, study_2015_SI$NO3_mgL)
t.test(study_2016_MC$NO3_mgL, study_2016_SI$NO3_mgL)
t.test(study_2017_MC$NO3_mgL, study_2017_SI$NO3_mgL)
t.test(study_2018_MC$NO3_mgL, study_2018_SI$NO3_mgL)

# turb sig higher in MC than SI all 4 years
t.test(study_2015_MC$Turb, study_2015_SI$Turb)
t.test(study_2016_MC$Turb, study_2016_SI$Turb)
t.test(study_2017_MC$Turb, study_2017_SI$Turb)
t.test(study_2018_MC$Turb, study_2018_SI$Turb)

# fdom sig higher in SI than MC in 2015, opposite other 3 years
t.test(study_2015_MC$FDOMqsu, study_2015_SI$FDOMqsu)
t.test(study_2016_MC$FDOMqsu, study_2016_SI$FDOMqsu)
t.test(study_2017_MC$FDOMqsu, study_2017_SI$FDOMqsu)
t.test(study_2018_MC$FDOMqsu, study_2018_SI$FDOMqsu)

# chl sig higher in MC than SI all 4 years
t.test(study_2015_MC$CHLugL, study_2015_SI$CHLugL)
t.test(study_2016_MC$CHLugL, study_2016_SI$CHLugL)
t.test(study_2017_MC$CHLugL, study_2017_SI$CHLugL)
t.test(study_2018_MC$CHLugL, study_2018_SI$CHLugL)

# bga sig higher in MC than SI all 4 years
t.test(study_2015_MC$BGAugL, study_2015_SI$BGAugL)
t.test(study_2016_MC$BGAugL, study_2016_SI$BGAugL)
t.test(study_2017_MC$BGAugL, study_2017_SI$BGAugL)
t.test(study_2018_MC$BGAugL, study_2018_SI$BGAugL)

##### comparing years #####

# bga MC: 2016, 2015, 2017, 2018
t.test(study_2015_MC$BGAugL, study_2016_MC$BGAugL) # 2015 > 2016
t.test(study_2015_MC$BGAugL, study_2017_MC$BGAugL) # 2017 > 2015
t.test(study_2015_MC$BGAugL, study_2018_MC$BGAugL) # 2018 > 2015
t.test(study_2016_MC$BGAugL, study_2017_MC$BGAugL) # 2017 > 2016
t.test(study_2016_MC$BGAugL, study_2018_MC$BGAugL) # 2018 > 2016
t.test(study_2017_MC$BGAugL, study_2018_MC$BGAugL) # 2018 > 2017

# bga SI: 2015/2017, 2016, 2018
t.test(study_2015_SI$BGAugL, study_2016_SI$BGAugL) # 2016 > 2015
t.test(study_2015_SI$BGAugL, study_2017_SI$BGAugL) # 2017 = 2015
t.test(study_2015_SI$BGAugL, study_2018_SI$BGAugL) # 2018 > 2015
t.test(study_2016_SI$BGAugL, study_2017_SI$BGAugL) # 2016 > 2017
t.test(study_2016_SI$BGAugL, study_2018_SI$BGAugL) # 2018 > 2016
t.test(study_2017_SI$BGAugL, study_2018_SI$BGAugL) # 2018 > 2017

# chl MC: 2016, 2018, 2017, 2015
t.test(study_2015_MC$CHLugL, study_2016_MC$CHLugL) # 2015 > 2016
t.test(study_2015_MC$CHLugL, study_2017_MC$CHLugL) # 2015 > 2017
t.test(study_2015_MC$CHLugL, study_2018_MC$CHLugL) # 2015 > 2018
t.test(study_2016_MC$CHLugL, study_2017_MC$CHLugL) # 2017 > 2016
t.test(study_2016_MC$CHLugL, study_2018_MC$CHLugL) # 2018 > 2016
t.test(study_2017_MC$CHLugL, study_2018_MC$CHLugL) # 2017 > 2018

# chl SI: 2015, 2017, 2018, 2016
t.test(study_2015_SI$CHLugL, study_2016_SI$CHLugL) # 2016 > 2015
t.test(study_2015_SI$CHLugL, study_2017_SI$CHLugL) # 2017 > 2015
t.test(study_2015_SI$CHLugL, study_2018_SI$CHLugL) # 2018 > 2015
t.test(study_2016_SI$CHLugL, study_2017_SI$CHLugL) # 2016 > 2017
t.test(study_2016_SI$CHLugL, study_2018_SI$CHLugL) # 2016 > 2018
t.test(study_2017_SI$CHLugL, study_2018_SI$CHLugL) # 2018 > 2017

# fdom MC: 2015, 2017, 2016/2018
t.test(study_2015_MC$FDOMqsu, study_2016_MC$FDOMqsu) # 2016 > 2015
t.test(study_2015_MC$FDOMqsu, study_2017_MC$FDOMqsu) # 2017 > 2015
t.test(study_2015_MC$FDOMqsu, study_2018_MC$FDOMqsu) # 2018 > 2015
t.test(study_2016_MC$FDOMqsu, study_2017_MC$FDOMqsu) # 2016 > 2017
t.test(study_2016_MC$FDOMqsu, study_2018_MC$FDOMqsu) # 2018 = 2016
t.test(study_2017_MC$FDOMqsu, study_2018_MC$FDOMqsu) # 2018 > 2017

# fdom SI: 2015, 2017, 2018, 2016
t.test(study_2015_SI$FDOMqsu, study_2016_SI$FDOMqsu) # 2016 > 2015
t.test(study_2015_SI$FDOMqsu, study_2017_SI$FDOMqsu) # 2017 > 2015
t.test(study_2015_SI$FDOMqsu, study_2018_SI$FDOMqsu) # 2018 > 2015
t.test(study_2016_SI$FDOMqsu, study_2017_SI$FDOMqsu) # 2016 > 2017
t.test(study_2016_SI$FDOMqsu, study_2018_SI$FDOMqsu) # 2016 > 2018
t.test(study_2017_SI$FDOMqsu, study_2018_SI$FDOMqsu) # 2018 > 2017

# nO3 MC: 2015, 2017, 2018, 2016
t.test(study_2015_MC$NO3_mgL, study_2016_MC$NO3_mgL) # 2016 > 2015
t.test(study_2015_MC$NO3_mgL, study_2017_MC$NO3_mgL) # 2017 > 2015
t.test(study_2015_MC$NO3_mgL, study_2018_MC$NO3_mgL) # 2018 > 2015
t.test(study_2016_MC$NO3_mgL, study_2017_MC$NO3_mgL) # 2016 > 2017
t.test(study_2016_MC$NO3_mgL, study_2018_MC$NO3_mgL) # 2016 > 2018
t.test(study_2017_MC$NO3_mgL, study_2018_MC$NO3_mgL) # 2018 > 2017

# nO3 SI: 2015, 2017, 2018, 2016
t.test(study_2015_SI$NO3_mgL, study_2016_SI$NO3_mgL) # 2016 > 2015
t.test(study_2015_SI$NO3_mgL, study_2017_SI$NO3_mgL) # 2017 > 2015
t.test(study_2015_SI$NO3_mgL, study_2018_SI$NO3_mgL) # 2018 > 2015
t.test(study_2016_SI$NO3_mgL, study_2017_SI$NO3_mgL) # 2016 > 2017
t.test(study_2016_SI$NO3_mgL, study_2018_SI$NO3_mgL) # 2016 > 2018
t.test(study_2017_SI$NO3_mgL, study_2018_SI$NO3_mgL) # 2018 > 2017

# turb MC: 2015, 2018, 2016, 2017
t.test(study_2015_MC$Turb, study_2016_MC$Turb) # 2016 > 2015
t.test(study_2015_MC$Turb, study_2017_MC$Turb) # 2017 > 2015
t.test(study_2015_MC$Turb, study_2018_MC$Turb) # 2018 > 2015
t.test(study_2016_MC$Turb, study_2017_MC$Turb) # 2017 > 2016
t.test(study_2016_MC$Turb, study_2018_MC$Turb) # 2016 > 2018
t.test(study_2017_MC$Turb, study_2018_MC$Turb) # 2017 > 2018

# turb SI: 2015, 2016, 2018, 2017
t.test(study_2015_SI$Turb, study_2016_SI$Turb) # 2016 > 2015
t.test(study_2015_SI$Turb, study_2017_SI$Turb) # 2017 > 2015
t.test(study_2015_SI$Turb, study_2018_SI$Turb) # 2018 > 2015
t.test(study_2016_SI$Turb, study_2017_SI$Turb) # 2017 > 2016
t.test(study_2016_SI$Turb, study_2018_SI$Turb) # 2018 > 2016
t.test(study_2017_SI$Turb, study_2018_SI$Turb) # 2017 > 2018


study_all <- WQ_hourly_discharge %>% 
  filter(date >= mdy("06-02-2015") & date <= mdy("10-21-2015") |
         date >= mdy("04-13-2016") & date <= mdy("10-31-2016") |
         date >= mdy("06-09-2017") & date <= mdy("10-26-2017") |
         date >= mdy("05-08-2018") & date <= mdy("10-12-2018") )

bga_box <- study_all %>% 
           drop_na() %>% 
           ggplot(aes(x = factor(year(dateTime)), y = BGAugL, col = site)) +
           geom_boxplot()

chl_box <- study_all %>% 
           drop_na() %>% 
           ggplot(aes(x = factor(year(dateTime)), y = CHLugL, col = site)) +
           geom_boxplot()

fdom_box <- study_all %>% 
            drop_na() %>% 
            ggplot(aes(x = factor(year(dateTime)), y = FDOMqsu, col = site)) +
            geom_boxplot()

no3_box <- study_all %>% 
           drop_na() %>% 
           ggplot(aes(x = factor(year(dateTime)), y = NO3_mgL, col = site)) +
           geom_boxplot()

turb_box <- study_all %>% 
            drop_na() %>% 
            ggplot(aes(x = factor(year(dateTime)), y = Turb, col = site)) +
            geom_boxplot()

ggarrange(bga_box, chl_box, fdom_box, no3_box, turb_box, 
          nrow = 3, ncol = 2, common.legend = T)

```

```{r hydrographs}
check_during_event <- function(date){
  for (i in 1:18){
    if(date >= date(eventsTurb$startDate[i]) &
       date <= date(eventsTurb$endDate[i])) {
      return(T)
       }
  }
  return(F)
}

WQ_hourly_discharge_startDates <- WQ_hourly_discharge %>% 
  mutate (startDate = date(dateTime) %in% date(eventsFinal2$startDate),
          during_event = sapply(date,check_during_event))

just_starts <- WQ_hourly_discharge_startDates %>% 
  filter(startDate, site == "MC", hour %in% c(0,1))

just_events <- WQ_hourly_discharge_startDates %>% filter(during_event) %>% 
  mutate(hourlyDischarge = case_when(date %in% just_starts$date ~ -1,
                                     TRUE ~ hourlyDischarge))

just_events$hourlyDischarge <- replace(just_events$hourlyDischarge, 
                                       which(just_events$hourlyDischarge < 0), NA)

library(patchwork)
precipHydro <- precip2 %>% filter(month(Date) %in% c(4, 5, 6, 7, 8, 9, 10)) %>% 
  mutate(date = Date, metPrecip = precip*2.54) 

png(filename = "hydrographsPrecip.png", width = 6, height = 4, 
    units = "in", res = 300, bg = "transparent")
WQ_hourly_discharge_startDates %>% filter(month(date) %in% c(4, 5, 6, 7, 8, 9, 10)) %>% 
  ggplot(aes(x = date, y = hourlyDischarge*0.02831685)) +
  scale_y_continuous(sec.axis = sec_axis(trans=~.*.003, name="Daily Precipitation (cm)")) +
  geom_bar(data = precipHydro, aes(x = Date, y = metPrecip/.003), stat = "identity",fill = "grey59") +
  geom_line(col = "black") +
  geom_line(data = just_events, aes(x = date, y = hourlyDischarge*0.02831685), col = "red") +
  xlab("Time") + ylab(expression(paste("Discharge (m"^{3} ,"s"^{-1}, ")"))) +
  geom_point(data = just_starts, 
             aes(x = date, y = hourlyDischarge*0.02831685),
             col = "blue", cex = 1) +
  geom_vline(xintercept = mdy("06-02-2015"), lty = 2) +
  geom_vline(xintercept = mdy("10-21-2015"), lty = 2) +
  geom_vline(xintercept = mdy("04-13-2016"), lty = 2) +
  geom_vline(xintercept = mdy("10-31-2016"), lty = 2) +
  geom_vline(xintercept = mdy("06-09-2017"), lty = 2) +
  geom_vline(xintercept = mdy("10-26-2017"), lty = 2) +
  geom_vline(xintercept = mdy("05-08-2018"), lty = 2) +
  geom_vline(xintercept = mdy("10-12-2018"), lty = 2) +
  facet_wrap(vars(as.factor(year(date))), scales = "free_x") +
  theme_bw() +
  theme(plot.background = element_rect(fill = "transparent",colour = NA),
        legend.background = element_rect(fill = "transparent"), legend.position = "bottom") 
dev.off()

png(filename = "hydrographsLong.png", width = 10, height = 2, 
    units = "in", res = 300, bg = "transparent")
WQ_hourly_discharge_startDates %>% filter(month(date) %in% c(4, 5, 6, 7, 8, 9, 10)) %>% 
  ggplot(aes(x = date, y = hourlyDischarge)) +
  geom_line(col = "black") +
  geom_line(data = just_events, aes(x = date, y = hourlyDischarge), col = "red") +
  xlab("Time") + ylab(expression(paste("Discharge (m"^{3} ,"s"^{-1}, ")"))) +
  geom_point(data = just_starts, 
             aes(x = date, y = hourlyDischarge),
             col = "blue", cex = 1) +
  geom_vline(xintercept = mdy("06-02-2015"), lty = 2) +
  geom_vline(xintercept = mdy("10-21-2015"), lty = 2) +
  geom_vline(xintercept = mdy("04-13-2016"), lty = 2) +
  geom_vline(xintercept = mdy("10-31-2016"), lty = 2) +
  geom_vline(xintercept = mdy("06-09-2017"), lty = 2) +
  geom_vline(xintercept = mdy("10-26-2017"), lty = 2) +
  geom_vline(xintercept = mdy("05-08-2018"), lty = 2) +
  geom_vline(xintercept = mdy("10-12-2018"), lty = 2) +
  facet_wrap(vars(as.factor(year(date))), scales = "free_x", nrow = 1) +
  theme_bw() +
  theme(plot.background = element_rect(fill = "transparent",colour = NA),
        legend.background = element_rect(fill = "transparent"), legend.position = "bottom") 
dev.off()



```

# High Discharge Event identification
Identifying high discharge events involves a combination of visual observation of the hydrographs and quantitative analysis. I decided to use a daily rate of change threshold to flag the beginning of an event, since this does a good job of capturing the starts of HDEs that are visible in the hydrograph. Through a process of testing different thresholds and observing the points captured, I settled on a threshold of 3000 cubic feet per second per day. This means that a day is flagged as the start of an event if the difference in average discharge between the day and the next day is greater than 3000 cubic feet per second, but the difference between the day and the previous day is less than this threshold.

First, we add a difference column to the daily discharge data (the difference between this day and the next day)
```{r adding discharge difference column}
# arrange by date
temp <- discharge %>% 
  rename(daily_discharge = X_00060_00003) %>% 
  arrange(Date) 


# calculate difference column
temp$diff <- lead(temp$daily_discharge, 1) - temp$daily_discharge
```

Then, we flag dates that are above the threshold and for which the previous date was below the threshold.
```{r flagging each start date}
# flag whenever the rate of change is above the threshold
threshold <- 3000
storm2 <- temp %>% 
  mutate(status = case_when(diff > threshold ~ "storm", T ~"base"))

# function to test whether previous day was over threshold
prevThresh <- function(date){
  prevDay <- storm2 %>% filter(Date == ((date)-1))
  #stat <- prevDay[["status"]]
  if(nrow(prevDay) == 0) {
    return (F)
  }
  if (prevDay$status == "storm"){
    return (T)
  } else {
    return (F)
  }
}

# flag where the threshold is crossed 
# (when a point is above the threshold and the previous point is not)
storm3 <- storm2 %>% 
  mutate(prevDay = sapply(Date, prevThresh)) %>% 
  mutate(crossed = !prevDay & status == "storm")

# get a df of just the storm starts
starts <- storm3 %>% filter(crossed == T & month(Date) %in% c(4, 5, 6, 7, 8, 9, 10))
```

# Missing Data

Next, we need to fill in some missing values for the constituent time series. For time series with values missing in just a single row, we simply use the average of the previous and next values.
```{r}
# SI
WQ_hourly_discharge[which(WQ_hourly_discharge$dateTime==ymd_hms("2015-09-03 12:05:00")),] %<>%
  mutate(NO3_mgL = reassign(WQ_hourly_discharge, "NO3_mgL", ymd_hms("2015-09-03 12:05:00")),
         Turb = reassign(WQ_hourly_discharge, "Turb", ymd_hms("2015-09-03 12:05:00")),
         CHLugL = reassign(WQ_hourly_discharge, "CHLugL", ymd_hms("2015-09-03 12:05:00")),
         FDOMqsu = reassign(WQ_hourly_discharge, "FDOMqsu", ymd_hms("2015-09-03 12:05:00")),
         BGAugL = reassign(WQ_hourly_discharge, "BGAugL", ymd_hms("2015-09-03 12:05:00")),
         Temp = reassign(WQ_hourly_discharge, "Temp", ymd_hms("2015-09-03 12:05:00")))

WQ_hourly_discharge[which(WQ_hourly_discharge$dateTime==ymd_hms("2016-08-18 13:28:00")),] %<>%
  mutate(NO3_mgL = reassign(WQ_hourly_discharge, "NO3_mgL", ymd_hms("2016-08-18 13:28:00")),
         Turb = reassign(WQ_hourly_discharge, "Turb", ymd_hms("2016-08-18 13:28:00")),
         CHLugL = reassign(WQ_hourly_discharge, "CHLugL", ymd_hms("2016-08-18 13:28:00")),
         FDOMqsu = reassign(WQ_hourly_discharge, "FDOMqsu", ymd_hms("2016-08-18 13:28:00")),
         BGAugL = reassign(WQ_hourly_discharge, "BGAugL", ymd_hms("2016-08-18 13:28:00")),
         Temp = reassign(WQ_hourly_discharge, "Temp", ymd_hms("2016-08-16 23:28:00")))

WQ_hourly_discharge[which(WQ_hourly_discharge$dateTime==ymd_hms("2016-08-16 23:28:00")),] %<>%
  mutate(Turb = reassign(WQ_hourly_discharge, "Turb", ymd_hms("2016-08-16 23:28:00")))



# main channel
WQ_hourly_discharge[which(WQ_hourly_discharge$dateTime==ymd_hms("2017-09-12 13:13:00")),] %<>%
  mutate(NO3_mgL = reassign(WQ_hourly_discharge, "NO3_mgL", ymd_hms("2017-09-12 13:13:00")),
         Turb = reassign(WQ_hourly_discharge, "Turb", ymd_hms("2017-09-12 13:13:00")),
         CHLugL = reassign(WQ_hourly_discharge, "CHLugL", ymd_hms("2017-09-12 13:13:00")),
         FDOMqsu = reassign(WQ_hourly_discharge, "FDOMqsu", ymd_hms("2017-09-12 13:13:00")),
         BGAugL = reassign(WQ_hourly_discharge, "BGAugL", ymd_hms("2017-09-12 13:13:00")),
         Temp = reassign(WQ_hourly_discharge, "Temp", ymd_hms("2017-09-12 13:13:00")))

WQ_hourly_discharge[which(WQ_hourly_discharge$dateTime==ymd_hms("2016-05-16 10:34:00")),] %<>%
  mutate(NO3_mgL = reassign(WQ_hourly_discharge, "NO3_mgL", ymd_hms("2016-05-16 10:34:00")),
         Turb = reassign(WQ_hourly_discharge, "Turb", ymd_hms("2016-05-16 10:34:00")),
         CHLugL = reassign(WQ_hourly_discharge, "CHLugL", ymd_hms("2016-05-16 10:34:00")),
         FDOMqsu = reassign(WQ_hourly_discharge, "FDOMqsu", ymd_hms("2016-05-16 10:34:00")),
         BGAugL = reassign(WQ_hourly_discharge, "BGAugL", ymd_hms("2016-05-16 10:34:00")),
         Temp = reassign(WQ_hourly_discharge, "Temp", ymd_hms("2016-05-16 10:34:00")))

WQ_hourly_discharge[which(WQ_hourly_discharge$dateTime==ymd_hms("2016-10-11 10:49:00")),] %<>%
  mutate(NO3_mgL = reassign(WQ_hourly_discharge, "NO3_mgL", ymd_hms("2016-10-11 10:49:00")),
         Turb = reassign(WQ_hourly_discharge, "Turb", ymd_hms("2016-10-11 10:49:00")),
         CHLugL = reassign(WQ_hourly_discharge, "CHLugL", ymd_hms("2016-10-11 10:49:00")),
         FDOMqsu = reassign(WQ_hourly_discharge, "FDOMqsu", ymd_hms("2016-10-11 10:49:00")),
         BGAugL = reassign(WQ_hourly_discharge, "BGAugL", ymd_hms("2016-10-11 10:49:00")),
         Temp = reassign(WQ_hourly_discharge, "Temp", ymd_hms("2016-10-11 10:49:00")))

WQ_hourly_discharge[which(WQ_hourly_discharge$dateTime==ymd_hms("2018-09-06 11:47:00")),] %<>%
  mutate(NO3_mgL = reassign(WQ_hourly_discharge, "NO3_mgL", ymd_hms("2018-09-06 11:47:00")),
         Turb = reassign(WQ_hourly_discharge, "Turb", ymd_hms("2018-09-06 11:47:00")),
         CHLugL = reassign(WQ_hourly_discharge, "CHLugL", ymd_hms("2018-09-06 11:47:00")),
         FDOMqsu = reassign(WQ_hourly_discharge, "FDOMqsu", ymd_hms("2018-09-06 11:47:00")),
         BGAugL = reassign(WQ_hourly_discharge, "BGAugL", ymd_hms("2018-09-06 11:47:00")),
         Temp = reassign(WQ_hourly_discharge, "Temp", ymd_hms("2018-09-06 11:47:00")))
```

# Generating event data
First, we need to create a data structure to hold all of the event data, including the event metrics and the hysteresis index, slope, and area for each constituent.
```{r}
# initialize the tibble
events <- tibble(startDate = Date(), endDate = Date(), 
                 length = numeric(), length_loop = numeric(),
                 dis_change = numeric(), dis_change_loop = numeric(),
                 max_dis = numeric(), avg_temp = numeric(),
                 HI_turb = numeric(), Slope_turb = numeric(), area_turb = numeric(),
                 HI_n03 = numeric(), Slope_n03 = numeric(), area_n03 = numeric(),
                 HI_chl = numeric(), Slope_chl = numeric(), area_chl = numeric(),
                 HI_fdom = numeric(), Slope_fdom = numeric(), area_fdom = numeric(),
                 HI_bga = numeric(), Slope_bga = numeric(), area_bga = numeric())

# adds a high-discharge event to the events tibble,
# given the data, start date, end date, and number of
# observations per day
add.event <- function(data, Site, start, end, freq){
  # subset the data to include just the event
  subset <- data %>% 
    filter(dateTime >= ymd_hm(start) & dateTime <= ymd_hm(end)) %>% 
    filter(site == "MC")
  # add rising or falling limb column
  peakDate <- (subset %>% slice_max(hourlyDischarge))[["dateTime"]]
  # RL if before peakDate, FL if after, peak if on peakDate
  subset <- subset %>% 
    mutate(limb = case_when(dateTime < peakDate ~ "RL",
                            dateTime > peakDate ~ "FL",
                            T ~ "peak"))
  turb <- storm_cq(data, Site, start, end, "Turb", 20, freq)
  n03 <- storm_cq(data, Site, start, end, "NO3_mgL", 20, freq)
  chl <- storm_cq(data, Site, start, end, "CHLugL", 20, freq)
  fdom <- storm_cq(data, Site, start, end, "FDOMqsu", 20, freq)
  bga <- storm_cq(data, Site, start, end, "BGAugL", 20, freq)
  return(events %>% add_row(startDate = ymd_hm(start), endDate = ymd_hm(end),
                     length = as.numeric(ymd_hm(end) - ymd_hm(start)), 
                     length_loop = storm.length.loop(subset),
                     dis_change = storm.dis.change(subset),
                     dis_change_loop = storm.dis.change.loop(subset),
                     max_dis = storm.max.dis(subset),
                     avg_temp = storm.avg.temp(subset),
                     HI_turb = turb[[2]], Slope_turb = turb[[4]], area_turb = turb[[3]],
                     HI_n03 = n03[[2]], Slope_n03 = n03[[4]], area_n03 = n03[[3]],
                     HI_chl = chl[[2]], Slope_chl = chl[[4]], area_chl = chl[[3]],
                     HI_fdom = fdom[[2]], Slope_fdom = fdom[[4]], area_fdom = fdom[[3]],
                     HI_bga = bga[[2]], Slope_bga = bga[[4]], area_bga = bga[[3]]))
}
```

Now, we can add all of the events to the initialized tibble, using the function defined above. This takes a long time to run, so we write the tibble to a csv once it's done so that we can read the csv back in instead of re-running it.
```{r, eval=F}
# add events -- this takes a long time to run 
events <- add.event(WQ_hourly_discharge, "MC", "2016-04-25-00-00", "2016-05-24-00-00", 12)
events <- add.event(WQ_hourly_discharge, "MC", "2016-07-13-00-00", "2016-08-10-00-00", 12)
events <- add.event(WQ_hourly_discharge, "MC", "2016-06-01-00-00", "2016-06-14-00-00", 12)
events <- add.event(WQ_hourly_discharge, "MC", "2016-08-10-00-00", "2016-09-05-00-00", 12)
events <- add.event(WQ_hourly_discharge, "MC", "2016-09-06-00-00", "2016-09-20-00-00", 12)
events <- add.event(WQ_hourly_discharge, "MC", "2016-09-20-00-00", "2016-10-15-00-00", 12)
events <- add.event(WQ_hourly_discharge, "MC", "2017-06-14-00-00", "2017-06-28-00-00", 12)
events <- add.event(WQ_hourly_discharge, "MC", "2017-07-19-00-00", "2017-07-25-00-00", 24)
events <- add.event(WQ_hourly_discharge, "MC", "2017-07-26-00-00", "2017-07-30-00-00", 24)
events <- add.event(WQ_hourly_discharge, "MC", "2017-08-06-00-00", "2017-08-13-00-00", 24)
events <- add.event(WQ_hourly_discharge, "MC", "2017-09-20-00-00", "2017-09-24-00-00", 24)
events <- add.event(WQ_hourly_discharge, "MC", "2018-06-15-00-00", "2018-07-31-00-00", 12)
events <- add.event(WQ_hourly_discharge, "MC", "2015-07-06-00-00", "2015-07-28-00-00", 24)
events <- add.event(WQ_hourly_discharge, "MC", "2015-08-18-00-00", "2015-09-04-00-00", 24)
events <- add.event(WQ_hourly_discharge, "MC", "2015-09-17-00-00", "2015-09-24-00-00", 24)
events <- add.event(WQ_hourly_discharge, "MC", "2015-09-09-00-00", "2015-09-17-00-00", 24)
events <- add.event(WQ_hourly_discharge, "MC", "2017-08-16-00-00", "2017-09-17-00-00", 24)
events <- add.event(WQ_hourly_discharge, "MC", "2018-09-03-00-00", "2018-09-15-00-00", 12)

# get rid of the duplicate rows
# (for some reason, some of the events were added multiple times)
eventsUnique <- distinct(events)

# add precip and antecedent discharge data 
eventsUnique2 <- eventsUnique %>% 
  mutate(precip = mapply(totalPrecip, startDate, endDate, MoreArgs = list(precip2)),
         precipYear = mapply(yearTotalPrecip, startDate, MoreArgs = list(precip2)),
         prev_month_dis = mapply(prev.month.discharge, startDate, MoreArgs = list(unitDischarge,900)))

# write data to csv
write_csv(eventsUnique2, "data/events.csv")
```

```{r}
eventsUnique2 <- read_csv("data/events.csv")

```

Finally, we reshape the data so that there are 4 rows for each event (one for each constituent). This makes plotting easier, so that each constituent doesn't have to be plotted separately. We also add columns for the direction of hysteresis (clockwise or counterclockwise) and slope (flushing or dilution). In Aguilera and Melack (2018), they classify clockwise hysteresis as HI >= 0.05, counterclockwise hysteresis as HI <= -0.05, flushing as $\Delta C$ >= 0.1, and dilution as $\Delta C$ <= 0.1. For now, I'm going to use 0.05 as the limit for flushing and dilution.

```{r reshaping the data}
eventsReshaped <- eventsUnique2 %>% 
  gather(key, value, -c(startDate, endDate, length, length_loop,
                        dis_change, dis_change_loop, max_dis, avg_temp, 
                        precip, precipYear, prev_month_dis)) %>% 
  separate(key, c("measure", "var"), "_") %>% 
  tidyr::spread(measure, value)

# add columns for directions of hysteresis, slope, and complexity
eventsFinal <- eventsReshaped %>% 
  mutate(H_dir = case_when(HI >= 0.05 ~ "clockwise",
                           HI <= -0.05 ~ "counter-clockwise", T ~ "none"),
         slope_dir = case_when(Slope >= 0.05 ~ "flushing",
                               Slope <= -0.05 ~ "dilution", T ~ "constant"),
         complex = case_when(area >= abs(2*HI) ~ "complex",
                             T ~ "not complex"))
# one event (starting 9/9/15) had a value >1 for fdom slope because of negative
# values at the beginning-- might be best to take out the whole fdom row for
# this event because the data jumps from negatives to 100 (the issue we've been 
# looking at for the time series)
eventsFinal2 <- eventsFinal[-c(70,87,88),]

```

# Data exploration and visualization

Now that we have a dataframe of all the events, we can visualize the data in various ways.

## Event summary
First, we can summarize the event metrics. Here is a heatmap showing when the events occurred over time (some data wrangling is required first to get a full dataframe of event counts):
```{r}
# use subset with just turbidity so events aren't repeated
eventsTurb <- eventsFinal2 %>% filter(var == "turb")

# discrete month column
eventsTurb2 <- eventsTurb %>% 
  mutate(month = case_when(month(startDate) == 4 ~ "April",
                           month(startDate) == 6 ~ "June",
                           month(startDate) == 7 ~ "July",
                           month(startDate) == 8 ~ "August",
                           month(startDate) == 9 ~ "September"))
eventsTurb2$month <- factor(eventsTurb2$month,
                            levels = c("April", "June", "July", "August", "September"))

# get number of events starting in each month of each year
eventsCount <- (eventsTurb2 %>% 
  mutate(year = as.factor(year(startDate))))[,c(1, 19, 20)] %>% 
  group_by(year, month) %>% tally() 

# make a dataframe of all possible combinations of month and year
allMonthYear <- eventsCount %>% expand(month, year) 
allMonthYear <- allMonthYear %>% as_tibble() %>% 
  add_row(month = "May", year = as.factor(2015)) %>% 
   add_row(month = "May", year = as.factor(2016)) %>% 
   add_row(month = "May", year = as.factor(2017)) %>% 
   add_row(month = "May", year = as.factor(2018))

# make full events count dataframe
eventsCount2 <- full_join(eventsCount, allMonthYear) %>% 
  mutate(n = case_when(!is.na(n) ~ as.character(n),
                       T ~ "0"))
eventsCount2$month <- factor(eventsCount2$month,
                            levels = c("April", "May", "June", "July", "August", "September"))

# figure: heatmap
png(filename = "Figures/monthly_event_count_grid.png", width = 6, height = 5.5,
    units = "in", res = 300, bg = "transparent")
eventsCount2 %>% 
  ggplot(aes(x = year, y = month, fill = n)) +
  geom_tile(color = "black") +
  scale_fill_manual(values = c("white", brewer.pal(3, "Paired")[1], brewer.pal(3, "Paired")[2]),
                    name = "Number of Events") +
  theme_bw() + xlab("Year") + ylab("Month") +
  theme(plot.background = element_rect(fill = "transparent",colour = NA),
        legend.background = element_rect(fill = "transparent"), legend.position = "bottom") 
dev.off()  
```

We can also make histograms to see the distributions of the event metrics:
```{r}
png(filename = "Figures/event_hists.png", width = 5, height = 5,
    units = "in", res = 300, bg = "transparent")
par(mfrow = c(3,3), bg = "transparent", mgp=c(2.5,1,0), mar = c(5,2,1,2))
hist(eventsTurb$avg_temp, breaks = c(14, 17, 20, 23, 26, 29),
     xlab = "Average temperature (C)", xlim = c(10, 30),
     ylab = "Number of events", main = "")
hist(eventsTurb$length, breaks = 5,
     xlab = "Length (days)", 
     ylab = "Number of events", main = "")
hist(eventsTurb$dis_change, breaks = c(0, 16000, 32000, 48000, 64000, 80000),
     xlab = "Change in discharge (ft^3/s)", 
     ylab = "Number of events", main = "")
hist(eventsTurb$max_dis, breaks = 5,
     xlab = "Peak discharge (ft^3/s)", 
     ylab = "Number of events", main = "")
hist(eventsTurb$precip, breaks = 5,
     xlab = "Cumulative precipitation (in)", 
     ylab = "Number of events", main = "")
hist(eventsTurb$precipYear, breaks = 5,
     xlab = "Precipitation since Jan 1 (in)", 
     ylab = "Number of events", main = "")
plot.new()
hist(eventsTurb$prev_month_dis, breaks = 5,
     xlab = "Previous month discharge", 
     ylab = "Number of events", main = "")
dev.off()
```

```{r proportions of events}
# get the total number of events classified as each slope/hi combo for MC
eventsRawFinal2 %>% filter(var == "bga") %>% group_by(H_dir, slope_dir) %>% count() 
eventsRawFinal2 %>% filter(var == "chl") %>% group_by(H_dir, slope_dir) %>% count()
eventsRawFinal2 %>% filter(var == "fdom") %>% group_by(H_dir, slope_dir) %>% count()
eventsRawFinal2 %>% filter(var == "n03") %>% group_by(H_dir, slope_dir) %>% count()
eventsRawFinal2 %>% filter(var == "turb") %>% group_by(H_dir, slope_dir) %>% count()

# get the total number of events classified as each slope/hi combo for SI
events_SI_raw_Final2 %>% filter(var == "bga") %>% group_by(H_dir, slope_dir) %>% count() #11
events_SI_raw_Final2 %>% filter(var == "chl") %>% group_by(H_dir, slope_dir) %>% count() #11
events_SI_raw_Final2 %>% filter(var == "fdom") %>% group_by(H_dir, slope_dir) %>% count() #11
events_SI_raw_Final2 %>% filter(var == "n03") %>% group_by(H_dir, slope_dir) %>% count() #11
events_SI_raw_Final2 %>% filter(var == "turb") %>% group_by(H_dir, slope_dir) %>% count() #11

# add a column for total number of events calculated for the constituent
eventsRawFinal2 <- eventsRawFinal2 %>% 
  mutate(varCount = case_when(var == "bga" ~ 16,
                              var == "chl" ~ 18,
                              var == "fdom" ~ 17,
                              var == "n03" ~ 18,
                              var == "turb" ~ 18))
```


## HI/slope summary
Next, we can look at the distributions of hysteresis and slope between constituents and events. 

First, we can make scatterplots with hysteresis index along the x-axis and slope along the y-axis. We can plot the data in this way both event-wise (a plot for each event with a dot for each constituent), or constituent-wise (a plot for each constituent with a dot for each event). 

```{r HI/slope scatterplots}
# New facet label names for constituents
var.labs <- c("Cyanobacteria", "Chlorophyll", "FDOM", "Nitrate", "Turbidity")
names(var.labs) <- c("bga", "chl", "fdom", "n03", "turb")

# event-wise plots
eventsFinal2 %>% 
  ggplot(aes(x = HI, y = Slope, col = var)) +
  geom_hline(yintercept=0, linetype="dashed", color = "black") +
  geom_vline(xintercept=0, linetype="dashed", color = "black") +
  geom_point() + xlim(-1, 1) + ylim(-1, 1) +
  facet_wrap(~as.factor(startDate)) +
  scale_color_brewer(palette = "Dark2", name = "Constituent", labels = var.labs) 

# constituent-wise plots
png(filename = "Figures/HI_slope_point.png", 
    width = 4, height = 5, units = "in", res = 300, bg = "transparent")
eventsFinal2 %>% 
  ggplot(aes(x = HI, y = Slope)) +
  geom_hline(yintercept=0, linetype="dashed", color = "red") +
  geom_vline(xintercept=0, linetype="dashed", color = "red") +
  geom_point() + xlim(-1, 1) + ylim(-1, 1) +
  facet_wrap(~var, labeller = labeller(var = var.labs), nrow = 3) +
  theme_bw() + xlab("Hysteresis Index") +
  theme(plot.background = element_rect(fill = "transparent",colour = NA))
dev.off()

png(filename = "Figures/HI_slope_point_v2.png", 
    width = 4, height = 5, units = "in", res = 300, bg = "transparent")
eventsV2Final2 %>% 
  ggplot(aes(x = HI, y = Slope)) +
  geom_hline(yintercept=0, linetype="dashed", color = "red") +
  geom_vline(xintercept=0, linetype="dashed", color = "red") +
  geom_point() + xlim(-1, 1) + ylim(-1, 1) +
  facet_wrap(~var, labeller = labeller(var = var.labs), nrow = 3) +
  theme_bw() + xlab("Hysteresis Index") +
  theme(plot.background = element_rect(fill = "transparent",colour = NA))
dev.off()

png(filename = "Figures/HI_slope_point_raw.png", 
    width = 4, height = 5, units = "in", res = 300, bg = "transparent")
eventsRawFinal2 %>% 
  ggplot(aes(x = HI, y = Slope)) +
  geom_hline(yintercept=0, linetype="dashed", color = "red") +
  geom_vline(xintercept=0, linetype="dashed", color = "red") +
  geom_point() + xlim(-1, 1) + ylim(-1, 1) +
  facet_wrap(~var, labeller = labeller(var = var.labs), nrow = 3) +
  theme_bw() + xlab("Hysteresis Index") +
  theme(plot.background = element_rect(fill = "transparent",colour = NA))
dev.off()


```

Next, we can visualize the categorical HI and slope data by making barplots
```{r HI/slope barplots}
# need to reorder the slope directions
eventsFinal2$slope_dir <- factor(eventsFinal2$slope_dir, 
                                 levels = c("flushing", "dilution", "constant"))

# dealing with complex hysteresis status
events_hyst_type <- eventsFinal2 %>% 
  mutate(hyst_type = case_when(complex == "complex" ~ "complex",
                               H_dir == "clockwise" ~ "clockwise",
                               H_dir == "counter-clockwise" ~ "counter-clockwise",
                               T ~ "no hysteresis"))
events_hyst_type$hyst_type <- factor(events_hyst_type$hyst_type, 
                                     levels = c("clockwise", "counter-clockwise", 
                                                "complex", "no hysteresis"))
# bar plot of hysteresis type
hiBar <- events_hyst_type %>% 
  ggplot(aes(x = var, fill = hyst_type)) +
  geom_bar(position = position_dodge(preserve = "single")) +
  scale_fill_brewer(palette = "Paired", name = "Hysteresis Type") +
  ylab("Number of Events") + xlab("Constituent") +
  scale_x_discrete(labels = var.labs) +
  scale_y_continuous(breaks = c(2, 4, 6, 8, 10, 12)) +
  theme_bw() + 
  theme(plot.background = element_rect(fill = "transparent", colour = NA),
        legend.background = element_rect(fill = "transparent")) 

# bar plot of slope type
slopeBar <- eventsFinal2 %>% 
  ggplot(aes(x = var, fill = slope_dir)) +
  geom_bar(position = "dodge") +
  scale_fill_brewer(palette = "Paired", name = "Slope Type") +
  ylab("Number of Events") + xlab("Constituent") +
  scale_x_discrete(labels = var.labs) +
  scale_y_continuous(breaks = c(2, 4, 6, 8, 10, 12)) +
  theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA),
                     legend.background = element_rect(fill = "transparent")) 

# figure with both bar plots
#png(filename = "HI_slope_bar.png", 
    #width = 6, height = 6, units = "in", res = 300, bg = "transparent")
ggarrange(hiBar, slopeBar, nrow = 2)
#dev.off()

# need to reorder the slope directions
eventsFinal2$slope_dir <- factor(eventsFinal2$slope_dir, 
                                 levels = c("flushing", "constant", "dilution"),
                                 labels = c("Flushing", "Constant", "Dilution"))
eventsV2Final2$slope_dir <- factor(eventsV2Final2$slope_dir, 
                                 levels = c("flushing", "constant", "dilution"),
                                 labels = c("Flushing", "Constant", "Dilution"))
eventsRawFinal2$slope_dir <- factor(eventsRawFinal2$slope_dir, 
                                 levels = c("flushing", "constant", "dilution"),
                                 labels = c("Flushing", "Constant", "Dilution"))
# need to reorder the hysteresis directions
eventsFinal2$H_dir <- factor(eventsFinal2$H_dir, 
                             levels = c("clockwise", "none", "counter-clockwise"),
                             labels = c("Clockwise", "No Hysteresis", "Counter-clockwise"))
eventsV2Final2$H_dir <- factor(eventsV2Final2$H_dir, 
                             levels = c("clockwise", "none", "counter-clockwise"),
                             labels = c("Clockwise", "No Hysteresis", "Counter-clockwise"))
eventsRawFinal2$H_dir <- factor(eventsRawFinal2$H_dir, 
                             levels = c("clockwise", "none", "counter-clockwise"),
                             labels = c("Clockwise", "No Hysteresis", "Counter-clockwise"))

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

png("bar_grid_MC_v2.png", width = 6, height = 5,
    units = "in", res = 300, bg = "transparent")
eventsV2Final2 %>% 
  ggplot(aes(x = var, y = (..count..)/nrow(eventsTurb), fill = var)) +
  geom_bar() +
  facet_grid(cols = vars(H_dir), rows = vars(slope_dir), drop = F) +
  theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA),
                     legend.background = element_rect(fill = "transparent"),
                     axis.text.x = element_blank(), axis.ticks.x = element_blank(),
                     legend.position = "bottom") +
  scale_fill_manual(values = cbPalette, labels = var.labs, name = "Constituent") +
  xlab("") + ylab("Proportion of events")
dev.off()

png("bar_grid_MC_raw.png", width = 6, height = 5,
    units = "in", res = 300, bg = "transparent")
eventsRawFinal2 %>% 
  ggplot(aes(x = var, y = (..count..)/18, fill = var)) +
  geom_bar() +
  facet_grid(cols = vars(H_dir), rows = vars(slope_dir), drop = F) +
  theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA),
                     legend.background = element_rect(fill = "transparent"),
                     axis.text.x = element_blank(), axis.ticks.x = element_blank(),
                     legend.position = "bottom") +
  scale_fill_manual(values = cbPalette, labels = var.labs, name = "Constituent") +
  xlab("") + ylab("Proportion of events") +
  ylim(c(0,0.5))
dev.off()

eventsRawFinal2 %>% 
  ggplot(aes(x = var, y = (..count..), fill = var)) +
  geom_bar() +
  facet_grid(cols = vars(H_dir), rows = vars(slope_dir), drop = F) +
  theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA),
                     legend.background = element_rect(fill = "transparent"),
                     axis.text.x = element_blank(), axis.ticks.x = element_blank(),
                     legend.position = "bottom") +
  scale_fill_manual(values = cbPalette, labels = var.labs, name = "Constituent") +
  xlab("") + ylab("Proportion of events") +
  ylim(c(0,0.6))

eventsRawFinal2 %>% 
  group_by(var, H_dir, slope_dir) %>% 
  count()
```

```{r classification without none/constant}
eventsRawFinal3 <- eventsRawFinal2 %>% 
  mutate(H_dir2 = case_when(HI > 0 ~ "Clockwise",
                            T ~ "Counter-clockwise"),
         slope_dir2 = case_when(Slope > 0 ~ "Flushing",
                                T ~ "Dilution"))

png("bar_grid_MC_4.png", width = 6, height = 5,
    units = "in", res = 300, bg = "transparent")
eventsRawFinal3 %>% 
  ggplot(aes(x = var, y = (..count..)/18, fill = var)) +
  geom_bar() +
  facet_grid(cols = vars(H_dir2), rows = vars(slope_dir2), drop = F) +
  theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA),
                     legend.background = element_rect(fill = "transparent"),
                     axis.text.x = element_blank(), axis.ticks.x = element_blank(),
                     legend.position = "bottom",
                     panel.background = element_rect(fill = "transparent")) +
  scale_fill_manual(values = cbPalette, labels = var.labs, name = "Constituent") +
  xlab("") + ylab("Proportion of events") +
  ylim(c(0,0.75))
dev.off()

png("bar_grid_SI_4.png", width = 6, height = 5,
    units = "in", res = 300, bg = "transparent")
events_SI_raw_Final3 %>% 
  ggplot(aes(x = var, y = (..count..)/11, fill = var)) +
  geom_bar() +
  facet_grid(cols = vars(H_dir2), rows = vars(slope_dir2), drop = F) +
  theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA),
                     legend.background = element_rect(fill = "transparent"),
                     axis.text.x = element_blank(), axis.ticks.x = element_blank(),
                     legend.position = "bottom",
                     panel.background = element_rect(fill = "transparent")) +
  scale_fill_manual(values = cbPalette, labels = var.labs, name = "Constituent") +
  xlab("") + ylab("Proportion of events") +
  ylim(c(0,0.75))
dev.off()


```


## Pairwise constituent comparisons
```{r}

```


# Relationships between event metrics and constituent responses

```{r single regressions function}
eventsFinalMetric <- eventsFinal2 %>% 
  mutate(dis_change = 0.02831685*dis_change,
         max_dis = 0.02831685*max_dis,
         precip = 2.54*precip,
         precipYear = 2.54*precipYear,
         prev_month_dis = 0.02831685*prev_month_dis)

eventsV2FinalMetric <- eventsV2Final2 %>% 
   mutate(dis_change = 0.02831685*dis_change,
         max_dis = 0.02831685*max_dis,
         precip = 2.54*precip,
         precipYear = 2.54*precipYear,
         prev_month_dis = 0.02831685*prev_month_dis)

eventsRawFinalMetric <- eventsRawFinal2 %>% 
   mutate(dis_change = 0.02831685*dis_change,
         max_dis = 0.02831685*max_dis,
         precip = 2.54*precip,
         precipYear = 2.54*precipYear,
         prev_month_dis = 0.02831685*prev_month_dis)

# function that returns the p value from a linear model output
lmp <- function (modelobject) {
    if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
    f <- summary(modelobject)$fstatistic
    p <- pf(f[1],f[2],f[3],lower.tail=F)
    attributes(p) <- NULL
    return(p)
}

# function to add a row of regression info to a tibble
add.regression.MC <- function(data, cons, met){
  subset <- data %>% filter(var == cons)
  metVec <- subset %>% pull(met)
  regHI <- lm(subset$HI ~ metVec)
  regSl <- lm(subset$Slope ~ metVec)
  
  return(regressions_MC %>% add_row(constituent = cons, metric = met,
                                    p_HI = lmp(regHI), 
                                    Rsqd_HI = summary(regHI)$r.squared,
                                    slope_HI = as.numeric(regHI$coefficients[2]),
                                    intercept_HI = as.numeric(regHI$coefficients[1]),
                                    p_Sl = lmp(regSl), 
                                    Rsqd_Sl = summary(regSl)$r.squared,
                                    slope_Sl = as.numeric(regSl$coefficients[2]),
                                    intercept_Sl = as.numeric(regSl$coefficients[1])))
  
}


# initialize table of regression info
regressions_MC <- tibble (constituent = character(), metric = character(),
                          p_HI = numeric(), Rsqd_HI = numeric(), 
                          slope_HI = numeric(), intercept_HI = numeric(),
                          p_Sl = numeric(), Rsqd_Sl = numeric(),
                          slope_Sl = numeric(), intercept_Sl = numeric())


# add all regressions
constituents <- c("bga", "chl", "fdom", "n03", "turb")
metrics <- c("length", "dis_change", "max_dis", "avg_temp", 
             "precip", "precipYear", "prev_month_dis")

for(c in constituents){
  for(m in metrics){
    regressions_MC <- add.regression.MC(eventsFinalMetric, c, m)
  }
}

write_csv(regressions_MC, "data/reg_MC_metric.csv")

##### v2 #####

# function to add a row of regression info to a tibble
add.regression.MC.V2 <- function(data, cons, met){
  subset <- data %>% filter(var == cons)
  metVec <- subset %>% pull(met)
  regHI <- lm(subset$HI ~ metVec)
  regSl <- lm(subset$Slope ~ metVec)
  
  return(regressions_MC_V2 %>% add_row(constituent = cons, metric = met,
                                    p_HI = lmp(regHI), 
                                    Rsqd_HI = summary(regHI)$r.squared,
                                    slope_HI = as.numeric(regHI$coefficients[2]),
                                    intercept_HI = as.numeric(regHI$coefficients[1]),
                                    p_Sl = lmp(regSl), 
                                    Rsqd_Sl = summary(regSl)$r.squared,
                                    slope_Sl = as.numeric(regSl$coefficients[2]),
                                    intercept_Sl = as.numeric(regSl$coefficients[1])))
  
}


# initialize table of regression info
regressions_MC_V2 <- tibble (constituent = character(), metric = character(),
                          p_HI = numeric(), Rsqd_HI = numeric(), 
                          slope_HI = numeric(), intercept_HI = numeric(),
                          p_Sl = numeric(), Rsqd_Sl = numeric(),
                          slope_Sl = numeric(), intercept_Sl = numeric())


# add all regressions
constituents <- c("bga", "chl", "fdom", "n03", "turb")
metrics <- c("length", "dis_change", "max_dis", "avg_temp", 
             "precip", "precipYear", "prev_month_dis")

for(c in constituents){
  for(m in metrics){
    regressions_MC_V2 <- add.regression.MC.V2(eventsV2FinalMetric, c, m)
  }
}

write_csv(regressions_MC_V2, "data/reg_MC_metric_v2.csv")


##### raw data version #####

# function to add a row of regression info to a tibble
add.regression.MC.raw <- function(data, cons, met){
  subset <- data %>% filter(var == cons)
  metVec <- subset %>% pull(met)
  regHI <- lm(subset$HI ~ metVec)
  regSl <- lm(subset$Slope ~ metVec)
  
  return(regressions_MC_raw %>% add_row(constituent = cons, metric = met,
                                    p_HI = lmp(regHI), 
                                    Rsqd_HI = summary(regHI)$r.squared,
                                    slope_HI = as.numeric(regHI$coefficients[2]),
                                    intercept_HI = as.numeric(regHI$coefficients[1]),
                                    p_Sl = lmp(regSl), 
                                    Rsqd_Sl = summary(regSl)$r.squared,
                                    slope_Sl = as.numeric(regSl$coefficients[2]),
                                    intercept_Sl = as.numeric(regSl$coefficients[1])))
  
}


# initialize table of regression info
regressions_MC_raw <- tibble (constituent = character(), metric = character(),
                          p_HI = numeric(), Rsqd_HI = numeric(), 
                          slope_HI = numeric(), intercept_HI = numeric(),
                          p_Sl = numeric(), Rsqd_Sl = numeric(),
                          slope_Sl = numeric(), intercept_Sl = numeric())


# add all regressions
constituents <- c("bga", "chl", "fdom", "n03", "turb")
metrics <- c("length", "dis_change", "max_dis", "avg_temp", 
             "precip", "precipYear", "prev_month_dis")

for(c in constituents){
  for(m in metrics){
    regressions_MC_raw <- add.regression.MC.raw(eventsRawFinalMetric, c, m)
  }
}

write_csv(regressions_MC_raw, "data/reg_MC_raw.csv")
```

To explore the relationships between event metrics and HI/slope responses for each constituent, we can create multiple linear regressions using the event metrics as predictor variables and HI and slope as response variables. Here are functions to do this:
```{r}
# multiple linear regression function for HI
reg_multi_HI <- function(variable){
  subset <- eventsFinal2 %>% 
    filter(var == variable)
  
  return((lm(subset$HI ~ 
             subset$dis_change + subset$length + subset$max_dis + subset$prev_month_dis +
             subset$avg_temp + subset$precip + subset$precipYear )))
}

# multiple linear regression function for slope
reg_multi_slope <- function(variable){
  subset <- eventsFinal2 %>% 
    filter(var == variable)
  
  return((lm(subset$Slope ~ 
             subset$dis_change + subset$length + subset$max_dis + subset$prev_month_dis +
             subset$avg_temp + subset$precip + subset$precipYear )))
}
```

To determine the best combination of predictor variables to model HI and slope for each constituent, we can conduct step AIC analyses:
```{r step AIC}
# ================ Turbidity ================
step(reg_multi_HI("turb"))
summary(lm(eventsTurb$HI ~ eventsTurb$length + eventsTurb$max_dis + 
             eventsTurb$precipYear + eventsTurb$precip))
step(reg_multi_slope("turb"))
summary(lm(eventsTurb$Slope ~ eventsTurb$prev_month_dis + 
             eventsTurb$avg_temp + eventsTurb$precip))

# ================ Nitrate =================
step(reg_multi_HI("n03"))
summary(lm(eventsN03$HI ~ eventsN03$avg_temp))
step(reg_multi_slope("n03"))
summary(lm(eventsN03$Slope ~ eventsN03$dis_change + eventsN03$precip +
             eventsN03$length + eventsN03$prev_month_dis))

# ================ FDOM =================
step(reg_multi_HI("fdom"))
summary(lm(eventsFDOM$HI ~ eventsFDOM$length))
step(reg_multi_slope("fdom"))
summary(lm(eventsFDOM$Slope ~ eventsFDOM$length))


# ================ Chlorophyll ==================
step(reg_multi_HI("chl"))
summary(lm(eventsChl$HI ~ 1))
step(reg_multi_slope("chl"))
summary(lm(eventsChl$Slope ~ 1))

# ================ BGA =================
step(reg_multi_HI("bga"))
summary(lm(eventsBGA$HI ~ 1))
step(reg_multi_slope("bga"))
summary(lm(eventsBGA$Slope ~ eventsBGA$prev_month_dis))
```

For the linear models that only include one predictor variable, we can plot the linear model:
```{r linear model plots}
#png(filename = "nO3_HI_temp.png", width = 4, height = 3,
   # units = "in", res = 300, bg = "transparent")
eventsN03 %>% ggplot(aes(x = avg_temp, y = HI)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
  geom_point() + geom_smooth(method = "lm", se = F) +
  xlab("Average Temperature") + ylab("Hysteresis Index") +
  theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA))
#dev.off()

#png(filename = "fdom_HI_length.png", width = 4, height = 3,
    #units = "in", res = 300, bg = "transparent")
eventsFDOM %>% ggplot(aes(x = length, y = HI)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
  geom_point() + geom_smooth(method = "lm", se = F) +
  xlab("Event length (days)") + ylab("Hysteresis Index") +
  theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA))
#dev.off()

#png(filename = "fdom_S_length.png", width = 4, height = 3,
    #units = "in", res = 300, bg = "transparent")
eventsFDOM %>% ggplot(aes(x = length, y = Slope)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
  geom_point() + geom_smooth(method = "lm", se = F) +
  xlab("Event length (days)") + ylab("Slope") +
  theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA))
#dev.off()
```

Since the ratio of the number of predictor variables to the sample size is too high, step AIC/ multiple regression may not be valid. Instead, we will report the statistics for each individual single-variable regression.
```{r single regressions}

```


To calculate the event metrics using a time lag, we will create a new WQ_hourly_discharge dataframe which uses the appropriate lag for each year
```{r time lag}
# average velocity for each year
yearly_velocity <- pool8_vel %>% 
  filter(month %in% c(4:10)) %>% 
  group_by(year) %>% 
  summarise(mean_vel = mean(VEL))

# function to calculate time lag for each year
year_lag <- function(Year){
  # get velocity from the relevant year
  vel <- yearly_velocity$mean_vel[which(yearly_velocity$year == Year)]
  # calculate time lag (43452.3 meters)
  seconds <- 43452.3/vel
  hours <- seconds/3600
  if(hours - floor(hours) < 0.5){
    return(floor(hours))
  }
  return(floor(hours) + 1)
}

MC_vel <- read_csv("data/Pool8_VEL_MC.csv")

MC_vel <- MC_vel %>% 
  mutate(date = mdy(DATE)) %>% 
  filter(year(date) > 2014, year(date) < 2019)

library(birk)

# function to calculate time lag for each year
event_lag <- function(startDate, endDate){
  vel <- MC_vel %>% drop_na()
  # get closest date in velocity df to start date
  closestStartDate <- vel$date[which.closest(vel$date, startDate)]
  # calculate average velocity 
  velStart <- vel %>% 
    filter(date == closestStartDate) %>% 
    summarize(meanVel = mean(VEL))
  meanVelStart <- velStart$meanVel
  # get closest date in velocity df to end date
  closestEndDate <- vel$date[which.closest(vel$date, endDate)]
  # calculate average velocity
  velEnd <- vel %>% 
    filter(date == closestEndDate) %>% 
    summarize(meanVel = mean(VEL))
  meanVelEnd <- velEnd$meanVel
  # calculate average of start and end date velocities
  meanVel <- mean(meanVelStart, meanVelEnd)
  # calculate time lag (43452.3 meters)
  seconds <- 43452.3/meanVel
  hours <- seconds/3600
  # return time lag rounded to the nearest hour
  if(hours - floor(hours) < 0.5){
    return(floor(hours))
  }
  return(floor(hours) + 1)
}

eventsLag <- eventsFinal2 %>% 
  mutate(lag = mapply(event_lag, date(startDate), date(endDate)))

```