#### impact of drought analysis

## Reiter et al.

## date: January 2017; mod June 2017; mod December 2017

## Asessment of the contribution of incentive programs specifically in the sac valley

# sac valley = colusa, butte, sutter and american CVJV basins

# required data:  "err2.csv";
##                 "ag_merge_scene_bybasin_precip2.csv";
##              "cvjv_basins_lookup.csv"; "inctdat_rev012018.csv"

#
## load in needed packages
library(gbm)
library(dismo)
library(ggplot2)
library(gamm4)

## get the sac valley only data isolated to run moodels. Try 20132015 and the then see if possible to split
## the flood curve by year
setwd("C:/water")

agbasin<-read.csv("ag_merge_scene_bybasin_precip2.csv", header=T)
summary(agbasin)
nrow(agbasin)

## basin name file

basin<-read.csv("cvjv_basins_lookup.csv", header = T)
names(basin)<-c("basin","basinname","region","area","perimeter")
agbasinn<-merge(agbasin, basin)
nrow(agbasinn)

## read in bias correction data for reference
bc<-read.csv("C:/water/err2.csv", header=T)

##
## drought indicator

agbasinn$drought<-ifelse(agbasinn$water.year%in%c(2001:2002,2007:2015),1,0)

agbasinn$zday = (agbasinn$yday-150)/100
agbasinn$id = c(1:nrow(agbasinn))

#split data by region

sac<-subset(agbasinn,agbasinn$basinname%in%c("Butte", "Colusa","Sutter","American"))
sj<-subset(agbasinn,agbasinn$basinname%in%c("San Joaquin"))

## split regional data by crops
##

## sac valley
sacrice<-subset(sac,sac$landcover.1=="Rice" & sac$prop.sampled>=0.5)


sacrice$nfloodbc<-ifelse(sacrice$year<2012,
                         round(sacrice$n.flooded-sacrice$n.flooded*.0297),
                         round(sacrice$n.flooded-sacrice$n.flooded*.0446))
sacrice$sprecip4 = sacrice$precip.4wk/100
sacrice$sprecip2 = sacrice$precip.2wk/100

sacrice1315<-subset(sacrice, sacrice$year>2011)
sacriceD<-subset(sacrice, sacrice$year<2013 & sacrice$drought==1)
sacriceND<-subset(sacrice, sacrice$year<2013 & sacrice$drought==0)


## sacrice

## 2013-2015 drought
sacricef0RD<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~s(zday, k=10), 
                   random=~(1|id), family=binomial, data=sacrice1315, weights=sacrice1315$prop.sampled)

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

#plots

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


## get estimates of toal ha 2014
## rice estimates from Dybala et al.(2017) - 169606 ha planted

Dagwaterricef0RD$ha14<-Dagwaterricef0RD$prop.flooded*169606
Dagwaterricef0RD$halcl14<-Dagwaterricef0RD$lcl*169606
Dagwaterricef0RD$haucl14<-Dagwaterricef0RD$ucl*169606

## 2013 - 216105 ha planted

Dagwaterricef0RD$ha13<-Dagwaterricef0RD$prop.flooded*216105
Dagwaterricef0RD$halcl13<-Dagwaterricef0RD$lcl*216105
Dagwaterricef0RD$haucl13<-Dagwaterricef0RD$ucl*216105


## read in tnc and whep data to add to plot

################ corrected with the decay

inctdat<-read.csv("inctdat_rev012018.csv")

## get only 1 year of data

inctdat1314<-subset(inctdat,inctdat$syear =="y1314")
inctdat1415<-subset(inctdat,inctdat$syear =="y1415")

png("sacricewater1314ha_inct.png", width=2000, height =1000, res=300)
ggplot(data = Dagwaterricef0RD,aes(x = yday, y = ha13)) + 
  geom_area(fill='grey70') +
  geom_bar(data=inctdat1314, aes(x=yday, y=rhaw, fill = program, alpha = 0.4), stat = "identity")+
  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('ha open water')+ ylim(0,100000)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(),
         axis.line = element_line(colour = "black"),legend.position="none")
dev.off()

