Outflow_RCP8.5_NoIFR.Rmd

---
title: "Output_flow"
author: "Gustavo Facincani Dourado"
date: "6/12/2020"
output: html_document
---

```{r}
library(reshape2)
library(lfstat)
library(scales)

Livneh_outflow_Mer <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/merced/historical/Livneh/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,11)]%>%
  filter(between(node, as.Date("1980-10-01"), as.Date("2010-09-30")))

colnames(Livneh_outflow_Mer) <- c("node", "Livneh (Historical)")
Livneh_outflow_Mer


CanESM2_outflow_Mer <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/merced/gcms/CanESM2_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,11)]

colnames(CanESM2_outflow_Mer) <- c("node", "CanESM2")
CanESM2_outflow_Mer

#CanESM2_MIF_Mer <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/merced/gcms/CanESM2_rcp85/InstreamFlowRequirement_Min Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,2)]

#CanESM2_MIF_Mer
#CanESM2_MIF_Mer

CNRM_outflow_Mer <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/merced/gcms/CNRM-CM5_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,11)]

colnames(CNRM_outflow_Mer) <- c("node", "CNRM-CM5")
CNRM_outflow_Mer

HadGEM2_outflow_Mer <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/merced/gcms/HadGEM2-ES_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,11)]

colnames(HadGEM2_outflow_Mer) <- c("node", "HadGEM2-ES")
HadGEM2_outflow_Mer

MIROC5_outflow_Mer <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/merced/gcms/MIROC5_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,11)]

colnames(MIROC5_outflow_Mer) <- c("node", "MIROC5")
MIROC5_outflow_Mer


```

```{r}
Livneh_Mer_outflow <- Livneh_outflow_Mer 
Livneh_Mer_outflow <- melt(Livneh_Mer_outflow, id = "node")
Livneh_Mer_outflow

Outflow_Mer <- cbind(CanESM2_outflow_Mer, CNRM_outflow_Mer[(-1)], 
                HadGEM2_outflow_Mer[(-1)], MIROC5_outflow_Mer[(-1)])
Outflow_Mer

```

```{r}
Outflow_Mer1 <- melt(Outflow_Mer, id = "node") 
Outflow_Mer1



Outflow_Mer2 <- rbind(Outflow_Mer1, Livneh_Mer_outflow) %>%
  mutate(Month = as.factor(format(as.Date(node, format = "%B"), "%b")),
         Year = water_year(node, origin = "usgs"),
         Scenario = as.factor("RCP 8.5")) %>%
group_by(Year, Month, variable) %>%
   dplyr::summarize(MonthlyTotal = sum(value)) 
#    group_by(Year, Scenario, variable) %>%
#  dplyr::mutate(AnnualTotal = sum(value))
Outflow_Mer2

Outflow_Mer2$variable <- factor(Outflow_Mer2$variable, levels = c("Livneh (Historical)", "HadGEM2-ES", "CanESM2", "MIROC5", "CNRM-CM5")) 
Outflow_Mer2 


Outflow_Mer3 <- Outflow_Mer2 %>%
  group_by(Year, variable) %>%
   dplyr::summarize(AnnualTotal = sum(MonthlyTotal)) 
Outflow_Mer3

```


```{r}
#RCP 8.5 Statistics
#Monthly
all_Outflow_monthly
All_Mon_Stats <- aggregate(MonthlyTotal ~ variable+Basin, all_Outflow_monthly, summary)
write.csv(All_Mon_Stats, "All_Monthly_Stats_Outflow_RCP8.5.csv")

all_Outflow_monthly
All_Mon_Stats <- aggregate(AnnualTotal ~ variable+Basin, all_Outflow_annual, summary)
write.csv(All_Mon_Stats, "All_Annual_Stats_Outflow_RCP8.5.csv")
```


