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

---
title: "Nontidal Wetlands"
author: "Mercury Program and Basin Planning Unit"
date: "2/4/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, include=FALSE}
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(matrixStats)

# Had issue with runtime:shiny resetting the WD of rproj so used getActiveProject()
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 to select
columns_select <- c("SourceRow", "SampleDate", "StationName", "SourceID", "MDL", "RL", "UnitMeHg", "Result", "ResultQualCode", "SiteType", "Comments", "Month", "Year", "DaysInMonth")
```



## BREW
```{r}
brew <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Source Data/02b_Nontidal Wetlands/USGS BREW_Aq MeHg Data.xlsx"), sheet = 2)
```

```{r}
# Initial cleaning of data & adding flow
brew_clean1 <- brew %>%
  clean_names() %>%
  mutate(sample_date = as.Date(sample_date, origin = "1899-12-30")) %>% 
  mutate(location_code = case_when(location_code == "inlet" ~ "Inflow",
                                location_code == "outlet" ~ "Outflow")) %>%
  mutate(res_qual_code = case_when(res_qual_code == "None" ~ NA))

# Add together particulate and filtered concentrations to get total unfiltered methylmercury concentrations
brew_unfiltered <- brew_clean1 %>% 
  group_by(sample_date, collection_time) %>% 
  summarise_at(vars(result), sum) %>% 
  ungroup %>% 
  rename(Result = result)

# Join unfiltered results back to table
brew_clean2 <- brew_unfiltered %>% 
  full_join(brew_clean1) %>% 
  mutate(Unique = paste(sample_date, collection_time, Result)) %>% 
  distinct(Unique, .keep_all=T)
```


```{r}
# Rename and select final columns
brew_final <- brew_clean2 %>% 
  dplyr::rename(
    SiteType = location_code,
    SampleDate = sample_date,
    StationName = station_code,
    ResultQualCode = res_qual_code,
    UnitMeHg = unit_name,
    SourceID = project_code,
    MDL = mdl,
    RL = rl,
    Comments = lab_result_comments
    ) %>%
  # Add missing columns
  add_column(SourceRow = 1:nrow(.)) %>% 
  mutate(Month = month(SampleDate),
         Year = year(SampleDate)) %>%
  mutate(DaysInMonth = days_in_month(Month)) %>%
  dplyr::select(all_of(columns_select))
```


## Twitchell Island
### Methylmercury
```{r}
twitchell_mehg <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Source Data/02b_Nontidal Wetlands/2008 CALFED Task 5.3a_Aq MeHg and Flow Data.xlsx"),
                    sheet = 2)
```

```{r}
twitchell_mehg_final <- twitchell_mehg %>%
  rowwise() %>% 
  mutate(`West Pond Outflow` = median(c(`West Pond Out 1`, `West Pond Out 2`), na.rm = TRUE),
         `East Pond Outflow` = median(c(`East Pond Out 1`, `East Pond Out 2`), na.rm = TRUE)) %>% 
  mutate(SampleDate = mdy(Date)) %>% 
  dplyr::select(SampleDate, `Source Water`, `West Pond Outflow`, `East Pond Outflow`) %>% 
  pivot_longer(cols = c(`West Pond Outflow`, `East Pond Outflow`, `Source Water`), names_to = "StationName", values_to = "Result") %>% 
  mutate(SiteType = case_when(StationName == "West Pond Outflow" | StationName == "East Pond Outflow" ~ "Outflow",
                              StationName == "Source Water" ~ "Inflow"),
         Month = month(SampleDate),
         Year = year(SampleDate),
         SourceID = "CALFED 2008 Task 5.3a (Twitchell Island)",
         MDL = NA,
         RL = NA,
         ResultQualCode = "=",
         UnitMeHg = "ng/L",
         Comments = NA_character_,
         SourceRow = 1:nrow(.)
         ) %>% 
  mutate(DaysInMonth = days_in_month(Month)) %>% 
  dplyr::select(all_of(columns_select))
```

### Flow
```{r}
twitchell_flow <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Source Data/02b_Nontidal Wetlands/2008 CALFED Task 5.3a_Aq MeHg and Flow Data.xlsx"),
                    sheet = 3)
