In samdchamberlain/salinity_study: Reproducible manuscript evaluating salinity impacts to Sacramento-San Joaquin Delta restored wetlands

^1^Department of Environmental Science, Policy, and Management, University of California, Berkeley, California, USA

#setwd to main package directory
knitr::opts_knit$set(root.dir = normalizePath("../"))

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>")

knitr::opts_chunk$set(echo = FALSE, warning=FALSE)

library(knitr)
library(tidyverse)
library(lubridate)
library(gridExtra)
library(scales)
library(zoo)
library(ggpubr)
library(entropy)
library(forecast)
library(stringr)
source('R/TR_functions.R'); source('R/MI_functions.R'); 
source('R/signficance_tests.R'); source('R/transformations.R')

Results

#load and process wetland sites into half-hourly, daily, and annual sums
source("R/peat6_processing.R") #wetland experiencing salinity rise
source("R/peat19_processing.R") #reference wetland w/o salinity rise

#create a growing month stamp
peat6_all$month <- month(peat6_all$datetime)
peat6_all$m_yr <- paste(peat6_all$year, peat6_all$month, sep = "")

#shuffle data for comparison w/ monthly gap magnitude retained
peat6_shuffled <- peat6_all %>%
  group_by(month) %>%
  mutate(Ta_shuf = sample(TA.x, length(TA.x), replace = F),
         wm_shuf = sample(wm, length(wm), replace = F))

#create data list by month
by_month <- peat6_shuffled %>%
  filter(year > 2011 & year < 2018) %>%
  group_by(m_yr) %>%
  nest()

clipped_monthly <- subset(peat6_monthly, year > 2011 & year < 2018)
clipped_monthly$T_shuf <- mi_timeseries(Ta_shuf, wm_shuf, data_list = by_month)

ggplot(clipped_monthly, aes(x = CH4_obs, y = T_shuf)) +
  stat_smooth(se = F) +
  geom_point(color = "red") +
  ylab(expression(bolditalic(MI)~"("~T[air]~","~F[CH4]~")")) +
  xlab("Total methane flux observations")

source("R/MI_functions.R")
source("R/signficance_tests.R")

## Processing for high salinity wetland

#create data list by month
by_month <- peat6_all %>%
  filter(year > 2011 & year < 2018) %>%
  group_by(m_yr) %>%
  nest()

#plot that needs to be dealt with later
clipped_monthly <- subset(peat6_monthly, year > 2011 & year < 2018)
clipped_monthly$date <- as.Date(clipped_monthly$datetime)
clipped_monthly$GPP_CH4_pear <- corr_timeseries(gpp_ANNnight*-1, wm, data_list = by_month, 
                                              method = 'pearson')
clipped_monthly$Tair_CH4_pear <- corr_timeseries(TA.x, wm, data_list = by_month, 
                                               method = 'pearson')

clipped_monthly$GPP_CH4_kendall <- corr_timeseries(gpp_ANNnight*-1, wm, data_list = by_month, 
                                              method = 'kendall')
clipped_monthly$Tair_CH4_kendall <- corr_timeseries(TA.x, wm, data_list = by_month, 
                                               method = 'kendall')

a <- ggplot(clipped_monthly) +
  geom_point(aes(x = date, y = Tair_CH4_pear), color = "red", alpha = 0.5) +
  geom_line(aes(x = date, y = ma(Tair_CH4_pear, 3)), color = "red", size = 1) +
  geom_point(aes(x = date, y = GPP_CH4_pear), color = "dark green", alpha = 0.5) +
  geom_line(aes(x = date, y = ma(GPP_CH4_pear, 3)), color = "dark green", size = 1) +  
  scale_x_date(date_breaks = "1 year", labels = date_format("%Y")) +
  ylab(expression("PCC (x,"~F[CH4]~")")) + xlab("") + ggtitle("Salinization wetland") +
  theme_bw()

b <- ggplot(clipped_monthly) +
  geom_point(aes(x = date, y = Tair_CH4_kendall), color = "red", alpha = 0.5) +
  geom_line(aes(x = date, y = ma(Tair_CH4_kendall, 3)), color = "red", size = 1) +
  geom_point(aes(x = date, y = GPP_CH4_kendall), color = "dark green", alpha = 0.5) +
  geom_line(aes(x = date, y = ma(GPP_CH4_kendall, 3)), color = "dark green", size = 1) +  
  scale_x_date(date_breaks = "1 year", labels = date_format("%Y")) +
  ylab(expression(tau~"(x,"~F[CH4]~")")) + xlab("") +
  theme_bw()