```{r}
#Merced annual outflow
Outflow_Mer2 %>%
  filter(variable == "Livneh") %>%

ggplot() + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_boxplot(aes( x = Month, y = MonthlyTotal, outlier.alpha = 0.3)) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
  #  scale_y_continuous(limits = c(0, NA),
   #                   expand = c(0, NA)) +
     scale_y_log10()+
 
  scale_x_discrete(limits=c("Oct","Nov","Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep")) +
  
  labs(title = "Merced River",
    subtitle = "CNRM-CM5",
    x = element_blank(),
       y = "Total River Outflow (mcm)") + #name of x axis
  theme(legend.title = element_blank(),
    legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -19),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))
#  png("Mer_outflow_CNRM-CM5.png", units ="in", width=8, height=5, res = 300)

```

```{r}
#Merced outflow
f <- function(x) {
  r <- quantile(x, probs = c(0, 0.25, 0.5, 0.75, 1))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  r
}
```

```{r}
ggplot(Outflow_Mer2) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_boxplot(aes( x = Month, y = MonthlyTotal, fill = variable), outlier.alpha = 0.3) + #plot monthly observed data in greenish blue
#  stat_summary(aes( x = Month, y = value, fill = variable),fun.data = f, outlier.alpha = 0.3, fun.y = sum, geom = "boxplot", stat = "identity")+
   #geom_line(aes(x = node, y = value)) +
  #scale_y_continuous(limits = c(0, NA),
   #                expand = c(0, NA),
    #               labels = scales::comma) +
  scale_y_log10(breaks=c(0,10,100,1000), labels=comma)+ 

    scale_x_discrete(limits=c("Oct","Nov","Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +

  labs(fill = "RCP 8.5",
       title = "Merced River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/month)") + #name of x axis
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -17),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
png("Mer_montlhy_outflow_GCMs_logscale.png", units ="in", width=8, height=5, res = 300)
```

```{r}
#annual aggreagate
ggplot(Outflow_Mer3) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_violin(draw_quantiles = c(.25, .5, .75), aes( x = variable, y = AnnualTotal, fill = variable)) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
    scale_y_continuous(limits = c(0, NA),
                      expand = c(0, NA),
                      labels = scales::comma) +
    scale_x_discrete(limits=c("Livneh (Historical)", "HadGEM2-ES", "CanESM2", "MIROC5", "CNRM-CM5"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +
# scale_y_log10(breaks=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000), labels=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000))+
  labs(subtitle = "RCP 8.5 Scenario",
       title = "Merced River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/year)") + #name of x axis
  theme(legend.position = "none",
        #legend.direction = "horizontal",
        #legend.box.margin = margin(t = -19),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
png("Mer_annual_outflow_GCMs.png", units ="in", width=8, height=5, res = 300)
```


```{r}
Livneh_Outflow_Tuo <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/tuolumne/historical/Livneh/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1, 7)]

Livneh_Outflow_Tuo <- Livneh_Outflow_Tuo %>%
# mutate(Livneh = rowSums(Livneh_Outflow_Tuo[c(2:5)])) %>%
      filter(between(node, as.Date("1980-10-01"), as.Date("2010-09-30")))

colnames(Livneh_Outflow_Tuo) <- c("node", "Livneh (Historical)")
Livneh_Outflow_Tuo

CanES_Outflow_Tuo <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/tuolumne/gcms/CanESM2_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1, 7)]

#CanES_Outflow_Tuo <- CanES_Outflow_Tuo %>%
# mutate(CanESM2 = rowSums(CanES_Outflow_Tuo[c(2:5)]))
colnames(CanES_Outflow_Tuo) <- c("node", "CanESM2")
CanES_Outflow_Tuo 

CNRM_Outflow_Tuo <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/tuolumne/gcms/CNRM-CM5_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1, 7)] 

#CNRM_Outflow_Tuo <- CNRM_Outflow_Tuo %>%
# mutate(`CNRM-CM5` = rowSums(CNRM_Outflow_Tuo[c(2:5)]))
colnames(CNRM_Outflow_Tuo) <- c("node", "CNRM-CM5")
CNRM_Outflow_Tuo

HadGEM_Outflow_Tuo <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/tuolumne/gcms/HadGEM2-ES_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1, 7)]

#HadGEM_Outflow_Tuo <- HadGEM_Outflow_Tuo %>%
# mutate(`HadGEM2-ES` = rowSums(HadGEM_Outflow_Tuo[c(2:5)]))
colnames(HadGEM_Outflow_Tuo) <- c("node", "HadGEM2-ES")
HadGEM_Outflow_Tuo

MIROC_Outflow_Tuo <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/tuolumne/gcms/MIROC5_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1, 7)]

#MIROC_Outflow_Tuo <- MIROC_Outflow_Tuo %>%
# mutate(MIROC5 = rowSums(MIROC_Outflow_Tuo[c(2:5)]))
colnames(MIROC_Outflow_Tuo) <- c("node", "MIROC5")
MIROC_Outflow_Tuo
```

