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

---
title: 'NPDES Allocation Comparison'
author: "Mercury Program and Basin Planning Unit"
date: "2/15/2023"
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(kableExtra) # better formatting of tables
library(ggbeeswarm) # for geom_beeswarm in ggplot - stops over plotting
library(janitor) # for clean_names()
library(shiny)

# Had issue trying to set WD with Shiny in R project, reset working directory of rproj
wd <- rstudioapi::getActiveProject()
setwd(wd)

source("R Functions/functions_estimate NDDNQ values.R")
source("R Functions/functions_QA data.R")
source("R Functions/ggplot_theme_border.r")
source("R Functions/ggplot_force_panelsizes.r")

file <- function(file_name){
  paste0(wd, "/Reeval_Allocations/04_Percent Allocations/", file_name)
}
```

## Load Data
```{r, eval=T, echo=T}
# Change chunck option to "eval=TRUE" if changes were made to "Tables 8.4a-h MeHg Inputs Per Subarea.xlsx"
# e.g. Mean concentrations were added to WWTP or MS4 facilities

# Idea is to read "Tables 8.4a-h MeHg Inputs Per Subarea.xlsx" and concatenate sheets into one dataframe
  # how to: https://readxl.tidyverse.org/articles/readxl-workflows.html
  # Will likely need to clean up spreadsheet and condense multiple header rows into one row & Add "NA" for missing values
  # Code will need to specify range to read for each sheet
sheets <- file("Tables 8.4a-h MeHg Inputs Per Subarea.xlsx") %>%
  excel_sheets() %>%
  set_names()

sheets <- grep("summary|north|south", sheets, ignore.case=T, value=T, invert=T) # do not select summary sheets or Yolo ByPass North & South


# specify Cell Ranges for sheets
#         Central Delta, Marsh Creek, MokelumneCosumnes Rivers, Sacramento River, San Joaquin River, West Delta, & Yolo Bypass - Combined
ranges <- list("A1:O20",    "A1:O10",                 "A1:O10",         "A1:O20",          "A1:O29",   "A1:O12",                "A1:O24")


# Concatenate sheets
allocations_npdes <- map2_df(
  sheets,
  ranges,
  ~ read_xlsx(file("Tables 8.4a-h MeHg Inputs Per Subarea.xlsx"), sheet = .x, range = .y, na=c("NA", "depends on tributary", "depends on discharge location", "depends on location", "Not in 2010 TMDL Staff Report")),
  .id = "Subarea"
) %>%
  clean_names() %>%
  filter(grepl("wwtp|ms4", me_hg_source, ignore.case=T)) %>%
  write_csv(file("output/WWTP & MS4 Allocation Complilation.csv"))
```

### Add Future Growth Allocations to assigned Facilities 
```{r, echo=T}
# Load Formatted Data
allocations_npdes <- read_csv(file("output/WWTP & MS4 Allocation Complilation.csv")) %>% 
  mutate(subarea_graph = case_when(grepl("Central", subarea) ~ "Central\nDelta",
                             grepl("Marsh", subarea) ~ "Marsh\nCreek",
                             grepl("Moke", subarea) ~ "Moke/Cos\nRivers",
                             grepl("Sac", subarea)  ~ "Sacramento\nRiver",
                             grepl("San", subarea)  ~ "San Joaquin\nRiver",
                             grepl("West", subarea) ~ "West\nDelta",
                             grepl("Yolo", subarea) ~ "YB\nCombined",
                             T ~ subarea),
         subarea_table = case_when(grepl("Central", subarea) ~ "Central Delta",
                             grepl("Marsh", subarea) ~ "Marsh Creek",
                             grepl("Moke", subarea) ~ "Moke/Cos Rivers",
                             grepl("Sac", subarea)  ~ "Sacramento River",
                             grepl("San", subarea)  ~ "San Joaquin River",
                             grepl("West", subarea) ~ "West Delta",
                             grepl("Yolo", subarea) ~ "YB Combined",
                             T ~ subarea)) %>% 
  arrange(me_hg_source, subarea, specific_source)

