# Date(s): 2/3/2022
# Project: Carli Creek Water Quality Assessment Project
# Created For: Pollutants, medians, OWQC, and Intern'l BMP database visulization. Also for visulaizing pollutant reductions compared to rainfall.
# Name: Chris Desiderati
###
setwd("C:/Users/irish/Documents/ACADEMIC/GRAD SCHOOL/Project_Carli/R")
getwd()
library(tidyverse) #activate tidyverse
library(lubridate) #activate lubridate to deal with time 
library(dplyr) #to transform data
library(reshape2)
library(cowplot)
library(egg)
#### Import Data ####
dta<-read.csv("Ecoli.csv") #import E. coli
class(dta$Sample.Pt.Descr)
dta<-rename(dta,Date=Collection.Date)#rename variable
dta$Date<-mdy(dta$Date) # convert factor to Date
class(dta$Date)
str(dta)
brks <- as.Date(c("2020-10-01","2020-12-01","2021-02-01","2021-04-01"))
#### E. coli #### 
# OWQC (site 5 only) + IBMPdb (bioretention only) for just E.coli
dta<-dta %>%
  mutate(IBMPdb_BR=ifelse(Sample.Pt.Descr=="3_PostTerr",158,NA),
         OWQC=ifelse(Sample.Pt.Descr=="5_DS"|Sample.Pt.Descr=="4_US",406,NA))#Create OWQC and IBMPdb BR median
dta<-dta %>% # Compute median by Monitoring Points
  group_by(Sample.Pt.Descr) %>%
  mutate(Medians=median(COLILERT))
dta1<-dta %>%
  filter(Sample.Pt.Descr=="3_PostTerr")# create special df for site 3 ribbon
g<-ggplot(data=dta,aes(x=Date))+
  geom_point(aes(y=COLILERT))+
  geom_line(aes(y=COLILERT))+
  geom_line(aes(y=Medians,linetype="Median"),color="steelblue",size=1)+
  geom_line(aes(y=OWQC,colour="OWQC"))+ #Oregon single-sample WQC
  geom_line(aes(y=IBMPdb_BR,colour="Intl BMP BR"))+ #Bioretention median
  geom_ribbon(data=dta1,aes(ymin=46.5,ymax=212),fill="green",alpha=0.1)+ #BR CFs
  theme_minimal()+
  scale_x_date(limits = as.Date(c('2020-10-05','2021-04-30')),
               date_labels="%b '%y",
               breaks=brks
  )+
  theme(
    axis.text.x = element_text(angle = 45, hjust=1),
    axis.title.x=element_blank(),
    legend.position="bottom",
    legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  labs(colour="Criterion",
       fill="Confidence Intervals",
       linetype="Statistic")+
  scale_linetype_manual(values = "dashed")+
  scale_color_manual(values=c("OWQC"="red",'Intl BMP BR'="green"))+
  ylab("E. coli, MPN/100 mL")+
  scale_y_log10()
g+facet_grid(cols=vars(Sample.Pt.Descr)) #Looks good
#### Nutrients ####
# IBMP (bioretention only) for Ammonia, Nitrate-Nitrite, and TP
library(reshape2) #to deal with reshaping data from 'wide' to 'long'
dtan<-read.csv("AmmNO3NO2TPIBMPdb.csv")
dtan<-rename(dtan,Date=Collection.Date)#rename variable
dtan$Date<-mdy(dtan$Date)
class(dtan$Date)

mdtan<-melt(dtan,id.vars=c("Event","Sample.Pt.Descr","Date"),measure.vars=c('Ammonia','Nitrate.Nitrite','Total.Phosphorous'))
str(mdtan)
levels(mdtan$variable)[levels(mdtan$variable)=="Nitrate.Nitrite"]<-"Nitrate-Nitrite"
levels(mdtan$variable)[levels(mdtan$variable)=="Total.Phosphorous"]<-"Total Phosphorous"
str(mdtan)
mdtan<-mdtan %>%
  mutate(NH3_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&variable=="Ammonia",0.05,NA),#IBMPdb NH3 medians
         NO3NO2_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&variable=="Nitrate-Nitrite",0.441,NA),#IBMPdb NO3NO2 medians
         TPO4_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&variable=="Total Phosphorous",0.240,NA)) #IBMPdb TPhos
e.s1<-mdtan %>%# Create special ribbon dataframes
  filter(Sample.Pt.Descr=="3_PostTerr",variable=='Ammonia')
e.s2<-mdtan %>%
  filter(Sample.Pt.Descr=="3_PostTerr",variable=='Nitrate-Nitrite')
e.s3<-mdtan %>%
  filter(Sample.Pt.Descr=="3_PostTerr",variable=='Total Phosphorous')
mdtan<-mdtan %>% # Compute median by Monitoring Points and variable
  group_by(Sample.Pt.Descr,variable)%>%
  mutate(Medians=median(value))
mdtan1<-mdtan%>%
  filter(variable=="Ammonia")
mdtan2<-mdtan%>%
  filter(variable=="Nitrate-Nitrite")
mdtan3<-mdtan%>%
  filter(variable=="Total Phosphorous")
