#### 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:2009,2013: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) Dagwaterricef0RD = data.frame(yday=seq(1,319,1), prop.flooded=0, n.flooded=0, n.sampled=0, prop.sampled=0) Dagwaterricef0RD$zday = (Dagwaterricef0ND$yday-150)/100 Dagwaterricef0RD$prop.flooded = predict(sacricef0RD$gam, Dagwaterricef0ND, type='response') predagwaterriceRD = as.data.frame(predict(sacricef0RD$gam, Dagwaterricef0ND, type='link', se.fit=T)) Dagwaterricef0RD$lcl = plogis(predagwaterriceND$fit-2*predagwaterriceND$se.fit) Dagwaterricef0RD$ucl = plogis(predagwaterriceND$fit+2*predagwaterriceND$se.fit) ## get estimates of toal ha 2014 ## basing rice estimate on dybaklaet al 2017 - 178,347ha Dagwaterricef0RD$ha14<-Dagwaterricef0RD$prop.flooded*178347 Dagwaterricef0RD$halcl14<-Dagwaterricef0RD$lcl*178347 Dagwaterricef0RD$haucl14<-Dagwaterricef0RD$ucl*178347 ## 2013 Dagwaterricef0RD$ha13<-Dagwaterricef0RD$prop.flooded*227242 Dagwaterricef0RD$halcl13<-Dagwaterricef0RD$lcl*227242 Dagwaterricef0RD$haucl13<-Dagwaterricef0RD$ucl*227242 ### inctdat<-read.csv("inctdat_rev052018.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$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_may2018.csv") write.csv(ptotal2summary,"ptotal2summary_may2018.csv") ### ### not corrected minimum estimates ### 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_may2018.csv") write.csv(ptotal2summary,"ptotal2summaryM_may2018.csv") ############################################# inctdat<-read.csv("inctdat_rev052018.csv") ## get the prportion contributed on the days when active inctdatA<-subset(inctdat, inctdat$ha>0) ## average per day proportion 2014 BR fall ## get bird returns in fall inctdatfallBR<-subset(inctdatA, inctdatA$yday<=122) summary(inctdatfallBR) nrow(inctdatfallBR) ## get total in fall Dagwaterricef0RDfallBR<-subset(Dagwaterricef0RD, Dagwaterricef0RD$yday>61 & Dagwaterricef0RD$yday<=122) summary(Dagwaterricef0RDfallBR) m1<-merge(inctdatfallBR,Dagwaterricef0RDfallBR, by = "yday") nrow(m1) m1$ptotal<- m1$ha/m1$ha14 m1$ptotal2<-ifelse(m1$ptotal>1,1,m1$ptotal) summary(m1$ptotal2) ## get whep ## winter is 1 nov - 31 january ## get total in fall inctdatfwinWHEP<-subset(inctdatA, inctdatA$yday>122 & inctdatA$yday<=214 ) summary(inctdatfwinWHEP) nrow(inctdatfwinWHEP) Dagwaterricef0RDwinWHEP<-subset(Dagwaterricef0RD, Dagwaterricef0RD$yday>122 & Dagwaterricef0RD$yday<=214) summary(Dagwaterricef0RDwinWHEP) m2<-merge(inctdatfwinWHEP,Dagwaterricef0RDwinWHEP, by = "yday") nrow(m2) m2$ptotal<-ifelse(m2$syear=="y1314", m2$ha/m2$ha13,m2$ha/m2$ha14) m2$ptotal2<-ifelse(m2$ptotal>1,1,m2$ptotal) summary(m2$ptotal2) ## get whep & BR spring ## spring = feb1 onwards inctdatfspWHEP<-subset(inctdatA, inctdatA$yday>214 & inctdatA$program =="WHEP") summary(inctdatfspWHEP) nrow(inctdatfspWHEP) Dagwaterricef0RDspWHEP<-subset(Dagwaterricef0RD, Dagwaterricef0RD$yday>214 ) summary(Dagwaterricef0RDspWHEP) m3<-merge(inctdatfspWHEP,Dagwaterricef0RDspWHEP, by = "yday") nrow(m3) m3$ptotal<-ifelse(m3$syear=="y1314", m3$ha/m3$ha13,m3$ha/m3$ha14) summary(m3$ptotal) m3$ptotal2<-ifelse(m3$ptotal>1,1,m3$ptotal) summary(m3$ptotal2) ## BR spring inctdatfspBR<-subset(inctdatA, inctdatA$yday>214 & inctdatA$program =="BR") summary(inctdatfspBR) nrow(inctdatfspBR) Dagwaterricef0RDspBR<-subset(Dagwaterricef0RD, Dagwaterricef0RD$yday>214 ) summary(Dagwaterricef0RDspBR) m4<-merge(inctdatfspBR,Dagwaterricef0RDspBR, by = "yday") nrow(m4) m4$ptotal<-ifelse(m4$syear=="y1314", m4$ha/m4$ha13,m4$ha/m4$ha14) summary(m4$ptotal) m4$ptotal2<-ifelse(m4$ptotal>1,1,m4$ptotal) summary(m4$ptotal2) #### create stacked bar plot incplot<-read.csv("incplot.csv",header=T) names(incplot) incplot$seas <- factor(incplot$seas, levels = c("Fall", "Winter", "Spring")) incplot$ip <- factor(incplot$ip, levels = c("NI", "WHEP", "BR")) png("incplot_apr2018.png", height = 3.5, width = 3.5, units ="in", res=300) ggplot() + geom_bar(aes(y = p, x = seas, fill = ip), data = incplot, stat="identity")+ scale_fill_manual(name = "", values = c('gray','black','blue'), labels = c('No Incentive','WHEP','Bird Returns'))+ theme(legend.position = "top", legend.text=element_text(size=9), axis.title.y =element_text(size=12),axis.title.x =element_text(size=12))+ xlab("Season")+ ylab("Proportion Total Open Water in Rice") dev.off()