EDGETransportComparison.Rmd

title: "Compare scenarios Transport"
output:
  pdf_document: default
  html_document:
    df_print: paged
classoption: landscape
knitr::opts_chunk$set(dev = 'pdf')
require(ggplot2)
require(moinput)
require(rmndt)
require(quitte)
library(lucode)
library(magpie)
library(quitte)
library(cowplot)
# Set RDS files path
EJmode_all = readRDS("EJmode_all.RDS")
EJroad_all = readRDS("EJroad_all.RDS")
fleet_all = readRDS("fleet_all.RDS")
salescomp_all = readRDS("salescomp_all.RDS")
ESmodecap_all = readRDS("ESmodecap_all.RDS")
CO2km_int_newsales_all = readRDS("CO2km_int_newsales_all.RDS")
EJpass_all = readRDS("EJfuelsPass_all.RDS")
EJfrgt_all = readRDS("EJfuelsFrgt_all.RDS")
emidem_all = readRDS("emidem_all.RDS")

setConfig(forcecache=T)

cols <- c("NG" = "#d11141",
          "Liquids" = "#8c8c8c",
          "Hybrid Liquids" = "#ffc425",
          "Hybrid Electric" = "#f37735",
          "BEV" = "#00b159",
          "Electricity" = "#00b159",
          "FCEV" = "#00aedb",
          "Hydrogen" = "#00aedb",
          "Biodiesel" = "#66a182",
          "Synfuel" = "orchid",
          "Oil" = "#2e4057",
          "fuel price pkm" = "#edae49",
          "Operating costs registration and insurance" = "#8d96a3",
          "Operating costs maintenance" = "#00798c",
          "Capital cost" = "#d1495b",
          "International Aviation" = "#9acd32",
          "Domestic Aviation" = "#7cfc00",
          "Bus" = "#32cd32",
          "Passenger Rail" = "#2e8b57",
          "Freight Rail" = "#ee4000",
          "Trucks" = "#ff6a6a",
          "International Shipping" = "#cd2626",
          "Domestic Shipping" = "#ff4040",
          "Shipping" = "#ff4040",
          "Truck" = "#ff7f50",
          "Trucks (<3.5t)" = "#ff7f50",
          "Trucks (3.5t-16)" = "#8b0000",
          "Trucks (>16)" = "#fa8072",
          "Motorbikes" = "#1874cd", #"dodgerblue3",
          "Small Cars" = "#87cefa",
          "Large Cars" = "#6495ed",
          "Van" = " 	#40e0d0",
          "LDV" = "#00bfff",
          "Non motorized" = "#da70d6",
          "Freight"="#ff0000",
          "Freight (Inland)" = "#cd5555",
          "Pass non LDV" = "#6b8e23",
          "Pass" = "#66cdaa",
          "Pass non LDV (Domestic)" = "#54ff9f",
          "refined liquids enduse" = "#8c8c8c",
          "FE|Transport|Hydrogen" = "#00aedb",
          "FE|Transport|NG" = "#d11141",
          "FE|Transport|Liquids" = "#8c8c8c",
          "FE|Transport|Electricity" = "#00b159",
          "FE|Transport" = "#1e90ff",
          "FE|Buildings" = "#d2b48c",
          "FE|Industry" = "#919191",
          "Electricity_push" = "#00b159",
          "ElecEra" = "#00b159",
          "ElecEraWise" = "#68c6a4",
          "HydrHype" = "#00aedb",
          "HydrHypeWise" = "#o3878f",
          "Hydrogen_push" = "#00aedb",
          "Smart_lifestyles_Electricity_push" = "#68c6a4",
          # "Smart_lyfestiles_Electricity_push" = "#03878f", ##maybe "#o3878f"
          "Conservative_liquids" = "#113245",
          "ConvCase" = "#113245",
          "ConvCaseWise" = "#d11141",
          "Emi|CO2|Transport|Demand" = "#113245",
          "Emi|CO2|Industry|Gross" = "#919191",
          "Emi|CO2|Buildings|Direct" = "#d2b48c",
          "Emi|CO2|Energy|Supply|Gross" = "#f2b531",
          "Emi|CO2|CDR|BECCS" = "#ed5958",
          "Emi|CO2|Land-Use Change" = "#66a182",
          "Cons. + Synfuels" = "orchid",
          "Ctax_Conservative" = "#d11141")

legend_ord_modes <- c("Freight Rail", "Truck", "Shipping", "International Shipping", "Domestic Shipping",  "Trucks",
                "Motorbikes", "Small Cars", "Large Cars", "Van",
                "International Aviation", "Domestic Aviation","Bus", "Passenger Rail",
                "Freight", "LDV", "Pass non LDV", "Freight (Inland)", "Pass non LDV (Domestic)", "Non motorized")

