---
title: "Vulnerability of Estuarine Systems in the Contiguous United States to Water Quality Change Under Future Climate and Land-Use"
output: html_document
author: "Lise Montefiore"
date: "2022-10-19"

mainfont: Times New Roman
fontsize: 12pt
header-includes:
  - \usepackage{titling}
  - \pretitle{\begin{flushleft}\small\bfseries}
  - \posttitle{\end{flushleft}}  
  - \preauthor{\begin{flushleft}\Large}
  - \postauthor{\end{flushleft}}  
  - \predate{\begin{flushleft}\large}
  - \postdate{\end{flushleft}}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



 <br>
 
 
 
 
 
#### 1.  Overview

A vulnerability framework was developed to identify estuarine systems at higher risk of water quality change under future climate and land-use. The vulnerability was developed for 112 estuarine systems across the contiguous US and includes several factors capturing biological, physical, ecological, and social aspects of estuarine vulnerability to water quality change. for the analysis,historical and projected TN and TP yields were computed using STEPLgrid model (Montefiore & Nelson, 2022) under the climate scenarios RCPs 4.5 and 8.5 and the land-use land-cover (LULC) scenarios A1B, A2, B1, and B2. 

This document aims to reproduce Figures 3, 5, and 7 in the main manuscript the manuscript and Figure S6 in the supporting information.

#### 2.  Load the R packages and set your working directory 

```{r Load R Packages,message=FALSE, warning=FALSE}
rm(list=ls())
library(dplyr)
library(ggplot2)
library(patchwork)
library(reshape)
library(wesanderson)
setwd("/Volumes/Seagate Backup Plus Drive/chapter_3_share/manuscript/Dryad/")
```

#### 3.  Annual or seasonal TN and TP yields from the averaging ensemble climate models

```{r Annual or Seasonal projected TN and TP yields,warning=FALSE,message=FALSE}
nut<- read.csv("nutrient_yields.csv")

### Create TN plot
tn<- nut %>% filter(nutrient=="TN")
levels(tn$period)[levels(tn$period)=="fall"] <- "Fall"
levels(tn$period)[levels(tn$period)=="spring"] <- "Spring"
levels(tn$period)[levels(tn$period)=="summer"] <- "Summer"
levels(tn$period)[levels(tn$period)=="winter"] <- "Winter"
levels(tn$land_use)[levels(tn$land_use)=="NLCD 2016"] <- "NLCD"

tn[tn == "annual"] <- "Annual"
tn[tn == "fall"] <- "Fall"
tn[tn == "winter"] <- "Winter"
tn[tn == "spring"] <- "Spring"
tn[tn == "summer"] <- "Summer"
tn[tn == "NLCD 2016"] <- "NLCD"

tn2<- tn
tn2$land_use <- factor(tn2$land_use,
                 levels = c('NLCD','A1B', "A2", "B1", "B2"),ordered = TRUE)

## Filter annual period for the annual plot 
tn2<- tn2 %>% filter(period %in% "Annual") 
p1<- ggplot(tn2, aes(x=land_use, y=Yield, fill=climate)) + 
  geom_boxplot(outlier.size=0.1) +
  theme_bw(base_size = 10) +
  facet_grid( period ~ region , scales="free")+
  scale_fill_manual(values=c("#ED936B", "#7DC0A6", "#909FC7"), 
                    name = "", labels = c( "Historical",  "RCP 4.5", "RCP 8.5"))+
  labs(y = bquote('TN Yield'~(kg/km^2)), x="LULC")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme(legend.position = "none") +ggtitle("(a) TN")

######Create TP plot ############
tp<- nut %>% filter(nutrient=="TP")
tp[tp == "annual"] <- "Annual"
tp[tp == "fall"] <- "Fall"
tp[tp == "winter"] <- "Winter"
tp[tp == "spring"] <- "Spring"
tp[tp == "summer"] <- "Summer"
tp[tp == "NLCD 2016"] <- "NLCD"

tp2<- tp
tp2$land_use <- factor(tp2$land_use,
                       levels = c('NLCD','A1B', "A2", "B1", "B2"),ordered = TRUE)

## Filter annual period for the annual plot 
tp2<- tp2 %>% filter(period %in% "Annual") 

p2<- ggplot(tp2, aes(x=land_use, y=Yield, fill=climate)) + 
  geom_boxplot(outlier.size=0.1) +
  theme_bw(base_size = 10) + 
  facet_grid( period ~ region  , scales="free")+
  scale_fill_manual(values=c("#ED936B", "#7DC0A6", "#909FC7"), 
                    name = "", labels = c( "Historical",  "RCP 4.5", "RCP 8.5"))+
  labs(y = bquote('TP Yield'~(kg/km^2)), x="LULC")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  guides(fill=guide_legend(title="Climate Scenario")) +ggtitle("(b) TP")

p3<- p1 / p2
p3
# save the figure in your directory 
#ggsave(device = "png", filename = "annual_nutrient_yields", plot = p3, width=15, heigh=7.5)
```

#### 4.  Projected annual yields for the 20 climate ensemble models

The 10th  and 90th percentiles of the nutrient yields were computed using the 20 climate models outputs and were compared to the outputs corresponding to the multi-model ensemble mean. 

