#### impact of drought analysis

### wetlands

### modified from Nov 2016 file to consoder modeling bias corrected values
## bias correction derived from comparison of clasifiction with inow water and nonwater points
## correction separate for freshwater emergent wetland, corn and rice

## matt reiter january 31, 2017

## updated may 2017

## updated december 2017

## new data set now being used.

## addded in data to distinguish between public and private wetlands

## This code is for modeling wetland and ag individually by year set
## It includes code for making plots shown in the manuscript


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)

### make public provate = 1,0

wetwater1seas$lc<-ifelse(wetwater1seas$landcover==10,1,0)


### 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))

### make public provate = 1,0

wetwater1seasD$lc<-ifelse(wetwater1seasD$landcover==10,1,0)


f1D<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), 
           random=~(1|water.year) , family=binomial, data=wetwater1seasD, weights=wetwater1seasD$prop.sampled)

f2D<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10) +sprecip2, 
           random=~(1|water.year), family=binomial, data=wetwater1seasD, weights=wetwater1seasD$prop.sampled)

f3D<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10) +sprecip4, 
           random=~(1|water.year) , family=binomial, data=wetwater1seasD, weights=wetwater1seasD$prop.sampled)

f4D<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10) +lc, 
           random=~(1|water.year), family=binomial, data=wetwater1seasD, weights=wetwater1seasD$prop.sampled)

f5D<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10) +lc+sprecip2, 
           random=~(1|water.year), family=binomial, data=wetwater1seasD, weights=wetwater1seasD$prop.sampled)

f6D<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)  +lc+sprecip4, 
           random=~(1|water.year) , family=binomial, data=wetwater1seasD, weights=wetwater1seasD$prop.sampled)
summary(f1D$gam)
summary(f2D$gam)
summary(f3D$gam)
summary(f4D$gam)
summary(f5D$gam)
summary(f6D$gam)


save(f1D, file ="f1Dnoid.rda")
save(f2D, file ="f2Dnoid.rda")
save(f3D, file ="f3Dnoid.rda")

save(f4D, file ="f4Dnoid.rda")
save(f5D, file ="f5Dnoid.rda")
save(f6D, file ="f6Dnoid.rda")

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)

### make public private = 1,0

wetwater1seasND$lc<-ifelse(wetwater1seasND$landcover==10,1,0)


f1ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), 
            random=~(1|water.year) , family=binomial, data=wetwater1seasND, weights=wetwater1seasND$prop.sampled)

f2ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+sprecip2, 
            random=~(1|water.year) , family=binomial, data=wetwater1seasND, weights=wetwater1seasND$prop.sampled)

f3ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+sprecip4, 
            random=~(1|water.year) , family=binomial, data=wetwater1seasND, weights=wetwater1seasND$prop.sampled)
f4ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc, 
            random=~(1|water.year) , family=binomial, data=wetwater1seasND, weights=wetwater1seasND$prop.sampled)

f5ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc+sprecip2, 
            random=~(1|water.year) , family=binomial, data=wetwater1seasND, weights=wetwater1seasND$prop.sampled)

f6ND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc+sprecip4, 
            random=~(1|water.year) , family=binomial, data=wetwater1seasND, weights=wetwater1seasND$prop.sampled)

summary(f1ND$gam)
summary(f2ND$gam)
summary(f3ND$gam)
summary(f4ND$gam)
summary(f5ND$gam)
summary(f6ND$gam)


save(f1ND, file ="f1NDnoid.rda")
save(f2ND, file ="f2NDnoid.rda")
save(f3ND, file ="f3NDnoid.rda")
save(f4ND, file ="f4NDnoid.rda")
save(f5ND, file ="f5NDnoid.rda")
save(f6ND, file ="f6NDnoid.rda")


## data for plots from the top 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)

### make public provate = 1,0

wetwater2seas$lc<-ifelse(wetwater2seas$landcover==10,1,0)


f1r<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), 
            random=~(1|water.year) , family=binomial, data=wetwater2seas, weights=wetwater2seas$prop.sampled)