allocations_npdes %>% 
  kbl(align='lcccccl', caption='WWTP & MS4 Individual WLAs') %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)

allocations_facility_sum <- allocations_npdes %>% 
  mutate(specific_source = sub(" \\(allowable future growth)", "" , specific_source)) %>% # remove " (allowable future growth)" to add facility loads together
  group_by(me_hg_source, subarea_graph, subarea_table, specific_source) %>% 
  summarise(across(ends_with("_g_yr"), ~sum(.x, na.rm=T)),
            .groups="drop") %>% 
  replace (.==0, NA_real_) %>%  # Replace 0's with NA so 0 values aren't graphed
  mutate(dmcp_review_compliance = ifelse(dmcp_review_median_annual_me_hg_load_g_yr <= dmcp_review_me_hg_allocation_g_yr, 1, 0),
         tmdl_2010_compliance = ifelse(tmdl_2010_staff_report_average_annual_me_hg_load_g_yr <= tmdl_2010_staff_report_me_hg_allocation_g_yr, 1, 0),
         dmcp_review_average_2010reduction = ifelse(!is.na(tmdl_2010_staff_report_average_annual_me_hg_load_g_yr),
                                                    ifelse(dmcp_review_average_annual_me_hg_load_g_yr <= tmdl_2010_staff_report_average_annual_me_hg_load_g_yr, 1, 0),
                                                    NA_real_),
         dmcp_review_average_2010compliance = ifelse(!is.na(tmdl_2010_staff_report_me_hg_allocation_g_yr),
                                                     ifelse(dmcp_review_average_annual_me_hg_load_g_yr <= tmdl_2010_staff_report_me_hg_allocation_g_yr, 1, 0),
                                                     NA_real_),
         dmcp_review_average_2010either = ifelse(dmcp_review_average_2010reduction == 1 | dmcp_review_average_2010compliance == 1, 1, 0)) %>% 
  mutate(#WWTP Name Changes
         # specific_source=sub(" WWTP", "", specific_source),
         specific_source=sub("City of | Vocational Institution| CSD", "", specific_source),
         specific_source=sub("Regional County", "Reg Co", specific_source),
         specific_source=sub("Sacramento", "Sac", specific_source),
         #MS4 Name Changes
         specific_source=sub("San Joaquin", "SJ", specific_source),
         specific_source=sub("Community Services District", "CSD", specific_source),
         specific_source=sub("Sac Stormwater Quality Partnership", "SSQP", specific_source))

writexl::write_xlsx(allocations_facility_sum, path=file('output/allocations_facility_summary.xlsx'))

