# Code for producing snow-related metrics at Seeley snowshoe hare site and related figures
# 
# Author: jared.oyler
###############################################################################

#Source in required functions
source("/nfshome/jared.oyler/ecl_helios_workspace/snowshoe_hare/swat_snow.R", echo=FALSE)
source("/nfshome/jared.oyler/ecl_helios_workspace/snowshoe_hare/snowhare_funcs.R", echo=FALSE)

#Load met data and model SWE using newly calibrated params
setwd('/projects/daymet2/hare_sites/hare_wxTopo')
metdata <- load.metdata("metdata_seeley.csv")
snw.params <- to.params(c(1,0,0.05337819,0.00000000,3.20078580)) #newly calibrated params
#snw.params <- to.params(c(1,0,0.66942390,0.04473933,2.01484801))
metdata$SWE <- swat_snow(metdata,snw.params)


##########################################
#Snow metrics for 2009 - 2012
##########################################

swe.seasons <- calc.hare.snow.seasons(metdata)

swe.len.2010 <- swe.seasons$swe.len[swe.seasons$hare.yr == 2009]
swe.len.2011 <- swe.seasons$swe.len[swe.seasons$hare.yr == 2010]
swe.len.2012 <- swe.seasons$swe.len[swe.seasons$hare.yr == 2011]

swe.start.2010 <- swe.seasons$start.yday[swe.seasons$hare.yr == 2009]
swe.start.2011 <- swe.seasons$start.yday[swe.seasons$hare.yr == 2010]
swe.start.2012 <- swe.seasons$start.yday[swe.seasons$hare.yr == 2011]

swe.end.2010 <- swe.seasons$end.yday[swe.seasons$hare.yr == 2009]
swe.end.2011 <- swe.seasons$end.yday[swe.seasons$hare.yr == 2010]
swe.end.2012 <- swe.seasons$end.yday[swe.seasons$hare.yr == 2011]

max.swe.2010<-max(metdata$SWE[metdata$year.hare==2009])
max.swe.2011<-max(metdata$SWE[metdata$year.hare==2010])
max.swe.2012<-max(metdata$SWE[metdata$year.hare==2011])

####################################################################################
#Load and apply climate projection deltas and run snow model for projections
####################################################################################

#Limit met data to the baseline 1980s
metdata<-metdata[metdata$year >= 1970 & metdata$year <= 1999,]

#Directory where projection delta files are located
setwd('/projects/daymet2/hare_sites/hare_wxTopo/cmip5_deltas')

d.rcp45.2040s <- load.all.deltas("RCP45","20312060")
d.rcp45.2080s <- load.all.deltas("RCP45","20712100")
d.rcp85.2040s <- load.all.deltas("RCP85","20312060")
d.rcp85.2080s <- load.all.deltas("RCP85","20712100")

met.base <- list(base=metdata)
met.rcp45.2040s <- apply.all.deltas(metdata,d.rcp45.2040s)
met.rcp45.2080s <- apply.all.deltas(metdata,d.rcp45.2080s)
met.rcp85.2040s <- apply.all.deltas(metdata,d.rcp85.2040s)
met.rcp85.2080s <- apply.all.deltas(metdata,d.rcp85.2080s)

met.rcp45.2040s <- snow.all.deltas(met.rcp45.2040s,snw.params)
met.rcp45.2080s <- snow.all.deltas(met.rcp45.2080s,snw.params)
met.rcp85.2040s <- snow.all.deltas(met.rcp85.2040s,snw.params)
met.rcp85.2080s <- snow.all.deltas(met.rcp85.2080s,snw.params)

####################################################################################
#SNOW SEASON LENGTH AND PLOT
####################################################################################
swe.len.base <- swe.season.all.deltas(met.base,"1980s","swe.len")
swe.len.rcp45.2040s <- swe.season.all.deltas(met.rcp45.2040s,"RCP4.5\n2040s","swe.len")
swe.len.rcp45.2080s <- swe.season.all.deltas(met.rcp45.2080s,"RCP4.5\n2080s","swe.len")
swe.len.rcp85.2040s <- swe.season.all.deltas(met.rcp85.2040s,"RCP8.5\n2040s","swe.len")
swe.len.rcp85.2080s <- swe.season.all.deltas(met.rcp85.2080s,"RCP8.5\n2080s","swe.len")
window.swe.len <- append(swe.len.base$group.name,c(swe.len.rcp45.2040s$group.name,swe.len.rcp85.2040s$group.name,swe.len.rcp45.2080s$group.name,swe.len.rcp85.2080s$group.name))
window.swe.len <- factor(window.swe.len,levels=c("1980s","RCP4.5\n2040s","RCP8.5\n2040s","RCP4.5\n2080s","RCP8.5\n2080s"))
swe.len <- append(swe.len.base$swe.len,c(swe.len.rcp45.2040s$swe.len,swe.len.rcp85.2040s$swe.len,swe.len.rcp45.2080s$swe.len,swe.len.rcp85.2080s$swe.len))

