This report is automatically generated with the R package knitr (version 1.39) .

---
title: "Tidal Wetlands"
author: "Mercury Program and Basin Planning Unit"
date: "2/23/2022"
output:
  html_document:
    code_folding: show
    toc: TRUE
    toc_float: TRUE
    toc_depth: 3
runtime: shiny
assets:
  css:
    - "http://fonts.googleapis.com/css?family=Raleway:300"
    - "http://fonts.googleapis.com/css?family=Oxygen"
---

<style>
body{
  font-family: 'Oxygen', sans-serif;
  font-size: 16px;
  line-height: 24px;
}

h1,h2,h3,h4 {
  font-family: 'Raleway', sans-serif;
}

.container { width: 1250px; }
h3 {
  background-color: #D4DAEC;
  text-indent: 50px; 
}
h4 {
  text-indent: 75px;
  margin-top: 35px;
  margin-bottom: 5px;
}
</style>

```{r setup, echo=TRUE, warning=FALSE, message=FALSE, fig.width=9.5}
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE, fig.width=9.5)
```


```{r Libraries, echo=FALSE}
# library(here)
# library(lubridate)
# library(readxl)
library(janitor) # for clean_names()
# library(openxlsx) # for convertToDate
# library(shiny)
# library(naniar)
# library(survival) # for neardate()
# library(ggplot2)
# library(plotly)
# library(hrbrthemes)
# library(Hmisc) # for %nin%

# Had issue trying to set WD with Shiny in R project, reset working directory of rproj
wd <- rstudioapi::getActiveProject()
source(paste0(wd, "/R Functions/functions_estimate NDDNQ values.R"))
source(paste0(wd, "/R Functions/functions_QA data.R"))
```

# Load and Tidy Data

Define some variables for repeat use later
```{r}
columns_mehg_select <- c("SampleDateTime", "StationName", "SampleCode", "SourceID", "MDL", "RL", "UnitMeHg", "SiteType", "Result", "ResultQualCode", "LabComments", "QAComments", "Month")

columns_flow_select <- c("SampleDateTime", "StationName", "SourceID", "Flow", "UnitFlow", "SiteType", "Comments", "Month")
```


## All Wetlands MeHg
```{r}
all_mehg <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/DWR Tidal Wetlands_Aq MeHg and Flow Data.xlsx"), 
                       col_types = c("text", "text", "text", "date", "guess", "text", "numeric", "guess", "numeric","numeric", "text", "text", "text", "guess"),
                       sheet = 2)
```

Common cleaning code
```{r}
all_mehg_clean <- all_mehg %>% 
  clean_names() %>%
  filter(analyte == "MeHg- total") %>% 
  mutate(SampleTime = format(collection_time_pst, "%H:%M:%S")
    ) %>% 
  mutate(SampleDateTime = as.POSIXct(paste(sample_date, SampleTime)), format="%d/%m/%Y %H.%M.%S") %>% 
  dplyr::rename(
    StationName = station_name,
    SampleCode = sample_code,
    LabComments = lab_comments,
    QAComments = mme_comments,
    Result = result,
    RL = rl,
    MDL = mdl,
    UnitMeHg = units
    ) %>%
  # Add missing columns
  add_column(
    SourceID = "DWR Tidal Wetlands",
    ResultQualCode = NA_character_,
    ) %>% 
  mutate(Month = month(SampleDateTime)) 
# %>% 
  # mutate(SiteType = NA_character_) %>% 
  #          case_when(grepl("flood", SiteType, ignore.case = TRUE) 
  #                    ~ "Flood",
  #                    grepl("ebb", SiteType, ignore.case = TRUE)
  #                    ~ "Ebb")) %>% 
  # dplyr::select(all_of(columns_mehg_select))
```


## Westervelt
Load flow data
```{r}
west_flow <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/DWR Tidal Wetlands_Aq MeHg and Flow Data.xlsx"),
                    sheet = 3)
```

Clean flow data
```{r}
west_flow_clean <- west_flow %>% 
  clean_names() %>%
  dplyr::rename(
    SampleDateTime = date_and_time,
    Comments = comments,
    SiteType = tidal_cycle) %>% 
  add_column(
    StationName = "Westervelt", 
    UnitFlow = "cfs",
    SourceID = "DWR Tidal Wetlands") %>% 
  mutate(Month = month(SampleDateTime),
         Flow = case_when(!is.na(raw_flow) ~ raw_flow,
                          TRUE ~ interpolated_flow)) %>% 
  dplyr::select(all_of(columns_flow_select)) %>% 
  mutate(SiteType = case_when(SiteType == "Slack" ~ NA_character_,
                              TRUE ~ SiteType)) %>% 
  arrange(SampleDateTime) %>%
  fill(SiteType, .direction = "down") %>% 
  filter(!is.na(Flow))
```

Select methylmercury data & join with SiteType from flow
```{r}
west_mehg_final <- all_mehg_clean %>% 
  filter(StationName == "Westervelt") %>% 
  left_join(west_flow_clean) %>% 
  dplyr::select(all_of(columns_mehg_select))
```