c <- ggplot(subset(peat6_daily, year > 2011 & year < 2018), aes(x=as.Date(datetime), y=Sal)) +
  geom_line(color = "blue") +
  scale_y_continuous(limits = c(0, 7)) +
  scale_x_date(date_breaks = "1 year", labels = date_format("%Y")) +
  ylab("Salinity (psu)") + xlab("Year") +
  theme_bw()

#Multiple plots aligned
iA <- ggplotGrob(a)
iB <- ggplotGrob(b)
iC <- ggplotGrob(c)

maxWidth = grid::unit.pmax(iA$widths[2:5], iB$widths[2:5], iC$widths[2:5])
iA$widths[2:5] <- as.list(maxWidth); iB$widths[2:5] <- as.list(maxWidth); 
iC$widths[2:5] <- as.list(maxWidth); 

grid.arrange(iA, iB, iC, ncol=1)

time_series <- theme_bw() +
  theme(panel.grid.minor=element_blank(), panel.grid.major.y = element_blank(),
        axis.title.x = element_blank(), axis.title.y = element_text(size=9))

both_sites <- subset(rbind(peat19_daily, peat6_daily), year > 2011)
both_sites$date <- as.Date(both_sites$datetime)

a <- ggplot(both_sites, aes(x=date, y=mgCH4, color=site)) +
  geom_line(size = 0.5) + 
  scale_x_date(date_breaks = "1 year", labels = date_format("%Y")) +
  ylab(expression(F[CH4]~" "~"("~mg~C~m^{-2}~d^{-1}~")")) + 
  time_series +
  theme(legend.title=element_blank(), legend.margin = margin(0, 0, 0, 0),
        legend.position=c(.85, .80),
        legend.key.size =  unit(0.2, "in"), legend.text = element_text(size=8),
        legend.background = element_rect(fill = "transparent"))

b <- ggplot(both_sites, aes(x=date, y=gGEP, color=site)) +
  geom_line(size = 0.5) +
  scale_x_date(date_breaks = "1 year", labels = date_format("%Y")) +
  ylab(expression("GEP "~"("~g~C~m^{-2}~d^{-1}~")")) +
  time_series + theme(legend.position = "none")

c <- ggplot(both_sites, aes(x=date, y=gCO2, color=site)) +
  geom_line(size = 0.5) +
  ylab(expression("NEE "~"("~g~C~m^{-2}~d^{-1}~")")) +
  scale_x_date(date_breaks = "1 year", labels = date_format("%Y")) +
  time_series + theme(legend.position = "none")

d <- ggplot(both_sites, aes(x=date, y=Sal, color=site)) +
  geom_line() + 
  ylab("Salinity (psu)") +
  scale_x_date(date_breaks = "1 year", labels = date_format("%Y")) +
  time_series + theme(legend.position = "none")

#Multiple plots aligned
iA <- ggplotGrob(a)
iB <- ggplotGrob(b)
iC <- ggplotGrob(c)
iD <- ggplotGrob(d)


maxWidth = grid::unit.pmax(iA$widths[2:5], iB$widths[2:5], iC$widths[2:5], iD$widths[2:5])
iA$widths[2:5] <- as.list(maxWidth); iB$widths[2:5] <- as.list(maxWidth); 
iC$widths[2:5] <- as.list(maxWidth); iD$widths[2:5] <- as.list(maxWidth); 

grid.arrange(iA, iB, iC, iD, ncol=1)

#apply anomaly filter
anom_data <- peat6_all %>%
  arrange(year, time) %>%
  mutate(CH4_anom = anom_filter(wm_gf),
         TA_anom = anom_filter(TA.x),
         GPP_anom = anom_filter(gpp_ANNnight))

anom_data <- arrange(anom_data, datetime)
anom_data$CH4na_anom <- ifelse(is.na(anom_data$wm), NA, anom_data$CH4_anom) # re-insert missing values

# Focus only on high v. low salinity years without additional disturbance
# other disturbance includes the 2014 insect outbreak and 2017 flushing events
anom_data <- subset(anom_data, year %in% c(2012, 2013, 2015, 2016))

#create data list by annual growing season
by_year <- anom_data %>%
  filter(DOY > 100 & DOY < 300) %>%
  group_by(year) %>%
  nest()

#calculate lagged transfer entropy time series (up to 2 day lag at half-hourly timestep)
# for TA -> CH4
tr_wm_ta <- tr_timeseries(TA_anom, CH4na_anom, data_list = by_year, type = "observed", 
                           lags = 120, bins = 10, normalize = T)

mc_wm_ta <- tr_timeseries(TA_anom, CH4na_anom, data_list = by_year, type = "shuffled", 
                           lags = 120, bins = 10, runs = 100, alpha = 0.05, normalize = T)
