#' Scenario biomass and emissions
#'
#' Function \code{scenarioLiveCarbonStock} estimates the amount of living biomass per time step.
#' Functions \code{scenarioDeadWood} and \code{scenarioProducts} estimate the dead and product biomass generated in each time step.
#' Finally, function \code{scenarioEmissions} estimates the carbon source and sink effects for each time step of the simulated scenario.
#'
#' @param scenario The result of function \code{\link{IFNscenario}}.
#' @param productDestination Data frame of product destination by species and diameter class
#' @param variable Either 'biomass' (default) or 'volume'
#' @param ... Additional parameters to \code{\link{IFNproducts}}.
#'
#' @name scenarioEmissions
#'
#' @return Functions \code{scenarioLiveCarbonStock} and \code{scenarioDeadWood} return a list of two data frames with
#' biomass (dry weight or CO2, depending on \code{as.CO2}). The first data data frame contains biomass values in Mg for the target area, whereas the second contains
#' biomass values in Mg/ha.
#'
#' Function \code{scenarioProducts} can estimate both biomass or volumes (by calling function \code{\link{IFNproducts}}). Hence, and compared to the previous functions, the two data frames can correspond to volumes (m3 and m3/ha, respectively) instead of biomass values.
#'
#' Function \code{scenarioEmissions} returns a list with the following elements:
#' \itemize{
#' \item{\code{initialDead}: Initial values of dead matter components (in Mg/ha).}
#' \item{\code{initialProducts}: Initial values of products (in Mg/ha).}
#' \item{\code{liveWoodChange}: Changes in live biomass components (in Mg/ha) per time step.}
#' \item{\code{deadWoodGeneration}: Generation of dead biomass components (in Mg/ha) per time step.}
#' \item{\code{productWoodGeneration}: Generation of products (in Mg/ha) per time step.}
#' \item{\code{carbonBalance}: List of two frames with carbon balance details (including net emissions) in Mg/yr or Mg/ha/yr, respectively.}
#' \item{\code{finalDead}: Final values of dead matter components (in Mg/ha).}
#' \item{\code{finalProducts}: Final values of products (in Mg/ha).}
#' }
#' @export
#'
scenarioProducts<-function(scenario, productDestination, variable = "biomass",
...) {
variable = match.arg(variable, c("biomass", "volume"))
products = names(productDestination)[-c(1:4)]
if(variable=="biomass") destination = c("stumps","coarse_slash","medium_slash", "fine_slash", products)
else destination = c("slash", products)
treeDataCut = scenario$treeDataCut
plotArea = scenario$plotDataSuppl$Area
totalArea = sum(plotArea, na.rm=T)
names(plotArea) = row.names(scenario$plotDataSuppl)
y = IFNproducts(x= treeDataCut,
productDestination = productDestination,
variable = variable,
...)
y$Area = plotArea[y$ID]
if(variable=="biomass") {
for(var in destination) y[[var]] = (y[[var]]/1000)*(y$Area) # Transform to Mg
} else {
for(var in destination) y[[var]] = (y[[var]])*(y$Area) # Transform to m3
}
res = y %>% group_by(Step) %>%
summarize_at(destination, sum, na.rm=T)
res = as.data.frame(res)
row.names(res) = res$Step
res = res[,-1]
res_ha = res
res_ha[,destination] = res_ha[,destination]/totalArea
if(variable=="biomass") return(list(Mg = res, Mg_ha = res_ha))
return(list(m3 = res, m3_ha = res_ha))
}
#' @rdname scenarioEmissions
#' @param as.CO2 Flag to indicate output as kg of CO2 / ha instead of kg of dry weight / ha. Percentage of carbon per dry weight biomass by species are
#' taken from Montero et al. (2005) (in turn, from Ibáñez et al. 2002).
#' @param area Either 'Atlantic' or 'Mediterranean' to specify allometric equations specific to the area (for Pinus pinaster)
scenarioDeadWood<-function(scenario, as.CO2 = FALSE, area = NA) {
treeDataDead = scenario$treeDataDead
plotArea = scenario$plotDataSuppl$Area
names(plotArea) = row.names(scenario$plotDataSuppl)
totalArea = sum(plotArea)
destination = c("roots_dead", "coarse_dead", "medium_dead", "fine_dead")
translateToDestiny<-function(Bs) {
res = data.frame(matrix(0, nrow=nrow(Bs), ncol = length(destination)))
names(res)<- destination
row.names(res)<-1:nrow(res)
#Set missing values to zero
Bs$Roots[is.na(Bs$Roots)] = 0
Bs$Needles[is.na(Bs$Needles)] = 0
Bs$Leaves[is.na(Bs$Leaves)] = 0
Bs$Bark[is.na(Bs$Bark)] = 0
Bs$Branches[is.na(Bs$Branches)] = 0
Bs$Stem[is.na(Bs$Stem)] = 0
#assign
res$roots_dead =Bs$Roots
res$fine_dead = Bs$Leaves + Bs$Needles + Bs$Bark
res$medium_dead = Bs$Branches
res$coarse_dead = Bs$Stem
return(res)
}
isSteps = "Step" %in% names(treeDataDead)
if(isSteps) {
steps = unique(treeDataDead$Step)
numSteps = length(steps)
res = data.frame(matrix(0, nrow = numSteps, ncol = length(destination)))
names(res) =destination
row.names(res) = steps
for(i in 1:length(steps)) {
Bs = IFNbiomass(treeDataDead[treeDataDead$Step==steps[i],],
as.CO2 = as.CO2, area = area)
Bs$plotArea = plotArea[Bs$ID]
Bs$Roots = Bs$Roots*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Stem = Bs$Stem*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Branches = Bs$Branches*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Leaves = Bs$Leaves*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Needles = Bs$Needles*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Bark = Bs$Bark*Bs$plotArea/1000 # Translate from kg/ha to Mg
prodi = translateToDestiny(Bs)
res[i,destination] = colSums(prodi, na.rm=T)
}
return(list(Mg = res, Mg_ha = res/totalArea))
} else {
Bs = IFNbiomass(treeDataCut, DBHclasses = DBHclasslimits,
as.CO2 = as.CO2, area = area)
Bs$plotArea = plotArea[Bs$ID]
Bs$Roots = Bs$Roots*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Stem = Bs$Stem*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Branches = Bs$Branches*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Leaves = Bs$Leaves*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Needles = Bs$Needles*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Bark = Bs$Bark*Bs$plotArea/1000 # Translate from kg/ha to Mg
prod = translateToDestiny(Bs)
res = colSums(prod, na.rm=T)
return(list(Mg = res, Mg_ha = res/totalArea))
}
}
#' @rdname scenarioEmissions
scenarioLiveCarbonStock<-function(scenario, as.CO2 = FALSE, area = NA) {
treeDataSequence = scenario$treeDataSequence
plotArea = scenario$plotDataSuppl$Area
names(plotArea) = row.names(scenario$plotDataSuppl)
totalArea = sum(plotArea)
destination = c("roots_live", "coarse_live", "medium_live", "fine_live")
translateToDestiny<-function(Bs) {
res = data.frame(matrix(0, nrow=nrow(Bs), ncol = length(destination)))
names(res)<- destination
row.names(res)<-1:nrow(res)
#Set missing values to zero
Bs$Roots[is.na(Bs$Roots)] = 0
Bs$Needles[is.na(Bs$Needles)] = 0
Bs$Leaves[is.na(Bs$Leaves)] = 0
Bs$Bark[is.na(Bs$Bark)] = 0
Bs$Branches[is.na(Bs$Branches)] = 0
Bs$Stem[is.na(Bs$Stem)] = 0
#assign
res$roots_live =Bs$Roots
res$fine_live = Bs$Leaves + Bs$Needles + Bs$Bark
res$medium_live = Bs$Branches
res$coarse_live = Bs$Stem
return(res)
}
isSteps = "Step" %in% names(treeDataSequence)
steps = unique(treeDataSequence$Step)
numSteps = length(steps)
res = data.frame(matrix(0, nrow = numSteps, ncol = length(destination)))
names(res) =destination
row.names(res) = steps
for(i in 1:length(steps)) {
Bs = IFNbiomass(treeDataSequence[treeDataSequence$Step==steps[i],],
as.CO2 = as.CO2, area = area)
Bs$plotArea = plotArea[Bs$ID]
Bs$Roots = Bs$Roots*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Stem = Bs$Stem*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Branches = Bs$Branches*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Leaves = Bs$Leaves*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Needles = Bs$Needles*Bs$plotArea/1000 # Translate from kg/ha to Mg
Bs$Bark = Bs$Bark*Bs$plotArea/1000 # Translate from kg/ha to Mg
prodi = translateToDestiny(Bs)
res[i,destination] = colSums(prodi, na.rm=T)
}
return(list(Mg = res, Mg_ha = res/totalArea))
}
#' @rdname scenarioEmissions
#' @param initialDead Initial values of dead matter components (in Mg/ha)
#' @param initialProducts Initial values of product stocks (in Mg/ha)
#' @param deadHalfLives Half-live (in yrs) of dead matter components.
#' @param productHalfLives Half-live (in yrs) of products.
scenarioEmissions<-function(scenario, productDestination,
initialDead = c("roots_dead" = 5, "coarse_dead" = 3, "medium_dead" = 3, "fine_dead" = 1),
initialProducts = c("trituration" = 0.5, "poles" = 0.1, "boards" = 0.5, "construction" = 0.5, "firewood" = 0.1, "veneer" = 0.1),
deadHalfLives = c("roots_dead" = 20, "coarse_dead" = 10, "medium_dead" = 5, "fine_dead" = 2),
productHalfLives = c("trituration" = 2, "poles" = 25, "boards" = 15, "construction" = 35, "firewood" = 2, "veneer" = 25),
as.CO2 = FALSE, area = NA) {
products = names(productDestination)[-c(1:4)]
slash = c("stumps", "coarse_slash", "medium_slash", "fine_slash")
dead = c("roots_dead", "coarse_dead", "medium_dead", "fine_dead")
totalArea = scenario$totalArea
steps = scenario$steps
cat(paste0("Calculating live carbon stocks...\n"))
scenLiveStock = scenarioLiveCarbonStock(scenario = scenario,
as.CO2 = as.CO2, area = area)
scenStockChange = scenLiveStock$Mg_ha[2:(steps+1),] - scenLiveStock$Mg_ha[1:steps,]
cat(paste0("Calculating dead wood...\n"))
scenDeadWood = scenarioDeadWood(scenario = scenario,
as.CO2 = as.CO2, area = area)
cat(paste0("Calculating products...\n"))
scenProd = scenarioProducts(scenario = scenario,
productDestination = productDestination,
as.CO2 = as.CO2, area = area)
slashProd = scenProd$Mg_ha[,slash]
deadWoodGeneration = scenDeadWood$Mg_ha[,dead]
for(i in 1:nrow(slashProd)) {
rn = row.names(slashProd)[i]
if(rn %in% row.names(deadWoodGeneration)) {
deadWoodGeneration[rn,] = deadWoodGeneration[rn,] + slashProd[rn,]
}
}
productWoodGeneration = data.frame(matrix(0, nrow=nrow(deadWoodGeneration), ncol=length(products)))
names(productWoodGeneration) = products
for(i in 1:nrow(scenProd$Mg_ha)) {
rn = row.names(scenProd$Mg_ha)[i]
if(rn %in% row.names(productWoodGeneration)) {
productWoodGeneration[rn,] = scenProd$Mg_ha[rn,products]
}
}
deadComp = initialDead
prodComp = initialProducts
deadK = log(2)/deadHalfLives
productK = log(2)/productHalfLives
decompositionEmissions = rep(0, steps*10)
productEndlifeEmissions = rep(0, steps*10)
livestockVariation = rep(0, steps*10)
deadstockGeneration = rep(0, steps*10)
productstockGeneration = rep(0, steps*10)
NPPSink = rep(0, steps*10)
yearvec = 1:(steps*10)
stepvec = rep(0, steps*10)
for(step in 1:steps) {
for(i in 1:10) {
t = ((step-1)*10)+i
stepvec[t] = step
NPPSink[t] = sum(scenStockChange[step,]/10) + sum(as.numeric(deadWoodGeneration[step, ])/10) + sum(as.numeric(productWoodGeneration[step, ])/10)
livestockVariation[t] = sum(scenStockChange[step,]/10)
deadstockGeneration[t] = sum(as.numeric(deadWoodGeneration[step,])/10)
productstockGeneration[t] = sum(as.numeric(productWoodGeneration[step,])/10)
# if(smoothSink) {
# if((i<=5) && (step>1)) {
# p = (5+i)/10
# NPPSink[t] = NPPSink[t]*p + sum(scenStockChange[step-1,]/10)*(1-p)
# } else if((i>=6) && (step<steps)) {
# p = (i-6)/10
# NPPSink[t] = NPPSink[t]*(1-p) + sum(scenStockChange[step+1,]/10)*p
# }
# }
decompositionEmissions[t] = sum((1 - exp(-deadK))*deadComp)
productEndlifeEmissions[t] = sum((1 - exp(-productK))*prodComp)
deadComp = exp(-deadK)*deadComp + (as.numeric(deadWoodGeneration[step, ])/10)
prodComp = exp(-productK)*prodComp + (as.numeric(productWoodGeneration[step,])/10)
}
}
res_ha_yr = data.frame(Step = stepvec,
Years = yearvec,
LiveWoodChange = livestockVariation,
DeadWoodGeneration = deadstockGeneration,
ProductWoodGeneration = productstockGeneration,
NPPSink = NPPSink,
DecompositionEmissions = decompositionEmissions,
ProductEndlifeEmissions = productEndlifeEmissions,
EmissionNetEffect = decompositionEmissions + productEndlifeEmissions - NPPSink)
res_yr = data.frame(Step = stepvec,
Years = yearvec,
LiveWoodChange = livestockVariation*totalArea,
DeadWoodGeneration = deadstockGeneration*totalArea,
ProductWoodGeneration = productstockGeneration*totalArea,
NPPSink = NPPSink*totalArea,
DecompositionEmissions = decompositionEmissions*totalArea,
ProductEndlifeEmissions = productEndlifeEmissions*totalArea,
EmissionNetEffect = (decompositionEmissions + productEndlifeEmissions - NPPSink)*totalArea)
carbonBalance = list(Mg_yr = res_yr, Mg_ha_yr = res_ha_yr)
finalDead = deadComp
finalProducts = prodComp
return(list(initialDead = initialDead, initialProducts = initialProducts,
liveWoodChange = scenStockChange,
deadWoodGeneration = deadWoodGeneration,
productWoodGeneration = productWoodGeneration,
carbonBalance = carbonBalance,
finalDead = finalDead, finalProducts = finalProducts))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.