###################################################################################################################
# This R script injests the various trait and stand-level measurements, as well as climate data, and computes site-level summaries. 
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)
require(reshape2) # for casting
require(plyr) # for ddply 
require(SDMTools) # for wt.mean and wt.sd
setwd("")

plot.cords <- read.csv('data/plot_details/plot_coordinates.csv', header=T)
plot.dims <- read.csv('data/plot_details/plot_dimensions.csv', header=T)
lai <- read.csv('data/plot_summaries/lai_by_plot.csv', header=T)
stems <- read.csv('data/stem_survey/tree_C_npp_sap.csv', header=T)
phys <- read.csv('data/stem_traits/stem_physiology_data.csv', header=T)
treeC.stems <- read.csv('data/stem_survey/tree_C_npp_sap.csv', header=T)
treeC <- read.csv('data/plot_summaries/treeC_by_plot.csv', header=T)
treeC.by.sp.subplot <- read.csv('data/subplot_summaries/subplot_sp_treeC.csv', header=T)
sap <- read.csv('data/plot_summaries/sapwood_by_plot.csv', header=T)
sg <- read.csv('data/stem_traits/stem_sg_procesed.csv', header=T)
leaf <- read.csv('data/leaf_traits/leaf_processed.csv', header = T)
rwc <- read.csv('data/stem_traits/stem_rwc_data.csv', header=T)
rw.clim.resp <- read.csv('outputs/rw_clim_relations/rw_clim_plot_frac_pos_resp_max.csv', header = T)
rw.clim.cor <- read.csv('outputs/rw_clim_relations/rw_clim_plot_top_cors.csv', header = T)
tree.age <- read.csv('data/treerings/tree_ages.csv', header = T)
clim.gy <- read.csv('outputs/climate_sum_annual_avgs_by_site.csv', header=T)

# This is the data frame into all the plot-level summaries will reside
plot.data <- subset(plot.cords, subplot == 1)
plot.data <- plot.data[,-c(3:5)]
plot.data$plot.id <- paste(plot.data$strata, plot.data$plot, sep='.')

# 
plot.data$sp.pch <- sort(rep(c(8,15,17),5))

#--------------------------------------------------------------------------------------------------------------
# set factors, etc
#--------------------------------------------------------------------------------------------------------------
spoi <- c('juoc','pipo','abgr')
strata <- c('WJ','PP','GF')
plots <- 1:5
subplots <- 1:4
species <- unique(stems$species)
yoi <- 1964:2013
n.yoi <- length(yoi)

stems$strata <- factor(stems$strata, levels=c('WJ','PP','GF'))
stems$plot <- factor(stems$plot, levels=1:5)
stems$plot.id <- paste(stems$strata, stems$plot, sep = '.')
stems$plot.id <- factor(stems$plot.id, levels = unique(stems$plot.id))
stems$subplot <- factor(stems$subplot)
stems$subplot.sp.id <- paste(stems$strata, stems$plot, stems$subplot, stems$species, sep='.')
stems$subplot.id <- paste(stems$strata, stems$plot, stems$subplot, sep = '.')
stems$subplot.id <- factor(stems$subplot.id, levels = unique(stems$subplot.id))
stems$sp.frac <- rep(NA, nrow(stems))

phys$strata <- factor(phys$strata, levels=c('WJ','PP','GF'))
phys$plot <- factor(phys$plot, levels=1:5)
phys$plot.id <- paste(phys$strata, phys$plot, sep = '.')
phys$plot.id <- factor(phys$plot.id, levels = unique(phys$plot.id))
phys$subplot <- factor(phys$subplot)
phys$subplot.sp.id <- paste(phys$strata, phys$plot, phys$subplot, phys$species, sep='.')
phys$subplot.id <- paste(phys$strata, phys$plot, phys$subplot, sep = '.')
phys$subplot.id <- factor(phys$subplot.id, levels = unique(phys$subplot.id))
phys$sp.frac <- rep(NA, nrow(phys))

