R/GrowingStock.R

In pik-piam/magpie4: MAgPIE outputs R package for MAgPIE version 4.x

#' @title GrowingStock
#' @description reads woody growing stock out of a MAgPIE gdx file
#'
#' @export
#'
#' @param gdx GDX file
#' @param file a file name the output should be written to using write.magpie
#' @param level Level of regional aggregation; "cell", "reg" (regional), "glo" (global), "regglo" (regional and global) or any secdforest aggregation level defined in superAggregate
#' @param indicator If the reported numbers are relative (mio m3/ha) or absolute (mio. m3). Default is relative.
#' @details Growing stock for producing woody materials consist of growing stock from plantations (forestry), secondary and primary forest as well as other land (natveg)
#' @return Growing stock in m3 per ha
#' @author Abhijeet Mishra
#' @importFrom gdx readGDX out
#' @importFrom magclass clean_magpie dimSums collapseNames setYears write.magpie
#' @importFrom luscale superAggregate
#' @examples
#'
#' \dontrun{
#' x <- GrowingStock(gdx)
#' }
#'
GrowingStock <- function(gdx, file = NULL, level = "regglo", indicator = "relative") {
  if (level == "regglo") {
    ac_sub <- readGDX(gdx, "ac")
    #ac_sub <- ac_sub[-1] ## Dont count ac0
    
    wood_density <- 0.6 ## tDM/m3
    ## Multiple sources for this number
    ## Check Table 2.8.1 in 2013 Revised Supplementary Methods and Good Practice Guidance Arising from the Kyoto Protocol
    
    ## Read timber yield, this is already upscale in the model after calibration to FAO GS.
    ## Divide by density to convert from tDM/ha to m3/ha
    pm_timber_yield <- readGDX(gdx, "pm_timber_yield") / wood_density ### mio. ha * tDM per ha / tDM per m3 = mio. m3
    
    ## Land information - cluster level
    land_plantations       <- collapseNames(readGDX(gdx, "ov32_land", "ov_land_fore", select = list(type = "level"))[, , "plant"][, , ac_sub])
    land_afforest          <- collapseNames(dimSums(readGDX(gdx, "ov32_land", "ov_land_fore", select = list(type = "level"))[, , "plant", invert = T][, , ac_sub], dim = "type32"))
    land_secdforest        <- collapseNames(readGDX(gdx, "ov35_secdforest", select = list(type = "level"))[, , ac_sub])
    land_primforest        <- setNames(collapseNames(readGDX(gdx, "ov_land", select = list(type = "level"))[, , "primforest"]), "acx")
    land_other             <- collapseNames(readGDX(gdx, "ov35_other", select = list(type = "level"))[, , ac_sub])
    land_natfor            <- land_secdforest
    land_natfor[, , "acx"] <- land_natfor[, , "acx"] + land_primforest
    land_forest_total      <- land_plantations + land_afforest + land_natfor
    
    ## Only use ac_sub
    land_plantations       <- land_plantations[, , ac_sub]
    land_afforest          <- land_afforest[, , ac_sub]
    land_secdforest        <- land_secdforest[, , ac_sub]
    land_natfor            <- land_natfor[, , ac_sub]
    land_other             <- land_other[, , ac_sub]
    land_forest_total      <- land_forest_total[, , ac_sub]
    pm_timber_yield        <- pm_timber_yield[, , ac_sub]
    
    ## Yields x land = standing growing stocks (mio. m3)

    growing_stock_main <- superAggregate(
      data = mbind(
        setNames(dimSums(land_plantations * pm_timber_yield[, , "forestry"]             , dim = "ac"), "plantations"),
        setNames(dimSums(land_afforest    * pm_timber_yield[, , "secdforest"]           , dim = "ac"), "afforestation"),
        setNames(dimSums(land_secdforest  * pm_timber_yield[, , "secdforest"]           , dim = "ac"), "secdforest"),
        setNames(dimSums(land_primforest  * pm_timber_yield[, , "primforest"][, , "acx"], dim = 3)   , "primforest"),
        setNames(dimSums(land_other       * pm_timber_yield[, , "other"]                , dim = "ac"), "other")
      ),
      aggr_type = "sum", level = level
    )
    
    ## Create GS for natural and total forests, as this is absolute amount still it can be added up
    growing_stock_natfor <- setNames(growing_stock_main[, , "secdforest"] + growing_stock_main[, , "primforest"], "natfor")
    
    growing_stock_forest_total <- setNames(growing_stock_natfor + growing_stock_main[, , "plantations"] + growing_stock_main[, , "afforestation"], "forest")
    
    a <- mbind(growing_stock_forest_total, growing_stock_natfor, growing_stock_main)
    
    ## In case relative numbers are needed, divide by aggregated land information
    
    if (indicator == "relative") {
      growing_stock_forest_total  <- growing_stock_forest_total              / superAggregate(data = dimSums(land_forest_total, dim = "ac"), aggr_type = "sum", level = level)
      growing_stock_natfor        <- growing_stock_natfor                    / superAggregate(data = dimSums(land_natfor, dim = "ac"),       aggr_type = "sum", level = level)
      growing_stock_plantations   <- growing_stock_main[, , "plantations"]   / superAggregate(data = dimSums(land_plantations, dim = "ac"),  aggr_type = "sum", level = level)
      growing_stock_afforestation <- growing_stock_main[, , "afforestation"] / superAggregate(data = dimSums(land_afforest, dim = "ac"),     aggr_type = "sum", level = level)
      growing_stock_primforest    <- growing_stock_main[, , "primforest"]    / superAggregate(data = dimSums(land_primforest, dim = 3),      aggr_type = "sum", level = level)
      growing_stock_secdforest    <- growing_stock_main[, , "secdforest"]    / superAggregate(data = dimSums(land_secdforest, dim = "ac"),   aggr_type = "sum", level = level)
      growing_stock_other         <- growing_stock_main[, , "other"]         / superAggregate(data = dimSums(land_other, dim = "ac"),        aggr_type = "sum", level = level)
      growing_stock_main          <- mbind(growing_stock_plantations, growing_stock_afforestation, growing_stock_secdforest, growing_stock_primforest, growing_stock_other)
      a <- mbind(growing_stock_forest_total, growing_stock_natfor, growing_stock_main)
    }
  } else {
    message("ERROR - wrong regions")
    a <- NULL
  }
  
  out(a, file)
}