###################################################################################################################
# This R script creates figure 6 from Berner and Law (2015)
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)
require(reshape2)
require(plotrix)
require(nlme)
setwd("")
source('fun_fit_models.R')

plot.data <- read.csv('data/plot_summaries/plot_master_data.csv', header = T)

#--------------------------------------------------------------------------------------------
# Correlations and regression models
#--------------------------------------------------------------------------------------------
ht.lm <- lm(ht.m.wavg~cmi.gy.avg, subset(plot.data, plot.id != 'GF.1'))
#ht.lm <- lm(ht.m.wavg~cmi.gy.avg, plot.data)
ht.lm.summary <- summary(ht.lm)
ht.lm.summary

ht.power <- fit.power(x = plot.data$cmi.gy.avg, y = plot.data$ht.m.wavg)
ht.power$coefs

sg.lm <- lm(sg.g.cm3.wavg~cmi.gy.avg, plot.data)
sg.lm.summary <- summary(sg.lm)
sg.lm.summary

lai2sap.lm <- lm(lai.2.sap.wavg~cmi.gy.avg, plot.data)
lai2sap.lm.summary <- summary(lai2sap.lm)
lai2sap.lm.summary

#--------------------------------------------------------------------------------------------
# Plot 
#--------------------------------------------------------------------------------------------
pdf('figures/fig6_anatomy_by_cmi.pdf', 4,8)
#jpeg('figures/fig6_anatomy_by_cmi.jpeg', 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))


# wood density
plotCI(plot.data$cmi.gy.avg, plot.data$sg.g.cm3.wavg, uiw=plot.data$cmi.gy.se, err='x',  
       xlim=c(-1300,250), ylim=c(0.36,0.51),xlab = '', ylab='',cex = my.cex.axis,  pch=plot.data$sp.pch, scol='darkgrey', xaxt='n', yaxt='n')
plotCI(plot.data$cmi.gy.avg, plot.data$sg.g.cm3.wavg, uiw=plot.data$sg.g.cm3.wse, err='y', add=T,
       xlim=c(-1300,250), ylim=c(0.36,0.51),xlab = '', ylab='',cex = my.cex.axis,  pch=plot.data$sp.pch, scol='darkgrey', xaxt='n', yaxt='n')
mtext(side = 2, line = 2.5, cex = my.pch.cex, expression(paste("Stem wood density (g cm" ^-3, ')')))
axis(1, at = seq(-1250,250,500), labels = T, cex.axis=my.cex.axis)
axis(2, at = seq(0.40,0.55,0.05), labels = T, cex.axis=my.cex.axis)
text(-1300,0.505, '(a)', cex=1.3, font=2)
text(100, 0.505, expression('r' ^2~'= 0.36'),cex=my.cex.axis)
lines(plot.data$cmi.gy.avg, fitted(sg.lm), lwd=1.5, col='gray80')
sg.lm.summary

# tree height
plotCI(plot.data$cmi.gy.avg, plot.data$ht.m.wavg, uiw=plot.data$cmi.gy.se, err='x',  
       xlim=c(-1300,250), ylim=c(0,50),xlab = '', ylab='',cex = my.cex.axis,  pch=plot.data$sp.pch, scol='darkgrey', xaxt='n', yaxt='n')
plotCI(plot.data$cmi.gy.avg, plot.data$ht.m.wavg, uiw=plot.data$ht.m.wse, err='y', add=T,  
       xlim=c(-1300,250), ylim=c(0,50),xlab = '', ylab='',cex = my.cex.axis,  pch=plot.data$sp.pch, scol='darkgrey', xaxt='n', yaxt='n')
mtext(side = 2, line = 2.5, cex = my.pch.cex, "Max tree height (m)")
text(-1300,48.5, '(b)', cex=1.3, font=2)
axis(2, at = seq(0,50,10), labels = T, cex.axis=my.cex.axis)
axis(1, at = seq(-1250,250,500), labels = T, cex.axis=my.cex.axis)
text(100, 48.5, expression('r' ^2~'= 0.69'), cex=my.cex.axis)
#lines(plot.data$cmi.gy.avg[-9], fitted(ht.lm), col='gray80')
lines(ht.power$preds$x, ht.power$preds$modeled, lwd=1.5, col='gray80')


# leaf to sapwood area ratio 
plotCI(plot.data$cmi.gy.avg, plot.data$lai.2.sap.wavg, uiw=plot.data$cmi.gy.se, err='x',  
       xlim=c(-1300,250), ylim=c(0,0.23),xlab = '', ylab='',cex = my.cex.axis,  pch=plot.data$sp.pch, scol='darkgrey', xaxt='n', yaxt='n')
plotCI(plot.data$cmi.gy.avg, plot.data$lai.2.sap.wavg, uiw=plot.data$lai.2.sap.wse, err='y', add=T,  
       xlim=c(0,0.3), ylim=c(0,0.23),xlab = '', ylab='',cex = my.cex.axis,  pch=plot.data$sp.pch, scol='darkgrey', xaxt='n', yaxt='n')
mtext(side = 2, line = 2.5, cex = my.pch.cex, expression(paste('Leaf : sapwood area (m' ^2, ' cm ' ^-2,')')))
text(-1300,0.22, '(c)', cex=1.3, font=2)
#legend(x = -1300, y = 0.35, legend = c('Western juniper','Ponderosa pine','Grand fir'), pch = unique(plot.data$sp.pch), bty='n', ncol=1, cex=1.3)
legend('bottomright', legend = c('Western juniper','Ponderosa pine','Grand fir'), pch = unique(plot.data$sp.pch), bty='n', ncol=1, cex=1.3)
axis(1, at = seq(-1250,250,500), labels = T, cex.axis=my.cex.axis)
axis(2, at = seq(0,0.4,0.1), labels = T, cex.axis=my.cex.axis)
text(100, 0.22, expression('r' ^2~'= 0.76'),cex=my.cex.axis)
lines(plot.data$cmi.gy.avg, fitted(lai2sap.lm), lwd=1.5, col='gray80')
lai2sap.lm.summary
mtext(side = 1, line = 3, cex = my.pch.cex, my.xlab)

dev.off()