leaf <- subset(leaf, aci.sample=='No')
leaf$strata <- factor(leaf$strata, levels=c('WJ','PP','GF'))
leaf$plot <- factor(leaf$plot, levels=1:5)
leaf$plot.id <- paste(leaf$strata, leaf$plot, sep = '.')
leaf$plot.id <- factor(leaf$plot.id, levels = unique(leaf$plot.id))
leaf$subplot <- factor(leaf$subplot)
leaf$subplot.sp.id <- paste(leaf$strata, leaf$plot, leaf$subplot, leaf$species, sep='.')
leaf$subplot.id <- paste(leaf$strata, leaf$plot, leaf$subplot, sep = '.')
leaf$subplot.id <- factor(leaf$subplot.id, levels = unique(leaf$subplot.id))
leaf$sp.frac <- rep(NA, nrow(leaf))

sg$strata <- factor(sg$strata, levels=c('WJ','PP','GF'))
sg$plot <- factor(sg$plot, levels=1:5)
sg$plot.id <- paste(sg$strata, sg$plot, sep = '.')
sg$plot.id <- factor(sg$plot.id, levels = unique(sg$plot.id))
sg$subplot <- factor(sg$subplot)
sg$subplot.sp.id <- paste(sg$strata, sg$plot, sg$subplot, sg$species, sep='.')
sg$subplot.id <- paste(sg$strata, sg$plot, sg$subplot, sep = '.')
sg$subplot.id <- factor(sg$subplot.id, levels = unique(sg$subplot.id))
sg$sp.frac <- rep(NA, nrow(sg))

rwc$strata <- factor(rwc$strata, levels=c('WJ','PP','GF'))
rwc$plot <- factor(rwc$plot, levels=1:5)
rwc$plot.id <- paste(rwc$strata, rwc$plot, sep = '.')
rwc$plot.id <- factor(rwc$plot.id, levels = unique(rwc$plot.id))
rwc$subplot <- factor(rwc$subplot)
rwc$subplot.sp.id <- paste(rwc$strata, rwc$plot, rwc$subplot, rwc$species, sep='.')
rwc$subplot.id <- paste(rwc$strata, rwc$plot, rwc$subplot, sep = '.')
rwc$subplot.id <- factor(rwc$subplot.id, levels = unique(rwc$subplot.id))
rwc$sp.frac <- rep(NA, nrow(rwc))

head(tree.age)
tree.age$strata <- toupper(tree.age$strata)
tree.age$strata <- factor(tree.age$strata, levels=c('WJ','PP','GF'))
tree.age$plot <- factor(tree.age$plot, levels=1:5)
tree.age$plot.id <- paste(tree.age$strata, tree.age$plot, sep = '.')
tree.age$plot.id <- factor(tree.age$plot.id, levels = unique(tree.age$plot.id))
tree.age$subplot <- factor(tree.age$subplot)
tree.age$subplot.sp.id <- paste(tree.age$strata, tree.age$plot, tree.age$subplot, tree.age$species, sep='.')
tree.age$subplot.id <- paste(tree.age$strata, tree.age$plot, tree.age$subplot, sep = '.')
tree.age$subplot.id <- factor(tree.age$subplot.id, levels = unique(tree.age$subplot.id))
tree.age$sp.frac <- rep(NA, nrow(tree.age))
tree.age <- tree.age[which(tree.age$species%in%spoi),]

treeC.by.sp.subplot$subplot.sp.id <- paste(treeC.by.sp.subplot$strata, treeC.by.sp.subplot$plot, treeC.by.sp.subplot$subplot, treeC.by.sp.subplot$species, sep='.')
treeC.by.sp.subplot$subplot.sp.id <- factor(treeC.by.sp.subplot$subplot.sp.id, levels = unique(treeC.by.sp.subplot$subplot.sp.id))

#--------------------------------------------------------------------------------------------------------------
# Add spcies biomass weight to data frames
#--------------------------------------------------------------------------------------------------------------
subplot.sp.id <- unique(stems$subplot.sp.id)
for (i in subplot.sp.id){
  sp.frac <- subset(treeC.by.sp.subplot, subplot.sp.id==i)
  
  stem.index <- which(stems$subplot.sp.id%in%i)
  stems$sp.frac[stem.index] <- rep(sp.frac$sp.frac, length(stem.index)) 
  
  phys.index <- which(phys$subplot.sp.id%in%i)
  phys$sp.frac[phys.index] <- rep(sp.frac$sp.frac, length(phys.index)) 
  
  leaf.index <- which(leaf$subplot.sp.id%in%i)
  leaf$sp.frac[leaf.index] <- rep(sp.frac$sp.frac, length(leaf.index)) 
  
  sg.index <- which(sg$subplot.sp.id%in%i)
  sg$sp.frac[sg.index] <- rep(sp.frac$sp.frac, length(sg.index)) 

  rwc.index <- which(rwc$subplot.sp.id%in%i)
  rwc$sp.frac[rwc.index] <- rep(sp.frac$sp.frac, length(rwc.index)) 
}

