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

#------------------------------------------------------------------------------
# lattice error bar
# custom function for adding error bands to a lattice xyplot
# from: http://www.r-bloggers.com/confidence-bands-with-lattice-and-r/

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 = 'gray70', border = FALSE,
                ...)
}
#------------------------------------------------------------------------------


clim <- read.csv('data/plot_summaries/climate_by_plot.csv', header=T)
clim <- subset(clim, year>=1963)
clim$month <- rep(1:12, nrow(clim)/12)
clim$strata <- factor(clim$strata, levels=c('WJ','PP','GF'), labels = c('Western juniper','Ponderosa pine','Grand fir'))
head(clim)
yoi <- 1964:2013
n.yoi <- length(yoi)

# ----------------------------------------------------------------------
# calculate growing-year values
# ----------------------------------------------------------------------
index.drop.onepar <- (c(1:9,(((n.yoi+1)*12)-2):((n.yoi+1)*12)))
index.drop.full <- c()
for (i in seq(0,nrow(clim)-612,612)){
  index.drop.full <- c(index.drop.full, index.drop.onepar+i)
}
clim.gy <- clim[-index.drop.full,]
clim.gy$year.gy <- rep(sort(rep(yoi,12)),15)
drop.plots <- c('WJ.4','PP.4','GF.3')
clim.gy <- clim.gy[which(clim.gy$plot.id%in%drop.plots==F),]

cmi <- subset(clim.gy, clim.par=='cmi')
cmi.gy <- ddply(cmi, .(plot,strata,year.gy), plyr::summarize, cmi.gy.mm=sum(values))

cmi.gy.avg <- ddply(cmi.gy, .(strata,year.gy), plyr::summarize, cmi.gy.avg=mean(cmi.gy.mm, na.rm=T), 
                    cmi.gy.sd=sd(cmi.gy.mm,na.rm=T), n.obs=length(cmi.gy.mm))
cmi.gy.avg$cmi.gy.se <- with(cmi.gy.avg,cmi.gy.sd/sqrt(n.obs))
cmi.gy.avg$cmi.upper <- with(cmi.gy.avg, cmi.gy.avg+cmi.gy.se)
cmi.gy.avg$cmi.lower <- with(cmi.gy.avg, cmi.gy.avg-cmi.gy.se)
#cmi.gy.avg$strata <- factor(cmi.gy.avg$strata, levels = c('juoc','pipo','abgr'), labels = c('Western juniper','Ponderosa pine','Grand fir'))


write.table(cmi.gy.avg, 'outputs/climate_cmigy_timeseries_by_strata.csv', sep = ',', row.names = F, col.names = T)

# ----------------------------------------------------------------------
# plot CMI time series 
# ----------------------------------------------------------------------

#pdf('figures/fig2_cmigy_timeseries.pdf', 4.95,6)

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

clim.ts <- xyplot(cmi.gy.avg~year.gy|strata, cmi.gy.avg, groups = strata, ylab=my.ylab,xlab='Year', layout=c(1,3), 
       upper = cmi.gy.avg$cmi.upper, lower = cmi.gy.avg$cmi.lower, 
       scales=list(y=list(at=seq(-1500,1000,500), alternating=F), x=list(at=seq(1965,2010,5), rot=90)),
       panel = function(x, y, ...){
         panel.superpose(x, y, panel.groups = my.panel.bands, type='l',...)
         panel.xyplot(x,y, type='l', lwd=1, col='black')
         panel.abline(v=c(1973,1977,1992,1994,2001,2003), lty = 3, lwd=0.7)
       }
)

#------------------------------------------------------------------------------

rw <- read.csv('data/treerings/treering_master_data.csv', header = T)
rw <- subset(rw, confidence == 1)
rw <- subset(rw, reliable <= 1964)
rw <- subset(rw, year >= 1964)
rw$plot.id <- toupper(paste(rw$strata, rw$plot, sep='.'))
drop.plots <- c('WJ.4','PP.4','GF.3')
rw <- rw[which(rw$plot.id%in%drop.plots==F),]


rw.by.yr <- ddply(rw, .(species,year), plyr::summarize, rw.avg=mean(rw.mm, na.rm=T), 
                  rw.sd=sd(rw.mm,na.rm=T), n.obs=length(rw.mm))
rw.by.yr$rw.se <- with(rw.by.yr,rw.sd/sqrt(n.obs))
rw.by.yr$rw.upper <- with(rw.by.yr, rw.avg+rw.se)
rw.by.yr$rw.lower <- with(rw.by.yr, rw.avg-rw.se)
rw.by.yr$species <- factor(rw.by.yr$species, levels = c('juoc','pipo','abgr'), labels = c('Western juniper','Ponderosa pine','Grand fir'))




rw.ts <- xyplot(rw.avg~year|species, data=rw.by.yr, groups = species, ylab='Annual tree ring-width (mm)',xlab='Year', layout=c(1,3), 
       upper = rw.by.yr$rw.upper, lower = rw.by.yr$rw.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', ...)
         panel.xyplot(x,y, type='l', lwd=1, col='black')
         panel.abline(v=c(1973,1977,1992,1994,2001,2003), lty = 3, lwd=0.7)
       }
)

pdf('figures/fig2_cmi_rw_timeseries.pdf', width = 6.65, height=5)
jpeg('figures/fig2_cmi_rw_timeseries.jpg', width =18, height =15, units='cm', res=400)
sb <- trellis.par.get("strip.background") 
sb[["col"]][1] <- "lightgray"
trellis.par.set("strip.background", sb) 
print(clim.ts, position=c(0,0,0.5,1), more=T)
print(rw.ts, position=c(0.5,0,1,1))
dev.off()