###################################################################################################################
# This R script summarizes trait measurements by species
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)
require(latticeExtra)
require(plyr)
require(lme4)
require(lsmeans)

#setwd("C:/research/projects/cascade_tree_traits/")
setwd("C:/Users/bernerl/Google Drive/research/cascade_tree_traits/")
#setwd("C:/Users/Logan/Google Drive/research/cascade_tree_traits/")
sg <- read.csv('data/stem_traits/stem_sg_procesed.csv', header=T)
leaf <- read.csv('data/leaf_traits/leaf_processed.csv', header = T)
stems <- read.csv('data/stem_survey/stem_survey_data.csv', header = T)
npp <- read.csv('data/stem_survey/tree_C_npp_sap.csv', header = T)
sap <- read.csv('data/stem_traits/stem_sapwood_processed.csv', header = T)
rwc <- read.csv('data/stem_traits/stem_rwc_data.csv', header = T)
spoi <- c('juoc','pipo','abgr')
drop.sites <- c('WJ.4','PP.4','GF.3')
#-----------------------------------------------------------------------------------------------------
# Factorize and all that jazz
#-----------------------------------------------------------------------------------------------------
stems$strata <- factor(stems$strata, levels=c('WJ','PP','GF'))
stems$plot <- factor(stems$plot, levels=1:5)
stems$subplot <- factor(stems$subplot)
stems$plot.id <- with(stems, paste(strata,plot,sep='.'))
stems <- stems[-which(stems$plot.id%in%drop.sites),]
stems <- stems[which(stems$species%in%spoi),]

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$subplot <- factor(leaf$subplot)
leaf$plot.id <- with(leaf, paste(strata,plot,sep='.'))
leaf <- leaf[-which(leaf$plot.id%in%drop.sites),]
leaf <- leaf[which(leaf$species%in%spoi),]

sap$strata <- factor(sap$strata, levels=c('WJ','PP','GF'))
sap$plot <- factor(sap$plot, levels=1:5)
sap$subplot <- factor(sap$subplot)
sap$plot.id <- with(sap, paste(strata,plot,sep='.'))
sap <- sap[-which(sap$plot.id%in%drop.sites),]
sap <- sap[which(sap$species%in%spoi),]

rwc$strata <- factor(rwc$strata, levels=c('WJ','PP','GF'))
rwc$plot <- factor(rwc$plot, levels=1:5)
rwc$subplot <- factor(rwc$subplot)
rwc$plot.id <- with(rwc, paste(strata,plot,sep='.'))
rwc <- rwc[-which(rwc$plot.id%in%drop.sites),]
rwc <- rwc[which(rwc$species%in%spoi),]

sg$strata <- factor(sg$strata, levels=c('WJ','PP','GF'))
sg$plot <- factor(sg$plot, levels=1:5)
sg$subplot <- factor(sg$subplot)
sg$plot.id <- with(sg, paste(strata,plot,sep='.'))
sg <- sg[-which(sg$plot.id%in%drop.sites),]
sg <- sg[which(sg$species%in%spoi),]

npp$strata <- factor(npp$strata, levels=c('WJ','PP','GF'))
npp$plot <- factor(npp$plot, levels=1:5)
npp$subplot <- factor(npp$subplot)
npp$plot.id <- with(npp, paste(strata,plot,sep='.'))
npp <- npp[-which(npp$plot.id%in%drop.sites),]
npp <- npp[which(npp$species%in%spoi),]


#--------------------------------------------------------------------------------------------------------------
# SLA
#--------------------------------------------------------------------------------------------------------------
sla.subplot.avg <- ddply(leaf, .(strata,plot,subplot,species), plyr::summarize, plot.id=unique(plot.id),
                            sla.cm2.gC.avg.sub = mean(sla.cm2.gC), sla.cm2.gC.sd.sub = sd(sla.cm2.gC))
sla.plot.avg <- ddply(sla.subplot.avg, .(strata,plot,species), plyr::summarize,
                         sla.cm2.gC.avg = round(mean(sla.cm2.gC.avg.sub, na.rm = T),2), 
                         sla.cm2.gC.sd=round(sd(sla.cm2.gC.sd.sub, na.rm = T),2))
sla.sp.avg <- ddply(sla.plot.avg, .(species), plyr::summarize,
                       sla.cm2.gC.avg = round(mean(sla.cm2.gC.avg, na.rm = T)), 
                       sla.cm2.gC.sd=round(sd(sla.cm2.gC.sd, na.rm = T)))
#sla.sp.avg$sla.cm2.gC.se <- sla.sp.avg$sla.cm2.gC.sd/sqrt(4)
sla.sp.avg$n <- ddply(leaf, .(species), plyr::summarize, n=length(sla.cm2.gC))[,2]

