###################################################################################################################
# This R script computes stem wood density and sapwood area for indidivual trees from laboratory measurements for multiple species.
# The script must be run before estimating agb, npp, and sapwood area for all trees in the study.
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)

setwd("")
sp <- c('juoc','pipo','abgr')
sg <- read.csv('stem_traits/stem_sg_sapwood_data.csv', header=T)
stems <- read.csv('stem_survey/stem_survey_data.csv', header=T, fill=T)
phys <- read.csv('stem_traits/stem_physiology_data.csv', header=T)

strata <- unique(sg$strata)
plot <- unique(sg$plot)
subplot <- 1:5

# convert to factors
stems$strata <- factor(stems$strata)
stems$plot <- factor(stems$plot)
stems$subplot <- factor(stems$subplot)

# output data frame for tree stem specific gravity
tree.sg.df <- as.data.frame(matrix(NA, ncol=9, nrow=0))
colnames(tree.sg.df) <- c('strata','plot','subplot','tree.id','id','species','pithed','sg.g.cm3')
  

bark.mmnts <- phys[which(is.na(phys$bark.mm)==F),]
bark.mmnts$dbh.cm <- stems$dbh.cm[stems$id%in%bark.mmnts$id]
bark.mmnts$ht.total.m <- stems$ht.total.m[stems$id%in%bark.mmnts$id]
#write.table(bark.mmnts, "C:/Users/bernerl/Google Drive/Fieldwork Summer 2015/data/Berner_bark_measurements.csv", sep = ',', row.names = F, col.names = T)

xyplot(bark.mm~dbh.cm|species, bark.mmnts, type=c('p','r'))

#---------------------------------------------------------------------------------------
# calculate area-weighted stem specific gravity for each tree
for (i in strata){
  for (j in plot){
    for (k in subplot){
      trees.subplot <- subset(sg, strata == i & plot == j & subplot == k)
      for (l in unique(trees.subplot$tree.id)){
        tree <- subset(trees.subplot, tree.id==l)
        if (tree$pithed[1] == 'Yes'){
          area.total <- pi*sum(tree$length.mm)^2 # total tree basal area (w/o bark)
          area.rings <- pi*cumsum(rev(tree$length.mm))^2 # total area within each segment ; reversed because mmnts made from bark to pith
          area.segs <- area.rings - c(0,area.rings[-length(area.rings)]) # area of each segment
          tree.sg <- sum(tree$sg.g.cm3*(area.segs/area.total))
        
          } else {
            tree.phys <- subset(phys, strata == i & plot == j & subplot == k & tree.id==l)
            tree.stem <- subset(stems, strata == i & plot == j & subplot == k & tree.id==l)
            tree.r.mm <- tree.stem$dbh.cm*10/2-tree.phys$bark.mm # stem radius, excluding bark
            area.total <- pi*tree.r.mm^2
            area.rings <- pi*(tree.r.mm-cumsum(tree$length.mm))^2
            area.segs <- c(area.total, area.rings)-c(area.rings,0)
            #tree.sg.avg <- mean(tree$sg.g.cm3)
            tree.sg.avg <- tree$sg.g.cm3[length(tree$sg.g.cm3)]
            tree.sg <- sum(tree$sg.g.cm3*(area.segs[1:length(tree$sg.g.cm3)]/area.total)) + tree.sg.avg*area.segs[length(tree$sg.g.cm3)+1]/area.total
        }
        # store tree sg to data frame
        tree.sg.df <- rbind(tree.sg.df, data.frame(strata=i,plot=j,subplot=k,tree.id=l,id=tree$id[1], species=unique(tree$species),
                                                   pithed=unique(tree$pithed),sg.g.cm3=tree.sg))
      }
    }
  }
}


# compare sg of pithed and unpithed trees
bwplot(sg.g.cm3~pithed|species, tree.sg.df)
abgr.sg <- subset(tree.sg.df, species=='abgr')
pipo.sg <- subset(tree.sg.df, species=='pipo')
juoc.sg <- subset(tree.sg.df, species=='juoc')
t.test(sg.g.cm3~pithed,abgr.sg)
t.test(sg.g.cm3~pithed,pipo.sg)
t.test(sg.g.cm3~pithed,juoc.sg)
#tree.sg.df <- rbind(tree.sg.df, subset(abgr.sg, pithed=='No'))

