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

---
title: '03_Wastewater Treatment Plants'
subtitle: 'NPDES Wastload Allocations'
author: "Mercury Program and Basin Planning Unit"
date: "2/1/2023"
output:
  html_document:
    code_folding: show
    toc: TRUE
    toc_float: TRUE
    toc_depth: 3
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)
```


## Load Packages and Functions
```{r}
library(scales) # used for percent()
library(kableExtra) # better formatting of tables


AqImplGoal <- 0.059

# Issue with setting 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"))


# Read Data set
npdesStatus <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Source Data/03a_Municipal WWTPs (NPDES)/Data/03_NPDES Facility Info.xlsx"), guess_max = 300000, sheet = "Final Table") %>% 
  distinct(Discharger, Facility, Status)

wwtpSummary <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Source Data/03a_Municipal WWTPs (NPDES)/Output/03a_WWTP MeHg Annual Loads_2023-04-26.xlsx"), guess_max = 300000, sheet = "WWTP") %>% 
  filter(!is.na(Subarea)) %>%  # removes summary statistic rows
  left_join(npdesStatus, by=c("Discharger", "Facility")) %>%  # Adds status column
  mutate(effective_MeHg_median = if_else(is.na(final_MeHg_median), MeHg_median, final_MeHg_median),  # Creates col with effective median meHg conc (needed for facilities that have more than 1 discharge location)
         uses_ND_median = if_else(is.na(uses_ND_median), if_else(effective_MeHg_median < ND_median, 1, 0), uses_ND_median), # Adds whether ND value is used for facilities that have more than 1 discharge location
         effective_Flow_median = MedianLoad_grams * 1/365 * 1/effective_MeHg_median * 1000000000 * 0.2642 * 1/1000000) %>%  # Calculates effective median flow (MGD) based on calculated median load (needed for facilities that have more than 1 discharge location)
  arrange(Subarea, Discharger)

# Get column class from wwtpSummary
wwtpSummary_colTypes <- paste(map_chr(sapply(wwtpSummary, class), ~str_sub(., 1,1)), collapse = "")


indSummary1 <- readxl::read_xlsx(paste0(wd, "/Reeval_Source_Analysis/Source Data/03a_Municipal WWTPs (NPDES)/Output/03a_WWTP MeHg Annual Loads_2023-04-26.xlsx"), guess_max = 300000, sheet = "Ind") %>% 
  filter(!is.na(Subarea)) %>%  # removes summary statistic rows
  left_join(npdesStatus, by=c("Discharger", "Facility")) %>%  # Adds status column
  mutate(effective_MeHg_median = if_else(is.na(final_MeHg_median), MeHg_median, final_MeHg_median),  # Creates col with effective median meHg conc (needed for facilities that have more than 1 discharge location)
         uses_ND_median = if_else(is.na(uses_ND_median), if_else(effective_MeHg_median < ND_median, 1, 0), uses_ND_median), # Adds whether ND value is used for facilities that have more than 1 discharge location
         effective_Flow_median = MedianLoad_grams * 1/365 * 1/effective_MeHg_median * 1000000000 * 0.2642 * 1/1000000) %>%  # Calcualtes effective median flow (MGD) based on calculated median load (needed for facilities that have more than 1 discharge location)
  arrange(Subarea, Discharger)

# Asssign columns in indSummary the same class as wwtpSummary (otherwise columns with NA in indSummary won't match corresponding col in wwtpSummary and bind_rows will give an error)
indSummary <- type_convert(indSummary1, wwtpSummary_colTypes, guess_integer=T)


npdesSummary <- wwtpSummary %>%
  bind_rows(indSummary)

subareas <- npdesSummary %>% distinct(Subarea) %>% pull(Subarea)
# [1] "Central Delta"     "Marsh Creek"       "Sacramento"        "San Joaquin"       "West Delta"        "Yolo Bypass North"