#--------------------------------------------------------------------------------------------------------------
# Foliar N
#--------------------------------------------------------------------------------------------------------------
leaf$leaf.N.pcnt <- leaf$leaf.N.mg.g/10
leaf.N.subplot.avg <- ddply(leaf, .(strata,plot,subplot,species), plyr::summarize, plot.id=unique(plot.id),
                             leaf.N.pcnt.avg.sub = mean(leaf.N.pcnt), leaf.N.pcnt.sd.sub = sd(leaf.N.pcnt))
leaf.N.plot.avg <- ddply(leaf.N.subplot.avg, .(strata,plot,species), plyr::summarize,
                         leaf.N.pcnt.avg = round(mean(leaf.N.pcnt.avg.sub, na.rm = T),2), 
                         leaf.N.pcnt.sd=round(sd(leaf.N.pcnt.sd.sub, na.rm = T),2))
leaf.N.sp.avg <- ddply(leaf.N.plot.avg, .(species), plyr::summarize,
                       leaf.N.pcnt.avg = round(mean(leaf.N.pcnt.avg, na.rm = T),2), 
                       leaf.N.pcnt.sd=round(sd(leaf.N.pcnt.sd, na.rm = T),2))
#leaf.N.sp.avg$leaf.N.pcnt.se <- leaf.N.sp.avg$leaf.N.pcnt.sd/sqrt(4)
leaf.N.sp.avg$n <- ddply(leaf, .(species), plyr::summarize, n=length(leaf.N.pcnt))[,2]
#--------------------------------------------------------------------------------------------------------------
# Foliar C
#--------------------------------------------------------------------------------------------------------------
leaf$leaf.C.pcnt <- leaf$leaf.C.g.g
leaf.C.subplot.avg <- ddply(leaf, .(strata,plot,subplot,species), plyr::summarize, plot.id=unique(plot.id),
                            leaf.C.pcnt.avg.sub = mean(leaf.C.pcnt), leaf.C.pcnt.sd.sub = sd(leaf.C.pcnt))
leaf.C.plot.avg <- ddply(leaf.C.subplot.avg, .(strata,plot,species), plyr::summarize,
                         leaf.C.pcnt.avg = mean(leaf.C.pcnt.avg.sub, na.rm = T), 
                         leaf.C.pcnt.sd=sd(leaf.C.pcnt.sd.sub, na.rm = T))
leaf.C.sp.avg <- ddply(leaf.C.plot.avg, .(species), plyr::summarize,
                       leaf.C.pcnt.avg = round(mean(leaf.C.pcnt.avg, na.rm = T),1), 
                       leaf.C.pcnt.sd=round(sd(leaf.C.pcnt.sd, na.rm = T),1))
#leaf.C.sp.avg$leaf.C.pcnt.se <- round(leaf.C.sp.avg$leaf.C.pcnt.sd/sqrt(4),1)
leaf.C.sp.avg$n <- ddply(leaf, .(species), plyr::summarize, n=length(leaf.C.pcnt))[,2]
#--------------------------------------------------------------------------------------------------------------
# Foliar CN
#--------------------------------------------------------------------------------------------------------------
leaf.CN.subplot.avg <- ddply(leaf, .(strata,plot,subplot,species), plyr::summarize,
                            C.N.wt.ratio.avg.sub = mean(C.N.wt.ratio), C.N.wt.ratio.sd.sub = sd(C.N.wt.ratio))
leaf.CN.plot.avg <- ddply(leaf.CN.subplot.avg, .(strata,plot,species), plyr::summarize,
                         C.N.wt.ratio.avg = mean(C.N.wt.ratio.avg.sub, na.rm = T), 
                         C.N.wt.ratio.sd=sd(C.N.wt.ratio.sd.sub, na.rm = T))
leaf.CN.sp.avg <- ddply(leaf.CN.plot.avg, .(species), plyr::summarize,
                       C.N.wt.ratio.avg = round(mean(C.N.wt.ratio.avg, na.rm = T),1), 
                       C.N.wt.ratio.sd=round(sd(C.N.wt.ratio.sd, na.rm = T),1))
#leaf.CN.sp.avg$C.N.wt.ratio.se <- round(leaf.CN.sp.avg$C.N.wt.ratio.sd/sqrt(4),1)
leaf.CN.sp.avg$n <- ddply(leaf, .(species), plyr::summarize, n=length(C.N.wt.ratio))[,2]
#--------------------------------------------------------------------------------------------------------------
# Stem specific gravity
#--------------------------------------------------------------------------------------------------------------
sg.subplot.avg <- ddply(sg, .(strata,plot,subplot,species), plyr::summarize, 
                         sg.g.cm3.avg = mean(sg.g.cm3, na.rm = T), sg.g.cm3.sd = sd(sg.g.cm3,na.rm = T))