# use only trees that were pithed
#tree.sg.df <- subset(tree.sg.df, pithed == "Yes")
tree.sg.df <- tree.sg.df[-which(tree.sg.df$species=='juoc' & tree.sg.df$pithed=='No'),] # exclude juniper that were not pithed
tree.sg.df <- na.omit(tree.sg.df)
nrow(tree.sg.df)
table(tree.sg.df$pithed)

bwplot(sg.g.cm3~species, tree.sg.df, notch=T, ylim=c(0.3,.6))
#tree.sg.df.no.pith <- subset(tree.sg.df, pithed == "No")

write.csv(tree.sg.df, "stem_traits/stem_sg_procesed.csv")
#write.csv(tree.sg.df.no.pith, 'stem_traits/stem_sg_no_pith.csv')

#---------------------------------------------------------------------------------------
# calculate species-specific SG summaries
#---------------------------------------------------------------------------------------
sp.sg <- aggregate(sg.g.cm3~species, tree.sg.df, function(x){round(mean(x),3)})
sp.sg$sd <- aggregate(sg.g.cm3~species, tree.sg.df, function(x){round(sd(x),3)})[,2]
sp.sg$n <- aggregate(sg.g.cm3~species, tree.sg.df, function(x){round(length(x),3)})[,2]
sp.sg$se <- round(sp.sg$sd/sqrt(sp.sg$n),3)
write.table(sp.sg, 'species_summaries/sg_by_species.csv', sep=',', row.names=F, col.names=T)

#---------------------------------------------------------------------------------------
# calculate sapwood area
#---------------------------------------------------------------------------------------
sapwood <- subset(sg, segment == "SW")
sapwood$strata <- factor(sapwood$strata, levels=c('WJ','PP','GF'))
sapwood.tree <- aggregate(length.mm~species+tree.id+subplot+plot+strata, sapwood, FUN=sum)
sapwood.tree$id <- unique(sapwood$id)
sapwood.tree$dbh.cm <- stems$dbh.cm[stems$id%in%sapwood.tree$id]
sapwood.tree$sapwood.2.ba <- rep(NA, nrow(sapwood.tree))

# model bark thickness ~ dbh for each species
bark.mmnts <- phys[which(is.na(phys$bark.mm)==F),]
bark.mmnts$dbh.cm <- stems$dbh.cm[stems$id%in%bark.mmnts$id]
bark.mmnts$ht.total.m <- stems$ht.total.m[stems$id%in%bark.mmnts$id]

for (i in sapwood.tree$id){
  i.sap <- sapwood.tree[sapwood.tree$id == i,]
  i.stem <- stems[stems$id == i,]
  i.bark <- bark.mmnts[bark.mmnts$id ==i,]
  i.sp <- i.stem$species
  
  if (length(i.bark$bark.mm)==0){ # if bark measurements are lacking, then model bark diameter 
    bark.sp <- subset(bark.mmnts, species == i.sp)
    bark.lm <- lm(bark.mm~dbh.cm, bark.sp)
    bark.cm <- predict(bark.lm, newdata=data.frame(dbh.cm=i.stem$dbh.cm))/10
  } else {bark.cm <- i.bark$bark.mm/10}
  
  dbh.wobark.cm <- i.stem$dbh.cm-bark.cm*2 # tree dbh after removing bark
  ba.wobark.cm2 <- pi*(dbh.wobark.cm/2)^2 # tree basal area after removing bark
  sap.cm2 <- ba.wobark.cm2- (pi*((dbh.wobark.cm-(i.sap$length.mm/10)*2)/2)^2) # tree sapwood area
  sap.2.ba <- sap.cm2/ba.wobark.cm2
  sapwood.tree$sapwood.cm2[sapwood.tree$id == i] <- sap.cm2
  sapwood.tree$ba.cm2[sapwood.tree$id == i] <- ba.wobark.cm2
  sapwood.tree$sapwood.2.ba[sapwood.tree$id == i] <- sap.2.ba
}

#bwplot(sapwood.2.ba~species, sapwood.tree)
write.csv(sapwood.tree, 'stem_traits/stem_sapwood_processed.csv')

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