```

```{r}
# Note: when converting from DOY to day, assuming January 1st is day 1 (based on spreadsheet, where DOY 296 in 2003 is October 23)
twitchell_flow_clean <- twitchell_flow %>% 
  clean_names() %>% 
  mutate(SampleDate = as.Date(doy, origin = "2002-12-31")) %>%
  mutate(Month = month(SampleDate)) %>% 
  rename(percent_lost = percent_lost_to_evap_or_groundwater,
         Year = year) %>% 
  dplyr::select(SampleDate, Year, Month, total_out_cfs, total_in_cfs, percent_lost)

# Median of flow in, flow out, and percent lost (data for October, November, December; this flow represents winter months)
twitchell_winter_inflow_cfs <- median(twitchell_flow_clean$total_in_cfs)
twitchell_winter_outflow_cfs <- median(twitchell_flow_clean$total_out_cfs)
twitchell_winter_percentlost <- median(twitchell_flow_clean$percent_lost)

twitchell_summer_inflow_cfs <- 0.25 * twitchell_winter_inflow_cfs
twitchell_summer_outflow_cfs <- 0.186 * twitchell_winter_outflow_cfs

```

# Visualize 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)]
```


## BREW
```{r}
hist(brew_final$Result)
```

```{r}
(brew_mehg_plot <- ggplot(
  data=brew_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, 6), breaks = seq(0, 6, by = 1)) +
  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"))
  )

ggsave(filename = paste0(wd, '/Reeval_Source_Analysis/Source Data/02b_Nontidal Wetlands/Output/Figures/Fig 6.9 BREW MeHg.jpg'),
       plot = brew_mehg_plot,
       width = 9,
       height = 6,
       units = "in",
       dpi = 125)
```


## Twitchell
```{r}
hist(twitchell_mehg_final$Result)
```

```{r}
(twitchell_mehg_plot <- ggplot(
  data=twitchell_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, 6), breaks = seq(0, 6, by = 1)) +
   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"))
   )

ggsave(filename = paste0(wd, '/Reeval_Source_Analysis/Source Data/02b_Nontidal Wetlands/Output/Figures/Fig 6.10 Twitchell MeHg.jpg'),
       plot = twitchell_mehg_plot,
       width = 9,
       height = 6,
       units = "in",
       dpi = 125)
```

# Analyze Data: Nontidal Wetland Load Rate
## BREW
```{r}
brew_total_volume <- 80 * 1233000 # acre feet times conversion factor to get L

brew_total_area <- 87.7 * 4046.86 # acres times conversion factor to get m2 

brew_vwm <- c(12, 1:3)

# cfs * conversion factor to get L/day
brew_vw_flow <- 0.5 * 2446571.52

brew_wm <- c(9:11, 4:5)
brew_w_flow <- 0.125 * 2446571.52

brew_dm <- c(6, 7, 8)
brew_d_flow <- 0
```

```{r}
brew_outflow_mm <- brew_final %>%
  filter(SiteType == "Outflow") %>% 
  group_by(Month) %>%
  summarise(MonthlyOutflowMeHgMedian = median(Result))

brew_inflow_mm <- brew_final %>%
  filter(SiteType == "Inflow") %>% 
  group_by(Month) %>%
  summarise(MonthlyInflowMeHgMedian = median(Result))

brew_mm <- brew_outflow_mm %>%
  left_join(brew_inflow_mm) %>% 
  fill(MonthlyInflowMeHgMedian, .direction = "down") %>% 
  add_row(Month = c(6:8)) %>% 
  mutate(DaysInMonth = days_in_month(Month)) %>% 
  # There is no inflow MeHg data for May so Board staff assumed the concentration was the same as April inflow water concentration
  mutate(Flow = case_when(Month %in% brew_vwm ~ brew_vw_flow,
                          Month %in% brew_wm ~ brew_w_flow,
                          Month %in% brew_dm ~ brew_d_flow)) %>% 
  mutate(Volume = Flow * DaysInMonth) %>% 
  mutate(MonthlyInflowLoad = Volume * MonthlyInflowMeHgMedian,
         MonthlyOutflowLoad = Volume * MonthlyOutflowMeHgMedian) %>% 
  # Assuming there is no flow (and therefore no load) in or out in summer months
  mutate(MonthlyInflowLoad = case_when(is.na(MonthlyInflowLoad) ~ 0,
                                       TRUE ~ MonthlyInflowLoad),
         MonthlyOutflowLoad = case_when(is.na(MonthlyOutflowLoad) ~ 0,
                                        TRUE ~ MonthlyOutflowLoad)) %>%
  mutate(MonthlyNetLoad = MonthlyOutflowLoad - MonthlyInflowLoad) %>% # in ng/month
  mutate(DailyNetLoadRate = MonthlyNetLoad / brew_total_area / DaysInMonth) %>%  # in ng/m2/day
  arrange(Month)

brew_annual_load = sum(brew_mm$MonthlyNetLoad)
brew_annual_load_rate = brew_annual_load / brew_total_area / 365
# 0.299 ng/m2/day
```