allocations_facility_sum %>% 
  kbl(align='lcccccl', caption='WWTP & MS4 Total WLAs (Summed with Future Growth)') %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```

### DMCP Review & 2010 Allocation comparision
```{r}
allocation_comparison <- allocations_facility_sum %>%
  select(specific_source, me_hg_source, subarea_table, tmdl_2010_staff_report_average_annual_me_hg_load_g_yr, tmdl_2010_staff_report_me_hg_allocation_g_yr, dmcp_review_median_annual_me_hg_load_g_yr, dmcp_review_me_hg_allocation_g_yr
         ) %>%
  
  rename(Facility = specific_source,
         `Facility Type` = me_hg_source,
         Subarea = subarea_table,
         `2010 TMDL Avg Annual Load (g/yr)` = tmdl_2010_staff_report_average_annual_me_hg_load_g_yr,
         `2010 TMDL Allocation (g/yr)` = tmdl_2010_staff_report_me_hg_allocation_g_yr,
         `DMCP Review Median Annual Load (g/yr)` = dmcp_review_median_annual_me_hg_load_g_yr,
         `DMCP Review Allocation (g/yr)` = dmcp_review_me_hg_allocation_g_yr
         ) %>%
  
  # filter(!is.na(`2010 TMDL Allocation (g/yr)`)) %>% 
  
  mutate(`2010 Load Reduction (g/yr)`= ifelse( `2010 TMDL Avg Annual Load (g/yr)` - `2010 TMDL Allocation (g/yr)` > 0, `2010 TMDL Avg Annual Load (g/yr)` - `2010 TMDL Allocation (g/yr)`, 0)) %>%
  
  mutate(`DMCP Review Load Reduction (g/yr)` = ifelse(`DMCP Review Median Annual Load (g/yr)` - `DMCP Review Allocation (g/yr)` > 0, `DMCP Review Median Annual Load (g/yr)` - `DMCP Review Allocation (g/yr)`, 0)) %>%
  
    # mutate(`DMCP Review Load Reduction (g/yr)` = ifelse(!is.na(`2010 TMDL Allocation (g/yr)`), ifelse(`DMCP Review Median Annual Load (g/yr)` - `DMCP Review Allocation (g/yr)` > 0, `DMCP Review Median Annual Load (g/yr)` - `DMCP Review Allocation (g/yr)`, 0), NA_real_)) %>%
  
  arrange(desc(`Facility Type`), `DMCP Review Load Reduction (g/yr)`, Facility) %>% 
  
  add_row(Facility = "Cumulative Total", summarize(., across(contains("(g/yr)"), ~sum(.x, na.rm=T)))) %>%
  
  mutate(`2010 TMDL % Reduction Needed of Load to Meet Allocation` = round((`2010 TMDL Avg Annual Load (g/yr)` - `2010 TMDL Allocation (g/yr)`)/`2010 TMDL Avg Annual Load (g/yr)`, 2),
         `DMCP Review % Reduction Needed of Load to Meet Allocation` = round((`DMCP Review Median Annual Load (g/yr)` - `DMCP Review Allocation (g/yr)`)/`DMCP Review Median Annual Load (g/yr)`, 2))
    

writexl::write_xlsx(allocation_comparison, path=file('output/allocation_comparison.xlsx'))

allocation_comparison %>% 
  kbl(align='lcccccl', caption='WWTP & MS4 Total WLAs (Summed with Future Growth)') %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)

 
```


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

used_colors <- c("DMCP Review\nProposed WLA" = colors_graphs[9], "DMCP Review\nMedian Load"=colors_graphs[3],
                 "2010 TMDL\nWLA"=colors_graphs[9], "2010 TMDL\nAvg Load"=colors_graphs[11], "DMCP Review\nAvg Load"=colors_graphs[3])

base_theme <- theme_set(theme_light() +
                          theme(text = element_text(size=14, family="sans"),
                                axis.title=element_text(size=12, face="bold"),
                                axis.title.x = element_blank(),
                                axis.text.x = element_text(size=10.5, angle = 45, hjust = 1),
                                panel.grid.major.y=element_line(size=1),
                                panel.grid.minor=element_line(size=0.5),
                                legend.key.height=unit(1.5, "cm"),
                                legend.text = element_text(size = 14),
                                legend.title = element_blank()
                                )
                        )
```

## Idea to Replace Long labels with Alphabetic
```{r}
# https://stackoverflow.com/questions/30614105/create-abbreviated-legends-manually-for-long-x-labels-in-ggplot2
```



## DMCP Review Load & Allocation Comparison
```{r}
dmcp_allocations <- allocations_facility_sum %>% 
  select(me_hg_source:dmcp_review_me_hg_allocation_g_yr) %>%
  rename(`DMCP Review\nMedian Load` = dmcp_review_median_annual_me_hg_load_g_yr, 
         `DMCP Review\nProposed WLA` = dmcp_review_me_hg_allocation_g_yr) %>% 
  pivot_longer(cols=c("DMCP Review\nMedian Load", "DMCP Review\nProposed WLA"),
               names_to="DMCP Review",
               values_to="Methylmercury Load (g/yr)")
```