legend_ord_fuels <- c("BEV", "Electricity", "Hybrid Electric", "FCEV", "Hydrogen", "Hybrid Liquids", "Liquids", "Oil", "Biodiesel", "Synfuel", "NG")

legend_ord_costs <- c("fuel price pkm", "Operating costs registration and insurance", "Operating costs maintenance", "Capital cost")

legend_ord_emissions <- c("Emi|CO2|Industry|Gross", "Emi|CO2|Buildings|Direct", "Emi|CO2|Transport|Demand", "Emi|CO2|Energy|Supply|Gross", "Emi|CO2|Land-Use Change","Emi|CO2|CDR|BECCS")

legend_ord = c(legend_ord_modes, legend_ord_fuels, legend_ord_costs)

## mapping for scenario names
mapping_scens = data.table(scenario = c("Budg1100_ConvCase", "Budg1100_ElecEra", "Budg1100_HydrHype", "Budg1100_SynSurge", "NDC_ConvCase", "Budg1100_ConvCaseWise", "Budg1100_ElecEraWise", "Budg1100_HydrHypeWise", "Budg1100_SynSurgeWise", "NDC_ConvCaseWise"), scen_name = c("ConvCase", "ElecEra", "HydrHype", "SynSurge", "Baseline", "ConvCaseWise", "ElecEraWise", "HydrHypeWise", "SynSurgeWise", "BaselineWise"))

regionplot = "EUR"