## North Lindsey
Select methylmercury data
```{r}
linds_mehg <- all_mehg_clean %>% 
  filter(StationName == "N Lindsey Slough Tidal Wetland")
```

Load flow data
```{r}
linds_flow <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/DWR Tidal Wetlands_Aq MeHg and Flow Data.xlsx"),
                    sheet = 4)
```

Clean flow data
```{r}
linds_flow_clean <- linds_flow %>% 
  clean_names() %>%
  # Filter only for good data with Qualities 1, 2, and 70 (based on key in RAW data column D)
  filter(qual == 1 | qual == 2 | qual == 70) %>% 
  mutate(Comments = paste("Quality Code ", qual),
         SiteType = case_when(point < 0 ~ "Flood",
                              point > 0 ~ "Ebb",
                              point == 0 ~ NA_character_)) %>%
  dplyr::rename(
    Flow = point,
    SampleDateTime = date) %>% 
  add_column(
    StationName = "North Lindsey", 
    UnitFlow = "cfs",
    SourceID = "DWR Tidal Wetlands") %>% 
  mutate(Month = month(SampleDateTime)) %>% 
  dplyr::select(all_of(columns_flow_select)) %>%
  arrange(SampleDateTime) %>%
  fill(SiteType, .direction = "down")
```

```{r}
# Join both tables, use fill to append site type based on closest SampleDateTime
linds_combined <- linds_flow_clean %>%
  full_join(linds_mehg, by = c("SampleDateTime", "StationName", "SourceID", "Month")) %>%
  mutate(SiteTypeFilled = SiteType) %>% 
  arrange(SampleDateTime) %>% 
  fill(SiteTypeFilled, .direction = "up") %>% 
  select(-SiteType) %>% 
  rename(SiteType = SiteTypeFilled)

# Manually make sure that SiteType matches that of the next closest reported time
linds_check <- linds_combined %>% 
  dplyr::select(SampleDateTime, Flow, Result, SiteType, SiteType)
# It does in every scenario

# Join site types back to final tables
linds_mehg_final <- linds_mehg %>% 
  left_join(linds_combined) %>% 
  dplyr::select(all_of(columns_mehg_select))

linds_flow_clean2 <- linds_flow_clean %>% 
  dplyr::select(-SiteType) %>% 
  left_join(linds_combined) %>% 
  dplyr::select(all_of(columns_flow_select))
```


## Yolo Bypass
Select methylmercury data
```{r}
yolo_mehg <- all_mehg_clean %>% 
  filter(StationName == "YWA Tidal Wetland")
```

Load flow data
```{r}
yolo_flow <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/DWR Tidal Wetlands_Aq MeHg and Flow Data.xlsx"),
                    sheet = 5)
```

Clean flow data
```{r}
yolo_flow_clean <- yolo_flow %>% 
  clean_names() %>%
  # Filter only for good data with Qualities 1, 2, and 70 (based on key in RAW data column D)
  filter(qual == 1 | qual == 2 | qual == 70) %>%
  mutate(Comments = paste("Quality Code ", qual),
         SiteType = case_when(point < 0 ~ "Flood",
                              point > 0 ~ "Ebb",
                              point == 0 ~ NA_character_)) %>% 
  dplyr::rename(
    Flow = point,
    SampleDateTime = date) %>% 
  add_column(
    StationName = "YWA Tidal Wetland", 
    UnitFlow = "cfs",
    SourceID = "DWR Tidal Wetlands") %>% 
  mutate(Month = month(SampleDateTime)) %>% 
  dplyr::select(all_of(columns_flow_select)) %>%
  arrange(SampleDateTime) %>%
  fill(SiteType, .direction = "down")
```

```{r}
# Join both tables, use fill to append site type based on closest SampleDateTime
yolo_combined <- yolo_flow_clean %>%
  full_join(yolo_mehg, by = c("SampleDateTime", "StationName", "SourceID", "Month")) %>%
  mutate(SiteTypeFilled = SiteType) %>% 
  arrange(SampleDateTime) %>% 
  fill(SiteTypeFilled, .direction = "up") %>% 
  select(-SiteType) %>% 
  rename(SiteType = SiteTypeFilled)

# Manually make sure that SiteType matches that of the next closest reported time
yolo_check <- yolo_combined %>% 
  dplyr::select(SampleDateTime, Flow, Result, SiteType, SiteType)
# It does in every scenario

# Join site types back to final tables
yolo_mehg_final <- yolo_mehg %>% 
  left_join(yolo_combined) %>% 
  dplyr::select(all_of(columns_mehg_select))

yolo_flow_clean2 <- yolo_flow_clean %>% 
  dplyr::select(-SiteType) %>% 
  left_join(yolo_combined) %>% 
  dplyr::select(all_of(columns_flow_select))
```



## Clean flow data
### Westervelt
Visualize flow data
```{r}
west_plot <- west_flow_clean %>%
  ggplot(aes(x = SampleDateTime, y = Flow)) +
    geom_area(fill="#69b3a2", alpha=0.5) +
    geom_line(color="#69b3a2") +
    ylab("Flow (cfs)")

# Turn it interactive with ggplotly
ggplotly(west_plot)
```