f2r<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+sprecip2, 
            random=~(1|water.year) , family=binomial, data=wetwater2seas, weights=wetwater2seas$prop.sampled)

f3r<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+sprecip4, 
            random=~(1|water.year) , family=binomial, data=wetwater2seas, weights=wetwater2seas$prop.sampled)
f4r<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc, 
            random=~(1|water.year) , family=binomial, data=wetwater2seas, weights=wetwater2seas$prop.sampled)

f5r<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc+sprecip2, 
            random=~(1|water.year) , family=binomial, data=wetwater2seas, weights=wetwater2seas$prop.sampled)

f6r<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc+sprecip4, 
            random=~(1|water.year) , family=binomial, data=wetwater2seas, weights=wetwater2seas$prop.sampled)

summary(f1r$gam)
summary(f2r$gam)
summary(f3r$gam)
summary(f4r$gam)
summary(f5r$gam)
summary(f6r$gam)


save(f1r, file ="f1rnoid.rda")
save(f2r, file ="f2rnoid.rda")
save(f3r, file ="f3rnoid.rda")
save(f4r, file ="f4rnoid.rda")
save(f5r, file ="f5rnoid.rda")
save(f6r, file ="f6rnoid.rda")

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))

png("p111DNDeseas_dec17noid_noleg.png", width=2000, height =1000, res=300)
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")
dev.off()



#######################################################################

## 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)

### make public private = 1,0

wetwater1semi$lc<-ifelse(wetwater1semi$landcover==20,1,0)

### 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))

### make public private = 1,0

wetwater1semiD$lc<-ifelse(wetwater1semiD$landcover==20,1,0)


f1sD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), 
            random=~(1|water.year) , family=binomial, data=wetwater1semiD, weights=wetwater1semiD$prop.sampled)

f2sD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10) +sprecip2, 
            random=~(1|water.year), family=binomial, data=wetwater1semiD, weights=wetwater1semiD$prop.sampled)

f3sD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10) +sprecip4, 
            random=~(1|water.year) , family=binomial, data=wetwater1semiD, weights=wetwater1semiD$prop.sampled)

f4sD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10) +lc, 
            random=~(1|water.year), family=binomial, data=wetwater1semiD, weights=wetwater1semiD$prop.sampled)

f5sD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10) +lc+sprecip2, 
            random=~(1|water.year), family=binomial, data=wetwater1semiD, weights=wetwater1semiD$prop.sampled)

f6sD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)  +lc+sprecip4, 
            random=~(1|water.year) , family=binomial, data=wetwater1semiD, weights=wetwater1semiD$prop.sampled)
summary(f1sD$gam)
summary(f2sD$gam)
summary(f3sD$gam)
summary(f4sD$gam)
summary(f5sD$gam)
summary(f6sD$gam)


save(f1sD, file ="f1sDnoid.rda")
save(f2sD, file ="f2sDnoid.rda")
save(f3sD, file ="f3sDnoid.rda")

save(f4sD, file ="f4sDnoid.rda")
save(f5sD, file ="f5sDnoid.rda")
save(f6sD, file ="f6sDnoid.rda")

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)
### make public provate = 1,0

wetwater1semiND$lc<-ifelse(wetwater1semiND$landcover==20,1,0)


f1sND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), 
             random=~(1|water.year) , family=binomial, data=wetwater1semiND, weights=wetwater1semiND$prop.sampled)

f2sND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+sprecip2, 
             random=~(1|water.year) , family=binomial, data=wetwater1semiND, weights=wetwater1semiND$prop.sampled)

f3sND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+sprecip4, 
             random=~(1|water.year) , family=binomial, data=wetwater1semiND, weights=wetwater1semiND$prop.sampled)
f4sND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc, 
             random=~(1|water.year) , family=binomial, data=wetwater1semiND, weights=wetwater1semiND$prop.sampled)

f5sND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc+sprecip2, 
             random=~(1|water.year) , family=binomial, data=wetwater1semiND, weights=wetwater1semiND$prop.sampled)