## Twitchell
```{r}
# "The West Pond is 27,534.8 square meters and the East Pond is 26,472.9 square meters, or approximately 2.7 hectares." (pg. 3)
twitchell_total_area <- 27534.8 + 26472.9 # in m2 

twitchell_wintermonths <- c(11:12, 1:4)
twitchell_summermonths <- c(5:10)

# cfs * conversion factor to get L/day
twitchell_winter_inflow <- twitchell_winter_inflow_cfs * 2446571.52
twitchell_summer_inflow <- twitchell_summer_inflow_cfs * 2446571.52

twitchell_winter_outflow <- twitchell_winter_outflow_cfs * 2446571.52
twitchell_summer_outflow <- twitchell_summer_outflow_cfs * 2446571.52
```

```{r}
# for missing months, assume the median is the same as the median of the "sandwiching" months' medians
twitchell_outflow_mm <- twitchell_mehg_final %>%
  filter(SiteType == "Outflow") %>% 
  group_by(Month) %>%
  summarise(MonthlyOutflowMeHgMedian = median(Result))

t_outflow_jan <- twitchell_outflow_mm %>% 
  filter(Month == 12 | Month == 2)

t_outflow_oct <- twitchell_outflow_mm %>% 
  filter(Month == 9 | Month == 11)

# inflow
twitchell_inflow_mm <- twitchell_mehg_final %>%
  filter(SiteType == "Inflow") %>% 
  group_by(Month) %>%
  summarise(MonthlyInflowMeHgMedian = median(Result))

t_inflow_jan <- twitchell_inflow_mm %>% 
  filter(Month == 12 | Month == 2)

t_inflow_oct <- twitchell_inflow_mm %>% 
  filter(Month == 9 | Month == 11)

# join
twitchell_mm <- twitchell_outflow_mm %>%
  left_join(twitchell_inflow_mm) %>% 
  add_row(Month = c(1, 10)) %>% 
  rowwise() %>%
  mutate(MonthlyOutflowMeHgMedian = 
           case_when(Month == 1 ~ median(t_outflow_jan$MonthlyOutflowMeHgMedian),
                     Month == 10 ~ median(t_outflow_oct$MonthlyOutflowMeHgMedian),
                     TRUE ~ MonthlyOutflowMeHgMedian),
         MonthlyInflowMeHgMedian = 
           case_when(Month == 1 ~ median(t_inflow_jan$MonthlyInflowMeHgMedian),
                     Month == 10 ~ median(t_inflow_oct$MonthlyInflowMeHgMedian),
                     TRUE ~ MonthlyInflowMeHgMedian)) %>%
  arrange(Month) %>% 
  mutate(DaysInMonth = days_in_month(Month)) %>% 
  # There is no inflow MeHg data for May so Board staff assumed the concentration was the same as April inflow water concentration
  mutate(Inflow = 
           case_when(Month %in% twitchell_wintermonths ~ twitchell_winter_inflow,
                     Month %in% twitchell_summermonths ~ twitchell_summer_inflow),
         Outflow = 
           case_when(Month %in% twitchell_wintermonths ~ twitchell_winter_outflow,
                     Month %in% twitchell_summermonths ~ twitchell_summer_outflow)) %>% 
  mutate(InflowVolume = Inflow * DaysInMonth,
         OutflowVolume = Outflow * DaysInMonth) %>% 
  mutate(MonthlyInflowLoad = InflowVolume * MonthlyInflowMeHgMedian,
         MonthlyOutflowLoad = OutflowVolume * MonthlyOutflowMeHgMedian) %>% 
  mutate(MonthlyNetLoad = MonthlyOutflowLoad - MonthlyInflowLoad) %>% # in ng/month
  mutate(DailyNetLoadRate = MonthlyNetLoad / twitchell_total_area / DaysInMonth) # in ng/m2/day

twitchell_annual_load = sum(twitchell_mm$MonthlyNetLoad)
twitchell_annual_load_rate = twitchell_annual_load / twitchell_total_area / 365
# 5.479 ng/m2/day
```