### WWTP Compliance
```{r}
allocation_wwtp <- dmcp_allocations %>% 
  filter(grepl("wwtp", me_hg_source, ignore.case=T))

allocation_wwtp %>% 
  kbl(align='lcccccl', caption='WWTP Total WLAs (Summed with Future Growth)') %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)

# wwtp_wla_plot <- allocation_wwtp %>% 
#   ggplot() +
#   scale_fill_manual(values = used_colors[unique(dmcp_allocations$`DMCP Review`)])
 
zoom_value <- 0.65 

# Graph All 
wwtp_wla_all <- allocation_wwtp %>% 
  ggplot() +
  scale_fill_manual(values = used_colors[unique(dmcp_allocations$`DMCP Review`)]) +
  # Gray Rectangle Annotation to highlight Zoom Region 
  annotate(geom="rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=zoom_value, fill=alpha("grey", 0.3)) +
  
  # Full Scale of all data points
  geom_point(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=`DMCP Review`),
             size=3, colour="black", pch=21) +
  
  # Axis padding 
  scale_y_continuous(expand = expansion(mult = c(0.01, 0.03))) + 
  facet_grid(~ subarea_graph, scales="free", space="free") +
  
  # add minor grid lines and adjust space between plots
  theme(axis.text.x = element_text(angle = 59),
        panel.border = theme_border(type=c("right", "left"), colour="grey52", size = 0.6),
        panel.spacing.x = unit(0, "lines"))


# Graph Zoom of WWTPs less than < 1 g/yr
wwtp_wla_zoom <- allocation_wwtp %>% 
  ggplot() +
  scale_fill_manual(values = used_colors[unique(dmcp_allocations$`DMCP Review`)]) +
  # Gray Rectangle Annotation to highlight Zoom Region 
  annotate(geom="rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=zoom_value, fill=alpha("grey88", 0.3)) +
  
  # Zoom in on data points less than 1 g/yr
  # geom_point(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=`DMCP Review`),
  #            size=3, colour="black", pch=21) +
  geom_beeswarm(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=`DMCP Review`),
                size=3, colour="black", pch=21,
                cex=3, method="swarm") +
  
  # Axis padding 
  scale_y_continuous(limits=c(0,zoom_value), expand = expansion(mult = c(0.01, 0))) +
  facet_grid(~ subarea_graph, scales="free", space="free") +
  
  # add minor grid lines and add space between plots
  theme(axis.text.x = element_text(angle = 59),
        panel.border = theme_border(type=c("right", "left"), colour="grey52", size = 0.6),
        panel.spacing.x = unit(0, "lines"))



# View Plots
wwtp_wla_all
wwtp_wla_zoom

# Save Graphs
ggsave(filename = file("output/Fig 8.1a - DMCP Review Load & Allocation Comparison_WWTP_full.png"),
       plot = wwtp_wla_all,
       width = 9,
       height = 5.75,
       units = "in",
       dpi = 300)

ggsave(filename = file("output/Fig 8.1b - DMCP Review Load & Allocation Comparison_WWTP.zoom.png"),
       plot = wwtp_wla_zoom,
       width = 9,
       height = 5.75,
       units = "in",
       dpi = 300)
```