#--------------------------------------------------------------------------------------------------------------
# Stand C and ANPP
#--------------------------------------------------------------------------------------------------------------
plot.data$tree.kgCm2.2014.avg <- treeC$tree.kgCm2.2014.avg
plot.data$tree.kgCm2.2014.sd <- treeC$tree.kgCm2.2014.sd
plot.data$tree.kgCm2.2014.se <- treeC$tree.kgCm2.2014.sd/sqrt(4)

plot.data$tree.anpp.gCm2yr.avg <- treeC$tree.anpp.kgCm2yr.avg*1000
plot.data$tree.anpp.gCm2yr.sd <- treeC$tree.anpp.kgCm2yr.sd*1000
plot.data$tree.anpp.gCm2yr.se <- treeC$tree.anpp.kgCm2yr.sd/sqrt(4)*1000

#--------------------------------------------------------------------------------------------------------------
# Stand LAI
#--------------------------------------------------------------------------------------------------------------
plot.data$lai.hc.avg <- lai$lai.hc.avg
plot.data$lai.hc.sd <- lai$lai.hc.sd
plot.data$lai.hc.se <- lai$lai.hc.sd/sqrt(4)

#--------------------------------------------------------------------------------------------------------------
# Stand basal area (m2/ha)
#--------------------------------------------------------------------------------------------------------------
treeC.stems$strata <- factor(treeC.stems$strata,levels = c('WJ','PP','GF'))
treeC.stems$ba.cm2 <- 3.14159*(treeC.stems$dbh.cm/2)^2
stem.ba <- ddply(treeC.stems, .(strata, plot, subplot, size.class), plyr::summarize, ba.cm2.sum=sum(ba.cm2))
stem.ba$ba.m2 <- stem.ba$ba.cm2.sum/10^4

subplot.sizeclass.ba.m2ha <- as.data.frame(matrix(NA, nrow=0, ncol=5))
colnames(subplot.sizeclass.ba.m2ha) <- c('size.class','subplot','plot','strata','ba.m2ha')
subplot.sizeclass.ba.m2ha$strata <- factor(subplot.sizeclass.ba.m2ha$strata, levels=c('WJ','PP','GF'))

for (i in strata){
  for (j in plots){
    dims <- subset(plot.dims, strata == i & plot == j)
    ba.subplot <- subset(stem.ba, strata == i & plot == j)
    sapling.ba <- subset(ba.subplot, size.class=='sapling')
    mature.ba <- subset(ba.subplot, size.class=='mature')
    sapling.ba$ba.m2ha <-sapling.ba$ba.m2/dims$sapling.plot.area.m2*10^4
    mature.ba$ba.m2ha <-mature.ba$ba.m2/dims$stem.plot.area.m2*10^4
    subplot.sizeclass.ba.m2ha <- rbind(subplot.sizeclass.ba.m2ha, sapling.ba, mature.ba)
  }
}
subplot.ba.m2ha <- aggregate(ba.m2ha~subplot+plot+strata, subplot.sizeclass.ba.m2ha, sum)

plot.ba.m2ha  <- ddply(subplot.ba.m2ha, .(strata, plot), plyr::summarize, ba.m2ha.avg=round(mean(ba.m2ha),1),
                       ba.m2ha.sd=round(sd(ba.m2ha),1))
plot.ba.m2ha$ba.m2ha.se <- round(plot.ba.m2ha$ba.m2ha.sd/sqrt(4),1)

plot.data <- cbind(plot.data, plot.ba.m2ha)

