#' @title carbonstock
#' @description reads carbon stocks out of a MAgPIE gdx file
#'
#' @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 sum_cpool sum over carbon pool dimension (default = TRUE)
#' @param sum_land sum over land type dimension (default = TRUE)
#' @param subcategories NULL or vector of strings. If NULL, no subcategories are returned.
#' Meaningful options are "crop, "forestry" and "other"
#' @param stockType carbon stock type (default = "actual").
#' Options: "actual", "previousLandPattern" and "previousCarbonDensity".
#' @details carbon pools consist of vegetation carbon (vegc), litter carbon (litc) and soil carbon (soilc)
#' @return carbon stocks in MtC
#' @export
#' @author Florian Humpenoeder
#'
#' @importFrom magclass dimSums collapseNames write.magpie
#' @importFrom luscale superAggregate
#' @importFrom utils head
#' @examples
#' \dontrun{
#' x <- carbonstock(gdx)
#' }
carbonstock <- function(gdx, file = NULL, level = "cell", sum_cpool = TRUE,
sum_land = TRUE, subcategories = NULL, stockType = "actual") {
# read in carbon stocks
a <- readGDX(gdx, "ov_carbon_stock", select = list(type = "level"), react = "silent")
names(dimnames(a))[1] <- "j"
if (length(getDim(stockType, a)) > 0) {
a <- collapseNames(a[, , stockType])
}
dyn_som <- !is.null(readGDX(gdx, "ov59_som_pool", react = "silent"))
.testDiff <- function(x = NULL, xCheck = NULL, accuracy = 1e-6, message = NULL){
diff <- abs(x - xCheck)
if (any(diff > accuracy)){
where <- where(diff > accuracy)$true
warning(paste(message, "\n", paste(where$individual, collapse = " ")), call. = FALSE)
}
}
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"))
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"
croparea_carbon_stock <- dimSums(croparea_land, dim = 3) * fm_carbon_density
names(dimnames(croparea_carbon_stock))[[3]] <- "ag_pools"
croparea_carbon_stock <- add_dimension(croparea_carbon_stock, dim = 3.1, add = "land", "area")
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_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)
if (!"soilc" %in% getNames(crop_carbon_stock, dim = "ag_pools")) {
names(dimnames(crop_carbon_stock))[[3]] <- "land.c_pools"
if (dyn_som) {
# 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(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"))
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
.testDiff(ov59_som_pool, ov59_som_pool_check, 1e-3, "differences in ov59_som_pool detected")
# split crop som pool based with crop (area, fallow, treecover) 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 area, fallow, treecover and if available scm
i59_scm_target <- readGDX(gdx, "i59_scm_target", react = "silent")
i59_cratio_scm <- readGDX(gdx, "i59_scm_ratio", "i59_cratio_scm", react = "silent")
if (is.null(i59_scm_target)) {
i59_scm_target <- 0
i59_cratio_scm <- 1
}
zz <- ov59_som_pool_intermediate
zz[, , ] <- 0
zz[, , "area"] <- (dimSums(croparea_land * readGDX(gdx, "i59_cratio"), dim = 3) +
dimSums(croparea_land * i59_scm_target *
readGDX(gdx, "i59_cratio") * (i59_cratio_scm - 1), 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"] <- croptree_land * readGDX(gdx, "i59_cratio_treecover") *
readGDX(gdx, "f59_topsoilc_density")[, getYears(zz), ]
.testDiff(dimSums(zz, dim = 3), collapseNames(ov59_som_target[, , "crop"]),
1e-3, "differences in ov59_som_target detected")
ov59_som_target <- zz
# recalculate ov59_som_pool for area, fallow and treecover
ov59_som_pool <- ov59_som_target * i59_lossrate + ov59_som_pool_intermediate
# derive top soil density for area, fallow and treecover
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")
} else {
soilc_croparea <- dimSums(croparea_land, dim = 3) *
(readGDX(gdx, "i59_topsoilc_density")[, getYears(crop_carbon_stock), ] +
readGDX(gdx, "i59_subsoilc_density")[, getYears(crop_carbon_stock), ])
soilc_croparea <- add_dimension(soilc_croparea, dim = 3.1, add = "land", "area")
soilc_fallow <- fallow_land * (readGDX(gdx, "i59_topsoilc_density")[, getYears(crop_carbon_stock), ] +
readGDX(gdx, "i59_subsoilc_density")[, getYears(crop_carbon_stock), ])
soilc_fallow <- add_dimension(soilc_fallow, dim = 3.1, add = "land", "fallow")
soilc_croptree <- croptree_land *
collapseNames(readGDX(gdx, "fm_carbon_density")[, getYears(crop_carbon_stock), "secdforest"][, , "soilc"])
soilc_croptree <- add_dimension(soilc_croptree, dim = 3.1, add = "land", "treecover")
soilc <- mbind(soilc_croparea, soilc_fallow, soilc_croptree)
soilc <- add_dimension(soilc, dim = 3.2, add = "c_pools", nm = "soilc")
}
crop_carbon_stock <- mbind(crop_carbon_stock, soilc)
}
# check
.testDiff(dimSums(crop_carbon_stock, dim = 3.1), collapseNames(a[, , "crop"]),
1e-3, "Differences in crop land detected!")
getNames(crop_carbon_stock, dim = 1) <- paste("crop", getNames(crop_carbon_stock, dim = 1), sep = "_")
crop <- crop_carbon_stock
} else {
a1 <- a[, , "crop"]
getNames(a1, dim = 1) <- "area"
a2 <- a1
a2[, , ] <- 0
getNames(a2, dim = 1) <- "fallow"
a3 <- a1
a3[, , ] <- 0
getNames(a3, dim = 1) <- "treecover"
crop <- mbind(a1, a2, a3)
getNames(crop, dim = 1) <- paste("crop", getNames(crop, dim = 1), sep = "_")
}
} else {
crop <- a[, , "crop"]
}
if ("past" %in% subcategories) {
warning("There are no subcatgories for pasture. Returning total pasture area")
past <- a[, , "past"]
} else {
past <- a[, , "past"]
}
if ("forestry" %in% subcategories) {
# calculate detailed forestry land module carbon stock: aff, ndc, plant
ov32_land <- readGDX(gdx, "ov32_land", select = list(type = "level"))
if (!is.null(ov32_land)) {
p32_carbon_density_ac <- readGDX(gdx, "p32_carbon_density_ac", react = "silent")
if (is.null(p32_carbon_density_ac)) p32_carbon_density_ac <- readGDX(gdx, "pm_carbon_density_ac")
ov32_carbon_stock <- ov32_land * p32_carbon_density_ac
ov32_carbon_stock <- dimSums(ov32_carbon_stock, dim = "ac")
if (!"soilc" %in% getNames(ov32_carbon_stock, dim = "ag_pools")) {
names(dimnames(ov32_carbon_stock))[[3]] <- "land.c_pools"
if (dyn_som) {
ov59_som_pool <- readGDX(gdx, "ov59_som_pool", select = list(type = "level"))
ov_land <- readGDX(gdx, "ov_land", select = list(type = "level"))
ov59_som_pool <- round(ov59_som_pool, 15)
ov_land <- round(ov_land, 15)
top <- ov59_som_pool / ov_land
top[is.na(top)] <- 0
top[is.infinite(top)] <- 0
sub <- readGDX(gdx, "i59_subsoilc_density")[, getYears(top), ]
soilc <- dimSums(ov32_land, dim = 3.2) * (collapseNames(top[, , "forestry"]) + sub)
soilc <- add_dimension(soilc, dim = 3.2, add = "c_pools", nm = "soilc")
} else {
soilc <- dimSums(ov32_land, dim = 3.2) *
collapseNames(readGDX(gdx, "fm_carbon_density")[, getYears(ov32_land), "forestry"])[, , "soilc"]
}
ov32_carbon_stock <- mbind(ov32_carbon_stock, soilc)
}
# check
.testDiff(dimSums(ov32_carbon_stock, dim = 3.1),
collapseNames(a[, , "forestry"]),
1e-3, "Differences in ov32_carbon_stock detected!")
getNames(ov32_carbon_stock, dim = 1) <- paste("forestry", getNames(ov32_carbon_stock, dim = 1), sep = "_")
forestry <- ov32_carbon_stock
}
} else {
forestry <- a[, , "forestry"]
}
if ("primforest" %in% subcategories) {
warning("There are no subcatgories for primforest Returning total primforest area")
primforest <- a[, , "primforest"]
} else {
primforest <- a[, , "primforest"]
}
if ("secdforest" %in% subcategories) {
warning("There are no subcatgories for secdforest. Returning total secdforest area")
secdforest <- a[, , "secdforest"]
} else {
secdforest <- a[, , "secdforest"]
}
if ("urban" %in% subcategories) {
warning("There are no subcatgories for urban land. Returning total urban area")
urban <- a[, , "urban"]
} else {
urban <- a[, , "urban"]
}
if ("other" %in% subcategories) {
# calculate detailed other land carbon stock: othernat and youngsecdf
ov_land_other <- readGDX(gdx, "ov_land_other", select = list(type = "level"), react = "silent")
if (!is.null(ov_land_other)) {
p35_carbon_density_other <- readGDX(gdx, "p35_carbon_density_other")
other_carbon_stock <- ov_land_other * p35_carbon_density_other
other_carbon_stock <- dimSums(other_carbon_stock, dim = "ac")
if (!"soilc" %in% getNames(other_carbon_stock, dim = "ag_pools")) {
names(dimnames(other_carbon_stock))[[3]] <- "land.c_pools"
if (dyn_som) {
ov59_som_pool <- readGDX(gdx, "ov59_som_pool", select = list(type = "level"))
ov_land <- readGDX(gdx, "ov_land", select = list(type = "level"))
ov59_som_pool <- round(ov59_som_pool, 15)
ov_land <- round(ov_land, 15)
top <- ov59_som_pool / ov_land
top[is.na(top)] <- 0
top[is.infinite(top)] <- 0
sub <- readGDX(gdx, "i59_subsoilc_density")[, getYears(top), ]
soilc <- dimSums(ov_land_other, dim = "ac") * (collapseNames(top[, , "other"]) + sub)
soilc <- add_dimension(soilc, dim = 3.2, add = "c_pools", nm = "soilc")
} else {
soilc1 <- dimSums(ov_land_other[, , "othernat"], dim = "ac") *
collapseNames(readGDX(gdx, "fm_carbon_density")[, getYears(other_carbon_stock), "other"])[, , "soilc"]
soilc2 <- dimSums(ov_land_other[, , "youngsecdf"], dim = "ac") *
collapseNames(
readGDX(gdx, "fm_carbon_density")[, getYears(other_carbon_stock), "secdforest"])[, , "soilc"]
soilc <- mbind(soilc1, soilc2)
}
other_carbon_stock <- mbind(other_carbon_stock, soilc)
}
# check
.testDiff(dimSums(other_carbon_stock, dim = 3.1), collapseNames(a[, , "other"]),
1e-3, "Differences in other land carbon stock detected!")
getNames(other_carbon_stock, dim = 1) <- paste("other", getNames(other_carbon_stock, dim = 1), sep = "_")
other <- other_carbon_stock
} else {
a1 <- a[, , "other"]
getNames(a1, dim = 1) <- "othernat"
a2 <- a1
a2[, , ] <- 0
getNames(a2, dim = 1) <- "youngsecdf"
other <- mbind(a1, a2)
getNames(other, dim = 1) <- paste("other", getNames(other, dim = 1), sep = "_")
}
} else {
other <- a[, , "other"]
}
a <- mbind(crop, past, forestry, primforest, secdforest, urban, other)
}
# sum over land pools
if (sum_land) a <- dimSums(a, dim = "land")
# sum over carbon pools
if (sum_cpool) a <- dimSums(a, dim = "c_pools")
# aggregate over regions
if (level != "cell") a <- superAggregate(a, aggr_type = "sum", level = level, na.rm = FALSE)
out(a, file)
}
# nolint end
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.