###################################################################################################################
# This R script computes 10-yr growth increments for each tree, develops linear models relating both tree height
# and sapwood area to DBH, and then esitmates AGB, ANPP, and sapwood area for every tree in the study.
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################

rm(list=ls())
require(lattice)
require(plyr)

setwd("")

#-------------------------------------------------------------------------
# read in data files
sap <- read.csv('data/stem_traits/stem_sapwood_processed.csv', header=T, fill=T)
phys <- read.csv('data/stem_traits/stem_physiology_data.csv', header=T, fill=T)
stems <- read.csv('data/stem_survey/stem_survey_data.csv', header=T, fill=T)
sg <- read.csv('data/stem_traits/stem_sg_procesed.csv', header=T, fill=T)
allom <- read.csv('data/terrapnw_data/BiomassEquations_ECascades.csv', fill=T)
rw <- read.csv('data/treerings/treering_master_data.csv', header=T, fill=T)

species <- unique(stems$species)
n.trees <- nrow(stems)
strata <- unique(stems$strata)
plots <- 1:5

#-------------------------------------------------------------------------
# Calculate wood density for each site
#-------------------------------------------------------------------------
sp.sg.avg <- aggregate(sg.g.cm3~species, sg, mean)

#-------------------------------------------------------------------------
# Estimate 10 year mean growth increment for each tree by grouping trees into quantiles and applying
# measured increment to non-measured increments, like done by Hudiburg et al. ()
#-----------------------------------------------------------------------
stems$dbh.change.10yr.cm <- rep(NA, nrow(stems))
for (i in strata){
  for (j in 1:length(plots)){
    plot.stems <- subset(stems, strata==i & plot==j)
    plot.sp <- unique(plot.stems$species)
    for (k in plot.sp){
      sp.stems <- subset(plot.stems, species==k) # get species stem in plot
      sp.stems.cored <- sp.stems[tolower(sp.stems$id)%in%rw$id,] # id the trees that were cored
      # if no cores for a species were taken on a plot, then raise a flag and apply growth increment from other species 
      if (nrow(sp.stems.cored)==0){
        next()
      } else {
        sp.rw <- rw[rw$id%in%tolower(sp.stems.cored$id),] # get ringwidth data
        # for the subdominant species, group all trees together into the same size class
        if (nrow(sp.stems.cored)<=5){
          sp.stems$quartile <- 1
          sp.stems.cored$quartile <- 1
        } else {
          brks <- quantile(sp.stems$dbh.cm, probs = seq(0,1,1/3)) # identify dbh quartile (actually 1/3s)
          sp.stems <- within(sp.stems, quartile <- cut(dbh.cm, breaks = brks, labels = 1:3, include.lowest = TRUE)) # add quartile category 
          sp.stems.cored <- within(sp.stems.cored, quartile <- cut(dbh.cm, breaks = brks, labels = 1:3, include.lowest = TRUE)) # add quartile category 
        }
        sp.rw.10yrs <- subset(sp.rw, year>=2004) # subset ringwidth data to the past 10 years
        sp.rw.10yr.sum <- aggregate(rw.mm~id, sp.rw.10yrs, sum) # calculate total radial growth over time period
        sp.rw.10yr.sum$dbh.dif.10yr.cm <- (sp.rw.10yr.sum$rw.mm*2)/10 # convert change in radius to change in diameter, converting also from mm to cm
        sp.rw.10yr.sum$quartile <- sp.stems.cored$quartile
        sp.rw.10yr.sum.quart.avg <- aggregate(dbh.dif.10yr.cm~quartile, sp.rw.10yr.sum, mean, na.rm=T)
        sp.stems <- merge(sp.stems,sp.rw.10yr.sum.quart.avg, by = "quartile")
        stems$dbh.change.10yr.cm[stems$id%in%sp.stems$id] <- sp.stems$dbh.dif.10yr.cm
      }
    }
  }
}

#stems[is.na(stems$dbh.change.10yr.cm),]
stems$dbh.2004.cm <- stems$dbh.cm-stems$dbh.change.10yr.cm

#-------------------------------------------------------------------------
# Model tree ht ~ dbh | species using a Weibull-type function described in Yang et al. (1978) and Huang et al. (1992) both in CJFR. 
#-------------------------------------------------------------------------
sp.oi <- c('juoc','pipo','abgr')