#--------------------------------------------------------------------------------------------------------------
# tree density
#--------------------------------------------------------------------------------------------------------------
stem.trees <- subset(stems, dbh.cm>=10) # only take trees > 10 cm 
plot.dims.expanded <- plot.dims[sort(rep(seq(1,nrow(plot.dims)),4)),]
tree.dens.subplot <- ddply(stem.trees, .(strata,plot,subplot), plyr::summarize, n.trees=length(dbh.cm))
tree.dens.subplot$trees.ha <- tree.dens.subplot$n.trees/plot.dims.expanded$stem.plot.area.m2*10000 # tree/ha
tree.dens.plot <- ddply(tree.dens.subplot, .(strata,plot), plyr::summarize, 
                        trees.ha.avg=mean(trees.ha), trees.ha.sd=sd(trees.ha))
tree.dens.plot$trees.ha.se <- tree.dens.plot$trees.ha.sd/sqrt(4)
plot.data <- cbind(plot.data, tree.dens.plot[,3:5])

#--------------------------------------------------------------------------------------------------------------
# Mean AGB per tree
#--------------------------------------------------------------------------------------------------------------

tree.agb.subplot.avg <- ddply(stem.trees, .(strata,plot,subplot), plyr::summarize,tree.agb.kgC.avg = mean(tree.kgC.2014))
tree.agb.plot.avg <- ddply(tree.agb.subplot.avg, .(strata,plot), plyr::summarize, 
                      tree.avg.agb.kgC.avg = round(mean(tree.agb.kgC.avg),2), tree.avg.agb.cm.sd = round(sd(tree.agb.kgC.avg),2))
tree.agb.plot.avg$tree.avg.agb.kgC.se <- tree.agb.plot.avg$tree.avg.agb.cm.sd/sqrt(4)
plot.data <- cbind(plot.data, tree.agb.plot.avg[,3:5])

#--------------------------------------------------------------------------------------------------------------
# Max tree height
#--------------------------------------------------------------------------------------------------------------
# ht.subplot.wavg <- ddply(stems, "subplot.id", plyr::summarize, plot.id=unique(plot.id),
#                          ht.m.wavg = wt.mean(ht.total.m, sp.frac), ht.m.wsd = wt.sd(ht.total.m, sp.frac))
# ht.plot.wavg <- ddply(ht.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),
#                       ht.m.wavg = round(mean(ht.m.wavg),2), ht.m.wsd = round(sd(ht.m.wsd),2))
# ht.plot.wavg$ht.m.wse <- ht.plot.wavg$ht.m.wsd/sqrt(4)
# plot.data <- cbind(plot.data, ht.plot.wavg[,2:4])

ht.subplot.wavg <- ddply(stems, "subplot.id", plyr::summarize, plot.id=unique(plot.id),ht.m.max = max(ht.total.m))
ht.plot.wavg <- ddply(ht.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),
                      ht.m.wavg = round(mean(ht.m.max),2), ht.m.wsd = round(sd(ht.m.max),2))
ht.plot.wavg$ht.m.wse <- ht.plot.wavg$ht.m.wsd/sqrt(4)
plot.data <- cbind(plot.data, ht.plot.wavg[,2:4])

# ht.subplot.wavg <- ddply(stems, "subplot.id", plyr::summarize, plot.id=unique(plot.id),ht.m.max = mean(ht.total.m))
# ht.plot.wavg <- ddply(ht.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),
#                       ht.m.wavg = round(mean(ht.m.max),2), ht.m.wsd = round(sd(ht.m.max),2))
# ht.plot.wavg$ht.m.wse <- ht.plot.wavg$ht.m.wsd/sqrt(4)
# plot.data <- cbind(plot.data, ht.plot.wavg[,2:4])

#--------------------------------------------------------------------------------------------------------------
# Sapwood area 
#--------------------------------------------------------------------------------------------------------------
plot.data$sapwood.cm2.m2.wavg <- sap$sapwood.cm2m2
plot.data$sapwood.cm2.m2.wsd <- sap$sap.sd
plot.data$sapwood.cm2.m2.wse <- sap$sap.sd/sqrt(4)

plot.data$lai.2.sap.wavg <- plot.data$lai.hc.avg/plot.data$sapwood.cm2.m2.wavg
plot.data$lai.2.sap.wse <- with(plot.data, sqrt((sapwood.cm2.m2.wse/sapwood.cm2.m2.wavg)^2+(lai.hc.se/lai.hc.avg)^2)*lai.2.sap.wavg)