Remove incomplete tidal cycles
```{r}
west_flow_id <- unique_tidal_cycle(west_flow_clean)

west_flow_gaps <- identify_tidal_gaps(west_flow_id, 16)

# Match identified gaps from gapsts function with tide ID
west_gaps_id <- west_flow_gaps %>%
  filter(Gap == TRUE) %>%
  distinct(TideID)

# Remove these periods, and the first and last periods, from the flow data
west_flow_final <- west_flow_gaps %>%
  filter(TideID != "Flood1") %>%
  filter(TideID != "Ebb1738") %>%
  filter(TideID %are not% west_gaps_id$TideID)
```


### North Lindsey
Visualize flow data
```{r}
linds_plot <- linds_flow_clean %>% 
  
  ggplot(aes(x = SampleDateTime, y = Flow)) +
    geom_area(fill="#69b3a2", alpha=0.5) +
    geom_line(color="#69b3a2") +
    ylab("Flow (cfs)")

# Turn it interactive with ggplotly
ggplotly(linds_plot)
```

Remove incomplete tidal cycles
```{r}
linds_flow_id <- unique_tidal_cycle(linds_flow_clean2)

linds_flow_gaps <- identify_tidal_gaps(linds_flow_id, 16)

# Match identified gaps from gapsts function with tide ID
linds_gaps_id <- linds_flow_gaps %>%
  filter(Gap == TRUE) %>%
  distinct(TideID)

# Remove these periods, and the first and last periods, from the flow data
linds_flow_final <- linds_flow_gaps %>%
  filter(TideID != "Flood1") %>%
  filter(TideID != "Ebb1794") %>%
  filter(TideID %are not% linds_gaps_id$TideID)
```


### Yolo Bypass
Visualize flow data
```{r}
yolo_plot <- yolo_flow_clean %>%  # yolo_flow_final %>%
  ggplot(aes(x = SampleDateTime, y = Flow)) +
    geom_area(fill="#69b3a2", alpha=0.5) +
    geom_line(color="#69b3a2") +
    ylab("Flow (cfs)")

# Turn it interactive with ggplotly
ggplotly(yolo_plot)
```

Remove incomplete tidal cycles
```{r}
yolo_flow_id <- unique_tidal_cycle(yolo_flow_clean2)

yolo_flow_gaps <- identify_tidal_gaps(yolo_flow_id, 16)

# Match identified gaps from gapsts function with tide ID
yolo_gaps_id <- yolo_flow_gaps %>%
  filter(Gap == TRUE) %>%
  distinct(TideID)

# Remove these periods, and the first and last periods, from the flow data
yolo_flow_final <- yolo_flow_gaps %>%
  filter(TideID != "Flood1") %>%
  filter(TideID != "Ebb1865") %>%
  filter(TideID %are not% yolo_gaps_id$TideID)
```


# Visualize Methylmercury Data
```{r}
colors_graphs <- RColorBrewer::brewer.pal(n = 11, name = 'RdYlBu')
# [1] "#A50026" "#D73027" "#F46D43" "#FDAE61" "#FEE090" "#FFFFBF" "#E0F3F8" "#ABD9E9" "#74ADD1" "#4575B4" "#313695"

colors_graphs_2 <- colors_graphs[c(9,3)]
```

## Westervelt
```{r}
hist(west_mehg_final$Result)
```

```{r}
west_mehg_plot <- ggplot(data=west_mehg_final,
                         aes(x=month(Month, label=T, abbr=T),
                            y=Result,
                            color = SiteType)) +
   geom_point(size = 3) + 
   scale_y_continuous(expand=c(0,0), limits = c(0, 1.5), breaks = seq(0, 1.5, by = 0.25)) +
   scale_x_discrete(limits = as.character(month(1:12, label=T, abbr=T))) +
   scale_color_manual(name = "Site Type",
                      values = colors_graphs_2,
                     labels = function(x) str_wrap(x, width = 15)) +
   ylab("MeHg Result (ng/L)") +
   xlab("Month") +
   theme_light() + 
   theme(
     text = element_text(size=14, family="sans"),
     axis.title=element_text(size=14, face="bold"), 
     axis.text.x = element_text(angle = 0, hjust = 0.5),     
     panel.grid.major.y=element_line(size=1),
     panel.grid.minor=element_blank(),
     legend.text = element_text(size = 14), 
     legend.title = element_text(size = 14, face="bold"),
     legend.key.height=unit(1.0, "cm"))
west_mehg_plot

ggsave(filename = paste0(wd, '/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/Output/Fig 6.25 Westervelt MeHg.jpg'),
       plot = west_mehg_plot,
       width = 9,
       height = 6,
       units = "in",
       dpi = 125)
```


## North Lindsey
```{r}
hist(linds_mehg_final$Result)
```