f1<-ggplot(data = mdtan1,aes(Date))+
  geom_point(aes(y=value))+
  geom_line(aes(y=value))+
  geom_line(aes(y=Medians,linetype="Medians"),color="steelblue", size=1)+
  geom_line(aes(y=NH3_BR,color="BR"))+ #Bioretention NH3 median 
  geom_ribbon(data=e.s1,aes(ymin=0.05,ymax=0.06),fill='green',alpha=0.2)+
  scale_linetype_manual(values=2)+
  scale_color_manual(values="Green")+
  scale_x_date(limits = as.Date(c('2020-10-05','2021-04-30')),
               date_labels="%b '%y",
               breaks=brks
  )+
  scale_y_continuous(limits=c(0,2),breaks=c(0,0.5,1,1.5,2),sec.axis=sec_axis(~.,name="Ammonia",breaks=NULL,labels=NULL))+
  theme_minimal()+
  theme(axis.text.x=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y.left=element_blank(),
        legend.position="none",
        )
f1.<-f1+facet_grid(cols=vars(Sample.Pt.Descr))
f1.
f2<-ggplot(data = mdtan2,aes(Date))+
  geom_point(aes(y=value))+
  geom_line(aes(y=value))+
  geom_line(aes(y=Medians,linetype="Medians"),color="steelblue", size=1)+
  geom_line(aes(y=NO3NO2_BR,color="BR"))+ #NO3NO2
  geom_ribbon(data=e.s2,aes(ymin=0.380,ymax=0.507),fill='green',alpha=0.2)+
  scale_linetype_manual(values=2)+
  scale_color_manual(values="Green")+
  scale_x_date(limits = as.Date(c('2020-10-05','2021-04-30')),
               date_labels="%b '%y",
               breaks=brks
  )+
  scale_y_continuous(limits=c(0,2.2),breaks=c(0,0.5,1,1.5,2),sec.axis=sec_axis(~.,name="Nitrate+Nitrite",breaks=NULL,labels=NULL))+
  theme_minimal()+
  theme(axis.text.x=element_blank(),
        axis.title.x=element_blank(),
        legend.position="none"
  )+
  ylab("Concentration, mg (N or P)/L")
