Dedicated emissions reporting for climate assessment

R/reportEmiForClimateAssessment.R

@@ -0,0 +1,160 @@
+#' Reports emissions & air pollutant values from GDX for climate assessment in between Nash iterations before some of
+#' the energy system variables are defined. The report contains only a subset of reported emissions, with the main
+#' difference being that Emi|CO2|Energy and Industrial Processes is calculated from the difference between total and
+#' LUC emissions. Only global values from this function should be used, as it also skips the subtraction
+#' of certain non-regional sources, such as bunkers, from the regional information
+#'
+#' @param gdx a GDX as created by readGDX, or the file name of a gdx
+#' @param output a magpie object containing all needed variables generated by other report*.R functions
+#' @param regionSubsetList a list containing regions to create report variables region
+#' aggregations. If NULL (default value) only the global region aggregation "GLO" will
+#' be created.
+#' @param t temporal resolution of the reporting, default: t=c(seq(2005,2060,5),seq(2070,2110,10),2130,2150)
+#'
+#' @author Gabriel Abrahao, Tonn Rüter
+#' @examples
+#' \dontrun{
+#' reportEmiForClimateAssessment(gdx)
+#' }
+#' @export
+#'
+#' @importFrom gdx readGDX
+#' @importFrom magclass mbind mselect getItems getYears getSets new.magpie
+reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = NULL,
+                                          t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150)) {
+  # NOTE: This function is a copy of reportEmi with unnecessary parts removed
+  # Read Data from GDX ----
+
+  ####### get realisations #########
+  module2realisation <- readGDX(gdx, "module2realisation")
+  rownames(module2realisation) <- module2realisation$modules
+
+  # unit conversion parameters needed
+  sm_c_2_co2 <- readGDX(gdx, "sm_c_2_co2")
+  GtC_2_MtCO2 <- sm_c_2_co2 * 1000 # conversion of GtC to MtCO2
+  MtN2_to_ktN2O <- 44 / 28 * 1000 # conversion from MtN to ktN2O
+
+  # sets required
+  emiMac2sector <- readGDX(gdx, "emiMac2sector") # mapping of MAC sectors to emissions sectors and gases
+  macSector2emiMkt <- readGDX(gdx, "macSector2emiMkt") # mapping of MAC sectors to emissions markets
+
+  ### emissions variables from REMIND (see definitions in core/equations.gms)
+  # total GHG emissions
+  # total emissions by gas
+  vm_emiAllMkt <- readGDX(gdx, "vm_emiAllMkt", field = "l", restore_zeros = FALSE)[, t, ]
+  # total energy emissions from pe2se and se2fe conversions
+  vm_emiTeDetailMkt <- readGDX(gdx, c("vm_emiTeDetailMkt", "v_emiTeDetailMkt"), field = "l", restore_zeros = FALSE)[, t, ]
+  # total energy emissions in REMIND
+  # Emissions from MACs (currently: all emissions outside of energy CO2 emissions)
+  vm_emiMacSector <- readGDX(gdx, "vm_emiMacSector", field = "l", restore_zeros = FALSE)[, t, ]
+  # exogenous emissions (SO2, BC, OC)
+  # F-Gases
+  vm_emiFgas <- readGDX(gdx, "vm_emiFgas", field = "l", restore_zeros = FALSE)[, t, ]
+
+  #### Markets for total CO2 emissions
+  sel_vm_emiAllMkt_co2 <- if (getSets(vm_emiAllMkt)[[3]] == "emiTe") {
+    mselect(vm_emiAllMkt, emiTe = "co2")
+  } else {
+    mselect(vm_emiAllMkt, all_enty = "co2")
+  }
+
+  out <- mbind(
+    # Basic Total CO2 emissions
+    setNames(dimSums(sel_vm_emiAllMkt_co2, dim = 3) * GtC_2_MtCO2, "Emi|CO2 (Mt CO2/yr)"),
+    # Basic Land-use change CO2 emissions
+    setNames(dimSums(vm_emiMacSector[, , "co2luc"], dim = 3) * GtC_2_MtCO2, "Emi|CO2|+|Land-Use Change (Mt CO2/yr)")
+  )
+  # Basic Get Energy and Industrial Processes from the difference between total and Land-use change emissions, for
+  # climate assessment
+  out <- mbind(
+    out,
+    setNames(
+      out[, , "Emi|CO2 (Mt CO2/yr)"] - out[, , "Emi|CO2|+|Land-Use Change (Mt CO2/yr)"],
+      "Emi|CO2|Energy and Industrial Processes (Mt CO2/yr)"
+    )
+  )
+
+  #  Basic Non-CO2 GHG Emissions
+  # join mac mapping of sectors and markets
+  mac.map <- right_join(emiMac2sector, macSector2emiMkt, by = "all_enty")
+  # create magpie array with MAC emissions per sector, market and gas
+  # MAC emissions comprise non-energy CO2, CH4, N2O emissions
+  emiMAC <- new.magpie(
+    getItems(vm_emiMacSector, dim = "all_regi"),
+    getYears(vm_emiMacSector),
+    paste(
+      mac.map$all_enty,
+      mac.map$emi_sectors,
+      mac.map$all_emiMkt,
+      mac.map[[if ("emiAll" %in% names(mac.map)) "emiAll" else "all_enty1"]],
+      sep = "."
+    )
+  )
+  # rename dimensions for sake of understanding
+  getSets(emiMAC) <- c("region", "year", "macsector", "sector", "emiMkt", "gas")
+  emiMAC[, , mac.map$all_enty] <- vm_emiMacSector[, , mac.map$all_enty]
+
+  ### Basic Energy CH4 and N2O emissions, by markets -------
+  sel_vm_emiTeDetailMkt_ch4 <- if (getSets(vm_emiTeDetailMkt)[[6]] == "emiAll") {
+    mselect(vm_emiTeDetailMkt, emiAll = "ch4")
+  } else {
+    mselect(vm_emiTeDetailMkt, all_enty2 = "ch4")
+  }
+
+  sel_vm_emiTeDetailMkt_n2o <- if (getSets(vm_emiTeDetailMkt)[[6]] == "emiAll") {
+    mselect(vm_emiTeDetailMkt, emiAll = "n2o")
+  } else {
+    mselect(vm_emiTeDetailMkt, all_enty2 = "n2o")
+  }
+
+  out <- mbind(
+    out,
+    # total CH4 emissions
+    setNames(
+      dimSums(mselect(emiMAC, gas = "ch4"), dim = 3) + dimSums(sel_vm_emiTeDetailMkt_ch4, dim = 3),
+      "Emi|CH4 (Mt CH4/yr)"
+    ),
+    # total N2O emissions
+    setNames(
+      (dimSums(mselect(emiMAC, gas = "n2o"), dim = 3) + dimSums(sel_vm_emiTeDetailMkt_n2o, dim = 3)) * MtN2_to_ktN2O,
+      "Emi|N2O (kt N2O/yr)"
+    )
+  )
+
+  out <- mbind(
+    out,
+    ### Basic PFCs
+    setNames(vm_emiFgas[, , "emiFgasCF4"],      "Emi|CF4 (kt CF4/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasC2F6"],     "Emi|C2F6 (kt C2F6/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasC6F14"],    "Emi|C6F14 (kt C6F14/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasHFC"],      "Emi|HFC (kt HFC134a-equiv/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasHFC125"],   "Emi|HFC|HFC125 (kt HFC125/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasHFC134a"],  "Emi|HFC|HFC134a (kt HFC134a/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasHFC143a"],  "Emi|HFC|HFC143a (kt HFC143a/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasHFC227ea"], "Emi|HFC|HFC227ea (kt HFC227ea/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasHFC23"],    "Emi|HFC|HFC23 (kt HFC23/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasHFC245fa"], "Emi|HFC|HFC245fa (kt HFC245fa/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasHFC32"],    "Emi|HFC|HFC32 (kt HFC32/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasHFC43-10"], "Emi|HFC|HFC43-10 (kt HFC43-10/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasPFC"],      "Emi|PFC (kt CF4-equiv/yr)"),
+    setNames(vm_emiFgas[, , "emiFgasSF6"],      "Emi|SF6 (kt SF6/yr)"),
+    ### Basic F-Gases (Mt CO2eq)
+    setNames(vm_emiFgas[, , "emiFgasTotal"],    "Emi|GHG|+|F-Gases (Mt CO2eq/yr)")
+  )
+
+  # Add global values
+  out <- mbind(out, dimSums(out, dim = 1))
+  # add other region aggregations
+  if (!is.null(regionSubsetList)) {
+    out <- mbind(out, calc_regionSubset_sums(out, regionSubsetList))
+  }
+
+  # Run air pollution report
+  message("reportEmiForClimateAssessment executes reportEmiAirPol")
+  pollutants <- reportEmiAirPol(gdx)
+  # Combine both reports
+  out <- mbind(out, pollutants[getItems(out, dim = "all_regi"), getItems(out, dim = "tall"), ])
+
+  getSets(out)[3] <- "variable"
+  return(out)
+}