### MS4 Compliance Graph
```{r}
allocation_ms4 <- dmcp_allocations %>% 
  filter(grepl("ms4", me_hg_source, ignore.case=T))

allocation_ms4 %>% 
  kbl(align='lcccccl', caption='WWTP Total WLAs (Summed with Future Growth)') %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)


allocation_ms4 %>% 
  ggplot() +
  # geom_point(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=`DMCP Review`),
  #            size=3, colour="black", pch=21) +
  geom_beeswarm(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=`DMCP Review`),
                size=3, colour="black", pch=21,
                cex=3.5, method="swarm") +
  scale_fill_manual(values = used_colors[unique(dmcp_allocations$`DMCP Review`)]) +
  scale_y_continuous(expand = expansion(mult = c(0.01, 0.025))) + 
  facet_grid(. ~ subarea_graph, scales = "free", space = "free") +
  force_panelsizes(cols = c(.45, .25, .45, .5, 0.9, .35, .5)) +
  theme(axis.text.x = element_text(angle = 52),
        panel.border = theme_border(type=c("right", "left"), colour="grey52", size = 0.6),
        panel.spacing.x = unit(0, "lines"))

ggsave(filename = file("output/Fig 8.3 - DMCP Review Load & Allocation Comparison_MS4.png"),
       # plot = predictionGraph,
       width = 9,
       height = 5.75,
       units = "in",
       dpi = 300)
```



## Median Load & Allocation Compared to Mean Load & Allocation in 2010 TMDL Staff Report
```{r}
dmcp_2010_comparison <- allocations_facility_sum %>% 
  filter(!is.na(tmdl_2010_staff_report_me_hg_allocation_g_yr)) %>% # filter for facilities that were in 2010 TMDL
  select(me_hg_source:specific_source, dmcp_review_median_annual_me_hg_load_g_yr, dmcp_review_me_hg_allocation_g_yr,
         tmdl_2010_staff_report_average_annual_me_hg_load_g_yr, tmdl_2010_staff_report_me_hg_allocation_g_yr) %>%
  rename(`DMCP Review\nMedian Load` = dmcp_review_median_annual_me_hg_load_g_yr,
         `DMCP Review\nWLA` = dmcp_review_me_hg_allocation_g_yr,
         `2010 TMDL\nAvg Load` = tmdl_2010_staff_report_average_annual_me_hg_load_g_yr,
         `2010 TMDL\nWLA` = tmdl_2010_staff_report_me_hg_allocation_g_yr) %>% 
  pivot_longer(cols=c("DMCP Review\nMedian Load", "DMCP Review\nWLA", "2010 TMDL\nAvg Load", "2010 TMDL\nWLA"),
               names_to="Comparison",
               values_to="Methylmercury Load (g/yr)")
```

### WWTP Comparison
```{r}
compare_wwtp <- dmcp_2010_comparison %>% 
  filter(grepl("wwtp", me_hg_source, ignore.case=T))

compare_wwtp_plot <- compare_wwtp %>% 
  ggplot() +
  scale_fill_manual(values = used_colors[unique(compare_wwtp$`Comparison`)])

zoom_value <- 2.7

# Graph All 
compare_wwtp_all <- compare_wwtp_plot + 
  # Gray Rectangle Annotation to highlight Zoom Region 
  annotate(geom="rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=zoom_value, fill=alpha("grey", 0.3)) +
  
  # Full Scale of all data points
  geom_point(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=`Comparison`),
             size=3, colour="black", pch=21) +
  
  # Axis padding 
  scale_y_continuous(expand = expansion(mult = c(0.01, 0.03))) + 
  facet_grid(~ subarea_graph, scales="free", space="free") +
  
  # add minor grid lines and adjust space between plots
  theme(axis.text.x = element_text(angle = 56),
        panel.border = theme_border(type=c("right", "left"), colour="grey52", size = 0.6),
        panel.spacing.x = unit(0, "lines"))


# Graph Zoom of WWTPs less than < 1 g/yr
compare_wwtp_zoom <- compare_wwtp_plot + 
  # Gray Rectangle Annotation to highlight Zoom Region 
  annotate(geom="rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=zoom_value, fill=alpha("grey88", 0.3)) +
  
  # Zoom in on data points less than 1 g/yr
  # geom_point(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=`Comparison`),
  #            size=3, colour="black", pch=21) +
  geom_beeswarm(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=Comparison),
                size=3, colour="black", pch=21,
                cex=3.4, method="swarm") +
  
  # Axis padding 
  scale_y_continuous(limits=c(0,zoom_value), expand = expansion(mult = c(0.01, 0))) +
  facet_grid(~ subarea_graph, scales="free", space="free") +
  
  # add minor grid lines and add space between plots
  theme(axis.text.x = element_text(angle = 56),
        axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0)), # adds space after y-axis title so plot gridlines line up with compare_wwtp_all
        panel.border = theme_border(type=c("right", "left"), colour="grey52", size = 0.6),
        panel.spacing.x = unit(0, "lines"))

# View Plots
compare_wwtp_all
compare_wwtp_zoom

# Save Graphs
ggsave(filename = file("output/Fig 8.Xa - Compare DMCP Review Median to 2010 TMDL_WWTP_full.png"),
       plot = compare_wwtp_all,
       width = 9,
       height = 5.75,
       units = "in",
       dpi = 300)

ggsave(filename = file("output/Fig 8.Xb - Compare DMCP Review Median to 2010 TMDL_WWTP.zoom.png"),
       plot = compare_wwtp_zoom,
       width = 9,
       height = 5.75,
       units = "in",
       dpi = 300)
```