## Export Tables of Rates (for allocations)
```{r}
# setwd("R:/RB5/R5SSections/TMDL Basin NPS Delta/Units/Mercury Metals TMDL/TMDL_Delta_Review/Revised Staff Report/Reeval_Allocations")
# all_rates <- list('Twitchell' = twitchell_mm, 'Brew' = brew_mm)
# write.xlsx(all_rates, file = 'Gross Flux Calculations/Nontidal Wetland Rates.xlsx')
# setwd(wd)

# Added this code to save output in output folder
writexl::write_xlsx(list('Twitchell' = twitchell_mm, 'Brew' = brew_mm),
paste0(wd, '/Reeval_Source_Analysis/Source Data/02b_Nontidal Wetlands/Output/Nontidal Wetland Rates.xlsx'))
```

## Combine load from both wetlands
```{r}
twitchell_monthly_rates <- twitchell_mm %>% 
  dplyr::select(Month, DailyNetLoadRate) %>% 
  rename(TwitchellRate = DailyNetLoadRate)

brew_monthly_rates <- brew_mm %>% 
  dplyr::select(Month, DailyNetLoadRate) %>% 
  rename(BREWRate = DailyNetLoadRate)

all_rates <- twitchell_monthly_rates %>% 
  full_join(brew_monthly_rates) %>% 
  mutate(AllRates = (median(c(TwitchellRate, BREWRate))))

# Code below copies the table to clipboard for pasting into excel
# write.table(all_rates, "clipboard", sep="\t")
writexl::write_xlsx(all_rates, paste0(wd, '/Reeval_Source_Analysis/Source Data/02b_Nontidal Wetlands/Output/all_rates.xlsx'))

total_annual_load = brew_annual_load + twitchell_annual_load

total_area = brew_total_area + twitchell_total_area

total_annual_load_rate = total_annual_load / total_area / 365
```



*****CODE BELOW NOT USED IN FINAL REPORT, COMMENT OUT IF YOU DON'T WANT IT RUN*****
# Other Options Considered for BREW (See Appendix)
## Option 1
Separating year 1 and 2
Year 1
```{r}
# Total wetland volume (average in acre feet) provided in BREW Report, converting to L
total_volume <- 80 * 1233000

# Fill in September over 16 days
fill_year1 <- brew_final %>% filter(Year == 2014 & Month == 9 & SiteType == "Inflow")

# L * ng/L = Total inflow load in ng
fill_load_year1 <- total_volume * median(fill_year1$Result)

# Drain in May over 15 days, outflow sampling occurred right before draw-down in late April
outflow_year1 <- brew_final %>% filter(SampleDate == "2015-04-29")

# Total outflow load in ng
outflow_load_year1 <- total_volume * median(outflow_year1$Result)

# Total net load in ng/year (only on discharge time per year)
total_dry_year <- (outflow_load_year1 - fill_load_year1)
```

