.managePlot<-function(treeData,
type = "irregular",
thinning = "below", thinningBA = 10, thinningHB = NA, thinningPerc = 40,
finalMeanDBH = 20, finalPerc = "40-60-100",
plantingSpecies = NA, plantingDBH = 5, plantingHeight = 5, plantingDensity = 200,
finalPreviousStage = 0,
removeEmpty = TRUE, verbose = TRUE) {
BAtotal = as.numeric(basalArea(treeData))
HB = hartBeckingIndex(treeData)
if(is.na(HB)) {
print(treeData)
stop("NA Hart-becking index")
}
meanDBH = sum(treeData$N * treeData$DBH)/sum(treeData$N)
treeDataOut = treeData
treeDataCut = treeData
treeDataCut$N = 0
nextPrescription = list(type = type, thinning = thinning, thinningBA = thinningBA, thinningHB = thinningHB, thinningPerc = thinningPerc,
finalMeanDBH = finalMeanDBH, finalPerc = finalPerc,
plantingSpecies = plantingSpecies, plantingDBH = plantingDBH, plantingHeight = plantingHeight, plantingDensity = plantingDensity,
finalPreviousStage = finalPreviousStage)
action = "none"
if(type == "irregular") action = "thinning"
else if(type == "regular") {
if(verbose) cat(paste0(" mean DBH: ", round(meanDBH,1), " threshold ", finalMeanDBH))
if((meanDBH>finalMeanDBH)|| (finalPreviousStage>0)) {
action = "finalcut"
} else {
action = "thinning"
}
}
planting = FALSE
if(verbose) cat(paste0(" action: ", action))
#Irregular management
if(action=="thinning") {
thin = FALSE
if(!is.na(thinningBA)) {
if(verbose) cat(paste0(" BA: ", round(BAtotal,1), " threshold ", round(thinningBA)))
if(BAtotal > thinningBA) thin = TRUE
}
else if(!is.na(thinningHB)) {
if(verbose) cat(paste0(" HB: ", round(HB,1), " threshold ", round(thinningHB)))
if(HB < thinningHB) thin = TRUE
}
if(thin) {
BA2remove = BAtotal*(thinningPerc/100)
Nremoved = rep(0, nrow(treeData))
BAremoved = 0
if(verbose) cat(paste0(", type: ",thinning,", BA to extract: ", round(BA2remove,1)))
if(thinning %in% c("below", "above")) {
o = order(treeData$DBH, decreasing = (thinning=="above"))
cohort = 1
while(BA2remove > 0) {
r = (treeData$DBH[o[cohort]]/200)
n = treeData$N[o[cohort]]
BAcohort = pi*(r^2)*n
if(BAcohort > BA2remove) {
Nremoved[o[cohort]] = BA2remove/(pi*(r^2))
BA2remove = 0
} else {
Nremoved[o[cohort]] = n
BA2remove = BA2remove - BAcohort
}
cohort = cohort + 1
}
}
else if (thinning%in% c("below-intermediate", "above-intermediate")) { #Cut half of target BA from below/above and the other half using trees from all diameters (keep current distribution)
# Remove half as intermediate
sec = pi*(treeData$DBH/200)^2
propN = treeData$N/sum(treeData$N)
HalfBA2remove = BA2remove/2
while(HalfBA2remove > 0) {
if(sum(sec*propN)> HalfBA2remove) { # Correct to avoid cutting above existences
propN = propN*(HalfBA2remove/sum(sec*propN))
}
Nremoved = Nremoved + propN
HalfBA2remove = HalfBA2remove - sum(sec*propN)
}
BA2remove = BA2remove/2
#Remove remaining as below/above
o = order(treeData$DBH, decreasing = (thinning=="above-intermediate"))
cohort = 1
while(BA2remove > 0) {
r = (treeData$DBH[o[cohort]]/200)
n = treeData$N[o[cohort]] - Nremoved[o[cohort]]
BAcohort = pi*(r^2)*n
if(BAcohort > BA2remove) {
Nremoved[o[cohort]] = Nremoved[o[cohort]] + BA2remove/(pi*(r^2))
BA2remove = 0
} else {
Nremoved[o[cohort]] = Nremoved[o[cohort]]+n
BA2remove = BA2remove - BAcohort
}
cohort = cohort + 1
}
}
else if (thinning=="intermediate") { #Cut trees from all diameters (keep current distribution)
sec = pi*(treeData$DBH/200)^2
propN = treeData$N/sum(treeData$N)
while(BA2remove > 0) {
if(sum(sec*propN)> BA2remove) { # Correct to avoid cutting above existences
propN = propN*(BA2remove/sum(sec*propN))
}
Nremoved = Nremoved + propN
BA2remove = BA2remove - sum(sec*propN)
}
}
else {
s = strsplit(thinning,"/")[[1]]
# print(s)
ncl = length(s)
breaks = rep(0, ncl)
props = rep(0, ncl)
for(i in 1:ncl) {
s2 = strsplit(s[i],"-")[[1]]
breaks[i] = as.numeric(s2[1])
props[i] = as.numeric(s2[2])
}
breaks[ncl] = max(breaks[ncl], max(treeData$DBH)+1) # To include largest trees
breaks = c(0, breaks)
# print(breaks)
props = props/sum(props)
cl = as.numeric(cut(treeData$DBH, breaks))
propNmat = matrix(0, nrow(treeData), ncol = ncl)
sec = pi*(treeData$DBH/200)^2
BAcl = rep(0, ncl)
for(i in 1:ncl) {
# print(sum(cl==i))
propNmat[cl==i,i] = treeData$N[cl==i]/sum(treeData$N[cl==i])
BAcl[i] = sum(sec[cl==i]*treeData$N[cl==i])
}
# print(propNmat)
BA2removecl = BA2remove*props
# print(cbind(BA2removecl, BAcl))
BA2removecl = pmin(BA2removecl, BAcl) # Do not allow to cut more than is available for a given class
for(i in 1:ncl) {
propN = propNmat[,i]
while(BA2removecl[i] > 0) {
if(sum(sec*propN)> BA2removecl[i]) { # Correct to avoid cutting above existences
propN = propN*(BA2removecl[i]/sum(sec*propN))
}
Nremoved = Nremoved + propN
BA2removecl[i] = BA2removecl[i] - sum(sec*propN)
}
}
}
treeDataOut$N = treeDataOut$N - Nremoved
treeDataCut$N = treeDataCut$N + Nremoved
} else {
if(verbose) cat(paste0(" aborted"))
treeDataOut = treeData
treeDataCut = treeData[-(1:nrow(treeData)),]
}
}
else if(action=="finalcut") {
#Proportions in different stages
if(verbose) cat(paste0(", previous stage: ",finalPreviousStage))
finalPreviousStage = finalPreviousStage + 1
props = as.numeric(strsplit(finalPerc,"-")[[1]])/100
BA2remove = BAtotal*props[finalPreviousStage]
if(verbose) cat(paste0(", % to extract: ", round(props[finalPreviousStage]*100),", BA to extract: ", round(BA2remove,1)))
Nremoved = rep(0, nrow(treeData))
BAremoved = 0
#Cut from above in final cuts
o = order(treeData$DBH, decreasing = TRUE)
cohort = 1
while((BA2remove > 0) && (cohort<=length(o))) {
# print(BA2remove)
d = treeData$DBH[o[cohort]]
r = (d/200)
n = treeData$N[o[cohort]]
if((d > 12.5) || (!is.na(plantingSpecies))) { # Do not cut regeneration if there is no planting
BAcohort = pi*(r^2)*n
if(BAcohort > BA2remove) {
Nremoved[o[cohort]] = BA2remove/(pi*(r^2))
BA2remove = 0
} else {
Nremoved[o[cohort]] = n
BA2remove = BA2remove - BAcohort
}
}
cohort = cohort + 1
}
treeDataOut$N = treeDataOut$N - Nremoved
treeDataCut$N = treeDataCut$N + Nremoved
if(finalPreviousStage==length(props)) {
finalPreviousStage = 0
if(!is.na(plantingSpecies)) planting = TRUE
}
nextPrescription$finalPreviousStage = finalPreviousStage
if(verbose) cat(paste0(", final stage: ",finalPreviousStage))
}
if(removeEmpty) {
treeDataOut = treeDataOut[treeDataOut$N>0.000001, ]
treeDataCut = treeDataCut[treeDataCut$N>0.000001, ]
treeDataOut = treeDataOut[!is.na(treeDataOut$N),]
treeDataCut = treeDataCut[!is.na(treeDataCut$N),]
}
if(planting) {
action = paste0(action, " + planting")
if(verbose) cat(paste0(", planting"))
nr = nrow(treeDataOut)+1
treeDataOut[nr,] = treeDataCut[1,]
treeDataOut[nr,"Species"] = plantingSpecies
if("Name" %in% names(treeDataOut)) treeDataOut[nr,"Name"] = speciesNames(plantingSpecies)
treeDataOut[nr,"DBH"] = plantingDBH
treeDataOut[nr,"H"] = plantingHeight
treeDataOut[nr,"N"] = plantingDensity
}
if(verbose) cat(".\n")
# print(head(treeDataOut))
# print(head(treeDataCut))
return(list(action = action, treeDataOut = treeDataOut, treeDataCut = treeDataCut, nextPrescription = nextPrescription))
}
#' Management module
#'
#' Application of sylvicultural action as specified in prescription data
#'
#' @param treeData A data frame with tree records in rows and columns 'ID', 'Species', 'DBH' (in cm) and 'N' (ha-1)
#' @param prescriptionData A data frame with plots in rows (row names should be equal to 'ID' in tree data) and prescription options in columns:
#' \itemize{
#' \item{\code{type}: Either 'regular' or 'irregular'.}
#' \item{\code{thinning}: How trees are selected for cutting. Either:
#' \itemize{
#' \item{'above' (i.e. thinning from above).}
#' \item{'below' (i.e. thinning from below).}
#' \item{'intermediate' (i.e. according to the current diameter distribution).}
#' \item{A string specifying diameter breaks and percentages (see details).}
#' }
#' }
#' \item{\code{thinningBA}: Threshold for basal area. Thinning will be applied if the basal area of the stand is larger than the speficied value.}
#' \item{\code{thinningHB}: Threshold for Hart-Becking index (if thinningBA = NA). Thinning will be applied if the Hart-Becking index of the stand is smaller than the speficied threshold value.}
#' \item{\code{thinningPerc}: Percentage of initial basal area to be cut.}
#' \item{\code{finalMeanDBH}: For regular management, threshold of mean diameter (in cm) separating thinnings from final cuts.}
#' \item{\code{finalPrec}: For regular management, string with percentages of basal area to be removed in each final cut ('40-60-100' means a three final cuts,
#' where percentages of basal area to be removed each time is 40\%, 60\% and 100\%).}
#' \item{\code{plantingSpecies}: Tree species code to be planted after the final cut.}
#' \item{\code{plantingDBH}: Tree DBH to be planted after the final cut (in cm).}
#' \item{\code{plantingHeight}: Tree height to be planted after the final cut (in m).}
#' \item{\code{plantingDensity}: Tree density to be planted after the final cut (in ind·ha-1).}
#' \item{\code{finalPreviousStage}: Previous stage in final cut (0 - means thinning interventions, 1 - means first final cut, 2 - means second final cut, and so on).}
#' }
#' @param verbose A flag to indicate console output explaining the silvicultural operations done in each plot.
#'
#' @return A list with the following elements:
#' \itemize{
#' \item{\code{action}: A string of the silvicultural action performed (\code{'none'} for no action).}
#' \item{\code{treeDataOut}: Tree data frame after the application of management actions.}
#' \item{\code{treeDataCut}: Tree data frame with the trees that were removed from the input tree data.}
#' \item{\code{nextPrescription}: Data frame as the input prescription data, but where the stage of final cuts may have been modified for regular management sequential actions.}
#' }
#'
#' @details Final cuts in regular management are done from above (i.e. always starting with the largest trees). Previous cuts should be done from below, but it is left
#' to the user.
#'
#' When specifying thinning distribution in a string '15-60/30-30/40-10' means that diameters that 60\% of basal area to be cut is from trees between 0 to 15 cm,
#' 30\% from trees between 15 and 30 cm and 10\% from trees between 30 and 40 cm.
#'
#'
#' @examples
#' # Load example tree data (two forest plots)
#' data(exampleTreeData)
#' head(exampleTreeData)
#'
#' data(examplePrescriptionData)
#' examplePrescriptionData
#'
#' IFNmanagement(exampleTreeData, examplePrescriptionData)
IFNmanagement<-function(treeData, prescriptionData, verbose = FALSE) {
IDs = row.names(prescriptionData)
# print(IDs)
man = vector("list", length(IDs))
names(man) = IDs
ncut = 0
nout = 0
for(i in 1:length(IDs)) {
tdi = treeData[treeData$ID==IDs[i],]
if(verbose) {
cat(paste0(" Plot ", IDs[i],": "))
}
pdi = as.list(prescriptionData[IDs[i],])
man[[i]] = .managePlot(tdi,
type=pdi$type,
thinning = pdi$thinning, thinningBA = pdi$thinningBA, thinningHB = pdi$thinningHB, thinningPerc = pdi$thinningPerc,
finalMeanDBH = pdi$finalMeanDBH, finalPerc = pdi$finalPerc,
plantingSpecies = pdi$plantingSpecies, plantingDBH = pdi$plantingDBH, plantingHeight = pdi$plantingHeight, plantingDensity = pdi$plantingDensity,
finalPreviousStage = pdi$finalPreviousStage,
verbose = verbose)
ncut = ncut + nrow(man[[i]]$treeDataCut)
nout = nout + nrow(man[[i]]$treeDataOut)
}
treeDataOut = as.data.frame(matrix(nrow = nout, ncol = ncol(treeData)))
names(treeDataOut) = names(treeData)
treeDataCut = as.data.frame(matrix(nrow = ncut, ncol = ncol(treeData)))
names(treeDataCut) = names(treeData)
cntOut = 0
cntCut = 0
nextPrescription = prescriptionData
action = character(0)
for(i in 1:length(IDs)) {
tdo = man[[i]]$treeDataOut
tdc = man[[i]]$treeDataCut
if(nrow(tdc)>0) {
treeDataCut[(cntCut+1):(cntCut + nrow(tdc)),] = tdc
cntCut = cntCut + nrow(tdc)
}
if(nrow(tdo)>0) {
treeDataOut[(cntOut+1):(cntOut + nrow(tdo)),] = tdo
cntOut = cntOut + nrow(tdo)
}
nextPrescription[i,] = man[[i]]$nextPrescription
action = c(action, man[[i]]$action)
}
names(action) = IDs
if(nout>0) row.names(treeDataOut) = 1:nout
if(ncut>0) row.names(treeDataCut) = 1:ncut
return(list(action = action, treeDataOut = treeDataOut, treeDataCut = treeDataCut, nextPrescription = nextPrescription))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.