f6sND<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc+sprecip4, 
             random=~(1|water.year) , family=binomial, data=wetwater1semiND, weights=wetwater1semiND$prop.sampled)

summary(f1sND$gam)
summary(f2sND$gam)
summary(f3sND$gam)
summary(f4sND$gam)
summary(f5sND$gam)
summary(f6sND$gam)


save(f1sND, file ="f1sNDnoid.rda")
save(f2sND, file ="f2sNDnoid.rda")
save(f3sND, file ="f3sNDnoid.rda")
save(f4sND, file ="f4sNDnoid.rda")
save(f5sND, file ="f5sNDnoid.rda")
save(f6sND, file ="f6sNDnoid.rda")


## 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 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)

wetwater2semi<-subset(wetwater2,wetwater2$landcover==2&wetwater2$prop.sampled>=0.5|wetwater2$landcover==20&wetwater2$prop.sampled>=0.5)

### make public provate = 1,0

wetwater2semi$lc<-ifelse(wetwater2semi$landcover==20,1,0)


f1sr<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), 
            random=~(1|water.year) , family=binomial, data=wetwater2semi, weights=wetwater2semi$prop.sampled)

f2sr<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+sprecip2, 
            random=~(1|water.year) , family=binomial, data=wetwater2semi, weights=wetwater2semi$prop.sampled)

f3sr<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+sprecip4, 
            random=~(1|water.year) , family=binomial, data=wetwater2semi, weights=wetwater2semi$prop.sampled)
f4sr<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc, 
            random=~(1|water.year) , family=binomial, data=wetwater2semi, weights=wetwater2semi$prop.sampled)

f5sr<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc+sprecip2, 
            random=~(1|water.year) , family=binomial, data=wetwater2semi, weights=wetwater2semi$prop.sampled)

f6sr<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10)+lc+sprecip4, 
            random=~(1|water.year) , family=binomial, data=wetwater2semi, weights=wetwater2semi$prop.sampled)


summary(f1sr$gam)
summary(f2sr$gam)
summary(f3sr$gam)
summary(f4sr$gam)
summary(f5sr$gam)
summary(f6sr$gam)

names(f1sr$gam)
summary(f1sr$gam)$r.sq


save(f1sr, file ="f1srnoid.rda")
save(f2sr, file ="f2srnoid.rda")
save(f3sr, file ="f3srnoid.rda")
save(f4sr, file ="f4srnoid.rda")
save(f5sr, file ="f5srnoid.rda")
save(f6sr, file ="f6srnoid.rda")

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))

png("p111DNDesemi_dec17noid_noleg.png", width=2000, height =1000, res=300)
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")
dev.off()


#######

# combine 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)



########################

## make plots of the effect of public versus private lands

## used model with just LC in it

## create data

load(file = "C:/water/f4srnoid.rda")
load(file = "C:/water/f4sDnoid.rda")
load(file = "C:/water/f4Dnoid.rda")
load(file = "C:/water/f4rnoid.rda")
load(file = "C:/water/f4NDnoid.rda")
load(file = "C:/water/f4NDnoid.rda")
summary(f4r$gam)

newdatO = data.frame(yday=rep(seq(1,319,1),2), lc = c(rep(1,319),rep(0,319)),prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0)
newdatO$zday = (newdatO$yday-150)/100
newdatO$prop.flooded = predict(f4sr$gam, newdatO, type='response')
pred = as.data.frame(predict(f4sr$gam, newdatO, type='link', se.fit=T))
newdatO$lcl = plogis(pred$fit-2*pred$se.fit)
newdatO$ucl = plogis(pred$fit+2*pred$se.fit)

newdatO$own<-ifelse(newdatO$lc==1,"pc","pr")

png("f4srnoidPP.png", width=2000,height=1000, res=300)
ggplot(data = newdatO,aes(x = yday, y = prop.flooded, ymin=lcl, ymax=ucl,fill = own)) + 
  geom_ribbon(alpha = 0.4) + 
  scale_fill_manual(name = "Legend",
                    values = c("black","grey"),
                    labels = c('Public', 'Private'))+
  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'))+
  xlab(NULL) + ylab('Proportion open water')+
  theme_bw()
