R/timber.area.r

In nuaquilue/QLDM: Quebec Landscape Dynamic Model

#' Sustainable clear-cut area
#'
#' Calculates sustained-yield levels of clear cuts for each management unit in area
#' 
#' @param land A \code{landscape} data frame with forest stand records in rows
#' @param params A list of default parameters generated by the function \code{default.params()} or 
#' a customized list of model parameters
#' 
#' @return A data frame with the number of cells to clear cut per mgmt unit with actual names of management units
#' 
#' @export
#' 
#' @examples
#' data(landscape)
#' params = default.params()
#' cc.area = timber.area(landscape, params) 
#' 

timber.area <- function(land, params){  

  cat("  Timber supply even-aged stands in area", "\n" )
             
  # Initialize empty vector for the clear cut cells 
  cc.cells <- numeric(0)
  n.cc.cells<- numeric(0)
  
  # Name of the management units.
  land2 <- land[!is.na(land$mgmt.unit),]
  land2$even[land2$tsfire==0 & land2$spp!="NonFor"] <- 1 ## but only for forest spp. f
  
  units <- sort(unique(land2$mgmt.unit[!is.na(land2$mgmt.unit)]))
  
  # Harvest rates have to be calculated separately for each management unit:
  unit=units[24] # for testing  unit="2662"
  for(unit in units){  
    # cat(unit, "\n")
    # Separate locations that can be harvested (included) from those that cannot due to environmental or 
    # social constraints (excluded).
    # Some excluded areas are identified directly on the map based on local knowledge.
    # We need to consider excluded cells in some calculations because they contribute to
    # biodiversity objectives (even if they cannot be harvested).
    # Differentiate also between young and mature stands.
    
    s.inc <- length(land2$cell.id[land2$mgmt.unit == unit &  land2$age >= 0 & is.na(land2$exclus)])
    s.ex  <- length(land2$cell.id[land2$mgmt.unit == unit &  land2$age >= 0 & !is.na(land2$exclus)])
    s.inc.mat <- length(land2$cell.id[land2$mgmt.unit == unit &  land2$age >= land2$age.matu & is.na(land2$exclus)])
    s.ex.mat  <- length(land2$cell.id[land2$mgmt.unit == unit &  land2$age >= land2$age.matu & !is.na(land2$exclus)])
    
    # categories of burned area - young (cannot be salvaged) vs mature (can be salvaged)
    
    s.inc.burnt     <- length(land2$cell.id[land2$mgmt.unit == unit &  land2$tsfire==0 & is.na(land2$exclus)])
    s.inc.mat.burnt <- length(land2$cell.id[land2$mgmt.unit == unit &  (land2$age >= land2$age.matu) &  land2$tsfire==0 & is.na(land2$exclus)])
    s.inc.kill     <- length(land2$cell.id[land2$mgmt.unit == unit &  land2$age >= 0 &  land2$tssbw %in% c(0,5) & is.na(land2$exclus)])
    s.inc.mat.kill <- length(land2$cell.id[land2$mgmt.unit == unit &  (land2$age >= land2$age.matu) &  land2$tssbw %in% c(0,5) & is.na(land2$exclus)])
    
    #print(paste("tordeuse",s.inc.kill,s.inc.mat.kill))
    
    # Extract the portion that is managed through even-aged silviculture (clearcutting) based 
    # on species dominance. Some species are mostly managed through even aged silviculture (EPN,
    # SAB, PET, others), the rest through unevenaged silviculture.

   # even <- land2$mgmt.unit == unit & land2$SppGrp %in% c("EPN", "PET", "SAB", "OthCB", "OthCT", "OthDB") & is.na(land2$exclus) & land2$rndm<=0.95
  #  sum(even) 
  #  even[land2$mgmt.unit == unit & land2$SppGrp %in% c("BOJ", "ERS", "OthDT")& is.na(land2$exclus) & land2$rndm>0.95] <- 1

    
    land.ea <- land2[land2$mgmt.unit == unit & land2$even==1,] 
    dim(land.ea)
    # Get the area managed under an even-aged regime
    s.ea <- length(land.ea$cell.id)   
    
    # Area in mature (old) forests that should be maintained in the FMUs in order to meet the conservation target
    target.old.ha  <- params$target.old.pct * (s.inc + s.ex)
    target.old.ha.ea <- max(0, target.old.ha - s.ex.mat)
    target.old.pct.ea <- target.old.ha.ea/s.ea        
    
    # Subset of harvestable (mature even-aged) cells
    land.rec <- land.ea[land.ea$age >= land.ea$age.matu,]
    s.mat <- nrow(land.rec)
    
    #### Determine the sustained yield level
    
    # Number of strata corresponding to the number of different ages of maturity present in
    # each FMU. Only one stratum is used in the current version
    strates <- sort((unique(land.ea$age.matu)))
      # table(land.ea$age.matu)
    
    # Calculation of the expected abundance of harvestable stands during future planning periods, 
    # as the stands that are currently young will age and become harvestable
    recoltable <- matrix(0,length(strates), params$hor.plan)
    recoltable2 <- matrix(0,length(strates), params$hor.plan)
    for(j in 1:length(strates)) { # j=1
      age.mat.stra <- strates[j]
      TSD_strate <- land.ea$age[land.ea$age.matu==strates[j]]
      # maximum theoretical harvestable area per period for each stratum
      recoltable2[j,] <- length(TSD_strate)/(age.mat.stra/5) * (1:params$hor.plan)   
      # Determine the period when maturity will be reached for the different age classes
      for(per in 0:(params$hor.plan-1))  # per=0  
        recoltable[j,per+1] <- sum(TSD_strate >= (age.mat.stra-(per*5)))
      for(per in (age.mat.stra/5):params$hor.plan)
        recoltable[j,per] <- recoltable2[j,per]
    }
    
    # Total harvestable area, all strata combined, minus what has to be kept to satisfy 
    # the old forest target
    recoltable.s <- colSums(recoltable)
    recoltable.s1 <- pmax(0,recoltable.s-target.old.ha.ea)
    recoltable.s2 <- recoltable.s1/(1:params$hor.plan)
    # a priori reduction in maximal allowable harvest level to buffer fire impacts    
    recoltable.s3 <- recoltable.s2 #* a.priori
    # Number of cells to harvest (sustained yield level) 
    n.cc.cells.UA <-  max(0, round(min(recoltable.s3)*1))
    n.cc.cells <- c(n.cc.cells,n.cc.cells.UA)
  }
  
  ## Return number of cells to cut per mgmt unit (with actual names of management units)
  return(data.frame(mgmt.unit=units, x=n.cc.cells))  
  
}