```{r}
Livneh_Tuo_Outflow <- Livneh_Outflow_Tuo
Livneh_Tuo_Outflow <- melt(Livneh_Tuo_Outflow, id = "node")
Livneh_Tuo_Outflow

Outflow_Tuo <- cbind(CanES_Outflow_Tuo, CNRM_Outflow_Tuo[-1], 
                HadGEM_Outflow_Tuo[-1], MIROC_Outflow_Tuo[-1])
Outflow_Tuo 

Outflow_Tuo2 <- melt(Outflow_Tuo, id = "node")
Outflow_Tuo2
```

```{r}

Outflow_Tuo2 <- rbind(Outflow_Tuo2, Livneh_Tuo_Outflow) %>%
  mutate(Month = as.factor(format(as.Date(node, format = "%B"), "%b")),
                  Year = water_year(node, origin = "usgs"),
         Scenario = as.factor("RCP 8.5")) %>%
group_by(Year, Month, variable) %>%
   dplyr::summarize(MonthlyTotal = sum(value))
Outflow_Tuo2

Outflow_Tuo2$variable <- factor(Outflow_Tuo2$variable, levels = c("Livneh (Historical)", "HadGEM2-ES", "CanESM2", "MIROC5", "CNRM-CM5")) 

Outflow_Tuo2

Outflow_Tuo3 <- Outflow_Tuo2 %>%
  group_by(Year, variable) %>%
   dplyr::summarize(AnnualTotal = sum(MonthlyTotal)) 
Outflow_Tuo3


```

```{r}
#Merced annual generation
Outflow_Tuo2 %>%
  filter(variable == "CanESM2") %>%

ggplot() + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_boxplot(aes( x = Month, y = MonthlyTotal), outlier.alpha = 0.3) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
#    scale_y_continuous(limits = c(0, NA),
 #                     expand = c(0, NA)) +
   scale_y_log10()+
 
  scale_x_discrete(limits=c("Oct","Nov","Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep")) +
  
  labs(title = "Stanislaus River",
    subtitle = "CanESM2",
    x = element_blank(),
       y = "Total River Outflow (mcm)") + #name of x axis
  theme(legend.title = element_blank(),
    legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -19),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))
```

```{r}
#Merced annual generation
#Gen_Mer2 %>%
#  filter(variable == "CanESM2") %>%

ggplot(Outflow_Tuo2) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_boxplot(aes( x = Month, y = MonthlyTotal, fill = variable), outlier.alpha = 0.3) + #plot monthly observed data in greenish blue
#   geom_line(aes(x = node, y = value)) +
  #  scale_y_continuous(limits = c(0, NA),
 #                    expand = c(0, NA),
#                     labels = comma) +
   scale_y_log10(breaks=c(0, 10,100,1000), labels=scales::comma)+


    scale_x_discrete(limits=c("Oct","Nov","Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +

  labs(fill = "RCP 8.5",
       title = "Tuolumne River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/month)") + #name of x axis
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -17),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
  png("Tuo_outflow_GCMs_logscale.png", units ="in", width=8, height=5, res = 300)

```

```{r}
#annual aggreagate
ggplot(Outflow_Tuo3) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_violin(draw_quantiles = c(.25, .5, .75), aes( x = variable, y = AnnualTotal, fill = variable)) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
    scale_y_continuous(limits = c(0, NA),
                      expand = c(0, NA),
                      labels = scales::comma) +
#   scale_y_log10(breaks=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000), labels=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000))+
    scale_x_discrete(limits=c("Livneh (Historical)", "HadGEM2-ES", "CanESM2", "MIROC5", "CNRM-CM5"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +

  labs(subtitle = "RCP 8.5 Scenario",
       title = "Tuolumne River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/year)") + #name of x axis
  theme(legend.position = "none",
        #legend.direction = "horizontal",
        #legend.box.margin = margin(t = -19),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
png("Tuo_Outflow_GCMs2.png", units ="in", width=8, height=5, res = 300)
```


