#### impact of drought analysis ## Reiter et al. ## date: January 2017; mod June 2017; mod December 2017 ## This code is for modeling wetland and agriculture types individually by year set ## It includes code for making plots shown in the manuscript ## required data: "wetlands_merge_scene_dates_ownership.csv"; "err2.csv"; ## "Ag_merge_precip3.csv" setwd("C:/water") wetwater<-read.csv("wetlands_merge_scene_dates_ownership.csv", header=T) ## read in bias correction data for reference bc<-read.csv("C:/water/err2.csv", header=T) library(gamm4) library(ggplot2) library(gridExtra) library(cowplot) #library(gam) ## add in drought indicator wetwater$drought<-ifelse(wetwater$water.year%in%c(2001:2002,2007:2015),1,0) ## only 2000-2011 data wetwater1<- subset(wetwater,wetwater$year<2012) ## drought year 2000-2011 wetwater1d<-subset(wetwater1, wetwater1$drought==1) ## non-drought year 2000-2011 wetwater1nd<-subset(wetwater1, wetwater1$drought==0) ## 2013-2015 data - extreme drought wetwater2<-subset(wetwater, wetwater$year>2012) ## ####### grab only seasonal wetlands wetwater1seas<-subset(wetwater1,wetwater1$landcover==1 & wetwater1$prop.sampled>=0.5|wetwater1$landcover==10 & wetwater1$prop.sampled>=0.5) wetwater1seasD<-subset(wetwater1d,wetwater1d$landcover==1 & wetwater1d$prop.sampled>=0.5|wetwater1d$landcover==10 & wetwater1d$prop.sampled>=0.5) wetwater1seasND<-subset(wetwater1nd,wetwater1nd$landcover==1 & wetwater1nd$prop.sampled>=0.5|wetwater1nd$landcover==10 & wetwater1nd$prop.sampled>=0.5) ## standardize some covariates wetwater1seas$zday = (wetwater1seas$yday-150)/100 wetwater1seas$id = c(1:nrow(wetwater1seas)) wetwater1seas$sprecip4 = wetwater1seas$precip.4wk/100 wetwater1seas$sprecip2 = wetwater1seas$precip.2wk/100 ## correct for classification bias between sensors wetwater1seas$nfloodbc<-ifelse(wetwater1seas$year<2012, round(wetwater1seas$n.flooded-wetwater1seas$n.flooded*-.0127), round(wetwater1seas$n.flooded-wetwater1seas$n.flooded*-.1065)) ## some issueswith bias correction where the corrected number of flooded pixels is greater than the number sampled ## line below makes the correctd value equal to the number sampled. wetwater1seas$nfloodbc<-ifelse(wetwater1seas$nfloodbc>wetwater1seas$n.sampled,wetwater1seas$n.sampled,wetwater1seas$nfloodbc) ### now looking at modeling only the drought years 2000-2011 wetwater1seasD$zday = (wetwater1seasD$yday-150)/100 wetwater1seasD$id = c(1:nrow(wetwater1seasD)) wetwater1seasD$sprecip4 = wetwater1seasD$precip.4wk/100 wetwater1seasD$sprecip2 = wetwater1seasD$precip.2wk/100 ## correct for classification bias between sensors wetwater1seasD$nfloodbc<-ifelse(wetwater1seasD$year<2012, round(wetwater1seasD$n.flooded-wetwater1seasD$n.flooded*-0.0127), round(wetwater1seasD$n.flooded-wetwater1seasD$n.flooded*-.1065)) f1D<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=wetwater1seasD, weights=wetwater1seasD$prop.sampled) newdatf1D = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) newdatf1D$zday = (newdatf1D$yday-150)/100 newdatf1D$prop.flooded = predict(f1D$gam, newdatf1D, type='response') predD = as.data.frame(predict(f1D$gam, newdatf1D, type='link', se.fit=T)) newdatf1D$lcl = plogis(predD$fit-2*predD$se.fit) newdatf1D$ucl = plogis(predD$fit+2*predD$se.fit) ## the non-drought years 2000-2011 wetwater1seasND$zday = (wetwater1seasND$yday-150)/100 wetwater1seasND$id = c(1:nrow(wetwater1seasND)) wetwater1seasND$sprecip4 = wetwater1seasND$precip.4wk/100 wetwater1seasND$sprecip2 = wetwater1seasND$precip.2wk/100 ## correct for classification bias between sensors wetwater1seasND$nfloodbc<-ifelse(wetwater1seasND$year<2012, round(wetwater1seasND$n.flooded-wetwater1seasND$n.flooded*-0.0127), round(wetwater1seasND$n.flooded-wetwater1seasND$n.flooded*-0.1065)) ## some issueswith bias correction where the corrected number of flooded pixels is greater than the number sampled ## line below makes the correctd value equal to the number sampled. wetwater1seasND$nfloodbc<-ifelse(wetwater1seasND$nfloodbc>wetwater1seasND$n.sampled,wetwater1seasND$n.sampled,wetwater1seasND$nfloodbc) f1ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=wetwater1seasND, weights=wetwater1seasND$prop.sampled) ## data for plots from the fitted model newdatf1ND = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) newdatf1ND$zday = (newdatf1ND$yday-150)/100 newdatf1ND$prop.flooded = predict(f1ND$gam, newdatf1ND, type='response') predND = as.data.frame(predict(f1ND$gam, newdatf1ND, type='link', se.fit=T)) newdatf1ND$lcl = plogis(predND$fit-2*predND$se.fit) newdatf1ND$ucl = plogis(predND$fit+2*predND$se.fit) ############################################ ## fitting gam to 2013-2015 data wetwater2$zday = (wetwater2$yday-150)/100 wetwater2$id = c(1:nrow(wetwater2)) wetwater2$sprecip4 = wetwater2$precip.4wk/100 wetwater2$sprecip2 = wetwater2$precip.2wk/100 ## correct for classification bias between sensors wetwater2$nfloodbc<-ifelse(wetwater2$year<2012, round(wetwater2$n.flooded-wetwater2$n.flooded*-.0127), round(wetwater2$n.flooded-wetwater2$n.flooded*-.1065)) ## some issueswith bias correction where the corrected number of flooded pixels is greater than the number sampled ## line below makes the correctd value equal to the number sampled. wetwater2$nfloodbc<-ifelse(wetwater2$nfloodbc>wetwater2$n.sampled,wetwater2$n.sampled,wetwater2$nfloodbc) wetwater2seas<-subset(wetwater2,wetwater2$landcover==1&wetwater2$prop.sampled>=0.5|wetwater2$landcover==10&wetwater2$prop.sampled>=0.5) f1r<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=wetwater2seas, weights=wetwater2seas$prop.sampled) newdatf2 = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) newdatf2$zday = (newdatf2$yday-150)/100 newdatf2$prop.flooded = predict(f1r$gam, newdatf2, type='response') pred = as.data.frame(predict(f1r$gam, newdatf2, type='link', se.fit=T)) newdatf2$lcl = plogis(pred$fit-2*pred$se.fit) newdatf2$ucl = plogis(pred$fit+2*pred$se.fit) #########################################3 ## combine all three model plots newdatDNDe<-data.frame(rbind(newdatf1ND, newdatf1D, newdatf2)) newdatDNDe$g<-c(rep("non-drought",319), rep("drought",319), rep("critical",319)) seas<-ggplot(data = newdatDNDe,aes(x = yday, y = prop.flooded, ymin=lcl, ymax=ucl,fill = g)) + geom_ribbon(alpha = 0.4) + scale_fill_manual(name = "", values = c('black','green',"blue"))+ #labels = c('2013-2015','Drought: 2000-2011','Non-Drought: 2000-2011'))+ scale_x_continuous(limits=c(1,319), expand=c(0,0), breaks=c(1,32,62,93,124,154,185,216,244,275,305), labels=c('Jul','Aug','Sep','Oct','Nov','Dec','Jan','Feb','Mar','Apr','May'))+ scale_y_continuous(labels=scaleFUN)+ xlab(NULL) + ylab('Proportion Open Water')+ theme(legend.position="none")+ ggtitle("Seasonal") print(seas) ####################################################################### ## analysis for semi-permanent wetlands ####### grab only semional wetlands wetwater1semi<-subset(wetwater1,wetwater1$landcover==2 & wetwater1$prop.sampled>=0.5|wetwater1$landcover==20 & wetwater1$prop.sampled>=0.5) wetwater1semiD<-subset(wetwater1d,wetwater1d$landcover==2 & wetwater1d$prop.sampled>=0.5|wetwater1d$landcover==20 & wetwater1d$prop.sampled>=0.5) wetwater1semiND<-subset(wetwater1nd,wetwater1nd$landcover==2 & wetwater1nd$prop.sampled>=0.5|wetwater1nd$landcover==20 & wetwater1nd$prop.sampled>=0.5) ## standardize some covariates wetwater1semi$zday = (wetwater1semi$yday-150)/100 wetwater1semi$id = c(1:nrow(wetwater1semi)) wetwater1semi$sprecip4 = wetwater1semi$precip.4wk/100 wetwater1semi$sprecip2 = wetwater1semi$precip.2wk/100 ## correct for classification bias between sensors wetwater1semi$nfloodbc<-ifelse(wetwater1semi$year<2012, round(wetwater1semi$n.flooded-wetwater1semi$n.flooded*-.0127), round(wetwater1semi$n.flooded-wetwater1semi$n.flooded*-.1065)) ## some issueswith bias correction where the corrected number of flooded pixels is greater than the number sampled ## line below makes the correctd value equal to the number sampled. wetwater1semi$nfloodbc<-ifelse(wetwater1semi$nfloodbc>wetwater1semi$n.sampled,wetwater1semi$n.sampled,wetwater1semi$nfloodbc) ### now looking at modeling only the drought years 2000-2011 wetwater1semiD$zday = (wetwater1semiD$yday-150)/100 wetwater1semiD$id = c(1:nrow(wetwater1semiD)) wetwater1semiD$sprecip4 = wetwater1semiD$precip.4wk/100 wetwater1semiD$sprecip2 = wetwater1semiD$precip.2wk/100 ## correct for classification bias between sensors wetwater1semiD$nfloodbc<-ifelse(wetwater1semiD$year<2012, round(wetwater1semiD$n.flooded-wetwater1semiD$n.flooded*-0.0127), round(wetwater1semiD$n.flooded-wetwater1semiD$n.flooded*-.1065)) f1sD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=wetwater1semiD, weights=wetwater1semiD$prop.sampled) newdatf1sD = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) newdatf1sD$zday = (newdatf1sD$yday-150)/100 newdatf1sD$prop.flooded = predict(f1sD$gam, newdatf1sD, type='response') predD = as.data.frame(predict(f1sD$gam, newdatf1sD, type='link', se.fit=T)) newdatf1sD$lcl = plogis(predD$fit-2*predD$se.fit) newdatf1sD$ucl = plogis(predD$fit+2*predD$se.fit) ## the non-drought years 2000-2011 wetwater1semiND$zday = (wetwater1semiND$yday-150)/100 wetwater1semiND$id = c(1:nrow(wetwater1semiND)) wetwater1semiND$sprecip4 = wetwater1semiND$precip.4wk/100 wetwater1semiND$sprecip2 = wetwater1semiND$precip.2wk/100 ## correct for classification bias between sensors wetwater1semiND$nfloodbc<-ifelse(wetwater1semiND$year<2012, round(wetwater1semiND$n.flooded-wetwater1semiND$n.flooded*-0.0127), round(wetwater1semiND$n.flooded-wetwater1semiND$n.flooded*-0.1065)) ## some issueswith bias correction where the corrected number of flooded pixels is greater than the number sampled ## line below makes the correctd value equal to the number sampled. wetwater1semiND$nfloodbc<-ifelse(wetwater1semiND$nfloodbc>wetwater1semiND$n.sampled,wetwater1semiND$n.sampled,wetwater1semiND$nfloodbc) f1sND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=wetwater1semiND, weights=wetwater1semiND$prop.sampled) ## data for plots from the top model newdatf1sND = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) newdatf1sND$zday = (newdatf1sND$yday-150)/100 newdatf1sND$prop.flooded = predict(f1sND$gam, newdatf1sND, type='response') predND = as.data.frame(predict(f1sND$gam, newdatf1sND, type='link', se.fit=T)) newdatf1sND$lcl = plogis(predND$fit-2*predND$se.fit) newdatf1sND$ucl = plogis(predND$fit+2*predND$se.fit) ############################################ ## fitting gam to 2013-2015 data wetwater2$zday = (wetwater2$yday-150)/100 wetwater2$id = c(1:nrow(wetwater2)) wetwater2$sprecip4 = wetwater2$precip.4wk/100 wetwater2$sprecip2 = wetwater2$precip.2wk/100 ## correct for classification bias between sensors wetwater2$nfloodbc<-ifelse(wetwater2$year<2012, round(wetwater2$n.flooded-wetwater2$n.flooded*-.0127), round(wetwater2$n.flooded-wetwater2$n.flooded*-.1065)) ## some issues with bias correction where the corrected number of flooded pixels is greater than the number sampled ## line below makes the correctd value equal to the number sampled. wetwater2$nfloodbc<-ifelse(wetwater2$nfloodbc>wetwater2$n.sampled,wetwater2$n.sampled,wetwater2$nfloodbc) wetwater2semi<-subset(wetwater2,wetwater2$landcover==2&wetwater2$prop.sampled>=0.5|wetwater2$landcover==20&wetwater2$prop.sampled>=0.5) f1sr<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=wetwater2semi, weights=wetwater2semi$prop.sampled) newdatf2 = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) newdatf2$zday = (newdatf2$yday-150)/100 newdatf2$prop.flooded = predict(f1sr$gam, newdatf2, type='response') pred = as.data.frame(predict(f1sr$gam, newdatf2, type='link', se.fit=T)) newdatf2$lcl = plogis(pred$fit-2*pred$se.fit) newdatf2$ucl = plogis(pred$fit+2*pred$se.fit) #########################################3 ## combine all three model plots newdatDNDe<-data.frame(rbind(newdatf1sND, newdatf1sD, newdatf2)) newdatDNDe$g<-c(rep("non-drought",319), rep("drought",319), rep("critical",319)) semi<-ggplot(data = newdatDNDe,aes(x = yday, y = prop.flooded, ymin=lcl, ymax=ucl,fill = g)) + geom_ribbon(alpha = 0.4) + scale_fill_manual(name = "", values = c('black','green',"blue"))+ #labels = c('2013-2015','Drought: 2000-2011','Non-Drought: 2000-2011'))+ scale_x_continuous(limits=c(1,319), expand=c(0,0), breaks=c(1,32,62,93,124,154,185,216,244,275,305), labels=c('Jul','Aug','Sep','Oct','Nov','Dec','Jan','Feb','Mar','Apr','May'))+ scale_y_continuous(labels=scaleFUN)+ xlab(NULL) + ylab('Proportion Open Water')+ theme(legend.position="none")+ ggtitle("Semi-permanent") ####### # combine seasonal and semipermanent into single plot wet1<- plot_grid(seas, semi, nrow=2, ncol=1, align="hv") save_plot("wet_all2.png", wet1, nrow=2, ncol=1) wet2<- plot_grid(seas, semi, nrow=1, ncol=2, align="hv") save_plot("wet_all2b.png", wet2, nrow=1, ncol=2) ########################################################### ### AGRICULTURE MODELS ## separate models by year type ############# setwd("C:/water") agwater<-read.csv("Ag_merge_precip3.csv", header=T) ## read in bias correction data for reference bc<-read.csv("C:/water/err2.csv", header=T) ## ## drought indicator agwater$drought<-ifelse(agwater$water.year%in%c(2001:2002,2007:2015),1,0) agwater$zday = (agwater$yday-150)/100 agwater$id = c(1:nrow(agwater)) ## agwater valley agwaterrice<-subset(agwater,agwater$landcover=="Rice" & agwater$prop.sampled>=0.5) agwatercorn<-subset(agwater,agwater$landcover=="Corn" & agwater$prop.sampled>=0.5) agwaterother<-subset(agwater,agwater$landcover=="Field Crop" & agwater$prop.sampled>=0.5| agwater$landcover=="Grains" & agwater$prop.sampled>=0.5| agwater$landcover=="Row Crop" & agwater$prop.sampled>=0.5) agwaterrice$nfloodbc<-ifelse(agwaterrice$year<2012, round(agwaterrice$n.flooded-agwaterrice$n.flooded*.0297), round(agwaterrice$n.flooded-agwaterrice$n.flooded*.0446)) agwaterrice$sprecip4 = agwaterrice$precip.4wk/100 agwaterrice$sprecip2 = agwaterrice$precip.2wk/100 agwatercorn$nfloodbc<-ifelse(agwatercorn$year<2012, round(agwatercorn$n.flooded-agwatercorn$n.flooded*-0.06522), round(agwatercorn$n.flooded-agwatercorn$n.flooded*0.05319)) agwatercorn$sprecip4 = agwatercorn$precip.4wk/100 agwatercorn$sprecip2 = agwatercorn$precip.2wk/100 agwaterother$nfloodbc<-ifelse(agwaterother$year<2012, round(agwaterother$n.flooded-agwaterother$n.flooded*0.003), round(agwaterother$n.flooded-agwaterother$n.flooded*0.001)) agwaterother$sprecip4 = agwaterother$precip.4wk/100 agwaterother$sprecip2 = agwaterother$precip.2wk/100 #################################################### ## split each crop and region combination into 2000-2011 non-drought, 2000-2011 drought and.. ## 2013-2015 which is sever drought agwaterrice1315<-subset(agwaterrice, agwaterrice$year>2011) agwaterriceD<-subset(agwaterrice, agwaterrice$year<2013 & agwaterrice$drought==1) agwaterriceND<-subset(agwaterrice, agwaterrice$year<2013 & agwaterrice$drought==0) agwatercorn1315<-subset(agwatercorn, agwatercorn$year>2011) agwatercornD<-subset(agwatercorn, agwatercorn$year<2013 & agwatercorn$drought==1) agwatercornND<-subset(agwatercorn, agwatercorn$year<2013 & agwatercorn$drought==0) agwaterother1315<-subset(agwaterother, agwaterother$year>2011) agwaterotherD<-subset(agwaterother, agwaterother$year<2013 & agwaterother$drought==1) agwaterotherND<-subset(agwaterother, agwaterother$year<2013 & agwaterother$drought==0) ######################################################################## ## ## fit models for each combination of region and cover type ## model set agriculture ## RICE # 2000-2011 drought agwaterricef1<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=agwaterriceD, weights=agwaterriceD$prop.sampled) ####### # 2000-2011 non-drought agwaterricef1ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=agwaterriceND, weights=agwaterriceND$prop.sampled) ## 2013-2015 drought agwaterricef1RD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=agwaterrice1315, weights=agwaterrice1315$prop.sampled) ## CORN # 2000-2011 drought agwatercornf1<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=agwatercornD, weights=agwatercornD$prop.sampled) ####### # 2000-2011 non-drought agwatercornf1ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=agwatercornND, weights=agwatercornND$prop.sampled) ## 2013-2015 drought agwatercornf1RD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=agwatercorn1315, weights=agwatercorn1315$prop.sampled) # OTHER CROPS # 2000-2011 drought agwaterotherf1<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=agwaterotherD, weights=agwaterotherD$prop.sampled) ####### # 2000-2011 non-drought agwaterotherf1ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=agwaterotherND, weights=agwaterotherND$prop.sampled) ## 2013-2015 drought agwaterotherf1RD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), random=~(1|water.year) , family=binomial, data=agwaterother1315, weights=agwaterother1315$prop.sampled) ########################################################## ## now we need to make plots of each crop of interest ## code to standardize decimal places scaleFUN <- function(x) sprintf("%.2f", x) ## data for plots from the top model Dagwaterricef1ND = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwaterricef1ND$zday = (Dagwaterricef1ND$yday-150)/100 Dagwaterricef1ND$prop.flooded = predict(agwaterricef1ND$gam, Dagwaterricef1ND, type='response') predagwaterriceND = as.data.frame(predict(agwaterricef1ND$gam, Dagwaterricef1ND, type='link', se.fit=T)) Dagwaterricef1ND$lcl = plogis(predagwaterriceND$fit-2*predagwaterriceND$se.fit) Dagwaterricef1ND$ucl = plogis(predagwaterriceND$fit+2*predagwaterriceND$se.fit) Dagwaterricef1RD = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwaterricef1RD$zday = (Dagwaterricef1RD$yday-150)/100 Dagwaterricef1RD$prop.flooded = predict(agwaterricef1RD$gam, Dagwaterricef1RD, type='response') predagwaterriceRD = as.data.frame(predict(agwaterricef1RD$gam, Dagwaterricef1RD, type='link', se.fit=T)) Dagwaterricef1RD$lcl = plogis(predagwaterriceRD$fit-2*predagwaterriceRD$se.fit) Dagwaterricef1RD$ucl = plogis(predagwaterriceRD$fit+2*predagwaterriceRD$se.fit) Dagwaterricef1D = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwaterricef1D$zday = (Dagwaterricef1D$yday-150)/100 Dagwaterricef1D$prop.flooded = predict(agwaterricef1$gam, Dagwaterricef1D, type='response') predagwaterriceD = as.data.frame(predict(agwaterricef1$gam, Dagwaterricef1D, type='link', se.fit=T)) Dagwaterricef1D$lcl = plogis(predagwaterriceD$fit-2*predagwaterriceD$se.fit) Dagwaterricef1D$ucl = plogis(predagwaterriceD$fit+2*predagwaterriceD$se.fit) newdatagwaterriceDND<-data.frame(rbind(Dagwaterricef1ND, Dagwaterricef1D, Dagwaterricef1RD)) newdatagwaterriceDND$g<-c(rep("non-drought",319), rep("drought",319), rep("critical",319)) ## RICE rice<-ggplot(data = newdatagwaterriceDND,aes(x = yday, y = prop.flooded, ymin=lcl, ymax=ucl,fill = g)) + geom_ribbon(alpha = 0.4) + scale_fill_manual(name = "", values = c('black','green',"blue"))+ #labels = c('Extreme Drought: 2013-2015','Drought: 2000-2011','Non-Drought: 2000-2011'))+ scale_x_continuous(limits=c(1,319), expand=c(0,0), breaks=c(1,32,62,93,124,154,185,216,244,275,305), labels=c('Jul','Aug','Sep','Oct','Nov','Dec','Jan','Feb','Mar','Apr','May'))+ scale_y_continuous(labels=scaleFUN)+ xlab(NULL) + ylab('Proportion Open Water')+ theme(legend.position="none")+ ggtitle("Rice") ## CORN Dagwatercornf1ND = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwatercornf1ND$zday = (Dagwatercornf1ND$yday-150)/100 Dagwatercornf1ND$prop.flooded = predict(agwatercornf1ND$gam, Dagwatercornf1ND, type='response') predagwatercornND = as.data.frame(predict(agwatercornf1ND$gam, Dagwatercornf1ND, type='link', se.fit=T)) Dagwatercornf1ND$lcl = plogis(predagwatercornND$fit-2*predagwatercornND$se.fit) Dagwatercornf1ND$ucl = plogis(predagwatercornND$fit+2*predagwatercornND$se.fit) Dagwatercornf1RD = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwatercornf1RD$zday = (Dagwatercornf1RD$yday-150)/100 Dagwatercornf1RD$prop.flooded = predict(agwatercornf1RD$gam, Dagwatercornf1RD, type='response') predagwatercornRD = as.data.frame(predict(agwatercornf1RD$gam, Dagwatercornf1RD, type='link', se.fit=T)) Dagwatercornf1RD$lcl = plogis(predagwatercornRD$fit-2*predagwatercornRD$se.fit) Dagwatercornf1RD$ucl = plogis(predagwatercornRD$fit+2*predagwatercornRD$se.fit) Dagwatercornf1D = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwatercornf1D$zday = (Dagwatercornf1D$yday-150)/100 Dagwatercornf1D$prop.flooded = predict(agwatercornf1$gam, Dagwatercornf1D, type='response') predagwatercornD = as.data.frame(predict(agwatercornf1$gam, Dagwatercornf1D, type='link', se.fit=T)) Dagwatercornf1D$lcl = plogis(predagwatercornD$fit-2*predagwatercornD$se.fit) Dagwatercornf1D$ucl = plogis(predagwatercornD$fit+2*predagwatercornD$se.fit) newdatagwatercornDND<-data.frame(rbind(Dagwatercornf1ND, Dagwatercornf1D, Dagwatercornf1RD)) newdatagwatercornDND$g<-c(rep("non-drought",319), rep("drought",319), rep("critical",319)) corn<-ggplot(data = newdatagwatercornDND,aes(x = yday, y = prop.flooded, ymin=lcl, ymax=ucl,fill = g)) + geom_ribbon(alpha = 0.4) + scale_fill_manual(name = "", values = c('black','green',"blue"))+ #labels = c('Extreme Drought: 2013-2015','Drought: 2000-2011','Non-Drought: 2000-2011'))+ scale_x_continuous(limits=c(1,319), expand=c(0,0), breaks=c(1,32,62,93,124,154,185,216,244,275,305), labels=c('Jul','Aug','Sep','Oct','Nov','Dec','Jan','Feb','Mar','Apr','May'))+ scale_y_continuous(labels=scaleFUN)+ xlab(NULL) + ylab('Proportion Open Water')+ theme(legend.position="none")+ ggtitle("Corn") ## OTHER CROPS Dagwaterotherf1ND = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwaterotherf1ND$zday = (Dagwaterotherf1ND$yday-150)/100 Dagwaterotherf1ND$prop.flooded = predict(agwaterotherf1ND$gam, Dagwaterotherf1ND, type='response') predagwaterotherND = as.data.frame(predict(agwaterotherf1ND$gam, Dagwaterotherf1ND, type='link', se.fit=T)) Dagwaterotherf1ND$lcl = plogis(predagwaterotherND$fit-2*predagwaterotherND$se.fit) Dagwaterotherf1ND$ucl = plogis(predagwaterotherND$fit+2*predagwaterotherND$se.fit) Dagwaterotherf1RD = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwaterotherf1RD$zday = (Dagwaterotherf1RD$yday-150)/100 Dagwaterotherf1RD$prop.flooded = predict(agwaterotherf1RD$gam, Dagwaterotherf1RD, type='response') predagwaterotherRD = as.data.frame(predict(agwaterotherf1RD$gam, Dagwaterotherf1RD, type='link', se.fit=T)) Dagwaterotherf1RD$lcl = plogis(predagwaterotherRD$fit-2*predagwaterotherRD$se.fit) Dagwaterotherf1RD$ucl = plogis(predagwaterotherRD$fit+2*predagwaterotherRD$se.fit) Dagwaterotherf1D = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwaterotherf1D$zday = (Dagwaterotherf1D$yday-150)/100 Dagwaterotherf1D$prop.flooded = predict(agwaterotherf1$gam, Dagwaterotherf1D, type='response') predagwaterotherD = as.data.frame(predict(agwaterotherf1$gam, Dagwaterotherf1D, type='link', se.fit=T)) Dagwaterotherf1D$lcl = plogis(predagwaterotherD$fit-2*predagwaterotherD$se.fit) Dagwaterotherf1D$ucl = plogis(predagwaterotherD$fit+2*predagwaterotherD$se.fit) newdatagwaterotherDND<-data.frame(rbind(Dagwaterotherf1ND, Dagwaterotherf1D, Dagwaterotherf1RD)) newdatagwaterotherDND$g<-c(rep("non-drought",319), rep("drought",319), rep("critical",319)) other<-ggplot(data = newdatagwaterotherDND,aes(x = yday, y = prop.flooded, ymin=lcl, ymax=ucl,fill = g)) + geom_ribbon(alpha = 0.4) + scale_fill_manual(name = "", values = c('black','green',"blue"))+ #labels = c('Extreme Drought: 2013-2015','Drought: 2000-2011','Non-Drought: 2000-2011'))+ scale_x_continuous(limits=c(1,319), expand=c(0,0), breaks=c(1,32,62,93,124,154,185,216,244,275,305), labels=c('Jul','Aug','Sep','Oct','Nov','Dec','Jan','Feb','Mar','Apr','May'))+ scale_y_continuous(labels=scaleFUN)+ xlab(NULL) + ylab('Proportion Open Water')+ theme(legend.position="none")+ ggtitle("Other crops") ## combine into a single plot ag1<- plot_grid(rice, corn, other, nrow=3, ncol=1, align="hv") save_plot("ag_all3.png", ag1, nrow=3, ncol=1) ag2<- plot_grid(rice, corn, other, nrow=1, ncol=3, align="hv") save_plot("ag_all3b.png", ag2, nrow=1, ncol=3)