REMIND/EDGE-T variable harmonization

.buildlibrary

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

.pre-commit-config.yaml

@@ -15,7 +15,7 @@ repos:
     -   id: mixed-line-ending
 
 -   repo: https://github.com/lorenzwalthert/precommit
-    rev: 7910e0323d7213f34275a7a562b9ef0fde8ce1b9  # frozen: v0.4.2
+    rev: bae853d82da476eee0e0a57960ee6b741a3b3fb7  # frozen: v0.4.3
     hooks:
     -   id: parsable-R
     -   id: deps-in-desc

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: 'reporttransport: Reporting package for edgeTransport'
-version: 0.0.12
-date-released: '2024-07-30'
+version: 0.1.0
+date-released: '2024-08-23'
 abstract: This package contains edgeTransport-specific routines to report model results.
   The main functionality is to generate transport reporting variables in MIF format
   from a given edgeTransport model run folder or REMIND input data.

DESCRIPTION

@@ -1,8 +1,8 @@
 Type: Package
 Package: reporttransport
 Title: Reporting package for edgeTransport
-Version: 0.0.12
-Date: 2024-07-30
+Version: 0.1.0
+Date: 2024-08-23
 Authors@R:
     person("Johanna", "Hoppe", , "[email protected]", role = c("aut", "cre"))
 Description: This package contains edgeTransport-specific routines to
@@ -21,7 +21,8 @@ Imports:
     quitte,
     remind2,
     rmndt,
-    utils
+    utils, 
+    tibble
 Suggests:
     knitr,
     sf,

NAMESPACE

@@ -3,6 +3,7 @@
 export(aggregateVariables)
 export(checkForNAsAndDups)
 export(convertToMIF)
+export(harmonizeOutput)
 export(prepareForREMIND)
 export(reportAggregatedCosts)
 export(reportAnalyticsVarSet)
@@ -31,10 +32,12 @@ importFrom(gdxrrw,rgdx.param)
 importFrom(gdxrrw,rgdx.scalar)
 importFrom(quitte,aggregate_map)
 importFrom(quitte,as.quitte)
+importFrom(quitte,read.quitte)
 importFrom(quitte,write.mif)
 importFrom(remind2,toolRegionSubsets)
 importFrom(rmndt,approx_dt)
 importFrom(rmndt,disaggregate_dt)
 importFrom(rmndt,magpie2dt)
+importFrom(tibble,tribble)
 importFrom(utils,packageVersion)
 importFrom(utils,write.csv)

R/convertToMIF.R