#--------------------------------------------------------------------------------------------------------------
# Stem specific gravity
#--------------------------------------------------------------------------------------------------------------
sg.subplot.wavg <- ddply(sg, "subplot.id", plyr::summarize, plot.id=unique(plot.id),
                         sg.g.cm3.wavg = wt.mean(sg.g.cm3, sp.frac), sg.g.cm3.wsd = wt.sd(sg.g.cm3, sp.frac))
sg.plot.wavg <- ddply(sg.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id), sg.g.cm3.wavg = round(mean(sg.g.cm3.wavg, na.rm = T),2))
sg.plot.wavg$sg.g.cm3.wsd <- ddply(sg.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),sg.g.cm3.wsd = round(sd(sg.g.cm3.wavg, na.rm = T),2))[,2]
sg.plot.wavg$sg.g.cm3.wse <- sg.plot.wavg$sg.g.cm3.wsd/sqrt(4)
plot.data <- cbind(plot.data, sg.plot.wavg[,2:4])

#--------------------------------------------------------------------------------------------------------------
# Stem relative water content
#--------------------------------------------------------------------------------------------------------------

rwc.subplot.wavg <- ddply(rwc, "subplot.id", plyr::summarize, plot.id=unique(plot.id),
                          rwc.pcnt.wavg = wt.mean(rwc.pcnt, sp.frac), rwc.pcnt.wsd = wt.sd(rwc.pcnt, sp.frac))
rwc.plot.wavg <- ddply(rwc.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id), rwc.pcnt.wavg = round(mean(rwc.pcnt.wavg, na.rm = T),2))
rwc.plot.wavg$rwc.pcnt.wsd <- ddply(rwc.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),rwc.pcnt.wsd = round(sd(rwc.pcnt.wavg, na.rm = T),2))[,2]
rwc.plot.wavg$rwc.pcnt.wse <- rwc.plot.wavg$rwc.pcnt.wsd/sqrt(4)
plot.data <- cbind(plot.data, rwc.plot.wavg[,2:4])


#--------------------------------------------------------------------------------------------------------------
# Stand age (90th quantile)
#--------------------------------------------------------------------------------------------------------------
age.q90 <- ddply(tree.age, .(strata,plot), plyr::summarize,  age.avg=mean(age), age.sd=sd(age), age.min=min(age),age.max=max(age), q90=quantile(age, 0.9), n=length(age)) # mean tree age by species
age.q90$plot.id <- with(age.q90, paste(strata,plot,sep='.'))
age.q90$q90avg <- rep(NA, nrow(age.q90))
age.q90$q90sd <- rep(NA, nrow(age.q90))
age.q90$plot.id <- paste(age.q90$strata, age.q90$plot, sep='.')

for (i in age.q90$plot.id){
  q90 <- subset(age.q90, plot.id==i)$q90
  plot.ages <- subset(tree.age, plot.id==i)
  oldies <- plot.ages[plot.ages$age>=q90,]
  age.q90$q90avg[cnt] <- round(mean(oldies$age))
  age.q90$q90sd[cnt] <- round(sd(oldies$age))
  cnt=cnt+1
}
age.q90
plot.data <- cbind(plot.data, data.frame(age.90q.avg=age.q90$q90avg))

#--------------------------------------------------------------------------------------------------------------
# Specific leaf area
#--------------------------------------------------------------------------------------------------------------
sla.subplot.wavg <- ddply(leaf, "subplot.id", plyr::summarize, plot.id=unique(plot.id),
                          sla.cm2.gC.wavg = wt.mean(sla.cm2.gC, sp.frac), sla.cm2.gC.wsd = wt.sd(sla.cm2.gC, sp.frac))
sla.plot.wavg <- ddply(sla.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),
                       sla.cm2.gC.wavg = round(mean(sla.cm2.gC.wavg),2))
sla.plot.wavg$sla.cm2.gC.wsd <- ddply(sla.subplot.wavg, "plot.id", plyr::summarize, sla.cm2.gC.wsd = round(sd(sla.cm2.gC.wavg),2))[,2]
sla.plot.wavg$sla.cm2.gC.wse <- sla.plot.wavg$sla.cm2.gC.wsd/sqrt(4)
plot.data <- cbind(plot.data, sla.plot.wavg[,2:4])

