###################################################################################################################
# This R script creates figure 3 from Berner and Law (2015)
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)
require(reshape2)
require(plotrix)
setwd("")
load('workspaces/forest_by_cmi_nls.RDATA')
plot.data <- read.csv('data/plot_summaries/plot_master_data.csv', header = T)

#--------------------------------------------------------------------------------------------
# Plot mean site tree C stock, NPP, and LAI by CMIgy
#--------------------------------------------------------------------------------------------
pdf('figures/fig3_site_treeC_anpp_lai_by_cmigy.pdf', 4,8)
#jpeg('figures/fig3_site_treeC_anpp_lai_by_cmigy.jpg', 4,8, units = 'in', res = 500)
my.xlab  <- expression('Climate moisture index (PPT-ET'[0]*' ; mm yr'^-1*')')
my.pch.cex=1
my.cex.axis=1.2
layout(matrix(c(1,2,3),3,1, byrow = TRUE))
par.top <- par(mar=c(4,4.75,1,1))
plotCI(plot.data$cmi.gy.avg, plot.data$lai.hc.avg, uiw=plot.data$lai.hc.se, err='y', xaxt='n', yaxt='n', 
       xlim=c(-1300,250), ylim=c(0,7.2),xlab = '', ylab='',cex = my.cex.axis, pch=plot.data$sp.pch, scol='darkgrey')
plotCI(plot.data$cmi.gy.avg, plot.data$lai.hc.avg, uiw=plot.data$cmi.gy.se, err='x', add=T, xaxt='n', yaxt='n',   
       xlim=c(-1300,250), ylim=c(0,7.2),xlab = '', ylab='',cex = my.cex.axis, pch=plot.data$sp.pch, scol='darkgrey')
lines(lai.cmi.nls$preds$x, lai.cmi.nls$preds$modeled, lwd=1.5, col='gray80')
mtext(side = 2, line = 2.9, cex = my.pch.cex, expression(paste('Leaf area index (m' ^2, ' m' ^-2, ')')))
legend(x = -400, y=2, legend = c('Western juniper','Ponderosa pine','Grand fir'), pch = unique(plot.data$sp.pch), bty='n', ncol=1, cex=1.3)
axis(2, seq(0,8,2), las=2, cex.axis=my.cex.axis)
axis(1, at = seq(-1250,250,500), cex.axis=my.cex.axis, labels = T)
text(-1250,6.5,"(a)", cex=1.3, font=2)
text(100,6.5, expression('r' ^2~'= 0.68'), cex=my.cex.axis)

plotCI(plot.data$cmi.gy.avg, plot.data$tree.anpp.gCm2yr.avg, uiw=plot.data$tree.anpp.gCm2yr.se, err='y', xaxt='n',  yaxt='n', 
       xlim=c(-1300,250), ylim=c(0,380),xlab = '', ylab='',cex = my.cex.axis, pch=plot.data$sp.pch, scol='darkgrey')
plotCI(plot.data$cmi.gy.avg, plot.data$tree.anpp.gCm2yr.avg, uiw=plot.data$cmi.gy.se, err='x', add=T, xaxt='n',  yaxt='n', 
       xlim=c(-1300,250), ylim=c(0,380),xlab = '', ylab='',cex = my.cex.axis, pch=plot.data$sp.pch, scol='darkgrey')
lines(npp.cmi.nls$preds$x, npp.cmi.nls$preds$modeled, lwd=1.5, col='gray80')
mtext(side = 2, line = 2.9, cex = my.pch.cex, expression(paste('Tree ANPP (g C m' ^-2, ' yr' ^-1,')')))
axis(1, at = seq(-1250,250,500), labels = T, cex.axis=my.cex.axis)
axis(2, seq(0,375,75), las=2, cex.axis=my.cex.axis)
text(-1250,350,"(b)", cex=1.3, font=2)
text(100,350, expression('r' ^2~'= 0.88'), cex=my.cex.axis)

plotCI(plot.data$cmi.gy.avg, plot.data$tree.kgCm2.2014.avg, uiw=plot.data$tree.kgCm2.2014.se, err='y',  yaxt='n', 
       xlim=c(-1300,250), ylim=c(0,25),xlab = '', ylab='',cex = my.cex.axis, pch=plot.data$sp.pch, scol='darkgrey', xaxt='n')
plotCI(plot.data$cmi.gy.avg, plot.data$tree.kgCm2.2014.avg, uiw=plot.data$cmi.gy.se, err='x', add=T,  yaxt='n', 
       xlim=c(-1300,250), ylim=c(0,25),xlab = '', ylab='',cex = my.cex.axis, pch=plot.data$sp.pch, scol='darkgrey', xaxt='n')
lines(agb.cmi.nls$preds$x, agb.cmi.nls$preds$modeled, lwd=1.5, col='gray80')
mtext(side = 2, line = 2.9, cex = my.pch.cex, expression(paste('Tree AGB (kg C m' ^-2, ')')))
text(-1250,23,"(c)", cex=1.3, font=2)
axis(1, at = seq(-1250,250,500), labels = T, cex.axis=my.cex.axis)
axis(2, seq(0,25,5), las=2, cex.axis=my.cex.axis)
mtext(side = 1, line=3, cex=1, my.xlab)
text(100,23, expression('r' ^2~'= 0.67'), cex=my.cex.axis)

dev.off()