```{r Nutrient yields variability, pvf, fig.height = 8, fig.width = 10,warning=FALSE,message=FALSE}
g<- read.csv("nutrient_yields_variability.csv")
v1<- g%>% filter(nutrient=="TN") %>% filter(climate=="RCP45")
p1<- ggplot(v1, aes(x=land_use, y=value, fill=variable)) + 
  geom_boxplot() + facet_wrap(~region, scales="free")+
  theme_bw(base_size = 10) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(y = bquote('TN Yield'~(kg/km^2)), x="LULC")+
  scale_fill_manual(values=c( "#DEA25B", "#75BDAD", "#E1DDBC"), name = "")+
  theme(legend.position="none")+
  ggtitle("(a) TN-RCP 4.5")

v2<- g%>% filter(nutrient=="TN") %>% filter(climate=="RCP85")
p2<- ggplot(v2, aes(x=land_use, y=value, fill=variable)) + 
  geom_boxplot() + facet_wrap(~region, scales="free")+
  theme_bw(base_size = 10) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(y = bquote('TN Yield'~(kg/km^2)), x="LULC")+
  scale_fill_manual(values=c( "#DEA25B", "#75BDAD", "#E1DDBC"), 
                    name = "", labels = c(bquote(10^th ~ percentile), bquote(90^th ~ percentile), "Multi-model mean"))+
  ggtitle("(b) TN-RCP 8.5")


v3<- g%>% filter(nutrient=="TP") %>% filter(climate=="RCP45")
p3<- ggplot(v3, aes(x=land_use, y=value, fill=variable)) + 
  geom_boxplot() + facet_wrap(~region, scales="free")+
  theme_bw(base_size = 10) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(y = bquote('TP Yield'~(kg/km^2)), x="LULC")+
  scale_fill_manual(values=c( "#DEA25B", "#75BDAD", "#E1DDBC"), 
                    name = "", labels = c("Minimum", "Maximum", "Multi-model mean"))+
  theme(legend.position="none")+
  ggtitle("(c) TP-RCP 4.5")


v4<- g%>% filter(nutrient=="TP") %>% filter(climate=="RCP85")
p4<- ggplot(v4, aes(x=land_use, y=value, fill=variable)) + 
  geom_boxplot() + facet_wrap(~region, scales="free")+
  theme_bw(base_size = 10) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  labs(y = bquote('TP Yield'~(kg/km^2)), x="LULC")+
  scale_fill_manual(values=c( "#DEA25B", "#75BDAD", "#E1DDBC"), 
                    name = "", labels = c("Minimum", "Maximum", "Multi-model mean"))+
  theme(legend.position="none")+
  ggtitle("(d) TP-RCP 8.5")

### Combine the code ########
pvf<- (p1 + p2) /(p3 +p4)
pvf
# save the figure in your directory 
#ggsave(device = "png", filename = "climate_variability", plot = pvf, width=20, heigh=12)
```

#### 5.  Nutrient exposure scores corresponding to each of the considered climate models (n=20) under the climate and LULC scenarios.

Cells in the heatmaps are colored by the number of times the exposure scores fell within the given ranges when calculated for each climate model output. The exposure score computed with the multi-model ensemble mean for each watershed is indicated by the symbol “+”. 

```{r Exposure score variability, p1, fig.height = 15, fig.width = 10,warning=FALSE,message=FALSE}
##### Plot exposure variability for TN under RCP4.5 and 8.5
dd<- read.csv("exposure_variability.csv")
### Choice nutrient "TN" or "TP"
nutrient<- "TN"
### Choice scenarios "rcp45" or "rcp85"
rcp<- "rcp45"
### Choice LULC scenarios "A1B" or "A2", "B1", "B2"
lulc<- "A1B"
### color palette
pal <- wes_palette("Zissou1", 20, type = "continuous")

### Filter and create the plot
dd1<- dd %>% filter(nutrient==nutrient) %>%
             filter(rcp==rcp) %>%
            filter(lulc==lulc)

p1<- ggplot(dd1, aes(x = variable, y = name, fill = value)) +
  geom_tile() + ylab("") + xlab("") + 
  geom_text(aes(x = variable, y = name, label=MM_label))+ 
  scale_fill_gradientn(na.value="#F8F8F8", name="Number of climate models",colours = pal)+ 
  theme(legend.position="none")+
  theme(axis.text.x = element_text(angle = 45, hjust=1), 
        axis.text=element_text(size=10), axis.title=element_text(size=10),
        plot.title = element_text(size = 10)) +
  ggtitle(paste("(a) ", lulc,"-" ,rcp, nutrient,  sep="")) 
p1
```

#### 6.  Historical and Projected annual average event precipitation (mm)

Annual average event precipitation (mm) for the historical (1990-2020; green lines) and future (2035-2065) periods under RCP 4.5 (red line) and 8.5 (blue line) scenarios across the study regions. The minimum and maximum outputs of the climate model ensemble are shown as red and blue shaded areas for the RCP 4.5 and 8.5 scenarios, respectively.

``` {r Climate Event, warning=FALSE, message=FALSE}
dc<- read.csv("average_event_ precipitation.csv")

pc<- ggplot(dc) + 
  geom_ribbon(aes(x = year, y = Mean, 
                  ymin = Min, ymax = Max, fill=climate), alpha=0.2)+
  geom_line(aes(x = year, y = Mean, color=climate)) + 
  facet_wrap(~region)+ 
  labs(y = "Annual Average Precipitation Event", x="")+
  scale_fill_manual(values=c("#66c2a5", "#fc8d62", "#8da0cb"), name="")+
  scale_color_manual(values=c("#66c2a5", "#fc8d62", "#8da0cb"), name="")+
  theme_bw(base_size = 10) +
  theme(panel.grid.minor = element_blank(), panel.grid.major.x = element_blank())
  pc
```