```{r}
linds_mehg_plot <- ggplot(data=linds_mehg_final,
                           aes(x=month(Month, label=T, abbr=T),
                               y=Result,
                               color = SiteType)) +
   geom_point(size = 3) + 
   scale_y_continuous(expand=c(0,0), limits = c(0, 1.5), breaks = seq(0, 1.5, by = 0.25)) +
   scale_x_discrete(limits = as.character(month(1:12, label=T, abbr=T))) +
   scale_color_manual(name = "Site Type",
                      values = colors_graphs_2,
                     labels = function(x) str_wrap(x, width = 15)) +
   ylab("MeHg Result (ng/L)") +
   xlab("Month") +
   theme_light() + 
   theme(
     text = element_text(size=14, family="sans"),
     axis.title=element_text(size=14, face="bold"), 
     axis.text.x = element_text(angle = 0, hjust = 0.5),     
     panel.grid.major.y=element_line(size=1),
     panel.grid.minor=element_blank(),
     legend.text = element_text(size = 14), 
     legend.title = element_text(size = 14, face="bold"),
     legend.key.height=unit(1.0, "cm"))
linds_mehg_plot

ggsave(filename = paste0(wd, '/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/Output/Fig 6.26 North Lindsey MeHg.jpg'),
       plot = linds_mehg_plot,
       width = 9,
       height = 6,
       units = "in",
       dpi = 125)
```


## Yolo Bypass
```{r}
hist(yolo_mehg_final$Result)
```

```{r}
yolo_mehg_plot <- ggplot(data=yolo_mehg_final,
                           aes(x=month(Month, label=T, abbr=T),
                               y=Result,
                               color = SiteType)) +
   geom_point(size = 3) + 
   scale_y_continuous(expand=c(0,0), limits = c(0, 1.5), breaks = seq(0, 1.5, by = 0.25)) +
   scale_x_discrete(limits = as.character(month(1:12, label=T, abbr=T))) +
   scale_color_manual(name = "Site Type",
                      values = colors_graphs_2,
                     labels = function(x) str_wrap(x, width = 15)) +
   ylab("MeHg Result (ng/L)") +
   xlab("Month") +
   theme_light() + 
   theme(
     text = element_text(size=14, family="sans"),
     axis.title=element_text(size=14, face="bold"), 
     axis.text.x = element_text(angle = 0, hjust = 0.5),     
     panel.grid.major.y=element_line(size=1),
     panel.grid.minor=element_blank(),
     legend.text = element_text(size = 14), 
     legend.title = element_text(size = 14, face="bold"),
     legend.key.height=unit(1.0, "cm"))
yolo_mehg_plot

ggsave(filename = paste0(wd, '/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/Output/Fig 6.27 Yolo Bypass MeHg.jpg'),
       plot = yolo_mehg_plot,
       width = 9,
       height = 6,
       units = "in",
       dpi = 125)
```



# Analyze Data: Tidal Wetland Load Rate
Note: Assuming Ebb is synonymous with Drain and Flood is Synonymous with Source

## Formulas for summarizing data
```{r}
sum_flow <- function(flow_data, sitetype) {
  # flow_data - a dataframe, containing flow data in 15 minute time steps
  # sitetype - a character string, for example "Ebb" or "Flood"
  flow_summarized <- flow_data %>% 
  filter(SiteType == sitetype) %>% 
  # Multiply by seconds in 15 minutes to get volume per 15 minute period
  mutate(VolumeFifteen = abs(Flow) * 900) %>% 
  group_by(Month) %>% 
  summarise(MonthlyVolumeFifteen = sum(VolumeFifteen), n = n()) %>% 
  mutate(RepVolumeFifteen = MonthlyVolumeFifteen/n, # now in cf/15 minutes
         DaysInMonth = days_in_month(Month)) %>% 
  # Change to L/month
  mutate(MonthlyVolume = RepVolumeFifteen * 48 * DaysInMonth * 28.317) %>% 
  dplyr::select(Month, MonthlyVolume)
  return(flow_summarized)
}

sum_mehg <- function(mehg_data, flow_data, sitetype) {
  # mehg_data - a dataframe, containing methylmercury data
  # flow_data - a dataframe, containing flow data summarized by sum_flow()
  # sitetype - a character string, for example "Ebb" or "Flood"
  mehg_summarized <- mehg_data %>% 
  filter(SiteType == sitetype) %>% 
  group_by(Month) %>% 
  summarise(RepMeHg = median(Result)) %>% 
  left_join(flow_data) %>% 
  # L cancels, now in ng/month
  mutate(Load = RepMeHg * MonthlyVolume) %>%
  dplyr::select(Month, Load)
  return(mehg_summarized)
}

sum_load <- function(ebb_data, flood_data, area){
  # ebb_data - a dataframe containing ebb data summarized by sum_mehg()
  # flood_data - a dataframe containing flood data summarized by sum_mehg()
  # area - numeric value in meters squared
  load_summary <- ebb_data %>% 
  full_join(flood_data) %>% 
  mutate(MonthlyNetLoad = (EbbLoad - FloodLoad)/area) %>% 
  mutate(DaysInMonth = days_in_month(Month)) %>%
  # Now in ng/m2/day
  mutate(DailyNetLoad = MonthlyNetLoad/DaysInMonth) %>% 
  return(load_summary)
}
```

