add SOC and Residues reporting

.buildlibrary

@@ -1,4 +1,4 @@
-ValidationKey: '45934350'
+ValidationKey: '4773840'
 AutocreateReadme: yes
 AcceptedWarnings:
 - 'Warning: package ''.*'' was built under R version'

CITATION.cff

@@ -2,8 +2,8 @@ cff-version: 1.2.0
 message: If you use this software, please cite it using the metadata from this file.
 type: software
 title: 'magpie4: MAgPIE outputs R package for MAgPIE version 4.x'
-version: 2.3.10
-date-released: '2024-06-11'
+version: 2.4.0
+date-released: '2024-06-17'
 abstract: Common output routines for extracting results from the MAgPIE framework
   (versions 4.x).
 authors:
@@ -65,7 +65,6 @@ authors:
 - family-names: Vartika
   given-names: Singh
 license: LGPL-3.0
-keywords: ~
 repository-code: https://github.com/pik-piam/magpie4
 doi: 10.5281/zenodo.1158582
 

DESCRIPTION

@@ -1,8 +1,8 @@
 Type: Package
 Package: magpie4
 Title: MAgPIE outputs R package for MAgPIE version 4.x
-Version: 2.3.10
-Date: 2024-06-11
+Version: 2.4.0
+Date: 2024-06-17
 Authors@R: c(
     person("Benjamin Leon", "Bodirsky", , "[email protected]", role = c("aut", "cre")),
     person("Florian", "Humpenoeder", , "[email protected]", role = "aut"),

NAMESPACE

@@ -39,8 +39,10 @@ export(PriceGHG)
 export(PrimSecdOtherLand)
 export(ResidueBiomass)
 export(ResidueUsage)
+export(Residues)
 export(RotationLength)
 export(SOM)
+export(SOM2)
 export(Seed)
 export(Timber)
 export(TimberDemand)
@@ -223,6 +225,7 @@ export(reportProductionNr)
 export(reportProtectedArea)
 export(reportProtein)
 export(reportRelativeHourlyLaborCosts)
+export(reportResidues)
 export(reportRotationLength)
 export(reportRuralDemandShares)
 export(reportSDG1)
@@ -233,6 +236,7 @@ export(reportSDG3)
 export(reportSDG6)
 export(reportSDG9)
 export(reportSOM)
+export(reportSOM2)
 export(reportTau)
 export(reportTc)
 export(reportTimber)

R/Residues.R

@@ -0,0 +1,228 @@
+#' @title Residues
+#' @description reads various crop residue (carbon) outputs out of a MAgPIE gdx file
+#' @export
+#'
+#' @param gdx         GDX file
+#' @param level       Level of regional aggregation; "reg" (regional), "glo" (global), "regglo" (regional and global)
+#' @param products    Selection of products (either "kcr" or "kres")
+#' @param waterAggr   Aggregate irrigated and non-irriagted production or not (boolean).
+#' @param output      Switch between different outputs: "biomass", "fieldBalance", "resDemand", all
+#'
+#' @return production as MAgPIE object (unit depends on attributes)
+#' @author Kristine Karstens, Michael Crawford
+#' @seealso \code{\link{ResidueBiomass}}
+#' @examples
+#'
+#'   \dontrun{
+#'     x <- Residues(gdx)
+#'   }
+#'   @importFrom madrat toolAggregate
+
+Residues <- function(gdx, level = "regglo", products = "kres", waterAggr = TRUE, output = "all"){
+
+  attr     <- "nr" # before dm
+  kres     <- readGDX(gdx, "kres")
+  kresfull <- c(kres, "res_nouse")
+  kcr      <- readGDX(gdx, "kcr")
+  kcr2kres <- readGDX(gdx, "kres_kcr")
+  missingKcr   <- setdiff(kcr, kcr2kres$kcr)
+  missingMap   <- data.frame(kres = rep("res_nouse", length(missingKcr)),
+                             kcr  = missingKcr)
+  kcr2kresfull <- rbind(kcr2kres, missingMap)
+  fmAttributes <- readGDX(gdx, "fm_attributes")[, , kres]
+
+
+
+  if (level %in% c("reg", "regglo", "glo")) {
+
+    biomassAgKcr   <- collapseNames(readGDX(gdx, "ov_res_biomass_ag")[, , "level"][, , attr])
+    biomassBgKcr   <- collapseNames(readGDX(gdx, "ov_res_biomass_bg")[, , "level"][, , attr])
+    recyclingKcr   <- collapseNames(readGDX(gdx, "ov18_res_ag_recycling")[, , "level"][, , attr])
+    burnKcr        <- collapseNames(readGDX(gdx, "ov_res_ag_burn", "ov18_res_ag_burn",
+                                            format = "first_found")[, , "level"][, , attr])
+    removalKcr     <- collapseNames(readGDX(gdx, "ov18_res_ag_removal")[, , "level"][, , attr])
+    feedKres       <- dimSums(collapseNames(readGDX(gdx, "ov_dem_feed")[, , "level"][, , kres] *
+                                              fmAttributes[, , kres][, , "c"], collapsedim = "type"), dim = 3.1)
+    materialKres   <- collapseNames(readGDX(gdx, "ov_dem_material")[, , "level"][, , kres] *
+                                      fmAttributes[, , kres][, , "c"], collapsedim = "type")
+    bioenergyKres  <- collapseNames(readGDX(gdx, "ov_dem_bioen")[, , "level"][, , kres] *
+                                      fmAttributes[, , kres][, , "c"], collapsedim = "type")
+
+    # When aggregating to kres no water specific results are possible
+    wSpecific <- any(grepl("^w$", getSets(biomassAgKcr)))
+    if (waterAggr == FALSE && "wSpecific != TRUE" ) {
+      note("No water-type specific output available.")
+    }
+
+    if (products == "kres"){
+
+      .aggregateKcr2Kres <- function(input) {
+        if (wSpecific) input <- dimSums(input, dim = "w")
+        madrat::toolAggregate(input, rel = kcr2kresfull, from = "kcr", to = "kres",
+                      dim = 3.1, partrel = TRUE)[, , kresfull]
+      }
+      biomassAgKres      <- .aggregateKcr2Kres(biomassAgKcr)
+      biomassBgKres      <- .aggregateKcr2Kres(biomassBgKcr)
+      recyclingKres      <- .aggregateKcr2Kres(recyclingKcr)
+      burnKres           <- .aggregateKcr2Kres(burnKcr)
+      removalKres        <- .aggregateKcr2Kres(removalKcr)
+
+      # ########### TESTING #################
+      # prod     <- collapseNames(readGDX(gdx, "ov_prod_reg")[, , "level"])
+      # prodKres <- prod[, , kres] * fmAttributes[, , kres][, , "c"]
+      # range(removalKres-prodKres)
+      #
+      # supplyKres     <- collapseNames(readGDX(gdx, "ov_supply")[, , "level"][, , kres]) * fmAttributes[, , kres][, , "c"]
+      # range(prodKres-supplyKres)
+      # where(abs(prodKres-supplyKres) > 1)$true
+      # # $individual
+      # # i     t       kall.attributes
+      # # LAM "LAM" "y1995" "res_cereals.c"
+      # # CAZ "CAZ" "y1995" "res_fibrous.c"
+      # # EUR "EUR" "y1995" "res_fibrous.c"
+      # # MEA "MEA" "y1995" "res_fibrous.c"
+      # # REF "REF" "y1995" "res_fibrous.c"
+      # # CAZ "CAZ" "y2000" "res_fibrous.c"
+      # # REF "REF" "y2000" "res_fibrous.c"
+      # #
+      # # $regions
+      # # [1] "LAM" "CAZ" "EUR" "MEA" "REF"
+      # #
+      # # $years
+      # # [1] "y1995" "y2000"
+      # #
+      # # $data
+      # # [1] "res_cereals.c" "res_fibrous.c"
+      #
+      #
+      # fbask     <- collapseNames(readGDX(gdx, "im_feed_baskets")[, , kres])
+      # bflow     <- collapseNames(readGDX(gdx, "fm_feed_balanceflow")[, , kres])
+      # kap     <- readGDX(gdx, "kap")
+      # time    <- getYears(prod)
+      #
+      # feedKres_real <-  dimSums(prod[, , kap] * fbask[, time, kres] + bflow[, time, kres], dim = 3.1) *
+      #   fmAttributes[, , kres][, , "c"]
+      #
+      # range(feedKres - feedKres_real)
+      # where(abs(feedKres - feedKres_real) > 1)$true
+      # # $individual
+      # # i     t       kall.attributes
+      # # MEA "MEA" "y1995" "res_cereals.c"
+      # # OAS "OAS" "y1995" "res_cereals.c"
+      # # SSA "SSA" "y1995" "res_cereals.c"
+      # # USA "USA" "y1995" "res_cereals.c"
+      # # IND "IND" "y1995" "res_fibrous.c"
+      # # LAM "LAM" "y1995" "res_fibrous.c"
+      # # USA "USA" "y1995" "res_fibrous.c"
+      # # USA "USA" "y2000" "res_fibrous.c"
+      # # USA "USA" "y2005" "res_fibrous.c"
+      # # REF "REF" "y2010" "res_fibrous.c"
+      # #
+      # # $regions
+      # # [1] "MEA" "OAS" "SSA" "USA" "IND" "LAM" "REF"
+      # #
+      # # $years
+      # # [1] "y1995" "y2000" "y2005" "y2010"
+      # #
+      # # $data
+      # # [1] "res_cereals.c" "res_fibrous.c"
+
+      ########### TESTING #################
+
+      # CHECK
+      # [x] FEED DMEAND OTHER fibrous -> basket equation
+      # [] check attribute equation
+
+      biomass      <- mbind(add_dimension(biomassAgKres, add = "plantpart", nm = "ag"),
+                            add_dimension(biomassBgKres, add = "plantpart", nm = "bg"))
+
+      fieldBalance <- mbind(add_dimension(recyclingKres, add = "fieldbalance", nm = "recycle"),
+                            add_dimension(burnKres,      add = "fieldbalance", nm = "burned"),
+                            add_dimension(removalKres,   add = "fieldbalance", nm = "removal"))
+
+      resDemand <- mbind(add_dimension(feedKres,      add = "usage", nm = "feed"),
+                         add_dimension(materialKres,  add = "usage", nm = "other_util"),
+                         add_dimension(bioenergyKres, add = "usage", nm = "bioenergy"))
+
+      check <- round(dimSums(removalKres, dim = 3) - dimSums(resDemand, dim = 3), 5)
+      if (any(check != 0)) {
+        warning("Sum over all residue removal options and residue removal in total are not matching.")
+      }
+
+      check <- round(biomassAgKres - dimSums(fieldBalance, dim = 3.1), 5)
+      if (any(check != 0)) {
+        warning("Sum over all residue usage options and ag-residue biomass in total are not matching.")
+      }
+
+    } else if (products == "kcr"){
+
+      biomass      <- mbind(add_dimension(biomassAgKcr, add = "plantpart", nm = "ag"),
+                            add_dimension(biomassBgKcr, add = "plantpart", nm = "bg"))
+
+      fieldBalance <- mbind(add_dimension(recyclingKcr, add = "fieldbalance", nm = "recycle"),
+                            add_dimension(burnKcr,      add = "fieldbalance", nm = "burned"),
+                            add_dimension(removalKcr,   add = "fieldbalance", nm = "removal"))
+
+      resDemand <- mbind(add_dimension(feedKres,      add = "usage", nm = "feed"),
+                         add_dimension(materialKres,  add = "usage", nm = "other_util"),
+                         add_dimension(bioenergyKres, add = "usage", nm = "bioenergy"))
+
+      if(wSpecific) {
+        weight <- dimSums(removalKcr[, , kcr2kres$kcr], dim = 3.2)
+      } else {
+        weight <- removalKcr[, , kcr2kres$kcr]
+      }
+      resDemand <- toolAggregate(resDemand, kcr2kres, weight = weight, from = "kres", to = "kcr", dim = 3.2)
+
+      getSets(resDemand, fulldim = FALSE)[3] <- "usage.kcr.attributes"
+      resDemand <- add_columns(resDemand, addnm = missingKcr, dim = 3.2)
+      resDemand[is.na(resDemand)] <- 0
+
+      check <- round(dimSums(removalKcr, dim = 3) - dimSums(resDemand, dim = 3), 5)
+      if(any(check != 0)) {
+        warning(paste0("Sum over all residue removal options and residue removal in total are not matching.
+                       Non-matching in the range of ", toString(range(check))))
+      }
+
+      check <- round(biomassAgKcr - dimSums(fieldBalance, dim = 3.1), 5)
+      if(any(check != 0)) {
+        warning(paste0("Sum over all residue usage options and ag-residue biomass in total are not matching.
+                       Non-matching in the range of ", toString(range(check))))
+      }
+    } else { stop(paste0("Product type ", products, " unknown.")) }
+
+    ### reg, regglo, glo aggregation
+    biomass      <- gdxAggregate(gdx, biomass,      to = level, absolute = TRUE)
+    fieldBalance <- gdxAggregate(gdx, fieldBalance, to = level, absolute = TRUE)
+    resDemand    <- gdxAggregate(gdx, resDemand,    to = level, absolute = TRUE)
+
+  } else {
+    stop("Level not available.")
+  }
+
+  if (output == "biomass") {
+    out <- biomass
+  } else if (output == "fieldBalance") {
+    out <- fieldBalance
+  } else if (output == "resDemand") {
+    out <- resDemand
+  } else if (output == "all") {
+    if (products == "kcr" && wSpecific) {
+      out <- mbind(collapseNames(dimSums(biomass, dim ="w")),
+                   collapseNames(dimSums(fieldBalance, dim ="w")),
+                   collapseNames(resDemand))
+    } else {
+      out <- mbind(collapseNames(biomass),
+                   collapseNames(fieldBalance),
+                   collapseNames(resDemand))
+    }
+  } else {
+    stop("Output type not available.")
+  }
+
+  if (products == "kcr" && waterAggr == TRUE && output %in% c("biomass", "fieldBalance")) {
+    out <- dimSums(out, dim = "w")
+  }
+
+  return(out)
+}

R/SOM2.R

@@ -0,0 +1,90 @@
+#' @title SOM2
+#' @description Calculates soil organic carbon stock size
+#'              based on a MAgPIE gdx file (for threepool realization)
+#'
+#' @export
+#'
+#' @param gdx GDX file
+#' @param type "stock" (default) for absoulte values,
+#'             "density" for per hectar values
+#' @param noncropAggr aggregate non cropland types to 'noncropland'
+#'                     (if FALSE all land types of pools59 will be reported)
+#' @param level Level of regional aggregation;
+#'                "reg" (regional),
+#'                "glo" (global),
+#'                "regglo" (regional and global)
+#'
+#' @return A MAgPIE object containing som values
+#' @author Kristine Karstens
+#' @examples
+#' \dontrun{
+#' x <- SOM2(gdx)
+#' }
+#'
+SOM2 <- function(gdx, type = "stock", level = "regglo", noncropAggr = TRUE) {
+
+  somStock        <- readGDX(gdx, "ov59_topsoilc_actualstate", select = list(type = "level"))
+  equStockCrop    <- readGDX(gdx, "ov59_topsoilc_crop_steadystate", select = list(type = "level"))
+  equStockNoncrop <- readGDX(gdx, "ov59_topsoilc_noncrop_steadystate", select = list(type = "level"))
+  years           <- getYears(somStock)
+  natStock        <- readGDX(gdx, "p59_topsoilc_naturalstate")[, years, ]
+  landuse         <- readGDX(gdx, "ov_land", select = list(type = "level"))
+  landuse         <- gdxAggregate(gdx, landuse, to = "reg", absolute = TRUE)
+
+  .aggregatePools <- function(x) {
+    return(dimSums(x, dim = "sPools59"))
+  }
+
+  somStock <- .aggregatePools (somStock)
+  equStock <- mbind(setNames(dimSums(equStockCrop, dim = 3), "crop"),
+                    .aggregatePools(equStockNoncrop))
+  natStock <- .aggregatePools(natStock)
+
+  if(noncropAggr == TRUE) {
+    .aggregateNoncrop <- function(x) {
+      nc59 <- readGDX(gdx, "noncropland59", types = "sets", react = "silent")
+      return(mbind(x[, , "crop"], setNames(dimSums(x[, , nc59], dim = 3.1), "natveg")))
+    }
+    somStock <- .aggregateNoncrop(somStock)
+    equStock <- .aggregateNoncrop(equStock)
+    natStock <- .aggregateNoncrop(natStock)
+    landuse  <- .aggregateNoncrop(landuse)
+  }
+
+  .addTotal <- function(x){
+    return(mbind(x, add_dimension(dimSums(x, dim = 3.1), add = "landuse", nm = "total")))
+  }
+
+  somStock <- .addTotal(somStock)
+  equStock <- .addTotal(equStock)
+  natStock <- .addTotal(natStock)
+  landuse  <- .addTotal(landuse)
+
+  if (type == "density") {
+
+    .calcDensity <- function(x, land) {
+       x <- x / land
+       x[is.infinite(x)] <- NA
+       return(x)
+    }
+
+    somStock <- .calcDensity(somStock, landuse)
+    equStock <- .calcDensity(equStock, landuse)
+    natStock <- .calcDensity(natStock, landuse)
+    somStock <- gdxAggregate(gdx, somStock, weight = landuse, to = level, absolute = FALSE)
+    equStock <- gdxAggregate(gdx, equStock, weight = landuse, to = level, absolute = FALSE)
+    natStock <- gdxAggregate(gdx, natStock, weight = landuse, to = level, absolute = FALSE)
+
+  } else {
+
+    somStock <- gdxAggregate(gdx, somStock, to = level, absolute = TRUE)
+    equStock <- gdxAggregate(gdx, equStock, to = level, absolute = TRUE)
+    natStock <- gdxAggregate(gdx, natStock, to = level, absolute = TRUE)
+  }
+
+  out <- mbind(add_dimension(somStock, dim = 3.1, add = "state", nm =  "actualstate"),
+               add_dimension(natStock, dim = 3.1, add = "state", nm = "naturalstate"),
+               add_dimension(equStock, dim = 3.1, add = "state", nm = "steadystate"))
+
+  return(out)
+}

R/getReport.R

@@ -128,6 +128,8 @@ getReport <- function(gdx, file = NULL, scenario = NULL, filter = c(1, 2, 7),
       "reportWaterAvailability(gdx)",
       "reportAAI(gdx)",
       "reportSOM(gdx)",
+      "reportSOM2(gdx)",
+      "reportResidues(gdx)",
       "reportGrowingStock(gdx, indicator='relative')",
       "reportGrowingStock(gdx, indicator='absolute')",
       "reportSDG1(gdx)",