png("/projects/daymet2/hare_sites/hare_wxTopo/figs/length_boxplt.png")
boxplot(swe.len~window.swe.len,col=c("royalblue1","orange","red","orange","red"),xlab="Time Period",ylab="Days",main="Modeled Main Snow Season Duration",axes=FALSE,frame.plot=TRUE)
axis(2)
axis(1,tick=0,at=c(1,2,3,4,5),labels=levels(window.swe.len))
means.len <- tapply(swe.len,window.swe.len,mean)
points(means.len,pch=18,cex=1.5)

abline(h=swe.len.2010)
text(5.5,swe.len.2010+3,"2010")
abline(h=swe.len.2011)
text(5.5,swe.len.2011+3,"2011")
abline(h=swe.len.2012)
text(5.5,swe.len.2012+3,"2012")
dev.off()

####################################################################################
#MAX SWE AND PLOT
####################################################################################
max.swe.1980s <- maxswe.all.deltas(met.base,"1980s")
max.swe.rcp45.2040s<-maxswe.all.deltas(met.rcp45.2040s,"RCP4.5\n2040s")
max.swe.rcp45.2080s<-maxswe.all.deltas(met.rcp45.2080s,"RCP4.5\n2080s")
max.swe.rcp85.2040s<-maxswe.all.deltas(met.rcp85.2040s,"RCP8.5\n2040s")
max.swe.rcp85.2080s<-maxswe.all.deltas(met.rcp85.2080s,"RCP8.5\n2080s")
window.max.swe <- append(max.swe.1980s$group.name,c(max.swe.rcp45.2040s$group.name,max.swe.rcp85.2040s$group.name,max.swe.rcp45.2080s$group.name,max.swe.rcp85.2080s$group.name))
window.max.swe <- factor(window.max.swe,levels=c("1980s","RCP4.5\n2040s","RCP8.5\n2040s","RCP4.5\n2080s","RCP8.5\n2080s"))
max.swe <- append(max.swe.1980s$result,c(max.swe.rcp45.2040s$result,max.swe.rcp85.2040s$result,max.swe.rcp45.2080s$result,max.swe.rcp85.2080s$result))

png("/projects/daymet2/hare_sites/hare_wxTopo/figs/swemax_boxplt.png")
boxplot(max.swe~window.max.swe,col=c("royalblue1","orange","red","orange","red"),xlab="Time Period",ylab="cm",main="Modeled Annual Maximum Snow Water Equivalent",axes=FALSE,frame.plot=TRUE)
axis(2)
axis(1,tick=0,at=c(1,2,3,4,5),labels=levels(window.max.swe))
means.max <- tapply(max.swe,window.max.swe,mean)
points(means.max,pch=18,cex=1.5)

abline(h=max.swe.2010)
text(5.5,max.swe.2010+1,"2010")
abline(h=max.swe.2011)
text(5.5,max.swe.2011+1,"2011")
abline(h=max.swe.2012)
text(5.5,max.swe.2012+1,"2012")
dev.off()


#########################################
#SNOW SEASON LENGTH BY MODEL BOXPLOT (SUPPL. FIG.)
########################################
swe.len.base <- swe.season.all.deltas(met.base,"1980s","swe.len")

swe.len.rcp45.2040s.bymodel <- swe.season.all.deltas.bymodel(met.rcp45.2040s,"swe.len","2040s.45")
swe.len.rcp85.2040s.bymodel <- swe.season.all.deltas.bymodel(met.rcp85.2040s,"swe.len","2040s.85")
swe.len.rcp45.2080s.bymodel <- swe.season.all.deltas.bymodel(met.rcp45.2080s,"swe.len","2080s.45")
swe.len.rcp85.2080s.bymodel <- swe.season.all.deltas.bymodel(met.rcp85.2080s,"swe.len","2080s.85")

swe.len.2040s.all.bymodel = list(swe.len=c(),model.name=c())
swe.len.2040s.all.bymodel$swe.len <- append(swe.len.base$swe.len,c(swe.len.rcp45.2040s.bymodel$swe.len,swe.len.rcp85.2040s.bymodel$swe.len,swe.len.rcp45.2080s.bymodel$swe.len,swe.len.rcp85.2080s.bymodel$swe.len)) 
swe.len.2040s.all.bymodel$model.name <- append(rep("1980s",length(swe.len.base$swe.len)),c(swe.len.rcp45.2040s.bymodel$model.name,swe.len.rcp85.2040s.bymodel$model.name,swe.len.rcp45.2080s.bymodel$model.name,swe.len.rcp85.2080s.bymodel$model.name)) 

cols<-append("royalblue1",c(rep("orange",19),rep("red",19),rep("orange",19),rep("red",19)))

box.at <- seq(19*4)
box.at <- append(c(-1),box.at)

