IFNdyn: IFN forest dynamics and scenario analysis

Documented in mergeScenarios

#' Merge scenario results
#'
#' The function merges the set of scenario results given as input. This is useful, for example, to obtain statistics
#' at the level of CCAA from simulation results at the province level.
#'
#' @param x list of scenario results. Each element is the result of calling \code{\link{IFNscenario}}
#'
#' @return An object of class '\code{IFNscenario}' with merged results
mergeScenarios<-function(x) {

  n = length(x)

  #Checkings
  spp = character(0)
  domspp = character(0)
  plotIDs = character(0)
  demandSpp = character(0)
  numSteps = x[[1]]$steps
  stepOffset = x[[1]]$stepOffset
  for(i in 1:n) {
    if(x[[i]]$steps != numSteps) stop("All simulations need to be of the same number of steps.")
    if(x[[i]]$stepOffset != stepOffset) stop("All simulations need to have the same step offset.")
    demandSpp = c(demandSpp, row.names(x[[i]]$extractedDemand))
    spp = c(spp, row.names(x[[i]]$volume$standingSpp))
    domspp = c(domspp, row.names(x[[i]]$volume$standingDomSpp))
    plotIDs = c(plotIDs, row.names(x[[i]]$volume$standingPlot))
  }

  demandSpp = sort(unique(demandSpp))
  spp = as.character(sort(as.numeric(unique(spp))))
  domspp = as.character(sort(as.numeric(unique(domspp))))
  plotIDs = as.character(sort(as.numeric(unique(plotIDs))))

  # init output vars
  plots = 0
  totalArea = 0
  treeDataNew = NULL
  treeDataSequence = NULL
  treeDataDead = NULL
  treeDataCut = NULL

  plotDataSuppl = data.frame(row.names=plotIDs)
  plotDataSuppl$Area = NA
  plotDataSuppl$DominantSpecies = NA
  plotDataSuppl$DemandRow = NA
  plotDataSuppl$PrescriptionRow = NA
  plotDataSuppl$finalPreviousStage = NA

  plotActions = matrix(NA, nrow=length(plotIDs), ncol = numSteps)
  rownames(plotActions) = plotIDs
  colnames(plotActions) = 1:numSteps

  extractedDemand = data.frame(matrix(0, nrow=length(demandSpp), ncol=numSteps), row.names = demandSpp)
  names(extractedDemand)[1:numSteps] = 1:numSteps
  extractedDemand$AnnualMean = 0
  extractedDemand$Demand = 0
  managedArea = data.frame(matrix(0, nrow=length(demandSpp), ncol=numSteps), row.names = demandSpp)
  names(managedArea)[1:numSteps] = 1:numSteps
  managedArea$ManageableArea = 0

  #Standing volume
  standingVolumeSpp = matrix(0, length(spp), numSteps+1)
  rownames(standingVolumeSpp) = spp
  colnames(standingVolumeSpp) = 0:numSteps

  standingVolumePlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingVolumePlot) = plotIDs
  colnames(standingVolumePlot) = 0:numSteps

  standingVolumeDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingVolumeDomSpp) = domspp
  colnames(standingVolumeDomSpp) = 0:numSteps


  #Dead volume
  deadVolumeSpp = matrix(0, length(spp), numSteps)
  rownames(deadVolumeSpp) = spp
  colnames(deadVolumeSpp) = 1:numSteps

  deadVolumePlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadVolumePlot) = plotIDs
  colnames(deadVolumePlot) = 1:numSteps

  deadVolumeDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadVolumeDomSpp) = domspp
  colnames(deadVolumeDomSpp) = 1:numSteps


  #Extracted volume
  extractedVolumeSpp = matrix(0, length(spp), numSteps)
  rownames(extractedVolumeSpp) = spp
  colnames(extractedVolumeSpp) = 1:numSteps

  extractedVolumePlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedVolumePlot) = plotIDs
  colnames(extractedVolumePlot) = 1:numSteps

  extractedVolumeDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedVolumeDomSpp) = domspp
  colnames(extractedVolumeDomSpp) = 1:numSteps

  #Standing density
  standingDensitySpp = matrix(0, length(spp), numSteps+1)
  rownames(standingDensitySpp) = spp
  colnames(standingDensitySpp) = 0:numSteps

  standingDensityPlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingDensityPlot) = plotIDs
  colnames(standingDensityPlot) = 0:numSteps

  standingDensityDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingDensityDomSpp) = domspp
  colnames(standingDensityDomSpp) = 0:numSteps

  #Dead density
  deadDensitySpp = matrix(0, length(spp), numSteps)
  rownames(deadDensitySpp) = spp
  colnames(deadDensitySpp) = 1:numSteps

  deadDensityPlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadDensityPlot) = plotIDs
  colnames(deadDensityPlot) = 1:numSteps

  deadDensityDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadDensityDomSpp) = domspp
  colnames(deadDensityDomSpp) = 1:numSteps


  #Extracted density
  extractedDensitySpp = matrix(0, length(spp), numSteps)
  rownames(extractedDensitySpp) = spp
  colnames(extractedDensitySpp) = 1:numSteps

  extractedDensityPlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedDensityPlot) = plotIDs
  colnames(extractedDensityPlot) = 1:numSteps

  extractedDensityDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedDensityDomSpp) = domspp
  colnames(extractedDensityDomSpp) = 1:numSteps


  #Standing root biomass
  standingRootBiomassSpp = matrix(0, length(spp), numSteps+1)
  rownames(standingRootBiomassSpp) = spp
  colnames(standingRootBiomassSpp) = 0:numSteps

  standingRootBiomassPlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingRootBiomassPlot) = plotIDs
  colnames(standingRootBiomassPlot) = 0:numSteps

  standingRootBiomassDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingRootBiomassDomSpp) = domspp
  colnames(standingRootBiomassDomSpp) = 0:numSteps

  #Standing stem biomass
  standingStemBiomassSpp = matrix(0, length(spp), numSteps+1)
  rownames(standingStemBiomassSpp) = spp
  colnames(standingStemBiomassSpp) = 0:numSteps

  standingStemBiomassPlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingStemBiomassPlot) = plotIDs
  colnames(standingStemBiomassPlot) = 0:numSteps

  standingStemBiomassDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingStemBiomassDomSpp) = domspp
  colnames(standingStemBiomassDomSpp) = 0:numSteps

  #Standing branch biomass
  standingBranchBiomassSpp = matrix(0, length(spp), numSteps+1)
  rownames(standingBranchBiomassSpp) = spp
  colnames(standingBranchBiomassSpp) = 0:numSteps

  standingBranchBiomassPlot = matrix(0, length(plotIDs), numSteps+1)
  rownames(standingBranchBiomassPlot) = plotIDs
  colnames(standingBranchBiomassPlot) = 0:numSteps

  standingBranchBiomassDomSpp = matrix(0, length(domspp), numSteps+1)
  rownames(standingBranchBiomassDomSpp) = domspp
  colnames(standingBranchBiomassDomSpp) = 0:numSteps

  #Dead root biomass
  deadRootBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(deadRootBiomassSpp) = spp
  colnames(deadRootBiomassSpp) = 1:numSteps

  deadRootBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadRootBiomassPlot) = plotIDs
  colnames(deadRootBiomassPlot) = 1:numSteps

  deadRootBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadRootBiomassDomSpp) = domspp
  colnames(deadRootBiomassDomSpp) = 1:numSteps

  #Dead stem biomass
  deadStemBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(deadStemBiomassSpp) = spp
  colnames(deadStemBiomassSpp) = 1:numSteps

  deadStemBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadStemBiomassPlot) = plotIDs
  colnames(deadStemBiomassPlot) = 1:numSteps

  deadStemBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadStemBiomassDomSpp) = domspp
  colnames(deadStemBiomassDomSpp) = 1:numSteps


  #Dead branch biomass
  deadBranchBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(deadBranchBiomassSpp) = spp
  colnames(deadBranchBiomassSpp) = 1:numSteps

  deadBranchBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(deadBranchBiomassPlot) = plotIDs
  colnames(deadBranchBiomassPlot) = 1:numSteps

  deadBranchBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(deadBranchBiomassDomSpp) = domspp
  colnames(deadBranchBiomassDomSpp) = 1:numSteps


  #Extracted root biomass
  extractedRootBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(extractedRootBiomassSpp) = spp
  colnames(extractedRootBiomassSpp) = 1:numSteps

  extractedRootBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedRootBiomassPlot) = plotIDs
  colnames(extractedRootBiomassPlot) = 1:numSteps

  extractedRootBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedRootBiomassDomSpp) = domspp
  colnames(extractedRootBiomassDomSpp) = 1:numSteps

  #Extracted stem biomass
  extractedStemBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(extractedStemBiomassSpp) = spp
  colnames(extractedStemBiomassSpp) = 1:numSteps

  extractedStemBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedStemBiomassPlot) = plotIDs
  colnames(extractedStemBiomassPlot) = 1:numSteps

  extractedStemBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedStemBiomassDomSpp) = domspp
  colnames(extractedStemBiomassDomSpp) = 1:numSteps

  #Extracted branch biomass
  extractedBranchBiomassSpp = matrix(0, length(spp), numSteps)
  rownames(extractedBranchBiomassSpp) = spp
  colnames(extractedBranchBiomassSpp) = 1:numSteps

  extractedBranchBiomassPlot = matrix(0, length(plotIDs), numSteps)
  rownames(extractedBranchBiomassPlot) = plotIDs
  colnames(extractedBranchBiomassPlot) = 1:numSteps

  extractedBranchBiomassDomSpp = matrix(0, length(domspp), numSteps)
  rownames(extractedBranchBiomassDomSpp) = domspp
  colnames(extractedBranchBiomassDomSpp) = 1:numSteps
  # Parse data
  for(i in 1:n) {
    plots = plots + x[[i]]$plots
    totalArea = totalArea + x[[i]]$totalArea

    pdsi = x[[i]]$plotDataSuppl
    plotDataSuppl[row.names(pdsi),] = pdsi
    pai = x[[i]]$plotActions
    plotActions[row.names(pai),] = pai

    dspp = row.names(x[[i]]$extractedDemand)
    extractedDemand[dspp,1:numSteps] = extractedDemand[dspp,1:numSteps] + x[[i]]$extractedDemand[,1:numSteps]
    extractedDemand[dspp, "AnnualMean"] = extractedDemand[dspp, "AnnualMean"] + x[[i]]$extractedDemand[,"AnnualMean"]
    extractedDemand[dspp, "Demand"] = extractedDemand[dspp, "Demand"] + x[[i]]$extractedDemand[,"Demand"]
    managedArea[dspp,1:numSteps] = managedArea[dspp,1:numSteps] + x[[i]]$managedArea[,1:numSteps]
    managedArea[dspp, "ManageableArea"] = managedArea[dspp, "ManageableArea"] + x[[i]]$managedArea[,"ManageableArea"]

    cspp = row.names(x[[i]]$volume$standingSpp)
    standingVolumeSpp[cspp, ] = standingVolumeSpp[cspp, ] + x[[i]]$volume$standingSpp
    standingDensitySpp[cspp, ] = standingDensitySpp[cspp, ] + x[[i]]$density$standingSpp
    standingRootBiomassSpp[cspp, ] = standingRootBiomassSpp[cspp, ] + x[[i]]$rootBiomass$standingSpp
    standingStemBiomassSpp[cspp, ] = standingStemBiomassSpp[cspp, ] + x[[i]]$stemBiomass$standingSpp
    standingBranchBiomassSpp[cspp, ] = standingBranchBiomassSpp[cspp, ] + x[[i]]$branchBiomass$standingSpp

    deadVolumeSpp[cspp, ] = deadVolumeSpp[cspp, ] + x[[i]]$volume$deadSpp
    deadDensitySpp[cspp, ] = deadDensitySpp[cspp, ] + x[[i]]$density$deadSpp
    deadRootBiomassSpp[cspp, ] = deadRootBiomassSpp[cspp, ] + x[[i]]$rootBiomass$deadSpp
    deadStemBiomassSpp[cspp, ] = deadStemBiomassSpp[cspp, ] + x[[i]]$stemBiomass$deadSpp
    deadBranchBiomassSpp[cspp, ] = deadBranchBiomassSpp[cspp, ] + x[[i]]$branchBiomass$deadSpp

    extractedVolumeSpp[cspp, ] = extractedVolumeSpp[cspp, ] + x[[i]]$volume$extractedSpp
    extractedDensitySpp[cspp, ] = extractedDensitySpp[cspp, ] + x[[i]]$density$extractedSpp
    extractedRootBiomassSpp[cspp, ] = extractedRootBiomassSpp[cspp, ] + x[[i]]$rootBiomass$extractedSpp
    extractedStemBiomassSpp[cspp, ] = extractedStemBiomassSpp[cspp, ] + x[[i]]$stemBiomass$extractedSpp
    extractedBranchBiomassSpp[cspp, ] = extractedBranchBiomassSpp[cspp, ] + x[[i]]$branchBiomass$extractedSpp

    cpl = row.names(x[[i]]$volume$standingPlot)
    standingVolumePlot[cpl, ] = standingVolumePlot[cpl, ] + x[[i]]$volume$standingPlot
    standingDensityPlot[cpl, ] = standingDensityPlot[cpl, ] + x[[i]]$density$standingPlot
    standingRootBiomassPlot[cpl, ] = standingRootBiomassPlot[cpl, ] + x[[i]]$rootBiomass$standingPlot
    standingStemBiomassPlot[cpl, ] = standingStemBiomassPlot[cpl, ] + x[[i]]$stemBiomass$standingPlot
    standingBranchBiomassPlot[cpl, ] = standingBranchBiomassPlot[cpl, ] + x[[i]]$branchBiomass$standingPlot

    deadVolumePlot[cpl, ] = deadVolumePlot[cpl, ] + x[[i]]$volume$deadPlot
    deadDensityPlot[cpl, ] = deadDensityPlot[cpl, ] + x[[i]]$density$deadPlot
    deadRootBiomassPlot[cpl, ] = deadRootBiomassPlot[cpl, ] + x[[i]]$rootBiomass$deadPlot
    deadStemBiomassPlot[cpl, ] = deadStemBiomassPlot[cpl, ] + x[[i]]$stemBiomass$deadPlot
    deadBranchBiomassPlot[cpl, ] = deadBranchBiomassPlot[cpl, ] + x[[i]]$branchBiomass$deadPlot

    extractedVolumePlot[cpl, ] = extractedVolumePlot[cpl, ] + x[[i]]$volume$extractedPlot
    extractedDensityPlot[cpl, ] = extractedDensityPlot[cpl, ] + x[[i]]$density$extractedPlot
    extractedRootBiomassPlot[cpl, ] = extractedRootBiomassPlot[cpl, ] + x[[i]]$rootBiomass$extractedPlot
    extractedStemBiomassPlot[cpl, ] = extractedStemBiomassPlot[cpl, ] + x[[i]]$stemBiomass$extractedPlot
    extractedBranchBiomassPlot[cpl, ] = extractedBranchBiomassPlot[cpl, ] + x[[i]]$branchBiomass$extractedPlot


    cdomspp = row.names(x[[i]]$volume$standingDomSpp)
    standingVolumeDomSpp[cdomspp, ] = standingVolumeDomSpp[cdomspp, ] + x[[i]]$volume$standingDomSpp
    standingDensityDomSpp[cdomspp, ] = standingDensityDomSpp[cdomspp, ] + x[[i]]$density$standingDomSpp
    standingRootBiomassDomSpp[cdomspp, ] = standingRootBiomassDomSpp[cdomspp, ] + x[[i]]$rootBiomass$standingDomSpp
    standingStemBiomassDomSpp[cdomspp, ] = standingStemBiomassDomSpp[cdomspp, ] + x[[i]]$stemBiomass$standingDomSpp
    standingBranchBiomassDomSpp[cdomspp, ] = standingBranchBiomassDomSpp[cdomspp, ] + x[[i]]$branchBiomass$standingDomSpp

    deadVolumeDomSpp[cdomspp, ] = deadVolumeDomSpp[cdomspp, ] + x[[i]]$volume$deadDomSpp
    deadDensityDomSpp[cdomspp, ] = deadDensityDomSpp[cdomspp, ] + x[[i]]$density$deadDomSpp
    deadRootBiomassDomSpp[cdomspp, ] = deadRootBiomassDomSpp[cdomspp, ] + x[[i]]$rootBiomass$deadDomSpp
    deadStemBiomassDomSpp[cdomspp, ] = deadStemBiomassDomSpp[cdomspp, ] + x[[i]]$stemBiomass$deadDomSpp
    deadBranchBiomassDomSpp[cdomspp, ] = deadBranchBiomassDomSpp[cdomspp, ] + x[[i]]$branchBiomass$deadDomSpp

    extractedVolumeDomSpp[cdomspp, ] = extractedVolumeDomSpp[cdomspp, ] + x[[i]]$volume$extractedDomSpp
    extractedDensityDomSpp[cdomspp, ] = extractedDensityDomSpp[cdomspp, ] + x[[i]]$density$extractedDomSpp
    extractedRootBiomassDomSpp[cdomspp, ] = extractedRootBiomassDomSpp[cdomspp, ] + x[[i]]$rootBiomass$extractedDomSpp
    extractedStemBiomassDomSpp[cdomspp, ] = extractedStemBiomassDomSpp[cdomspp, ] + x[[i]]$stemBiomass$extractedDomSpp
    extractedBranchBiomassDomSpp[cdomspp, ] = extractedBranchBiomassDomSpp[cdomspp, ] + x[[i]]$branchBiomass$extractedDomSpp

    if(!is.null(treeDataNew)) {
      treeDataNew = rbind(treeDataNew, x[[i]]$treeDataNew)
    } else {
      treeDataNew = x[[i]]$treeDataNew
    }
    if("treeDataSequence" %in% names(x[[i]])) {
      if(!is.null(treeDataSequence)) {
        treeDataSequence = rbind(treeDataSequence, x[[i]]$treeDataSequence)
      } else {
        treeDataSequence = x[[i]]$treeDataSequence
      }
    }
    if("treeDataDead" %in% names(x[[i]])) {
      if(!is.null(treeDataDead)) {
        treeDataDead = rbind(treeDataDead, x[[i]]$treeDataDead)
      } else {
        treeDataDead = x[[i]]$treeDataDead
      }
    }
    if("treeDataCut" %in% names(x[[i]])) {
      if(!is.null(treeDataCut)) {
        treeDataCut = rbind(treeDataCut, x[[i]]$treeDataCut)
      } else {
        treeDataCut = x[[i]]$treeDataCut
      }
    }
  }
  extractedDemand$Coverage = 100*extractedDemand$AnnualMean/extractedDemand$Demand



  #Calculate variation
  # variation_t = standing_t - standing_(t-1)
  variationVolumeSpp = standingVolumeSpp[,2:(numSteps+1), drop=FALSE] - standingVolumeSpp[,1:(numSteps), drop=FALSE]
  variationVolumePlot = standingVolumePlot[,2:(numSteps+1), drop=FALSE] - standingVolumePlot[,1:(numSteps), drop=FALSE]
  variationVolumeDomSpp = standingVolumeDomSpp[,2:(numSteps+1), drop=FALSE] - standingVolumeDomSpp[,1:(numSteps), drop=FALSE]

  variationDensitySpp = standingDensitySpp[,2:(numSteps+1), drop=FALSE] - standingDensitySpp[,1:(numSteps), drop=FALSE]
  variationDensityPlot = standingDensityPlot[,2:(numSteps+1), drop=FALSE] - standingDensityPlot[,1:(numSteps), drop=FALSE]
  variationDensityDomSpp = standingDensityDomSpp[,2:(numSteps+1), drop=FALSE] - standingDensityDomSpp[,1:(numSteps), drop=FALSE]

  variationRootBiomassSpp = standingRootBiomassSpp[,2:(numSteps+1), drop=FALSE] - standingRootBiomassSpp[,1:(numSteps), drop=FALSE]
  variationRootBiomassPlot = standingRootBiomassPlot[,2:(numSteps+1), drop=FALSE] - standingRootBiomassPlot[,1:(numSteps), drop=FALSE]
  variationRootBiomassDomSpp = standingRootBiomassDomSpp[,2:(numSteps+1), drop=FALSE] - standingRootBiomassDomSpp[,1:(numSteps), drop=FALSE]

  variationStemBiomassSpp = standingStemBiomassSpp[,2:(numSteps+1), drop=FALSE] - standingStemBiomassSpp[,1:(numSteps), drop=FALSE]
  variationStemBiomassPlot = standingStemBiomassPlot[,2:(numSteps+1), drop=FALSE] - standingStemBiomassPlot[,1:(numSteps), drop=FALSE]
  variationStemBiomassDomSpp = standingStemBiomassDomSpp[,2:(numSteps+1), drop=FALSE] - standingStemBiomassDomSpp[,1:(numSteps), drop=FALSE]

  variationBranchBiomassSpp = standingBranchBiomassSpp[,2:(numSteps+1), drop=FALSE] - standingBranchBiomassSpp[,1:(numSteps), drop=FALSE]
  variationBranchBiomassPlot = standingBranchBiomassPlot[,2:(numSteps+1), drop=FALSE] - standingBranchBiomassPlot[,1:(numSteps), drop=FALSE]
  variationBranchBiomassDomSpp = standingBranchBiomassDomSpp[,2:(numSteps+1), drop=FALSE] - standingBranchBiomassDomSpp[,1:(numSteps), drop=FALSE]

  #Calculate stand growth (tree growth + mortality + recruitment)
  # Stand growth_t = standing_t - postmanagement_t = standing_t - (standing_(t-1) - extracted_t) = variation_t + extracted_t
  standgrowthVolumeSpp = variationVolumeSpp+extractedVolumeSpp
  standgrowthVolumePlot = variationVolumePlot+extractedVolumePlot
  standgrowthVolumeDomSpp = variationVolumeDomSpp+extractedVolumeDomSpp

  standgrowthDensitySpp = variationDensitySpp+extractedDensitySpp
  standgrowthDensityPlot = variationDensityPlot+extractedDensityPlot
  standgrowthDensityDomSpp = variationDensityDomSpp+extractedDensityDomSpp

  standgrowthRootBiomassSpp = variationRootBiomassSpp+extractedRootBiomassSpp
  standgrowthRootBiomassPlot = variationRootBiomassPlot+extractedRootBiomassPlot
  standgrowthRootBiomassDomSpp = variationRootBiomassDomSpp+extractedRootBiomassDomSpp

  standgrowthStemBiomassSpp = variationStemBiomassSpp+extractedStemBiomassSpp
  standgrowthStemBiomassPlot = variationStemBiomassPlot+extractedStemBiomassPlot
  standgrowthStemBiomassDomSpp = variationStemBiomassDomSpp+extractedStemBiomassDomSpp

  standgrowthBranchBiomassSpp = variationBranchBiomassSpp+extractedBranchBiomassSpp
  standgrowthBranchBiomassPlot = variationBranchBiomassPlot+extractedBranchBiomassPlot
  standgrowthBranchBiomassDomSpp = variationBranchBiomassDomSpp+extractedBranchBiomassDomSpp

  #Stats per hectare
  volumeStats = data.frame(standing = colSums(standingVolumeSpp)/totalArea,
                           variation = c(NA,colSums(variationVolumeSpp)/totalArea),
                           standgrowth = c(NA,colSums(standgrowthVolumeSpp)/totalArea),
                           dead = c(NA,colSums(deadVolumeSpp)/totalArea),
                           extracted = c(NA,colSums(extractedVolumeSpp)/totalArea))
  densityStats = data.frame(standing = colSums(standingDensitySpp)/totalArea,
                            variation = c(NA,colSums(variationDensitySpp)/totalArea),
                            standgrowth = c(NA,colSums(standgrowthDensitySpp)/totalArea),
                            dead = c(NA,colSums(deadDensitySpp)/totalArea),
                            extracted = c(NA,colSums(extractedDensitySpp)/totalArea))
  rootBiomassStats = data.frame(standing = colSums(standingRootBiomassSpp)/totalArea,
                                variation = c(NA,colSums(variationRootBiomassSpp)/totalArea),
                                standgrowth = c(NA,colSums(standgrowthRootBiomassSpp)/totalArea),
                                dead = c(NA,colSums(deadRootBiomassSpp)/totalArea),
                                extracted = c(NA,colSums(extractedRootBiomassSpp)/totalArea))
  stemBiomassStats = data.frame(standing = colSums(standingStemBiomassSpp)/totalArea,
                                variation = c(NA,colSums(variationStemBiomassSpp)/totalArea),
                                standgrowth = c(NA,colSums(standgrowthStemBiomassSpp)/totalArea),
                                dead = c(NA,colSums(deadStemBiomassSpp)/totalArea),
                                extracted = c(NA,colSums(extractedStemBiomassSpp)/totalArea))
  branchBiomassStats = data.frame(standing = colSums(standingBranchBiomassSpp)/totalArea,
                                  variation = c(NA,colSums(variationBranchBiomassSpp)/totalArea),
                                  standgrowth = c(NA,colSums(standgrowthBranchBiomassSpp)/totalArea),
                                  dead = c(NA,colSums(deadBranchBiomassSpp)/totalArea),
                                  extracted = c(NA,colSums(extractedBranchBiomassSpp)/totalArea))

  volume = list(summary = volumeStats,
                standingSpp = standingVolumeSpp,
                standingPlot = standingVolumePlot,
                standingDomSpp = standingVolumeDomSpp,
                variationSpp = variationVolumeSpp,
                variationPlot = variationVolumePlot,
                variationDomSpp = variationVolumeDomSpp,
                standgrowthSpp = standgrowthVolumeSpp,
                standgrowthPlot = standgrowthVolumePlot,
                standgrowthDomSpp = standgrowthVolumeDomSpp,
                deadSpp = deadVolumeSpp,
                deadPlot = deadVolumePlot,
                deadDomSpp = deadVolumeDomSpp,
                extractedSpp = extractedVolumeSpp,
                extractedPlot = extractedVolumePlot,
                extractedDomSpp = extractedVolumeDomSpp)
  density = list(summary = densityStats,
                 standingSpp = standingDensitySpp,
                 standingPlot = standingDensityPlot,
                 standingDomSpp = standingDensityDomSpp,
                 variationSpp = variationDensitySpp,
                 variationPlot = variationDensityPlot,
                 variationDomSpp = variationDensityDomSpp,
                 standgrowthSpp = standgrowthDensitySpp,
                 standgrowthPlot = standgrowthDensityPlot,
                 standgrowthDomSpp = standgrowthDensityDomSpp,
                 deadSpp = deadDensitySpp,
                 deadPlot = deadDensityPlot,
                 deadDomSpp = deadDensityDomSpp,
                 extractedSpp = extractedDensitySpp,
                 extractedPlot = extractedDensityPlot,
                 extractedDomSpp = extractedDensityDomSpp)
  rootBiomass = list(summary = rootBiomassStats,
                     standingSpp = standingRootBiomassSpp,
                     standingPlot = standingRootBiomassPlot,
                     standingDomSpp = standingRootBiomassDomSpp,
                     variationSpp = variationRootBiomassSpp,
                     variationPlot = variationRootBiomassPlot,
                     variationDomSpp = variationRootBiomassDomSpp,
                     standgrowthSpp = standgrowthRootBiomassSpp,
                     standgrowthPlot = standgrowthRootBiomassPlot,
                     standgrowthDomSpp = standgrowthRootBiomassDomSpp,
                     deadSpp = deadRootBiomassSpp,
                     deadPlot = deadRootBiomassPlot,
                     deadDomSpp = deadRootBiomassDomSpp,
                     extractedSpp = extractedRootBiomassSpp,
                     extractedPlot = extractedRootBiomassPlot,
                     extractedDomSpp = extractedRootBiomassDomSpp)
  stemBiomass = list(summary = stemBiomassStats,
                     standingSpp = standingStemBiomassSpp,
                     standingPlot = standingStemBiomassPlot,
                     standingDomSpp = standingStemBiomassDomSpp,
                     variationSpp = variationStemBiomassSpp,
                     variationPlot = variationStemBiomassPlot,
                     variationDomSpp = variationStemBiomassDomSpp,
                     standgrowthSpp = standgrowthStemBiomassSpp,
                     standgrowthPlot = standgrowthStemBiomassPlot,
                     standgrowthDomSpp = standgrowthStemBiomassDomSpp,
                     deadSpp = deadStemBiomassSpp,
                     deadPlot = deadStemBiomassPlot,
                     deadDomSpp = deadStemBiomassDomSpp,
                     extractedSpp = extractedStemBiomassSpp,
                     extractedPlot = extractedStemBiomassPlot,
                     extractedDomSpp = extractedStemBiomassDomSpp)
  branchBiomass = list(summary = branchBiomassStats,
                       standingSpp = standingBranchBiomassSpp,
                       standingPlot = standingBranchBiomassPlot,
                       standingDomSpp = standingBranchBiomassDomSpp,
                       variationSpp = variationBranchBiomassSpp,
                       variationPlot = variationBranchBiomassPlot,
                       variationDomSpp = variationBranchBiomassDomSpp,
                       standgrowthSpp = standgrowthBranchBiomassSpp,
                       standgrowthPlot = standgrowthBranchBiomassPlot,
                       standgrowthDomSpp = standgrowthBranchBiomassDomSpp,
                       deadSpp = deadBranchBiomassSpp,
                       deadPlot = deadBranchBiomassPlot,
                       deadDomSpp = deadBranchBiomassDomSpp,
                       extractedSpp = extractedBranchBiomassSpp,
                       extractedPlot = extractedBranchBiomassPlot,
                       extractedDomSpp = extractedBranchBiomassDomSpp)

  res = list(plots = plots,
             steps = numSteps,
             stepOffset = stepOffset,
             totalArea = totalArea,
             plotDataSuppl = plotDataSuppl,
             plotActions = plotActions,
             extractedDemand=extractedDemand,
             managedArea = managedArea,
             volume=volume,
             density=density,
             rootBiomass = rootBiomass,
             stemBiomass = stemBiomass,
             branchBiomass = branchBiomass,
             treeDataNew = treeDataNew,
             treeDataSequence = treeDataSequence,
             treeDataDead = treeDataDead,
             treeDataCut = treeDataCut)
  class(res) <- "IFNscenario"
  return(res)
}