```{r}
Livneh_Outflow_Stn <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/stanislaus/historical/Livneh/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,3)] 

Livneh_Outflow_Stn <-Livneh_Outflow_Stn %>%
# mutate(Livneh = rowSums(Livneh_Outflow_Stn[c(2:11)])) %>%
      filter(between(node, as.Date("1980-10-01"), as.Date("2010-09-30")))

colnames(Livneh_Outflow_Stn) <- c("node", "Livneh (Historical)")
Livneh_Outflow_Stn

CanES_Outflow_Stn <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/stanislaus/gcms/CanESM2_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,3)]

#CanES_Outflow_Stn <- CanES_Outflow_Stn %>%
#   mutate(CanESM2 = rowSums(CanES_Outflow_Stn[c(2:11)]))
colnames(CanES_Outflow_Stn) <- c("node", "CanESM2")
CanES_Outflow_Stn

CNRM_Outflow_Stn <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/stanislaus/gcms/CNRM-CM5_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,3)]

#CNRM_Outflow_Stn <- CNRM_Outflow_Stn %>%
#   mutate(`CNRM-CM5` = rowSums(CNRM_Outflow_Stn[c(2:11)]))
colnames(CNRM_Outflow_Stn) <- c("node", "CNRM-CM5")
CNRM_Outflow_Stn


HadGEM_Outflow_Stn <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/stanislaus/gcms/HadGEM2-ES_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,3)]


#HadGEM_Outflow_Stn <- HadGEM_Outflow_Stn %>%
#   mutate(`HadGEM2-ES` = rowSums(HadGEM_Outflow_Stn[c(2:11)]))
colnames(HadGEM_Outflow_Stn) <- c("node", "HadGEM2-ES")
HadGEM_Outflow_Stn

MIROC_Outflow_Stn <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/stanislaus/gcms/MIROC5_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1, 3)]

#MIROC_Outflow_Stn <- MIROC_Outflow_Stn %>%
#   mutate(MIROC5 = rowSums(MIROC_Outflow_Stn[c(2:11)]))
colnames(MIROC_Outflow_Stn) <- c("node", "MIROC5")
MIROC_Outflow_Stn

```

```{r}
Livneh_Stn_Outflow <- Livneh_Outflow_Stn 
Livneh_Stn_Outflow <- melt(Livneh_Stn_Outflow, id = "node")
Livneh_Stn_Outflow


Outflow_Stn <- cbind(CanES_Outflow_Stn, CNRM_Outflow_Stn[-1], 
                HadGEM_Outflow_Stn[-1], MIROC_Outflow_Stn[-1])
Outflow_Stn
 

Outflow_Stn2 <- melt(Outflow_Stn, id = "node")
Outflow_Stn2
```

```{r}
Outflow_Stn2 <- rbind(Outflow_Stn2, Livneh_Stn_Outflow) %>%
  mutate(Month = as.factor(format(as.Date(node, format = "%B"), "%b")),
                           Year = water_year(node, origin = "usgs"),
         Scenario = as.factor("RCP 8.5")) %>%
group_by(Year, Month, variable) %>%
   dplyr::summarize(MonthlyTotal = sum(value))
Outflow_Stn2
  

Outflow_Stn2$variable <- factor(Outflow_Stn2$variable, levels = c("Livneh (Historical)", "HadGEM2-ES", "CanESM2", "MIROC5", "CNRM-CM5")) 
Outflow_Stn2

Outflow_Stn3 <- Outflow_Stn2 %>%
  group_by(Year, variable) %>%
   dplyr::summarize(AnnualTotal = sum(MonthlyTotal)) 

Outflow_Stn3
```

