bugfix emisCO2

.buildlibrary

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

CITATION.cff

@@ -2,7 +2,7 @@ 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.10.1
+version: 2.10.2
 date-released: '2024-09-05'
 abstract: Common output routines for extracting results from the MAgPIE framework
   (versions 4.x).

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: magpie4
 Title: MAgPIE outputs R package for MAgPIE version 4.x
-Version: 2.10.1
+Version: 2.10.2
 Date: 2024-09-05
 Authors@R: c(
     person("Benjamin Leon", "Bodirsky", , "[email protected]", role = c("aut", "cre")),

R/carbonstock.R

@@ -38,11 +38,12 @@ carbonstock <- function(gdx, file=NULL, level="cell", sum_cpool=TRUE, sum_land=T
   if (!is.null(subcategories)) {
     if ("crop" %in% subcategories) {
       croparea_land <- readGDX(gdx, "ov_area", select = list(type = "level"))
-      fallow_land <- readGDX(gdx, "ov_fallow", select = list(type = "level"), react = "silent")
-      croptree_land <- readGDX(gdx, "ov29_treecover", select = list(type = "level"), react = "silent")
-      if(!is.null(fallow_land) && !is.null(croptree_land)) {
+      fallow_land <- readGDX(gdx, "ov_fallow", select = list(type = "level"))
+      croptree_land <- readGDX(gdx, "ov_treecover", select = list(type = "level"))
+      if(sum(fallow_land) > 0 || sum(croptree_land) > 0) {
 
         p29_carbon_density_ac <- readGDX(gdx, "p29_carbon_density_ac")
+        croptree_land_ac <- readGDX(gdx, "ov29_treecover", select = list(type = "level"))
         fm_carbon_density <- collapseNames(readGDX(gdx, "fm_carbon_density")[,getYears(p29_carbon_density_ac),getNames(p29_carbon_density_ac,dim=2)][,,"crop"])
         names(dimnames(fm_carbon_density))[[3]] <- "ag_pools"
 
@@ -53,7 +54,7 @@ carbonstock <- function(gdx, file=NULL, level="cell", sum_cpool=TRUE, sum_land=T
         fallow_carbon_stock <- fallow_land * fm_carbon_density
         fallow_carbon_stock <- add_dimension(fallow_carbon_stock, dim = 3.1, add = "land", "fallow")
 
-        croptree_carbon_stock <- dimSums(croptree_land * p29_carbon_density_ac, dim = "ac")
+        croptree_carbon_stock <- dimSums(croptree_land_ac * p29_carbon_density_ac, dim = "ac")
         croptree_carbon_stock <- add_dimension(croptree_carbon_stock, dim = 3.1, add = "land", "treecover")
 
         crop_carbon_stock <- mbind(croparea_carbon_stock,fallow_carbon_stock,croptree_carbon_stock)
@@ -64,7 +65,7 @@ carbonstock <- function(gdx, file=NULL, level="cell", sum_cpool=TRUE, sum_land=T
             # compose crop areas
             crop <- mbind(add_dimension(dimSums(croparea_land, dim = 3), dim = 3.1, add = "land", "area"),
                           add_dimension(fallow_land, dim = 3.1, add = "land", "fallow"),
-                          add_dimension(dimSums(croptree_land, dim = "ac"), dim = 3.1, add = "land", "treecover"))
+                          add_dimension(croptree_land, dim = 3.1, add = "land", "treecover"))
 
             # recalculate ov59_som_pool for checking
             ov59_som_pool_check <- readGDX(gdx, "ov59_som_pool", select = list(type = "level"))
@@ -92,7 +93,7 @@ carbonstock <- function(gdx, file=NULL, level="cell", sum_cpool=TRUE, sum_land=T
             zz[,,] <- 0
             zz[,,"area"] <- dimSums(croparea_land * readGDX(gdx,"i59_cratio"), dim=3) * readGDX(gdx,"f59_topsoilc_density")[,getYears(zz),]
             zz[,,"fallow"] <- fallow_land * readGDX(gdx,"i59_cratio_fallow") * readGDX(gdx,"f59_topsoilc_density")[,getYears(zz),]
-            zz[,,"treecover"] <- dimSums(croptree_land, dim = "ac") * readGDX(gdx,"i59_cratio_treecover") * readGDX(gdx,"f59_topsoilc_density")[,getYears(zz),]
+            zz[,,"treecover"] <- croptree_land * readGDX(gdx,"i59_cratio_treecover") * readGDX(gdx,"f59_topsoilc_density")[,getYears(zz),]
             if(abs(sum(dimSums(zz,dim=3) - collapseNames(ov59_som_target[,,"crop"]))) > 1e-6) warning("differences in ov59_som_target detected")
             ov59_som_target <- zz
 

R/emisCO2.R

@@ -135,63 +135,217 @@ emisCO2 <- function(gdx, file = NULL, level = "cell", unit = "gas",
 
   ###
   # compose list of carbon-densities per land-use type, compartment
+
   composeDensities <- function() {
 
-    cs <- carbonstock(gdx, level = "cell", subcategories = c("crop", "other", "forestry"), sum_land = FALSE, sum_cpool = FALSE)
-    cs[cs < 1e-8] <- 0
-    a  <- land(gdx, level = "cell", subcategories = c("crop", "other", "forestry"))
-    a[a < 1e-8] <- 0
-    cd <- cs / a
-    cd[is.na(cd)] <- 0
-    cd[is.infinite(cd)] <- 0
-
-    cropArea <- cd[, , "crop_area"]
-    cropFallow <- cd[, , "crop_fallow"]
-    cropTreecover <- cd[, , "crop_treecover"]
-    pasture <- cd[, , "past"]
-    forestry <- cd[, , c("forestry_ndc", "forestry_aff", "forestry_plant")]
-    primforest <- cd[, , "primforest"]
-    secdforest <- cd[, , "secdforest"]
-    other <- cd[, , c("other_othernat", "other_youngsecdf")]
-    urban <- cd[, , "urban"]
-
-    forestryAC <- readGDX(gdx, "p32_carbon_density_ac")[, years, ]
-    getSets(forestryAC)["d3.1"] <- "land"
-    getNames(forestryAC, dim = 1) <- paste("forestry", getNames(forestryAC, dim = 1), sep = "_")
-    forestry <- add_dimension(forestry, dim = 3.2, add = "ac", getNames(forestryAC, dim = "ac"))
-    forestry[, , getNames(forestryAC)] <- forestryAC
-
-    secdforestAC <- readGDX(gdx, "pm_carbon_density_secdforest_ac", "pm_carbon_density_ac", format = "first_found")[, years, ]
-    secdforest <- add_dimension(secdforest, dim = 3.2, add = "ac", getNames(secdforestAC, dim = "ac"))
-    secdforest[, , getNames(secdforestAC)] <- secdforestAC
-
-    cropTreecoverAC <- readGDX(gdx, "p29_carbon_density_ac", react = "silent")
-    if (is.null(cropTreecoverAC)) {
-      cropTreecoverAC <- secdforestAC
-      cropTreecoverAC[, , ] <- 0
+    # currently, for ac types all the soil carbon densities are assumed equal
+    .addSOM <- function(x, soilc) {
+
+      getSets(x)["d3.3"] <- "c_pools"
+      x <- add_columns(x, addnm = "soilc", dim = "c_pools", fill = 1)
+      soilc <- x[, , "soilc"] * soilc
+      x <- mbind(x[, , "soilc", invert = TRUE], soilc)
+
+      return(x)
     }
-    cropTreecover <- add_dimension(cropTreecover, dim = 3.2, add = "ac", getNames(cropTreecoverAC, dim = "ac"))
-    cropTreecover[, , getNames(cropTreecoverAC)] <- cropTreecoverAC
 
-    otherAC <- readGDX(gdx, "p35_carbon_density_other", "pm_carbon_density_ac", format = "first_found")[, years, ]
-    if (getSets(otherAC)["d3.1"] == "othertype35") {
-      getSets(otherAC)["d3.1"] <- "land"
-      getNames(otherAC, dim = 1) <- paste("other", getNames(otherAC, dim = 1), sep = "_")
+    # --- non-age class carbon densities
+    carbonDensities <- readGDX(gdx, "fm_carbon_density")[, years, ]
+    croparea          <- carbonDensities[, , "crop"][,,c("vegc","litc")]
+    names(dimnames(croparea)) <- c("j.region","t","land.ag_pools")
+    getNames(croparea, dim = "land") <- "crop_area"
+    fallow            <- croparea
+    getNames(fallow, dim = "land") <- "crop_fallow"
+    pasture           <- carbonDensities[, , "past"]
+    urban             <- carbonDensities[, , "urban"]
+    primforest        <- carbonDensities[, , "primforest"]
+
+    # --- age class carbon densities, excl forestry
+
+    secdforest <- readGDX(gdx, "pm_carbon_density_secdforest_ac", "pm_carbon_density_ac", format = "first_found")[, years, ]
+    secdforest <- add_dimension(secdforest, dim = 3.1, add = "land", nm = "secdforest")
+
+    other <- readGDX(gdx, "ov_land_other", select = list(type = "level"), react = "silent")
+    if(!is.null(other)) {
+      other <- readGDX(gdx, "p35_carbon_density_other", react = "silent")
+      getSets(other)["d3.1"] <- "land"
+      other <- other[, years, ]
+      getNames(other,dim="land") <- paste("other", getNames(other,dim="land"), sep = "_")
+    } else {
+      other <- readGDX(gdx, "pm_carbon_density_ac")[, years, ]
+      other <- add_dimension(other, dim = 3.1, add = "land", nm = "other")
     }
-    other <- add_dimension(other, dim = 3.2, add = "ac", getNames(otherAC, dim = "ac"))
-    other[, , getNames(otherAC)] <- otherAC
 
-    densities <- list(cropArea      = cropArea,
-                  cropFallow    = cropFallow,
-                  cropTreecover = cropTreecover,
-                  pasture       = pasture,
-                  forestry      = forestry,
-                  primforest    = primforest,
-                  secdforest    = secdforest,
-                  other         = other,
-                  urban         = urban)
+    croptree <- readGDX(gdx, "p29_carbon_density_ac", react = "silent")
+    if(!is.null(croptree)) {
+      croptree <- add_dimension(croptree[, years, ], dim = 3.1, add = "land", nm = "crop_treecover")
+    } else {
+      croptree <- secdforest
+      croptree[,,] <- 0
+      getNames(croptree, dim = 1) <- "crop_treecover"
+    }
+
+    # --- forestry
+    forestry <- readGDX(gdx, "p32_carbon_density_ac", react = "silent")
+    getSets(forestry)["d3.1"] <- "land"
+    getNames(forestry, dim = 1) <- paste("forestry", getNames(forestry, dim = 1), sep = "_")
+
+    # --- SOM emissions
+    dynSom <- !is.null(readGDX(gdx, "ov59_som_pool", react = "silent"))
+    if (dynSom) {
+
+      #topsoilCarbonDensities <- readGDX(gdx, "f59_topsoilc_density")[, years, ]
+      subsoilCarbonDensities <- readGDX(gdx, "i59_subsoilc_density")[, years, ]
+
+      ov59_som_pool <- readGDX(gdx, "ov59_som_pool", select = list(type = "level"))
+      ov_land <- readGDX(gdx, "ov_land", select = list(type = "level"))
+      topsoilCarbonDensities <- ov59_som_pool / ov_land
+      topsoilCarbonDensities[is.na(topsoilCarbonDensities)] <- 0
+      topsoilCarbonDensities[is.infinite(topsoilCarbonDensities)] <- 0
+
+      # croparea[, , "soilc"]   <- topsoilCarbonDensities[,,"crop"] + subsoilCarbonDensities
+      # fallow[, , "soilc"]     <- topsoilCarbonDensities[,,"crop"] + subsoilCarbonDensities
+      pasture[, , "soilc"]    <- topsoilCarbonDensities[,,"past"] + subsoilCarbonDensities
+      urban[, , "soilc"]      <- topsoilCarbonDensities[,,"urban"] + subsoilCarbonDensities
+      primforest[, , "soilc"] <- topsoilCarbonDensities[,,"primforest"] + subsoilCarbonDensities
+
+      ### BEGIN croptreecover calculations
+
+      # compose crop areas
+      croparea_land <- readGDX(gdx, "ov_area", select = list(type = "level"))
+      fallow_land <- readGDX(gdx, "ov_fallow", select = list(type = "level"))
+      croptree_land <- readGDX(gdx, "ov_treecover", select = list(type = "level"))
+
+      if(sum(fallow_land) > 0 || sum(croptree_land) > 0) {
+        crop <- mbind(add_dimension(dimSums(croparea_land, dim = 3), dim = 3.1, add = "land", "crop_area"),
+                      add_dimension(fallow_land, dim = 3.1, add = "land", "crop_fallow"),
+                      add_dimension(croptree_land, dim = 3.1, add = "land", "crop_treecover"))
+        # recalculate ov59_som_pool for checking
+        ov59_som_pool_check <- readGDX(gdx, "ov59_som_pool", select = list(type = "level"))
+
+        ov59_som_target <- readGDX(gdx, "ov59_som_target", select = list(type = "level"))
+        i59_lossrate <- readGDX(gdx,"i59_lossrate")
+        p59_carbon_density <- readGDX(gdx,"p59_carbon_density")[,getYears(i59_lossrate),]
+        names(dimnames(p59_carbon_density))[[3]] <- "land_from"
+        names(dimnames(p59_carbon_density))[[2]] <- "t"
+        ov_lu_transitions <- readGDX(gdx, "ov_lu_transitions", select = list(type = "level"))
+        names(dimnames(ov_lu_transitions))[[3]] <- "land_from.land"
+        ov59_som_pool_intermediate <- (1 - i59_lossrate) * dimSums(p59_carbon_density * ov_lu_transitions, dim="land_from")
+
+        ov59_som_pool <- ov59_som_target * i59_lossrate + ov59_som_pool_intermediate
+
+        if(abs(sum(ov59_som_pool - ov59_som_pool_check)) > 1e-6) warning("differences in ov59_som_pool detected")
+
+        # split crop som pool based with crop (crop_area, crop_fallow, crop_tree) as weight
+        w <- crop / dimSums(crop, dim=3)
+        w[is.na(w)] <- 1/3
+        ov59_som_pool_intermediate <- collapseNames(ov59_som_pool_intermediate[,,"crop"]) * w
+
+        # recalculate ov59_som_target for crop_area, crop_fallow and crop_tree
+        zz <- ov59_som_pool_intermediate
+        zz[,,] <- 0
+        zz[,,"crop_area"] <- dimSums(croparea_land * readGDX(gdx,"i59_cratio"), dim=3) * readGDX(gdx,"f59_topsoilc_density")[,getYears(zz),]
+        zz[,,"crop_fallow"] <- fallow_land * readGDX(gdx,"i59_cratio_fallow") * readGDX(gdx,"f59_topsoilc_density")[,getYears(zz),]
+        zz[,,"crop_tree"] <- croptree_land * readGDX(gdx,"i59_cratio_treecover") * readGDX(gdx,"f59_topsoilc_density")[,getYears(zz),]
+        if(abs(sum(dimSums(zz,dim=3) - collapseNames(ov59_som_target[,,"crop"]))) > 1e-6) warning("differences in ov59_som_target detected")
+        ov59_som_target <- zz
+
+        # recalculate ov59_som_pool for crop_area, crop_fallow and crop_tree
+        ov59_som_pool <- ov59_som_target * i59_lossrate + ov59_som_pool_intermediate
+
+        # derive top soil density for crop_area, crop_fallow and crop_tree
+        top <- ov59_som_pool / crop
+        top[is.na(top)] <- 0
+        top[is.infinite(top)] <- 0
+        sub <- readGDX(gdx, "i59_subsoilc_density")[,getYears(top),]
+        soilc <- crop * (top + sub)
+        soilc <- add_dimension(soilc,dim=3.2,add="c_pools",nm = "soilc")
+
+        croparea    <- mbind(croparea, soilc[,,"crop_area"])
+        names(dimnames(croparea))[[3]] <- "land.c_pools"
+        fallow    <- mbind(fallow, soilc[,,"crop_fallow"])
+        names(dimnames(fallow))[[3]] <- "land.c_pools"
+        croptree    <- .addSOM(croptree, soilc[,,"crop_treecover"])
+      } else {
+        croparea    <- mbind(croparea, setNames(topsoilCarbonDensities[,,"crop"] + subsoilCarbonDensities,"crop_area.soilc"))
+        names(dimnames(croparea))[[3]] <- "land.c_pools"
+        getSets(fallow)["d3.2"] <- "c_pools"
+        fallow <- add_columns(fallow, "soilc", dim = 3.2, fill = 0)
+        croptree    <- .addSOM(croptree, 0)
+      }
+      ### END croptreecover calculations
+
+      secdforestSOM <- topsoilCarbonDensities[, , "secdforest"] + subsoilCarbonDensities
+      secdforest    <- .addSOM(secdforest, secdforestSOM)
+
+      if(all(c("other_othernat","other_youngsecdf") %in% getNames(other,dim="land"))) {
+        othernatSOM <- setNames(topsoilCarbonDensities[, , "other"] + subsoilCarbonDensities, "other_othernat")
+        youngsecdfSOM <- setNames(topsoilCarbonDensities[, , "other"] + subsoilCarbonDensities, "other_youngsecdf")
+        otherSOM <- mbind(othernatSOM,youngsecdfSOM)
+      } else {
+        otherSOM <- topsoilCarbonDensities[, , "other"] + subsoilCarbonDensities
+      }
+      other    <- .addSOM(other, otherSOM)
+
+      forestrySOM <- collapseNames(topsoilCarbonDensities[, , "forestry"]) + subsoilCarbonDensities
+      forestry    <- .addSOM(forestry, forestrySOM)
+
+    } else {
+
+      # --- cropland and pasture
+      topsoilCarbonDensities <- readGDX(gdx, "i59_topsoilc_density")[, years, ] # i59, not f59 as in dynamic realization
+      subsoilCarbonDensities <- readGDX(gdx, "i59_subsoilc_density")[, years, ]
+
+      cropareaSOM <- topsoilCarbonDensities + subsoilCarbonDensities
+      cropareaSOM <- add_dimension(cropareaSOM, dim = 3, add = "c_pools", nm = "soilc")
+      getSets(croparea)["d3.2"] <- "c_pools"
+      croparea <- add_columns(croparea, addnm = "soilc", dim = 3.2, fill = 0)
+      croparea[,,"soilc"] <- cropareaSOM
+
+      fallowSOM <- topsoilCarbonDensities + subsoilCarbonDensities
+      fallowSOM <- add_dimension(fallowSOM, dim = 3, add = "c_pools", nm = "soilc")
+      getSets(fallow)["d3.2"] <- "c_pools"
+      fallow <- add_columns(fallow, addnm = "soilc", dim = 3.2, fill = 0)
+      fallow[,,"soilc"] <- fallowSOM
+
+      # --- all other types
+      carbonDensities       <- readGDX(gdx, "fm_carbon_density")[, years, ]
+
+      pasture[, , "soilc"]    <- carbonDensities[, , "past"][, , "soilc"]
+      urban[, , "soilc"]      <- carbonDensities[, , "urban"][, , "soilc"]
+      primforest[, , "soilc"] <- carbonDensities[, , "primforest"][, , "soilc"]
+
+      croptreeSOM <- collapseDim(carbonDensities[, , "secdforest"][, , "soilc"], dim = "land")
+      croptree    <- .addSOM(croptree, croptreeSOM)
+
+      secdforestSOM <- collapseDim(carbonDensities[, , "secdforest"][, , "soilc"], dim = "land")
+      secdforest    <- .addSOM(secdforest, secdforestSOM)
+
+      if(all(c("other_othernat","other_youngsecdf") %in% getNames(other,dim="land"))) {
+        othernatSOM <- collapseDim(carbonDensities[, , "other"][, , "soilc"], dim = "land")
+        othernatSOM <- add_dimension(othernatSOM, dim = 3.1, add = "land", nm = "other_othernat")
+        youngsecdfSOM <- collapseDim(carbonDensities[, , "secdforest"][, , "soilc"], dim = "land")
+        youngsecdfSOM <- add_dimension(youngsecdfSOM, dim = 3.1, add = "land", nm = "other_youngsecdf")
+        otherSOM <- mbind(othernatSOM,youngsecdfSOM)
+      } else {
+        otherSOM <- collapseDim(carbonDensities[, , "other"][, , "soilc"], dim = "land")
+      }
+      other    <- .addSOM(other, otherSOM)
+
+      forestrySOM <- collapseDim(carbonDensities[, , "forestry"][, , "soilc"], dim = "land")
+      forestry    <- .addSOM(forestry, forestrySOM)
+
+    }
 
-    return(densities)
+    densities <- list(cropArea      = croparea,
+                      cropFallow    = fallow,
+                      cropTreecover = croptree,
+                      pasture       = pasture,
+                      forestry      = forestry,
+                      primforest    = primforest,
+                      secdforest    = secdforest,
+                      other         = other,
+                      urban         = urban)
   }
 
   ###

R/land.R

@@ -51,13 +51,11 @@ land <- function(gdx, file = NULL, level = "reg", types = NULL, subcategories =
     if (!is.null(subcategories)) {
       if ("crop" %in% subcategories) {
         croparea_land <- readGDX(gdx, "ov_area", select = list(type = "level"))
-        fallow_land <- readGDX(gdx, "ov_fallow", select = list(type = "level"), react = "silent")
-        croptree_land <- readGDX(gdx, "ov29_treecover", select = list(type = "level"), react = "silent")
-        if(is.null(fallow_land)) fallow_land <- new.magpie(getCells(croparea_land), getYears(croparea_land), NULL, fill = 0, sets = c("j.region","t","d3"))
-        if(is.null(croptree_land)) croptree_land <- new.magpie(getCells(croparea_land), getYears(croparea_land), NULL, fill = 0, sets = c("j.region","t","d3"))
+        fallow_land <- readGDX(gdx, "ov_fallow", select = list(type = "level"))
+        croptree_land <- readGDX(gdx, "ov_treecover", select = list(type = "level"))
         crop <- mbind(add_dimension(dimSums(croparea_land, dim = 3), dim = 3.1, add = "land", "area"),
                       add_dimension(fallow_land, dim = 3.1, add = "land", "fallow"),
-                      add_dimension(dimSums(croptree_land, dim = 3), dim = 3.1, add = "land", "treecover"))
+                      add_dimension(croptree_land, dim = 3.1, add = "land", "treecover"))
         getNames(crop,dim=1) <- paste("crop",getNames(crop,dim=1),sep="_")
 
         #crop <- add_dimension(crop, dim = 3.1, add = "land", "crop")