Year 2
```{r}
# Flow through starting in September, lasting 16 days in September
# Assuming flow through went halfway (16 days) into December
# Total days then are 16 (September) + 16 (December) + 31 (October) + 30 (November)
flowthrough_days <- 16 + 16 + 31 + 30

# 0.5 cfs of flow through flow, times conversion factor to get L/day
flowthrough_volume <- 0.5 * 2446571.52

flowthrough_inflow_year2 <- brew_final %>% 
  filter(Year == 2015 & Month >= 9 & Month <= 12 & SiteType == "Inflow") %>% 
  group_by(Month) %>% 
  mutate(MonthlyInflowResult = median(Result)) %>% 
  distinct(MonthlyInflowResult) %>%
  mutate(DaysInMonth = case_when(Month == 9 ~ 16,
                                 Month == 10 ~ 31,
                                 Month == 11 ~ 30,
                                 Month == 12 ~ 16)) %>% 
  mutate(MonthlyFlow = flowthrough_volume * DaysInMonth) %>% # in L per those numbers of days
  mutate(MonthlyInflowLoad = MonthlyFlow * MonthlyInflowResult) # in ng per those numbers of days

flowthrough_outflow_year2 <- brew_final %>%
  filter(Year == 2015 & Month >= 9 & Month <= 12 & SiteType == "Outflow") %>% 
  group_by(Month) %>% 
  mutate(MonthlyOutflowResult = median(Result)) %>% 
  distinct(MonthlyOutflowResult)

flowthrough_year2 <- flowthrough_inflow_year2 %>% 
  full_join(flowthrough_outflow_year2) %>% 
  mutate(MonthlyOutflowResult = case_when(Month == 10 ~ median(c(2.73043707, 0.17698460)), # assuming the October outflow concentration is represented by the median of September and November outflow concentration
                                        TRUE ~ MonthlyOutflowResult)) %>% 
  mutate(MonthlyOutflowLoad = MonthlyFlow * MonthlyOutflowResult) %>% 
  mutate(DailyNetLoad = (MonthlyOutflowLoad - MonthlyInflowLoad)) # in ng per those number of days

# Total flow through
flowthrough_load <- sum(flowthrough_year2$DailyNetLoad) # in ng per those number of days


# Outflow load in May 2016
outflow_year2 <- brew_final %>% filter(Year == 2016 & Month == 5)

# ng/day
outflow_load_year2 <- total_volume * median(outflow_year2$Result)

# ng/year in wet year
total_wet_year <- outflow_load_year2 + flowthrough_load

# ng/m2/day
option1_rate <- (median(c(total_wet_year, total_dry_year))) / total_area / 365
# 0.459 ng/m2/day
```


## Option 2
Pooling together both years
```{r}
# Monthly Median MeHg Concentration Difference
test_brew_outflow_mm <- brew_final %>%
  filter(SiteType == "Outflow") %>% 
  group_by(Year, Month) %>%
  summarise(MonthlyOutflowMeHgMedian = median(Result))

test_brew_inflow_mm <- brew_final %>%
  filter(SiteType == "Inflow") %>% 
  group_by(Year, Month) %>%
  summarise(MonthlyInflowMeHgMedian = median(Result))

test_brew_mm <- brew_outflow_mm %>%
  left_join(brew_inflow_mm) 
# %>%
  # There is no inflow MeHg data for May so we will not use this month
  # filter(!is.na(MonthlyInflowMeHgMedian))

# Monthly Flows

# test_brew_outflow_flow <- brew_final %>%
#   filter(SiteType == "Outflow") %>% 
#   group_by(Year, Month) %>%
#   summarise(MonthlyOutflowFlowMedian = median(Flow)) %>% 
#   mutate(DaysInMonth = days_in_month(Month)) %>% 
#   mutate(MonthlyOutflowVolume = (MonthlyOutflowFlowMedian * 86400 * DaysInMonth * 28.317))

# test_brew_inflow_flow <- brew_final %>%
#   filter(SiteType == "Inflow") %>% 
#   group_by(Month) %>%
#   summarise(MonthlyInflowFlowMedian = median(Flow)) %>% 
#   mutate(DaysInMonth = days_in_month(Month)) %>% 
#   mutate(MonthlyInflowVolume = (MonthlyInflowFlowMedian * 86400 * DaysInMonth * 28.317))

# test_brew_flow_mm <- brew_outflow_flow %>%
#   left_join(brew_inflow_flow) %>%
#   # There is no associated inflow MeHg data for May so we will not use this month
#   # filter(!is.na(MonthlyInflowFlowMedian)) %>% 
#   dplyr::select(Year, Month, MonthlyInflowVolume, MonthlyOutflowVolume)

# Total study wetland acreage converting to meters squared
test_brew_area <- 87.7 * 4046.86

# Calculate daily load (units now in ng/day)
# test_brew_load <- brew_flow_mm %>% 
#   full_join(brew_mm) %>% 
#   mutate(MonthlyOutflowLoad = MonthlyOutflowVolume * MonthlyOutflowMeHgMedian,
#          MonthlyInflowLoad = MonthlyInflowVolume * MonthlyInflowMeHgMedian) 
# %>%
#   # multiply by days in each month to get load for each month
#   # subtract outflow - inflow to get total monthly loads
#   mutate(MonthlyLoad = DailyOutflowLoad - DailyInflowLoad)

# test_annual_inflow <- sum(brew_load$MonthlyInflowLoad, na.rm = TRUE)/3
# 
# test_annual_outflow <- sum(brew_load$MonthlyOutflowLoad)/3 #divide by 3 years

# Calculate Rate (units now in ng/m2/day)

# option2_rate <- (annual_outflow - annual_inflow)/brew_acreage/365

# Comes out as 1.216 when using medians for everything (except the final calculation is a mean because median resulted in 0 because there is 0 flow for most months)
# Comes out as 1.117 when using means for everything
# Comes out as 0 when using geomeans for everything
```


