Add a reporting script and two templates for the edge_esm module and EDGE-T.

scripts/output/single/EDGETransportValidation.R

+require(rmarkdown)
+
+## run EDGE transport validation output if required
+if(cfg$gms$transport == "edge_esm"){
+  md_template <- "EDGETransportReport.Rmd"
+  file.copy(file.path("./scripts/output/single/notebook_templates", md_template), outputdir)
+  rmarkdown::render(path(outputdir, md_template), output_format="pdf_document")
+  if(cfg$gms$c_keep_iteration_gdxes == 1){
+    md_template <- "EDGETransportMultiIterationAnalysis.Rmd"
+    file.copy(file.path("./scripts/output/single/notebook_templates", md_template), outputdir)
+    rmarkdown::render(path(outputdir, md_template), output_format="pdf_document")
+  }
+}

scripts/output/single/notebook_templates/EDGETransportMultiIterationAnalysis.Rmd

+---
+title: "Analysis of the transport module for a run with resolved iteration domain"
+output: html_document
+---
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+require(ggplot2)
+require(lusweave)
+require(rmndt)
+
+
+knitr::opts_chunk$set(fig.width=12, fig.height=12) 
+
+folder <- "./"
+files <- list.files(path = folder, pattern = "fulldata_[0-9]+\\.gdx")
+
+year_toplot = 2050
+iter_toplot = 25
+maxiter = 100
+
+print(paste0("Year: ", year_toplot))
+print(paste0("Iteration: ", iter_toplot))
+
+```
+
+```{r, echo=FALSE}
+addyrs <- function(dt, yrcol="ttot"){
+    dt[, year := as.numeric(get(yrcol))][, (yrcol) := NULL]
+    return(dt)
+}
+
+get_trp_shares <- function(gdx){
+    tes <- readGDX(gdx, "teEs_dyn35")
+    
+    shares_data <- readgdx(gdx, "pm_shFeCes")[all_teEs %in% tes]
+    shares_data <- addyrs(shares_data)
+    
+    setnames(shares_data, "value", "tech_share")
+    return(shares_data)
+}
+
+get_fuel_prices <- function(gdx){
+  REMINDyears <- c(1990,
+                   seq(2005, 2060, by = 5),
+                   seq(2070, 2110, by = 10),
+                   2130, 2150)
+  ## report prices from REMIND gdx in 2005$/MJ
+  
+  tdptwyr2dpgj <- 31.71  #TerraDollar per TWyear to Dollar per GJ
+  CONV_2005USD_1990USD <- 0.67
+  startyear <- 2020
+  
+  ## load entries from the gdx
+  fety <- readGDX(gdx, c("entyFe", "fety"), format = "first_found")
+  
+  budget.m <- readGDX(gdx, name = "qm_budget", types = "equations", field = "m",
+                      format = "first_found")[, REMINDyears[REMINDyears >= startyear],]  # Alternative: calcPrice
+  
+  interpolate_first_timesteps <- function(obj){
+    ## interpolate values for 1990, 2005 and 2010
+    obj = time_interpolate(obj, c(1990, seq(2005, startyear, 5)),
+                           integrate_interpolated_years=T,
+                           extrapolation_type = "constant")
+    return(obj)
+  }
+  
+  budget.m <- interpolate_first_timesteps(budget.m)
+  budget.m <- lowpass(budget.m)
+  
+  
+  bal_eq <- "qm_balFeForCesAndEs"
+  febal.m <- readGDX(gdx, name = bal_eq, types = "equations",
+                      field = "m", format = "first_found")[, REMINDyears[REMINDyears>=startyear], fety]
+    
+  if(any(febal.m > 0)){
+    sprintf("Found positive marginals on %s. We correct this, but the issue should be adressed.", bal_eq)
+    febal.m[febal.m > 0] <- -1e-10
+  }
+    
+  febal.m <- interpolate_first_timesteps(febal.m)
+    
+  ## in some regions and time steps, 0 final energy demand for an entry could give problems
+  tmp <- magclass::setNames(lowpass(lowpass(febal.m[, , "fegat"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "fegat")
+    
+    
+  tmp <- mbind(tmp, magclass::setNames(lowpass(lowpass(febal.m[, , "feelt"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "feelt"))
+    
+  tmp <- mbind(tmp, magclass::setNames(lowpass(lowpass(febal.m[, , "feh2t"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "feh2t"))
+    
+  tmp <- mbind(tmp, magclass::setNames(lowpass(lowpass(febal.m[, , "fedie"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "fedie"))
+    
+  tmp <- mbind(tmp, magclass::setNames(lowpass(lowpass(febal.m[, , "fepet"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "fepet"))
+                              
+  tmp <- magpie2dt(tmp, regioncol = "all_regi", yearcol = "year", datacols = "all_enty")
+
+  setnames(tmp, old = "value", new = "fe_price")
+  return(tmp)
+}
+
+get_demand <- function(gdx){
+    ces_nodes <- readGDX(gdx, "ppfen_dyn35")
+    demand <- readgdx(gdx, "vm_cesIO")[all_in %in% ces_nodes]
+    demand <- addyrs(demand, "tall")
+    setnames(demand, "value", "CES_demand")
+    return(demand)
+}
+```
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+sorted_files <- paste0(folder, "fulldata_", 1:length(files), ".gdx")
+
+shares <- lapply(sorted_files, get_trp_shares)
+prices <- lapply(sorted_files, get_fuel_prices)
+demand <- lapply(sorted_files, get_demand)
+
+library(stringr)
+# lapply(files, function(fname){str_extract(fname, "[0-9]+")})
+
+for(fname in files){
+    idx <- as.numeric(str_extract(fname, "[0-9]+"))
+    shares[[idx]][, iter := idx]
+    prices[[idx]][, iter := idx]
+    demand[[idx]][, iter := idx]
+}
+
+shares <- rbindlist(shares)
+prices <- rbindlist(prices)
+demand <- rbindlist(demand)
+
+## we are interested only in years up to 2100
+shares = shares[year <= 2100]
+prices = prices[year <= 2100]
+demand = demand[year <= 2100]
+```
+
+##Prices and shares in the iteration domain
+
+```{r, echo=FALSE}
+toplot <- prices[year == year_toplot  & iter < maxiter]
+
+ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+    geom_line() +
+    labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Iteration Domain, in ", year_toplot)) +
+    facet_wrap(~ all_regi, ncol=2, scales = "free" )
+# +
+#   ylim(c(0,50))
+```
+
+```{r, echo=FALSE}
+# toplot <- prices[year %in% c(year_toplot-5,year_toplot,year_toplot+5)]
+# toplot <- toplot[,.(fe_price=mean(fe_price)), by=c("all_regi","all_enty","iter")]
+# 
+# ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line() +
+#     labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Iteration Domain, average of ", year_toplot-5, ", ", year_toplot, ", ", year_toplot+5)) +
+#     facet_wrap(~ all_regi, ncol=2, scales = "free")
+
+
+```
+
+
+```{r, echo=FALSE}
+# toplot <- prices[year %in% c(year_toplot-5,year_toplot)]
+# toplot <- toplot[,.(fe_price=mean(fe_price)), by=c("all_regi","all_enty","iter")]
+# 
+# ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line() +
+#     labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Iteration Domain, average of ", year_toplot-5, ", ", year_toplot)) +
+#     facet_wrap(~ all_regi, ncol=2, scales = "free")
+
+
+```
+
+
+
+```{r, echo=FALSE}
+# toplot <- prices[year %in% c(year_toplot-10,year_toplot-5,year_toplot)]
+# toplot <- toplot[,.(fe_price=mean(fe_price)), by=c("all_regi","all_enty","iter")]
+# 
+# ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line() +
+#     labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Iteration Domain, average of ", year_toplot-10, ", ", year_toplot-5, ", ", year_toplot)) +
+#     facet_wrap(~ all_regi, ncol=2, scales = "free")
+
+
+```
+
+
+```{r, echo=FALSE}
+# toplot <- prices[year == year_toplot & all_regi == "EUR"]
+# 
+# ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line(size = 2) +
+#     labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("European Fuel Prices in the Iteration Domain, in ", year_toplot)) +
+#    # facet_wrap(~ all_regi, ncol=2)+
+#   theme_light()+
+#   scale_color_brewer(palette = "Set1") 
+```
+
+
+```{r, echo =FALSE}
+# toplot <- shares[year == year_toplot & iter < maxiter]
+# 
+# ggplot(aes(iter, tech_share, colour=all_teEs, shares=all_teEs), data=toplot) +
+#     geom_line() +
+#     labs(x = "iteration #", y = "Tech Share", color = "LDV Tech Type", title=paste0("Tech Shares in the Iteration Domain, in ", year_toplot)) +
+#     facet_wrap(~ all_regi, ncol=2)
+```
+
+```{r, echo=FALSE}
+
+toplot <- demand[year == year_toplot & all_in %in% c("entrp_pass_sm", "entrp_pass_lo")  & iter < maxiter]
+
+ggplot(aes(iter, CES_demand, colour=all_in, shares=all_in), data=toplot) +
+    geom_line() +
+    labs(x = "iteration #", y = "CES node aggr. pass. demand (trillion pkm)", color = "CES category", title=paste0("Aggregated Passenger Tranport Demand in the Iteration Domain, in ", year_toplot)) +
+    facet_wrap(~ all_regi, ncol=2)
+```
+
+
+```{r, echo=FALSE}
+
+# toplot <- demand[year == year_toplot & all_in %in% c("entrp_pass_sm", "entrp_pass_lo") & all_regi == "EUR"]
+# 
+# ggplot(aes(iter, CES_demand, colour=all_in, shares=all_in), data=toplot) +
+#     geom_line(size = 2) +
+#     labs(x = "iteration #", y = "CES node aggr. pass. demand (trillion pkm)", color = "CES category", title=paste0("Europe: aggregated Passenger Tranport Demand in the Iteration Domain, in ", year_toplot)) +
+#    # facet_wrap(~ all_regi, ncol=2)+
+#   theme_light()
+```
+
+
+```{r, echo=FALSE}
+toplot <- demand[year == year_toplot & all_in %in% c("entrp_frgt_sm", "entrp_frgt_lo")  & iter < maxiter]
+ggplot(aes(iter, CES_demand, colour=all_in, shares=all_in), data=toplot) +
+    geom_line() +
+    labs(x = "iteration #", y = "CES node aggr. freight demand (trillion tkm)", color = "CES category", title= paste0("Aggregated Freight Transport Demand in the Iteration Domain, in ", year_toplot)) +
+    facet_wrap(~ all_regi, ncol=2)
+```
+
+##Prices in the time domain
+
+```{r, echo=FALSE}
+# toplot <- prices[iter == iter_toplot]
+# 
+# ggplot(aes(year, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line() +
+#     labs(x = "Year", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Time Domain, in iteration ", iter_toplot)) +
+#     facet_wrap(~ all_regi, ncol=2)
+# +
+#   ylim(c(0,50))
+```
+
+```{r, echo=FALSE}
+
+enties = c("fegat", "fepet", "fedie", "feelt", "feh2t")
+
+for (enty in enties) {
+  toplot <- prices[all_enty == enty  & iter < maxiter]
+  toplot[, iter := as.numeric(iter)]
+
+  plot=ggplot(aes(year, fe_price, group=interaction(iter, all_enty), color=iter), data=toplot) +
+    geom_line() +
+    labs(x = "Year", y = "US$ 2015/GJ", color = "Iteration number", title=paste0("Fuel Prices in the Time Domain, across iterations, for ", enty)) +
+    facet_wrap(~ all_regi, scales = "free", ncol = 2)
+
+  print(plot)
+}
+
+```
+
+```{r, echo =FALSE}
+
+
+enty_vals = unique(prices$all_enty)
+toplot <- prices[year>=2005 & iter < maxiter]
+toplot[, iter := as.numeric(iter)]
+toplot[, min := min(fe_price), by = c("all_regi", "all_enty", "year")]
+toplot[, max := max(fe_price), by = c("all_regi", "all_enty", "year")]
+
+plot= ggplot()+
+  geom_ribbon(data = toplot, aes(x=year, group = all_enty, fill = all_enty, ymin = min, ymax = max), alpha = 0.5)+
+  facet_wrap(~ all_regi, scales = "free", ncol = 2)+
+  theme_minimal() +
+  facet_wrap(~ all_regi,ncol = 2) +
+  geom_line(data = toplot[iter == max(iter)], aes(year, fe_price, group=all_enty, color = all_enty), linetype = "dashed")+
+    labs(x = "Year", y = "US$ 2015/GJ", fill = "Fuel Type", color = "Final iteration", title=paste0("Fuel Prices in the Time Domain, across iterations"))
+
+
+print(plot)
+```
+
+
+## Shares in the time domain
+
+```{r, echo =FALSE}
+
+
+in_vals = unique(shares$all_in)
+
+for (in_val in in_vals) {
+  toplot <- shares[year>=2005 & all_in == in_val  & iter < maxiter]
+  toplot[, iter := as.numeric(iter)]
+
+  plot=ggplot() +
+    geom_line(aes(year, tech_share, group=interaction(iter, all_teEs), color=all_teEs, alpha = iter, linetype ="Intermediate Iter." ), data=toplot) +
+    geom_line(data = toplot[iter == max(iter)], aes(year, tech_share, group=all_teEs, linetype ="Last Iter.", color ="Last Iter."), color="black")+
+    scale_linetype_manual("Iteration Status",values=c("Intermediate Iter." = 1, "Last Iter."=2))+
+    labs(x = "Year", y = "[-]", color = "Fuel Type", title=paste0("Tech shares in the Time Domain, across iterations, for ", in_val)) +
+    facet_wrap(~ all_regi, scales = "free", ncol = 2)+
+    theme_minimal() +
+    facet_wrap(~ all_regi,ncol = 2)
+
+  print(plot)
+  
+}
+```
+
+
+
+
+```{r, echo =FALSE}
+
+
+in_vals = unique(shares$all_in)
+
+for (in_val in in_vals) {
+  toplot <- shares[year>=2005 & all_in == in_val  & iter < maxiter]
+  toplot[, iter := as.numeric(iter)]
+  toplot[, min := min(tech_share), by = c("all_regi", "all_teEs", "year")]
+  toplot[, max := max(tech_share), by = c("all_regi", "all_teEs", "year")]
+
+ plot= ggplot()+
+  geom_ribbon(data = toplot, aes(x=year, group = all_teEs, fill = all_teEs, ymin = min, ymax = max), alpha = 0.5)+
+  facet_wrap(~ all_regi, scales = "free", ncol = 2)+
+  theme_minimal() +
+  facet_wrap(~ all_regi,ncol = 2) +
+  geom_line(data = toplot[iter == max(iter)], aes(year, tech_share, group=all_teEs, color = all_teEs), linetype = "dashed")+
+    labs(x = "Year", y = "[-]", color = "Fuel Type", title=paste0("Tech shares in the Time Domain, across iterations, for ", in_val))
+
+
+  print(plot)
+  
+}
+```
+
+
+
+```{r, echo =FALSE}
+
+
+in_vals = unique(shares$all_in)
+
+for (in_val in in_vals) {
+  toplot <- demand[year>=2005 & all_in == in_val  & iter < maxiter]
+  toplot[, iter := as.numeric(iter)]
+  if (in_val %in% c("entrp_pass_sm","entrp_pass_lo")) {
+    ylabel = "[billion pkm]"
+  } else {
+    ylabel = "[billion tkm]"
+  }
+
+ plot= ggplot()+
+   geom_line(data = toplot, aes(x=year, y = CES_demand, group = iter, color = all_in, alpha = iter, linetype ="Intermediate Iter." ))+
+   geom_line(data = toplot[iter == max(iter)], aes(year, CES_demand, group=all_in, linetype ="Last Iter.", color ="Last Iter."), color="black")+
+   facet_wrap(~ all_regi, scales = "free", ncol = 2)+
+   theme_minimal() +
+   facet_wrap(~ all_regi,ncol = 2) +
+   scale_linetype_manual("Iteration Status",values=c("Intermediate Iter." = 1, "Last Iter."=2))+
+    labs(x = "Year", y = ylabel, color = "Transport Type", title=paste0("CES node value in Time Domain, across iterations, for ", in_val))
+
+
+  print(plot)
+  
+}
+```
+## A quality measure of convergence
+
+For starters, let us implement something rather simple:
+
+$$Q{r,s}=\sum_{t}\frac{\sum_{i=m}^N [s_{t,i}-Avg_i(s_{t,i})]^2}{Avg^2_i(s_{t,i})}$$
+
+
+for a variable $s$ over times $t$ with iteration index $i$ (first m iterations being ignored) and where $Avg_i$ denotes the average over the iteration index dimension. The region index $r$ is omitted on the right hand side.
+
+
+```{r, echo =FALSE}
+quality <- function(arr, m=3){
+    # we ignore the first m iterations
+    sum((arr[-(1:m)] - ave(arr[-(1:m)]))^2/ave(arr[-(1:m)])^2)
+}
+
+shares[, quality := quality(.SD$tech_share), by=c("all_regi", "all_teEs", "year")]
+top <- head(shares[order(-quality), max(quality), by=c("all_regi", "all_teEs", "year")], 15)
+top
+
+```
+
+```{r, echo=FALSE}
+ggplot(aes(iter, tech_share), data=shares[all_regi == top[1]$all_regi & all_teEs == top[1]$all_teEs & year == top[1]$year]) +
+    geom_line() +
+    labs(x = "iteration #", y = "Tech Share", color = "year", title=paste(top[1]$all_te, "Shares in", top[1]$all_regi, "in the Iteration Domain, in ", top[1]$year))
+```
+

scripts/output/single/notebook_templates/EDGETransportReport.Rmd

+---
+title: "Compare scenarios Transport"
+output:
+  html_document:
+    df_print: paged
+---
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+require(ggplot2)
+require(moinput)
+require(data.table)
+require(dplyr)
+require(remind)
+require(gdxdt)
+require(gdx)
+require(rmndt)
+require(data.table)
+
+require(edgeTrpLib)
+```
+
+
+```{r, echo=FALSE, warning=FALSE}
+
+dem_shares <- list()
+intensity <- list()
+demand_km <- list()
+demand_ej <- list()
+sw_tech <- list()
+prices_FV <- list()
+
+datapath <- function(fname){
+  file.path("./input_EDGE/", fname)
+}
+
+mapspath <- function(fname, scenariopath=""){
+    file.path("../../modules/35_transport/edge_esm/input", fname)
+}
+
+## Load mappings
+EDGE2CESmap <- fread(mapspath("mapping_CESnodes_EDGE.csv"))
+
+REMIND2ISO_MAPPING <- fread("../../config/regionmappingH12.csv")[, .(iso = CountryCode,
+                                                                         region = RegionCode)]
+
+EDGE2teESmap <- fread(mapspath("mapping_EDGE_REMIND_transport_categories.csv"))
+
+years <- c(1990,
+           seq(2005, 2060, by = 5),
+           seq(2070, 2110, by = 10),
+           2130, 2150)
+
+REMINDyears <- c(1990,
+           seq(2005, 2060, by = 5),
+           seq(2070, 2110, by = 10),
+           2130, 2150)
+
+## include the paths to the scenarios you want to compare
+output_folders <- "./"
+
+for(output_folder in output_folders){
+  ## load gdx for fuel prices and demand
+  gdx = file.path(output_folder, "fulldata.gdx")
+  ## load policy scenario
+  load(file.path(output_folder, "config.Rdata"))
+  REMIND_scenario <- cfg$gms$cm_GDPscen
+  EDGE_scenario <- cfg$gms$cm_EDGEtr_scen
+  policy_scenario <- cfg$gms$c_expname
+  
+  scen <- paste0(REMIND_scenario, "-", EDGE_scenario, "-", policy_scenario)
+  
+  ## load demand
+  ES_demand = readREMINDdemand(gdx, REMIND2ISO_MAPPING, EDGE2teESmap, years)
+
+  ## load input data
+  int_dat <- readRDS(datapath("harmonized_intensities.RDS"))
+  nonfuel_costs <- readRDS(datapath("UCD_NEC_iso.RDS"))
+  sw_data <- readRDS(datapath("SW.RDS"))
+  vot_data <- readRDS(datapath("VOT_iso.RDS"))
+  logit_params <- readRDS(datapath("logit_exp.RDS"))
+  price_nonmot <- readRDS(datapath("price_nonmot.RDS"))
+  
+    ## FIXME: hotfix to make the (empty) vot_data$value_time_VS1 with the right column types. Probably there is another way to do that, did not look for it.
+  vot_data$value_time_VS1$iso = as.character(vot_data$value_time_VS1$iso)
+  vot_data$value_time_VS1$subsector_L1 = as.character(vot_data$value_time_VS1$subsector_L1)
+  vot_data$value_time_VS1$vehicle_type = as.character(vot_data$value_time_VS1$vehicle_type)
+  vot_data$value_time_VS1$year = as.numeric(vot_data$value_time_VS1$year)
+  vot_data$value_time_VS1$time_price = as.numeric(vot_data$value_time_VS1$time_price)
+  
+  ## calculate prices
+  REMIND_prices <- merge_prices(
+    gdx = gdx,
+    REMINDmapping = REMIND2ISO_MAPPING,
+    REMINDyears = REMINDyears,
+    intensity_data = int_dat,
+    nonfuel_costs = nonfuel_costs)
+
+  ## calculates logit
+  logit_data <- calculate_logit(
+    REMIND_prices[tot_price > 0],
+    REMIND2ISO_MAPPING,
+    vot_data = vot_data,
+    sw_data = sw_data,
+    logit_params = logit_params,
+    intensity_data = int_dat,
+    price_nonmot = price_nonmot)
+
+
+  shares <- logit_data[["share_list"]] ## shares of alternatives for each level of the logit function
+  mj_km_data <- logit_data[["mj_km_data"]] ## energy intensity at a technology level
+  prices_FV[[scen]] <- REMIND_prices[, EDGE_scenario := scen] ## prices at each level of the logit function, 1990USD/pkm
+
+  ## calculate energy intensity and FE demand at a REMIND-region level for the desired level of aggregation
+  res <- shares_intensity_and_demand(
+    logit_shares=shares,
+    MJ_km_base=mj_km_data,
+    REMIND2ISO_MAPPING=REMIND2ISO_MAPPING,
+    EDGE2CESmap=EDGE2CESmap,
+    REMINDyears=REMINDyears,
+    demand_input=ES_demand)
+
+  dem_shares[[scen]] <- res$demand[, EDGE_scenario := scen]
+  intensity[[scen]] <- res$demandI[, EDGE_scenario := scen]
+  demand_km[[scen]] <- res$demandF_plot_pkm[, EDGE_scenario := scen]
+  demand_ej[[scen]] <- res$demandF_plot_EJ[, EDGE_scenario := scen]
+  sw_tech[[scen]] <- sw_data$FV_final_SW[, EDGE_scenario := scen]
+}
+
+dem_shares <- rbindlist(dem_shares)
+intensity <- rbindlist(intensity)
+demand_km <- rbindlist(demand_km)
+demand_ej <- rbindlist(demand_ej)
+sw_tech <- rbindlist(sw_tech)
+prices_FV <- rbindlist(prices_FV)
+```
+
+
+
+```{r, echo=FALSE}
+## plot settings
+years_plot = c(2010,2015,2020,2025,2030,2040,2050) ## in bar charts, these are the time steps that are represented
+year_single = 2050
+region_plot = "NEU" ## in case is a region specific plot, this region is represented
+sector_plot ="trn_pass" ## in case is a sector specific plot, this sector is represented
+
+
+##conversion rate 2005->1990 USD
+CONV_2005USD_1990USD=0.67
+
+# print(paste0("Scenario: ", REMIND_scenario))
+print(paste0("Regional plots are about ",region_plot))
+print(paste0("Sectoral plots are about ",sector_plot))
+
+## maps
+cesmap <- data.table(CES_parent=c("_p_sm", "_p_lo", "_f_lo", "_f_sm"),
+                     CES_label=c("Passenger, Short-to-Medium Distances",
+                                 "Passenger, Long Distances",
+                                 "Freight, Long Distances",
+                                 "Freight, Short-to-Medium Distances"))
+
+EDGE_sectormap <- data.table(sector=c("trn_pass", "trn_freight", "trn_aviation_intl", "trn_shipping_intl"),
+                     CES_label=c("Passenger, Short-to-Medium Distances",
+                                 "Passenger, Long Distances",
+                                 "Freight, Long Distances",
+                                 "Freight, Short-to-Medium Distances"))
+
+```
+
+
+```{r, echo=FALSE}
+## aggregate demands to REMIND regions
+demandF_plot_EJ <- demand_ej[,c("EDGE_scenario", "sector","subsector_L3","subsector_L2",
+                                "subsector_L1","vehicle_type","technology", "iso","year","demand_EJ")]
+demandF_plot_pkm <- demand_km[,c("EDGE_scenario", "sector","subsector_L3","subsector_L2",
+                                 "subsector_L1","vehicle_type","technology","iso","year","demand_F")]
+
+demandF_plot_EJ=aggregate_dt(demandF_plot_EJ,REMIND2ISO_MAPPING,
+                             datacols = c("EDGE_scenario", "sector", "subsector_L3", "subsector_L2", "subsector_L1",
+                                          "vehicle_type", "technology"),
+                             valuecol = "demand_EJ")
+
+demandF_plot_pkm=aggregate_dt(demandF_plot_pkm,REMIND2ISO_MAPPING,
+                              datacol = c("EDGE_scenario", "sector","subsector_L3","subsector_L2",
+                                          "subsector_L1","vehicle_type","technology"),
+                              valuecol = "demand_F")
+## add GLO
+glo <- demandF_plot_pkm[,.(region="GLO", demand_F=sum(demand_F)),
+                        by=eval(names(demandF_plot_pkm)[2:9])]
+demandF_plot_pkm <- rbind(demandF_plot_pkm, glo)
+glo <- demandF_plot_EJ[,.(region="GLO", demand_EJ=sum(demand_EJ)),
+                        by=eval(names(demandF_plot_EJ)[2:9])]
+demandF_plot_EJ <- rbind(demandF_plot_EJ, glo)
+
+```
+
+## ES
+
+
+```{r, echo=FALSE}
+##chunk of code that plots the ES
+ES_modes_bar1=function(demandpkm){
+  #group by subsector_L3 and summarise the demand
+  df=demandpkm[, .(demand_F=sum(demand_F)
+                   ), by = c("EDGE_scenario", "region", "year","sector","subsector_L1")]
+  df[,demand_F:=demand_F   ## in millionkm
+     *1e-6                 ## in trillion km
+     ]
+  df=df[order(year)]
+  # #filter only 2020, 2050 and 2100
+  df=df[year %in% years_plot,]
+  #separate into passenger and freight categories
+  pass=c("trn_pass","trn_aviation_intl")
+  freight=c("trn_freight","trn_shipping_intl")
+  ## give proper names to the categories  
+  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
+  #plot
+  plot_p=ggplot()+
+    geom_bar(data=df%>%filter(sector %in% pass, year %in% years_plot, region == region_plot),
+             aes(x=year,y=demand_F,group=mode,fill=mode),position=position_stack(),stat="identity")+
+    facet_wrap(~EDGE_scenario)+
+    ggtitle("Energy Services Demand - Passenger Transport Modes")+
+    theme_light()+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1))+
+    scale_x_continuous(breaks=years_plot)+
+    xlab("Year")+
+    ylab("Energy Services Demand (trillion pkm)")+ 
+    guides(fill=guide_legend(title="Transport mode"))+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1),
+          axis.text = element_text(size=13),
+          title = element_text(size=13),
+          legend.text = element_text(size=13))
+  
+  plot_f=ggplot()+
+    geom_bar(data=df%>%filter(sector %in% freight,year %in% years_plot, region == region_plot),
+             aes(x=year,y=demand_F,group=mode,fill=mode),position=position_stack(),stat="identity")+
+    facet_wrap(~EDGE_scenario)+
+    ggtitle("Energy Services Demand - Freight Transport Modes")+
+    theme_light()+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1))+
+    scale_x_continuous(breaks=years_plot)+
+    xlab("Year")+
+    ylab("Energy Services Demand (trillion tkm)")+ 
+    guides(fill=guide_legend(title="Transport mode"))+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1),
+          axis.text = element_text(size=13),
+          title = element_text(size=13),
+          legend.text = element_text(size=13))
+  
+  plot=list(plot_p,plot_f)
+  return(plot)
+}
+
+p=ES_modes_bar1(demandpkm=demandF_plot_pkm)
+p[[1]]
+p[[2]]
+
+```
+
+```{r, echo=FALSE}
+##chunk of code that plots the ES
+ES_modes_bar=function(demandpkm){
+  demandpkm[technology == "LA-BEV", technology := "BEV"]
+  ## use proper non-fuel mode names
+  demandpkm[technology %in% c("Cycle_tmp_technology", "Walk_tmp_technology"), technology := "Human Powered"]
+  #group by subsector_L3 and summarise the demand
+  df=demandpkm[, .(demand_F=sum(demand_F)),
+               by = c("EDGE_scenario", "region", "year","sector","subsector_L1", "technology")]
+  df[,demand_F:=demand_F   ## in millionkm
+     *1e-6                 ## in trillion km
+     ]
+  df=df[order(year)]
+  #separate into passenger and freight categories
+  pass=c("trn_pass","trn_aviation_intl")
+  freight=c("trn_freight","trn_shipping_intl")
+  ## give proper names to the categories
+  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
+  #plot
+  plot_p=ggplot()+
+    geom_bar(data=df%>%filter(sector %in% pass,year %in% years_plot,region==region_plot),
+             aes(x=year,y=demand_F,group=technology,fill=technology),
+             position=position_stack(),stat="identity")+
+    ggtitle("Energy Services Demand by Technology, Passenger Transport, EUR")+
+    theme_light()+
+    facet_wrap(~EDGE_scenario)+
+    scale_x_continuous(breaks=years_plot)+
+    xlab("Year")+
+    ylab("Energy Services Demand (trillion pkm)")+ 
+    guides(fill=guide_legend(title="Technology"))+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1),
+          axis.text = element_text(size=13),
+          title = element_text(size=13),
+          legend.text = element_text(size=13))
+  
+  plot_f=ggplot()+
+    geom_bar(data=df%>%filter(sector %in% freight,year %in% years_plot, region==region_plot),
+             aes(x=year,y=demand_F, group=technology, fill=technology),
+             position=position_stack(),stat="identity")+
+    ggtitle(paste0("Energy Services Demand by Technology, Freight Transport, ", region_plot))+
+    facet_wrap(~EDGE_scenario) +
+    theme_light()+
+    scale_x_continuous(breaks=years_plot)+
+    xlab("Year")+
+    ylab("Energy Services Demand (trllion tkm)")+ 
+    guides(fill=guide_legend(title="Technology"))+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1),
+          axis.text = element_text(size=13),
+          title = element_text(size=13),
+          legend.text = element_text(size=13))
+  
+  plot=list(plot_p,plot_f)
+  return(plot)
+}
+
+p=ES_modes_bar(demandpkm=demandF_plot_pkm)
+p[[1]]
+p[[2]]
+
+```
+
+```{r, echo=FALSE}
+## plot ES for LDVs only divided by fuel
+ES_modes_LDV_bar=function(demandpkm){
+  demandpkm[technology == "LA-BEV", technology := "BEV"]
+  ## use proper non-fuel mode names
+  demandpkm[technology %in% c("Cycle_tmp_technology", "Walk_tmp_technology"), technology := "Human Powered"]
+  #group by subsector_L3 and summarise the demand
+  df=demandpkm[, .(demand_F=sum(demand_F)),
+               by = c("EDGE_scenario", "region", "year","sector","subsector_L1", "technology")]
+  df[,demand_F:=demand_F   ## in millionkm
+     *1e-6                 ## in trillion km
+     ]
+  df=df[order(year)]
+  ## give proper names to the categories
+  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
+  
+  
+  ## select order of facets
+  df$technology = factor(df$technology, levels=c("Liquids","Hybrid Liquids","NG","BEV","FCEV"))
+  
+  #plot
+  plot_LDV=ggplot()+
+    geom_bar(data=df%>%filter(year %in% years_plot,region==region_plot, mode %in% c("4W","2W")),
+             aes(x=year,y=demand_F,group=technology,fill=technology), alpha = 0.9,
+             position=position_stack(),stat="identity")+
+    ggtitle(paste0("Energy Services Demand by Technology, LDVs", region_plot))+
+    theme_light()+
+    facet_wrap(~EDGE_scenario)+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1))+
+    xlab("Year")+
+    ylab("Energy Services Demand (trillion pkm)")+ 
+    guides(fill=guide_legend(title="Technology"))+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1),
+          axis.text = element_text(size=13),
+          title = element_text(size=13),
+          legend.text = element_text(size=13),
+          strip.text.x = element_text(size = 13, color = "black"),
+          strip.background=element_rect(fill="white"))+
+    scale_x_continuous(breaks=years_plot)+
+    scale_fill_brewer(palette = "Set1")
+    
+  
+  return(plot_LDV)
+}
+
+p=ES_modes_LDV_bar(demandpkm=demandF_plot_pkm)
+
+p
+
+```
+
+## FE
+
+```{r, echo=FALSE}
+FE_modes_bar=function(demandEJ){
+  #group by subsector_L1 and summarise the demand
+  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
+              by = c("EDGE_scenario", "region", "year","subsector_L1","subsector_L3","sector")]
+  df=df[order(year)]
+  df=df[year %in% years_plot,]
+  #separate into passenger and freight categories
+  pass=c("trn_pass","trn_aviation_intl")
+  freight=c("trn_freight","trn_shipping_intl")
+  ## give proper names to the categories
+  df <- merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
+  #plot
+  plot_p=ggplot()+
+    geom_bar(data=df%>%filter(sector %in% pass, region==region_plot),
+             aes(x=year,y=demand_EJ,group=mode,fill=mode),position=position_stack(),stat="identity",color="black")+
+    facet_wrap(~EDGE_scenario)+
+    ylab("Energy (EJ)") +
+    ggtitle(paste0("Final Energy Demand by Mode, Passenger, ", region_plot))+
+    theme(axis.text.x = element_text(angle = 90))+
+    scale_x_continuous(breaks=years_plot)
+  
+  plot_f=ggplot()+
+    geom_bar(data=df%>%filter(sector %in% freight,region==region_plot),
+             aes(x=year,y=demand_EJ,group=mode,fill=mode),position=position_stack(),stat="identity",color="black")+
+    facet_wrap(~EDGE_scenario) +
+    ggtitle(paste0("Final Energy Demand by Mode, Freight, ", region_plot))+
+    ylab("Energy (EJ)") +
+    theme(axis.text.x = element_text(angle = 90))+
+    scale_x_continuous(breaks=years_plot)
+  
+  plot=list(plot_p,plot_f)
+  return(plot)
+}
+
+p=FE_modes_bar(demandEJ = demandF_plot_EJ)
+p[[1]]
+p[[2]]
+```
+
+```{r, echo=FALSE}
+## function that calculates FE split and splits out the liquids by source
+FE_modes_bar_oilcomponent=function(demandEJ, msect="trn_pass", region_plot){
+  if(msect == "trn_pass")
+    sector_display = "Passenger"
+  if(msect == "trn_freight")
+    sector_display = "Freight"
+    
+  #group by subsector_L1 and summarise the demand
+  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
+              by = c("EDGE_scenario", "region", "year","subsector_L1","subsector_L3","sector","technology")]
+  df=df[order(year)]
+  df=df[year %in% years_plot,]
+  ## give proper names to the categories
+  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
+  
+  df[,technology := ifelse(technology == "LA-BEV", "BEV", technology)]
+  df[,technology := ifelse(technology == "Electric", "El. Trains", technology)]
+  ## select order of facerts
+  df$technology = factor(df$technology, levels=c("BEV","FCEV","Hybrid Liquids", "El. Trains", "NG","Liquids", "Coal"))
+
+  #plot
+  plot_psm = ggplot()+
+    geom_bar(data=df%>%filter(sector == msect, region == region_plot),
+             aes(x=year,y=demand_EJ,group=technology,fill=technology),
+             position=position_stack(),stat="identity", alpha = 0.9)+
+    theme_light()+
+    ggtitle(paste0("Final Energy Demand by Tech, ", sector_display, ", ", region_plot))+
+    theme(axis.text.x = element_text(angle = 90),
+        strip.text.x = element_text(size = 13, color = "black"))+
+    scale_x_continuous(breaks=years_plot)+
+    scale_fill_brewer(palette = "Set1")+
+    xlab("Year")+
+    ylab("Final energy demand [EJ]")+ 
+    facet_wrap(~EDGE_scenario) +
+    guides(fill=guide_legend(title="Technology"))
+  
+    plot_psm_LDV = ggplot()+
+    geom_bar(data=df%>%filter(sector == msect, region == region_plot, mode == "4W"),
+             aes(x=year,y=demand_EJ,group=technology,fill=technology),
+             position=position_stack(),stat="identity", alpha = 0.9)+
+    theme_light()+
+    ggtitle(paste0("Final Energy Demand by Tech, LDVs, ", region_plot))+
+    theme(axis.text.x = element_text(angle = 90),
+          strip.text.x = element_text(size = 13, color = "black"),
+          strip.background=element_rect(fill="white"),
+          axis.text = element_text(size=13),
+          title = element_text(size=13),
+          legend.text = element_text(size=13))+
+    scale_x_continuous(breaks=years_plot)+
+    scale_fill_brewer(palette = "Set1")+
+    xlab("Year")+
+    ylab("Final energy demand [EJ]")+ 
+    facet_wrap(~EDGE_scenario) +
+    guides(fill=guide_legend(title="Technology"))  
+  
+    plot_list = list(plot_psm, plot_psm_LDV)
+  return(plot_list)
+}
+
+
+
+FE_modes_bar_oilcomponent(demandEJ = demandF_plot_EJ, msect="trn_pass", region=region_plot)
+FE_modes_bar_oilcomponent(demandEJ = demandF_plot_EJ, msect="trn_pass", region="GLO")
+
+```
+
+
+```{r, echo=FALSE}
+FE_modes_bar1=function(demandEJ, msect="trn_pass", region_plot){
+  if(msect == "trn_pass")
+    sector_display = "Passenger"
+  if(msect == "trn_freight")
+    sector_display = "Freight"
+
+  ## group by subsector_L1 and summarise the demand
+  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
+              by = c("EDGE_scenario", "region", "year","subsector_L1","sector")]
+  df=df[order(year) & year %in% years_plot]
+  ## give proper names to the categories
+  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
+  #plot
+  plot_p=ggplot()+
+    geom_bar(data=df%>%filter(sector == msect, region == region_plot),
+             aes(x=year,y=demand_EJ,group=mode,fill=mode),
+             position=position_stack(),stat="identity",color="black")+
+    facet_wrap(~EDGE_scenario) +
+    ylab("Energy (EJ)") +
+    ggtitle(paste0("Final Energy Demand by Mode, ", sector_display, ", ", region_plot))+
+    theme(axis.text.x = element_text(angle = 90))+
+    scale_x_continuous(breaks=years_plot)
+  
+  return(plot_p)
+}
+
+FE_modes_bar1(demandEJ = demandF_plot_EJ, region_plot = region_plot)
+```
+
+
+
+```{r, echo=FALSE}
+FE_modes_bar_oilVSelec=function(demandEJ, region_plot){
+  sector_display = "Total transport"
+  ## group by subsector_L1 and summarise the demand
+  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
+              by = c("EDGE_scenario", "region", "year","subsector_L1","sector","technology")]
+  df[, tech_plot := ifelse(technology %in% c("BEV","Electric"), "Electriciy", NA)]
+  df[, tech_plot := ifelse(technology %in% c("Liquids", "Hybrid Liquids"), "Liquids", tech_plot)]
+  df=df[!is.na(tech_plot),] ## only liquids and electric driven entries interesting
+  df=df[order(year) & year %in% years_plot]
+  ## give proper names to the categories
+  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
+  df[,short_names:=ifelse(mode %in% c("Buses","Rail Passenger","High Speed Rail"),"Other Passenger",NA)]
+  df[,short_names:=ifelse(mode %in% c("2W","4W","Three Wheelers"),"LDV",short_names)]
+  df[,short_names:=ifelse(mode %in% c("Domestic Aviation","International Aviation"),"Aviation",short_names)]
+  df[,short_names:=ifelse(mode %in% c("International Shipping","Domestic Shipping"),"Shipping",short_names)]
+  df[,short_names:=ifelse(mode %in% c("Road Freight","Rail Freight"),"Road and Rail Freight",short_names)]
+  
+  #plot
+  plot_p=ggplot()+
+    geom_bar(data=df%>%filter(region == region_plot, year ==2050),
+             aes(x=tech_plot,y=demand_EJ,group=short_names,fill=short_names),
+             position=position_stack(),stat="identity",alpha=0.95)+
+    facet_wrap(~EDGE_scenario) +
+    ylab("Energy (EJ)") +
+    ggtitle(paste0("Final Energy Demand by Mode in 2050, total transport ", region_plot))+
+    theme_light()+
+    theme(axis.text.x = element_text(angle = 90),
+          axis.title.x = element_blank(),
+          strip.text.x = element_text(size = 13, color = "black"),
+          strip.background=element_rect(fill="white"),
+          axis.text = element_text(size=13),
+          title = element_text(size=13),
+          legend.text = element_text(size=13))+
+    scale_fill_brewer(palette = "Set2")+
+    guides(fill=guide_legend(title="Transport mode"))
+  
+  return(plot_p)
+}
+
+FE_modes_bar_oilVSelec(demandEJ = demandF_plot_EJ, region_plot = region_plot)
+```
+
+
+## FE composition
+
+```{r, echo=FALSE}
+
+FE_modeshares_area=function(demandEJ){
+  #group by subsector_L3 and summarise the demand
+  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
+              by = c("EDGE_scenario", "region", "year","subsector_L3")]
+  #order by year
+  df=df[order(year)]
+  df=df[year>=2005,]
+  #plot
+  plot=ggplot()+
+    geom_area(data=df %>% filter(year <= max(years_plot), region == region_plot),
+              aes(x=year,y=demand_EJ,group=subsector_L3,fill=subsector_L3),
+              color="black")+
+    facet_wrap(~EDGE_scenario)+
+    ylab("Energy (EJ)") +
+    ggtitle("Final Energy Demand, Mode composition")+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1))
+  
+  return(plot)
+}
+
+p=FE_modeshares_area(demandEJ = demandF_plot_EJ)
+p
+## ggsave("FE_modeshares.png")
+```
+
+
+
+```{r, echo=FALSE}
+## fuel use by sector
+fuel_shares_area=function(demandEJ, msect="trn_pass", region_plot){
+  if(msect == "trn_pass")
+    sector_display = "Passenger"
+  if(msect == "trn_freight")
+    sector_display = "Freight"
+  ##group by sector and technology and summarise demand
+  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
+              by = c("EDGE_scenario", "region", "year","technology","sector")]
+
+  df=df[order(year) & year>=2005,]
+  #plot
+  plot1=ggplot()+
+    geom_area(data=df%>%filter(sector == msect, year <= max(years_plot), region == region_plot),
+              aes(x=year,y=demand_EJ,group=technology,fill=technology),position="fill")+
+    facet_wrap(~EDGE_scenario)+
+    ggtitle(paste0("Final Energy Demand, ", sector_display, ", Fuel Composition"))+
+    ylab("Share")
+    theme(axis.text.x = element_text(angle = 90, hjust = 1))
+    
+  return(plot1)
+}
+
+p=fuel_shares_area(demandEJ = demandF_plot_EJ, region_plot = region_plot)
+p
+
+```
+
+
+```{r, echo=FALSE}
+
+SW_trend_plot = function(FV_SW,sector_plot){
+  if (sector_plot == "trn_pass") {
+  FV_SW=FV_SW[iso=="DEU" & vehicle_type =="Large Car and SUV",]
+  } else if (sector_plot =="trn_freight"){
+  FV_SW=FV_SW[iso=="DEU" & vehicle_type =="Truck (16-32t)",]
+  } else if (sector_plot =="trn_aviation_intl"){
+  FV_SW=FV_SW[iso=="DEU" & subsector_L3 =="International Aviation",]
+  } else if (sector_plot =="trn_shipping_intl"){
+  FV_SW=FV_SW[iso=="DEU" & subsector_L3 =="International Ship",]
+  }
+  FV_SW[,type:=ifelse(technology=="Liquids", "Conventional ICE (Liquid fuels)",NA)]
+  FV_SW[,type:=ifelse(technology=="NG", "Natural Gas ICE",type)]
+  FV_SW[,type:=ifelse(technology=="BEV", "Alternative fuels: BEV",type)]
+  FV_SW[,type:=ifelse(technology=="FCEV", "Alternative fuels: FCEV",type)]
+  FV_SW[,type:=ifelse(technology=="Hybrid Liquids", "Unconventional ICE (Hybrid)",type)]
+  
+  p=ggplot()+
+    geom_line(data=FV_SW%>%filter(year>= min(years_plot), year<=max(years_plot)),aes(x=year,y=sw,group=type,color=type),alpha = 0.8,size=1.5)+
+    theme_light()+
+    facet_wrap(~EDGE_scenario)+
+    theme(axis.text.x = element_text(angle = 90, hjust = 1),
+          axis.text = element_text(size=13),
+          title = element_text(size=14),
+          legend.text = element_text(size=13))+
+    scale_x_continuous(breaks=years_plot)+
+    xlab("Year")+
+    ylab ("Preference factors tech. types [-]")+
+    ggtitle(paste0("Preference factors trend for tech. types for ", sector_plot, " [-]"))+
+    theme(strip.text.x = element_text(size=13,color="black"),
+          strip.background = element_rect(fill="white",color = "black"))+
+    scale_color_discrete(name="Technology type")
+  return(p)
+  }
+
+p=SW_trend_plot(FV_SW=sw_tech,sector_plot)
+p
+```
+
+