f2.<-f2+facet_grid(cols=vars(Sample.Pt.Descr))+theme(strip.text.x = element_blank())
f2.
f3<-ggplot(data = mdtan3,aes(Date))+
  geom_point(aes(y=value))+
  geom_line(aes(y=value))+
  geom_line(aes(y=Medians,linetype="Medians"),color="steelblue", size=1)+
  geom_line(aes(y=TPO4_BR,color="BR"))+ #TPO4
  geom_ribbon(data=e.s3,aes(ymin=0.190,ymax=0.270),fill='green',alpha=0.2)+
  scale_y_continuous(limits=c(0,0.4),breaks=c(0,0.1,0.2,0.3,0.4),sec.axis=sec_axis(~.,name="Total Phosphorous",breaks=NULL,labels=NULL))+
  theme_minimal()+
  scale_linetype_manual(values=2)+
  scale_color_manual(values="Green")+
  scale_x_date(limits = as.Date(c('2020-10-05','2021-04-30')),
               date_labels="%b '%y",
               breaks=brks
  )+
  theme(axis.text.x= element_text(angle = 45, hjust = 1),
        axis.title.x=element_blank(),
        axis.title.y.left=element_blank(),
        legend.position="bottom",
        legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  labs(color="Int'l BMP median",
       linetype="Statistic")+
  ylab("Concentration, mg (N or P)/L")
f3.<-f3+facet_grid(cols=vars(Sample.Pt.Descr))+theme(strip.text.x = element_blank())
f3.
ggarrange(f1.,f2.,f3., ncol=1)
#### All Metals ####
# OWQC (site 5 only) IBMP (bioretention only) for cadmium, copper, lead, zinc
dta2<-read.csv("CdZnPbCuWQC3.csv") #import csv with 4 metals
dta2<-rename(dta2,Date=Collection.Date)#rename variable
dta2$Date<-mdy(dta2$Date) # convert factor to Date
str(dta2)

dta2$Metal <- factor(dta2$Metal, levels=c("Zinc","Copper","Lead","Cadmium")) #re-prioritize metals
dta2<-dta2 %>% #Create IBMPdb metal medians for bioretention type
  mutate(TZn_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Zinc",12.8,NA),
         DZn_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Zinc",12.5,NA),
         TPb_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Lead",0.932,NA),
         DPb_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Lead",0.0739,NA),
         TCu_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Copper",7.12,NA),
         DCu_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Copper",7.54,NA),
         TCd_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Cadmium",0.0825,NA),
         DCd_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Cadmium",0.0668,NA))
f.s1<-dta2 %>%# Create special ribbon dataframes
  filter(Sample.Pt.Descr=="3_PostTerr",Metal=='Zinc')
f.s2<-dta2 %>%
  filter(Sample.Pt.Descr=="3_PostTerr",Metal=='Lead')
f.s3<-dta2 %>%
  filter(Sample.Pt.Descr=="3_PostTerr",Metal=='Copper')
f.s4<-dta2 %>%
  filter(Sample.Pt.Descr=="3_PostTerr",Metal=='Cadmium')

dta2<-dta2 %>% # Compute Total Metal medians Monitoring Points
  group_by(Sample.Pt.Descr,Metal)%>%
  mutate(Total.Median=median(Total))
dta2<-dta2 %>% # Compute Dissolved Metal medians Monitoring Points
  group_by(Sample.Pt.Descr,Metal)%>%
  mutate(Dissolved.Median=median(Diss))

l<-ggplot(data=dta2,aes(Date))+
  geom_point(aes(y=Total,fill="Total"),shape=21,color="transparent")+
  geom_line(aes(y=Total),color="red")+
  geom_point(aes(y=Diss,fill="Dissolved"),shape=21,color="transparent")+
  geom_line(aes(y=Diss),color="blue")+
  geom_line(aes(y=Acute.IWQC,linetype="Acute IWQC",color="Acute IWQC"))+
  geom_line(aes(y=Chronic.IWQC,linetype="Chronic IWQC",color="Chronic IWQC"))+
  geom_line(aes(y=Acute.WQC,linetype="Acute WQC",color="Acute WQC"))+
  geom_line(aes(y=Chronic.WQC,linetype="Chronic WQC",color="Chronic WQC"))+
  geom_line(data=dta2, aes(x=Date,y=TZn_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #Bioretention TZn median
  geom_ribbon(data=f.s1,aes(ymin=11,ymax=14),fill='red',alpha=0.2)+ #RP CFs
  geom_line(data=dta2, aes(x=Date,y=DZn_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #DZn
  geom_ribbon(data=f.s1,aes(ymin=9,ymax=13.8),fill='blue',alpha=0.2)+ #RP CFs
  geom_line(data=dta2, aes(x=Date,y=TPb_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #TPb
  geom_ribbon(data=f.s2,aes(ymin=0.723,ymax=1.07),fill='red',alpha=0.2)+ #RP CFs
  geom_line(data=dta2, aes(x=Date,y=DPb_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #DPb
  geom_ribbon(data=f.s2,aes(ymin=0.0506,ymax=0.0878),fill='blue',alpha=0.2)+ #RP CFs
  geom_line(data=dta2, aes(x=Date,y=TCu_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #TCu
  geom_ribbon(data=f.s3,aes(ymin=6.4,ymax=8.2),fill='red',alpha=0.2)+ #RP CFs
  geom_line(data=dta2, aes(x=Date,y=DCu_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #DCu
  geom_ribbon(data=f.s3,aes(ymin=6.5,ymax=8.4),fill='blue',alpha=0.2)+ #RP CFs
  geom_line(data=dta2, aes(x=Date,y=TCd_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #TCd
  geom_ribbon(data=f.s4,aes(ymin=0.0647,ymax=0.1),fill='red',alpha=0.2)+ #RP CFs
  geom_line(data=dta2, aes(x=Date,y=DCd_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #DCd
  geom_ribbon(data=f.s4,aes(ymin=0.0444,ymax=0.0885),fill='blue',alpha=0.2)+ #RP CFs
  theme_minimal()+
  scale_linetype_manual(breaks=c("Acute IWQC","Chronic IWQC","Acute WQC","Chronic WQC","Intl BMP BR"),values=c(1,2,1,2,1),labels=c("Acute IWQC","Chronic IWQC","Acute WQC","Chronic WQC","Intl BMP BR"))+
  scale_color_manual(breaks=c("Acute IWQC","Chronic IWQC","Acute WQC","Chronic WQC","Intl BMP BR"),values=c("salmon3","salmon3","darkgray","darkgray","green"))+
  scale_fill_manual(breaks=c("Total","Dissolved"),values=c("Total"="red","Dissolved"="blue"))+
  scale_x_date(limits = as.Date(c("2020-10-01","2021-04-30")),
               date_labels="%b '%y",
               date_breaks="1 month"
  )+
  theme(
    axis.text.x = element_text(angle = 45, hjust=1),
    axis.title.x=element_blank(),
    legend.box="horizontal",
    legend.position="bottom",
    legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  labs(colour="Criterion",
       linetype="Criterion",
       fill="Fraction")+
  ylab("Concentration, ug/L")
l+facet_grid(Metal~Sample.Pt.Descr, scales="free") 
#### All Metals - cadmium ####
dta3<-read.csv("CdZnPbCuWQC3.csv") #import csv with 4 metals
dta3<-rename(dta3,Date=Collection.Date)#rename variable
dta3$Date<-mdy(dta3$Date) # convert factor to Date
str(dta3)

dta3$Metal <- factor(dta3$Metal, levels=c("Zinc","Copper","Lead","Cadmium")) #re-prioritize metals
dta3<-dta3%>%
  filter(Metal=="Zinc"|Metal=="Lead"|Metal=="Copper") # Remove cadmium
dta3<-dta3 %>% #Create IBMPdb metal medians for bioretention type
  mutate(TZn_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Zinc",12.8,NA),
         DZn_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Zinc",12.5,NA),
         TPb_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Lead",0.932,NA),
         DPb_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Lead",0.0739,NA),
         TCu_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Copper",7.12,NA),
         DCu_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&Metal=="Copper",7.54,NA))

g.s1<-dta3 %>%# Create special ribbon dataframes
  filter(Sample.Pt.Descr=="3_PostTerr",Metal=='Zinc')
g.s2<-dta3 %>%
  filter(Sample.Pt.Descr=="3_PostTerr",Metal=='Lead')
g.s3<-dta3 %>%
  filter(Sample.Pt.Descr=="3_PostTerr",Metal=='Copper')

dta3<-dta3 %>% # Compute Total Metal medians Monitoring Points
  group_by(Sample.Pt.Descr,Metal)%>%
  mutate(Total.Median=median(Total))
dta3<-dta3 %>% # Compute Dissolved Metal medians Monitoring Points
  group_by(Sample.Pt.Descr,Metal)%>%
  mutate(Dissolved.Median=median(Diss))

m<-ggplot(data=dta3,aes(Date))+
  geom_point(aes(y=Total,fill="Total"),shape=21,color="transparent")+
  geom_line(aes(y=Total),color="red")+
  geom_point(aes(y=Diss,fill="Dissolved"),shape=21,color="transparent")+
  geom_line(aes(y=Diss),color="blue")+
  geom_line(aes(y=Acute.IWQC,linetype="Acute IWQC",color="Acute IWQC"))+
  geom_line(aes(y=Chronic.IWQC,linetype="Chronic IWQC",color="Chronic IWQC"))+
  geom_line(aes(y=Acute.WQC,linetype="Acute WQC",color="Acute WQC"))+
  geom_line(aes(y=Chronic.WQC,linetype="Chronic WQC",color="Chronic WQC"))+
  geom_line(data=dta3, aes(x=Date,y=TZn_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #Bioretention TZn median
  geom_ribbon(data=g.s1,aes(ymin=11,ymax=14),fill='red',alpha=0.2)+ #RP CFs
  geom_line(data=dta3, aes(x=Date,y=DZn_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #DZn
  geom_ribbon(data=g.s1,aes(ymin=9,ymax=13.8),fill='blue',alpha=0.2)+ #RP CFs
  geom_line(data=dta3, aes(x=Date,y=TPb_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #TPb
  geom_ribbon(data=g.s2,aes(ymin=0.723,ymax=1.07),fill='red',alpha=0.2)+ #RP CFs
  geom_line(data=dta3, aes(x=Date,y=DPb_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #DPb
  geom_ribbon(data=g.s2,aes(ymin=0.0506,ymax=0.0878),fill='blue',alpha=0.2)+ #RP CFs
  geom_line(data=dta3, aes(x=Date,y=TCu_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #TCu
  geom_ribbon(data=g.s3,aes(ymin=6.4,ymax=8.2),fill='red',alpha=0.2)+ #RP CFs
  geom_line(data=dta3, aes(x=Date,y=DCu_BR,linetype="Intl BMP BR",color="Intl BMP BR"))+ #DCu
  geom_ribbon(data=g.s3,aes(ymin=6.5,ymax=8.4),fill='blue',alpha=0.2)+ #RP CFs
  theme_minimal()+
  scale_linetype_manual(breaks=c("Acute IWQC","Chronic IWQC","Acute WQC","Chronic WQC","Intl BMP BR"),values=c(1,2,1,2,1),labels=c("Acute IWQC","Chronic IWQC","Acute WQC","Chronic WQC","Intl BMP BR"))+
  scale_color_manual(breaks=c("Acute IWQC","Chronic IWQC","Acute WQC","Chronic WQC","Intl BMP BR"),values=c("salmon3","salmon3","darkgray","darkgray","green"))+
  scale_fill_manual(breaks=c("Total","Dissolved"),values=c("Total"="red","Dissolved"="blue"))+
  scale_x_date(limits = as.Date(c("2020-10-01","2021-04-30")),
               date_labels="%b '%y",
               breaks=brks
  )+
  theme(
    axis.text.x = element_text(angle = 45, hjust=1),
    axis.title.x=element_blank(),
    legend.box="horizontal",
    legend.position="bottom",
    legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  labs(colour="Criterion",
       linetype="Criterion",
       fill="Fraction")+
  ylab("Concentration, ug/L")
m+facet_grid(Metal~Sample.Pt.Descr, scales="free")
#Medians only, no WQC or BMP medians#
m2<-ggplot(data=dta3,aes(Date))+
  geom_point(aes(y=Total,fill="Total"),shape=21,color="transparent")+
  geom_line(aes(y=Total),color="red")+
  geom_point(aes(y=Diss,fill="Dissolved"),shape=21,color="transparent")+
  geom_line(aes(y=Diss),color="blue")+
  geom_line(aes(y=Total.Median,color="Total",linetype="Total"),size=1)+
  geom_line(aes(y=Dissolved.Median,color="Dissolved",linetype="Dissolved"),size=1)+
  theme_minimal()+
  scale_linetype_manual(breaks=c("Total","Dissolved"),values=c(2,2),labels=c("Total","Dissolved"))+
  scale_color_manual(breaks=c("Total","Dissolved"),values=c("darkred","darkblue"))+
  scale_fill_manual(breaks=c("Total","Dissolved"),values=c("Total"="red","Dissolved"="blue"))+
  scale_x_date(limits = as.Date(c("2020-10-01","2021-04-30")),
               date_labels="%b '%y",
               breaks=brks
  )+
  theme(
    axis.text.x = element_text(angle = 45, hjust=1),
    axis.title.x=element_blank(),
    legend.box="horizontal",
    legend.position="bottom",
    legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  labs(colour="Median",
       linetype="Median",
       fill="Fraction")+
  ylab("Concentration, ug/L")
m2+facet_grid(Metal~Sample.Pt.Descr, scales="free")
#### Solids ####
#IBMP (bioretention only) + Willamette River guidance value medians
dta4<-read.csv("TSSTDSTSBMPdb.csv")
class(dta4$Collection.Date)
dta4<-rename(dta4,Date=Collection.Date)#rename variable
dta4$Date<-mdy(dta4$Date)
class(dta4$Date)


mdta4<-melt(dta4,id.vars=c("Event","Sample.Pt.Descr","Date"),measure.vars=c('TDS','TSS','TS'))
str(mdta4)
mdta4$variable <- factor(mdta4$variable, levels=c('TSS','TDS','TS')) #re-prioritize
mdta4<-mdta4 %>% # Create IBMPdb medians
  mutate(TSS_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&variable=="TSS",10,NA),
         TDS_BR=ifelse(Sample.Pt.Descr=="3_PostTerr"&variable=="TDS",210,NA),
         TDS_guid=ifelse(Sample.Pt.Descr=="5_DS"&variable=="TDS",100.0,NA)) #Will. R guidance value
h.s1<-mdta4 %>%
  filter(Sample.Pt.Descr=="3_PostTerr",variable=='TSS')# create special df for site 3/5 ribbon
h.s2<-mdta4 %>%
  filter(Sample.Pt.Descr=="3_PostTerr",variable=='TDS')
h.s3<-mdta4 %>%
  filter(Sample.Pt.Descr=="5_DS",variable=='TDS')
mdta4<-mdta4 %>% # Compute median by Monitoring Points and variable
  group_by(Sample.Pt.Descr,variable)%>%
  mutate(Median=median(value))

n<-ggplot(data = mdta4,aes(Date))+
  geom_point(aes(y=value))+
  geom_line(aes(y=value))+
  geom_line(aes(y=TSS_BR,color="BR"))+ #Bioretention TSS median
  geom_ribbon(data=h.s1,aes(ymin=8,ymax=11),fill='green',alpha=0.2)+ #RP CFs
  geom_line(aes(y=TDS_BR,color="BR"))+ #TDS
  geom_ribbon(data=h.s2,aes(ymin=175,ymax=298),fill='green',alpha=0.2)+ #RP CFs
  geom_line(aes(y=TDS_guid,color="Will. River Guidance"))+
  geom_line(aes(y=Median,linetype="Median"),color="steelblue",size=1)+
  theme_minimal()+
  scale_linetype_manual(breaks=c("Median"),values=c(2),label=c("Median"))+
  scale_color_manual(breaks=c("BR","Will. River Guidance"),values=c("Green","Red"))+
  scale_x_date(limits = as.Date(c("2020-10-05","2021-04-30")),
               date_labels="%b '%y",
               breaks=brks
  )+
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        axis.title.x=element_blank(),
        legend.position="bottom",
        legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  labs(linetype="Statistic",color="Criterion")+
  ylab("Concentration, mg/L")

n+facet_grid(variable~Sample.Pt.Descr, scales="free")

#### Concentration/Mass Reduction figures ####
df<-read.csv("DailyPrecip.csv") #import Station KCMO3 Precip Data
df2<-read.csv("PollRed2.csv") #import concentration/Mass reduction Data
df3<-read.csv("ClackFlowDaily.csv") #import Clackamas River flow
df4<-read.csv("Campaign.csv") #import event dates for identifying arrows in rain/flow figures

#Legend Shrinker function for figures#
SmLeg <- function(myPlot, pointSize = 1.5, textSize = 10, spaceLegend = 0.4) {
  myPlot +
    guides(shape = guide_legend(override.aes = list(size = pointSize)),
           color = guide_legend(override.aes = list(size = pointSize))) +
    theme(legend.title = element_text(size = textSize), 
          legend.text  = element_text(size = textSize),
          legend.key.size = unit(spaceLegend, "lines"))
}

Nuts<-df2 %>%
  filter(Pollutant=="Ammonia"|Pollutant=="Nitrate-Nitrite"|Pollutant=="Total Phosphorous")# Nuts df
Metals<-df2 %>%
  filter(Pollutant=="Total Cu"|Pollutant=="Total Pb"|Pollutant=="Total Zn"|Pollutant=="Diss Cu"|Pollutant=="Diss Pb"|Pollutant=="Diss Zn")# Metals df
Solids<-df2 %>% #Solids df
  filter(Pollutant=="Total Suspended Solids"|Pollutant=="Total Diss Solids"|Pollutant=="Total Solids")

####Rain Figure/panel####
df3$Date<-mdy(df3$Date)
df4$Date<-mdy(df4$Date)
df$Date<-mdy(df$Date)
df4<-df4 %>%
  left_join(df %>%
              select(Date,Daily.Total..in),
            by="Date")#add Precip to df4
df4<-df4 %>%
  left_join(df3 %>%
              select(Date,Daily.Volume.M.cf),
            by="Date")#Add Clack Flow to df4
class(df$Date)
df$Date<-ymd(df$Date)
str(df)
df <- df %>% # create new variables for summary statistics
  mutate(Year=as.factor(year(Date)), #Year, duh
          DOY=yday(Date),#DOY = Day of Year
          MOY=as.factor(months(Date)),#MOY = Month of Year
          MDay=(format(as.Date(Date),"%b-%d"))) #Month and day
dfM<-df %>% #Add historical monthly averages variable, excluding study period
  filter(Date<"2020-10-01")%>%
  group_by(MOY) %>%
  summarize(MAvg=mean(Daily.Total..in))
df<-df %>%
  left_join(dfM, by="MOY") #add MAvg to df table
dfD<-df %>% #Add historical daily averages variable, excluding study period
  filter(Date<"2020-10-01")%>%
  group_by(DOY) %>%
  summarize(DAvg=mean(Daily.Total..in))
df<-df %>%
  left_join(dfD, by="DOY") #add DAvg to df table

rain<-ggplot(data=df)+
  geom_line(aes(x=Date,y=Daily.Total..in,color="Study Period",linetype="Study Period"))+
  geom_line(aes(x=Date,y=DAvg, color="Historic Average, by day",linetype="Historic Average, by day"))+
  geom_line(aes(x=Date,y=MAvg, color="Historic Average, by month",linetype="Historic Average, by month"))+
  scale_x_date(limits = as.Date(c("2020-10-05","2021-04-15")),date_labels="%b '%y",breaks="months")+ #'limits' only shows study period of data
  theme_minimal()+
  scale_color_manual(breaks=c("Study Period","Historic Average, by day","Historic Average, by month"),values=c("black","red","blue"))+
  scale_linetype_manual(breaks=c("Study Period","Historic Average, by day","Historic Average, by month"),values=c(1,2,5))+
  theme(axis.title.x=element_blank(),
        axis.text.x = element_text(angle = 60,vjust=1,hjust=1),
        legend.title=element_blank(),
        legend.position=c(0.85,0.80),
        legend.direction="vertical",
        legend.box="vertical",
        legend.background = element_rect(fill="gray90", size=.5, linetype="dotted"))+
  ylab("Daily Precipitation, in.")+
  labs(color='',linetype='')
rain

rain1<-ggplot(data=df)+
  geom_line(aes(x=Date,y=Daily.Total..in,color="Study Period",linetype="Study Period"))+
  geom_line(aes(x=Date,y=DAvg, color="Historic Average, by day",linetype="Historic Average, by day"))+
  geom_line(aes(x=Date,y=MAvg, color="Historic Average, by month",linetype="Historic Average, by month"))+
  scale_x_date(limits = as.Date(c("2020-10-05","2021-04-15")),date_labels="%b '%y",breaks="months")+ #'limits' only shows study period of data
  theme_minimal()+
  scale_color_manual(breaks=c("Study Period","Historic Average, by day","Historic Average, by month"),values=c("black","red","blue"))+
  scale_linetype_manual(breaks=c("Study Period","Historic Average, by day","Historic Average, by month"),values=c(1,2,5))+
  theme(axis.title.x=element_blank(),
        axis.text.x = element_text(angle = 60,vjust=1,hjust=1),
        legend.title=element_blank(),
        legend.position=c(0.85,0.80),
        legend.direction="vertical",
        legend.box="vertical",
        legend.background = element_rect(fill="gray90", size=.5, linetype="dotted"))+
  ylab("Daily Precipitation, in.")+
  labs(color='',linetype='')+
  geom_segment(data=df4,aes(x=Date,xend=Date,y=Daily.Total..in+0.5,yend=Daily.Total..in+0.2),color="darkgreen",arrow=arrow(length=unit(0.2,"cm")))+
  geom_text(data=df4,aes(x=Date,y=Daily.Total..in+0.6,label=Event),color="darkgreen")

####Clackamas River Flow figure/panel####
str(df3)
df3$Date<-ymd(df3$Date)
str(df3)
df3 <- df3 %>% # create new variables for summary statistics
  mutate(Year=as.factor(year(Date)), #Year, duh
         DOY=yday(Date),#DOY = Day of Year
         MOY=as.factor(months(Date)),#MOY = Month of Year
         MDay=(format(as.Date(Date),"%b-%d"))) #Month and day
df3M<-df3 %>% #Add historical monthly averages variable, excluding study period
  filter(Date<"2020-10-01")%>%
  group_by(MOY) %>%
  summarize(MAvg=mean(Daily.Volume.M.cf))
df3<-df3 %>%
  left_join(df3M, by="MOY") #add MAvg to df table
df3D<-df3 %>% #Add historical daily averages variable, excluding study period
  filter(Date<"2020-10-01")%>%
  group_by(DOY)%>%
  summarise(DAvg=mean(Daily.Volume.M.cf))
df3<-df3 %>%
  left_join(df3D, by="DOY") #add DAvg to df table
Cflow<-ggplot(data=df3)+
  geom_line(aes(x=Date,y=Daily.Volume.M.cf,color="Study Period",linetype="Study Period"))+
  geom_line(aes(x=Date,y=DAvg, color="Historic Average, by day",linetype="Historic Average, by day"))+
  geom_line(aes(x=Date,y=MAvg, color="Historic Average, by month",linetype="Historic Average, by month"))+
  scale_x_date(limits = as.Date(c("2020-10-05","2021-04-15")),date_labels="%b '%y",breaks="months")+ #'limits' only shows study period of data
  theme_minimal()+
  scale_color_manual(breaks=c("Study Period","Historic Average, by day","Historic Average, by month"),values=c("black","red","blue"))+
  scale_linetype_manual(breaks=c("Study Period","Historic Average, by day","Historic Average, by month"),values=c(1,2,5))+
  theme(axis.title.x=element_blank(),
        axis.text.x = element_text(angle = 60,vjust=1,hjust=1),
        legend.title=element_blank(),
        legend.position=c(0.85,0.80),
        legend.direction="vertical",
        legend.box="vertical",
        legend.background = element_rect(fill="gray90", size=.5, linetype="dotted"))+
  ylab(expression(paste("Daily Total Volume",", ",10^{"6"}, " cf")))+
  labs(color='',linetype='')+
  geom_segment(data=df4,aes(x=Date,xend=Date,y=Daily.Volume.M.cf+400,yend=Daily.Volume.M.cf+100),color="darkgreen",arrow=arrow(length=unit(0.2,"cm")))+
  geom_text(data=df4,aes(x=Date,y=Daily.Volume.M.cf+500,label=Event),color="darkgreen")

Cflow
#### Study Site Flows ####
SS1<-read.csv("Site1.csv")
SS1$DateTime<-mdy_hm(SS1$DateTime)
SS3<-read.csv("Site3.csv")
SS3$DateTime<-mdy_hm(SS3$DateTime)
SS4<-read.csv("Site4.csv")
SS4$DateTime<-mdy_hm(SS4$DateTime)
SS13<-SS1%>%
  left_join(SS3,by="DateTime")
SS13<-SS13%>%
  rename(Site.1=Flow.x)%>%
  rename(Site.3=Flow.y)%>%
  rename(Date=DateTime)
str(SS13)
SS13$Date<-ymd_hms(SS13$Date)
df4<-read.csv("Campaign.csv")
df4$Date<-as.POSIXct(mdy(df4$Date), force_tz="")
str(df4)
df4<-df4 %>%
  left_join(SS13,
            by="Date")#Add Site Flows to df4
s<-ggplot(data=SS13)+
  geom_line(aes(x=Date,y=Site.1,color="1"))+
  geom_line(aes(x=Date,y=Site.3, color="3"))+
  scale_x_datetime(limits = as_datetime(c("2020-12-01","2021-02-01")),date_labels="%b %d",breaks="3 days")+ #'limits' only shows study period of data
  theme_minimal()+
  scale_color_manual(values=c("red","blue"))+
  theme(axis.title.x=element_blank(),
    axis.text.x = element_text(angle = 60,vjust=1,hjust=1),
    legend.position=c(0.85,0.80),
    legend.direction="horizontal",
    legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  ylim(0,7.5)+
  labs(color="Site")+
  ylab("Discharge, cfs")+
  geom_segment(data=df4,aes(x=Date,xend=Date,y=Site.1+1,yend=Site.1+0.5),color="darkgreen",arrow=arrow(length=unit(0.2,"cm")))+
  geom_text(data=df4,aes(x=Date,y=Site.1+1.5,label=Event),color="darkgreen")
SmLeg(s)
s1<-ggplot(data=SS13)+
  geom_line(aes(x=Date,y=Site.1,color="1"))+
  geom_line(aes(x=Date,y=Site.3, color="3"))+
  scale_x_datetime(limits = as_datetime(c("2020-12- 01","2021-02-01")),date_labels="%b %d '%y",breaks="1 weeks")+ #'limits' only shows study period of data
  theme_minimal()+
  scale_color_manual(values=c("red","blue"))+
  theme(axis.title.x=element_blank(),
        axis.text.x = element_text(angle = 60,vjust=1,hjust=1),
        legend.position=c(0.85,0.80),
        legend.direction="horizontal",
        legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  ylim(0,7.5)+
  labs(color="Site")+
  ylab("Discharge, cfs")
s1
#### Nutrients versus Rain ####
library(cowplot)
library(egg)
Nuts$Date<-mdy(Nuts$Date)
o<-ggplot(data=Nuts) +
  geom_col( aes(x=Date, y=Reduction,fill=Pollutant),position="dodge2") +
  geom_hline(aes(yintercept=0),linetype=2)+
  scale_x_date(limits = as.Date(c("2020-10-05","2021-04-15")),date_breaks="1 month")+
  theme_minimal()+
  scale_y_continuous(limits=c(-350,100))+  
  theme(
    axis.text.x=element_blank(),
    axis.title.x=element_blank(),
    legend.position = c(0.15,0.05),
    legend.direction="vertical",
    legend.box="vertical",
    legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  ylab("Percent Reduction")
o
og<-SmLeg(o)+facet_grid(rows=vars(Type),scales="free")
og
ggarrange(og, rain, ncol=1,heights=c(0.66,0.33))

Metalsc<-Metals %>% filter(Type=="Concentration")
Metalsc$Pollutant<- recode_factor(Metalsc$Pollutant,"Total Cu"="T Copper","Diss Cu"="D Copper","Total Pb"="T Lead","Diss Pb"="D Lead","Total Zn"="T Zinc","Diss Zn"="D Zinc") #rename
str(Metalsc)
Metalsc$Pollutant <- factor(Metalsc$Pollutant, levels=c("T Copper","D Copper","T Lead","D Lead","T Zinc","D Zinc")) #re-prioritize metals
levels(Metalsc$Pollutant)
str(Metalsc)
Metalsc$Date<-mdy(Metalsc$Date)
p<-ggplot(data=Metalsc) +
  geom_col( aes(x=Date, y=Reduction,fill=Pollutant),position="dodge2") +
  geom_hline(aes(yintercept=0),linetype=2)+
  scale_x_date(limits = as.Date(c("2020-10-05","2021-04-15")),date_breaks="1 month")+
  theme_minimal()+
  scale_y_continuous(limits=c(-350,100))+
  scale_fill_manual(values=c("darkred","red","darkblue","blue","darkgreen","green"))+
  theme(
    axis.text.x=element_blank(),
    axis.title.x=element_blank(),
    legend.position = c(0.5,0.2),
    legend.direction="horizontal",
    legend.box="vertical",
    legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  ylab("Percent Reduction \n by Concentration")
p
p<-SmLeg(p)+guides(fill=guide_legend(ncol=3))
p
ggarrange(p, rain, ncol=1,heights=c(0.66,0.33))
## by mass ##
Metalsm<-Metals %>% filter(Type=="Mass")
Metalsm$Pollutant<- recode_factor(Metalsm$Pollutant,"Total Cu"="T Copper","Diss Cu"="D Copper","Total Pb"="T Lead","Diss Pb"="D Lead","Total Zn"="T Zinc","Diss Zn"="D Zinc") #rename
Metalsm$Pollutant <- factor(Metalsm$Pollutant, levels=c("T Copper","D Copper","T Lead","D Lead","T Zinc","D Zinc")) #re-prioritize metals
levels(Metalsm$Pollutant)
Metalsm$Date<-mdy(Metalsm$Date)
q<-ggplot(data=Metalsm) +
  geom_col( aes(x=Date, y=Reduction,fill=Pollutant),position="dodge2") +
  geom_hline(aes(yintercept=0),linetype=2)+
  scale_x_date(limits = as.Date(c("2020-10-05","2021-04-15")),date_breaks="1 month")+
  theme_minimal()+
  scale_y_continuous(limits=c(-350,100))+
  scale_fill_manual(values=c("darkred","red","darkblue","blue","darkgreen","green"))+
  theme(
    axis.text.x=element_blank(),
    axis.title.x=element_blank(),
    legend.position = c(0.5,0.2),
    legend.direction="horizontal",
    legend.box="vertical",
    legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  ylab("Percent Reduction \n by Mass")
q
q<-SmLeg(q)+guides(fill=guide_legend(ncol=3))
q
ggarrange(q, rain, ncol=1,heights=c(0.66,0.33))

#### Solids versus Rain ####
Solids$Date<-mdy(Solids$Date)
Solids$Pollutant<- recode_factor(Solids$Pollutant,"Total Diss Solids"="Total Dissolved Solids") #rename
Solids$Pollutant <- factor(Solids$Pollutant, levels=c("Total Solids","Total Suspended Solids","Total Dissolved Solids"))
levels(Solids$Pollutant)
r<-ggplot(data=Solids)+
  geom_col( aes(x=Date, y=Reduction,fill=Pollutant),position="dodge2")+
  geom_hline(aes(yintercept=0),linetype=2)+
  scale_x_date(limits = as.Date(c("2020-10-05","2021-04-15")),date_breaks="1 month")+
  theme_minimal()+
  scale_fill_manual(values=c("black","darkred","red"))+
  scale_y_continuous(limits=c(-200,100))+
  theme(
    axis.text.x=element_blank(),
    axis.title.x=element_blank(),
    legend.position = c(0.15,0.15),
    legend.direction="vertical",
    legend.box="vertical",
    legend.background = element_rect(fill="gray90", size=.5, linetype="dotted")
  )+
  ylab("Percent Reduction")
r
r<-SmLeg(r)+facet_grid(rows=vars(Type))
r
ggarrange(r, rain, ncol=1,heights=c(0.66,0.33))