## Westervelt
```{r}
# Total study wetland area converting to meters squared
west_area <- 478 * 4047

# Summarize flow data
west_ebb_flow <- sum_flow(west_flow_final, "Ebb")

west_flood_flow <- sum_flow(west_flow_final, "Flood")
  
# Summarize methylmercury data & join with flow
west_ebb <- west_mehg_final %>% 
  sum_mehg(west_ebb_flow, "Ebb") %>% 
  rename(EbbLoad = Load)

west_flood <- west_mehg_final %>% 
  sum_mehg(west_flood_flow, "Flood") %>% 
  rename(FloodLoad = Load)

west_load <- sum_load(west_ebb, west_flood, west_area)
```


## North Lindsey
```{r}
# Total study wetland area converting to meters squared
linds_area <- 20.2 * 4047

# Summarize flow data
linds_ebb_flow <- sum_flow(linds_flow_final, "Ebb")

linds_flood_flow <- sum_flow(linds_flow_final, "Flood")
  
# Summarize methylmercury data & join with flow
linds_ebb <- linds_mehg_final %>% 
  sum_mehg(linds_ebb_flow, "Ebb") %>% 
  rename(EbbLoad = Load)

linds_flood <- linds_mehg_final %>% 
  sum_mehg(linds_flood_flow, "Flood") %>% 
  rename(FloodLoad = Load)

linds_load <- sum_load(linds_ebb, linds_flood, linds_area)
```


## Yolo Bypass
```{r}
# Total study wetland area converting to meters squared
yolo_area <- 728.5 * 4047

# Summarize flow data
yolo_ebb_flow <- sum_flow(yolo_flow_final, "Ebb")

yolo_flood_flow <- sum_flow(yolo_flow_final, "Flood")
  
# Summarize methylmercury data & join with flow
yolo_ebb <- yolo_mehg_final %>% 
  sum_mehg(yolo_ebb_flow, "Ebb") %>% 
  rename(EbbLoad = Load)

yolo_flood <- yolo_mehg_final %>% 
  sum_mehg(yolo_flood_flow, "Flood") %>% 
  rename(FloodLoad = Load)

yolo_load <- sum_load(yolo_ebb, yolo_flood, yolo_area)
```

## Export Tables of Rates (for allocations)
```{r}
all_rates <- list('Westervelt' = west_load, 'North Lindsey' = linds_load, 'Yolo Bypass' = yolo_load)

writexl::write_xlsx(all_rates, paste0(wd, "/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/Output/Tidal Wetlands Monthly Summary.xlsx"))
```


## All Tidal Wetlands Rate
```{r}
all_tidal <- bind_rows(west_load, linds_load, yolo_load)

all_tidal_rates <- all_tidal %>%
  group_by(Month) %>%
  summarise(MedianRate = median(DailyNetLoad))
  # %>% mutate(Month = month.abb[Month])

# Code below copies the table to clipboard for pasting into excel
# write.table(all_tidal_rates, "clipboard", sep="\t")
writexl::write_xlsx(all_tidal_rates, paste0(wd, "/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/Output/all_tidal_rates.xlsx"))
```

Using these rate(S) to calculate loads. See excel file "Tidal Wetlands_Load Calculations.xlsx" to see load calcualtions.


Below is the code to calculate the specific rates needed to evaluate the first 2 load calculation options.
```{r}
# Option 1
annual_median_rate <- median(all_tidal_rates$MedianRate)

# Option 2
tidal_warm_months <- all_tidal_rates %>% 
  filter(Month >= 3 & Month <= 9)

tidal_cool_months <- all_tidal_rates %>% 
  filter(Month >= 10 | Month <= 2)

warm_median_rate <- median(tidal_warm_months$MedianRate)

cool_median_rate <- median(tidal_cool_months$MedianRate)
```



Write to Excel to use in Source Load (04_Tidal Wetlands)
```{r}
tidal_wetlands_data <- list("Westervelt MeHg" = west_mehg_final, "North Lindsey MeHg" = linds_mehg_final, "Yolo Bypass MeHg" = yolo_mehg_final, "Westervelt Flow" = west_flow_final, "North Lindsey Flow" = linds_flow_final, "Yolo Bypass Flow" = yolo_flow_final)

writexl::write_xlsx(tidal_wetlands_data, paste0(wd, "/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/Output/Tidal Wetlands Clean Data.xlsx"))
```



