#' IFN scenario simulation engine
#'
#' Performs simulations of forest dynamics including management prescriptions (on a demand-basis) for a set of forest plots
#'
#' @param treeData A data frame with tree records in rows and columns 'ID', 'Species', 'DBH' (in cm), 'H (in m)' and 'N' (ha-1)
#' @param plotData A data frame with plot 'ID' in rows and columns 'SWHC' (mm), 'Rad' (MJ), 'Temp' (degrees C), 'Prec' (mm), 'PET' (mm), 'slope' (degrees), 'Managed' (boolean) and 'Estrato' (stratum code).
#' It can also include a column 'Area' (ha) with the surface corresponding to each forest plot. In this case, 'strata' will not be used.
#' @param plotDataDyn A data frame with dynamic plot data, with columns 'ID', 'Step', 'Temp' (degrees C), 'Prec' (mm) and 'PET' (mm)
#' @param prescription A data frame with plots in rows (row names should be equal to 'ID' in tree data) and prescription options in columns (see \code{\link{IFNmanagement}})
#' @param strata A data frame with surface areas per stratum (see \code{\link{estratosIFN3}}). Used to determine the surface corresponding to each plot.
#' @param demand A data frame with annual volume to be extracted by species (see \code{\link{defaultDemand}})
#' @param nb An object of class \code{listw} (see functions \code{\link{IFNknn}} and \code{\link{nb2listw}}).
#' @param numSteps Number of steps to simulate (1 step = 10 years)
#' @param demandOffset Initial demand offset (for simulations continuing previous ones).
#' @param stepOffset Offset of time steps to read dynamic plot data (for simulations continuing previous ones).
#' @param as.CO2 Boolean flag to indicate output as kg of CO2 instead of kg of dry weight (see \code{\link{IFNbiomass}}).
#' @param sequence Boolean flag to indicate that the sequence of tree lists by step should be returned.
#' If \code{sequence = TRUE}, an additional data frame is included in the output, called \code{treeDataSequence}.
#' This table includes a column 'Step' to indicate the simulation step.
#' @param dead.table Boolean flag to indicate that a table of trees dead during the simulation is also desired.
#' This table includes a column 'Step' to indicate at which step the trees died.
#' @param cut.table Boolean flag to indicate that a table of trees cut during the simulation is also desired.
#' This table includes a column 'Step' to indicate at which step the trees were cut.
#' @param verbose Boolean flag to indicate extra console output (reporting management interventions in each plot)
#' @param ... Additional options to be passed to function \code{\link{IFNdynamics}}
#'
#' @return A list with the following items:
#' \itemize{
#' \item{\code{plots}: number of IFN plots included in the simulation.}
#' \item{\code{totalArea}: total area (in ha) that the set of plots represent.}
#' \item{\code{plotDataSuppl}: A data frame with supplementary plot data, with columns:
#' \itemize{
#' \item{\code{Area}: Area (in ha) represented by each plot.}
#' \item{\code{DominantSpecies}: Code corresponding to the dominant species of each plot.}
#' \item{\code{DemandRow}: Demand species row index corresponding to each plot.}
#' \item{\code{PrescriptionRow}: Prescription row index corresponding to each plot.}
#' \item{\code{finalPreviousStage}: Previous stage of the final cuts (for plots with regular management).}
#' }
#' }
#' \item{\code{steps}: number of 10-yr steps simulated.}
#' \item{\code{extractedDemand}: a data frame with the volume (m3) extracted for each demand species and time step.
#' The three last columns include: \code{AnnualDemand} the mean annual extraction (simulated), \code{Demand} the annual demand (input), and \code{Coverage} the coverage of the demand (in percent), for each species.
#' The last row ('Total') reports statistics across species.}
#' \item{\code{managedArea}: a data frame with the area (ha) managed for each demand species and time step.
#' The last row ('Total') reports statistics across species.}
#' \item{\code{volume, density, rootBiomass, stemBiomass, branchBiomass}: simulation results by variable. Volume is in m3, density in individuals and biomass in Mg.
#' Each of these elements is in turn a list with the following elements:
#' \itemize{
#' \item{\code{summary}: data frame with a summary of results for each time step (standing, variation, dead and extracted). Here values are given per hectare.}
#' \item{\code{standingSpp}: volume, biomass or density of living trees, per species and time step.}
#' \item{\code{standingPlot}: volume, biomass or density of living trees, per plot and time step.}
#' \item{\code{standingDomSpp}: volume, biomass or density of living trees, per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#' \item{\code{variationSpp}: variation (i.e. difference between the current time step and the previous one) in volume, biomass or density of living trees, per species and time step.}
#' \item{\code{variationPlot}: variation (i.e. difference between the current time step and the previous one) in volume, biomass or density of living trees, per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#' \item{\code{variationDomSpp}: variation (i.e. difference between the current time step and the previous one) in volume, biomass or density of living trees, per plot and time step.}
#' \item{\code{standgrowthSpp}: stand growth in volume, biomass or density of living trees, per species and time step.}
#' \item{\code{standgrowthPlot}: stand growth in volume, biomass or density of living trees, per plot and time step.}
#' \item{\code{standgrowthDomSpp}: stand growth in volume, biomass or density of living trees, per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#' \item{\code{deadSpp}: volume, biomass or density of dead trees (i.e. dead during the corresponding time step), per species and time step.}
#' \item{\code{deadPlot}: volume, biomass or density of dead trees (i.e. dead during the corresponding time step), per plot and time step.}
#' \item{\code{deadDomSpp}: volume, biomass or density of dead trees (i.e. dead during the corresponding time step), per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#' \item{\code{extractedSpp}: volume, biomass or density of trees extracted (during the corresponding time step), per species and time step.}
#' \item{\code{extractedPlot}: volume, biomass or density of trees extracted (during the corresponding time step), per plot and time step.}
#' \item{\code{extractedDomSpp}: volume, biomass or density of trees extracted (during the corresponding time step), per dominant species (i.e. pooling all plots with the same dominant species) and time step.}
#' }
#' }
#' \item{\code{treeDataNew}: tree data frame at the end of the simulated period.}
#' \item{\code{treeDataSequence}: tree data frame for the whole sequence of steps (if \code{sequence = TRUE}).}
#' \item{\code{treeDataDead}: tree data frame with dead trees for every step (if \code{dead.table = TRUE}).}
#' \item{\code{treeDataCut}: tree data frame with cut trees for every step (if \code{cut.table = TRUE}).}
#' }
#'
#' @examples
#'
#' \dontrun{
#' library(IFNread)
#' data(defaultDemand)
#' data(defaultPrescription)
#' data(estratosIFN3)
#' data(plotDataIFN)
#'
#' accessfilesIFN3 = list.files("D:/Recerca/Datasets/IFN/Sources/IFN3/MSAccess/",
#' pattern = "*[.]accdb", full.names = T)
#'
#' prov = 27 # Lugo
#'
#' treeData = readPiesMayoresIFN3(accessfilesIFN3[prov])
#'
#' IFNscenario(treeData, plotDataIFN, defaultPrescription,
#' estratosIFN3[[prov]], defaultDemand[[prov]], numSteps = 3)
#' }
#'
IFNscenario<-function(treeData, plotData, prescription, strata, demand, nb = NULL, numSteps = 1,
plotDataDyn = NULL, demandOffset = NULL, stepOffset = 0,
as.CO2 = TRUE, sequence = FALSE, dead.table = FALSE, cut.table = FALSE,
verbose = FALSE, ...) {
scen = NULL
plotDataSuppl = NULL
if(inherits(treeData, "data.frame")) {
cat(paste("A. Data preparation...\n"))
#Check variable type
treeData$ID = as.character(treeData$ID)
treeData$Species = as.character(treeData$Species)
selDBH = is.na(treeData$DBH)
if(sum(selDBH)>0) {
cat(paste0(" Removing ", sum(selDBH), " records out of ", length(selDBH)," because of missing diameter.\n"))
treeData = treeData[!selDBH,]
}
selDBH = (treeData$DBH==0)
if(sum(selDBH)>0) {
cat(paste0(" Removing ", sum(selDBH), " records out of ", length(selDBH)," because of zero diameter.\n"))
treeData = treeData[!selDBH,]
}
selH = is.na(treeData$H)
if(sum(selH)>0) {
cat(paste0(" Removing ", sum(selH), " records out of ", length(selH)," because of missing height.\n"))
treeData = treeData[!selH,]
}
selH = (treeData$H==0)
if(sum(selH)>0) {
cat(paste0(" Removing ", sum(selH), " records out of ", length(selH)," because of zero height.\n"))
treeData = treeData[!selH,]
}
selSpecies = is.na(treeData$Species)
if(sum(selSpecies)>0) {
cat(paste0(" Removing ", sum(selSpecies), " records out of ",length(selSpecies)," because of missing species.\n"))
treeData = treeData[!selSpecies,]
}
selSpecies = !(treeData$Species %in% biomass_species_match$ID_TAXON)
if(sum(selSpecies)>0) {
wc = unique(treeData$Species[selSpecies])
cat(paste0(" Removing ", sum(selSpecies), " records out of ",length(selSpecies)," because of wrong species code(s): ",paste0(wc,collapse=","),".\n"))
treeData = treeData[!selSpecies,]
}
}
else if(inherits(treeData, "IFNscenario")) {
cat(paste("A. Data loading from previous scenario...\n"))
scen = treeData
plotDataSuppl = scen$plotDataSuppl
treeData = scen$treeDataNew
demandOffset = scen$demandOffset
}
#Check plotData
selSWHC = is.na(plotData$SWHC)
if(sum(selSWHC)>0) {
cat(paste0(" Removing ", sum(selSWHC), " plots from plotData because of missing 'SWHC'.\n"))
plotData = plotData[!selSWHC,]
}
selRad = is.na(plotData$Rad)
if(sum(selRad)>0) {
cat(paste0(" Removing ", sum(selRad), " plots from plotData because of missing 'Rad'.\n"))
plotData = plotData[!selRad,]
}
selTemp = is.na(plotData$Temp)
if(sum(selTemp)>0) {
cat(paste0(" Removing ", sum(selTemp), " plots from plotData because of missing 'Temp'.\n"))
plotData = plotData[!selTemp,]
}
selSlope = is.na(plotData$slope)
if(sum(selSlope)>0) {
cat(paste0(" Removing ", sum(selSlope), " plots from plotData because of missing 'slope'.\n"))
plotData = plotData[!selSlope,]
}
selEstrato = is.na(plotData$Estrato)
if(sum(selEstrato)>0) {
cat(paste0(" Removing ", sum(selEstrato), " plots from plotData because of missing 'estrato'.\n"))
plotData = plotData[!selEstrato,]
}
#Check plot codes
plotIDs = unique(treeData$ID)
sel = (plotIDs %in% row.names(plotData))
nmis = sum(!sel)
if(nmis>0) {
cat(paste(" Removing tree data for", nmis, "plots not found in 'plotData'.\n"))
plotIDs = plotIDs[sel]
treeData = treeData[treeData$ID %in% plotIDs,]
}
cat(paste(" Subsetting", length(plotIDs), "plots from 'plotData'.\n"))
plotData = plotData[row.names(plotData) %in% plotIDs,] # Subset records
plotData = plotData[plotIDs,] # Sort plot data records
if(!is.null(plotDataDyn)) {
cat(paste(" Subsetting", length(plotIDs), "plots from 'plotDataDyn'.\n"))
plotDataDyn = plotDataDyn[plotDataDyn$ID %in% plotIDs,] # Subset records
}
#If 'Area' is not defined in 'plotData', determine area (hectares) corresponding to each plot according to IFN strata definition
nPar = table(plotData$Estrato)
if(!("Area" %in% names(plotData))) {
cat(paste(" Calculating plot area from strata.\n"))
plotData$Area = NA
for(i in 1:length(plotIDs)) {
plotData[i, "Area"] = strata[plotData[i, "Estrato"],"Area"]/nPar[plotData[i, "Estrato"]]
}
} else {
cat(paste(" Using plot area found in 'plotData'.\n"))
}
selArea = is.na(plotData$Area)
if(sum(selArea)>0) {
cat(paste(" Removing data for", sum(selArea), "plots with missing plot area.\n"))
ids = row.names(plotData)[selArea]
plotData = plotData[!selArea, ]
treeData = treeData[!(treeData$ID %in% ids),]
plotIDs = plotIDs[!(plotIDs %in% ids)]
}
#For all plots, determine dominant species (may be missing)
BA = basalArea(treeData, bySpecies = T)
plotData$DominantSpecies = colnames(BA)[apply(BA, 1, which.max)]
if(!is.null(plotDataSuppl)) {
cat(paste(" Replacing dominant species from input scenario.\n"))
plotData[row.names(plotDataSuppl), "DominantSpecies"] = plotDataSuppl$DominantSpecies
}
#Determine prescription and demand rows for each forest plot
if(!("PrescriptionRow" %in% names(plotData))) {
cat(paste(" Initializing prescription row.\n"))
plotData$PrescriptionRow = NA
pr_s = strsplit(as.character(prescription$Espcodes), split=",")
getPrescriptionRow<-function(sp) {
for(i in 1:length(pr_s)) if(as.character(sp) %in% pr_s[[i]]) return(i)
return(NA)
}
for(i in 1:length(plotIDs)) {
plotData[i, "PrescriptionRow"] = getPrescriptionRow(plotData[i,"DominantSpecies"])
}
} else {
cat(paste(" Using prescription row found in 'plotData'.\n"))
}
if(!is.null(plotDataSuppl)) {
cat(paste(" Replacing prescription row from input scenario.\n"))
plotData[row.names(plotDataSuppl), "PrescriptionRow"] = plotDataSuppl$PrescriptionRow
}
if(!("DemandRow" %in% names(plotData))) {
cat(paste(" Initializing demand row.\n"))
plotData$DemandRow = NA
dem_s = strsplit(as.character(demand$Species), split=",")
getDemandRow<-function(sp) {
for(i in 1:length(dem_s)) if(as.character(sp) %in% dem_s[[i]]) return(i)
return(NA)
}
for(i in 1:length(plotIDs)) {
plotData[i, "DemandRow"] = getDemandRow(plotData[i,"DominantSpecies"])
}
} else {
cat(paste(" Using demand row found in 'plotData'.\n"))
}
if(!is.null(plotDataSuppl)) {
cat(paste(" Replacing demand row from input scenario.\n"))
plotData[row.names(plotDataSuppl), "DemandRow"] = plotDataSuppl$DemandRow
}
if(!("finalPreviousStage" %in% names(plotData))) {
cat(paste(" Initializing final previous stage.\n"))
plotData$finalPreviousStage = NA
for(i in 1:length(plotIDs)) {
if(plotData[i, "Managed"]) plotData[i, "finalPreviousStage"] = prescription[plotData$PrescriptionRow[i], "finalPreviousStage"]
}
} else {
cat(paste(" Using final previous stage found in 'plotData'.\n"))
}
if(!is.null(plotDataSuppl)) {
cat(paste(" Replacing previous stage from input scenario.\n"))
plotData[row.names(plotDataSuppl), "finalPreviousStage"] = plotDataSuppl$finalPreviousStage
}
#Initial demand offset for all species
if(!is.null(demandOffset)) {
cat(" Setting initial demand offset to input.\n")
if(length(demandOffset)!= nrow(demand)) stop("Wrong demand offset length.")
TenYrDemandOffset = demandOffset
} else {
cat(" Setting initial demand offset to zero.\n")
TenYrDemandOffset = rep(0, nrow(demand))
}
manTurnOff = 0
for(i in 1:length(plotIDs)) {
if(is.na(plotData[i, "PrescriptionRow"]) || is.na(plotData[i, "DemandRow"])) {
if(plotData[i, "Managed"]) {
plotData[i, "Managed"] = FALSE
manTurnOff = manTurnOff+1
}
}
}
if(manTurnOff>0) cat(paste(" Management for", manTurnOff, "plots has been turned off because of missing prescription/demand.\n"))
# print(table(plotData$DominantSpecies, plotData$Managed))
# print(table(plotData$DominantSpecies, plotData$finalPreviousStage))
#Divide treeData between managed and non-managed
managedIDs = plotIDs[plotData$Managed]
nonmanagedIDs = plotIDs[!plotData$Managed]
cat(paste0(" ",length(managedIDs), " managed plots / ", length(nonmanagedIDs), " non-managed plots / ", nrow(treeData)," initial tree records.\n"))
prescription$plantingSpecies = NA
prescription$plantingDBH = 5
prescription$plantingHeight = 5
prescription$plantingDensity = 200
#Init output matrices
plotActions = matrix(NA, length(plotIDs), numSteps)
rownames(plotActions) = plotIDs
colnames(plotActions) = 1:numSteps
extractedDemand = matrix(0, nrow(demand)+1, numSteps)
rownames(extractedDemand) = c(demand$Name,"Total")
colnames(extractedDemand) = 1:numSteps
managedArea = matrix(0, nrow(demand)+1, numSteps)
rownames(managedArea) = c(demand$Name,"Total")
colnames(managedArea) = 1:numSteps
domspp = as.character(sort(as.numeric(unique(plotData$DominantSpecies))))
spp = as.character(sort(as.numeric(unique(c(treeData$Species, domspp))))) ## Add dominant species if not currently present
#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
#Copying area data
if(!("Area" %in% names(treeData))) treeData$Area = NA
if(!("Provincia" %in% names(treeData))) treeData$Provincia = NA
cat(" Copying area data to tree table.\n")
for(i in 1:nrow(plotData)) {
id = row.names(plotData)[i]
treeData$Provincia[treeData$ID==id] = plotData[id, "Province"]
treeData$Area[treeData$ID==id] = plotData[id, "Area"]
}
#Initial height
cat(" Calculating initial height.\n")
treeData$H = IFNheight(treeData, plotDataIFN, verbose = verbose)
#Initial volume
cat(" Calculating initial volume.\n")
initialVolume = IFNvolume(treeData, verbose = verbose)
volspp = tapply(initialVolume$VCC*treeData$Area, initialVolume$Species, sum, na.rm=T)
standingVolumeSpp[names(volspp), 1] = as.numeric(volspp)
volplot = tapply(initialVolume$VCC*treeData$Area, initialVolume$ID, sum, na.rm=T)
standingVolumePlot[names(volplot), 1] = as.numeric(volplot)
voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
standingVolumeDomSpp[names(voldomspp), 1] = as.numeric(voldomspp)
#Initial density
cat(" Calculating initial density.\n")
densspp = tapply(treeData$N*treeData$Area, treeData$Species, sum, na.rm=T)
standingDensitySpp[names(volspp), 1] = as.numeric(densspp)
densplot = tapply(treeData$N*treeData$Area, treeData$ID, sum, na.rm=T)
standingDensityPlot[names(densplot), 1] = as.numeric(densplot)
densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
standingDensityDomSpp[names(densdomspp), 1] = as.numeric(densdomspp)
#Initial biomass (Mg)
cat(" Calculating initial biomass.\n")
initialBiomass = IFNbiomass(treeData, as.CO2 = as.CO2, verbose = verbose)
rootspp = 0.001*tapply(initialBiomass$Roots*treeData$Area, initialBiomass$Species, sum, na.rm=T)
stemspp = 0.001*tapply(initialBiomass$Stem*treeData$Area, initialBiomass$Species, sum, na.rm=T)
branchspp = 0.001*tapply(initialBiomass$Branches*treeData$Area, initialBiomass$Species, sum, na.rm=T)
standingRootBiomassSpp[names(rootspp), 1] = as.numeric(rootspp)
standingStemBiomassSpp[names(stemspp), 1] = as.numeric(stemspp)
standingBranchBiomassSpp[names(branchspp), 1] = as.numeric(branchspp)
rootplot = 0.001*tapply(initialBiomass$Roots*treeData$Area, initialBiomass$ID, sum, na.rm=T)
stemplot = 0.001*tapply(initialBiomass$Stem*treeData$Area, initialBiomass$ID, sum, na.rm=T)
branchplot = 0.001*tapply(initialBiomass$Branches*treeData$Area, initialBiomass$ID, sum, na.rm=T)
standingRootBiomassPlot[names(rootplot), 1] = as.numeric(rootplot)
standingStemBiomassPlot[names(stemplot), 1] = as.numeric(stemplot)
standingBranchBiomassPlot[names(branchplot), 1] = as.numeric(branchplot)
rootdomspp = tapply(rootplot, plotData[names(rootplot), "DominantSpecies"], sum, na.rm=T)
stemdomspp = tapply(stemplot, plotData[names(stemplot), "DominantSpecies"], sum, na.rm=T)
branchdomspp = tapply(branchplot, plotData[names(branchplot), "DominantSpecies"], sum, na.rm=T)
standingRootBiomassDomSpp[names(rootdomspp), 1] = as.numeric(rootdomspp)
standingStemBiomassDomSpp[names(stemdomspp), 1] = as.numeric(stemdomspp)
standingBranchBiomassDomSpp[names(branchdomspp), 1] = as.numeric(branchdomspp)
cat("\n") ## End of data preparation
cat(paste("B. Simulations...\n"))
# Store initial state as sequence
if(sequence) {
treeDataSeq = data.frame(Step = rep(0, nrow(treeData)),
ID = treeData$ID,
Species = treeData$Species,
DBH = treeData$DBH,
N = treeData$N,
H = treeData$H,
stringsAsFactors = F)
}
if(dead.table) treeDataDeadSeq = treeDataSeq[numeric(0),]
if(cut.table) treeDataCutSeq = treeDataSeq[numeric(0),]
for(step in 1:numSteps) { # START OF STEP LOOP
cat(paste0(" *** Step ", step, "/", numSteps," ***\n"))
#B.00 calculate NspNmatrix
NspNmatrix = NULL
if(!is.null(nb)) {
cat(paste0(" Calculating Nsp/N matrix.\n"))
Nmatrix = tapply(treeData[,"N"], list(treeData[,"ID"], treeData[,"Species"]),
FUN = sum, na.rm=T)
NspNmatrix = sweep(Nmatrix,1, rowSums(Nmatrix, na.rm=T), "/")
rm(Nmatrix)
}
#B.0 Replace dynamic plot data
if(!is.null(plotDataDyn)) {
cat(paste0(" Drawing step [",(step + stepOffset), "] data from plotDataDyn.\n"))
plotDataDynStep = plotDataDyn[plotDataDyn$Step == (step + stepOffset),, drop=FALSE]
if(nrow(plotDataDynStep)==0) stop(paste0("Step data [", (step + stepOffset), "] not defined in plotDataDyn."))
row.names(plotDataDynStep) = plotDataDynStep$ID
}
#B. Process dynamics of non-managed plots
treeDataNonManagedNew = NULL
treeDataNonManagedDead = NULL
if(length(nonmanagedIDs)>0) {
cat(paste(" Processing", length(nonmanagedIDs),"non-managed plots:"))
treeDataNonManaged = treeData[treeData$ID %in% nonmanagedIDs, , drop = FALSE]
plotDataNonManaged = plotData[nonmanagedIDs,, drop = FALSE]
if(!is.null(plotDataDyn)) {
plotDataNonManaged$Temp = plotDataDynStep[nonmanagedIDs,"Temp"]
plotDataNonManaged$Prec = plotDataDynStep[nonmanagedIDs,"Prec"]
plotDataNonManaged$PET = plotDataDynStep[nonmanagedIDs,"PET"]
}
cat(paste0(" dynamics"))
#B.1 Calculate dynamics
dyn = IFNdynamics(treeDataNonManaged, plotDataNonManaged,
nb = nb, NspNmatrix = NspNmatrix,
numSteps = 1,
dead.table = TRUE, ...)
treeDataNonManagedNew = dyn$out[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
treeDataNonManagedDead = dyn$dead[, c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
#B.2 Add provincia/area
treeDataNonManagedNew$Provincia = NA
treeDataNonManagedNew$Area = NA
for(i in 1:nrow(plotDataNonManaged)) {
id = row.names(plotDataNonManaged)[i]
treeDataNonManagedNew$Provincia[treeDataNonManagedNew$ID==id] = plotDataNonManaged[id, "Province"]
treeDataNonManagedNew$Area[treeDataNonManagedNew$ID==id] = plotDataNonManaged[id, "Area"]
treeDataNonManagedDead$Provincia[treeDataNonManagedDead$ID==id] = plotDataNonManaged[id, "Province"]
treeDataNonManagedDead$Area[treeDataNonManagedDead$ID==id] = plotDataNonManaged[id, "Area"]
}
#B.3 Calculate standing/dead volume and density
cat(paste(", volume"))
nmvolume = IFNvolume(treeDataNonManagedNew)
volspp = tapply(nmvolume$VCC*treeDataNonManagedNew$Area, nmvolume$Species, sum, na.rm=T)
standingVolumeSpp[names(volspp), step+1] = standingVolumeSpp[names(volspp), step+1] + as.numeric(volspp)
volplot = tapply(nmvolume$VCC*treeDataNonManagedNew$Area, nmvolume$ID, sum, na.rm=T)
standingVolumePlot[names(volplot), step+1] = standingVolumePlot[names(volplot), step+1] + as.numeric(volplot)
voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
standingVolumeDomSpp[names(voldomspp), step+1] = standingVolumeDomSpp[names(voldomspp), step+1] + as.numeric(voldomspp)
nmdeadvolume = IFNvolume(treeDataNonManagedDead)
volspp = tapply(nmdeadvolume$VCC*treeDataNonManagedDead$Area, nmdeadvolume$Species, sum, na.rm=T)
deadVolumeSpp[names(volspp), step] = deadVolumeSpp[names(volspp), step] + as.numeric(volspp)
volplot = tapply(nmdeadvolume$VCC*treeDataNonManagedDead$Area, nmdeadvolume$ID, sum, na.rm=T)
deadVolumePlot[names(volplot), step] = deadVolumePlot[names(volplot), step] + as.numeric(volplot)
voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
deadVolumeDomSpp[names(voldomspp), step] = deadVolumeDomSpp[names(voldomspp), step] + as.numeric(voldomspp)
cat(paste(", density"))
densspp = tapply(treeDataNonManagedNew$N*treeDataNonManagedNew$Area, nmvolume$Species, sum, na.rm=T)
standingDensitySpp[names(densspp), step+1] = standingDensitySpp[names(densspp), step+1] + as.numeric(densspp)
densplot = tapply(treeDataNonManagedNew$N*treeDataNonManagedNew$Area, nmvolume$ID, sum, na.rm=T)
standingDensityPlot[names(densplot), step+1] = standingDensityPlot[names(densplot), step+1] + as.numeric(densplot)
densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
standingDensityDomSpp[names(densdomspp), step+1] = standingDensityDomSpp[names(densdomspp), step+1] + as.numeric(densdomspp)
densspp = tapply(treeDataNonManagedDead$N*treeDataNonManagedDead$Area, nmdeadvolume$Species, sum, na.rm=T)
deadDensitySpp[names(densspp), step] = deadDensitySpp[names(densspp), step] +as.numeric(densspp)
densplot = tapply(treeDataNonManagedDead$N*treeDataNonManagedDead$Area, nmdeadvolume$ID, sum, na.rm=T)
deadDensityPlot[names(densplot), step] = deadDensityPlot[names(densplot), step] + as.numeric(densplot)
densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
deadDensityDomSpp[names(densdomspp), step] = deadDensityDomSpp[names(densdomspp), step] + as.numeric(densdomspp)
#B.4 Calculate standing/dead biomass
cat(paste(", biomass"))
nmlivebiomass = IFNbiomass(treeDataNonManagedNew, as.CO2 = as.CO2)
biomassspp = 0.001*tapply(nmlivebiomass$Roots*treeDataNonManagedNew$Area, nmlivebiomass$Species, sum, na.rm=T)
standingRootBiomassSpp[names(biomassspp), step+1] = standingRootBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(nmlivebiomass$Roots*treeDataNonManagedNew$Area, nmlivebiomass$ID, sum, na.rm=T)
standingRootBiomassPlot[names(biomassplot), step+1] = standingRootBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
standingRootBiomassDomSpp[names(biomassdomspp), step+1] = standingRootBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(nmlivebiomass$Stem*treeDataNonManagedNew$Area, nmlivebiomass$Species, sum, na.rm=T)
standingStemBiomassSpp[names(biomassspp), step+1] = standingStemBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(nmlivebiomass$Stem*treeDataNonManagedNew$Area, nmlivebiomass$ID, sum, na.rm=T)
standingStemBiomassPlot[names(biomassplot), step+1] = standingStemBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
standingStemBiomassDomSpp[names(biomassdomspp), step+1] = standingStemBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(nmlivebiomass$Branches*treeDataNonManagedNew$Area, nmlivebiomass$Species, sum, na.rm=T)
standingBranchBiomassSpp[names(biomassspp), step+1] = standingBranchBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(nmlivebiomass$Branches*treeDataNonManagedNew$Area, nmlivebiomass$ID, sum, na.rm=T)
standingBranchBiomassPlot[names(biomassplot), step+1] = standingBranchBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
standingBranchBiomassDomSpp[names(biomassdomspp), step+1] = standingBranchBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)
nmdeadbiomass = IFNbiomass(treeDataNonManagedDead, as.CO2 = as.CO2)
biomassspp = 0.001*tapply(nmdeadbiomass$Roots*treeDataNonManagedDead$Area, nmdeadbiomass$Species, sum, na.rm=T)
deadRootBiomassSpp[names(biomassspp), step] = deadRootBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(nmdeadbiomass$Roots*treeDataNonManagedDead$Area, nmdeadbiomass$ID, sum, na.rm=T)
deadRootBiomassPlot[names(biomassplot), step] = deadRootBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
deadRootBiomassDomSpp[names(biomassdomspp), step] = deadRootBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(nmdeadbiomass$Stem*treeDataNonManagedDead$Area, nmdeadbiomass$Species, sum, na.rm=T)
deadStemBiomassSpp[names(biomassspp), step] = deadStemBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(nmdeadbiomass$Stem*treeDataNonManagedDead$Area, nmdeadbiomass$ID, sum, na.rm=T)
deadStemBiomassPlot[names(biomassplot), step] = deadStemBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
deadStemBiomassDomSpp[names(biomassdomspp), step] = deadStemBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(nmdeadbiomass$Branches*treeDataNonManagedDead$Area, nmdeadbiomass$Species, sum, na.rm=T)
deadBranchBiomassSpp[names(biomassspp), step] = deadBranchBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(nmdeadbiomass$Branches*treeDataNonManagedDead$Area, nmdeadbiomass$ID, sum, na.rm=T)
deadBranchBiomassPlot[names(biomassplot), step] = deadBranchBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
deadBranchBiomassDomSpp[names(biomassdomspp), step] = deadBranchBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
cat(".\n") ## End of non-managed plots (B)
}
# print(head(treeDataNonManagedNew))
# print(head(treeDataNonManagedDead))
# stop()
#C. Process dynamics of managed plots
cat(paste(" Processing", length(managedIDs),"managed plots:\n"))
plotDataManaged = plotData[managedIDs,]
treeDataManaged = treeData[treeData$ID %in% managedIDs, ]
treeDataManagedNew = NULL
treeDataManagedDead = NULL
for(r in 1:nrow(demand)) {
sel = (plotDataManaged$DemandRow == r)
if(sum(sel)>0) {
cat(paste0(" [",r,"/",nrow(demand),"] - ", sum(sel), " plots dominated by ", demand$Name[r], ":"))
#C.0 Subset plots dominaged by species
plotDataManagedSel = plotDataManaged[sel,]
treeDataManagedSel = treeDataManaged[treeDataManaged$ID %in% row.names(plotDataManagedSel), ]
prescriptionDataManagedSel = prescription[plotDataManagedSel$PrescriptionRow,-c(1,2)]
# Copy final previous stage (from previous step) from plot data
prescriptionDataManagedSel$finalPreviousStage = plotDataManaged$finalPreviousStage[sel]
row.names(prescriptionDataManagedSel) = row.names(plotDataManagedSel)
#C.2 Replace dynamic plot data
if(!is.null(plotDataDyn)) {
manSelIDs = row.names(plotDataManagedSel)
plotDataManagedSel$Temp = plotDataDynStep[manSelIDs,"Temp"]
plotDataManagedSel$Prec = plotDataDynStep[manSelIDs,"Prec"]
plotDataManagedSel$PET = plotDataDynStep[manSelIDs,"PET"]
}
#C.3 Calculate Hart-Becking
HB = hartBeckingIndex(treeDataManagedSel)
plotDataManagedSel$HB = HB[row.names(plotDataManagedSel)]
#C.4 Set demand
TenYrDemand = TenYrDemandOffset[r]+demand$AnnualDemand[r]*10
managedAreaSel = 0
Extracted = 0
treeDataOut = NULL #Stores the result of management
treeDataCut = NULL #Stores the result of management
treated = 0
if(!verbose) cat(paste0(" management"))
regularPlots = prescriptionDataManagedSel$type=="regular"
if(!verbose) cat(paste0(" [regular ", sum(regularPlots), " irregular ", sum(!regularPlots),"]"))
# if(!verbose) cat(paste0(" management [Dem. ", round(TenYrDemand),"m3"))
#C.5 Process first plots with regular models and previous stage > 0
finalCutPlots = (regularPlots) & (prescriptionDataManagedSel$finalPreviousStage>0)
if(sum(finalCutPlots)>0) {
if(!verbose) cat(paste0(", final cuts ", sum(finalCutPlots)," "))
plotDataFC = plotDataManagedSel[finalCutPlots, ]
prescriptionDataFC = prescriptionDataManagedSel[finalCutPlots,]
prescriptionDataFC$plantingSpecies = plotDataFC$DominantSpecies
treeDataFC = treeDataManagedSel[treeDataManagedSel$ID %in% row.names(plotDataManagedSel)[finalCutPlots], ]
if(nrow(treeDataFC)>0) {
treeDataFCPost = IFNmanagement(treeDataFC, prescriptionDataFC, verbose = verbose)
#store action
action = treeDataFCPost$action
plotActions[names(action), step] = action
#add area to tree data cut
treeDataFCPost$treeDataCut$Area = NA
for(i in 1:nrow(plotDataFC)) {
id = row.names(plotDataFC)[i]
treeDataFCPost$treeDataCut$Area[treeDataFCPost$treeDataCut$ID==id] = plotDataFC[id, "Area"]
}
#Calculate extracted volume
FCvolume = IFNvolume(treeDataFCPost$treeDataCut)
VCCplot = tapply(FCvolume$VCC*treeDataFCPost$treeDataCut$Area, FCvolume$ID, "sum", na.rm=TRUE) #Remove NA (TODO check volume equations)
#Update extractedVolumeSpp
volspp = tapply(FCvolume$VCC*treeDataFCPost$treeDataCut$Area, FCvolume$Species, sum, na.rm=T)
extractedVolumeSpp[names(volspp), step] = extractedVolumeSpp[names(volspp), step] + as.numeric(volspp)
volplot = tapply(FCvolume$VCC*treeDataFCPost$treeDataCut$Area, FCvolume$ID, sum, na.rm=T)
extractedVolumePlot[names(volplot), step] = extractedVolumePlot[names(volplot), step] + as.numeric(volplot)
voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
extractedVolumeDomSpp[names(voldomspp), step] = extractedVolumeDomSpp[names(voldomspp), step] + as.numeric(voldomspp)
#Update extractedDensitySpp
densspp= tapply(treeDataFCPost$treeDataCut$N*treeDataFCPost$treeDataCut$Area, FCvolume$Species, sum, na.rm=T)
extractedDensitySpp[names(densspp), step] = extractedDensitySpp[names(densspp), step] + as.numeric(densspp)
densplot= tapply(treeDataFCPost$treeDataCut$N*treeDataFCPost$treeDataCut$Area, FCvolume$ID, sum, na.rm=T)
extractedDensityPlot[names(densplot), step] = extractedDensityPlot[names(densplot), step] + as.numeric(densplot)
densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
extractedDensityDomSpp[names(densdomspp), step] = extractedDensityDomSpp[names(densdomspp), step] + as.numeric(densdomspp)
#Update extracted biomass
FCbiomass = IFNbiomass(treeDataFCPost$treeDataCut, as.CO2 = as.CO2)
biomassspp = 0.001*tapply(FCbiomass$Roots*treeDataFCPost$treeDataCut$Area, FCbiomass$Species, sum, na.rm=T)
extractedRootBiomassSpp[names(biomassspp), step] = extractedRootBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(FCbiomass$Roots*treeDataFCPost$treeDataCut$Area, FCbiomass$ID, sum, na.rm=T)
extractedRootBiomassPlot[names(biomassplot), step] = extractedRootBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
extractedRootBiomassDomSpp[names(biomassdomspp), step] = extractedRootBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(FCbiomass$Stem*treeDataFCPost$treeDataCut$Area, FCbiomass$Species, sum, na.rm=T)
extractedStemBiomassSpp[names(biomassspp), step] = extractedStemBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(FCbiomass$Stem*treeDataFCPost$treeDataCut$Area, FCbiomass$ID, sum, na.rm=T)
extractedStemBiomassPlot[names(biomassplot), step] = extractedStemBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
extractedStemBiomassDomSpp[names(biomassdomspp), step] = extractedStemBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(FCbiomass$Branches*treeDataFCPost$treeDataCut$Area, FCbiomass$Species, sum, na.rm=T)
extractedBranchBiomassSpp[names(biomassspp), step] = extractedBranchBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(FCbiomass$Branches*treeDataFCPost$treeDataCut$Area, FCbiomass$ID, sum, na.rm=T)
extractedBranchBiomassPlot[names(biomassplot), step] = extractedBranchBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
extractedBranchBiomassDomSpp[names(biomassdomspp), step] = extractedBranchBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
#Add to extracted
Extracted = Extracted + sum(VCCplot, na.rm=T)
treated = treated + sum(finalCutPlots)
#Add to managed area
managedAreaSel = managedAreaSel + sum(plotDataFC$Area)
#save final previous stage for next step
prescriptionDataManagedSel$finalPreviousStage[finalCutPlots] = treeDataFCPost$nextPrescription$finalPreviousStage
plotDataManagedSel$finalPreviousStage[finalCutPlots] = prescriptionDataManagedSel$finalPreviousStage[finalCutPlots]
#Store result
if(!is.null(treeDataOut)) {
treeDataOut = rbind(treeDataOut, treeDataFCPost$treeDataOut)
} else {
treeDataOut = treeDataFCPost$treeDataOut
}
if(!is.null(treeDataCut)) {
treeDataCut = rbind(treeDataCut, treeDataFCPost$treeDataCut)
} else {
treeDataCut = treeDataFCPost$treeDataCut
}
}
}
if(sum(!finalCutPlots)>0) {
plotDataTh = plotDataManagedSel[!finalCutPlots, ]
prescriptionDataTh = prescriptionDataManagedSel[!finalCutPlots,]
treeDataTh = treeDataManagedSel[treeDataManagedSel$ID %in% row.names(plotDataManagedSel)[!finalCutPlots], ]
plotActions[row.names(plotDataTh), step] = "none"
#If demand has not been fulfilled continue
if((Extracted<TenYrDemand)) {
#Find difference between HB and thinningHB
diffHB = plotDataTh$HB-prescriptionDataTh$thinningHB
o = order(diffHB, decreasing=FALSE) #Sort from largest difference to more negative difference
i = 1
while((Extracted<TenYrDemand) && (i<=length(o))) {
plotDataOne = plotDataTh[o[i], ]
prescriptionDataOne = prescriptionDataTh[o[i],]
prescriptionDataOne$plantingSpecies = plotDataOne$DominantSpecies
treeDataOne = treeDataTh[treeDataTh$ID %in% row.names(plotDataOne), ]
#update treeDataTh
treeDataTh = treeDataTh[!(treeDataTh$ID %in% row.names(plotDataOne)), ]
if(nrow(treeDataOne)>0) {
treeDataOnePost = IFNmanagement(treeDataOne, prescriptionDataOne, verbose = verbose)
#save next prescription
prescriptionDataManagedSel$finalPreviousStage[!finalCutPlots][o[i]] = treeDataOnePost$nextPrescription$finalPreviousStage
plotDataManagedSel$finalPreviousStage[!finalCutPlots][o[i]] = prescriptionDataManagedSel$finalPreviousStage[!finalCutPlots][o[i]]
#store actions
action = treeDataOnePost$action
plotActions[names(action), step] = action
#Calculate extracted volume
if(nrow(treeDataOnePost$treeDataCut)>0) {
FCvolume = IFNvolume(treeDataOnePost$treeDataCut)
VCCarea = sum(FCvolume$VCC, na.rm=T)*as.numeric(plotDataOne$Area)
if(is.na(VCCarea)) {
stop("VCCarea = NA")
}
#Update extracted volume and density
volspp = tapply(FCvolume$VCC, FCvolume$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedVolumeSpp[names(volspp), step] = extractedVolumeSpp[names(volspp), step] + as.numeric(volspp)
volplot = tapply(FCvolume$VCC, FCvolume$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedVolumePlot[names(volplot), step] = extractedVolumePlot[names(volplot), step] + as.numeric(volplot)
voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
extractedVolumeDomSpp[names(voldomspp), step] = extractedVolumeDomSpp[names(voldomspp), step] + as.numeric(voldomspp)
densspp= tapply(treeDataOnePost$treeDataCut$N, FCvolume$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedDensitySpp[names(densspp), step] = extractedDensitySpp[names(densspp), step] + as.numeric(densspp)
densplot= tapply(treeDataOnePost$treeDataCut$N, FCvolume$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedDensityPlot[names(densplot), step] = extractedDensityPlot[names(densplot), step] + as.numeric(densplot)
densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
extractedDensityDomSpp[names(densdomspp), step] = extractedDensityDomSpp[names(densdomspp), step] + as.numeric(densdomspp)
#Update extracted biomass
FCbiomass = IFNbiomass(treeDataOnePost$treeDataCut, as.CO2 = as.CO2)
biomassspp = 0.001*tapply(FCbiomass$Roots, FCbiomass$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedRootBiomassSpp[names(biomassspp), step] = extractedRootBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(FCbiomass$Roots, FCbiomass$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedRootBiomassPlot[names(biomassplot), step] = extractedRootBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
extractedRootBiomassDomSpp[names(biomassdomspp), step] = extractedRootBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(FCbiomass$Stem, FCbiomass$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedStemBiomassSpp[names(biomassspp), step] = extractedStemBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(FCbiomass$Stem, FCbiomass$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedStemBiomassPlot[names(biomassplot), step] = extractedStemBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
extractedStemBiomassDomSpp[names(biomassdomspp), step] = extractedStemBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(FCbiomass$Branches, FCbiomass$Species, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedBranchBiomassSpp[names(biomassspp), step] = extractedBranchBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(FCbiomass$Branches, FCbiomass$ID, sum, na.rm=T)*as.numeric(plotDataOne$Area)
extractedBranchBiomassPlot[names(biomassplot), step] = extractedBranchBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
extractedBranchBiomassDomSpp[names(biomassdomspp), step] = extractedBranchBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
#Add to extracted
Extracted = Extracted + VCCarea
treated = treated+1
#Add to managed area
managedAreaSel = managedAreaSel + as.numeric(plotDataOne$Area)
}
#Store result
if(!is.null(treeDataOut)) {
treeDataOut = rbind(treeDataOut, treeDataOnePost$treeDataOut)
} else {
treeDataOut = treeDataOnePost$treeDataOut
}
if(!is.null(treeDataCut)) {
treeDataCut = rbind(treeDataCut, treeDataOnePost$treeDataCut)
} else {
treeDataCut = treeDataOnePost$treeDataCut
}
}
# print(treeDataOnePost$nextPrescription)
#Increase plot counter
i = i+1
}
}
#Add the remaining (non-intervention plots) to treeDataOut
if(nrow(treeDataTh)>0) {
if(!is.null(treeDataOut)) {
# print(head(treeDataTh))
# print(head(treeDataOut))
treeDataOut = rbind(treeDataOut, treeDataTh)
} else {
treeDataOut = treeDataTh
}
}
}
#Copy final previous stage to plotDataManaged
plotDataManaged$finalPreviousStage[sel] = plotDataManagedSel$finalPreviousStage # copy to plotDataManaged (for next step)
#Add to cut table (for output)
if(cut.table) {
# print("merging cut")
if(!is.null(treeDataCut) && (nrow(treeDataCut)>0)) {
treeDataCutStep = treeDataCut
treeDataCutStep$Step = step
if(is.null(treeDataCutSeq)) {
treeDataCutSeq = treeDataCutStep[,c("Step","ID","Species", "DBH","N", "H")]
} else {
# print(names(treeDataCut))
# print(names(treeDataCutSeq))
treeDataCutSeq = rbind(treeDataCutSeq,treeDataCutStep[, names(treeDataCutSeq), drop=FALSE])
}
}
}
# Uodate demand offset (for next step)
TenYrDemandOffset[r] = TenYrDemand - Extracted
# if(!verbose) {
# cat(paste0(", Extr. ", round(Extracted),"m3, Diff. ", round(TenYrDemandOffset[r]),"m3, Treated ", treated,"]"))
# } else {
# cat(paste0("\n Management [Dem. ", round(TenYrDemand),"m3"))
# cat(paste0(", Extr. ", round(Extracted),"m3, Diff. ", round(TenYrDemandOffset[r]),"m3, Treated ", treated,"]"))
# }
extractedDemand[r,step] = Extracted
extractedDemand[nrow(demand)+1,step] = extractedDemand[nrow(demand)+1,step]+Extracted
managedArea[r,step] = managedAreaSel
managedArea[nrow(demand)+1,step] = managedArea[nrow(demand)+1,step] + managedAreaSel
# return(list(treeDataOut, plotDataManagedSel))
#C.6 Run forest dynamics for all managed plots
cat(paste(", dynamics"))
treeDataManagedSelNew = NULL
treeDataManagedSelDead = NULL
idsIrregular = row.names(prescriptionDataManagedSel)[prescriptionDataManagedSel$type=="irregular"]
tdoirr= treeDataOut[treeDataOut$ID %in% idsIrregular,, drop=FALSE]
if(nrow(tdoirr)>0) {
tdoirr$OIFFIN = as.character(1:nrow(tdoirr))
##Add cut records with OIFFIN = 0 to calculate BAL_ext
if(!is.null(treeDataCut)) {
tdci = treeDataCut[treeDataCut$ID %in% idsIrregular,, drop=FALSE]
if(nrow(tdci)>0) {
tdci$OIFFIN = "0"
tdoirr = rbind(tdoirr,tdci)
}
}
dyn = IFNdynamics(tdoirr, plotDataManagedSel[idsIrregular, , drop=FALSE],
nb = nb, NspNmatrix = NspNmatrix,
numSteps = 1,
ingrowth = T,
dead.table = T,
...)
treeDataManagedSelNew = dyn$out[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
treeDataManagedSelDead = dyn$dead[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
}
idsRegular = row.names(prescriptionDataManagedSel)[prescriptionDataManagedSel$type!="irregular"]
tdor= treeDataOut[treeDataOut$ID %in% idsRegular,, drop=FALSE]
if(nrow(tdor)>0) {
tdor$OIFFIN = as.character(1:nrow(tdor))
##Add cut records with OIFFIN = 0 to calculate BAL_ext
if(!is.null(treeDataCut)) {
tdci = treeDataCut[treeDataCut$ID %in% idsRegular,, drop=FALSE]
if(nrow(tdci)>0) {
tdci$OIFFIN = "0"
tdor = rbind(tdor,tdci)
}
}
dyn = IFNdynamics(tdor, plotDataManagedSel[idsRegular, , drop=FALSE],
nb = nb, NspNmatrix = NspNmatrix,
numSteps = 1,
ingrowth = F,
dead.table = T,
...)
if(!is.null(treeDataManagedSelNew)) treeDataManagedSelNew = rbind(treeDataManagedSelNew, dyn$out[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE])
else treeDataManagedSelNew = dyn$out[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
if(!is.null(treeDataManagedSelDead)) treeDataManagedSelDead = rbind(treeDataManagedSelDead, dyn$dead[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE])
else treeDataManagedSelDead = dyn$dead[,c("ID", "Species", "Name", "DBH", "N", "H"), drop = FALSE]
}
if(nrow(treeDataManagedSelNew)>0) {
treeDataManagedSelNew$Provincia = plotDataManagedSel[treeDataManagedSelNew$ID, "Province"]
treeDataManagedSelNew$Area = plotDataManagedSel[treeDataManagedSelNew$ID, "Area"]
}
if(nrow(treeDataManagedSelDead)>0) {
treeDataManagedSelDead$Provincia = plotDataManagedSel[treeDataManagedSelDead$ID, "Province"]
treeDataManagedSelDead$Area = plotDataManagedSel[treeDataManagedSelDead$ID, "Area"]
}
#C.7 Calculate standing/dead volume and density
cat(paste(", volume"))
mvolume = IFNvolume(treeDataManagedSelNew)
volspp = tapply(mvolume$VCC*treeDataManagedSelNew$Area, mvolume$Species, sum, na.rm=T)
standingVolumeSpp[names(volspp), step+1] = standingVolumeSpp[names(volspp), step+1] + as.numeric(volspp)
volplot = tapply(mvolume$VCC*treeDataManagedSelNew$Area, mvolume$ID, sum, na.rm=T)
standingVolumePlot[names(volplot), step+1] = standingVolumePlot[names(volplot), step+1] + as.numeric(volplot)
voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
standingVolumeDomSpp[names(voldomspp),step+ 1] = standingVolumeDomSpp[names(voldomspp), step+1] + as.numeric(voldomspp)
mdeadvolume = IFNvolume(treeDataManagedSelDead)
volspp = tapply(mdeadvolume$VCC*treeDataManagedSelDead$Area, mdeadvolume$Species, sum, na.rm=T)
deadVolumeSpp[names(volspp), step] = deadVolumeSpp[names(volspp), step] + as.numeric(volspp)
volplot = tapply(mdeadvolume$VCC*treeDataManagedSelDead$Area, mdeadvolume$ID, sum, na.rm=T)
deadVolumePlot[names(volplot), step] = deadVolumePlot[names(volplot), step] + as.numeric(volplot)
voldomspp = tapply(volplot, plotData[names(volplot), "DominantSpecies"], sum, na.rm=T)
deadVolumeDomSpp[names(voldomspp), step] = deadVolumeDomSpp[names(voldomspp), step] + as.numeric(voldomspp)
cat(paste(", density"))
densspp = tapply(treeDataManagedSelNew$N*treeDataManagedSelNew$Area, mvolume$Species, sum, na.rm=T)
standingDensitySpp[names(densspp), step+1] = standingDensitySpp[names(densspp), step+1] + as.numeric(densspp)
densplot =tapply(treeDataManagedSelNew$N*treeDataManagedSelNew$Area, mvolume$ID, sum, na.rm=T)
standingDensityPlot[names(densplot), step+1] = standingVolumePlot[names(densplot), step+1] + as.numeric(densplot)
densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
standingDensityDomSpp[names(densdomspp),step+ 1] = standingDensityDomSpp[names(densdomspp), step+1] + as.numeric(densdomspp)
densspp = tapply(treeDataManagedSelDead$N*treeDataManagedSelDead$Area, mdeadvolume$Species, sum, na.rm=T)
deadDensitySpp[names(densspp), step] = as.numeric(densspp)
densplot = tapply(treeDataManagedSelDead$N*treeDataManagedSelDead$Area, mdeadvolume$ID, sum, na.rm=T)
deadDensityPlot[names(densplot), step] = as.numeric(densplot)
densdomspp = tapply(densplot, plotData[names(densplot), "DominantSpecies"], sum, na.rm=T)
deadDensityDomSpp[names(densdomspp), step] = deadDensityDomSpp[names(densdomspp), step] + as.numeric(densdomspp)
#C.8 Calculate standing/dead biomass
cat(paste(", biomass"))
mlivebiomass = IFNbiomass(treeDataManagedSelNew, as.CO2 = as.CO2)
biomassspp = 0.001*tapply(mlivebiomass$Roots*treeDataManagedSelNew$Area, mlivebiomass$Species, sum, na.rm=T)
standingRootBiomassSpp[names(biomassspp), step+1] = standingRootBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(mlivebiomass$Roots*treeDataManagedSelNew$Area, mlivebiomass$ID, sum, na.rm=T)
standingRootBiomassPlot[names(biomassplot), step+1] = standingRootBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
standingRootBiomassDomSpp[names(biomassdomspp), step+1] = standingRootBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(mlivebiomass$Stem*treeDataManagedSelNew$Area, mlivebiomass$Species, sum, na.rm=T)
standingStemBiomassSpp[names(biomassspp), step+1] = standingStemBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(mlivebiomass$Stem*treeDataManagedSelNew$Area, mlivebiomass$ID, sum, na.rm=T)
standingStemBiomassPlot[names(biomassplot), step+1] = standingStemBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
standingStemBiomassDomSpp[names(biomassdomspp), step+1] = standingStemBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(mlivebiomass$Branches*treeDataManagedSelNew$Area, mlivebiomass$Species, sum, na.rm=T)
standingBranchBiomassSpp[names(biomassspp), step+1] = standingBranchBiomassSpp[names(biomassspp), step+1] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(mlivebiomass$Branches*treeDataManagedSelNew$Area, mlivebiomass$ID, sum, na.rm=T)
standingBranchBiomassPlot[names(biomassplot), step+1] = standingBranchBiomassPlot[names(biomassplot), step+1] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
standingBranchBiomassDomSpp[names(biomassdomspp), step+1] = standingBranchBiomassDomSpp[names(biomassdomspp), step+1] + as.numeric(biomassdomspp)
mdeadbiomass = IFNbiomass(treeDataManagedSelDead, as.CO2 = as.CO2)
biomassspp = 0.001*tapply(mdeadbiomass$Roots*treeDataManagedSelDead$Area, mdeadbiomass$Species, sum, na.rm=T)
deadRootBiomassSpp[names(biomassspp), step] = deadRootBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(mdeadbiomass$Roots*treeDataManagedSelDead$Area, mdeadbiomass$ID, sum, na.rm=T)
deadRootBiomassPlot[names(biomassplot), step] = deadRootBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
deadRootBiomassDomSpp[names(biomassdomspp), step] = deadRootBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(mdeadbiomass$Stem*treeDataManagedSelDead$Area, mdeadbiomass$Species, sum, na.rm=T)
deadStemBiomassSpp[names(biomassspp), step] = deadStemBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(mdeadbiomass$Stem*treeDataManagedSelDead$Area, mdeadbiomass$ID, sum, na.rm=T)
deadStemBiomassPlot[names(biomassplot), step] = deadStemBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
deadStemBiomassDomSpp[names(biomassdomspp), step] = deadStemBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
biomassspp = 0.001*tapply(mdeadbiomass$Branches*treeDataManagedSelDead$Area, mdeadbiomass$Species, sum, na.rm=T)
deadBranchBiomassSpp[names(biomassspp), step] = deadBranchBiomassSpp[names(biomassspp), step] + as.numeric(biomassspp)
biomassplot = 0.001*tapply(mdeadbiomass$Branches*treeDataManagedSelDead$Area, mdeadbiomass$ID, sum, na.rm=T)
deadBranchBiomassPlot[names(biomassplot), step] = deadBranchBiomassPlot[names(biomassplot), step] + as.numeric(biomassplot)
biomassdomspp = tapply(biomassplot, plotData[names(biomassplot), "DominantSpecies"], sum, na.rm=T)
deadBranchBiomassDomSpp[names(biomassdomspp), step] = deadBranchBiomassDomSpp[names(biomassdomspp), step] + as.numeric(biomassdomspp)
cat(paste(".\n"))
if(nrow(treeDataManagedSelDead)>0) {
if(is.null(treeDataManagedDead)) {
treeDataManagedDead = treeDataManagedSelDead
} else {
treeDataManagedDead = rbind(treeDataManagedDead,treeDataManagedSelDead)
}
}
if(nrow(treeDataManagedSel)>0) {
if(is.null(treeDataManagedNew)){
treeDataManagedNew = treeDataManagedSelNew
} else {
treeDataManagedNew = rbind(treeDataManagedNew, treeDataManagedSelNew)
}
}
} else {
cat(paste0(" [",r,"/",nrow(demand),"] - ", sum(sel), " plots dominated by ", demand$Name[r], "\n"))
}
}
#Copy final previous stage to plotData
plotData[managedIDs, "finalPreviousStage"] = plotDataManaged$finalPreviousStage
cat("\n") ## End of managed plots (C)
# print(head(treeDataManagedNew))
# print(head(treeDataManagedDead))
## Assemble tree data for next step
treeData = rbind(treeDataNonManagedNew, treeDataManagedNew)
treeDataDead = rbind(treeDataNonManagedDead, treeDataManagedDead)
if(sequence) {
treeDataStep = treeData
treeDataStep$Step = step
treeDataSeq = rbind(treeDataSeq, treeDataStep[,names(treeDataSeq)])
}
if(dead.table) {
treeDataDead$Step = step
if(is.null(treeDataDeadSeq)) {
treeDataDeadSeq = treeDataDead[,c("Step","ID","Species", "DBH","N", "H"), drop=FALSE]
} else {
treeDataDeadSeq = rbind(treeDataDeadSeq,treeDataDead[,names(treeDataDeadSeq), drop=FALSE])
}
}
} # END OF STEP LOOP
# C. Assemble output
cat(paste("C. Post-processing...\n"))
cat(" Calculating coverage of demand.\n")
extractedDemand = as.data.frame(extractedDemand)
extractedDemand$AnnualMean = rowMeans(extractedDemand)/10
extractedDemand$Demand = c(demand$AnnualDemand, sum(demand$AnnualDemand))
extractedDemand$Coverage = 100*extractedDemand$AnnualMean/extractedDemand$Demand
cat(" Calculating managed area.\n")
managedArea = as.data.frame(managedArea)
md = tapply(plotData$Area[plotData$Managed], plotData$DemandRow[plotData$Managed], sum, na.rm =T)
managedArea$ManageableArea = 0
managedArea$ManageableArea[as.numeric(names(md))] = md
managedArea$ManageableArea[nrow(managedArea)] = sum(md)
totalArea = sum(plotData$Area)
#Calculate variation
# variation_t = standing_t - standing_(t-1)
cat(" Calculating variation of stands between steps.\n")
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)
cat(" Calculating stand growth between steps.\n")
# 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
cat(" Calculating statistics per hectare.\n")
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))
cat(" Assembling output.\n")
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)
plotDataSuppl = plotData[,c("Area", "DominantSpecies", "DemandRow", "PrescriptionRow", "finalPreviousStage")]
res = list(plots = nrow(plotData),
steps = numSteps,
stepOffset = stepOffset,
totalArea = totalArea,
plotDataSuppl = plotDataSuppl,
plotActions = plotActions,
extractedDemand=extractedDemand,
demandOffset = TenYrDemandOffset,
managedArea = managedArea,
volume=volume,
density=density,
rootBiomass = rootBiomass,
stemBiomass = stemBiomass,
branchBiomass = branchBiomass,
treeDataNew = treeData)
if(sequence) {
if((!is.null(treeDataSeq)) && (nrow(treeDataSeq)>0)) row.names(treeDataSeq) = 1:nrow(treeDataSeq)
res[["treeDataSequence"]] = treeDataSeq
}
if(dead.table) {
if((!is.null(treeDataDeadSeq)) && (nrow(treeDataDeadSeq)>0)) row.names(treeDataDeadSeq) = 1:nrow(treeDataDeadSeq)
res[["treeDataDead"]] = treeDataDeadSeq
}
if(cut.table) {
if((!is.null(treeDataCutSeq)) && (nrow(treeDataCutSeq)>0)) row.names(treeDataCutSeq) = 1:nrow(treeDataCutSeq)
res[["treeDataCut"]] = treeDataCutSeq
}
class(res) <- "IFNscenario"
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.