##
png("sacricewater1415ha_inct.png", width=2000, height =1000, res=300)
ggplot(data = Dagwaterricef0RD,aes(x = yday, y = ha14)) + 
  geom_area(fill='grey70') +
  geom_bar(data=inctdat1415, aes(x=yday, y=rhaw, fill = program, alpha = 0.4), stat = "identity")+
  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('ha open water')+ ylim(0,100000)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
              panel.background = element_blank(),
         axis.line = element_line(colour = "black"),legend.position="none")
dev.off()



### caculate he proportion total by program and day
# first whep

inctdat1314whep<-subset(inctdat1314, inctdat1314$program=="WHEP")
ptotal1314whep<-inctdat1314whep$rhaw/Dagwaterricef0RD$ha13
total1314whep<-sum(inctdat1314whep$rhaw)
inctdat1415whep<-subset(inctdat1415, inctdat1415$program=="WHEP")
ptotal1415whep<-inctdat1415whep$rhaw/Dagwaterricef0RD$ha14
total1415whep<-sum(inctdat1415whep$rhaw)
## then BR

inctdat1314br<-subset(inctdat1314, inctdat1314$program=="BR")
ptotal1314br<-inctdat1314br$rhaw/Dagwaterricef0RD$ha13
total1314br<-sum(inctdat1314br$rhaw)
inctdat1415br<-subset(inctdat1415, inctdat1415$program=="BR")
ptotal1415br<-inctdat1415br$rhaw/Dagwaterricef0RD$ha14
total1415br<-sum(inctdat1415br$rhaw)

## combine into table

ptotal<-data.frame(cbind(ptotal1314whep,ptotal1415whep,ptotal1314br,ptotal1415br))
str(ptotal)
names(ptotal)

## get summary stats

s1<-summary(ptotal$ptotal1314whep)
s2<-summary(ptotal$ptotal1415whep)
s3<-summary(ptotal$ptotal1314br)
s4<-summary(ptotal$ptotal1415br)

ptotalsummary<-data.frame(rbind(s1,s2,s3,s4))
ptotalsummary$yr<-c("y1314","y1415","y1314","y1415")
ptotalsummary$pr<-c("WHEP","WHEP","BR","BR")
ptotalsummary$yrpr<-c("y1314w","y1415w","y1314b","y1415b")
ptotalsummary$sd<-c(sd(ptotal$ptotal1314whep),sd(ptotal$ptotal1415whep),
                    sd(ptotal$ptotal1314br),sd(ptotal$ptotal1415br))
ptotalsummary$sum<-(c(total1314whep, total1415whep,total1314br,total1415br))

ptotalsummary$n<-c(nrow(ptotal),nrow(ptotal),
                    nrow(ptotal),nrow(ptotal))## issue with some days actually exceeding what we expect based on the 2 year average ..hmmm


## remove those months
### caculate he proportion total by program and day
# first whep

inctdat1314whep<-subset(inctdat1314, inctdat1314$program=="WHEP")
ptotal1314whep<-data.frame(inctdat1314whep$rhaw/Dagwaterricef0RD$ha13)
names(ptotal1314whep)<-c("ptotal1314w")
ptotal1314whep$i<-ifelse(ptotal1314whep$ptotal1314w<=1&ptotal1314whep$ptotal1314w>0,"lt1","gt1")
ptotal1314whep2<-subset(ptotal1314whep,ptotal1314whep$i=="lt1")

inctdat1415whep<-subset(inctdat1415, inctdat1415$program=="WHEP")
ptotal1415whep<-data.frame(inctdat1415whep$rhaw/Dagwaterricef0RD$ha14)
names(ptotal1415whep)<-c("ptotal1415w")
ptotal1415whep$i<-ifelse(ptotal1415whep$ptotal1415w<=1&ptotal1415whep$ptotal1415w>0,"lt1","gt1")
ptotal1415whep2<-subset(ptotal1415whep,ptotal1415whep$i=="lt1")
names(ptotal1415whep2)
## then BR

inctdat1314br<-subset(inctdat1314, inctdat1314$program=="BR")
ptotal1314br<-data.frame(inctdat1314br$rhaw/Dagwaterricef0RD$ha13)
names(ptotal1314br)<-c("ptotal1314b")
ptotal1314br$i<-ifelse(ptotal1314br$ptotal1314b<=1&ptotal1314br$ptotal1314b>0,"lt1","gt1")
ptotal1314br2<-subset(ptotal1314br,ptotal1314br$i=="lt1")