*****CODE BELOW NOT USED IN FINAL REPORT, COMMENT OUT IF YOU DON'T WANT IT RUN*****
# Other Options Considered (See Appendix)
## Means instead of medians
```{r}
test_sum_mehg <- function(mehg_data, flow_data, sitetype) {
  mehg_summarized <- mehg_data %>% 
  filter(SiteType == sitetype) %>% 
  group_by(Month) %>% 
  summarise(RepMeHg = mean(Result)) %>% 
  left_join(flow_data) %>% 
  # L cancels, now in ng/month
  mutate(Load = RepMeHg * MonthlyVolume) %>%
  dplyr::select(Month, Load)
  return(mehg_summarized)
}

## Westervelt
# Summarize flow data
test_west_ebb_flow <- sum_flow(west_flow_final, "Ebb")

test_west_flood_flow <- sum_flow(west_flow_final, "Flood")
  
# Summarize methylmercury data & join with flow
test_west_ebb <- west_mehg_final %>% 
  test_sum_mehg(test_west_ebb_flow, "Ebb") %>% 
  rename(EbbLoad = Load)

test_west_flood <- west_mehg_final %>% 
  test_sum_mehg(test_west_flood_flow, "Flood") %>% 
  rename(FloodLoad = Load)

test_west_load <- sum_load(test_west_ebb, test_west_flood, west_area)

## North Lindsey
# Summarize flow data
test_linds_ebb_flow <- sum_flow(linds_flow_final, "Ebb")

test_linds_flood_flow <- sum_flow(linds_flow_final, "Flood")
  
# Summarize methylmercury data & join with flow
test_linds_ebb <- linds_mehg_final %>% 
  test_sum_mehg(test_linds_ebb_flow, "Ebb") %>% 
  rename(EbbLoad = Load)

test_linds_flood <- linds_mehg_final %>% 
  test_sum_mehg(test_linds_flood_flow, "Flood") %>% 
  rename(FloodLoad = Load)

test_linds_load <- sum_load(test_linds_ebb, test_linds_flood, linds_area)


## Yolo Bypass
# Summarize flow data
test_yolo_ebb_flow <- sum_flow(yolo_flow_final, "Ebb")

test_yolo_flood_flow <- sum_flow(yolo_flow_final, "Flood")
  
# Summarize methylmercury data & join with flow
test_yolo_ebb <- yolo_mehg_final %>% 
  test_sum_mehg(test_yolo_ebb_flow, "Ebb") %>% 
  rename(EbbLoad = Load)

test_yolo_flood <- yolo_mehg_final %>% 
  test_sum_mehg(test_yolo_flood_flow, "Flood") %>% 
  rename(FloodLoad = Load)

test_yolo_load <- sum_load(test_yolo_ebb, test_yolo_flood, yolo_area)


## All Tidal Wetlands Rate
test_all_tidal <- bind_rows(test_west_load, test_linds_load, test_yolo_load)

test_all_tidal_rates <- test_all_tidal %>% 
  group_by(Month) %>% 
  summarise(MeanRate = mean(DailyNetLoad)) 
  # %>% mutate(Month = month.abb[Month])

# Code below copies the table to clipboard for pasting into excel
# write.table(all_tidal_rates, "clipboard", sep="\t")
writexl::write_xlsx(test_all_tidal_rates, paste0(wd, "/Reeval_Source_Analysis/Loss Data/04_Tidal Wetlands/Output/test_all_tidal_rates.xlsx"))


# Option 1
test_annual_mean_rate <- mean(test_all_tidal_rates$MeanRate)

# Option 2
test_tidal_warm_months <- test_all_tidal_rates %>% 
  filter(Month >= 3 & Month <= 9)

test_tidal_cool_months <- test_all_tidal_rates %>% 
  filter(Month >= 10 | Month <= 2)

test_warm_mean_rate <- mean(test_tidal_warm_months$MeanRate)

test_cool_mean_rate <- mean(test_tidal_cool_months$MeanRate)
```