```{r}
Outflow_Stn2 %>%
  filter(variable == "CanESM2") %>%

ggplot() + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_boxplot(aes( x = Month, y = value, outlier.alpha = 0.3)) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
#    scale_y_continuous(limits = c(0, NA),
#                      expand = c(0, NA)) +
   scale_y_log10()+
 
  scale_x_discrete(limits=c("Oct","Nov","Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep")) +
  
  labs(title = "Stanislaus River",
    subtitle = "CanESM2",
    x = element_blank(),
       y = "Total River Outflow (mcm)") + #name of x axis
  theme(legend.title = element_blank(),
    legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -19),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))
#png("Stn_outflow_CNRM-CM5.png", units ="in", width=8, height=5, res = 300)
```

```{r}
ggplot(Outflow_Stn2) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_boxplot(aes( x = Month, y = MonthlyTotal, fill = variable), outlier.alpha = 0.3) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
    scale_y_continuous(limits = c(0, NA),
                      expand = c(0, NA),
                     labels = comma) +
#     scale_x_discrete(limits=c("Oct","Nov","Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +
#   scale_y_log10(breaks=c(0, 10,100,1000), labels=comma)+

  labs(fill = "RCP 8.5",
       title = "Stanislaus River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/month)") + #name of x axis
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -17),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
  png("Stn_outflow_GCMs.png", units ="in", width=8, height=5, res = 300)

```

```{r}
#annual aggreagate
ggplot(Outflow_Stn3) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_violin(draw_quantiles = c(.25, .5, .75), aes( x = variable, y = AnnualTotal, fill = variable)) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
    scale_y_continuous(limits = c(0, NA),
                      expand = c(0, NA),
                      labels = comma) +
#   scale_y_log10(breaks=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000), labels=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000))+
    scale_x_discrete(limits=c("Livneh (Historical)", "HadGEM2-ES", "CanESM2", "MIROC5", "CNRM-CM5"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +

  labs(subtitle = "RCP 8.5 Scenario",
       title = "Stanislaus River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/year)") + #name of x axis
  theme(legend.position = "none",
        #legend.direction = "horizontal",
        #legend.box.margin = margin(t = -19),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
png("Stn_outflow_GCMs2.png", units ="in", width=8, height=5, res = 300)
```

```{r}
Livneh_Outflow_USJ <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/upper_san_joaquin/historical/Livneh/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,3)] 

Livneh_Outflow_USJ <-Livneh_Outflow_USJ %>%
# mutate(Livneh = rowSums(Livneh_Outflow_Stn[c(2:11)])) %>%
      filter(between(node, as.Date("1980-10-01"), as.Date("2010-09-30")))

colnames(Livneh_Outflow_USJ) <- c("node", "Livneh (Historical)")
Livneh_Outflow_USJ

CanES_Outflow_USJ <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/upper_san_joaquin/gcms/CanESM2_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,3)]

#CanES_Outflow_Stn <- CanES_Outflow_Stn %>%
#   mutate(CanESM2 = rowSums(CanES_Outflow_Stn[c(2:11)]))
colnames(CanES_Outflow_USJ) <- c("node", "CanESM2")
CanES_Outflow_USJ

CNRM_Outflow_USJ <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/upper_san_joaquin/gcms/CNRM-CM5_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,3)]

#CNRM_Outflow_Stn <- CNRM_Outflow_Stn %>%
#   mutate(`CNRM-CM5` = rowSums(CNRM_Outflow_Stn[c(2:11)]))
colnames(CNRM_Outflow_USJ) <- c("node", "CNRM-CM5")
CNRM_Outflow_USJ


HadGEM_Outflow_USJ <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/upper_san_joaquin/gcms/HadGEM2-ES_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1,3)]


#HadGEM_Outflow_Stn <- HadGEM_Outflow_Stn %>%
#   mutate(`HadGEM2-ES` = rowSums(HadGEM_Outflow_Stn[c(2:11)]))
colnames(HadGEM_Outflow_USJ) <- c("node", "HadGEM2-ES")
HadGEM_Outflow_USJ

MIROC_Outflow_USJ <- read_csv("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/Pywr models/results/Binary IFRs x Prices/upper_san_joaquin/gcms/MIROC5_rcp85/Output_Flow_mcm.csv", col_types = cols(node = col_date(), .default = col_double()))[-c(1:3), c(1, 3)]

#MIROC_Outflow_Stn <- MIROC_Outflow_Stn %>%
#   mutate(MIROC5 = rowSums(MIROC_Outflow_Stn[c(2:11)]))
colnames(MIROC_Outflow_USJ) <- c("node", "MIROC5")
MIROC_Outflow_USJ

```