dev.off()



#############

#r-squared and AIC from model

## seasonal drought

load(file = "C:/water/f1Dnoid.rda")
load(file = "C:/water/f2Dnoid.rda")
load(file = "C:/water/f3Dnoid.rda")
load(file = "C:/water/f4Dnoid.rda")
load(file = "C:/water/f5Dnoid.rda")
load(file = "C:/water/f6Dnoid.rda")

## make AIC table with adj Rsq values

aic<-c(AIC(f1D$mer),AIC(f2D$mer),AIC(f3D$mer),AIC(f4D$mer),AIC(f5D$mer),AIC(f6D$mer))
mod<-c("d1","d2","d3","d4","d5","d6")
aictabD<-data.frame(aic,mod)
aictabD$daic<-aictabD$aic-(min(aictabD$aic))
aictabD$wcalc<-exp(-0.5*aictabD$daic)
aictabD$w<-aictabD$wcalc/sum(aictabD$wcalc)
aictabD$rsq<-c(summary(f1D$gam)$r.sq,summary(f2D$gam)$r.sq,summary(f3D$gam)$r.sq,
               summary(f4D$gam)$r.sq,summary(f5D$gam)$r.sq,summary(f6D$gam)$r.sq)

write.csv(aictabD,"aictabD.csv")

## seasonal no drought

load(file = "C:/water/f1NDnoid.rda")
load(file = "C:/water/f2NDnoid.rda")
load(file = "C:/water/f3NDnoid.rda")
load(file = "C:/water/f4NDnoid.rda")
load(file = "C:/water/f5NDnoid.rda")
load(file = "C:/water/f6NDnoid.rda")


## make AIC table with adj Rsq values

aic<-c(AIC(f1ND$mer),AIC(f2ND$mer),AIC(f3ND$mer),AIC(f4ND$mer),AIC(f5ND$mer),AIC(f6ND$mer))
mod<-c("nd1","nd2","nd3","nd4","nd5","nd6")
aictabND<-data.frame(aic,mod)
aictabND$daic<-aictabND$aic-(min(aictabND$aic))
aictabND$wcalc<-exp(-0.5*aictabND$daic)
aictabND$w<-aictabND$wcalc/sum(aictabND$wcalc)
aictabND$rsq<-c(summary(f1ND$gam)$r.sq,summary(f2ND$gam)$r.sq,summary(f3ND$gam)$r.sq,
               summary(f4ND$gam)$r.sq,summary(f5ND$gam)$r.sq,summary(f6ND$gam)$r.sq)

write.csv(aictabND,"aictabND.csv")




## recent drought



load(file = "C:/water/f1rnoid.rda")
load(file = "C:/water/f2rnoid.rda")
load(file = "C:/water/f3rnoid.rda")
load(file = "C:/water/f4rnoid.rda")
load(file = "C:/water/f5rnoid.rda")
load(file = "C:/water/f6rnoid.rda")


## make AIC table with adj Rsq values

aic<-c(AIC(f1r$mer),AIC(f2r$mer),AIC(f3r$mer),AIC(f4r$mer),AIC(f5r$mer),AIC(f6r$mer))
mod<-c("r1","r2","r3","r4","r5","r6")
aictabr<-data.frame(aic,mod)
aictabr$daic<-aictabr$aic-(min(aictabr$aic))
aictabr$wcalc<-exp(-0.5*aictabr$daic)
aictabr$w<-aictabr$wcalc/sum(aictabr$wcalc)
aictabr$rsq<-c(summary(f1r$gam)$r.sq,summary(f2r$gam)$r.sq,summary(f3r$gam)$r.sq,
                summary(f4r$gam)$r.sq,summary(f5r$gam)$r.sq,summary(f6r$gam)$r.sq)


write.csv(aictabr,"aictabr.csv")