## Seasonal rates
```{r}
# Westervelt
# Monthly Median MeHg Concentration Difference
west_ebb_mm <- west_mehg_final %>%
  filter(SiteType == "Ebb") %>%
  group_by(Month) %>%
  summarise(MonthlyEbbMeHgMedian = median(Result))

west_flood_mm <- west_mehg_final %>%
  filter(SiteType == "Flood") %>%
  group_by(Month) %>%
  summarise(MonthlyFloodMeHgMedian = median(Result))

west_mm <- west_ebb_mm %>%
  left_join(west_flood_mm)


# Monthly Flows

#Original "west_ebb_flow" - only contains "Day" column and does not join with "west_flood_flow" below so copied "west_flood_flow" generation code

# west_ebb_flow <- west_flow_final %>%
#   filter(SiteType == "Ebb") %>%
#   # multiply by seconds in 15 minutes to get total volume per 15 minute time step (in cubic feet)
#   mutate(TotalFlow = Flow * 900) %>%
#   mutate(DaysInMonth = days_in_month(Month)) %>%
#   mutate(Day = day(SampleDateTime)) %>%
#   group_by(Day) %>%
#   summarise

# corrected using west_flood_flow as an example since they are merged later
west_ebb_flow <- west_flow_final %>%
  filter(SiteType == "Ebb") %>%
  group_by(Month) %>%
  summarise(MonthlyFloodFlowMedian = median(Flow, na.rm = TRUE)) %>%
  mutate(DaysInMonth = days_in_month(Month)) %>%
  # Absolute value of flood flow
  mutate(MonthlyEbbFlow = abs(MonthlyFloodFlowMedian * 28.317 * DaysInMonth))
  # Represents a daily flow on any given day in that month

# total flood volume for the study period in cubic feet, divided by periods of 15 minutes so it's ebb volume per 15 minutes
# west_ebb_volume <- sum(west_ebb_flow$TotalFlow)/nrow(west_ebb_flow)

# Remove na's in median, missing flow data
west_flood_flow <- west_flow_final %>%
  filter(SiteType == "Flood") %>%
  group_by(Month) %>%
  summarise(MonthlyFloodFlowMedian = median(Flow, na.rm = TRUE)) %>%
  mutate(DaysInMonth = days_in_month(Month)) %>%
  # Absolute value of flood flow
  mutate(MonthlyFloodFlow = abs(MonthlyFloodFlowMedian * 28.317 * DaysInMonth))
  # Represents a daily flow on any given day in that month

west_flow_mm <- west_ebb_flow %>%
  left_join(west_flood_flow) %>%
  dplyr::select(Month, MonthlyEbbFlow, MonthlyFloodFlow)


# North Lindsey
# Monthly Median MeHg Concentration Difference
linds_ebb_mm <- linds_mehg_final %>%
  filter(SiteType == "Ebb") %>%
  group_by(Month) %>%
  summarise(MonthlyEbbMeHgMedian = median(Result))

linds_flood_mm <- linds_mehg_final %>%
  filter(SiteType == "Flood") %>%
  group_by(Month) %>%
  summarise(MonthlyFloodMeHgMedian = median(Result))

linds_mm <- linds_ebb_mm %>%
  left_join(linds_flood_mm)


# Monthly Flows
linds_ebb_flow <- linds_flow_final %>%
  filter(SiteType == "Ebb") %>%
  group_by(Month) %>%
  summarise(MonthlyEbbFlowMedian = median(Flow, na.rm = TRUE)) %>%
  mutate(DaysInMonth = days_in_month(Month)) %>%
  mutate(MonthlyEbbFlow = (MonthlyEbbFlowMedian * 28.317 * DaysInMonth))

# Remove na's in median, missing flow data
linds_flood_flow <- linds_flow_final %>%
  filter(SiteType == "Flood") %>%
  group_by(Month) %>%
  summarise(MonthlyFloodFlowMedian = median(Flow, na.rm = TRUE)) %>%
  mutate(DaysInMonth = days_in_month(Month)) %>%
  # Absolute value of flood flow
  mutate(MonthlyFloodFlow = abs(MonthlyFloodFlowMedian * 28.317 * DaysInMonth))
  # Represents a daily flow on any given day in that month

linds_flow_mm <- linds_ebb_flow %>%
  left_join(linds_flood_flow) %>%
  dplyr::select(Month, MonthlyEbbFlow, MonthlyFloodFlow)


# Yolo
# Monthly Median MeHg Concentration Difference
yolo_ebb_mm <- yolo_mehg_final %>%
  filter(SiteType == "Ebb") %>%
  group_by(Month) %>%
  summarise(MonthlyEbbMeHgMedian = median(Result))

yolo_flood_mm <- yolo_mehg_final %>%
  filter(SiteType == "Flood") %>%
  group_by(Month) %>%
  summarise(MonthlyFloodMeHgMedian = median(Result))

yolo_mm <- yolo_ebb_mm %>%
  left_join(yolo_flood_mm)


# Monthly Flows
yolo_ebb_flow <- yolo_flow_final %>%
  filter(SiteType == "Ebb") %>%
  group_by(Month) %>%
  summarise(MonthlyEbbFlowMedian = median(Flow, na.rm = TRUE)) %>%
  mutate(DaysInMonth = days_in_month(Month)) %>%
  mutate(MonthlyEbbFlow = (MonthlyEbbFlowMedian * 28.317 * DaysInMonth))

# Remove na's in median, missing flow data
yolo_flood_flow <- yolo_flow_final %>%
  filter(SiteType == "Flood") %>%
  group_by(Month) %>%
  summarise(MonthlyFloodFlowMedian = median(Flow, na.rm = TRUE)) %>%
  mutate(DaysInMonth = days_in_month(Month)) %>%
  # Absolute value of flood flow
  mutate(MonthlyFloodFlow = abs(MonthlyFloodFlowMedian * 28.317 * DaysInMonth))
# Represents a daily flow on any given day in that month

yolo_flow_mm <- yolo_ebb_flow %>%
  left_join(yolo_flood_flow) %>%
  dplyr::select(Month, MonthlyEbbFlow, MonthlyFloodFlow)


# All
# Join all 3 tables
tidal_wetlands <- yolo_flow_mm %>%
  full_join(yolo_mm) %>%
  full_join(linds_flow_mm) %>%
  full_join(linds_mm) %>%
  full_join(west_flow_mm) %>%
  full_join(west_mm) %>%
  # Average 3 wetland monthly averages & daily flows
  group_by(Month) %>%
  summarise(MonthlyFloodFlow = sum(MonthlyFloodFlow, na.rm = TRUE),
            MonthlyEbbFlow = sum(MonthlyEbbFlow, na.rm = TRUE),
            MedianFloodMeHg = median(MonthlyFloodMeHgMedian, na.rm = TRUE),
            MedianEbbMeHg = median(MonthlyEbbMeHgMedian, na.rm = TRUE)) %>%
  ungroup()

# Calculate daily loads, in ng/month
tidal_loads <- tidal_wetlands %>%
  mutate(MonthlyEbbLoad = MonthlyEbbFlow * MedianEbbMeHg,
         MonthlyFloodLoad = MonthlyFloodFlow * MedianFloodMeHg) %>%
  mutate(MonthlyLoad = MonthlyEbbLoad - MonthlyFloodLoad)

# Calculate total area
tidal_area <- west_area + linds_area + yolo_area

# Calculate Rate (units now in ng/m2/day)
tidal_rate <- (sum(tidal_loads$MonthlyLoad)/tidal_area)/365

# Warm months
tidal_warm_months <- tidal_loads %>%
  filter(Month >= 3 & Month <= 9)

tidal_warm_rate <- (sum(tidal_warm_months$MonthlyLoad)/tidal_area)/365

# Cool months
tidal_cool_months <- tidal_loads %>%
  filter(Month >= 10 | Month <= 2)

tidal_cool_rate <- (sum(tidal_cool_months$MonthlyLoad)/tidal_area)/365
# In both cases the rate is higher in warm months than cool months.
```


