R/timber.volume.r

In nuaquilue/QLDM: Quebec Landscape Dynamic Model

#' Sustainable clear-cut volume
#'
#' Calculates sustained-yield levels of clear cuts for each management unit in volume
#' 
#' @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
#' @param time.step Number of years of each time step
#' 
#' @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.vol = timber.volume(landscape, params) 
#' 

timber.volume <- function(land, params, time.step=5){  

  cat("  Timber supply even-aged stands in volume", "\n" )
             
  # Initialize empty vector for the clear cut cells 
  cc.cells <- initial.volume<- numeric(0)
  
  # Name of the management units.
  land2 <- land[!is.na(land$mgmt.unit),]
  units <- as.character(sort(unique(land2$mgmt.unit[!is.na(land2$mgmt.unit)])))

  # Harvest rates have to be calculated separately for each management unit:
  for(unit in units){  
    
    # 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$age >= 0 &  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 & 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)])
    
    # 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.
    land.ea1 <- land2[land2$mgmt.unit == unit &  
                      land2$spp %in% c("EPN", "PET", "SAB", "OTH.RES.N", "OTH.RES.S", "OTH.FEU.N") & is.na(land2$exclus),]
    land.ea2 <- land2[land2$mgmt.unit == unit &  
                      land2$spp %in% c("BOJ", "ERS", "OTH.FEU.S") & is.na(land2$exclus),]
    land.ea3 <- land.ea2[runif(nrow(land.ea2))<0.05,]
    land.ea4 <- land.ea1[runif(nrow(land.ea1))<0.95,]
    land.ea <- rbind(land.ea4,land.ea3)
    
    # 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: alternative approach to compute in volume
    strates <- sort((unique(land.ea$age.matu)))
    recoltable <- matrix(0,length(strates), params$hor.plan)
    recoltable2 <- matrix(0,length(strates), params$hor.plan)
    for(j in 1:length(strates)){ 
      
      # Maturity age of the starte
      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
      # VOLUME: It is assumed that communities are exactly at maturity
      land.ea2 <- land.ea[land.ea$age.matu == age.mat.stra,]
      land.ea2$age <- land.ea2$age.matu
      vol.max <- sum(stand.volume(land.ea2))*400
      recoltable2[j,] <- vol.max/(age.mat.stra/time.step) * (1:params$hor.plan)   
      
      # We go back to the initial age for the following calculations
      land.ea2 <- land.ea[land.ea$age.matu == age.mat.stra,]
      
      # Determine the period when maturity will be reached for the different age classes
      for(per in 0:(params$hor.plan-1)){
        # on calcule le volume des peuplements matures 
        if(sum(land.ea2$age>=land.ea2$age.matu)>0){
          vol.act <- sum(stand.volume(land.ea2[land.ea2$age>=land.ea2$age.matu,]))*400
        }
        else{
          vol.act=0
        }
        recoltable[j,per+1] <- vol.act
        # pour chaque periode, on update l'?ge des peuplements pour la p?riopde suivante
        land.ea2$age <- land.ea2$age + time.step
      } 
      for(per in (age.mat.stra/time.step):params$hor.plan)
        recoltable[j,per] <- recoltable2[j,per]
    }
    
    recoltable.s <- colSums(recoltable, na.rm=T)
    recoltable.s1 <- recoltable.s# pmax(0,recoltable.s-target.old.ha.ea)
    recoltable.s2 <- recoltable.s1/(1:params$hor.plan)
    recoltable.s3 <- recoltable.s2 #* a.priori
    vol.UA <-  max(0, round(min(recoltable.s3)*1))
    initial.volume <- c(initial.volume,vol.UA)
  }
  
  ## Return volume to cut per mgmt unit 
  return(data.frame(mgmt.unit=units, x=initial.volume))
  
}