```{r}
Livneh_USJ_Outflow <- Livneh_Outflow_USJ 
Livneh_USJ_Outflow <- melt(Livneh_USJ_Outflow, id = "node")
Livneh_USJ_Outflow


Outflow_USJ <- cbind(CanES_Outflow_USJ, CNRM_Outflow_USJ[-1], 
                HadGEM_Outflow_USJ[-1], MIROC_Outflow_USJ[-1])
Outflow_USJ
 

Outflow_USJ2 <- melt(Outflow_USJ, id = "node")
Outflow_USJ2
```

```{r}
Outflow_USJ2 <- rbind(Outflow_USJ2, Livneh_USJ_Outflow) %>%
  mutate(Month = as.factor(format(as.Date(node, format = "%B"), "%b")),
                           Year = water_year(node, origin = "usgs"),
         Scenario = as.factor("RCP 8.5")) %>%
group_by(Year, Month, variable) %>%
   dplyr::summarize(MonthlyTotal = sum(value))
Outflow_USJ2
  

Outflow_USJ2$variable <- factor(Outflow_USJ2$variable, levels = c("Livneh (Historical)", "HadGEM2-ES", "CanESM2", "MIROC5", "CNRM-CM5")) 
Outflow_USJ2

Outflow_USJ3 <- Outflow_USJ2 %>%
  group_by(Year, variable) %>%
   dplyr::summarize(AnnualTotal = sum(MonthlyTotal)) 

Outflow_USJ3
```

```{r}
ggplot(Outflow_USJ2) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_boxplot(aes( x = Month, y = MonthlyTotal, fill = variable), outlier.alpha = 0.3) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
#    scale_y_continuous(limits = c(0, NA),
#                      expand = c(0, NA),
#                     labels = comma) +
     scale_x_discrete(limits=c("Oct","Nov","Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +
   scale_y_log10(breaks=c(0,10,100,1000), labels=comma)+

  labs(fill = "RCP 8.5",
       title = "Upper San Joaquin River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/month)") + #name of x axis
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -17),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
  png("USJ_outflow_GCMs_logscale.png", units ="in", width=8, height=5, res = 300)

```

```{r}
ggplot(Outflow_USJ2) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_boxplot(aes( x = Month, y = MonthlyTotal, fill = variable), outlier.alpha = 0.3) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
    scale_y_continuous(limits = c(0, NA),
                      expand = c(0, NA),
                     labels = comma) +
     scale_x_discrete(limits=c("Oct","Nov","Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +
#   scale_y_log10(breaks=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000), labels=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000))+

  labs(fill = "RCP 8.5",
       title = "Upper San Joaquin River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/month)") + #name of x axis
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -17),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
  png("USJ_outflow_GCMs.png", units ="in", width=8, height=5, res = 300)

```

```{r}
#annual aggreagate
ggplot(Outflow_USJ3) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_violin(draw_quantiles = c(.25, .5, .75), aes( x = variable, y = AnnualTotal, fill = variable)) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
    scale_y_continuous(limits = c(0, NA),
                      expand = c(0, NA),
                      labels = comma) +
#   scale_y_log10(breaks=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000), labels=c(-10000000,-1000000,-100000,-10000, -1000, -100, -10, 0.1,1,10,100,1000))+
    scale_x_discrete(limits=c("Livneh (Historical)", "HadGEM2-ES", "CanESM2", "MIROC5", "CNRM-CM5"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +

  labs(subtitle = "RCP 8.5 Scenario",
       title = "Upper San Joaquin River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/year)") + #name of x axis
  theme(legend.position = "none",
        #legend.direction = "horizontal",
        #legend.box.margin = margin(t = -19),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
png("USJ_outflow_GCMs2.png", units ="in", width=8, height=5, res = 300)
```

