###################################################################################################################
# This R script computes AGB, ANPP, and sapwood area for every subplot and plot in the study.
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)
require(dplyr)
setwd("")

strata <- c('WJ','PP','GF')
plots <- 1:5
subplots <- 1:4

plot.dims <- read.csv('data/plot_details/plot_dimensions.csv', header=T)
stems <- read.csv('data/stem_survey/tree_C_npp_sap.csv', header=T)
sg <- read.csv('data/stem_traits/stem_sg_procesed.csv', header=T)
phys <- read.csv('data/stem_traits/stem_physiology_data.csv', header=T, fill=T)
lai <- read.csv('data/lai2200/lai_processed.csv', header = T)
leaf <- read.csv('data/leaf_traits/leaf_processed.csv', header = T)

#--------------------------------------------------------------------------------------------------------------
# set factors, build data frames
#--------------------------------------------------------------------------------------------------------------
stems$strata <- factor(stems$strata, levels=c('WJ','PP','GF'))
sg$strata <- factor(sg$strata, levels=c('WJ','PP','GF'))
sg$plot <- factor(sg$plot)
sg$subplot <- factor(sg$subplot)
species <- unique(stems$species)
sg.sp <- subset(sg, strata == "WJ" & species == "juoc" | strata == "PP" & species == "pipo" | strata == "GF" & species == "abgr")

#--------------------------------------------------------------------------------------------------------------
# calculate average leaf longevity and sla for each species
#--------------------------------------------------------------------------------------------------------------
sp.leaf.traits <- ddply(phys, .(species), plyr::summarize, ll.yr.avg = mean(ll.yr, na.rm = T))
sp.leaf.traits$ll.yr.avg[sp.leaf.traits$species == 'juoc'] <- 5.6 # leaf longevity for juniperous Runyon et al. (1994) Eco. Apps.
sp.leaf.traits$ll.yr.avg[sp.leaf.traits$species == 'cade'] <- 8.9 # leaf longevity for thuja plicata from Harlow et al. (2005) 
sp.leaf.traits$sla.cm2.gC.avg <-  ddply(leaf, .(species), plyr::summarize, sla.cm3.gC.avg = mean(sla.cm2.gC, na.rm = T))$sla.cm3.gC.avg
sp.leaf.traits$leaf.C.g.g.avg <-  ddply(leaf, .(species), plyr::summarize, leaf.C.g.g.avg = mean(leaf.C.g.g, na.rm = T))$leaf.C.g.g.avg

#--------------------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------------------
# PLOT/SUBPLOT TREE CARBON STOCK CALCULATIONS
#--------------------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------------------

#--------------------------------------------------------------------------------------------------------------
# # caluclate tree C stocks (kgC/m2) and ANPP (kgC/m2/yr) for both mature trees and saplings by species and subplot  (calc seperately because of differences in plot size)  
#--------------------------------------------------------------------------------------------------------------
subplot.agb.sp.size <- aggregate(tree.kgC.2014~species+size.class+subplot+plot+strata, stems, sum)
subplot.anpp.sp.size <- aggregate(tree.10yr.avg.wood.npp.kgCyr~species+size.class+subplot+plot+strata, stems, sum)

subplot.agb.sp.size.kgCm2 <- as.data.frame(matrix(NA, nrow=0, ncol=9))
colnames(subplot.agb.sp.size.kgCm2) <- c('species','size.class','subplot','plot','strata',
                                         'tree.kgC.2014','tree.kgCm2.2014','tree.wood.anpp.kgCyr.10yr.avg','tree.wood.anpp.kgCm2yr.10yr.avg')
# i='GF'
# j=1
for (i in strata){
  for (j in plots){
    dims <- subset(plot.dims, strata == i & plot == j)
    agb <-  subset(subplot.agb.sp.size, strata == i & plot == j)
    agb.sapling <- subset(agb, size.class=='sapling')
    agb.mature <- subset(agb, size.class=='mature')
    agb.sapling$tree.kgCm2.2014 <- agb.sapling$tree.kgC.2014/dims$sapling.plot.area.m2
    agb.mature$tree.kgCm2.2014 <- agb.mature$tree.kgC.2014/dims$stem.plot.area.m2
    
    npp <-  subset(subplot.anpp.sp.size, strata == i & plot == j)
    npp.sapling <- subset(npp, size.class=='sapling')
    npp.mature <- subset(npp, size.class=='mature')
    npp.sapling$tree.wood.anpp.kgCm2yr.10yr.avg <- npp.sapling$tree.10yr.avg.wood.npp.kgCyr/dims$sapling.plot.area.m2
    npp.mature$tree.wood.anpp.kgCm2yr.10yr.avg <- npp.mature$tree.10yr.avg.wood.npp.kgCyr/dims$stem.plot.area.m2
    subplot.agb.sp.size.kgCm2 <- rbind(subplot.agb.sp.size.kgCm2, cbind(agb.sapling, npp.sapling[,6:7]), cbind(agb.mature, npp.mature[,6:7]))
  }
}