### MS4 Comparison
```{r}
compare_ms4 <- dmcp_2010_comparison %>% 
  filter(grepl("ms4", me_hg_source, ignore.case=T))

# Look into geom_dotplot if overplotting is an issue
compare_ms4 %>% 
  ggplot() +
  # geom_point(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=Comparison),
  #            size=3, colour="black", pch=21) +
  geom_beeswarm(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=Comparison),
                size=3, colour="black", pch=21,
                cex=3.4, method="swarm") +
  
  scale_fill_manual(values = used_colors[unique(dmcp_2010_comparison$Comparison)]) +
  scale_y_continuous(expand = expansion(mult = c(0.01, 0.025))) + 
  facet_grid(. ~ subarea_graph, scales = "free", space = "free") +
  force_panelsizes(cols = c(.5, .2, .35, .6, .25, .25)) +
  theme(axis.text.x = element_text(angle = 61),
        panel.border = theme_border(type=c("right", "left"), colour="grey52", size=0.6),
        panel.spacing.x = unit(0, "lines"))

ggsave(filename = file("output/Fig 8.X - Compare DMCP Review Median to 2010 TMDL_MS4.png"),
       # plot = predictionGraph,
       width = 9,
       height = 5.75,
       units = "in",
       dpi = 300)
```



## Mean Load Compared to Mean Load & Allocation in 2010 TMDL Staff Report
```{r}
dmcp_2010_comparison <- allocations_facility_sum %>% 
  filter(!is.na(tmdl_2010_staff_report_me_hg_allocation_g_yr)) %>% # filter out facilities that were in 2010 TMDL
  select(me_hg_source:specific_source, dmcp_review_average_annual_me_hg_load_g_yr:tmdl_2010_staff_report_me_hg_allocation_g_yr) %>%
  rename(`DMCP Review\nAvg Load` = dmcp_review_average_annual_me_hg_load_g_yr,
         `2010 TMDL\nAvg Load` = tmdl_2010_staff_report_average_annual_me_hg_load_g_yr,
         `2010 TMDL\nWLA` = tmdl_2010_staff_report_me_hg_allocation_g_yr) %>% 
  pivot_longer(cols=c("DMCP Review\nAvg Load", "2010 TMDL\nAvg Load", "2010 TMDL\nWLA"),
               names_to="Comparison",
               values_to="Methylmercury Load (g/yr)")
```