inctdat1415br<-subset(inctdat1415, inctdat1415$program=="BR")
ptotal1415br<-data.frame(inctdat1415br$rhaw/Dagwaterricef0RD$ha14)
names(ptotal1415br)<-c("ptotal1415b")
ptotal1415br$i<-ifelse(ptotal1415br$ptotal1415b<=1.0&ptotal1415br$ptotal1415b>0,"lt1","gt1")
ptotal1415br2<-subset(ptotal1415br,ptotal1415br$i=="lt1")
summary(ptotal1415br2)


## get summary stats

s12<-summary(ptotal1314whep2$ptotal1314w)
s22<-summary(ptotal1415whep2$ptotal1415w)
s32<-summary(ptotal1314br2$ptotal1314b)
s42<-summary(ptotal1415br2$ptotal1415b)

ptotal2summary<-data.frame(rbind(s12,s22,s32,s42))
ptotal2summary$yr<-c("y1314","y1415","y1314","y1415")
ptotal2summary$pr<-c("WHEP","WHEP","BR","BR")
ptotal2summary$yrpr<-c("y1314w","y1415w","y1314b","y1415b")
ptotal2summary$sd<-c(sd(ptotal1314whep2$ptotal1314w),sd(ptotal1415whep2$ptotal1415w),
                    sd(ptotal1314br2$ptotal1314b),sd(ptotal1415br2$ptotal1415b))

ptotal2summary$n<-c(nrow(ptotal1314whep2$ptotal1314w),nrow(ptotal1415whep2$ptotal1415w),
                     nrow(ptotal1314br2$ptotal1314b),nrow(ptotal1415br2$ptotal1415b))


# writet data to tabel

write.csv(ptotalsummary,"ptotalsummary.csv")
write.csv(ptotal2summary,"ptotal2summary.csv")

###

## read in tnc data to add to plot

### not corrected minimum estimates

inctdat<-read.csv("inctdat_rev012018.csv")


## make TNC/whep file that matches the planing cycle for each year
## need to add in WHEP acres

## get only 1 year of data

inctdat1314<-subset(inctdat,inctdat$syear =="y1314")
inctdat1415<-subset(inctdat,inctdat$syear =="y1415")


### caculate he proportion total by program and day
# first whep

inctdat1314whep<-subset(inctdat1314, inctdat1314$program=="WHEP")
ptotal1314whep<-inctdat1314whep$ha/Dagwaterricef0RD$ha13
total1314whep<-sum(inctdat1314whep$ha)
inctdat1415whep<-subset(inctdat1415, inctdat1415$program=="WHEP")
ptotal1415whep<-inctdat1415whep$ha/Dagwaterricef0RD$ha14
total1415whep<-sum(inctdat1415whep$ha)

## then BR

inctdat1314br<-subset(inctdat1314, inctdat1314$program=="BR")
ptotal1314br<-inctdat1314br$ha/Dagwaterricef0RD$ha13
total1314br<-sum(inctdat1314br$ha)
inctdat1415br<-subset(inctdat1415, inctdat1415$program=="BR")
ptotal1415br<-inctdat1415br$ha/Dagwaterricef0RD$ha14
total1415br<-sum(inctdat1415br$ha)

## combine into table

ptotal<-data.frame(cbind(ptotal1314whep,ptotal1415whep,ptotal1314br,ptotal1415br))
str(ptotal)
names(ptotal)

## get summary stats

s1<-summary(ptotal$ptotal1314whep)
s2<-summary(ptotal$ptotal1415whep)
s3<-summary(ptotal$ptotal1314br)
s4<-summary(ptotal$ptotal1415br)

ptotalsummary<-data.frame(rbind(s1,s2,s3,s4))
ptotalsummary$yr<-c("y1314","y1415","y1314","y1415")
ptotalsummary$pr<-c("WHEP","WHEP","BR","BR")
ptotalsummary$yrpr<-c("y1314w","y1415w","y1314b","y1415b")
ptotalsummary$sd<-c(sd(ptotal$ptotal1314whep),sd(ptotal$ptotal1415whep),
                    sd(ptotal$ptotal1314br),sd(ptotal$ptotal1415br))

