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"