# Taken from "04_Percent Allocations/Tables 8.4a-h MeHg Inputs Per Subarea.xlsx"
percentAllocations <- c("Central Delta" = .2439, "Marsh Creek" = .2181, "Sacramento" = .6392, "San Joaquin" = .3806, "West Delta" = .6935, "Yolo Bypass North" = .0679)

names(percentAllocations) <- subareas  # added this because names(percentAllocations) %in% subareas is False for Central Delta, reason unknown. True for all other subareas though
```

## Subarea Unassigned Future Growth Allocations (Recreating Table 8.3a)
```{r}
unassignedNewWWTP <- wwtpSummary %>% 
  group_by(Subarea) %>% 
  summarize(`Existing Eff Vol (MGD)` = sum(effective_Flow_median, na.rm=T),
            `Projected Vol for New WWTPS (MGD)` = `Existing Eff Vol (MGD)` * .25, # 25% projected population growth from 2020 to 2060
            `Unassigned MeHg Load for New WWTPS (g/yr)`  = (`Projected Vol for New WWTPS (MGD)` * 1/0.2642 * 1000000 * 365) * (AqImplGoal * 1/1000000000) # (MGD * 1 L/0.2642 gal * 1,000,000 gal/1 Mgal * 365d/yr) * (ng/L * g/1,000,000,000 ng)
            ) %>% 
  add_row(Subarea = "Total", summarize(., across(where(is.numeric), ~sum(.x, na.rm=T))))

unassignedNewWWTP %>% 
  kbl(align='lcc', caption="Existing Municipal WWTP Volume Discharged to Each Subarea & Projected Increase") %>%
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```

## Facilities with MeHg Conc. Less Than Aq Impl. Goal (Recreating Table 8.3b)
```{r}
dilutionWWTP <- wwtpSummary %>%
  select(Subarea, Discharger, Facility, Status, effective_MeHg_median, uses_ND_median, effective_Flow_median, MedianLoad_grams) %>% 
  filter(effective_MeHg_median <= AqImplGoal) %>% 
  mutate(`Projected Effluent Volume for Future Growth (MGD)` = if_else(Status != "Closing", ifelse(grepl("lathrop", Discharger, ignore.case=T), 6.0-effective_Flow_median, effective_Flow_median * 0.5), NA_real_), # accounting for a 50% Flow increase & giving lathrop allocation up to 6 MGD from Anti-Deg analysis
         `MeHg Allocation for Future Growth (g/yr)` = (`Projected Effluent Volume for Future Growth (MGD)` * 1/0.2642 * 1000000 * 365) * (effective_MeHg_median * 1/1000000000) # (MGD * 1 L/0.2642 gal * 1,000,000 gal/1 Mgal * 365d/yr) * (ng/L * g/1,000,000,000 ng)
  )

dilutionWWTP %>% 
  kbl(align='lcc', caption="Existing Municipal WWTP MeHg Conc. Less Than Aq Implementation Goal") %>%
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```


## WWTP Allocations by Subarea (WWTPs in Table 8.4)
```{r}
npdesLoads <- npdesSummary %>%
  select(Subarea, Discharger, Facility, Status, FacilityType, effective_MeHg_median, effective_Flow_median, MedianLoad_grams)


