R/SummaryAGWP2.r
In GabrielaLopez/BiomasaFP: Tools for analysing data downloaded from ForestPlots.net
SummaryAGWP2<-function (xdataset, AGBEquation, dbh ="D4",rec.meth=0,height.data=NULL,param.type="Best",palm.eq=TRUE){
#Use AGBEquation to set Chv14 argument
if(all.equal(AGBEquation,AGBChv14)==TRUE){Chv14=TRUE}else{Chv14=FALSE}
AGBData <- AGBEquation (xdataset, dbh,height.data=height.data,param.type=param.type)
recs<-unobs.recs(AGBData)
test<-recs[[1]]
cns<-c()
for(i in 1:ncol(test)){
cns<-c(cns,rep(i,nrow(test)))
}
plts<-rep(dimnames(test)[[1]],ncol(test))
unobs.recs<-as.vector(test)
LS.unobs.recs<-as.vector(recs[[2]])
LS.dead<-as.vector(recs[[3]])
unobs.dat<-data.frame("PlotViewID"=plts,"Census.No"=cns+1,unobs.recs,LS.unobs.recs,LS.dead,stringsAsFactors=F)
if(palm.eq==TRUE){
AGBData[AGBData$Monocot==1,]$AGBind<-GoodmanPalm(AGBData[AGBData$Monocot==1,],dbh=dbh)
AGBData[AGBData$Monocot==1,]$AGBDead<-GoodmanPalm(AGBData[AGBData$Monocot==1,],dbh=paste(dbh,"_D",sep=""))
}
IndAL <- aggregate (Alive/PlotArea ~ PlotViewID + Census.No, data = AGBData, FUN=sum )
AGBAlive <-aggregate (AGBind/PlotArea ~ PlotViewID + Census.No, data = AGBData, FUN=sum )
#AGBData$Census.prev<-AGBData[match(paste(AGBData$TreeID,AGBData$Census.No-1),paste(AGBData$TreeID,AGBData$Census.No)),"Census.Mean.Date"]
#AGBData$Delta.time<-AGBData$Census.Mean.Date-AGBData$Census.prev
AGBData$Census.prev<-AGBData[match(paste(AGBData$PlotViewID,AGBData$Census.No-1),paste(AGBData$PlotViewID,AGBData$Census.No)),"Census.Mean.Date"]
AGBData$Delta.time<-AGBData$Census.Mean.Date-AGBData$Census.prev
#Estimate AGB of recruits at 100mm
Height100<-height.mod(100,AGBData$a_par,AGBData$b_par,AGBData$c_par)
if(Chv14==F){
AGBData$AGB.100<-(0.0509*AGBData$WD * ((100/10)^2)* Height100)/1000
}else{
AGBData$AGB.100<-(0.0673*(AGBData$WD * ((100/10)^2)* Height100)^0.976)/1000
}
if(palm.eq==TRUE){
if(length(AGBData[AGBData$Monocot==1,]$AGB.100)>0){
AGBData[AGBData$Monocot==1,]$AGB.100<-exp(-3.3488+(2.7483*log(10)))/1000
}
}
AGBData$AGWPRec<-AGBData$AGBind-AGBData$AGB.100
#Calculate growth of surviving trees
#AGBData$AGB.Prev<-AGBData[match(paste(AGBData$TreeID,AGBData$Census.No-1),paste(AGBData$TreeID,AGBData$Census.No)),"AGBind"]
AGBData$AGB.Prev<-AGBData[match(paste(AGBData$TreeID,AGBData$PlotViewID,AGBData$Census.No-1),paste(AGBData$TreeID,AGBData$PlotViewID,AGBData$Census.No)),"AGBind"]
AGBData$Delta.AGB<-AGBData$AGBind-AGBData$AGB.Prev
AGBData$Survive2<-ifelse(AGBData$Census.No>1 & AGBData$Recruit==0 & AGBData$Alive==1 & AGBData$Snapped==0,1,0)
#Includes snapped trees - for stem dynamics
AGBData$Survive<-ifelse(AGBData$Census.No>1 & AGBData$Recruit==0 & AGBData$Alive==1 ,1,0)
AGB.surv<-aggregate(Delta.AGB/PlotArea ~ PlotViewID + Census.No, data = AGBData[AGBData$Survive2==1,], FUN=sum )
IndSurv<-aggregate(Survive/PlotArea ~ PlotViewID + Census.No, data = AGBData, FUN=sum )
#Calculate mean WD per census
WD<-aggregate(WD ~ PlotViewID + Census.No, data = AGBData[AGBData$Alive==1,], FUN=function(x)mean(x,na.rm=T))
# Find Recruits
Recruits<- AGBData[AGBData$Recruit==1,]
if(rec.meth==100){
AGBRec <- aggregate (AGWPRec/PlotArea ~ PlotViewID + Census.No, data = Recruits, FUN=sum ) # NOTE: divides by PlotArea, so you get value per ha !!
}else{
AGBRec <- aggregate (AGBind/PlotArea ~ PlotViewID + Census.No, data = Recruits, FUN=sum )
}
colnames(AGBRec) <- c('PlotViewID','Census.No','AGBRec.PlotArea')
IndRec <- aggregate (Recruit/PlotArea ~ PlotViewID + Census.No, data = Recruits, FUN=sum )
#merge recruit information
Recs<- merge(AGBRec, IndRec, by = c('PlotViewID','Census.No'), all.x=TRUE)
# Dead stems get only stems that are dead and have an AGB_D, to count only dead trees ones
DeadTrees <-AGBData[AGBData$Dead==1 & !is.na(AGBData$D1_D),]
#Match in size class growth rate
IndDead <- aggregate (Dead/PlotArea ~ PlotViewID + Census.No, data = DeadTrees, FUN=sum )
AGBDe <-aggregate (AGBDead/PlotArea ~ PlotViewID + Census.No, data = AGBData, FUN=sum )
#merge Dead trees
Deads <- merge(AGBDe, IndDead, by = c('PlotViewID','Census.No'),all.x=TRUE)
# Unobserved growth of dead trees
growth.rate<-SizeClassGrowth(xdataset,dbh=dbh)
dead2<-merge(DeadTrees,growth.rate,by="PlotViewID",all.x=T)
dead2$LS.dead<-unobs.dat[match(paste(dead2$PlotViewID,dead2$Census.No),paste(unobs.dat$PlotViewID,unobs.dat$Census.No)),"LS.dead"]
dead2$DBH.death<-ifelse(dead2[,paste(dbh,"_D",sep="")]<200,dead2$LS.dead*dead2$Class1,
ifelse(dead2[,paste(dbh,"_D",sep="")]<400,dead2$LS.dead*dead2$Class2,
dead2$LS.dead*dead2$Class3))
dead2$DBH.death<-dead2$DBH.death+dead2[,paste(dbh,"_D",sep="")]
#Height at death
dead2$Height.dead<-height.mod(dead2$DBH.death,dead2$a_par,dead2$b_par,dead2$c_par)
#AGB at death
if(Chv14==F){
dead2$AGB.death2<-(0.0509*dead2$WD * ((dead2$DBH.death/10)^2)* dead2$Height.dead)/1000
}else{
dead2$AGB.death2<-(0.0673*(dead2$WD * ((dead2$DBH.death/10)^2)* dead2$Height.dead)^0.976)/1000
}
dead2<-dead2[!is.na(dead2$Monocot),]
if(palm.eq==TRUE){
dead2$AGB.death2[dead2$Monocot==1]<-GoodmanPalm(dead2[dead2$Monocot==1,],paste(dbh,"_D",sep=""))
}
dead2$Unobs.dead<-dead2$AGB.death2-dead2$AGBDead
unobs.dead<-aggregate(Unobs.dead/PlotArea~PlotViewID + Census.No, data = dead2, FUN=sum )
# merge all summaries
mergeAGBAlive <- merge (AGBAlive, IndAL, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeB <- merge(mergeAGBAlive, Recs, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeC <- merge (mergeB, Deads, by= c('PlotViewID','Census.No'),all.x=TRUE)
mergeD<-merge (mergeC,AGB.surv, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeE<-merge (mergeD,IndSurv, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeF<-merge(mergeE,unobs.dead, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeF<-merge(mergeF,WD,by=c('PlotViewID','Census.No'),all.x=TRUE)
#Correct names of merge F to be legal
names(mergeF)<-sub("/",".",names(mergeF)) # NOTE: .PlotArea replaces /PlotArea, so actually values per ha (plotarea weighted)
# the block of code hereafter prepares tmp, which will be merged to mergeF in the end
# tmp adds only gains from unobserved recruits, all other biomass variables are calculated above
# problem : the following line tries to identify a unique set of abc parameters per plotviewid
# solution : remove TreeID-specific optimal parameters (a_Best_tree,...) from AGBData and add plotview-specific optimal parameters (a_Best,...) at this point
if (all.equal(param.type, "Best_tree") == TRUE){ # CHANGED: new line added
AGBData$a_par <- NULL ; AGBData$b_par <- NULL ; AGBData$c_par <- NULL # CHANGED: new line added
parameters <- unique(height.data[, c("PlotViewID", "Census.No", "a_Best", "b_Best", "c_Best")]) # CHANGED: new line added
names(parameters) <- c("PlotViewID", "Census.No", "a_par", "b_par", "c_par") # CHANGED: new line added
AGBData <- merge(AGBData,parameters,all.x=TRUE) } # CHANGED: new line added
#Estimate growth of unobserved recruits - 1st merge growth rate data with allometric parameter data
unobs.dat2<-unique(AGBData[!is.na(AGBData$Delta.time),c("PlotViewID","Census.No","a_par","b_par","c_par","Delta.time")])
unobs.dat2<-merge(unobs.dat2,growth.rate,by="PlotViewID",all.x=T)
#Then merge with existing data to get number of recruits and deaths
tmp<-merge(mergeF,unobs.dat2,by = c('PlotViewID','Census.No'),all.x=TRUE)
#Get number of stems in previous census
tmp$PrevStems<-tmp[match(paste(tmp$PlotViewID,tmp$Census.No-1),paste(tmp$PlotViewID,tmp$Census.No)),"Alive.PlotArea"]
#Estimate recruitment and mortality rates
stems<-tmp$PrevStems
death<-tmp$Dead.PlotArea
recruits<-tmp$Recruit.PlotArea
interval<-tmp$Delta.time
tmp$rec.rate<-recruits/stems/interval
tmp$mort.rate<-death/stems/interval
#tmp$rec.rate <- replace(tmp$rec.rate,is.na(tmp$rec.rate),0) # CHANGE: added, otherwise function doesn't work as crashes on plotviews with 0 recruitment
#tmp$mort.rate <- replace(tmp$mort.rate,is.na(tmp$mort.rate),0)
#Get time weighted mean of rec and mort rates for each plot
a<-split(tmp,f=tmp$PlotViewID)
rec1<-unlist(lapply(a,function(x)weighted.mean(x$rec.rate,x$Delta.time,na.rm=T)))
mort1<-unlist(lapply(a,function(x)weighted.mean(x$mort.rate,x$Delta.time,na.rm=T)))
rec.means<-data.frame("PlotViewID"=names(rec1),"Rec.mean"=as.numeric(rec1),stringsAsFactors=F)
mort.means<-data.frame("PlotViewID"=names(mort1),"Mort.mean"=as.numeric(mort1),stringsAsFactors=F)
plot.means<-merge(rec.means,mort.means,by="PlotViewID",all.x=T)
tmp<-merge(tmp,plot.means,by="PlotViewID",all.x=T)
tmp$LS.rec<-unobs.dat[match(paste(tmp$PlotViewID,tmp$Census.No),paste(unobs.dat$PlotViewID,unobs.dat$Census.No)),"LS.unobs.recs"]
#Growth of unobserved recruits
unobs.D<-tmp$Class1*tmp$LS.rec
#D at death
D.death<-100+unobs.D
#Get height at death
height.death<-height.mod(D.death,tmp$a_par,tmp$b_par,tmp$c_par)
#Height at point of recruitment
height.rec<-height.mod(100,tmp$a_par,tmp$b_par,tmp$c_par)
#WD_mean<-tmp$MeanWD # Commented out - mean across all censuses
WD_mean<-tmp$WD # Use mean in plot and census
#Estimate AGB at death
if(Chv14==F){
AGB.death<-(0.0509*WD_mean * ((D.death/10)^2)* height.death)/1000
}else{
AGB.death<-(0.0673*(WD_mean * ((D.death/10)^2)* height.death)^0.976)/1000
}
#Estimate AGB at point of recruitment
if(Chv14==F){
AGB.rec<-(0.0509*WD_mean * ((100/10)^2)* height.rec)/1000
}else{
AGB.rec<-(0.0673*(WD_mean * ((100/10)^2)* height.rec)^0.976)/1000
}
#Get number of unobserved recruits
tmp$unobs.rec<-unobs.dat[match(paste(tmp$PlotViewID,tmp$Census.No),paste(unobs.dat$PlotViewID,unobs.dat$Census.No)),"unobs.recs"]
#Get growth of unobs recruits
#Option for recruits starting at 100
if(rec.meth==100){
unobs.growth<-(tmp$unobs.rec*(AGB.death-AGB.rec))
}else{
unobs.growth<-(tmp$unobs.rec*AGB.death)
}
tmp$UnobsGrowth<-unobs.growth
tmp$UnobsRec<-tmp$unobs.rec
#Simplyfy tmp and merge with MergeF
tmp<-tmp[,c("PlotViewID","Census.No","rec.rate","mort.rate","UnobsGrowth","UnobsRec","Delta.time")]
mergeG<-merge (mergeF,tmp, by = c('PlotViewID','Census.No'),all.x=TRUE)
#Create output file
SummaryB <- mergeG
SummaryB <- SummaryB[order(SummaryB$PlotViewID, SummaryB$Census.No, decreasing=FALSE), ]
#Add AGWP per ha
SummaryB$AGWP.PlotArea<-with(SummaryB,ifelse(is.na(Delta.AGB.PlotArea),0,Delta.AGB.PlotArea)+ifelse(is.na(AGBRec.PlotArea),0,AGBRec.PlotArea)+ifelse(is.na(Unobs.dead.PlotArea),0,Unobs.dead.PlotArea)+ifelse(is.na(UnobsGrowth),0,UnobsGrowth)) # NOTE: .PlotArea means plotarea weighted so per ha
SummaryB$AGWP.PlotArea.Year<-SummaryB$AGWP.PlotArea/SummaryB$Delta.time
#Neaten up column names
names(SummaryB)<-c("PlotViewID","Census.No","AGB.ha","Stems.ha",
"AGWPrec.ha","Recruit.ha","AGBmort.ha","Mortality.ha","AGWPsurv.ha","SurvivingStems.ha","UnobsAGWPmort.ha","Mean.WD", # NOTE: from mergeF
"Recruitment.stem.year","Mortality.stem.year","UnobsAGWPrec.ha","UnobsRecruits.ha","CensusInterval", # NOTE: from tmp, all per ha as PlotArea weighted !
"AGWP.ha","AGWP.ha.year")
#Extract recruitment of monocots in each census
if(sum(Recruits$Monocot)>0 &!is.na(sum(Recruits$Monocot))){
if(rec.meth==100){
palm.rec<-aggregate(AGWPRec/PlotArea ~ PlotViewID+Census.No,data=Recruits[Recruits$Monocot==1,],FUN=sum)
}else{
palm.rec<-aggregate(AGBind/PlotArea ~ PlotViewID+Census.No,data=Recruits[Recruits$Monocot==1,],FUN=sum)
}
names(palm.rec)[3]<-"MonocotRecruits.ha"
SummaryB<-merge(SummaryB,palm.rec,by = c('PlotViewID','Census.No'),all.x=T)
}
if(sum(DeadTrees$Monocot)>0 &!is.na(sum(DeadTrees$Monocot))){
#Mortality of monocots in each census
palm.dead <-aggregate (AGBDead/PlotArea ~ PlotViewID + Census.No, data = AGBData[AGBData$Monocot==1,], FUN=sum )
names(palm.dead)[3]<-"MonocotMortality.ha"
SummaryB<-merge(SummaryB,palm.dead,by = c('PlotViewID','Census.No'),all.x=T)
}
SummaryB
}
GabrielaLopez/BiomasaFP documentation built on June 23, 2020, 9:12 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
SummaryAGWP2<-function (xdataset, AGBEquation, dbh ="D4",rec.meth=0,height.data=NULL,param.type="Best",palm.eq=TRUE){
#Use AGBEquation to set Chv14 argument
if(all.equal(AGBEquation,AGBChv14)==TRUE){Chv14=TRUE}else{Chv14=FALSE}
AGBData <- AGBEquation (xdataset, dbh,height.data=height.data,param.type=param.type)
recs<-unobs.recs(AGBData)
test<-recs[[1]]
cns<-c()
for(i in 1:ncol(test)){
cns<-c(cns,rep(i,nrow(test)))
}
plts<-rep(dimnames(test)[[1]],ncol(test))
unobs.recs<-as.vector(test)
LS.unobs.recs<-as.vector(recs[[2]])
LS.dead<-as.vector(recs[[3]])
unobs.dat<-data.frame("PlotViewID"=plts,"Census.No"=cns+1,unobs.recs,LS.unobs.recs,LS.dead,stringsAsFactors=F)
if(palm.eq==TRUE){
AGBData[AGBData$Monocot==1,]$AGBind<-GoodmanPalm(AGBData[AGBData$Monocot==1,],dbh=dbh)
AGBData[AGBData$Monocot==1,]$AGBDead<-GoodmanPalm(AGBData[AGBData$Monocot==1,],dbh=paste(dbh,"_D",sep=""))
}
IndAL <- aggregate (Alive/PlotArea ~ PlotViewID + Census.No, data = AGBData, FUN=sum )
AGBAlive <-aggregate (AGBind/PlotArea ~ PlotViewID + Census.No, data = AGBData, FUN=sum )
#AGBData$Census.prev<-AGBData[match(paste(AGBData$TreeID,AGBData$Census.No-1),paste(AGBData$TreeID,AGBData$Census.No)),"Census.Mean.Date"]
#AGBData$Delta.time<-AGBData$Census.Mean.Date-AGBData$Census.prev
AGBData$Census.prev<-AGBData[match(paste(AGBData$PlotViewID,AGBData$Census.No-1),paste(AGBData$PlotViewID,AGBData$Census.No)),"Census.Mean.Date"]
AGBData$Delta.time<-AGBData$Census.Mean.Date-AGBData$Census.prev
#Estimate AGB of recruits at 100mm
Height100<-height.mod(100,AGBData$a_par,AGBData$b_par,AGBData$c_par)
if(Chv14==F){
AGBData$AGB.100<-(0.0509*AGBData$WD * ((100/10)^2)* Height100)/1000
}else{
AGBData$AGB.100<-(0.0673*(AGBData$WD * ((100/10)^2)* Height100)^0.976)/1000
}
if(palm.eq==TRUE){
if(length(AGBData[AGBData$Monocot==1,]$AGB.100)>0){
AGBData[AGBData$Monocot==1,]$AGB.100<-exp(-3.3488+(2.7483*log(10)))/1000
}
}
AGBData$AGWPRec<-AGBData$AGBind-AGBData$AGB.100
#Calculate growth of surviving trees
#AGBData$AGB.Prev<-AGBData[match(paste(AGBData$TreeID,AGBData$Census.No-1),paste(AGBData$TreeID,AGBData$Census.No)),"AGBind"]
AGBData$AGB.Prev<-AGBData[match(paste(AGBData$TreeID,AGBData$PlotViewID,AGBData$Census.No-1),paste(AGBData$TreeID,AGBData$PlotViewID,AGBData$Census.No)),"AGBind"]
AGBData$Delta.AGB<-AGBData$AGBind-AGBData$AGB.Prev
AGBData$Survive2<-ifelse(AGBData$Census.No>1 & AGBData$Recruit==0 & AGBData$Alive==1 & AGBData$Snapped==0,1,0)
#Includes snapped trees - for stem dynamics
AGBData$Survive<-ifelse(AGBData$Census.No>1 & AGBData$Recruit==0 & AGBData$Alive==1 ,1,0)
AGB.surv<-aggregate(Delta.AGB/PlotArea ~ PlotViewID + Census.No, data = AGBData[AGBData$Survive2==1,], FUN=sum )
IndSurv<-aggregate(Survive/PlotArea ~ PlotViewID + Census.No, data = AGBData, FUN=sum )
#Calculate mean WD per census
WD<-aggregate(WD ~ PlotViewID + Census.No, data = AGBData[AGBData$Alive==1,], FUN=function(x)mean(x,na.rm=T))
# Find Recruits
Recruits<- AGBData[AGBData$Recruit==1,]
if(rec.meth==100){
AGBRec <- aggregate (AGWPRec/PlotArea ~ PlotViewID + Census.No, data = Recruits, FUN=sum ) # NOTE: divides by PlotArea, so you get value per ha !!
}else{
AGBRec <- aggregate (AGBind/PlotArea ~ PlotViewID + Census.No, data = Recruits, FUN=sum )
}
colnames(AGBRec) <- c('PlotViewID','Census.No','AGBRec.PlotArea')
IndRec <- aggregate (Recruit/PlotArea ~ PlotViewID + Census.No, data = Recruits, FUN=sum )
#merge recruit information
Recs<- merge(AGBRec, IndRec, by = c('PlotViewID','Census.No'), all.x=TRUE)
# Dead stems get only stems that are dead and have an AGB_D, to count only dead trees ones
DeadTrees <-AGBData[AGBData$Dead==1 & !is.na(AGBData$D1_D),]
#Match in size class growth rate
IndDead <- aggregate (Dead/PlotArea ~ PlotViewID + Census.No, data = DeadTrees, FUN=sum )
AGBDe <-aggregate (AGBDead/PlotArea ~ PlotViewID + Census.No, data = AGBData, FUN=sum )
#merge Dead trees
Deads <- merge(AGBDe, IndDead, by = c('PlotViewID','Census.No'),all.x=TRUE)
# Unobserved growth of dead trees
growth.rate<-SizeClassGrowth(xdataset,dbh=dbh)
dead2<-merge(DeadTrees,growth.rate,by="PlotViewID",all.x=T)
dead2$LS.dead<-unobs.dat[match(paste(dead2$PlotViewID,dead2$Census.No),paste(unobs.dat$PlotViewID,unobs.dat$Census.No)),"LS.dead"]
dead2$DBH.death<-ifelse(dead2[,paste(dbh,"_D",sep="")]<200,dead2$LS.dead*dead2$Class1,
ifelse(dead2[,paste(dbh,"_D",sep="")]<400,dead2$LS.dead*dead2$Class2,
dead2$LS.dead*dead2$Class3))
dead2$DBH.death<-dead2$DBH.death+dead2[,paste(dbh,"_D",sep="")]
#Height at death
dead2$Height.dead<-height.mod(dead2$DBH.death,dead2$a_par,dead2$b_par,dead2$c_par)
#AGB at death
if(Chv14==F){
dead2$AGB.death2<-(0.0509*dead2$WD * ((dead2$DBH.death/10)^2)* dead2$Height.dead)/1000
}else{
dead2$AGB.death2<-(0.0673*(dead2$WD * ((dead2$DBH.death/10)^2)* dead2$Height.dead)^0.976)/1000
}
dead2<-dead2[!is.na(dead2$Monocot),]
if(palm.eq==TRUE){
dead2$AGB.death2[dead2$Monocot==1]<-GoodmanPalm(dead2[dead2$Monocot==1,],paste(dbh,"_D",sep=""))
}
dead2$Unobs.dead<-dead2$AGB.death2-dead2$AGBDead
unobs.dead<-aggregate(Unobs.dead/PlotArea~PlotViewID + Census.No, data = dead2, FUN=sum )
# merge all summaries
mergeAGBAlive <- merge (AGBAlive, IndAL, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeB <- merge(mergeAGBAlive, Recs, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeC <- merge (mergeB, Deads, by= c('PlotViewID','Census.No'),all.x=TRUE)
mergeD<-merge (mergeC,AGB.surv, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeE<-merge (mergeD,IndSurv, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeF<-merge(mergeE,unobs.dead, by = c('PlotViewID','Census.No'),all.x=TRUE)
mergeF<-merge(mergeF,WD,by=c('PlotViewID','Census.No'),all.x=TRUE)
#Correct names of merge F to be legal
names(mergeF)<-sub("/",".",names(mergeF)) # NOTE: .PlotArea replaces /PlotArea, so actually values per ha (plotarea weighted)
# the block of code hereafter prepares tmp, which will be merged to mergeF in the end
# tmp adds only gains from unobserved recruits, all other biomass variables are calculated above
# problem : the following line tries to identify a unique set of abc parameters per plotviewid
# solution : remove TreeID-specific optimal parameters (a_Best_tree,...) from AGBData and add plotview-specific optimal parameters (a_Best,...) at this point
if (all.equal(param.type, "Best_tree") == TRUE){ # CHANGED: new line added
AGBData$a_par <- NULL ; AGBData$b_par <- NULL ; AGBData$c_par <- NULL # CHANGED: new line added
parameters <- unique(height.data[, c("PlotViewID", "Census.No", "a_Best", "b_Best", "c_Best")]) # CHANGED: new line added
names(parameters) <- c("PlotViewID", "Census.No", "a_par", "b_par", "c_par") # CHANGED: new line added
AGBData <- merge(AGBData,parameters,all.x=TRUE) } # CHANGED: new line added
#Estimate growth of unobserved recruits - 1st merge growth rate data with allometric parameter data
unobs.dat2<-unique(AGBData[!is.na(AGBData$Delta.time),c("PlotViewID","Census.No","a_par","b_par","c_par","Delta.time")])
unobs.dat2<-merge(unobs.dat2,growth.rate,by="PlotViewID",all.x=T)
#Then merge with existing data to get number of recruits and deaths
tmp<-merge(mergeF,unobs.dat2,by = c('PlotViewID','Census.No'),all.x=TRUE)
#Get number of stems in previous census
tmp$PrevStems<-tmp[match(paste(tmp$PlotViewID,tmp$Census.No-1),paste(tmp$PlotViewID,tmp$Census.No)),"Alive.PlotArea"]
#Estimate recruitment and mortality rates
stems<-tmp$PrevStems
death<-tmp$Dead.PlotArea
recruits<-tmp$Recruit.PlotArea
interval<-tmp$Delta.time
tmp$rec.rate<-recruits/stems/interval
tmp$mort.rate<-death/stems/interval
#tmp$rec.rate <- replace(tmp$rec.rate,is.na(tmp$rec.rate),0) # CHANGE: added, otherwise function doesn't work as crashes on plotviews with 0 recruitment
#tmp$mort.rate <- replace(tmp$mort.rate,is.na(tmp$mort.rate),0)
#Get time weighted mean of rec and mort rates for each plot
a<-split(tmp,f=tmp$PlotViewID)
rec1<-unlist(lapply(a,function(x)weighted.mean(x$rec.rate,x$Delta.time,na.rm=T)))
mort1<-unlist(lapply(a,function(x)weighted.mean(x$mort.rate,x$Delta.time,na.rm=T)))
rec.means<-data.frame("PlotViewID"=names(rec1),"Rec.mean"=as.numeric(rec1),stringsAsFactors=F)
mort.means<-data.frame("PlotViewID"=names(mort1),"Mort.mean"=as.numeric(mort1),stringsAsFactors=F)
plot.means<-merge(rec.means,mort.means,by="PlotViewID",all.x=T)
tmp<-merge(tmp,plot.means,by="PlotViewID",all.x=T)
tmp$LS.rec<-unobs.dat[match(paste(tmp$PlotViewID,tmp$Census.No),paste(unobs.dat$PlotViewID,unobs.dat$Census.No)),"LS.unobs.recs"]
#Growth of unobserved recruits
unobs.D<-tmp$Class1*tmp$LS.rec
#D at death
D.death<-100+unobs.D
#Get height at death
height.death<-height.mod(D.death,tmp$a_par,tmp$b_par,tmp$c_par)
#Height at point of recruitment
height.rec<-height.mod(100,tmp$a_par,tmp$b_par,tmp$c_par)
#WD_mean<-tmp$MeanWD # Commented out - mean across all censuses
WD_mean<-tmp$WD # Use mean in plot and census
#Estimate AGB at death
if(Chv14==F){
AGB.death<-(0.0509*WD_mean * ((D.death/10)^2)* height.death)/1000
}else{
AGB.death<-(0.0673*(WD_mean * ((D.death/10)^2)* height.death)^0.976)/1000
}
#Estimate AGB at point of recruitment
if(Chv14==F){
AGB.rec<-(0.0509*WD_mean * ((100/10)^2)* height.rec)/1000
}else{
AGB.rec<-(0.0673*(WD_mean * ((100/10)^2)* height.rec)^0.976)/1000
}
#Get number of unobserved recruits
tmp$unobs.rec<-unobs.dat[match(paste(tmp$PlotViewID,tmp$Census.No),paste(unobs.dat$PlotViewID,unobs.dat$Census.No)),"unobs.recs"]
#Get growth of unobs recruits
#Option for recruits starting at 100
if(rec.meth==100){
unobs.growth<-(tmp$unobs.rec*(AGB.death-AGB.rec))
}else{
unobs.growth<-(tmp$unobs.rec*AGB.death)
}
tmp$UnobsGrowth<-unobs.growth
tmp$UnobsRec<-tmp$unobs.rec
#Simplyfy tmp and merge with MergeF
tmp<-tmp[,c("PlotViewID","Census.No","rec.rate","mort.rate","UnobsGrowth","UnobsRec","Delta.time")]
mergeG<-merge (mergeF,tmp, by = c('PlotViewID','Census.No'),all.x=TRUE)
#Create output file
SummaryB <- mergeG
SummaryB <- SummaryB[order(SummaryB$PlotViewID, SummaryB$Census.No, decreasing=FALSE), ]
#Add AGWP per ha
SummaryB$AGWP.PlotArea<-with(SummaryB,ifelse(is.na(Delta.AGB.PlotArea),0,Delta.AGB.PlotArea)+ifelse(is.na(AGBRec.PlotArea),0,AGBRec.PlotArea)+ifelse(is.na(Unobs.dead.PlotArea),0,Unobs.dead.PlotArea)+ifelse(is.na(UnobsGrowth),0,UnobsGrowth)) # NOTE: .PlotArea means plotarea weighted so per ha
SummaryB$AGWP.PlotArea.Year<-SummaryB$AGWP.PlotArea/SummaryB$Delta.time
#Neaten up column names
names(SummaryB)<-c("PlotViewID","Census.No","AGB.ha","Stems.ha",
"AGWPrec.ha","Recruit.ha","AGBmort.ha","Mortality.ha","AGWPsurv.ha","SurvivingStems.ha","UnobsAGWPmort.ha","Mean.WD", # NOTE: from mergeF
"Recruitment.stem.year","Mortality.stem.year","UnobsAGWPrec.ha","UnobsRecruits.ha","CensusInterval", # NOTE: from tmp, all per ha as PlotArea weighted !
"AGWP.ha","AGWP.ha.year")
#Extract recruitment of monocots in each census
if(sum(Recruits$Monocot)>0 &!is.na(sum(Recruits$Monocot))){
if(rec.meth==100){
palm.rec<-aggregate(AGWPRec/PlotArea ~ PlotViewID+Census.No,data=Recruits[Recruits$Monocot==1,],FUN=sum)
}else{
palm.rec<-aggregate(AGBind/PlotArea ~ PlotViewID+Census.No,data=Recruits[Recruits$Monocot==1,],FUN=sum)
}
names(palm.rec)[3]<-"MonocotRecruits.ha"
SummaryB<-merge(SummaryB,palm.rec,by = c('PlotViewID','Census.No'),all.x=T)
}
if(sum(DeadTrees$Monocot)>0 &!is.na(sum(DeadTrees$Monocot))){
#Mortality of monocots in each census
palm.dead <-aggregate (AGBDead/PlotArea ~ PlotViewID + Census.No, data = AGBData[AGBData$Monocot==1,], FUN=sum )
names(palm.dead)[3]<-"MonocotMortality.ha"
SummaryB<-merge(SummaryB,palm.dead,by = c('PlotViewID','Census.No'),all.x=T)
}
SummaryB
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.