R/IFNmanagement.R

Defines functions IFNmanagement .managePlot

Documented in IFNmanagement

.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))
}
miquelcaceres/IFNdyn documentation built on Feb. 1, 2021, 10:55 a.m.