sg.plot.avg <- ddply(sg.subplot.avg, .(strata,plot,species), plyr::summarize, 
                     sg.g.cm3.avg = round(mean(sg.g.cm3.avg, na.rm = T),2))
sg.sp.avg <- ddply(sg.plot.avg, .(species), plyr::summarize, 
                   sg.g.cm3.avg = round(mean(sg.g.cm3.avg, na.rm = T),2))
sg.sp.avg$sg.g.cm3.sd <- ddply(sg.plot.avg, .(species), plyr::summarize, sg.g.cm3.sd = round(sd(sg.g.cm3.avg,na.rm = T),2))[,2]
#sg.sp.avg$sg.g.cm3.se <- sg.sp.avg$sg.g.cm3.sd/sqrt(4)
sg.sp.avg$n <- ddply(sg, .(species), plyr::summarize, n=length(sg.g.cm3))[,2]
#--------------------------------------------------------------------------------------------------------------
# Max height
#--------------------------------------------------------------------------------------------------------------
ht.subplot.avg <- ddply(stems, .(strata,plot,subplot,species), plyr::summarize, 
                        ht.max.m = max(ht.total.m, na.rm = T))
ht.plot.avg <- ddply(ht.subplot.avg, .(strata,plot,species), plyr::summarize, 
                     ht.max.m.avg = round(mean(ht.max.m, na.rm = T),1))
ht.sp.avg <- ddply(ht.plot.avg, .(species), plyr::summarize, 
                   ht.max.m.avg = round(mean(ht.max.m.avg, na.rm = T),1))
ht.sp.avg$ht.max.m.sd <- ddply(ht.plot.avg, .(species), plyr::summarize, ht.max.m.sd = round(sd(ht.max.m.avg,na.rm = T),1))[,2]
#ht.sp.avg$ht.max.m.se <- ht.sp.avg$ht.max.m.sd/sqrt(4)
ht.sp.avg$n <- ddply(stems, .(species), plyr::summarize, n=length(ht.total.m))[,2]

#--------------------------------------------------------------------------------------------------------------
# Sapwood
#--------------------------------------------------------------------------------------------------------------
bwplot(sapwood.2.ba~species,sap, notch=T)
bwplot(water.wt.pcnt~species,rwc, notch=T,ylim=c(30,80))
sap.subplot.avg <- ddply(npp, .(strata,plot,subplot,species), plyr::summarize, 
                         sap.2.ba = mean(sap.2.ba.cm2.m2, na.rm = T))
sap.plot.avg <- ddply(sap.subplot.avg, .(strata,plot,species), plyr::summarize, 
                      sap.2.ba.avg = round(mean(sap.2.ba, na.rm = T),1))
sap.sp.avg <- ddply(sap.plot.avg, .(species), plyr::summarize, 
                    sap.kgC.yr.avg = round(mean(sap.kgC.yr.avg, na.rm = T),1))
sap.sp.avg$sap.kgC.yr.sd <- ddply(sap.plot.avg, .(species), plyr::summarize, sap.kgC.yr.sd = round(sd(sap.kgC.yr.avg,na.rm = T),1))[,2]
#sap.sp.avg$sap.kgC.yr.se <- sap.sp.avg$sap.kgC.yr.sd/sqrt(4)
sap.sp.avg$n <- ddply(sap, .(species), plyr::summarize, n=length(tree.10yr.avg.asap.kgCyr))[,2]
sap.sp.avg


#--------------------------------------------------------------------------------------------------------------
# Compile species averages
#--------------------------------------------------------------------------------------------------------------
fancy <- function(x){paste(x[,2],"±",x[,3]," (",x[,4],")",sep='')}
fancy2 <- function(x){paste(x[,2],"±",x[,3],sep='')}


traits.avgs <- as.data.frame(matrix(NA, nrow = 3,ncol = 8))
colnames(traits.avgs) <- c('Species','SG','H','NPP','SLA','C','N','CN')
traits.avgs$Species <- c('abgr','juoc','pipo')

traits.avgs$SG <- fancy(sg.sp.avg)
traits.avgs$H <- fancy(ht.sp.avg)
traits.avgs$NPP <- fancy(npp.sp.avg)
traits.avgs$SLA <- fancy2(sla.sp.avg)
traits.avgs$C <- fancy2(leaf.C.sp.avg)
traits.avgs$N <- fancy2(leaf.N.sp.avg)
traits.avgs$CN <- fancy(leaf.CN.sp.avg)

write.table(traits.avgs, 'outputs/trait_species_avg_sd.csv', sep = ',', col.names = T, row.names = F)