# for CH4 -> TA
tr_ta_wm <- tr_timeseries(CH4na_anom, TA_anom, data_list = by_year, type = "observed", 
                           lags = 120, bins = 10, normalize = T)

mc_ta_wm <- tr_timeseries(CH4na_anom, TA_anom, data_list = by_year, type = "shuffled", 
                           lags = 120, bins = 10, runs = 100, alpha = 0.05, normalize = T)

#unlist data to tidy dataframe
wm.ta_df <- tidy_list(tr_ta_wm); wm.ta_shuf <- tidy_list(mc_ta_wm)
ta.wm_df <- tidy_list(tr_wm_ta); ta.wm_shuf <- tidy_list(mc_wm_ta)

a <- ggplot() +
  geom_line(data = ta.wm_df, aes(x = hour, y = tr)) +
  geom_line(data = ta.wm_shuf, aes(x = hour, y = tr), color = "red") +
  facet_wrap(~year, nrow = 1) +
  scale_x_continuous(breaks = seq(0, 60, by = 12)) +
  scale_y_continuous(limits = c(0, 0.043)) +
  ylab(expression(bolditalic(TE)~"("~T[a] %->% F[CH4]~")"))

b <- ggplot() +
  geom_line(data = wm.ta_df, aes(x = hour, y = tr)) +
  geom_line(data = wm.ta_shuf, aes(x = hour, y = tr), color = "red") +
  facet_wrap(~year, nrow = 1) +
  scale_x_continuous(breaks = seq(0, 60, by = 12)) +
  scale_y_continuous(limits = c(0, 0.043)) +
  xlab("Lag (hours)")  + ylab(expression(bolditalic(TE)~"("~F[CH4] %->% T[a]~")"))

#Multiple plots aligned
iA <- ggplotGrob(a)
iB <- ggplotGrob(b)

maxWidth = grid::unit.pmax(iA$widths[2:5], iB$widths[2:5])
iA$widths[2:5] <- as.list(maxWidth); iB$widths[2:5] <- as.list(maxWidth); 

grid.arrange(iA, iB, ncol=1)

#calculate lagged transfer entropy time series (up to 2 day lag at half-hourly timestep)
# for TA -> CH4
tr_wm_gep <- tr_timeseries(GPP_anom, CH4na_anom, data_list = by_year, type = "observed", 
                           lags = 120, bins = 10, normalize = T)

mc_wm_gep <- tr_timeseries(GPP_anom, CH4na_anom, data_list = by_year, type = "shuffled", 
                           lags = 120, bins = 10, runs = 100, alpha = 0.05, normalize = T)
# for CH4 -> TA
tr_gep_wm <- tr_timeseries(CH4na_anom, GPP_anom, data_list = by_year, type = "observed", 
                           lags = 120, bins = 10, normalize = T)

mc_gep_wm <- tr_timeseries(CH4na_anom, GPP_anom, data_list = by_year, type = "shuffled", 
                           lags = 120, bins = 10, runs = 100, alpha = 0.05, normalize = T)

#unlist data to tidy dataframe
wm.gep_df <- tidy_list(tr_gep_wm); wm.gep_shuf <- tidy_list(mc_gep_wm)
gep.wm_df <- tidy_list(tr_wm_gep); gep.wm_shuf <- tidy_list(mc_wm_gep)

a <- ggplot() +
  geom_line(data = gep.wm_df, aes(x = hour, y = tr)) +
  geom_line(data = gep.wm_shuf, aes(x = hour, y = tr), color = "red") +
  facet_wrap(~year, nrow = 1) +
  scale_x_continuous(breaks = seq(0, 60, by = 12)) +
  scale_y_continuous(limits = c(0, 0.043)) +
  ylab(expression(bolditalic(TE)~"("~GEP %->% F[CH4]~")"))

b <- ggplot() +
  geom_line(data = wm.gep_df, aes(x = hour, y = tr)) +
  geom_line(data = wm.gep_shuf, aes(x = hour, y = tr), color = "red") +
  facet_wrap(~year, nrow = 1) +
  scale_x_continuous(breaks = seq(0, 60, by = 12)) +
  scale_y_continuous(limits = c(0, 0.043)) +
  xlab("Lag (hours)")  + ylab(expression(bolditalic(TE)~"("~F[CH4] %->% GEP~")"))

#Multiple plots aligned
iA <- ggplotGrob(a)
iB <- ggplotGrob(b)

maxWidth = grid::unit.pmax(iA$widths[2:5], iB$widths[2:5])
iA$widths[2:5] <- as.list(maxWidth); iB$widths[2:5] <- as.list(maxWidth); 

grid.arrange(iA, iB, ncol=1)

samdchamberlain/salinity_study documentation built on May 4, 2019, 1:23 p.m.