## Grouping all the data together
```{r}
# Join tables
# MeHg
all_tidal_wetlands_mehg <- west_mehg_final %>%
  full_join(linds_mehg_final) %>%
  full_join(yolo_mehg_final)
#Flow
all_tidal_wetlands_flow <- west_flow_final %>%
  full_join(linds_flow_final) %>%
  full_join(yolo_flow_final)

# Summarise flood
all_flood_mehg <- all_tidal_wetlands_mehg %>%
  filter(SiteType == "Flood") %>%
  group_by(Month) %>%
  summarise(MonthlyFloodMeHgMedian = median(Result))

all_flood_flow <- all_tidal_wetlands_flow %>%
  filter(SiteType == "Flood") %>%
  group_by(Month) %>%
  summarise(MonthlyFloodFlowMedian = median(Flow, na.rm = TRUE)) %>%
  mutate(DaysInMonth = days_in_month(Month)) %>%
  # Absolute value of flood flow
  mutate(MonthlyFloodFlow = abs(MonthlyFloodFlowMedian * 28.317 * DaysInMonth))

# Summarise ebb
all_ebb_mehg <- all_tidal_wetlands_mehg %>%
  filter(SiteType == "Ebb") %>%
  group_by(Month) %>%
  summarise(MonthlyEbbMeHgMedian = median(Result))

all_ebb_flow <- all_tidal_wetlands_flow %>%
  filter(SiteType == "Ebb") %>%
  group_by(Month) %>%
  summarise(MonthlyEbbFlowMedian = median(Flow, na.rm = TRUE)) %>%
  mutate(DaysInMonth = days_in_month(Month)) %>%
  # Absolute value of flood flow
  mutate(MonthlyEbbFlow = abs(MonthlyEbbFlowMedian * 28.317 * DaysInMonth))

all_joined <- all_flood_mehg %>%
  full_join(all_flood_flow) %>%
  full_join(all_ebb_mehg) %>%
  full_join(all_ebb_flow) %>%
  mutate(MonthlyEbbLoad = MonthlyEbbFlow * MonthlyEbbMeHgMedian,
         MonthlyFloodLoad = MonthlyFloodFlow * MonthlyFloodMeHgMedian) %>%
  mutate(MonthlyLoad = MonthlyEbbLoad - MonthlyFloodLoad)

# Warm months
all_tidal_warm_months <- all_joined %>%
  filter(Month >= 6 & Month <= 11)

all_tidal_warm_rate <- (sum(all_tidal_warm_months$MonthlyLoad)/tidal_area)/365

# Cool months
all_tidal_cool_months <- all_joined %>%
  filter(Month == 12 | Month <= 5)

tidal_cool_rate <- (sum(all_tidal_cool_months$MonthlyLoad)/tidal_area)/365

all_tidal_rate <- (sum(all_joined$MonthlyLoad)/tidal_area)
```
## Error: <text>:18:1: unexpected '<'
## 17: 
## 18: <
##     ^

The R session information (including the OS info, R version and all packages used):

    sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] rstudioapi_0.13    knitr_1.39         magrittr_2.0.3     tidyselect_1.1.2  
##  [5] R6_2.5.1           rlang_1.0.2        fastmap_1.1.0      fansi_1.0.3       
##  [9] stringr_1.4.0      highr_0.9          dplyr_1.0.9        tools_4.2.2       
## [13] xfun_0.31          utf8_1.2.2         DBI_1.1.2          cli_3.3.0         
## [17] htmltools_0.5.2    ellipsis_0.3.2     assertthat_0.2.1   yaml_2.3.5        
## [21] readxl_1.4.0       digest_0.6.29      tibble_3.1.7       lifecycle_1.0.1   
## [25] crayon_1.5.1       RColorBrewer_1.1-3 purrr_0.3.4        vctrs_0.4.1       
## [29] glue_1.6.2         evaluate_0.15      rmarkdown_2.14     stringi_1.7.6     
## [33] compiler_4.2.2     pillar_1.7.0       cellranger_1.1.0   generics_0.1.2    
## [37] writexl_1.4.0      pkgconfig_2.0.3
    Sys.time()
## [1] "2023-12-26 15:56:36 PST"