#--------------------------------------------------------------------------------
# calculate tree C and fraction of total tree C for each species in each subplot. Also calculate foliage npp 
#--------------------------------------------------------------------------------
# sum sapling and mature tree AGB to get total AGB for each sp in each subplot
subplot.agb.sp.kgCm2 <- aggregate(tree.kgCm2.2014~species+subplot+plot+strata, subplot.agb.sp.size.kgCm2, sum)
subplot.agb.sp.kgCm2$tree.wood.anpp.kgCm2yr <- aggregate(tree.wood.anpp.kgCm2yr.10yr.avg~species+subplot+plot+strata, subplot.agb.sp.size.kgCm2, sum)[,5]
subplot.agb.sp.kgCm2$subplot.id <- paste(subplot.agb.sp.kgCm2$strata, subplot.agb.sp.kgCm2$plot, subplot.agb.sp.kgCm2$subplot, sep='.')
subplot.agb.sp.kgCm2$subplot.id <- factor(subplot.agb.sp.kgCm2$subplot.id, levels=unique(subplot.agb.sp.kgCm2$subplot.id))
sp.subplot.frac <- NA

for (i in unique(subplot.agb.sp.kgCm2$subplot.id)){
  i.agb <- subset(subplot.agb.sp.kgCm2, subplot.id==i)
  subplot.sp.frac <- i.agb$tree.kgCm2.2014/sum(i.agb$tree.kgCm2.2014)
  sp.subplot.frac <- c(sp.subplot.frac, subplot.sp.frac)
}
subplot.agb.sp.kgCm2$sp.frac.C <- as.vector(na.omit(sp.subplot.frac))

write.table(subplot.agb.sp.kgCm2, 'data/subplot_summaries/treeC_by_subplot_sp.csv', sep=',', row.names=F, col.names=T)
#barchart(tree.kgCm2.2014~subplot.id, subplot.agb.sp.kgCm2, groups=species, stack=T, ylab='Aboveground carbon stock (kgC/m2)', auto.key=T)
#barchart(sp.frac~subplot.id, subplot.agb.sp.kgCm2, groups=species, stack=T, ylab='Species fraction', auto.key=T)

#--------------------------------------------------------------------------------
# calculate total tree C and ANPP for each subplot
#--------------------------------------------------------------------------------
subplot.agb.sp.kgCm2$ll.avg.yrs <- sp.leaf.traits$ll.yr.avg[match(subplot.agb.sp.kgCm2$species, sp.leaf.traits$species)]
subplot.agb.sp.kgCm2$sla.cm2.gC.avg <- sp.leaf.traits$sla.cm2.gC.avg[match(subplot.agb.sp.kgCm2$species, sp.leaf.traits$species)]
subplot.agb.sp.kgCm2$lai.hc <- lai$lai.hc[match(subplot.agb.sp.kgCm2$subplot.id, lai$subplot.id)]

subplot.tree.kgCm2 <- ddply(subplot.agb.sp.kgCm2, .(strata, plot, subplot), plyr::summarize, tree.wood.kgCm2.2014=sum(tree.kgCm2.2014),
                            tree.leaf.kgCm2=sum((lai.hc/(sla.cm2.gC.avg/10))*sp.frac.C),
                            tree.wood.anpp.kgCm2yr=sum(tree.wood.anpp.kgCm2yr), 
                            tree.leaf.anpp.kgCm2yr=sum((lai.hc/(sla.cm2.gC.avg/10)/ll.avg.yrs)*sp.frac.C))
subplot.tree.kgCm2$tree.kgCm2.2014 <- subplot.tree.kgCm2$tree.wood.kgCm2.2014+subplot.tree.kgCm2$tree.leaf.kgCm2
subplot.tree.kgCm2$tree.anpp.kgCm2yr <- subplot.tree.kgCm2$tree.wood.anpp.kgCm2yr+subplot.tree.kgCm2$tree.leaf.anpp.kgCm2yr
head(subplot.tree.kgCm2)
write.table(subplot.tree.kgCm2, 'data/subplot_summaries/treeC_by_subplot.csv', sep=',', row.names=F, col.names=T)

plot.biomass.prod <-  subplot.tree.kgCm2 %>% group_by(strata,plot) %>% summarise_each(funs(mean,sd))
write.table(plot.biomass.prod, 'data/plot_summaries/tree_biomass_prod_by_subplot.csv', sep=',', row.names=F, col.names=T)