```{r}
#Facet wraps of all basins

Outflow_Mer2$Basin <- "Merced River"
Outflow_Mer3$Basin <- "Merced River"
Outflow_Tuo2$Basin <- "Tuolumne River"
Outflow_Tuo3$Basin <- "Tuolumne River"
Outflow_Stn2$Basin <- "Stanislaus River"
Outflow_Stn3$Basin <- "Stanislaus River"
Outflow_USJ2$Basin <- "Upper San Joaquin River"
Outflow_USJ3$Basin <- "Upper San Joaquin River"


all_Outflow_monthly <- rbind(Outflow_Mer2, Outflow_Tuo2, Outflow_Stn2 ,Outflow_USJ2)

all_Outflow_monthly$Basin <- factor(all_Outflow_monthly$Basin, levels = c("Stanislaus River", "Tuolumne River", "Merced River", "Upper San Joaquin River"))
all_Outflow_monthly

all_Outflow_annual <- rbind(Outflow_Mer3, Outflow_Tuo3, Outflow_Stn3 ,Outflow_USJ3)

all_Outflow_annual$Basin <- factor(all_Outflow_annual$Basin, levels = c("Stanislaus River", "Tuolumne River", "Merced River", "Upper San Joaquin River"))
all_Outflow_annual

```

```{r}
all_Outflow_monthly$Month <- factor(all_Outflow_monthly$Month, levels = c("Oct", "Nov", "Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep"))
all_Outflow_monthly


All_Montlhy_Stats_Outflow <- aggregate(MonthlyTotal ~ Month+variable+Basin, all_Outflow_monthly, summary)
write.csv(All_Montlhy_Stats_Outflow, "All_Monthly_Stats_Outflow_mcm_RCP8.5.csv")

All_Annual_Stats_Outflow <- aggregate(AnnualTotal ~ variable+Basin, all_Outflow_annual, summary)
write.csv(All_Annual_Stats_Outflow, "All_Annual_Stats_Outflow_mcm_RCP8.5.csv")

all_Outflow_annual
```

```{r}
ggplot(all_Outflow_monthly) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_boxplot(aes( x = Month, y = MonthlyTotal, fill = variable), outlier.alpha = 0.3) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
    #scale_y_log10(breaks=c(0, 0.1,1,10,100,1000), labels=scales::comma)+
    scale_y_continuous(limits = c(0, NA),
                     expand = c(0, NA),
                 labels = scales::comma) +
    scale_x_discrete(limits=c("Oct","Nov","Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +

  labs(fill = "RCP 8.5",
       #title = "Upper San Joaquin River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/month)") + #name of x axis
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -17),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        panel.spacing.y = unit(0, "lines"))+
  facet_wrap(~Basin, scales = "free_y")+
png("AllBasins_outflow_GCMs_RCP8.5.png", units ="in", width=8, height=5, res = 300)
```

```{r}
#annual aggreagate
ggplot(all_Outflow_annual) + 
  theme_bw(base_size=12, base_family='Times New Roman') + #change font to Times New Roman, 12pt, Bold
    geom_violin(draw_quantiles = c(.25, .5, .75), aes( x = variable, y = AnnualTotal, fill = variable)) + #plot monthly observed data in greenish blue
   #geom_line(aes(x = node, y = value)) +
    scale_y_continuous(limits = c(0, NA),
                      expand = c(0, NA),
                      labels = scales::comma) +
    scale_x_discrete(limits=c("Livneh (Historical)", "HadGEM2-ES", "CanESM2", "MIROC5", "CNRM-CM5"))+
  scale_fill_manual(values = c("#5A5A5A", "#E69F00", "#F0E442", "#009E73", "#56B4E9")) +

  labs(fill = "RCP 8.5",
       #title = "Upper San Joaquin River",
    x = element_blank(), 
       y = "Total River Outflow (mcm/year)") + #name of x axis
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(t = -17),
         axis.ticks.x=element_blank(),
        axis.text.x=element_blank(),
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        strip.placement = "outside",
        strip.background = element_blank(),
        legend.text = element_text(size = 9),
       legend.key.size = unit(0.75,"line"),
        panel.spacing.y = unit(0, "lines"))+
    facet_wrap(~Basin, scales = "free_y")+
png("AllBasins_outflow_GCMs2_RCP8.5.png", units ="in", width=8, height=5, res = 300)
```