@@ -139,8 +139,6 @@ convertToMIF <- function(vars, GDPMER, helpers, scenario, model, gdx,  isTranspo
   toMIF[, model := model][, scenario := scenario]
   toMIF <- as.quitte(toMIF)
 
-  toMIF <- renameDuplicateVariables(toMIF)
-
   if (anyNA(toMIF)) stop("MIF output contains NAs.
                          Please check reportAndAggregatedMIF()")
   if (anyDuplicated(toMIF)) stop("MIF output contains duplicates.

R/harmonizeOutput.R

@@ -0,0 +1,104 @@
+#' Harmonize data to last REMIND iteration in coupled mode
+#'
+#' The two models are coupled via energy service demand, energy intensity and capital expenditure on CES node level.
+#' As REMIND converges the output between EDGE-T and REMIND still slightly differs.
+#' This function calculates harmonization factors that are applied on the EDGE-T energy intensity and final energy values.
+#' Furthermore, the harmonization factors indicate how well the two models converged and can be used for analysis purposes.
+#' @param REMINDoutput Output of a REMIND wihtoutPlus.mif to harmonize with
+#' @param edgetOutputDir edgeTransport output directory
+#' @param baseVarSet calculated base variable set including energy intensity and final energy on fleet level
+#' @param data list of raw model output data
+#'
+#' @returns The function returns a list with the harmonized energy intensity and final energy data
+#' @author Johanna Hoppe
+#' @importFrom quitte read.quitte
+#' @importFrom tibble tribble
+#' @import data.table
+#' @export
+
+harmonizeOutput <- function(REMINDoutput, edgetOutputDir, baseVarSet, data) {
+
+    ## Calculate FE harmonization factors that are then applied to the energy intensity
+    harmMap <- tribble(
+      ~variable,                                                  ~ all_teEs,
+      "FE|Transport|Freight|Short-Medium distance|Electricity",    "te_eselt_frgt_sm",
+      "FE|Transport|Freight|Short-Medium distance|Liquids",        "te_esdie_frgt_sm",
+      "FE|Transport|Freight|Short-Medium distance|Hydrogen",       "te_esh2t_frgt_sm",
+      "FE|Transport|Freight|Short-Medium distance|Gases",          "te_esgat_frgt_sm",
+      "FE|Transport|Freight|Long distance|Diesel Liquids",         "te_esdie_frgt_lo",
+      "FE|Transport|Pass|Short-Medium distance|Electricity",       "te_eselt_pass_sm",
+      "FE|Transport|Pass|Short-Medium distance|Diesel Liquids",    "te_esdie_pass_sm",
+      "FE|Transport|Pass|Short-Medium distance|Petrol Liquids",    "te_espet_pass_sm",
+      "FE|Transport|Pass|Short-Medium distance|Hydrogen",          "te_esh2t_pass_sm",
+      "FE|Transport|Pass|Short-Medium distance|Gases",             "te_esgat_pass_sm",
+      "FE|Transport|Pass|Long distance|Diesel Liquids",            "te_esdie_pass_lo"
+    )
+    REMINDoutputToHarmonize <- REMINDoutput[variable %in% unique(harmMap$variable)]
+    if (!length(unique(REMINDoutputToHarmonize$variable)) == length(unique(harmMap$variable))) {
+      stop("Variables for EDGE-T harmonization are missing in the REMIND output. Please check reportEdgeTransport()")
+    }
+    demByTechMap <- unique(data$helpers$mapEdgeToREMIND[, c("all_teEs", "univocalName", "technology")])
+    demByTechMap <- demByTechMap[!is.na(all_teEs)]
+    demByTechMap[technology %in% c("BEV", "Electric"), technology := "Electricity"]
+    demByTechMap[technology == "FCEV", technology := "Hydrogen"]
+    EDGEToutputToHarmonize <- merge(baseVarSet$ext$fleetFEdemand[period %in% timeResReporting], demByTechMap, by = c("univocalName", "technology"), all.x = TRUE) # nolint: object_name_linter
+    EDGEToutputToHarmonize <- EDGEToutputToHarmonize[period %in% timeResReporting & !subsectorL1 %in% c("Walk", "Cycle")]
+    EDGEToutputToHarmonize[, variable := NULL][, unit := NULL]
+    EDGEToutputToHarmonize <- merge(EDGEToutputToHarmonize, harmMap, by = "all_teEs", all.x = TRUE)
+    EDGEToutputToHarmonize <- EDGEToutputToHarmonize[, .(value = sum(value)), by = c("region", "period", "variable", "all_teEs")]
+    setnames(EDGEToutputToHarmonize, "value", "edget")
+    setnames(REMINDoutputToHarmonize, "value", "remind")
+    harmFactors <- merge(EDGEToutputToHarmonize, REMINDoutputToHarmonize, by = c("region", "period", "variable"), all.x = TRUE)
+    harmFactors[, factor := ifelse(!remind < 1e-6 & !edget < 1e-6, remind/edget, 1)][, c("edget", "remind") := NULL]
+    harmFactors[, variable := paste0("Harmonization factor|", variable)][, unit := "-"][, scenario:= data$scenarioName][, model:= data$modelName]
+    storeData(edgetOutputDir, list(harmFactors = harmFactors))
+
+    #Hybrids are listed in energy intensity but are attributed partly to electricity and petrol in the reporting of final energy
+    #The liquids part of the final energy needs to be scaled differntly than the electricity part. This is not possible using the current approach.
+    #Hence, the energy intensity of hybrids is not scaled and the FE is scaled in postprocessing. This leads to a small inconsistency in the dataset
+    #of eneryg service demand, energy intensity and final energy
+
+    #Apply factors
+    harmonizedEnergyIntensity <- merge(baseVarSet$int$fleetEnergyIntensity[period %in% timeResReporting], 
+      data$helpers$mapEdgeToREMIND[, c("all_teEs", "univocalName", "technology")], by = c("univocalName", "technology"), all.x = TRUE)
+    harmonizedEnergyIntensity <- merge(harmonizedEnergyIntensity, harmFactors[, c("region", "period", "all_teEs", "factor")], 
+      by = c("region", "period", "all_teEs"), all.x = TRUE)
+    harmonizedEnergyIntensity[technology == "Hybrid electric", factor := 1]
+    harmonizedEnergyIntensity[, value := ifelse(!subsectorL1 %in% c("Walk", "Cycle"), value * factor, value)][, c("factor", "all_teEs") := NULL]
+    if (anyNA(harmonizedEnergyIntensity)) stop("Variable harmonization did not work. Please check harmonizeOutput() in reporttransport.")
+
+    ## Harmonize final energy
+    harmonizedFinalEnergy <- merge(baseVarSet$ext$fleetFEdemand[period %in% timeResReporting], 
+      demByTechMap, by = c("univocalName", "technology"), all.x = TRUE)
+    harmonizedFinalEnergy <- merge(harmonizedFinalEnergy, harmFactors[, c("region", "period", "all_teEs", "factor")], 
+      by = c("region", "period", "all_teEs"), all.x = TRUE)
+    harmonizedFinalEnergy[, value := ifelse(!subsectorL1 %in% c("Walk", "Cycle"), value * factor, value)]
+    FEcheck <- harmonizedFinalEnergy[!subsectorL1 %in% c("Walk", "Cycle")]
+    harmonizedFinalEnergy[, c("factor", "all_teEs") := NULL]    
+    if (anyNA(harmonizedFinalEnergy)) stop("Variable harmonization did not work. Please check harmonizeOutput() in reporttransport.") 
+
+    ## Check if harmonization worked and how big the deviation due to the hybrids is
+    # test harmonization
+    FEcheck <- FEcheck[, .(value = sum(value)), by = c("region", "period", "all_teEs")]
+    FEcheck <- merge(FEcheck, harmMap, by = "all_teEs", all.x = TRUE)
+    FEcheck <- merge(FEcheck, REMINDoutputToHarmonize, by = c("region", "period", "variable"), all.x = TRUE)
+    FEcheck[, diff := ifelse(!remind < 1e-6 & !value < 1e-6, (value - remind) / value, 0)]
+    if (max(FEcheck$diff > 1e-5)) stop("FE harmonization did not work. Please check harmonizeOutput()")
+
+    # test the deviation (as if the dataset would be consistent)
+    calculatedFinalEnergy <- reportFinalEnergy(harmonizedEnergyIntensity, baseVarSet$ext$fleetESdemand, data$scenSpecLoadFactor, data$hybridElecShare, data$helpers)
+    test <- calculatedFinalEnergy[period %in% timeResReporting & !subsectorL1 %in% c("Walk", "Cycle")]
+    test[, variable := NULL][, unit := NULL]
+    test <- merge(test, demByTechMap, by = c("univocalName", "technology"), all.x = TRUE)
+    test <- merge(test, harmMap, by = "all_teEs", all.x = TRUE)
+    test <- test[, .(value = sum(value)), by = c("region", "period", "variable")]
+    test <- merge(test, REMINDoutputToHarmonize, by = c("region", "period", "variable"), all.x = TRUE)
+    test[, diff := ifelse(!remind < 1e-6 & !value < 1e-6, (value - remind) / value, 0)] 
+
+    harmonizedVars <- list(
+      harmonizedEnergyIntensity = harmonizedEnergyIntensity,
+      harmonizedFinalEnergy = harmonizedFinalEnergy
+    )
+
+    return(harmonizedVars)
+}

R/reportEdgeTransport.R

@@ -30,17 +30,19 @@
 #' @param isAnalyticsReported Switch for activating reporting of model analytics data
 #' @param isREMINDinputReported Switch for activating reporting of REMIND input data
 #' @param isStored Switch for activating data storage and creating the transport.MIF file
+#' @param isHarmonized Switch for activating ouput harmonization with REMIND
 #'
 #' @returns The function either returns the REMINDinputData if isREMINDinputReported is
 #'          enabled or the transport data in MIF format
 #' @author Johanna Hoppe
 #' @importFrom quitte write.mif
+#' @importFrom tibble tribble
 #' @import data.table
 #' @export
 
 reportEdgeTransport <- function(folderPath = file.path(".", "EDGE-T"), data = NULL, isTransportReported = TRUE,
                                 isTransportExtendedReported = FALSE, isAnalyticsReported = FALSE,
-                                isREMINDinputReported = FALSE, isStored = TRUE, ...) {
+                                isREMINDinputReported = FALSE, isStored = TRUE, isHarmonized = FALSE, ...) {
 
   # If you want to change timeResReporting to timesteps outside the modeleled timesteps,
   # please add an interpolation step
@@ -115,11 +117,19 @@ reportEdgeTransport <- function(folderPath = file.path(".", "EDGE-T"), data = NU
     if (is.null(data$scenSpecPrefTrends)) data$scenSpecPrefTrends <- readRDS(file.path(folderPath, "2_InputDataPolicy", "scenSpecPrefTrends.RDS"))
     if (is.null(data$initialIncoCosts)) data$initialIncoCosts <- readRDS(file.path(folderPath, "2_InputDataPolicy", "initialIncoCosts.RDS"))
   }
+
   #########################################################################
   ## Report output variables
   #########################################################################
   # Base variable set that is needed to report REMIND input data and additional detailed transport data
   baseVarSet <- reportBaseVarSet(data = data, timeResReporting = timeResReporting)
+  if (isHarmonized) {
+    REMINDoutput <- as.data.table(read.quitte(data$remindReportingFile))
+    harmonizedVars <- harmonizeOutput(REMINDoutput, edgetOutputDir, baseVarSet, data)
+    baseVarSet$int$fleetEnergyIntensity <- harmonizedVars$harmonizedEnergyIntensity
+    baseVarSet$ext$fleetFEdemand <- harmonizedVars$harmonizedFinalEnergy
+  }
+
   reporting <- baseVarSet
   outputVars <- baseVarSet
 
@@ -158,6 +168,7 @@ reportEdgeTransport <- function(folderPath = file.path(".", "EDGE-T"), data = NU
                               gdx                         = data$gdxPath,
                               isTransportExtendedReported = isTransportExtendedReported)
 
+    if (isHarmonized) reporting <- reporting[!variable %in% REMINDoutput$variable]
     if (isStored) write.mif(reporting, file.path(folderPath, "Transport.mif"))
   }
 

R/storeData.R

@@ -44,8 +44,10 @@ storeData <- function(outputFolder, varsList = NULL, ...) {
                        "fleetVehNumbersIterations",
                        "endogenousCostsIterations",
                        "endogenousCosts",
+                       "sectorESdemand",
                        "ESdemandFVsalesLevel",
-                       "fleetVehiclesPerTech")) subfolder <- "4_Output"
+                       "fleetVehiclesPerTech",
+                       "harmFactors")) subfolder <- "4_Output"
     if (varName %in% c("REMINDinputData")) subfolder <- "5_REMINDinputData"
     if (is.null(subfolder)) stop(paste0("No subfolder assigned to ", varName))
 