png("/projects/daymet2/hare_sites/hare_wxTopo/figs/length_boxplt_all.png",width=183,height=100,units="mm",res=300,pointsize=10)
boxplot(swe.len.2040s.all.bymodel$swe.len~swe.len.2040s.all.bymodel$model.name,at=box.at,col=cols,axes=FALSE,frame.plot=TRUE,xlab="Climate Model Ensemble Member",ylab="Days",outcex=.75,main="Modeled Main Snow Season Duration") #medlwd=1
axis(2)
axis(1,labels=FALSE,at=seq(76))
mtext(side=1,rep(LETTERS[1:19],4),at=seq(19*4),line=.5,cex=.6)
means <- tapply(swe.len.2040s.all.bymodel$swe.len,swe.len.2040s.all.bymodel$model.name,mean)
points(box.at,means,pch=18,cex=0.75)
abline(v=.5,lty=3)
abline(v=19.5,lty=3)
abline(v=38.5,lty=3)
abline(v=38.5+19,lty=3)
mtext(side=1,"1980s",at=-1.5,line=.2,cex=0.8)
text(6+.25,13,"RCP4.5 2040s",cex=0.9)
text(25+.25,13,"RCP8.5 2040s",cex=0.9)
text(44+.25,13,"RCP4.5 2080s",cex=0.9)
text(63+.25,13,"RCP8.5 2080s",cex=0.9)
segments(c(1.5-1),c(means.len[2]),c(20.5-1),c(means.len[2]),lwd=1)
segments(c(20.5-1),c(means.len[3]),c(39.5-1),c(means.len[3]),lwd=1)
segments(c(39.5-1),c(means.len[4]),c(58.5-1),c(means.len[4]),lwd=1)
segments(c(58.5-1),c(means.len[5]),c(78-1),c(means.len[5]),lwd=1)
mtext(side=3,"A",at=-1,line=.2,font=2)
dev.off()

#########################################
#MAX SWE BY MODEL BOXPLOT (SUPPL. FIG.)
########################################
max.swe.1980s <- maxswe.all.deltas(met.base,"1980s")

swe.max.rcp45.2040s.bymodel <- maxswe.all.deltas.bymodel(met.rcp45.2040s,"2040s.45")
swe.max.rcp85.2040s.bymodel <- maxswe.all.deltas.bymodel(met.rcp85.2040s,"2040s.85")
swe.max.rcp45.2080s.bymodel <- maxswe.all.deltas.bymodel(met.rcp45.2080s,"2080s.45")
swe.max.rcp85.2080s.bymodel <- maxswe.all.deltas.bymodel(met.rcp85.2080s,"2080s.85")

swe.max.all.bymodel = list(swe.max=c(),model.name=c())
swe.max.all.bymodel$swe.max <- append(max.swe.1980s$result,c(swe.max.rcp45.2040s.bymodel$swe.max,swe.max.rcp85.2040s.bymodel$swe.max,swe.max.rcp45.2080s.bymodel$swe.max,swe.max.rcp85.2080s.bymodel$swe.max)) 
swe.max.all.bymodel$model.name <- append(rep("1980s",length(max.swe.1980s$result)),c(swe.max.rcp45.2040s.bymodel$model.name,swe.max.rcp85.2040s.bymodel$model.name,swe.max.rcp45.2080s.bymodel$model.name,swe.max.rcp85.2080s.bymodel$model.name)) 

cols<-append("royalblue1",c(rep("orange",19),rep("red",19),rep("orange",19),rep("red",19)))

box.at <- seq(19*4)
box.at <- append(c(-1),box.at)
png("/projects/daymet2/hare_sites/hare_wxTopo/figs/maxswe_boxplt_all.png",width=183,height=100,units="mm",res=300,pointsize=10)
boxplot(swe.max.all.bymodel$swe.max~swe.max.all.bymodel$model.name,at=box.at,col=cols,axes=FALSE,frame.plot=TRUE,xlab="Climate Model Ensemble Member",ylab="cm",outcex=.75,main="Modeled Annual Maximum Snow Water Equivalent") #medlwd=1
axis(2)
axis(1,labels=FALSE,at=seq(76))
mtext(side=1,rep(LETTERS[1:19],4),at=seq(19*4),line=.5,cex=.6)
means <- tapply(swe.max.all.bymodel$swe.max,swe.max.all.bymodel$model.name,mean)
points(box.at,means,pch=18,cex=0.75)
abline(v=.5,lty=3)
abline(v=19.5,lty=3)
abline(v=38.5,lty=3)
abline(v=38.5+19,lty=3)
mtext(side=1,"1980s",at=-1.5,line=.2,cex=0.8)
text(6+.25,0,"RCP4.5 2040s",cex=0.9)
text(25+.25,0,"RCP8.5 2040s",cex=0.9)
text(44+.25,0,"RCP4.5 2080s",cex=0.9)
text(63+.25,0,"RCP8.5 2080s",cex=0.9)
segments(c(1.5-1),c(means.max[2]),c(20.5-1),c(means.max[2]),lwd=1)
segments(c(20.5-1),c(means.max[3]),c(39.5-1),c(means.max[3]),lwd=1)
segments(c(39.5-1),c(means.max[4]),c(58.5-1),c(means.max[4]),lwd=1)
segments(c(58.5-1),c(means.max[5]),c(78-1),c(means.max[5]),lwd=1)
mtext(side=3,"B",at=-1,line=.2,font=2)
dev.off()