###################################################################################################################
# This R script creates figure S2 from Berner and Law (2015)
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(plyr)
setwd("")

#----------------------------------------------------------------------------------
my.panel.bands <- function(x, y, upper, lower, fill, col,
                           subscripts, ..., font, fontface)
{
  upper <- upper[subscripts]
  lower <- lower[subscripts]
  panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
                col = 'gray60', border = FALSE,
                ...)
}
#----------------------------------------------------------------------------------

strata.water.gy <- read.csv('outputs/climate_water_growingyear_timeseries_by_strata.csv', header = T)
jpeg('figures/figS1_cmigy_anom_timeseries.jpg', 4.95,6, units = 'in', res = 300)

sb <- trellis.par.get("strip.background") 
sb[["col"]][1] <- "lightgray"
trellis.par.set("strip.background", sb) 
my.ylab <- expression('Climate moisture index anomaly (PPT-ET'[0]*' ; mm yr'^-1*')')

# plot anomalies
strata.water.gy$strata <- factor(strata.water.gy$strata, levels=c('WJ','PP','GF'), labels = c('Western juniper','Ponderosa pine','Grand fir'))
strata.cmi.gy <- subset(strata.water.gy, clim.par == 'cmi')

xyplot(clim.anom.avg~year.gy|strata, strata.cmi.gy, groups=strata, layout=c(1,3), ylab = my.ylab, xlab= 'Growing year',
       upper = strata.cmi.gy$clim.anom.upper, lower = strata.cmi.gy$clim.anom.lower, 
       scales=list(x=list(at=seq(1965,2010,5), rot=90)),
       panel = function(x, y, ...){
         panel.superpose(x, y, panel.groups = my.panel.bands, type='l', col='gray',...)
         panel.xyplot(x, y, type=c('b'), cex=0.6, lty=1,lwd=1,col='black', ...)
         panel.abline(v=c(1985,1994), lty = 3, lwd=0.7)
         panel.abline(h=c(0), lty = 2, lwd=0.7)}
)

dev.off()