# build df to store regression coefs
ht.nlms <- as.data.frame(matrix(NA, ncol=10,nrow=length(sp.oi)))
colnames(ht.nlms) <- c('species','a','a.se','b','b.se','c','c.se','r2','rmse','n')

cnt=1
for (i in sp.oi){
  sp.dat <- subset(stems, species == i)
  my.nlm <- nls(ht.total.m~1.3+exp(a+b/(dbh.cm+c)), start=list(a=3, b=-10, c=0),sp.dat) # Weibull-type function of Yang et al. (1978)
  #my.nlm <- nls(ht.total.m~1.3+a*(1-exp(-b*dbh.cm^c)), start=list(a=30, b=0, c=1),sp.dat)
  my.nlm.summary <- summary(my.nlm)
  
  resids <- my.nlm.summary$residuals
  resid.ss <- sum((sp.dat$ht.total.m-predict(my.nlm))^2)
  total.ss <- (nrow(sp.dat)-1)*var(sp.dat$ht.total.m)
  r2 <- round(1-resid.ss/total.ss,2)
  rmse  <- round(sqrt(mean(resids^2)),2)
  ht.nlms$species[cnt] <- i
  ht.nlms$a[cnt] <- round(my.nlm.summary$coefficients[1,1],3)
  ht.nlms$a.se[cnt] <- round(my.nlm.summary$coefficients[1,2],3)
  ht.nlms$b[cnt] <- round(my.nlm.summary$coefficients[2,1],3)
  ht.nlms$b.se[cnt] <- round(my.nlm.summary$coefficients[2,2],3)
  ht.nlms$c[cnt] <- round(my.nlm.summary$coefficients[3,1],3)
  ht.nlms$c.se[cnt] <- round(my.nlm.summary$coefficients[3,2],3)
  ht.nlms$r2[cnt] <- r2
  ht.nlms$rmse[cnt] <- rmse
  ht.nlms$n[cnt] <- nrow(sp.dat)
  assign(paste('ht',i,'nlm', sep='.'), value=my.nlm)
  cnt=cnt+1
}
ht.nlms

# Estimate height for all trees (for obs vs pred reasons) 
for (i in 1:n.trees){
  i.tree <- stems[i,]
  i.sp <- i.tree$species
  if (i.sp=='abpr'){
    i.nlm <- get("ht.abgr.nlm")
  } else if (i.sp == 'cade'){
    i.nlm <- get("ht.abgr.nlm")
  } else if (i.sp == 'pico'){
    i.nlm <- get("ht.pipo.nlm")
  } else {
    i.nlm <- get(paste("ht",i.sp,'nlm', sep='.'))
  }
  stems$ht.total.2014.m.nlm[i] <- predict(i.nlm, newdata=data.frame(dbh.cm=i.tree$dbh.cm))
  stems$ht.total.2004.m.nlm[i] <- predict(i.nlm, newdata=data.frame(dbh.cm=i.tree$dbh.2004.cm)) # estimate tree height in 2004 based on dbh growth rates
}
stems$ht.total.2004.m.corrected <- with(stems, ht.total.m-(ht.total.2014.m.nlm-ht.total.2004.m.nlm))
# plot ht~dbh relationships
pdf('figures/species_ht_dbh_scatterplots.pdf', 8,6)
xyplot(ht.total.m~dbh.cm|species, stems, type='p', ylab='Tree height (m)', xlab='Tree DBH (cm)')
dev.off()

# plot predicted vs observed tree heights
pdf('figures/species_ht_dbh_scatterplots.pdf', 8,6)
xyplot(ht.total.2004.m.nlm~ht.total.m|species, stems, ylab='Predicted tree height (m)', xlab='Measured tree height (m)',
       panel = function(x,y,...) {
         panel.abline(a=0, b = 1)
         panel.xyplot(x,y)
       })
dev.off()

# output ht by dbh regression stats
write.table(ht.nlms, 'outputs/species_ht_dbh_allometry.csv', sep=',', col.names=T, row.names=F)


#-------------------------------------------------------------------------
# Model tree AGB and NPP for each tree (and component)
#-----------------------------------------------------------------------
stems$size.class[stems$dbh.cm<10] <- 'sapling'
stems$size.class[stems$dbh.cm>=10] <- 'mature'

component.C <- as.data.frame(matrix(NA,nrow=0,ncol=8))
colnames(component.C) <- c('strata','plot','subplot','tree.id','species','component','agb.gC.2014','agb.gC.2004')

