From c1fba3e9def0728c784f301c11f37b7f99b3cae0 Mon Sep 17 00:00:00 2001
From: Marianna Rottoli <marianna.rottoli@mail.polimi.it>
Date: Tue, 31 Mar 2020 17:05:30 +0200
Subject: [PATCH] Updated EDGE single output script. Added EDGE comparison
output script.
---
scripts/output/comparison/EDGEcomparison.R | 207 +++++++++++++++++
.../EDGETransportComparison.Rmd | 218 ++++++++++++++++++
.../EDGETransportReport.Rmd | 63 +++--
3 files changed, 474 insertions(+), 14 deletions(-)
create mode 100644 scripts/output/comparison/EDGEcomparison.R
create mode 100644 scripts/output/comparison/notebook_templates/EDGETransportComparison.Rmd
diff --git a/scripts/output/comparison/EDGEcomparison.R b/scripts/output/comparison/EDGEcomparison.R
new file mode 100644
index 0000000..2038639
--- /dev/null
+++ b/scripts/output/comparison/EDGEcomparison.R
@@ -0,0 +1,207 @@
+# | (C) 2006-2019 Potsdam Institute for Climate Impact Research (PIK)
+# | authors, and contributors see CITATION.cff file. This file is part
+# | of REMIND and licensed under AGPL-3.0-or-later. Under Section 7 of
+# | AGPL-3.0, you are granted additional permissions described in the
+# | REMIND License Exception, version 1.0 (see LICENSE file).
+# | Contact: remind@pik-potsdam.de
+
+require(rmarkdown)
+require(lucode)
+require(quitte)
+require(data.table)
+require(rmndt)
+
+
+if(!exists("source_include")) {
+ ## Define arguments that can be read from command line
+ readArgs("outputdirs")
+}
+
+## Check if the output is EDGE based, otherwise remove the output directory from the list of compared output folders
+for (outputdir in outputdirs) {
+ load(file.path(outputdir, "config.Rdata"))
+ if(cfg$gms$transport != "edge_esm"){
+ print(paste0("The scenario ", outputdir, " is not EDGE-based and will be excluded from the reporting"))
+ outputdirs = outputdirs[outputdirs != outputdir]
+ }
+}
+
+
+gdx_name = "fulldata.gdx"
+
+
+emi_all = NULL
+salescomp_all = NULL
+fleet_all = NULL
+EJLDV_all = NULL
+EJmode_all = NULL
+
+
+
+
+scenNames <- getScenNames(outputdirs)
+EDGEdata_path <- path(outputdirs, paste("EDGE-T/"))
+gdx_path <- path(outputdirs,gdx_name)
+scenNames <- getScenNames(outputdirs)
+
+names(gdx_path) <- scenNames
+names(EDGEdata_path) <- scenNames
+
+REMIND2ISO_MAPPING <- fread("config/regionmappingH12.csv")[, .(iso = CountryCode, region = RegionCode)]
+
+
+
+SalesFun = function(shares_LDV, newcomp, sharesVS1){
+ ## I need the total demand for each region to get the average composition in Europe (sales are on a country level)
+ ## First I calculate the total demand for new sales using the shares on FV level (in newcomp) and on VS1 level
+ newcomp = merge(newcomp, sharesVS1[,.(shareVS1 = share, iso, year, vehicle_type, subsector_L1)], all.x=TRUE, by = c("iso", "year", "vehicle_type", "subsector_L1"))
+ newcomp[, newdem := totdem*sharetech_new*shareVS1]
+ newcomp = newcomp[,.(value = sum(newdem)), by = c("iso", "year", "subsector_L1")]
+ ## I have to interpolate in time the sales nto to loose the sales composition annual values
+ newcomp=approx_dt(dt=newcomp, unique(shares_LDV$year),
+ xcol= "year",
+ ycol = "value",
+ idxcols=c("iso","subsector_L1"),
+ extrapolate=T)
+
+ setnames(newcomp, new = "newdem", old = "value")
+
+ ## I calculate the sales composition (disrespective to the vehicle type)
+ shares_LDV = unique(shares_LDV[,c("iso","year", "technology", "shareFS1")])
+ shares_LDV <- shares_LDV[,.(shareFS1=sum(shareFS1)),by=c("iso","technology","year")]
+
+ ## I calculate the weighted regional sales (depending on the total volume of sales per country in each region)
+ shares_LDV = merge(shares_LDV, newcomp)
+ shares_LDV = merge(shares_LDV, REMIND2ISO_MAPPING, by = "iso")
+ shares_LDV[, demfuel := shareFS1*newdem, by = c("year", "iso", "technology")]
+ shares_LDV = shares_LDV[, .(demfuel = sum(demfuel)), by = c("year", "region", "technology")]
+ shares_LDV[, shareFS1 := demfuel/sum(demfuel), by = c("year", "region")]
+
+ ## plot features
+ shares_LDV[, technology := factor(technology, levels = c("BEV", "Hybrid Electric", "FCEV", "Hybrid Liquids", "Liquids", "NG"))]
+
+ return(shares_LDV)
+}
+
+
+fleetFun = function(vintcomp, newcomp, sharesVS1){
+ vintcomp = vintcomp[,.(totdem, iso, subsector_L1, year, technology,vehicle_type, sector, sharetech_vint)]
+ newcomp = newcomp[,.(iso, subsector_L1, year, technology,vehicle_type, sector, sharetech_new)]
+
+ allfleet = merge(newcomp, vintcomp, all =TRUE, by = c("iso", "sector", "subsector_L1", "vehicle_type", "technology", "year"))
+ allfleet = merge(allfleet, sharesVS1[,.(shareVS1 = share, iso, year, vehicle_type, subsector_L1)], all.x=TRUE, by = c("iso", "year", "vehicle_type", "subsector_L1"))
+ allfleet[,vintdem:=totdem*sharetech_vint*shareVS1]
+ allfleet[,newdem:=totdem*sharetech_new*shareVS1]
+ allfleet=melt(allfleet, id.vars = c("iso", "sector", "subsector_L1", "vehicle_type", "technology",
+ "year"), measure.vars = c("vintdem", "newdem"))
+ allfleet[,alpha:=ifelse(variable == "vintdem", 0, 1)]
+
+ load_factor = 1.5
+ annual_mileage = 15000
+ allfleet = allfleet[,.(value = sum(value/load_factor/annual_mileage)), by = c("iso", "technology", "variable", "year")]
+
+ allfleet = merge(allfleet, REMIND2ISO_MAPPING, by = "iso")
+ allfleet = allfleet[,.(value = sum(value)), by = c("region", "technology", "variable", "year")]
+ allfleet[,alphaval := ifelse(variable =="vintdem", 1,0)]
+ allfleet[, technology := factor(technology, levels = c("BEV", "Hybrid Electric", "FCEV", "Hybrid Liquids", "Liquids", "NG"))]
+
+ return(allfleet)
+}
+
+
+EJLDVFun <- function(demandEJ){
+ demandEJ = demandEJ[subsector_L1 == "trn_pass_road_LDV_4W",]
+ demandEJ <- demandEJ[, c("sector", "subsector_L3", "subsector_L2", "subsector_L1", "vehicle_type", "technology", "iso", "year", "demand_EJ")]
+
+ demandEJ = merge(demandEJ, REMIND2ISO_MAPPING, by = "iso")
+ demandEJ[technology == "Hybrid Liquids", technology := "Liquids"]
+ demandEJ[technology == "FCEV", technology := "Hydrogen"]
+ demandEJ[technology == "BEV", technology := "Electricity"]
+ demandEJ = demandEJ[, .(demand_EJ = sum(demand_EJ)), by = c("region", "year", "technology")]
+ return(demandEJ)
+
+}
+
+
+EJmodeFun = function(demandEJ){
+
+ demandEJ[, aggr_mode := ifelse(subsector_L2 == "trn_pass_road_LDV", "LDV", NA)]
+ demandEJ[, aggr_mode := ifelse(subsector_L3 %in% c("Passenger Rail", "HSR", "International Aviation", "Domestic Aviation"), "Pass non LDV", aggr_mode)]
+ demandEJ[, aggr_mode := ifelse(subsector_L2 %in% c("trn_pass_road_bus", "Bus"), "Pass non LDV", aggr_mode)]
+ demandEJ[, aggr_mode := ifelse(is.na(aggr_mode), "Freight", aggr_mode)]
+ demandEJ[, veh := ifelse(grepl("Large|SUV|Midsize|Multipurpose Vehicle|Van|3W Rural", vehicle_type), "Large Cars", NA)]
+ demandEJ[, veh := ifelse(grepl("Subcompact|Compact|Mini|Three-Wheeler", vehicle_type), "Small Cars", veh)]
+ demandEJ[, veh := ifelse(grepl("Motorcycle|Moped|Scooter", vehicle_type), "Motorbikes", veh)]
+ demandEJ[, veh := ifelse(grepl("bus|Bus", vehicle_type), "Bus", veh)]
+ demandEJ[, veh := ifelse(grepl("Truck", vehicle_type) & vehicle_type != "Light Truck and SUV", "Truck", veh)]
+ demandEJ[, veh := ifelse(grepl("Freight Rail_tmp_vehicletype", vehicle_type), "Freight Rail", veh)]
+ demandEJ[, veh := ifelse(grepl("Passenger Rail|HSR", vehicle_type), "Passenger Rail", veh)]
+ demandEJ[, veh := ifelse(subsector_L3 == "Domestic Ship", "Domestic Shipping", veh)]
+ demandEJ[, veh := ifelse(subsector_L3 == "International Ship", "International Shipping", veh)]
+ demandEJ[, veh := ifelse(subsector_L3 == "Domestic Aviation", subsector_L3, veh)]
+ demandEJ[, veh := ifelse(subsector_L3 == "International Aviation", subsector_L3, veh)]
+ demandEJ[, veh := ifelse(is.na(veh), vehicle_type, veh)]
+ demandEJ = demandEJ[,.(demand_EJ = sum(demand_EJ)), by = c("iso", "year", "aggr_mode", "veh")]
+
+ demandEJ[, vehicle_type_plot := factor(veh, levels = c("LDV","Freight Rail", "Truck","Domestic Shipping", "International Shipping",
+ "Motorbikes", "Small Cars", "Large Cars", "Van",
+ "Domestic Aviation", "International Aviation", "Bus", "Passenger Rail",
+ "Freight", "Freight (Inland)", "Pass non LDV", "Pass non LDV (Domestic)"))]
+
+ demandEJ = merge(demandEJ, REMIND2ISO_MAPPING, by = "iso")
+ demandEJ = demandEJ[,.(demand_EJ= sum(demand_EJ)), by = c("region", "year", "vehicle_type_plot", "aggr_mode")]
+
+
+ return(demandEJ)
+}
+
+
+for (outputdir in outputdirs) {
+ ## load mif file
+ name_mif = list.files(path = outputdir, pattern = "REMIND_generic", full.names = F)
+ name_mif = name_mif[!grepl("withoutPlu", name_mif)]
+ miffile <- as.data.table(read.quitte(paste0(outputdir, "/", name_mif)))
+
+ ## load RDS files
+ sharesVS1 = readRDS(paste0(outputdir, "/EDGE-T/" ,"shares.RDS"))[["VS1_shares"]]
+ newcomp = readRDS(paste0(outputdir, "/EDGE-T/" ,"newcomp.RDS"))
+ vintcomp = readRDS(paste0(outputdir, "/EDGE-T/" ,"vintcomp.RDS"))
+ shares_LDV = readRDS(paste0(outputdir, "/EDGE-T/" ,"annual_sales.RDS"))
+ demandEJ = readRDS(paste0(outputdir, "/EDGE-T/" , "demandF_plot_EJ.RDS"))
+
+ ## calculate sales
+ salescomp = SalesFun(shares_LDV, newcomp, sharesVS1)
+ ## calculate fleet compositons
+ fleet = fleetFun(vintcomp, newcomp, sharesVS1)
+ ## calculate EJ from LDVs by technology
+ EJLDV = EJLDVFun(demandEJ)
+ ## calculate FE demand by mode
+ EJmode = EJmodeFun(demandEJ)
+
+ ## add scenario dimension to the results
+ fleet[, scenario := as.character(unique(miffile$scenario))]
+ salescomp[, scenario := unique(miffile$scenario)]
+ EJLDV[, scenario := as.character(unique(miffile$scenario))]
+ EJmode[, scenario := as.character(unique(miffile$scenario))]
+
+ ## rbind scenarios
+ salescomp_all = rbind(salescomp_all, salescomp)
+ fleet_all = rbind(fleet_all, fleet)
+ EJLDV_all = rbind(EJLDV_all, EJLDV)
+ EJmode_all = rbind(EJmode_all, EJmode)
+
+}
+
+dir.create("output/comparerunEDGE", showWarnings = FALSE)
+
+outdir = "output/comparerunEDGE/"
+md_template = "EDGETransportComparison.Rmd"
+
+saveRDS(EJmode_all, "output/comparerunEDGE/EJmode_all.RDS")
+saveRDS(salescomp_all, "output/comparerunEDGE/salescomp_all.RDS")
+saveRDS(fleet_all, "output/comparerunEDGE/fleet_all.RDS")
+saveRDS(EJLDV_all, "output/comparerunEDGE/EJLDV_all.RDS")
+
+file.copy(file.path("./scripts/output/comparison/notebook_templates", md_template), outdir)
+rmarkdown::render(path(outdir, md_template), output_format="pdf_document")
+
diff --git a/scripts/output/comparison/notebook_templates/EDGETransportComparison.Rmd b/scripts/output/comparison/notebook_templates/EDGETransportComparison.Rmd
new file mode 100644
index 0000000..2fbe9e6
--- /dev/null
+++ b/scripts/output/comparison/notebook_templates/EDGETransportComparison.Rmd
@@ -0,0 +1,218 @@
+---
+title: "Compare scenarios Transport"
+output:
+ html_document:
+ df_print: paged
+---
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+require(ggplot2)
+require(moinput)
+require(rmndt)
+require(quitte)
+library(lucode)
+library(magpie)
+library(quitte)
+```
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+# Set mif path
+EJmode_all = readRDS("EJmode_all.RDS")
+EJLDV_all = readRDS("EJLDV_all.RDS")
+fleet_all = readRDS("fleet_all.RDS")
+salescomp_all = readRDS("salescomp_all.RDS")
+
+iso_plot = c("DEU", "USA", "IND", "CHN")
+
+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",
+ "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",
+ "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",
+ "Hydrogen_push" = "#00aedb",
+ "Smart_lifestyles_Electricity_push" = "#68c6a4",
+ # "Smart_lyfestiles_Electricity_push" = "#03878f", ##maybe "#o3878f"
+ "Conservative_liquids" = "#113245",
+ "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", "Trucks (<3.5t)", "Trucks (3.5t-16)", "Trucks (>16)", "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", "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)
+
+
+
+
+```
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+
+```
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+## 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")+
+ guides(fill = guide_legend(reverse=TRUE))+
+ theme_minimal()+
+ facet_grid(year~region)+
+ theme(axis.text.x = element_text(angle = 90, size=8, vjust=0.5, hjust=1),
+ title = element_text(size=8),
+ axis.line = element_line(size = 0.5, colour = "grey"),
+ # legend.text = element_text(size=8),
+ strip.text = element_text(size=8),
+ 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)
+## Sales composition
+
+
+```
+
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+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 = 8),
+ axis.text.y = element_text(size=8),
+ axis.line = element_line(size = 0.5, colour = "grey"),
+ axis.title = element_text(size = 8),
+ title = element_text(size = 8),
+ # legend.text = element_text(8),
+ # legend.title = element_text(size = 8),
+ strip.text = element_text(size=8),
+ strip.background = element_rect(color = "grey"))+
+ labs(x = "", y = "[%]", title = "Market share of new LDV sales")
+ return(plot)
+}
+
+salescompf(salescomp_all)
+
+```
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+EJLDVpf = function(dt){
+ dt[, technology := factor(technology, levels = legend_ord)]
+ plot = ggplot()+
+ geom_area(data = dt, aes(x=year, y=demand_EJ, group = technology, fill = technology), color="black",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)+
+ theme(axis.text.x = element_text(angle = 90, size = 8, vjust=0.5, hjust=1),
+ axis.text.y = element_text(size = 8),
+ axis.title = element_text(size = 8),
+ axis.line = element_line(size = 0.5, colour = "grey"),
+ title = element_text(size = 8),
+ # legend.text = element_text(size =8),
+ # legend.title = element_text(size = 8),
+ strip.text = element_text(size=8),
+ strip.background = element_rect(color = "grey"))
+ return(plot)
+}
+
+EJLDVpf(EJLDV_all)
+
+```
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+EJmodeFUn = function(dt){
+ 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", 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)+
+ theme(axis.text.x = element_text(angle = 90, size = 8, vjust=0.5, hjust=1),
+ axis.text.y = element_text(size=8),
+ axis.title = element_text(size = 8),
+ title = element_text(size = 8),
+ # legend.text = element_text(size = 8),
+ # legend.title = element_text(size =8),
+ strip.text = element_text(size=8),
+ axis.line = element_line(size = 0.5, colour = "grey"),
+ strip.background = element_rect(color = "grey"))
+ return(plot)
+}
+
+EJmodeFUn(EJmode_all)
+
+```
\ No newline at end of file
diff --git a/scripts/output/single/notebook_templates/EDGETransportReport.Rmd b/scripts/output/single/notebook_templates/EDGETransportReport.Rmd
index 6874f33..6dbf3cd 100644
--- a/scripts/output/single/notebook_templates/EDGETransportReport.Rmd
+++ b/scripts/output/single/notebook_templates/EDGETransportReport.Rmd
@@ -205,26 +205,61 @@ intcompplotf(EF_shares, shares$FV_shares, shares$VS1_shares)
# Sales of LDVs
```{r, echo=FALSE, warning=FALSE}
-salesplot = function(annual_sales){
- annual_sales = unique(annual_sales[,c("iso","year", "technology", "shareFS1")])
- annual_sales <- annual_sales[,.(shareFS1=sum(shareFS1)),by=c("iso","technology","year")]
-
- p = ggplot()+
- geom_bar(data = annual_sales[year <= 2100 & year >= 2020 & iso == iso_plot], aes(x=as.numeric(as.character(year)),y=shareFS1, group = technology, fill = technology), position = position_stack(), stat = "identity")+
+salesplot= function(shares_LDV, newcomp, sharesVS1){
+ ## I need the total demand for each region to get the average composition in Europe (sales are on a country level)
+ ## First I calculate the total demand for new sales using the shares on FV level (in newcomp) and on VS1 level
+ newcomp = merge(newcomp, sharesVS1[,.(shareVS1 = share, iso, year, vehicle_type, subsector_L1)], all.x=TRUE, by = c("iso", "year", "vehicle_type", "subsector_L1"))
+ newcomp[, newdem := totdem*sharetech_new*shareVS1]
+ newcomp = newcomp[,.(value = sum(newdem)), by = c("iso", "year", "subsector_L1")]
+ ## I have to interpolate in time the sales nto to loose the sales composition annual values
+ newcomp=approx_dt(dt=newcomp, unique(shares_LDV$year),
+ xcol= "year",
+ ycol = "value",
+ idxcols=c("iso","subsector_L1"),
+ extrapolate=T)
+
+ setnames(newcomp, new = "newdem", old = "value")
+
+ ## I calculate the sales composition (disrespective to the vehicle type)
+ shares_LDV = unique(shares_LDV[,c("iso","year", "technology", "shareFS1")])
+ shares_LDV <- shares_LDV[,.(shareFS1=sum(shareFS1)),by=c("iso","technology","year")]
+
+ ## I calculate the weighted regional sales (depending on the total volume of sales per country in each region)
+ shares_LDV = merge(shares_LDV, newcomp)
+ shares_LDV = merge(shares_LDV, REMIND2ISO_MAPPING, by = "iso")
+ shares_LDV[, demfuel := shareFS1*newdem, by = c("year", "iso", "technology")]
+ shares_LDV = shares_LDV[, .(demfuel = sum(demfuel)), by = c("year", "region", "technology")]
+ shares_LDV[, shareFS1 := demfuel/sum(demfuel), by = c("year", "region")]
+
+ ## plot features
+ shares_LDV[, technology := factor(technology, levels = c("BEV", "Hybrid Electric", "FCEV", "Hybrid Liquids", "Liquids", "NG"))]
+
+
+ plot = ggplot()+
+ geom_bar(data = shares_LDV, aes(x=as.numeric(as.character(year)),y=shareFS1, group = technology, fill = technology), position = position_stack(), stat = "identity")+
theme_minimal()+
+ facet_grid(~ region )+
scale_fill_manual("Technology", values = cols)+
expand_limits(y = c(0,1))+
- scale_x_continuous(breaks = c(2020, 2030, 2050, 2100))+
- theme(axis.text.x = element_text(angle = 90, vjust = +0.1),
- strip.background = element_rect(color = "grey"),
- legend.position = "none")+
- labs(x = "", y = "Market share on LDVs [%]", title = paste0("Sales composition, example of ", iso_plot))
-
- return(p)
+ scale_x_continuous(breaks = c(2015,2030,2050, 2100))+
+#+
+ theme(
+axis.text.x = element_text(angle = 90, vjust=0.5, hjust=1, size = 12),
+ axis.text.y = element_text(size=8),
+ axis.line = element_line(size = 0.5, colour = "grey"),
+ axis.title = element_text(size = 12),
+ title = element_text(size = 12),
+ # legend.text = element_text(12),
+ #legend.title = element_text(size = 12),
+ strip.text = element_text(size=12),
+ strip.background = element_rect(color = "grey")
+)+
+ labs(x = "", y = "[%]", title = "Market share of new LDV sales")
+ return(plot)
}
+salesplot(annual_sales, newcomp, shares$VS1_shares)
-salesplot(annual_sales)
```
# Final energy demand
--
GitLab