## Vintages
vintcomparisonpf = function(dt){
  dt = dt[year %in% c(2015, 2050, 2100)]
  plot = ggplot()+
    geom_bar(data = dt,
             aes(x=scenario, y=value, group=interaction(variable, technology),
                 fill = technology, width=.75), alpha = 0.5, position="stack", stat = "identity", width = 0.5)+
    geom_bar(data = dt,
             aes(x=scenario, y=value, group=interaction(variable, technology),
                 fill = technology, alpha = factor(alphaval),  width=.75), position="stack", stat = "identity", width = 0.5, color = "black", size=0.05)+
    guides(fill = guide_legend(reverse=TRUE))+
    theme_minimal()+
    facet_grid(year~region)+
    theme(axis.text.x = element_text(angle = 90, size=14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size=14),
          axis.title.y = element_text(size=14),
          title = element_text(size=14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          legend.text = element_text(size=14),
          strip.text = element_text(size=14),
          strip.background = element_rect(color = "grey"))+
    scale_alpha_discrete(breaks = c(1,0), name = "Status", labels = c("Vintages","New additions")) +
    guides(linetype=FALSE,
           fill=guide_legend(reverse=FALSE, title="Transport mode"))+
    scale_fill_manual(values = cols)+
    labs(y = "[million Veh]", x="", title = "LDV fleet")
    return(plot)
}

vintcomparisonpf(fleet_all)

salescompf = function(dt){

  plot = ggplot()+
    geom_bar(data = dt, aes(x=as.numeric(as.character(year)),y=shareFS1, group = technology, fill = technology), position = position_stack(), stat = "identity")+
    theme_minimal()+
    facet_grid(region ~ scenario)+
    scale_fill_manual("Technology", values = cols)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2015,2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90, vjust=0.5, hjust=1, size = 14),
          axis.text.y = element_text(size = 14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"))+
    labs(x = "", y = "[%]", title = "Market share of new LDV sales")
  return(plot)
}

salescompf(salescomp_all)

EJroadpf = function(dt){
  dt[, technology := factor(technology, levels = legend_ord)]
  dt = dt[year >= 2020]
  plotLDV = ggplot()+
    geom_area(data = dt[subsector_L1 == "trn_pass_road_LDV_4W"], aes(x=year, y=demand_EJ, group = technology, fill = technology), color = "black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "LDV Final Energy demand")+
    theme_minimal()+
    facet_grid(scenario~region)+
    scale_fill_manual("Technology", values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90, size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"))


   plotBus = ggplot()+
    geom_area(data = dt[subsector_L1 %in% c("trn_pass_road_bus_tmp_subsector_L1", "Bus_tmp_subsector_L1")], aes(x=year, y=demand_EJ, group = technology, fill = technology), color = "black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "Buses Final Energy demand")+
    theme_minimal()+
    facet_grid(scenario~region)+
    scale_fill_manual("Technology", values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90, size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"))


   plotTruck = ggplot()+
    geom_area(data = dt[subsector_L1 %in% c("trn_freight_road_tmp_subsector_L1")], aes(x=year, y=demand_EJ, group = technology, fill = technology), color = "black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "Trucks Final Energy demand")+
    theme_minimal()+
    facet_grid(scenario~region)+
    scale_fill_manual("Technology", values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90, size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"))

  return(plotlist = list(plotLDV = plotLDV, plotBus = plotBus, plotTruck = plotTruck))
}

EJroadpf(EJroad_all)

EJmodepf = function(dt){
  dt = dt[year >= 2020]
  plot = ggplot()+
    geom_area(data = dt, aes(x=year, y=demand_EJ, group = interaction(vehicle_type_plot,aggr_mode), fill = vehicle_type_plot), color = "black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "Total transport final energy demand")+
    theme_minimal()+
    facet_grid(scenario~region)+
    scale_fill_manual("Vehicle Type",values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020,2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90, size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          strip.background = element_rect(color = "grey"))
  return(plot)
}

EJmodepf(EJmode_all)

ESmodecappf = function(dt){
  dt[, vehicle_type_plot := factor(vehicle_type_plot, levels = legend_ord)]
  plot_frgt = ggplot()+
    geom_area(data = dt[mode == "freight" & year >= 2020], aes(x=year, y=cap_dem, group = vehicle_type_plot, fill = vehicle_type_plot), color="black", size=0.05, position= position_stack())+
    labs(x = "", y = "Energy Services demand [tkm/cap]")+
    theme_minimal()+
    facet_grid(scenario~region)+
    scale_fill_manual("Vehicle Type",values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020,2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90,  size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"),
          axis.line = element_line(size = 0.5, colour = "grey"))


  plot_pass = ggplot()+
    geom_area(data = dt[mode == "pass" & year >= 2020], aes(x=year, y=cap_dem, group = vehicle_type_plot, fill = vehicle_type_plot), color="black", size=0.05, position= position_stack())+
    labs(x = "", y = "Energy Services demand [pkm/cap]")+
    theme_minimal()+
    facet_grid(scenario~region)+
    scale_fill_manual("Vehicle Type",values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020,2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90,  size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"),
          axis.line = element_line(size = 0.5, colour = "grey"))


    return(list(plot_pass = plot_pass, plot_frgt = plot_frgt))
}


ESmodecappf(ESmodecap_all)

CO2km_int_newsalespf = function(dt){
  dt = dt[!is.na(gCO2_km_ave)]
  plot = ggplot()+
    geom_line(data = dt[year >= 2020], aes(x = year, y = gCO2_km_ave, group = scenario, color = scenario))+
    labs(title = expression(paste(CO["2"], " intensity of LDVs new additions")), y = expression(paste("[", gCO["2"], "/km]")), x = "")+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030, 2050, 2100))+
    theme_minimal()+
    facet_grid(~region)+
    theme(axis.text.x = element_text(angle = 90,  size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"),
          axis.line = element_line(size = 0.5, colour = "grey"))+
    guides(linetype = FALSE)
  return(plot)
}

CO2km_int_newsalespf(CO2km_int_newsales_all)
## passenger by fuel
EJfuels_pf = function(dt_p, dt_f){
  dt_p = dt_p[year >= 2020]
  dt_p = dt_p[, .(demand_EJ = sum(demand_EJ)), by = c("subtech", "year", "region", "scenario")]
  plotp = ggplot()+
    geom_area(data = dt_p, aes(x=year, y=demand_EJ, group = subtech, fill = subtech), color="black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "Passenger transport FE demand by fuel")+
    theme_minimal()+
    facet_grid(scenario~region)+
    scale_fill_manual("Technology",values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90,  size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"),
          axis.line = element_line(size = 0.5, colour = "grey"))

    dt_f = dt_f[year >= 2020]

  plotf_lo = ggplot()+
    geom_area(data = dt_f[sector == "trn_shipping_intl"], aes(x=year, y=demand_EJ, group = subtech, fill = subtech), color="black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "International freight FE demand by fuel")+
    theme_minimal()+
    facet_grid(scenario~region)+
    scale_fill_manual("Technology",values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90,  size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"),
          axis.line = element_line(size = 0.5, colour = "grey"))

  plotf_sm = ggplot()+
    geom_area(data = dt_f[sector == "trn_freight"], aes(x=year, y=demand_EJ, group = subtech, fill = subtech), color="black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "Short-medium freight FE demand by fuel")+
    theme_minimal()+
    facet_grid(scenario~region)+
    scale_fill_manual("Technology",values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90,  size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"),
          axis.line = element_line(size = 0.5, colour = "grey"))

  plot = list(plotf_lo = plotf_lo, plotf_sm = plotf_sm, plotp = plotp)
  return(plot)
}

EJfuels_pf(dt_p = EJpass_all, dt_f = EJfrgt_all)
emidem_pf = function(dt){
 dt[, scenario := as.character(scenario)]
  plot = ggplot()+
    geom_line(data = dt, aes(x = year, y = value, group = scenario, color = scenario))+
    labs(x = "", y = "CO2 emissions [Mt/CO2]", title = "Emissions from transport demand")+
    theme_minimal()+
    facet_grid(~region)+
    theme(axis.text.x = element_text(angle = 90,  size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"),
          axis.line = element_line(size = 0.5, colour = "grey"))

  return(plot)
}

emidem_pf(emidem_all)

vintcomparison_regi_pf = function(dt, mapping_scens, rp){
  dt = dt[year %in% c(2020, 2030, 2050) & region == rp]
  ## apply more esthetic scenario names
  dt = merge(dt, mapping_scens, by = "scenario")
  ## if scenario name not available in the mapping, apply scenario extended name (e.g. NPi, Budget950... are not recorded in the mapping)
  dt[is.na(scen_name), scen_name := scenario]
  ## define the maximum on the y-axis depending on the maximum value across scenarios
  dt[, sum := sum(value), by = c("year", "scenario")]
  dt[, maxval := max(sum)]
  ## select a random scenario for 2020 and attribute the label "Historical"
  p1 = ggplot()+
    geom_bar(data = dt[year == 2020 & scenario == unique(dt$scenario)[1]][, scen_name := "Historical"],
             aes(x=scen_name, y=value, group= technology,
                 fill = technology, width=.75), position="stack", stat = "identity", width = 0.5, alpha = 0.9)+
    theme_minimal()+
    ylim(0,unique(dt[, maxval]))+
    facet_wrap(~ year, nrow = 1)+
    theme(axis.text.x = element_text(angle = 90, size=14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size=14),
          axis.title.y = element_text(size=14),
          title = element_text(size=14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          legend.text = element_text(size=14),
          strip.text = element_text(size=14),
          strip.background = element_rect(color = "grey"),
          legend.position = "none")+
    scale_fill_manual(values = cols)+
    labs(y = "[million Veh]", x="")

  p2 = ggplot()+
    geom_bar(data = dt[year != 2020],
             aes(x=scen_name, y=value, group=interaction(variable, technology),
                 fill = technology, width=.75), position="stack", stat = "identity", width = 0.5, alpha = 0.9)+
    theme_minimal()+
    ylim(0,unique(dt[, maxval]))+
    facet_wrap(~ year, nrow = 1)+
    theme(axis.text.x = element_text(angle = 90, size=14, vjust=0.5, hjust=1),
          axis.text.y = element_blank(),
          axis.title.y = element_text(size=14),
          title = element_text(size=14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          legend.text = element_text(size=14),
          strip.text = element_text(size=14),
          strip.background = element_rect(color = "grey"))+
    guides(fill=guide_legend(title="Powertrain"))+
    scale_fill_manual(values = cols)+
    labs(y = "", x="")

  plot = plot_grid(p1, p2, align = "h", ncol = 2, rel_widths = c(0.15,0.85))

  return(plot)
}

p = vintcomparison_regi_pf(fleet_all, mapping_scens, rp = regionplot)

p

aspect_ratio <- 1.5
height <- 6
ggsave("pvint.png", p, dpi=500, height = height , width = height * aspect_ratio)
salescom_regi_pf = function(dt, rp){

  plot = ggplot()+
    geom_area(data = dt[region == rp], aes(x=as.numeric(as.character(year)), y = shareFS1, group = technology, fill = technology), position = position_fill())+
    theme_minimal()+
    facet_wrap( ~ scenario, nrow = 1)+
    scale_fill_manual("Technology", values = cols)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2015,2030,2050, 2100))+
    scale_y_continuous(labels = scales::percent)+
    theme(axis.text.x = element_text(angle = 90, vjust=0.5, hjust=1, size = 14),
          axis.text.y = element_text(size = 14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"))+
    labs(x = "", y = "[%]", title = "Market share of new LDV sales")
  return(plot)
}

p = salescom_regi_pf(salescomp_all, rp = regionplot)

p

aspect_ratio <- 2
height <- 5
ggsave("psales.png", p, dpi=500, height = height , width = height * aspect_ratio)


CO2km_int_regi_newsalespf = function(dt, rp){
  dt = dt[!is.na(gCO2_km_ave)]
  if (rp == "EUR"){
   ## add historical values
  historical_values = data.table(year = c(2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018), emi = c(159, 157, 145, 140, 137, 132, 128, 124, 120, 119, 119, 120))

  targets = data.table(name = c("2021 target", "2025 target", "2030 target"), value = c(95, 95*(1-0.15), 95*(1-0.37)))

   plot = ggplot()+
    geom_line(data = dt[year >= 2020 & region == rp], aes(x = year, y = gCO2_km_ave, group = scenario, color = scenario))+
    geom_point(data = historical_values, aes(x = year, y = emi), color = "grey20")+
    geom_hline(data = targets, aes(yintercept = value, linetype = name), color = "grey20", size=0.1)+
    geom_text(data = targets, aes(y = value+5, x = c(2025, 2030, 2035), label = name), size = 5)+
    labs(title = expression(paste(CO["2"], " intensity of LDVs new additions")), y = expression(paste("[", gCO["2"], "/km]")), x = "")+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030, 2050, 2100))+
    theme_minimal()+
    theme(axis.text.x = element_text(angle = 90,  size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"),
          axis.line = element_line(size = 0.5, colour = "grey"))+
    guides(linetype = FALSE)

  } else {
  ## historical values are not available
  plot = ggplot()+
    geom_line(data = dt[year >= 2020 & region == rp], aes(x = year, y = gCO2_km_ave, group = scenario, color = scenario))+
    labs(title = expression(paste(CO["2"], " intensity of LDVs new additions")), y = expression(paste("[", gCO["2"], "/km]")), x = "")+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030, 2050, 2100))+
    theme_minimal()+
    theme(axis.text.x = element_text(angle = 90,  size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"),
          axis.line = element_line(size = 0.5, colour = "grey"))+
    guides(linetype = FALSE)

   }
  return(plot)
}

p = CO2km_int_regi_newsalespf(CO2km_int_newsales_all, rp = regionplot)

p

aspect_ratio <- 1.5
height <- 6
ggsave("pCO2int.png", p, dpi=500, height = height , width = height * aspect_ratio)
EJroad_regi_pf = function(dt, rp){
  dt[, technology := factor(technology, levels = legend_ord)]
  dt = dt[year >= 2020]
  plotLDV = ggplot()+
    geom_area(data = dt[subsector_L1 == "trn_pass_road_LDV_4W" & region == rp], aes(x=year, y=demand_EJ, group = technology, fill = technology), color = "black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "LDV Final Energy demand")+
    theme_minimal()+
    facet_wrap(~scenario, nrow = 1)+
    scale_fill_manual("Technology", values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90, size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"))


   plotBus = ggplot()+
    geom_area(data = dt[subsector_L1 %in% c("trn_pass_road_bus_tmp_subsector_L1", "Bus_tmp_subsector_L1") & region == rp], aes(x=year, y=demand_EJ, group = technology, fill = technology), color = "black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "Buses Final Energy demand")+
    theme_minimal()+
    facet_wrap(~scenario, nrow = 1)+
    scale_fill_manual("Technology", values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90, size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"))


   plotTruck = ggplot()+
    geom_area(data = dt[subsector_L1 %in% c("trn_freight_road_tmp_subsector_L1") & region == rp], aes(x=year, y=demand_EJ, group = technology, fill = technology), color = "black", size=0.05, position= position_stack())+
    labs(x = "", y = "[EJ]", title = "Trucks Final Energy demand")+
    theme_minimal()+
    facet_wrap(~scenario, nrow = 1)+
    scale_fill_manual("Technology", values = cols, breaks=legend_ord)+
    expand_limits(y = c(0,1))+
    scale_x_continuous(breaks = c(2020, 2030,2050, 2100))+
    theme(axis.text.x = element_text(angle = 90, size = 14, vjust=0.5, hjust=1),
          axis.text.y = element_text(size = 14),
          axis.title = element_text(size = 14),
          axis.line = element_line(size = 0.5, colour = "grey"),
          title = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          strip.background = element_rect(color = "grey"))

  return(plotlist = list(plotLDV = plotLDV, plotBus = plotBus, plotTruck = plotTruck))
}

plist = EJroad_regi_pf(EJroad_all, rp = regionplot)

plist

pLDV = plist[["plotLDV"]]

pBus = plist[["plotBus"]]

pTruck = plist[["plotTruck"]]

aspect_ratio <- 1.5
height <- 6
ggsave("pLDV.png", pLDV, dpi=500, height = height , width = height * aspect_ratio)
ggsave("pBus.png", pBus, dpi=500, height = height , width = height * aspect_ratio)
ggsave("pTruck.png", pTruck, dpi=500, height = height , width = height * aspect_ratio)