@@ -59,7 +61,7 @@ storeData <- function(outputFolder, varsList = NULL, ...) {
 
   storeCSV <- function(varName, vars, outputFolder, subfolder = NULL) {
     if (is.null(subfolder)) subfolder <- allocateFile(varName)
-    write.csv(vars[[varName]], file.path(outputFolder, subfolder, paste0(varName, ".csv")))
+    write.csv(vars[[varName]], file.path(outputFolder, subfolder, paste0(varName, ".csv")), row.names = FALSE)
   }
   vars <- list()
   if (!is.null(varsList)) vars <- varsList

README.md

@@ -1,6 +1,6 @@
 # Reporting package for edgeTransport
 
-R package **reporttransport**, version **0.0.12**
+R package **reporttransport**, version **0.1.0**
 
 [![CRAN status](https://www.r-pkg.org/badges/version/reporttransport)](https://cran.r-project.org/package=reporttransport)  [![R build status](https://github.com/pik-piam/reporttransport/workflows/check/badge.svg)](https://github.com/pik-piam/reporttransport/actions) [![codecov](https://codecov.io/gh/pik-piam/reporttransport/branch/master/graph/badge.svg)](https://app.codecov.io/gh/pik-piam/reporttransport) [![r-universe](https://pik-piam.r-universe.dev/badges/reporttransport)](https://pik-piam.r-universe.dev/builds)
 
@@ -41,7 +41,7 @@ In case of questions / problems please contact Johanna Hoppe <johanna.hoppe@pik-
 
 To cite package **reporttransport** in publications use:
 
-Hoppe J (2024). _reporttransport: Reporting package for edgeTransport_. R package version 0.0.12, <https://github.com/pik-piam/reporttransport>.
+Hoppe J (2024). _reporttransport: Reporting package for edgeTransport_. R package version 0.1.0, <https://github.com/pik-piam/reporttransport>.
 
 A BibTeX entry for LaTeX users is
 
@@ -50,7 +50,7 @@ A BibTeX entry for LaTeX users is
   title = {reporttransport: Reporting package for edgeTransport},
   author = {Johanna Hoppe},
   year = {2024},
-  note = {R package version 0.0.12},
+  note = {R package version 0.1.0},
   url = {https://github.com/pik-piam/reporttransport},
 }
 ```

inst/compareScenarios/cs_01_energy_demand.Rmd

@@ -22,10 +22,10 @@ walk(items, showLinePlots, data = data)
 ```{r}
 tot <- "FE|Transport edge"
 items <- c(
-  "FE|Transport edge|Electricity",
-  "FE|Transport edge|Hydrogen",
-  "FE|Transport edge|Gases",
-  "FE|Transport edge|Liquids"
+  "FE|Transport|Electricity",
+  "FE|Transport|Hydrogen",
+  "FE|Transport|Gases",
+  "FE|Transport|Liquids"
   )
 showAreaAndBarPlots(data, items, tot, orderVars = "user")
 showAreaAndBarPlots(data, items, tot, fill = TRUE, orderVars = "user")
@@ -65,12 +65,12 @@ walk(items, showLinePlots, data = data)
 
 ### Passenger 
 ```{r}
-tot <- "FE|Transport edge|Pass"
+tot <- "FE|Transport|Pass"
 items <- c(
-  "FE|Transport edge|Pass|Electricity",
-  "FE|Transport edge|Pass|Hydrogen",
-  "FE|Transport edge|Pass|Gases",
-  "FE|Transport edge|Pass|Liquids"
+  "FE|Transport|Pass|Electricity",
+  "FE|Transport|Pass|Hydrogen",
+  "FE|Transport|Pass|Gases",
+  "FE|Transport|Pass|Liquids"
   )
 showAreaAndBarPlots(data, items, tot, orderVars = "user")
 showAreaAndBarPlots(data, items, tot, fill = TRUE, orderVars = "user")
@@ -190,12 +190,12 @@ walk(items, showLinePlots, data = data)
 
 ### Freight 
 ```{r}
-tot <- "FE|Transport edge|Freight"
+tot <- "FE|Transport|Freight"
 items <- c(
-  "FE|Transport edge|Freight|Electricity",
-  "FE|Transport edge|Freight|Hydrogen",
-  "FE|Transport edge|Freight|Gases",
-  "FE|Transport edge|Freight|Liquids"
+  "FE|Transport|Freight|Electricity",
+  "FE|Transport|Freight|Hydrogen",
+  "FE|Transport|Freight|Gases",
+  "FE|Transport|Freight|Liquids"
   )
 showAreaAndBarPlots(data, items, tot, orderVars = "user")
 showAreaAndBarPlots(data, items, tot, fill = TRUE, orderVars = "user")
@@ -382,7 +382,7 @@ walk(items, showLinePlots, data = data)
 ### Passenger with bunkers
 ```{r}
 tot_wbunk <- "FE|Transport|Pass with bunkers"
-tot_wobunk <- "FE|Transport edge|Pass"
+tot_wobunk <- "FE|Transport|Pass"
 items <- c(
     "FE|Transport|Bunkers|Pass|International Aviation",
     "FE|Transport|Pass|Domestic Aviation",
@@ -403,7 +403,7 @@ walk(c("FE|Transport|Pass|Aviation", items[1:5], "FE|Transport|Pass|Road|LDV", i
 ### Freight
 ```{r}
 tot_wbunk <- "FE|Transport|Freight with bunkers"
-tot_wobunk <- "FE|Transport edge|Freight"
+tot_wobunk <- "FE|Transport|Freight"
 items <- c(
       "FE|Transport|Bunkers|Freight|International Shipping",
       "FE|Transport|Freight|Domestic Shipping",

inst/compareScenarios/cs_02_energy_services.Rmd

@@ -2,7 +2,7 @@
 ## Passenger with bunkers
 ```{r}
 tot_wbunk <- "ES|Transport|Pass with bunkers"
-tot_wobunk <- "ES|Transport edge|Pass"
+tot_wobunk <- "ES|Transport|Pass"
 items <- c(
     "ES|Transport|Bunkers|Pass|International Aviation",
     "ES|Transport|Pass|Domestic Aviation",
@@ -25,7 +25,7 @@ walk(c( "ES|Transport|Pass|Aviation", tot_wobunk, items[1:5], "ES|Transport|Pass
 ## Passenger per Capita
 ```{r}
 items <- c(
-    "ES|Transport edge|Pass pCap",
+    "ES|Transport|Pass pCap",
     "ES|Transport|Pass|Aviation pCap",
     "ES|Transport|Bunkers|Pass|International Aviation pCap",
     "ES|Transport|Pass|Domestic Aviation pCap",
@@ -123,7 +123,7 @@ walk(items, showLinePlots, data=data)
 ```{r}
 
 items <- c(
-      "ES|Transport edge|Freight pCap",
+      "ES|Transport|Freight pCap",
       "ES|Transport|Bunkers|Freight|International Shipping pCap",
       "ES|Transport|Freight|Domestic Shipping pCap",
       "ES|Transport|Freight|Rail pCap",

inst/compareScenarios/preprocessing.R

@@ -108,7 +108,7 @@ showLinePlotsByVariable <- function(
 pCapVariables <- tribble(
   ~variable, ~newUnit, ~conversionFactor,
   "GDP|PPP", "kUS$2005", 1e6,
-  "ES|Transport edge|Pass", "pkm/yr", 1e9,
+  "ES|Transport|Pass", "pkm/yr", 1e9,
   "ES|Transport|Pass|Aviation", "pkm/yr", 1e9,
   "ES|Transport|Bunkers|Pass|International Aviation", "pkm/yr", 1e9,
   "ES|Transport|Pass|Domestic Aviation", "pkm/yr", 1e9,
@@ -117,7 +117,7 @@ pCapVariables <- tribble(
   "ES|Transport|Pass|Non-motorized|Cycle", "pkm/yr", 1e9,
   "ES|Transport|Pass|Rail|non-HSR", "pkm/yr", 1e9,
   "ES|Transport|Pass|Rail|HSR", "pkm/yr", 1e9,
-  "ES|Transport edge|Freight", "tkm/yr", 1e9,
+  "ES|Transport|Freight", "tkm/yr", 1e9,
   "ES|Transport|Bunkers|Freight|International Shipping", "tkm/yr", 1e9,
   "ES|Transport|Freight|Road", "tkm/yr", 1e9,
   "ES|Transport|Freight|Domestic Shipping", "tkm/yr", 1e9,