# i=3
# j=1
for (i in 1:n.trees){
  i.tree <- stems[i,]
  allom.sp <- subset(allom, species == as.character(i.tree$species))
  
  #-----------------------------------------------
  # seedlings (dbh < 10 cm)
  if (i.tree$dbh.cm < 10){
    # 2014 sizes
    allom.sapling <- subset(allom.sp, eqn.type=="SB")
    eqn.expression <- as.character(gsub('Y', i.tree$ht.total.m, as.character(gsub('X', i.tree$dbh.cm, as.character(allom.sapling$eqn)))))
    agb.2014 <- eval(parse(text=eqn.expression))
    agb.gC.2014 <- agb.2014*0.51
    
    # 2004 sizes
    eqn.expression <- as.character(gsub('Y', i.tree$ht.total.2004.m.corrected, as.character(gsub('X', i.tree$dbh.2004.cm, as.character(allom.sapling$eqn)))))
    agb.2004 <- eval(parse(text=eqn.expression))
    agb.gC.2004 <- agb.2004*0.51

    # store seedling data
    component.C <- rbind(component.C, data.frame(strata=i.tree$strata, plot=i.tree$plot, subplot=i.tree$subplot,
                                             tree.id=i.tree$tree.id, species=i.tree$species, component='SB', 
                                             agb.gC.2014=agb.gC.2014, agb.gC.2004=agb.gC.2004))
  } 
  
  #-----------------------------------------------
  # non-seedlings
  else {
    allom.sp <- subset(allom.sp, eqn.type != 'SB')
    
    for (j in 1:nrow(allom.sp)){
      # model tree C for 2014 tree size
      eqn.expression <- as.character(gsub('Y', i.tree$ht.total.m, as.character(gsub('X', i.tree$dbh.cm, as.character(allom.sp$eqn[j])))))
      agb.2014 <- eval(parse(text=eqn.expression))
      # convert bole volume from m3 to g DM
      if (allom.sp$eqn.type[j] == "BV"){
        sp.sg <- subset(sp.sg.avg, species==i.tree$species)
        agb.2014 <- agb.2014*(10^6)*sp.sg$sg.g.cm3
      }
      # conver dry matter to carbon, assuming 51%C for wood and 48%C for foliage
      if (allom.sp$eqn.type[j] == "FB"){
        agb.gC.2014 <- agb.2014*0.48 # convert g to kg and assume DM is 48% C 
      } else {
        agb.gC.2014 <- agb.2014*0.52 # convert g to kg and assume DM is 52% C - Lamlon and Savidge (2003) Biomass and Bioenergy
      }

      # model tree C for 2004 tree size
      eqn.expression <- as.character(gsub('Y', i.tree$ht.total.2004.m.corrected, as.character(gsub('X', i.tree$dbh.2004.cm, as.character(allom.sp$eqn[j])))))
      agb.2004 <- eval(parse(text=eqn.expression))
      # convert bole volume from m3 to g DM
      if (allom.sp$eqn.type[j] == "BV"){
        sp.sg <- subset(sp.sg.avg, species==i.tree$species)
        agb.2004 <- agb.2004*(10^6)*sp.sg$sg.g.cm3
      }
      # conver dry matter to carbon, assuming 51%C for wood and 48%C for foliage
      if (allom.sp$eqn.type[j] == "FB"){
        agb.gC.2004 <- agb.2004*0.48 # convert g to kg and assume DM is 48% C 
      } else {
        agb.gC.2004 <- agb.2004*0.52 # convert g to kg and assume DM is 52% C 
      }      
        
        # store non-seedling data 
        component.C <- rbind(component.C, data.frame(strata=i.tree$strata, plot=i.tree$plot, subplot=i.tree$subplot,
                                             tree.id=i.tree$tree.id, species=i.tree$species, component=allom.sp$eqn.type[j], 
                                             agb.gC.2014=agb.gC.2014, agb.gC.2004=agb.gC.2004))
    }
  }
}


#-------------------------------------------------------------------------
# Aggregate component C mass to total tree C for 2014 and 2004  
#-------------------------------------------------------------------------
component.C$strata <- factor(component.C$strata, levels=c('WJ','PP','GF'))
tree.C <- aggregate(cbind(agb.gC.2014,agb.gC.2004)~tree.id+subplot+plot+strata, subset(component.C, component != "FB"), sum, na.action=na.pass)
stems$tree.kgC.2014 <- tree.C$agb.gC.2014/1000
stems$tree.kgC.2004 <- tree.C$agb.gC.2004/1000