#--------------------------------------------------------------------------------
# calculate total tree C and ANPP for each species in each plot
#--------------------------------------------------------------------------------
plot.agb.sp.kgCm2 <- aggregate(tree.kgCm2.2014~species+plot+strata, subplot.agb.sp.kgCm2, sum) # total AGB summed across all subplots (ie. not avg of subplots)
plot.agb.sp.kgCm2$tree.anpp.kgCm2yr <- aggregate(tree.wood.anpp.kgCm2yr~species+plot+strata, subplot.agb.sp.kgCm2, sum)[,4]
colnames(plot.agb.sp.kgCm2) <- c('species','plot','strata','tree.kgCm2.2014.total','tree.anpp.kgCm2yr.total')
plot.agb.sp.kgCm2$plot.id <- paste(plot.agb.sp.kgCm2$strata, plot.agb.sp.kgCm2$plot, sep='.')
sp.plot.frac <- NA
#i='GF.4'
for (i in unique(plot.agb.sp.kgCm2$plot.id)){
  i.agb <- subset(plot.agb.sp.kgCm2, plot.id==i)
  sp.plot.frac <- c(sp.plot.frac, round(i.agb$tree.kgCm2.2014.total/sum(i.agb$tree.kgCm2.2014.total),3))
}
plot.agb.sp.kgCm2$sp.frac <- na.omit(sp.plot.frac)
plot.agb.sp.kgCm2$tree.kgCm2.2014.avg <- aggregate(tree.kgCm2.2014~species+plot+strata, subplot.agb.sp.kgCm2, mean)[,4]
plot.agb.sp.kgCm2$tree.kgCm2.2014.sd <- aggregate(tree.kgCm2.2014~species+plot+strata, subplot.agb.sp.kgCm2, sd)[,4]

write.table(plot.agb.sp.kgCm2, 'data/plot_summaries/treeC_by_plot_sp.csv', sep=',', row.names=F, col.names=T)

#--------------------------------------------------------------------------------------------------------------
# calculate mean/sd tree C for each plot 
#--------------------------------------------------------------------------------------------------------------
plot.tree.kgCm2.avg <- aggregate(tree.kgCm2.2014~plot+strata, subplot.tree.kgCm2, mean)
plot.tree.kgCm2.avg$sd <- aggregate(tree.kgCm2.2014~plot+strata, subplot.tree.kgCm2, sd)[,3]
colnames(plot.tree.kgCm2.avg) <- c('plot','strata','tree.kgCm2.2014.avg','tree.kgCm2.2014.sd')
plot.tree.kgCm2.avg$tree.anpp.kgCm2yr.avg <- aggregate(tree.anpp.kgCm2yr~plot+strata, subplot.tree.kgCm2, mean)[,3]
plot.tree.kgCm2.avg$tree.anpp.kgCm2yr.sd <- aggregate(tree.anpp.kgCm2yr~plot+strata, subplot.tree.kgCm2, sd)[,3]
plot.data <- plot.tree.kgCm2.avg
write.table(plot.tree.kgCm2.avg, 'data/plot_summaries/treeC_by_plot.csv', sep=',', row.names=F, col.names=T)

#--------------------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------------------
# SAPWOOD CALCULATIONS
#--------------------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------------------

#--------------------------------------------------------------------------------------------------------------
# calculate total sapwood area (cm2/m2) for each subplot
#--------------------------------------------------------------------------------------------------------------
plot.sap <- aggregate(sapwood.cm2~size.class+subplot+plot+strata, stems, sum)
plot.sap.cm2m2 <- as.data.frame(matrix(NA, nrow=0, ncol=5))
colnames(plot.sap.cm2m2) <- c('size.class','subplot','plot','strata','sapwood.cm2m2')

# normalize subplot sapwood area by plot area
for (i in strata){
  for (j in plots){
    dims <- subset(plot.dims, strata == i & plot == j)
    sap <-  subset(plot.sap, strata == i & plot == j)
    sap.sapling <- subset(sap, size.class=='sapling')
    sap.mature <- subset(sap, size.class=='mature')
    sap.sapling$sapwood.cm2m2 <-sap.sapling$sapwood.cm2/dims$sapling.plot.area.m2
    sap.mature$sapwood.cm2m2 <-sap.mature$sapwood.cm2/dims$stem.plot.area.m2
    plot.sap.cm2m2 <- rbind(plot.sap.cm2m2, sap.sapling, sap.mature)
  }
}

subplot.sap.cm2m2 <- aggregate(sapwood.cm2m2~subplot+plot+strata, plot.sap.cm2m2, sum)
write.table(subplot.sap.cm2m2, 'data/subplot_summaries/sapwood_by_subplot.csv', sep=',', row.names=F, col.names=T)

#--------------------------------------------------------------------------------------------------------------
# calculate average sapwood area (cm2/m2) for each plot
#--------------------------------------------------------------------------------------------------------------
plot.sap.cm2m2.avg <- aggregate(sapwood.cm2m2~plot+strata, subplot.sap.cm2m2, mean)
plot.sap.cm2m2.avg$sap.sd <- aggregate(sapwood.cm2m2~plot+strata, subplot.sap.cm2m2, sd)[,3]
write.table(plot.sap.cm2m2.avg, 'data/plot_summaries/sapwood_by_plot.csv', sep=',', row.names=F, col.names=T)

#--------------------------------------------------------------------------------------------------------------
# END SCRIPT
#--------------------------------------------------------------------------------------------------------------