#--------------------------------------------------------------------------------------------------------------
# Foliar C
#--------------------------------------------------------------------------------------------------------------
leaf.C.subplot.wavg <- ddply(leaf, "subplot.id", plyr::summarize, plot.id=unique(plot.id),
                             leaf.C.g.g.wavg.sub = wt.mean(leaf.C.g.g, sp.frac), leaf.C.g.g.wsd.sub = wt.sd(leaf.C.g.g, sp.frac))
leaf.C.plot.wavg <- ddply(leaf.C.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),
                          leaf.C.g.g.wavg = round(mean(leaf.C.g.g.wavg.sub, na.rm = T),2), leaf.C.g.g.wsd=round(sd(leaf.C.g.g.wsd.sub, na.rm = T),2))
leaf.C.plot.wavg$leaf.C.g.g.wse <- leaf.C.plot.wavg$leaf.C.g.g.wsd/sqrt(4)
plot.data <- cbind(plot.data, leaf.C.plot.wavg[,2:4])

#--------------------------------------------------------------------------------------------------------------
# Foliar N
#--------------------------------------------------------------------------------------------------------------
#leaf <- leaf[-which(is.na(leaf$leaf.N.mg.g)),]
leaf.N.subplot.wavg <- ddply(leaf, "subplot.id", plyr::summarize, plot.id=unique(plot.id),
                             leaf.N.mg.g.wavg.sub = wt.mean(leaf.N.mg.g, sp.frac), leaf.N.mg.g.wsd.sub = wt.sd(leaf.N.mg.g, sp.frac))
leaf.N.plot.wavg <- ddply(leaf.N.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),
                          leaf.N.mg.g.wavg = round(mean(leaf.N.mg.g.wavg.sub, na.rm = T),2), leaf.N.mg.g.wsd=round(sd(leaf.N.mg.g.wsd.sub, na.rm = T),2))
leaf.N.plot.wavg$leaf.N.mg.g.wse <- leaf.N.plot.wavg$leaf.N.mg.g.wsd/sqrt(4)
plot.data <- cbind(plot.data, leaf.N.plot.wavg[,2:4])

#--------------------------------------------------------------------------------------------------------------
# Foliar C:N
#--------------------------------------------------------------------------------------------------------------
#leaf <- leaf[-which(is.na(leaf$leaf.N.mg.g)),]
leaf.CN.subplot.wavg <- ddply(leaf, "subplot.id", plyr::summarize, plot.id=unique(plot.id),
                              leaf.CN.g.g.wavg.sub = wt.mean(C.N.wt.ratio, sp.frac), leaf.N.mg.g.wsd.sub = wt.sd(C.N.wt.ratio, sp.frac))
leaf.CN.plot.wavg <- ddply(leaf.CN.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),
                           leaf.CN.g.g.wavg = round(mean(leaf.CN.g.g.wavg.sub, na.rm = T),2), leaf.CN.g.g.wsd=round(sd(leaf.CN.g.g.wavg.sub, na.rm = T),2))
leaf.CN.plot.wavg$leaf.CN.g.g.wse <- leaf.CN.plot.wavg$leaf.CN.g.g.wsd/sqrt(4)
plot.data <- cbind(plot.data, leaf.CN.plot.wavg[,2:4])

#--------------------------------------------------------------------------------------------------------------
# leaf longevity
#--------------------------------------------------------------------------------------------------------------

leaf.ll.subplot.wavg <- ddply(phys, "subplot.id", plyr::summarize, plot.id=unique(plot.id),
                              leaf.ll.yr.wavg.sub = wt.mean(ll.yr, sp.frac), leaf.ll.yr.wsd.sub = wt.sd(ll.yr, sp.frac))
leaf.ll.plot.wavg <- ddply(leaf.ll.subplot.wavg, "plot.id", plyr::summarize, plot.id=unique(plot.id),
                           leaf.ll.yr.wavg = round(mean(leaf.ll.yr.wavg.sub, na.rm = T),2), leaf.ll.yr.wsd=round(sd(leaf.ll.yr.wavg.sub, na.rm = T),2))
leaf.ll.plot.wavg$leaf.ll.yr.wse <- leaf.ll.plot.wavg$leaf.ll.yr.wsd/sqrt(4)
plot.data <- cbind(plot.data, leaf.ll.plot.wavg[,2:4])