ptotalsummary$sum<-(c(total1314whep, total1415whep,total1314br,total1415br))
ptotalsummary$n<-c(nrow(ptotal),nrow(ptotal),
                   nrow(ptotal),nrow(ptotal))## issue with some days actually exceeding what we expect based on the 2 year average ..hmmm


## remove those months
### caculate he proportion total by program and day
# first whep

inctdat1314whep<-subset(inctdat1314, inctdat1314$program=="WHEP")
ptotal1314whep<-data.frame(inctdat1314whep$ha/Dagwaterricef0RD$ha13)
names(ptotal1314whep)<-c("ptotal1314w")
ptotal1314whep$i<-ifelse(ptotal1314whep$ptotal1314w<=1&ptotal1314whep$ptotal1314w>0,"lt1","gt1")
ptotal1314whep2<-subset(ptotal1314whep,ptotal1314whep$i=="lt1")
inctdat1415whep<-subset(inctdat1415, inctdat1415$program=="WHEP")
ptotal1415whep<-data.frame(inctdat1415whep$ha/Dagwaterricef0RD$ha14)
names(ptotal1415whep)<-c("ptotal1415w")
ptotal1415whep$i<-ifelse(ptotal1415whep$ptotal1415w<=1&ptotal1415whep$ptotal1415w>0,"lt1","gt1")
ptotal1415whep2<-subset(ptotal1415whep,ptotal1415whep$i=="lt1")
names(ptotal1415whep2)

## then BR

inctdat1314br<-subset(inctdat1314, inctdat1314$program=="BR")
ptotal1314br<-data.frame(inctdat1314br$ha/Dagwaterricef0RD$ha13)
names(ptotal1314br)<-c("ptotal1314b")
ptotal1314br$i<-ifelse(ptotal1314br$ptotal1314b<=1&ptotal1314br$ptotal1314b>0,"lt1","gt1")
ptotal1314br2<-subset(ptotal1314br,ptotal1314br$i=="lt1")

inctdat1415br<-subset(inctdat1415, inctdat1415$program=="BR")
ptotal1415br<-data.frame(inctdat1415br$ha/Dagwaterricef0RD$ha14)
names(ptotal1415br)<-c("ptotal1415b")
ptotal1415br$i<-ifelse(ptotal1415br$ptotal1415b<=1.0&ptotal1415br$ptotal1415b>0,"lt1","gt1")
ptotal1415br2<-subset(ptotal1415br,ptotal1415br$i=="lt1")
summary(ptotal1415br2)

## get summary stats

s12<-summary(ptotal1314whep2$ptotal1314w)
s22<-summary(ptotal1415whep2$ptotal1415w)
s32<-summary(ptotal1314br2$ptotal1314b)
s42<-summary(ptotal1415br2$ptotal1415b)

ptotal2summary<-data.frame(rbind(s12,s22,s32,s42))
ptotal2summary$yr<-c("y1314","y1415","y1314","y1415")
ptotal2summary$pr<-c("WHEP","WHEP","BR","BR")
ptotal2summary$yrpr<-c("y1314w","y1415w","y1314b","y1415b")
ptotal2summary$sd<-c(sd(ptotal1314whep2$ptotal1314w),sd(ptotal1415whep2$ptotal1415w),
                     sd(ptotal1314br2$ptotal1314b),sd(ptotal1415br2$ptotal1415b))

ptotal2summary$n<-c(nrow(ptotal1314whep2$ptotal1314w),nrow(ptotal1415whep2$ptotal1415w),
                    nrow(ptotal1314br2$ptotal1314b),nrow(ptotal1415br2$ptotal1415b))

haobs<-data.frame(c(sum(inctdat1314whep$ha),sum(inctdat1415whep$ha),
                    sum(inctdat1314whep$ha),sum(inctdat1415whep$ha)))

# write data to table

write.csv(ptotalsummary,"ptotalsummaryM.csv")
write.csv(ptotal2summary,"ptotal2summaryM.csv")