#leaf.C <- aggregate(cbind(agb.gC.2014,agb.gC.2004)~tree.id+subplot+plot+strata, subset(component.C, component == 'FB'), sum, na.action=na.pass)
# stems$leaf.kgC.2014 <- rep(NA, nrow(stems))
# stems$leaf.kgC.2004 <- rep(NA, nrow(stems))
# stems$leaf.kgC.2014[match(leaf.C$tree.id, stems$tree.id)] <- leaf.C$agb.gC.2014/1000
# stems$leaf.kgC.2004[match(leaf.C$tree.id, stems$tree.id)] <- leaf.C$agb.gC.2004/1000
# stems$leaf.kgC.2004.2014.avg <- rowMeans(cbind(stems$leaf.kgC.2014, stems$leaf.kgC.2004))

#-------------------------------------------------------------------------
# Calculate 10-yr mean NPP for both wood (add foliage NPP later)
#-------------------------------------------------------------------------
stems$tree.10yr.avg.wood.npp.kgCyr <- (stems$tree.kgC.2014-stems$tree.kgC.2004)/10

#stems$tree.10yr.avg.anpp.kgCyr <- with(stems, tree.10yr.avg.foliage.npp.kgCyr+tree.10yr.avg.wood.npp.kgCyr)
# pdf('figures/tree_npp_by_dbh.pdf', 9,6)
# xyplot(tree.10yr.avg.anpp.kgCyr~dbh.cm|species, stems, xlab='tree dbh (cm)', ylab='10-year mean tree anpp (gC/yr)')
# dev.off()

#-------------------------------------------------------------------------
# Model tree sapwood area ~ dhb using linear regression
#-------------------------------------------------------------------------
xyplot(sapwood.cm2~dbh.cm|species, sap, type=c('p','r'))
sapwood.lms <- as.data.frame(matrix(NA, ncol=8,nrow=length(species)))
colnames(sapwood.lms) <- c('species','slope','slope_se','intercept','p','r2','df','rmse')

cnt=1
i='juoc'
for (i in species){
  sp.dat <- subset(sap, species == i)
  my.lm <- lm(sapwood.cm2~0+dbh.cm, sp.dat)
  my.lm.smry <- summary(my.lm)
  resids <- my.lm.smry$residuals
  rmse  <- sqrt(mean(resids^2))
  
  sapwood.lms$species[cnt] <- i
  sapwood.lms$slope[cnt] <- round(my.lm.smry$coefficients[1],3)                                                       
  sapwood.lms$slope_se[cnt] <- round(my.lm.smry$coefficients[2],3)                                                       
  #sapwood.lms$intercept[cnt] <- round(my.lm.smry$coefficients[1],1)                                                       
  #sapwood.lms$p[cnt] <-  round(my.lm.smry[[4]][[8]],3)                                                     
  sapwood.lms$p[cnt] <-  round(my.lm.smry[[4]][[4]],3)
  sapwood.lms$r2[cnt] <- round(my.lm.smry$r.squared,2)                                   
  sapwood.lms$df[cnt] <- paste(my.lm.smry$df[1],my.lm.smry$df[2],sep=',')
  sapwood.lms$rmse[cnt] <- rmse
  assign(paste('sapwood',i,'lm', sep='.'), value=my.lm)
  cnt=cnt+1
}

# apply sapwood model to each tree
for (i in 1:n.trees){
  i.tree <- stems[i,]
  i.sp <- i.tree$species
  i.lm <- get(paste("sapwood",i.sp,'lm', sep='.'))
  stems$sapwood.cm2[i] <- predict(i.lm, newdata=data.frame(dbh.cm=i.tree$dbh.cm))
}

# output sapwood by dbh regression stats
write.table(sapwood.lms,'outputs/species_sapwood_by_dbh_lm.csv', sep=',', col.names=T, row.names=F)

#-------------------------------------------------------------------------
# output results
#-------------------------------------------------------------------------
write.csv(stems, 'data/stem_survey/tree_C_npp_sap.csv')
#-------------------------------------------------------------------------

#-------------------------------------------------------------------------
# end script
#-------------------------------------------------------------------------
