###################################################################################################################
# This R script corrects LAI-2200 measurements for interception by wood and for shoot clumping.
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)
setwd('')

# equation for wood area index pulled from C:\research\data\Master Field data 1999_2006\Master data\Master LAI.xls
WAI <- function(BA){2.061*(1-exp(-0.0062 * BA))} # where BASAL AREA is in m^2ha-1.

lai <- read.csv('lai2200/lai_data.csv', header=T)
treeC <- read.csv('stem_survey/tree_C_npp_sap.csv', header=T)
plot.dims <- read.csv('plot_details/plot_dimensions.csv', header=T)
treeC.subplot.sp <- read.csv('subplot_summaries/treeC_by_subplot_sp.csv', header=T)

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

# factorize
lai$strata <- factor(lai$strata, levels=c('WJ','PP','GF'))
lai$plot <- factor(lai$plot)
lai$subplot <- factor(lai$subplot)
lai$subplot.id <- factor(lai$subplot.id)

treeC$strata <- factor(treeC$strata, levels=c('WJ','PP','GF'))

#--------------------------------------------------------------------------
# calculate wood area index for each subplot based on basal area
#--------------------------------------------------------------------------
treeC$ba.cm2 <- 3.14159*(treeC$dbh.cm/2)^2
stem.ba <- aggregate(ba.cm2~size.class+subplot+plot+strata, treeC, sum)
stem.ba$ba.m2 <- stem.ba$ba.cm2/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)
lai$wai <- round(WAI(subplot.ba.m2ha$ba.m2ha),2)



#--------------------------------------------------------------------------
# Calculate needle to shoot clumping factor for each subplot (gamma correction), 
#   - weighing corrections for each subplot by the relative basal area of each species.   
#   - Pipo used also for Pico (Law et al. 2001, Tree Phys)
#   - Abam correction used for Abgr and Abpr (Frazer et al. 2000, CJFR)
#   - Thpl correction (cross-site mean = 1.01 +- 0.034) used for Cade (Frazer et al. 2000, CJFR)
#--------------------------------------------------------------------------
gamma.factors <- data.frame(species=c('abgr','abpr','cade','juoc','pico','pipo'), gamma=c(2.35,2.35,1.01,1,1.29,1.29)) 
lai$gamma <- rep(NA, nrow(lai))

for (i in unique(lai$subplot.id)){
  i.treeC.subplot.sp <- subset(treeC.subplot.sp, subplot.id==i)
  i.gamma <- gamma.factors[gamma.factors$species%in%i.treeC.subplot.sp$species,]
  subplot.gamma <- sum(i.gamma$gamma*i.treeC.subplot.sp$sp.frac)
  lai$gamma[lai$subplot.id%in%i] <- round(subplot.gamma,2)
}

#--------------------------------------------------------------------------
# Calculate needle to clumping at levels above shoots for each subplot (omega correction)
#--------------------------------------------------------------------------
#lai$omega <- rep(0.88, nrow(lai))# mean omega from 181 plots measured using TRAC as part of EPA and ORCA projects (0.88 +- 0.11 SD)
# NOTE:: The LAI-2200 automatically accounts for the apparant clumping factor (ACF), which is equivalent to the omega correction.
# Therefore, there is no need to correct for omega.

#--------------------------------------------------------------------------
# Calculate half-surface area LAI correced for clumping and light interception by wood as per Chen (1996) and Law et al. (2008)
#--------------------------------------------------------------------------
#lai$lai.hc <- with(lai, (lai*gamma/omega-wai))
lai$lai.hc <- with(lai, (lai*gamma-wai))
write.table(lai, 'lai2200/lai_processed.csv', sep=',', row.names=F, col.names=T)
write.table(lai, 'subplot_summaries/lai_by_subplot.csv', sep=',', row.names=F, col.names=T)

#--------------------------------------------------------------------------
# Calculate site mean/sd 
#--------------------------------------------------------------------------
# Measurements on several subplots were suprisingly different than other measurements made on the plots; these mmnts were omitted
# These measurements were dropped
lai <- lai[-which(lai$subplot.id%in%c('GF.1.1','GF.2.4','GF.5.4')),]

lai.plot.stats <-aggregate(lai.hc~plot+strata, lai, mean)
lai.plot.stats$lai.hc <- round(lai.plot.stats$lai.hc, 2)
colnames(lai.plot.stats)[3] <- "lai.hc.avg"
lai.plot.stats$lai.hc.sd <- round(aggregate(lai.hc~plot+strata, lai, sd)[,3],2)
lai.plot.stats$n <- aggregate(lai.hc~plot+strata, lai, function(x,...){length(x,...)})[,3]

write.table(lai.plot.stats, 'plot_summaries/lai_by_plot.csv', sep=',', row.names=F, col.names=T)

#--------------------------------------------------------------------------
# End script
#--------------------------------------------------------------------------