#--------------------------------------------------------------------------------------------------------------
# Max fraction of population responding to spei
#--------------------------------------------------------------------------------------------------------------

rw.clim.resp$cor.pcnt <- rw.clim.resp$pos.cor.pcnt.avg
rw.clim.resp$cor.pcnt.sd <- rw.clim.resp$pos.cor.pcnt.sd
plot.data$rw.cmi.max.pos.cor.pcnt <- rep(NA, nrow(plot.data))
plot.data$rw.cmi.max.pos.cor.pcnt.sd <- rep(NA, nrow(plot.data))
rw.pcnt.resp.cmi <- subset(rw.clim.resp, clim.par=='cmi')
rw.pcnt.resp.cmi$plot.id <- paste(toupper(rw.pcnt.resp.cmi$strata), rw.pcnt.resp.cmi$plot, sep = '.')
plot.data$rw.cmi.max.pos.cor.pcnt[which(plot.data$plot.id%in%rw.pcnt.resp.cmi$plot.id)] <- rw.pcnt.resp.cmi$cor.pcnt
plot.data$rw.cmi.max.pos.cor.pcnt.sd[which(plot.data$plot.id%in%rw.pcnt.resp.cmi$plot.id)] <- rw.pcnt.resp.cmi$cor.pcnt.sd

#--------------------------------------------------------------------------------------------------------------
# Mean RWI-CMI correlation 
#--------------------------------------------------------------------------------------------------------------
rw.cor.cmi <- subset(rw.clim.cor, clim.par=='cmi')
rw.cor.cmi$plot.id <- paste(toupper(rw.cor.cmi$strata), rw.cor.cmi$plot, sep = '.')
plot.data$rw.cmi.cor.avg[which(plot.data$plot.id%in%rw.cor.cmi$plot.id)] <- rw.cor.cmi$cor.avg
plot.data$rw.cmi.cor.sd[which(plot.data$plot.id%in%rw.cor.cmi$plot.id)] <- rw.cor.cmi$cor.sd
plot.data$rw.cmi.cor.se[which(plot.data$plot.id%in%rw.cor.cmi$plot.id)] <- rw.cor.cmi$cor.sd/sqrt(20)

#--------------------------------------------------------------------------------------------------------------
# climate data
#--------------------------------------------------------------------------------------------------------------
head(clim.gy)
tmean <- subset(clim.gy, clim.par=='tmean')
plot.data$tmean.gy.avg <- tmean$clim.avg
plot.data$tmean.gy.sd <-  tmean$clim.sd
plot.data$tmean.gy.se <-  tmean$clim.sd/sqrt(4)

ppt <- subset(clim.gy, clim.par=='ppt')
plot.data$ppt.gy.avg <- ppt$clim.avg
plot.data$ppt.gy.sd <-  ppt$clim.sd
plot.data$ppt.gy.se <-  ppt$clim.sd/sqrt(4)

vpd <- subset(clim.gy, clim.par=='vpd')
plot.data$vpd.gy.avg <- vpd$clim.avg
plot.data$vpd.gy.sd <-  vpd$clim.sd
plot.data$vpd.gy.se <-  vpd$clim.sd/sqrt(4)

eto <- subset(clim.gy, clim.par=='eto')
plot.data$eto.gy.avg <- eto$clim.avg
plot.data$eto.gy.sd <-  eto$clim.sd
plot.data$eto.gy.se <-  eto$clim.sd/sqrt(4)

cmi <- subset(clim.gy, clim.par=='cmi')
plot.data$cmi.gy.avg <- cmi$clim.avg
plot.data$cmi.gy.sd <-  cmi$clim.sd
plot.data$cmi.gy.se <-  cmi$clim.sd/sqrt(4)


#--------------------------------------------------------------------------------------------------------------
# output plot summaries 
#--------------------------------------------------------------------------------------------------------------
write.table(plot.data, 'data/plot_summaries/plot_master_data_all_sites.csv', sep = ',', row.names = F, col.names = T)

drop.plots <- c('WJ.4','PP.4','GF.3') # drop sites with young trees 
plot.data <- plot.data[which(plot.data$plot.id%in%drop.plots==F),]

write.table(plot.data, 'data/plot_summaries/plot_master_data.csv', sep = ',', row.names = F, col.names = T)

