From e8e709193751723ce46abf2c20f8c3535c48cccf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tonn=20R=C3=BCter?= Date: Fri, 14 Jun 2024 13:55:18 +0200 Subject: [PATCH 1/6] First draft of CO2 emissions reporting --- R/reportEmiForClimateAssessment.R | 248 ++++++++++++++++++++++++++++++ 1 file changed, 248 insertions(+) create mode 100644 R/reportEmiForClimateAssessment.R diff --git a/R/reportEmiForClimateAssessment.R b/R/reportEmiForClimateAssessment.R new file mode 100644 index 00000000..5d7f8d39 --- /dev/null +++ b/R/reportEmiForClimateAssessment.R @@ -0,0 +1,248 @@ +#' Reports emissions & air pollutant values from GDX for climate assessment in between Nash iterations +#' +#' @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, Felix Schreyer, Simón Moreno Leiva +#' @examples +#' \dontrun{ +#' reportEmiForClimateAssessment(gdx) +#' } +#' @export +#' +#' @importFrom gdx readGDX +#' @importFrom magclass mbind mselect mselect<- getItems getRegions getYears getSets dimSums new.magpie + +require(gdx) +require(magclass) +require(dplyr) +source("./R/reportFE.R") +source("./R/reportEmiAirPol.R") + +reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = NULL, + t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150)) { + # emissions calculation requires information from other reporting functions + # if (is.null(output)) { + # message("reportEmi executes reportFE") + # output <- mbind(output, reportFE(gdx, regionSubsetList = regionSubsetList, t = t)) + # } + + # if (is.null(output)) { + # # Error in `.dimextract()`: + # # ! subscript out of bounds ("GLO") + # # Run `rlang::last_trace()` to see where the error occurred. + # # > rlang::last_trace() + # # + # # Error in `.dimextract()`: + # # ! subscript out of bounds ("GLO") + # # --- + # # Backtrace: + # # ▆ + # # 1. └─global reportEmiForClimateAssessment(postsolve_gdx_path) + # # 2. ├─magclass::mbind(...) + # # 3. └─global reportEmiAirPol(gdx, regionSubsetList = regionSubsetList, t = t) + # # 4. ├─methods (local) `[<-`(`*tmp*`, "GLO", , , value = ``) at ./R/reportEmiAirPol.R:114:7 + # # 5. └─magclass (local) `[<-`(`*tmp*`, "GLO", , , value = ``) + # # 6. └─magclass (local) .local(x, i, j, ..., value = value) + # # 7. └─magclass:::.dimextract(x, i, 1, pmatch = pmatch) + # message("reportEmiForClimateAssessment executes reportAirPol") + # output <- mbind(output, reportEmiAirPol(gdx, regionSubsetList = regionSubsetList, t = t)) + # } + + # intialize varibles used in dplyr operations + + # 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 + + # other parameters required + + # switches relevant for emissions reporting + + # 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 + + # additional sets needed + # SE carriers for solids, liquids, gases + + + # SE carriers by origin + + if (is.null(entySEbio <- readGDX(gdx, "entySEbio", react = "silent"))) + entySEbio <- c("sesobio", "seliqbio", "segabio") + + if (is.null(entySEsyn <- readGDX(gdx, "entySEsyn", react = "silent")) || + (length(entySEbio) == length(entySEsyn) && all(entySEbio == entySEsyn))) { + entySEsyn <- c("seliqsyn", "segasyn") + } + + + ### 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) + # total energy emissions in REMIND + # emissions from MAC curves (non-energy 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, ] + # Emissions from MACs (currently: all emissions outside of energy CO2 emissions) + vm_emiMacSector <- readGDX(gdx, "vm_emiMacSector", field = "l", restore_zeros = FALSE) + # energy extraction energy-related CO2 emissions + v_emiEnFuelEx <- readGDX(gdx, "v_emiEnFuelEx", field = "l", restore_zeros = FALSE) + # set to zero if not available + if (is.null(v_emiEnFuelEx)) { + v_emiEnFuelEx <- vm_emiMacSector[, , "co2"] * 0 + } + + # START OF BASICMODE + # if basicmode is TRUE, we need only a subset of reported emissions. + # This is meant to be used for the climate assessment, which has to be run + # inside REMIND before some of the energy system variables are defined + # + # Most of the code here comes from the normal function. The only real difference is + # 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 + # } else { + out <- NULL + + #### 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") + } + + # Basic Total CO2 emissions + out <- mbind( + out, + setNames(dimSums(sel_vm_emiAllMkt_co2, dim = 3) * GtC_2_MtCO2, "Emi|CO2 (Mt CO2/yr)") + ) + # Basic Land-use change CO2 emissions + out <- mbind( + out, + 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( + getRegions(vm_emiMacSector), + 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)" + ) + ) + + ### Basic PFCs ---- + out <- mbind( + out, + 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) + out <- mbind( + out, + # F-Gases + setNames(vm_emiFgas[, , "emiFgasTotal"], + "Emi|GHG|+|F-Gases (Mt CO2eq/yr)" + ) + ) + + # add global values + out <- mbind(out, dimSums(out, dim = 1)) + + # out <- mbind(out, dimSums(out, dim = 1)) + + getSets(out)[3] <- "variable" + return(out) +} + +postsolve_gdx_path <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/fulldata_postsolve.gdx" +emiForCA <- reportEmiForClimateAssessment(postsolve_gdx_path) + + + +emiForCA From db60b5a2fd4c4af467fc180dd9a05b6718d82e7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tonn=20R=C3=BCter?= Date: Tue, 18 Jun 2024 14:24:31 +0200 Subject: [PATCH 2/6] Prototype of seleton reportEmi works, comparison was difficult --- R/reportEmiForClimateAssessment.R | 107 +++++++++++++++++++++++++++--- 1 file changed, 96 insertions(+), 11 deletions(-) diff --git a/R/reportEmiForClimateAssessment.R b/R/reportEmiForClimateAssessment.R index 5d7f8d39..c414f016 100644 --- a/R/reportEmiForClimateAssessment.R +++ b/R/reportEmiForClimateAssessment.R @@ -17,13 +17,10 @@ #' #' @importFrom gdx readGDX #' @importFrom magclass mbind mselect mselect<- getItems getRegions getYears getSets dimSums new.magpie - -require(gdx) -require(magclass) -require(dplyr) -source("./R/reportFE.R") -source("./R/reportEmiAirPol.R") - +library(gdx) +library(magclass) +library(dplyr) +library(remind2) reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = NULL, t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150)) { # emissions calculation requires information from other reporting functions @@ -237,12 +234,100 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = # out <- mbind(out, dimSums(out, dim = 1)) getSets(out)[3] <- "variable" + message("reportEmiForClimateAssessment about to return") return(out) } +# +# REFERENCE data (sans air pollution, right??) +# +library(quitte) +mif_path <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/reportemi_basicmode_fulldata_postsolve.mif" +# ref <- as.quitte(mif_path) +ref <- as.quitte(mif_path) %>% arrange(.data$variable, .data$region, .data$period) +unique(ref$region) # LAM OAS SSA EUR NEU MEA REF CAZ CHA IND JPN USA GLO +unique(ref$variable) +ref + +# +# EMISSIONS 1: Using the reduced version of reportEmi +# postsolve_gdx_path <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/fulldata_postsolve.gdx" -emiForCA <- reportEmiForClimateAssessment(postsolve_gdx_path) - - +# postsolve_gdx_path <- "fulldata_postsolve.gdx" +# emiForCA <- reportEmiForClimateAssessment(postsolve_gdx_path) +# getRegions(emiForCA) +emiForCa <- reportEmiForClimateAssessment("fulldata_postsolve.gdx") + +# For the life of me I cannot figure out why theres region "dummy" instead of "GLO" in +# the data returned to me. Gabriel assured me that when he reads THE SAME FILE, no +# "dummy" region values occur. So work around this for development purposes.. +regions_sans_dummy <- getItems(emiForCa, dim = "all_regi") +regions_sans_dummy[regions_sans_dummy == "dummy"] <- "GLO" +# regions_sans_dummy +getItems(emiForCa, dim = "all_regi") <- regions_sans_dummy + +# In order to be able to compare emiForCa with ref, we need to sort the data! HOWEVER, +# during the as.quitte step, apparently the LEVELS of som factor columns are sorted +# differently in ref than in emiForCa. So we need to sort the levels of the factor +# mifForCa <- mifForCa %>% +foo <- as.quitte(emiForCa) +unique(foo$region) +foo$region[foo$region != "(Missing)"] +levels(foo$region)[-length(levels(foo$region))] +levels(ref$region) +levels(foo$variable) + +mifForCa <- as.quitte(emiForCa) %>% + mutate( + variable = factor(variable, levels = sort(levels(variable))), + region = factor(region, levels = levels(region)[levels(region) != "(Missing)"]), + unit = factor(unit, levels = sort(levels(unit)))) %>% + arrange(.data$variable, .data$region, .data$period) +# mifForCa <- mifForCa %>% arrange(.data$region) %>% arrange(.data$variable, .data$period) +# mifForCa <- as.quitte(emiForCa) %>% arrange(.data$variable, .data$region, .data$period) +# levels(mifForCa$variable) <- +# mifForCa <- as.quitte(emiForCA) +# mifForCa <- mifForCa[order(mifForCa$variable, mifForCa$region, mifForCa$period), ] +# rownames(mifForCa) <- NULL +# rownames(mifForCa) <- 1:nrow(mifForCa) +# mifForCa <- mifForCa %>% arrange(.data$variable, .data$region, .data$period) +# mifForCa %>% arrange(.data$variable) + +unique(mifForCa$region) # LAM OAS SSA EUR NEU MEA REF CAZ CHA IND JPN USA GLO +unique(mifForCa$variable) + +# +# VERIFY the similarity between reference and newly created data sets +# +# Check if the data sets share the same variables, regions, scenarios and models +# as an indication that we're not completely off. The set difference of ref +# and emiForCa should be empty in case they share exactly the same variables etc. +# Gives [1, 5], i.e. all elements in the first set that are not in the second set +setdiff(c(1, 2, 3, 5), c(2, 3, 4)) +variable_difference <- setdiff(unique(ref$region), unique(mifForCa$region)) +region_difference <- setdiff(unique(ref$region), unique(mifForCa$region)) +scenario_difference <- setdiff(unique(ref$scenario), unique(mifForCa$scenario)) +model_difference <- setdiff(unique(ref$model), unique(mifForCa$model)) +# This should show all zeros +c(length(variable_difference), length(region_difference), length(scenario_difference), length(model_difference)) + +# Finally compare ALL VALUES. Should print TRUE in current config +all.equal(mifForCa, ref) + +# +# MBIND to results of reportEmiAirPol. MOVE THIS TO reportEmiForClimateAssessment.R !!! +# +# emiAirPol <- reportEmiAirPol("fulldata_postsolve.gdx") +emiAirPol <- reportEmiAirPol(postsolve_gdx_path) +# Convert to tibble, check for differences in factors.. +mifAirPol <- as.quitte(emiAirPol) +variable_difference <- setdiff(unique(mifForCa$region), unique(mifAirPol$region)) +region_difference <- setdiff(unique(mifForCa$region), unique(mifAirPol$region)) +scenario_difference <- setdiff(unique(mifForCa$scenario), unique(mifAirPol$scenario)) +model_difference <- setdiff(unique(mifForCa$model), unique(mifAirPol$model)) +c(length(variable_difference), length(region_difference), length(scenario_difference), length(model_difference)) + +# Check dimnames(emiForCa) to see which dim are available +# out <- mbind(emiForCa, emiAirPol[getRegions(emiForCA), getYears(emiForCa), ]) +out <- mbind(emiForCa, emiAirPol[getItems(emiForCa, dim = "all_regi"), getItems(emiForCa, dim = "tall"), ]) -emiForCA From f9057bfe948c28915397b59bacd7f656d5a123e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tonn=20R=C3=BCter?= Date: Tue, 18 Jun 2024 14:43:23 +0200 Subject: [PATCH 3/6] Test options for emiAirPol call with regards to mbind --- R/reportEmiForClimateAssessment.R | 62 ++++++++----------------------- 1 file changed, 15 insertions(+), 47 deletions(-) diff --git a/R/reportEmiForClimateAssessment.R b/R/reportEmiForClimateAssessment.R index c414f016..662ac2a6 100644 --- a/R/reportEmiForClimateAssessment.R +++ b/R/reportEmiForClimateAssessment.R @@ -17,77 +17,33 @@ #' #' @importFrom gdx readGDX #' @importFrom magclass mbind mselect mselect<- getItems getRegions getYears getSets dimSums new.magpie -library(gdx) -library(magclass) -library(dplyr) -library(remind2) reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = NULL, t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150)) { - # emissions calculation requires information from other reporting functions - # if (is.null(output)) { - # message("reportEmi executes reportFE") - # output <- mbind(output, reportFE(gdx, regionSubsetList = regionSubsetList, t = t)) - # } - - # if (is.null(output)) { - # # Error in `.dimextract()`: - # # ! subscript out of bounds ("GLO") - # # Run `rlang::last_trace()` to see where the error occurred. - # # > rlang::last_trace() - # # - # # Error in `.dimextract()`: - # # ! subscript out of bounds ("GLO") - # # --- - # # Backtrace: - # # ▆ - # # 1. └─global reportEmiForClimateAssessment(postsolve_gdx_path) - # # 2. ├─magclass::mbind(...) - # # 3. └─global reportEmiAirPol(gdx, regionSubsetList = regionSubsetList, t = t) - # # 4. ├─methods (local) `[<-`(`*tmp*`, "GLO", , , value = ``) at ./R/reportEmiAirPol.R:114:7 - # # 5. └─magclass (local) `[<-`(`*tmp*`, "GLO", , , value = ``) - # # 6. └─magclass (local) .local(x, i, j, ..., value = value) - # # 7. └─magclass:::.dimextract(x, i, 1, pmatch = pmatch) - # message("reportEmiForClimateAssessment executes reportAirPol") - # output <- mbind(output, reportEmiAirPol(gdx, regionSubsetList = regionSubsetList, t = t)) - # } - - # intialize varibles used in dplyr operations - + # 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 - # other parameters required - - # switches relevant for emissions reporting - # 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 - # additional sets needed - # SE carriers for solids, liquids, gases - - # SE carriers by origin - if (is.null(entySEbio <- readGDX(gdx, "entySEbio", react = "silent"))) entySEbio <- c("sesobio", "seliqbio", "segabio") if (is.null(entySEsyn <- readGDX(gdx, "entySEsyn", react = "silent")) || - (length(entySEbio) == length(entySEsyn) && all(entySEbio == entySEsyn))) { + (length(entySEbio) == length(entySEsyn) && all(entySEbio == entySEsyn))) { entySEsyn <- c("seliqsyn", "segasyn") } - ### emissions variables from REMIND (see definitions in core/equations.gms) # total GHG emissions # total emissions by gas @@ -232,7 +188,10 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = out <- mbind(out, dimSums(out, dim = 1)) # out <- mbind(out, dimSums(out, dim = 1)) + # message("reportEmiForClimateAssessment executes reportEmiAirPol") + # emiAirPol <- mbind(output, reportFE(gdx, regionSubsetList = regionSubsetList, t = t)) + out <- mbind(emiForCa, emiAirPol[getItems(emiForCa, dim = "all_regi"), getItems(emiForCa, dim = "tall"), ]) getSets(out)[3] <- "variable" message("reportEmiForClimateAssessment about to return") return(out) @@ -319,6 +278,13 @@ all.equal(mifForCa, ref) # # emiAirPol <- reportEmiAirPol("fulldata_postsolve.gdx") emiAirPol <- reportEmiAirPol(postsolve_gdx_path) +# emiAirPol <- reportEmiAirPol(postsolve_gdx_path, regionSubsetList = c("GLO")) +# emiAirPol <- reportEmiAirPol( +# postsolve_gdx_path, +# regionSubsetList = getItems(emiForCa, dim = "all_regi"), +# t = getItems(emiForCa, dim = "tall") +# ) + # Convert to tibble, check for differences in factors.. mifAirPol <- as.quitte(emiAirPol) variable_difference <- setdiff(unique(mifForCa$region), unique(mifAirPol$region)) @@ -330,4 +296,6 @@ c(length(variable_difference), length(region_difference), length(scenario_differ # Check dimnames(emiForCa) to see which dim are available # out <- mbind(emiForCa, emiAirPol[getRegions(emiForCA), getYears(emiForCa), ]) out <- mbind(emiForCa, emiAirPol[getItems(emiForCa, dim = "all_regi"), getItems(emiForCa, dim = "tall"), ]) - +# For some reason, providing t = getItems(emiForCa, dim = "tall") to reportEmiAirPol still returns years that are +# not in emiForCa. So we still need to subset the years, regions +# out <- mbind(emiForCa, emiAirPol) From 23c52299302ffae0990c20bd5594b107884c3d95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tonn=20R=C3=BCter?= Date: Wed, 19 Jun 2024 15:34:51 +0200 Subject: [PATCH 4/6] Finally figured out the dimSums issue, function works now --- R/reportEmiForClimateAssessment.R | 64 ++++++++++++++++++++----------- 1 file changed, 41 insertions(+), 23 deletions(-) diff --git a/R/reportEmiForClimateAssessment.R b/R/reportEmiForClimateAssessment.R index 662ac2a6..b46e3dcc 100644 --- a/R/reportEmiForClimateAssessment.R +++ b/R/reportEmiForClimateAssessment.R @@ -16,12 +16,19 @@ #' @export #' #' @importFrom gdx readGDX -#' @importFrom magclass mbind mselect mselect<- getItems getRegions getYears getSets dimSums new.magpie +#' @importFrom magclass mbind mselect mselect<- getItems getRegions getYears getSets new.magpie +# library(remind2) +source("R/dimSums.R") +source("R/reportEmiAirPol.R") +library(gdx) +# library(magclass) +library(dplyr) 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 ---- - + gdx <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/fulldata_postsolve.gdx" + t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150) ####### get realisations ######### module2realisation <- readGDX(gdx, "module2realisation") rownames(module2realisation) <- module2realisation$modules @@ -36,13 +43,13 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = macSector2emiMkt <- readGDX(gdx, "macSector2emiMkt") # mapping of MAC sectors to emissions markets # SE carriers by origin - if (is.null(entySEbio <- readGDX(gdx, "entySEbio", react = "silent"))) - entySEbio <- c("sesobio", "seliqbio", "segabio") + # if (is.null(entySEbio <- readGDX(gdx, "entySEbio", react = "silent"))) + # entySEbio <- c("sesobio", "seliqbio", "segabio") - if (is.null(entySEsyn <- readGDX(gdx, "entySEsyn", react = "silent")) || - (length(entySEbio) == length(entySEsyn) && all(entySEbio == entySEsyn))) { - entySEsyn <- c("seliqsyn", "segasyn") - } + # if (is.null(entySEsyn <- readGDX(gdx, "entySEsyn", react = "silent")) || + # (length(entySEbio) == length(entySEsyn) && all(entySEbio == entySEsyn))) { + # entySEsyn <- c("seliqsyn", "segasyn") + # } ### emissions variables from REMIND (see definitions in core/equations.gms) # total GHG emissions @@ -51,19 +58,11 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = # total energy emissions from pe2se and se2fe conversions vm_emiTeDetailMkt <- readGDX(gdx, c("vm_emiTeDetailMkt", "v_emiTeDetailMkt"), field = "l", restore_zeros = FALSE) # total energy emissions in REMIND - # emissions from MAC curves (non-energy emissions) + # 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, ] - # Emissions from MACs (currently: all emissions outside of energy CO2 emissions) - vm_emiMacSector <- readGDX(gdx, "vm_emiMacSector", field = "l", restore_zeros = FALSE) - # energy extraction energy-related CO2 emissions - v_emiEnFuelEx <- readGDX(gdx, "v_emiEnFuelEx", field = "l", restore_zeros = FALSE) - # set to zero if not available - if (is.null(v_emiEnFuelEx)) { - v_emiEnFuelEx <- vm_emiMacSector[, , "co2"] * 0 - } # START OF BASICMODE # if basicmode is TRUE, we need only a subset of reported emissions. @@ -81,8 +80,10 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = #### Markets for total CO2 emissions sel_vm_emiAllMkt_co2 <- if (getSets(vm_emiAllMkt)[[3]] == "emiTe") { + print("emiTe") mselect(vm_emiAllMkt, emiTe = "co2") } else { + print("all_enty") mselect(vm_emiAllMkt, all_enty = "co2") } @@ -174,7 +175,6 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = setNames(vm_emiFgas[, , "emiFgasPFC"], "Emi|PFC (kt CF4-equiv/yr)"), setNames(vm_emiFgas[, , "emiFgasSF6"], "Emi|SF6 (kt SF6/yr)") ) - ### Basic F-Gases (Mt CO2eq) out <- mbind( out, @@ -184,19 +184,33 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = ) ) - # add global values + # 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 # out <- mbind(out, dimSums(out, dim = 1)) - # message("reportEmiForClimateAssessment executes reportEmiAirPol") - # emiAirPol <- mbind(output, reportFE(gdx, regionSubsetList = regionSubsetList, t = t)) + message("reportEmiForClimateAssessment executes reportEmiAirPol") + pollutants <- reportEmiAirPol(postsolve_gdx_path) + - out <- mbind(emiForCa, emiAirPol[getItems(emiForCa, dim = "all_regi"), getItems(emiForCa, dim = "tall"), ]) + # Combine both reports + out <- mbind(out, pollutants[getItems(out, dim = "all_regi"), getItems(out, dim = "tall"), ]) + getSets(out)[3] <- "variable" message("reportEmiForClimateAssessment about to return") return(out) } +library(quitte) +postsolve_gdx_path <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/fulldata_postsolve.gdx" +emiForCa <- reportEmiForClimateAssessment(postsolve_gdx_path) +emiForCa +foo <- as.quitte(emiForCa) + # # REFERENCE data (sans air pollution, right??) # @@ -216,6 +230,8 @@ postsolve_gdx_path <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/full # emiForCA <- reportEmiForClimateAssessment(postsolve_gdx_path) # getRegions(emiForCA) emiForCa <- reportEmiForClimateAssessment("fulldata_postsolve.gdx") +# Check those values, they seem of. E.g. +# emiForCa["GLO", , "Emi|CO2 (Mt CO2/yr)"] # For the life of me I cannot figure out why theres region "dummy" instead of "GLO" in # the data returned to me. Gabriel assured me that when he reads THE SAME FILE, no @@ -262,7 +278,7 @@ unique(mifForCa$variable) # as an indication that we're not completely off. The set difference of ref # and emiForCa should be empty in case they share exactly the same variables etc. # Gives [1, 5], i.e. all elements in the first set that are not in the second set -setdiff(c(1, 2, 3, 5), c(2, 3, 4)) +setdiff(c(1, 2, 3, 5), c(2, 3, 4)) variable_difference <- setdiff(unique(ref$region), unique(mifForCa$region)) region_difference <- setdiff(unique(ref$region), unique(mifForCa$region)) scenario_difference <- setdiff(unique(ref$scenario), unique(mifForCa$scenario)) @@ -276,6 +292,7 @@ all.equal(mifForCa, ref) # # MBIND to results of reportEmiAirPol. MOVE THIS TO reportEmiForClimateAssessment.R !!! # +# No matter if regionSubsetList or t is provided, the size and runtimes of mifAirPol stays the same # emiAirPol <- reportEmiAirPol("fulldata_postsolve.gdx") emiAirPol <- reportEmiAirPol(postsolve_gdx_path) # emiAirPol <- reportEmiAirPol(postsolve_gdx_path, regionSubsetList = c("GLO")) @@ -299,3 +316,4 @@ out <- mbind(emiForCa, emiAirPol[getItems(emiForCa, dim = "all_regi"), getItems( # For some reason, providing t = getItems(emiForCa, dim = "tall") to reportEmiAirPol still returns years that are # not in emiForCa. So we still need to subset the years, regions # out <- mbind(emiForCa, emiAirPol) +mifOut <- as.quitte(out) \ No newline at end of file From 364cfa71af7f3384d99849016d5ed18f57164164 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tonn=20R=C3=BCter?= Date: Wed, 19 Jun 2024 17:07:13 +0200 Subject: [PATCH 5/6] Cleaned up reportEmiForClimateAssessment --- R/reportEmiForClimateAssessment.R | 78 +++++++++---------------------- 1 file changed, 23 insertions(+), 55 deletions(-) diff --git a/R/reportEmiForClimateAssessment.R b/R/reportEmiForClimateAssessment.R index b46e3dcc..9196e048 100644 --- a/R/reportEmiForClimateAssessment.R +++ b/R/reportEmiForClimateAssessment.R @@ -1,4 +1,15 @@ #' Reports emissions & air pollutant values from GDX for climate assessment in between Nash iterations +#' # START OF BASICMODE +# if basicmode is TRUE, we need only a subset of reported emissions. +# This is meant to be used for the climate assessment, which has to be run +# inside REMIND before some of the energy system variables are defined +# +# Most of the code here comes from the normal function. The only real difference is +# 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 @@ -16,7 +27,7 @@ #' @export #' #' @importFrom gdx readGDX -#' @importFrom magclass mbind mselect mselect<- getItems getRegions getYears getSets new.magpie +#' @importFrom magclass mbind mselect mselect<- getItems getYears getSets new.magpie # library(remind2) source("R/dimSums.R") source("R/reportEmiAirPol.R") @@ -27,8 +38,7 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = 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 ---- - gdx <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/fulldata_postsolve.gdx" - t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150) + ####### get realisations ######### module2realisation <- readGDX(gdx, "module2realisation") rownames(module2realisation) <- module2realisation$modules @@ -42,21 +52,12 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = emiMac2sector <- readGDX(gdx, "emiMac2sector") # mapping of MAC sectors to emissions sectors and gases macSector2emiMkt <- readGDX(gdx, "macSector2emiMkt") # mapping of MAC sectors to emissions markets - # SE carriers by origin - # if (is.null(entySEbio <- readGDX(gdx, "entySEbio", react = "silent"))) - # entySEbio <- c("sesobio", "seliqbio", "segabio") - - # if (is.null(entySEsyn <- readGDX(gdx, "entySEsyn", react = "silent")) || - # (length(entySEbio) == length(entySEsyn) && all(entySEbio == entySEsyn))) { - # entySEsyn <- c("seliqsyn", "segasyn") - # } - ### 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) + 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, ] @@ -64,40 +65,19 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = # F-Gases vm_emiFgas <- readGDX(gdx, "vm_emiFgas", field = "l", restore_zeros = FALSE)[, t, ] - # START OF BASICMODE - # if basicmode is TRUE, we need only a subset of reported emissions. - # This is meant to be used for the climate assessment, which has to be run - # inside REMIND before some of the energy system variables are defined - # - # Most of the code here comes from the normal function. The only real difference is - # 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 - # } else { - out <- NULL - #### Markets for total CO2 emissions sel_vm_emiAllMkt_co2 <- if (getSets(vm_emiAllMkt)[[3]] == "emiTe") { - print("emiTe") mselect(vm_emiAllMkt, emiTe = "co2") } else { - print("all_enty") mselect(vm_emiAllMkt, all_enty = "co2") } - # Basic Total CO2 emissions out <- mbind( - out, - setNames(dimSums(sel_vm_emiAllMkt_co2, dim = 3) * GtC_2_MtCO2, "Emi|CO2 (Mt CO2/yr)") - ) - # Basic Land-use change CO2 emissions - out <- mbind( - out, + # 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( @@ -109,14 +89,12 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = ) # 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( - getRegions(vm_emiMacSector), + getItems(vm_emiMacSector, dim = "all_regi"), getYears(vm_emiMacSector), paste( mac.map$all_enty, @@ -157,9 +135,9 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = ) ) - ### Basic PFCs ---- 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)"), @@ -173,15 +151,9 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = 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) - out <- mbind( - out, - # F-Gases - setNames(vm_emiFgas[, , "emiFgasTotal"], - "Emi|GHG|+|F-Gases (Mt CO2eq/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 @@ -192,16 +164,12 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = } # Run air pollution report - # out <- mbind(out, dimSums(out, dim = 1)) message("reportEmiForClimateAssessment executes reportEmiAirPol") pollutants <- reportEmiAirPol(postsolve_gdx_path) - - # Combine both reports out <- mbind(out, pollutants[getItems(out, dim = "all_regi"), getItems(out, dim = "tall"), ]) - + getSets(out)[3] <- "variable" - message("reportEmiForClimateAssessment about to return") return(out) } From fee3d27c807a04a34d63395c6ef69fa608bf4f99 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tonn=20R=C3=BCter?= Date: Wed, 19 Jun 2024 17:17:10 +0200 Subject: [PATCH 6/6] Removed dev code, added documentation --- R/reportEmiForClimateAssessment.R | 149 +++--------------------------- 1 file changed, 11 insertions(+), 138 deletions(-) diff --git a/R/reportEmiForClimateAssessment.R b/R/reportEmiForClimateAssessment.R index 9196e048..a4ccf60c 100644 --- a/R/reportEmiForClimateAssessment.R +++ b/R/reportEmiForClimateAssessment.R @@ -1,25 +1,17 @@ -#' Reports emissions & air pollutant values from GDX for climate assessment in between Nash iterations -#' # START OF BASICMODE -# if basicmode is TRUE, we need only a subset of reported emissions. -# This is meant to be used for the climate assessment, which has to be run -# inside REMIND before some of the energy system variables are defined -# -# Most of the code here comes from the normal function. The only real difference is -# 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 +#' 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) +#' @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, Felix Schreyer, Simón Moreno Leiva +#' @author Gabriel Abrahao, Tonn Rüter #' @examples #' \dontrun{ #' reportEmiForClimateAssessment(gdx) @@ -27,13 +19,7 @@ #' @export #' #' @importFrom gdx readGDX -#' @importFrom magclass mbind mselect mselect<- getItems getYears getSets new.magpie -# library(remind2) -source("R/dimSums.R") -source("R/reportEmiAirPol.R") -library(gdx) -# library(magclass) -library(dplyr) +#' @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 @@ -137,7 +123,7 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = out <- mbind( out, - ### Basic PFCs + ### 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)"), @@ -165,123 +151,10 @@ reportEmiForClimateAssessment <- function(gdx, output = NULL, regionSubsetList = # Run air pollution report message("reportEmiForClimateAssessment executes reportEmiAirPol") - pollutants <- reportEmiAirPol(postsolve_gdx_path) + 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) -} - -library(quitte) -postsolve_gdx_path <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/fulldata_postsolve.gdx" -emiForCa <- reportEmiForClimateAssessment(postsolve_gdx_path) -emiForCa -foo <- as.quitte(emiForCa) - -# -# REFERENCE data (sans air pollution, right??) -# -library(quitte) -mif_path <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/reportemi_basicmode_fulldata_postsolve.mif" -# ref <- as.quitte(mif_path) -ref <- as.quitte(mif_path) %>% arrange(.data$variable, .data$region, .data$period) -unique(ref$region) # LAM OAS SSA EUR NEU MEA REF CAZ CHA IND JPN USA GLO -unique(ref$variable) -ref - -# -# EMISSIONS 1: Using the reduced version of reportEmi -# -postsolve_gdx_path <- "/p/projects/piam/abrahao/tonns_stuff/reportemi_stuff/fulldata_postsolve.gdx" -# postsolve_gdx_path <- "fulldata_postsolve.gdx" -# emiForCA <- reportEmiForClimateAssessment(postsolve_gdx_path) -# getRegions(emiForCA) -emiForCa <- reportEmiForClimateAssessment("fulldata_postsolve.gdx") -# Check those values, they seem of. E.g. -# emiForCa["GLO", , "Emi|CO2 (Mt CO2/yr)"] - -# For the life of me I cannot figure out why theres region "dummy" instead of "GLO" in -# the data returned to me. Gabriel assured me that when he reads THE SAME FILE, no -# "dummy" region values occur. So work around this for development purposes.. -regions_sans_dummy <- getItems(emiForCa, dim = "all_regi") -regions_sans_dummy[regions_sans_dummy == "dummy"] <- "GLO" -# regions_sans_dummy -getItems(emiForCa, dim = "all_regi") <- regions_sans_dummy - -# In order to be able to compare emiForCa with ref, we need to sort the data! HOWEVER, -# during the as.quitte step, apparently the LEVELS of som factor columns are sorted -# differently in ref than in emiForCa. So we need to sort the levels of the factor -# mifForCa <- mifForCa %>% -foo <- as.quitte(emiForCa) -unique(foo$region) -foo$region[foo$region != "(Missing)"] -levels(foo$region)[-length(levels(foo$region))] -levels(ref$region) -levels(foo$variable) - -mifForCa <- as.quitte(emiForCa) %>% - mutate( - variable = factor(variable, levels = sort(levels(variable))), - region = factor(region, levels = levels(region)[levels(region) != "(Missing)"]), - unit = factor(unit, levels = sort(levels(unit)))) %>% - arrange(.data$variable, .data$region, .data$period) -# mifForCa <- mifForCa %>% arrange(.data$region) %>% arrange(.data$variable, .data$period) -# mifForCa <- as.quitte(emiForCa) %>% arrange(.data$variable, .data$region, .data$period) -# levels(mifForCa$variable) <- -# mifForCa <- as.quitte(emiForCA) -# mifForCa <- mifForCa[order(mifForCa$variable, mifForCa$region, mifForCa$period), ] -# rownames(mifForCa) <- NULL -# rownames(mifForCa) <- 1:nrow(mifForCa) -# mifForCa <- mifForCa %>% arrange(.data$variable, .data$region, .data$period) -# mifForCa %>% arrange(.data$variable) - -unique(mifForCa$region) # LAM OAS SSA EUR NEU MEA REF CAZ CHA IND JPN USA GLO -unique(mifForCa$variable) - -# -# VERIFY the similarity between reference and newly created data sets -# -# Check if the data sets share the same variables, regions, scenarios and models -# as an indication that we're not completely off. The set difference of ref -# and emiForCa should be empty in case they share exactly the same variables etc. -# Gives [1, 5], i.e. all elements in the first set that are not in the second set -setdiff(c(1, 2, 3, 5), c(2, 3, 4)) -variable_difference <- setdiff(unique(ref$region), unique(mifForCa$region)) -region_difference <- setdiff(unique(ref$region), unique(mifForCa$region)) -scenario_difference <- setdiff(unique(ref$scenario), unique(mifForCa$scenario)) -model_difference <- setdiff(unique(ref$model), unique(mifForCa$model)) -# This should show all zeros -c(length(variable_difference), length(region_difference), length(scenario_difference), length(model_difference)) - -# Finally compare ALL VALUES. Should print TRUE in current config -all.equal(mifForCa, ref) - -# -# MBIND to results of reportEmiAirPol. MOVE THIS TO reportEmiForClimateAssessment.R !!! -# -# No matter if regionSubsetList or t is provided, the size and runtimes of mifAirPol stays the same -# emiAirPol <- reportEmiAirPol("fulldata_postsolve.gdx") -emiAirPol <- reportEmiAirPol(postsolve_gdx_path) -# emiAirPol <- reportEmiAirPol(postsolve_gdx_path, regionSubsetList = c("GLO")) -# emiAirPol <- reportEmiAirPol( -# postsolve_gdx_path, -# regionSubsetList = getItems(emiForCa, dim = "all_regi"), -# t = getItems(emiForCa, dim = "tall") -# ) - -# Convert to tibble, check for differences in factors.. -mifAirPol <- as.quitte(emiAirPol) -variable_difference <- setdiff(unique(mifForCa$region), unique(mifAirPol$region)) -region_difference <- setdiff(unique(mifForCa$region), unique(mifAirPol$region)) -scenario_difference <- setdiff(unique(mifForCa$scenario), unique(mifAirPol$scenario)) -model_difference <- setdiff(unique(mifForCa$model), unique(mifAirPol$model)) -c(length(variable_difference), length(region_difference), length(scenario_difference), length(model_difference)) - -# Check dimnames(emiForCa) to see which dim are available -# out <- mbind(emiForCa, emiAirPol[getRegions(emiForCA), getYears(emiForCa), ]) -out <- mbind(emiForCa, emiAirPol[getItems(emiForCa, dim = "all_regi"), getItems(emiForCa, dim = "tall"), ]) -# For some reason, providing t = getItems(emiForCa, dim = "tall") to reportEmiAirPol still returns years that are -# not in emiForCa. So we still need to subset the years, regions -# out <- mbind(emiForCa, emiAirPol) -mifOut <- as.quitte(out) \ No newline at end of file +} \ No newline at end of file