### WWTP Comparison
```{r}
compare_wwtp <- dmcp_2010_comparison %>% 
  filter(grepl("wwtp", me_hg_source, ignore.case=T))

compare_wwtp_plot <- compare_wwtp %>% 
  ggplot() +
  scale_fill_manual(values = used_colors[unique(compare_wwtp$`Comparison`)])

zoom_value <- 2.7

# Graph All 
compare_wwtp_all <- compare_wwtp_plot + 
  # Gray Rectangle Annotation to highlight Zoom Region 
  annotate(geom="rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=zoom_value, fill=alpha("grey", 0.3)) +
  
  # Full Scale of all data points
  geom_point(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=`Comparison`),
             size=3, colour="black", pch=21) +
  
  # Axis padding 
  scale_y_continuous(expand = expansion(mult = c(0.01, 0.03))) + 
  facet_grid(~ subarea_graph, scales="free", space="free") +
  
  # add minor grid lines and adjust space between plots
  theme(axis.text.x = element_text(angle = 56),
        panel.border = theme_border(type=c("right", "left"), colour="grey52", size = 0.6),
        panel.spacing.x = unit(0, "lines"))


# Graph Zoom of WWTPs less than < 1 g/yr
compare_wwtp_zoom <- compare_wwtp_plot + 
  # Gray Rectangle Annotation to highlight Zoom Region 
  annotate(geom="rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=zoom_value, fill=alpha("grey88", 0.3)) +
  
  # Zoom in on data points less than 1 g/yr
  # geom_point(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=`Comparison`),
  #            size=3, colour="black", pch=21) +
  geom_beeswarm(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=Comparison),
                size=3, colour="black", pch=21,
                cex=3.4, method="swarm") +
  
  # Axis padding 
  scale_y_continuous(limits=c(0,zoom_value), expand = expansion(mult = c(0.01, 0))) +
  facet_grid(~ subarea_graph, scales="free", space="free") +
  
  # add minor grid lines and add space between plots
  theme(axis.text.x = element_text(angle = 56),
        axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0)), # adds space after y-axis title so plot gridlines line up with compare_wwtp_all
        panel.border = theme_border(type=c("right", "left"), colour="grey52", size = 0.6),
        panel.spacing.x = unit(0, "lines"))

# View Plots
compare_wwtp_all
compare_wwtp_zoom

# Save Graphs
ggsave(filename = file("output/Fig 8.2a - Compare DMCP Review Avg to 2010 TMDL_WWTP_full.png"),
       plot = compare_wwtp_all,
       width = 9,
       height = 5.75,
       units = "in",
       dpi = 300)

ggsave(filename = file("output/Fig 8.2b - Compare DMCP Review Avg to 2010 TMDL_WWTP.zoom.png"),
       plot = compare_wwtp_zoom,
       width = 9,
       height = 5.75,
       units = "in",
       dpi = 300)
```


### MS4 Comparison
```{r}
compare_ms4 <- dmcp_2010_comparison %>% 
  filter(grepl("ms4", me_hg_source, ignore.case=T))

# Look into geom_dotplot if overplotting is an issue
compare_ms4 %>% 
  ggplot() +
  # geom_point(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=Comparison),
  #            size=3, colour="black", pch=21) +
  geom_beeswarm(aes(x=specific_source, y=`Methylmercury Load (g/yr)`, fill=Comparison),
                size=3, colour="black", pch=21,
                cex=3.4, method="swarm") +
  
  scale_fill_manual(values = used_colors[unique(dmcp_2010_comparison$Comparison)]) +
  scale_y_continuous(expand = expansion(mult = c(0.01, 0.025))) + 
  facet_grid(. ~ subarea_graph, scales = "free", space = "free") +
  force_panelsizes(cols = c(.5, .2, .35, .6, .25, .25)) +
  theme(axis.text.x = element_text(angle = 61),
        panel.border = theme_border(type=c("right", "left"), colour="grey52", size=0.6),
        panel.spacing.x = unit(0, "lines"))

ggsave(filename = file("output/Fig 8.4 - Compare DMCP Review Avg to 2010 TMDL_MS4.png"),
       # plot = predictionGraph,
       width = 9,
       height = 5.75,
       units = "in",
       dpi = 300)
```
## Error: <text>:17:1: unexpected '<'
## 16: ---
## 17: <
##     ^

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 14:14:11 PST"