## Option 3
Seasonal Rates
```{r}
# brew_warm_months <- brew_load %>% 
#   filter(Month >= 3 & Month <= 9)
# 
# brew_warm_rate <- mean(brew_warm_months$DailyLoad)/brew_acreage
# # When using all means, 2.964
# # When using all geomeans, 0
# # When using all medians, 3.07
# 
# brew_cool_months <- brew_load %>% 
#   filter(Month >= 10 | Month <= 2)
# 
# brew_cool_rate <- mean(brew_cool_months$DailyLoad)/brew_acreage
# # When using all means, 0.008
# # When using all geomeans, 0
# # When using medians, 0.08

# In both cases the rate is higher in warm months than cool months.
```
## 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     
## 
## other attached packages:
##  [1] lubridate_1.8.0    plotly_4.10.0      readxl_1.4.0       actuar_3.2-2      
##  [5] NADA_1.6-1.1       forcats_0.5.1      stringr_1.4.0      dplyr_1.0.9       
##  [9] purrr_0.3.4        readr_2.1.2        tidyr_1.2.0        tibble_3.1.7      
## [13] ggplot2_3.3.6      tidyverse_1.3.1    fitdistrplus_1.1-8 survival_3.4-0    
## [17] MASS_7.3-58.1     
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-45    assertthat_0.2.1   digest_0.6.29      utf8_1.2.2        
##  [5] R6_2.5.1           cellranger_1.1.0   backports_1.4.1    reprex_2.0.1      
##  [9] evaluate_0.15      httr_1.4.3         highr_0.9          pillar_1.7.0      
## [13] rlang_1.0.2        lazyeval_0.2.2     rstudioapi_0.13    data.table_1.14.2 
## [17] Matrix_1.5-1       rmarkdown_2.14     labeling_0.4.2     splines_4.2.2     
## [21] htmlwidgets_1.5.4  munsell_0.5.0      broom_0.8.0        compiler_4.2.2    
## [25] modelr_0.1.8       xfun_0.31          pkgconfig_2.0.3    htmltools_0.5.2   
## [29] tidyselect_1.1.2   viridisLite_0.4.0  fansi_1.0.3        crayon_1.5.1      
## [33] tzdb_0.3.0         dbplyr_2.2.0       withr_2.5.0        grid_4.2.2        
## [37] jsonlite_1.8.0     gtable_0.3.0       lifecycle_1.0.1    DBI_1.1.2         
## [41] magrittr_2.0.3     scales_1.2.0       writexl_1.4.0      cli_3.3.0         
## [45] stringi_1.7.6      farver_2.1.0       fs_1.5.2           xml2_1.3.3        
## [49] ellipsis_0.3.2     generics_0.1.2     vctrs_0.4.1        expint_0.1-7      
## [53] RColorBrewer_1.1-3 tools_4.2.2        glue_1.6.2         hms_1.1.1         
## [57] yaml_2.3.5         fastmap_1.1.0      colorspace_2.0-3   rvest_1.0.2       
## [61] knitr_1.39         haven_2.5.0
    Sys.time()
## [1] "2023-12-27 10:04:25 PST"