npdesAllocations <- list()
for(subarea in subareas){
  allowWWTPgrowth <- unassignedNewWWTP %>% 
    filter(grepl(subarea, Subarea)) %>%
    mutate(Discharger = "Unassigned NPDES Facility Allocation",
           Facility = "(allowable future growth)",
           `MeHg Conc Used to Calc Allocation` = AqImplGoal) %>%
    rename(`MGD Used to Calc Allocation` = `Projected Vol for New WWTPS (MGD)`,
           `MeHg Waste Load Allocation` = `Unassigned MeHg Load for New WWTPS (g/yr)`) %>% 
    select(Subarea, Discharger, Facility, `MeHg Conc Used to Calc Allocation`, `MGD Used to Calc Allocation`, `MeHg Waste Load Allocation`)
  
  
  npdesSubarea <- npdesLoads %>%
    filter(grepl(subarea, Subarea)) %>% 
    mutate(
      `Percent Allocation` = ifelse(effective_MeHg_median > AqImplGoal, ifelse(AqImplGoal/effective_MeHg_median > percentAllocations[subarea], AqImplGoal/effective_MeHg_median, percentAllocations[subarea]), 1),
      `MeHg Conc Used to Calc Allocation` = effective_MeHg_median * `Percent Allocation`,
      `MGD Used to Calc Allocation`       = effective_Flow_median,
      `MeHg Waste Load Allocation`        = (`MeHg Conc Used to Calc Allocation` * 1/1000000000) * (`MGD Used to Calc Allocation` * 1/0.2642 * 1000000 * 365),  # (ng/L * g/1,000,000,000 ng) * (MGD * 1 L/0.2642 gal * 1,000,000 gal/1 Mgal * 365d/yr)
           )

  dischargers <- paste(npdesSubarea$Discharger, npdesSubarea$Facility)
  for(discharger in dischargers){
    
    npdesDischarger <- npdesSubarea %>% 
      filter(grepl(discharger, paste(Discharger,Facility)))
    
    if(npdesDischarger$Status != "Closing" & npdesDischarger$effective_MeHg_median <= AqImplGoal){
      
      npdesDischargerGrowth <- npdesDischarger %>%
        mutate(FacilityType                  = "(allowable future growth)",
               `MGD Used to Calc Allocation` = ifelse(grepl("lathrop", Discharger, ignore.case=T), 6.0-effective_Flow_median, effective_Flow_median * 0.5), # Allow Lathrop to increase up 6 MGD based on their Antidegradation Analysis, otherwise allow 50% flow increase for WWTPs if MeHg conc is less than AqImplGoal
               `MeHg Waste Load Allocation`  = (`MeHg Conc Used to Calc Allocation` * 1/1000000000) * (`MGD Used to Calc Allocation` * 1/0.2642 * 1000000 * 365))  # (ng/L * g/1,000,000,000 ng) * (MGD * 1 L/0.2642 gal * 1,000,000 gal/1 Mgal * 365d/yr)
      
      npdesSubarea <- npdesSubarea %>% 
        bind_rows(npdesDischargerGrowth)
    }
  }

  npdesSubarea <- npdesSubarea %>%    
    arrange(Discharger, Facility, desc(FacilityType)) %>% 
    bind_rows(allowWWTPgrowth)
    
  
  subareaTotals <- npdesSubarea %>%
  add_row(Facility = "Total", summarize(., across(contains(c('load')), ~sum(.x, na.rm=T)))) %>%
  slice_tail(n=1) 

  
  npdesAllocations[[subarea]] <- npdesSubarea %>%
    bind_rows(., subareaTotals) %>%
    rename(`Existing Median MeHg Conc`  = effective_MeHg_median,
           `Existing Median Flow (MGD)` = effective_Flow_median,
           `Existing Annual MeHg Load` = MedianLoad_grams) %>%
    mutate(`Percent Allocation` = percent(`Percent Allocation`, accuracy=0.01))
  
}
```


## Export to excel
```{r}
writexl::write_xlsx(list(`Unassigned Subarea`=unassignedNewWWTP,
                         `Dilution WWTP`=dilutionWWTP),
                    paste0(wd, "/Reeval_Allocations/03_NPDES WWTP WLAs/Output/NPDES WWTP Unassgnd & Future Growth WLAs_", today(), ".xlsx"))

writexl::write_xlsx(npdesAllocations,
                    paste0(wd, "/Reeval_Allocations/03_NPDES WWTP WLAs/Output/NPDES WWTP WLAs by Subarea_", today(), ".xlsx"))
```
## Error: <text>:19:1: unexpected '<'
## 18: 
## 19: <
##     ^

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] compiler_4.2.2 magrittr_2.0.3 tools_4.2.2    stringi_1.7.6  highr_0.9     
## [6] knitr_1.39     stringr_1.4.0  xfun_0.31      evaluate_0.15
    Sys.time()
## [1] "2023-12-26 13:35:29 PST"