From 63c406acbd7e83c6d0b09e4740af21504734bc4b Mon Sep 17 00:00:00 2001
From: Marianna Rottoli <marianna.rottoli@mail.polimi.it>
Date: Fri, 20 Mar 2020 11:16:24 +0100
Subject: [PATCH] Reporting script updates and small bugfixes.

---
 .../EDGETransportReport.Rmd                   | 77 +++++++++++++------
 1 file changed, 52 insertions(+), 25 deletions(-)

diff --git a/scripts/output/single/notebook_templates/EDGETransportReport.Rmd b/scripts/output/single/notebook_templates/EDGETransportReport.Rmd
index 4286f64..6874f33 100644
--- a/scripts/output/single/notebook_templates/EDGETransportReport.Rmd
+++ b/scripts/output/single/notebook_templates/EDGETransportReport.Rmd
@@ -18,7 +18,10 @@ require(magclass)
 require(quitte)
 require(ggpubr)
 require(gridExtra)
-require(edgeTrpLib)
+#require(edgeTrpLib)
+require(devtools)
+load_all("../../../../../../edgetrplib/")
+
 ```
 
 
@@ -52,10 +55,11 @@ demand_ej <- readRDS(datapath(fname = "demandF_plot_EJ.RDS")) ## detailed final
 vintcomp <- readRDS(datapath(fname = "vintcomp.RDS"))
 newcomp <- readRDS(datapath(fname = "newcomp.RDS"))
 shares <- readRDS(datapath(fname = "shares.RDS"))
-inco_tech <- readRDS(datapath(fname = "inco_costs.RDS"))
+pref <- readRDS(datapath(fname = "pref.RDS"))
 EF_shares <- readRDS(datapath(fname = "EF_shares.RDS"))
 annual_sales <- readRDS(datapath(fname = "annual_sales.RDS"))
 mj_km_data <- readRDS(datapath(fname = "mj_km_data.RDS"))
+stations <- readRDS(datapath(fname = "stations.RDS"))
 
 name_mif = list.files(pattern = "REMIND_generic", full.names = F)
 name_mif = name_mif[!grepl("withoutPlu", name_mif)]
@@ -122,22 +126,47 @@ p
 
 ```{r, echo=FALSE, warning=FALSE}
 
-p=ggplot()+
-  geom_bar(data = inco_tech[iso == iso_plot & subsector_L1 == "trn_pass_road_LDV_4W" & vehicle_type == "Large Car and SUV" & year<=2100 & year>=2010], aes(x = as.character(year), y = pinco, group = technology, fill = technology), position = position_stack(), stat = "identity")+
+plotinconv = function(inco_tech, iso_plot, vehicle_type){
+  p=ggplot()+
+  geom_bar(data = inco_tech[iso == iso_plot & subsector_L1 == "trn_pass_road_LDV_4W" & vehicle_type == vehicle_type & year<=2100 & year>=2010], aes(x = as.character(year), y = value, group = logit_type, fill = logit_type), position = position_stack(), stat = "identity")+
   facet_grid(~technology)+
-    theme_minimal()+
-  scale_fill_manual(values = cols)+
+  theme_minimal()+
+  # scale_fill_manual(values = cols)+
   expand_limits(y = c(0,0.8))+
-      scale_x_discrete(breaks = c(2015,2050,2100))+
+  scale_x_discrete(breaks = c(2015,2050,2100))+
   theme(axis.text.x = element_text(angle = 90, vjust = +0.1),
-        legend.position = "none",
         strip.background = element_rect(color = "grey"))+
   labs(x = "", y = "Inconvenience cost [$/pkm]", title = paste0("Example of ", iso_plot))
 
-p
+return(p)
+}
+
+plotinconv(inco_tech = pref$FV_final_pref[subsector_L1 == "trn_pass_road_LDV_4W"], iso_plot = "DEU", vehicle_type = "Large Car and SUV")
+plotinconv(inco_tech = pref$FV_final_pref[subsector_L1 == "trn_pass_road_LDV_4W"], iso_plot = "USA", vehicle_type = "Large Car")
+plotinconv(inco_tech = pref$FV_final_pref[subsector_L1 == "trn_pass_road_LDV_4W"], iso_plot = "JPN", vehicle_type = "Large Car")
+plotinconv(inco_tech = pref$FV_final_pref[subsector_L1 == "trn_pass_road_LDV_4W"], iso_plot = "CHN", vehicle_type = "Large Car")
+plotinconv(inco_tech = pref$FV_final_pref[subsector_L1 == "trn_pass_road_LDV_4W"], iso_plot = "IND", vehicle_type = "Large Car")
 
 ```
 
+# Share of stations
+
+```{r, echo=FALSE, warning=FALSE}
+
+stationsplot = function(stations, iso_plot){
+  p = ggplot()+
+    geom_line(data = stations[iso == iso_plot], aes(x= year, y = fracst))+
+    facet_grid(~technology)+
+    theme_minimal()+
+    labs(title = paste0("Stations for each fuel, ", iso_plot))
+
+  return(p)
+}
+
+stationsplot(stations, iso_plot)
+```
+
+
 # Endogenous intensity for Liquids
 
 ```{r, echo=FALSE, message=FALSE, warning=FALSE}
@@ -181,11 +210,11 @@ salesplot = function(annual_sales){
   annual_sales <- annual_sales[,.(shareFS1=sum(shareFS1)),by=c("iso","technology","year")]
 
   p = ggplot()+
-    geom_bar(data = annual_sales[year<=2050  & year>=2015 & iso == iso_plot], aes(x=as.numeric(as.character(year)),y=shareFS1, group = technology, fill = technology), position = position_stack(), stat = "identity")+
+    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")+
     theme_minimal()+
     scale_fill_manual("Technology", values = cols)+
     expand_limits(y = c(0,1))+
-    scale_x_continuous(breaks = c(2015,2030,2050))+
+    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")+
@@ -210,13 +239,11 @@ demandEJplotf = function(demandEJ){
   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(vehicle_type %in% c("Truck (0-1t)", "Truck (0-3.5t)", "Truck (0-2.7t)", "Truck (0-2t)"), "Trucks (<3.5t)", NA)]
-  demandEJ[, veh := ifelse(vehicle_type %in% c("Truck (16-32t)", "Truck (3.5-16t)", "Truck (6-15t)", "Truck (4.5-12t)", "Truck (2.7-4.5t)", "Truck (4.5-15t)"), "Trucks (3.5t-16)", veh)]
-  demandEJ[, veh := ifelse(vehicle_type %in% c("Truck (>15t)", "Truck (16-32t)", "Truck (>32t)" ), "Trucks (>16)", veh)]
-  demandEJ[, veh := ifelse(grepl("Large|SUV|Midsize|Multipurpose Vehicle|Van|3W Rural", vehicle_type), "Large Cars", veh)]
+  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)]
@@ -226,12 +253,12 @@ demandEJplotf = function(demandEJ){
   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", "Trucks (<3.5t)", "Trucks (3.5t-16)", "Truck (>12t)", "Trucks (>16)", "Trucks","Domestic Shipping", "International Shipping",
+  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)"))]
 
-  legend_ord <- c("Freight Rail", "Trucks (<3.5t)", "Trucks (3.5t-16)", "Truck (>12t)", "Trucks (>16)", "International Shipping","Domestic Shipping",  "Trucks",
+  legend_ord <- c("Freight Rail", "Truck", "International Shipping","Domestic Shipping",
                   "Motorbikes", "Small Cars", "Large Cars", "Van",
                   "International Aviation", "Domestic Aviation","Bus", "Passenger Rail",
                   "Freight", "LDV", "Pass non LDV", "Freight (Inland)", "Pass non LDV (Domestic)")
@@ -313,31 +340,30 @@ demandpkmplotf = function(demandpkm){
   demandpkm[, aggr_mode := ifelse(sector %in% c("trn_freight", "trn_shipping_intl"), "Freight", aggr_mode)]
   demandpkm[, aggr_mode := ifelse(sector %in% c("trn_aviation_intl"), "Pass. non LDV", aggr_mode)]
   demandpkm[, aggr_mode := ifelse(subsector_L2 %in% c("trn_pass_road_bus", "HSR_tmp_subsector_L2", "Passenger Rail_tmp_subsector_L2", "Cycle_tmp_subsector_L2", "Walk_tmp_subsector_L2", "Domestic Aviation_tmp_subsector_L2", "Bus") | subsector_L1 %in% c("trn_pass_road_LDV_2W"), "Pass. non LDV", aggr_mode)]
-
-  demandpkm[, veh := ifelse(vehicle_type %in% c("Truck (0-1t)", "Truck (0-3.5t)"), "Trucks (<3.5t)", "Trucks (<3.5t)")]
-  demandpkm[, veh := ifelse(vehicle_type %in% c("Truck (16-32t)", "Truck (3.5-16t)", "Truck (6-15t)"), "Trucks (3.5t-16)", veh)]
-  demandpkm[, veh := ifelse(vehicle_type %in% c("Truck (>15t)", "Truck (16-32t)", "Truck (>32t)" ), "Trucks (>16)", veh)]
-  demandpkm[, veh := ifelse(grepl("Large|SUV|Midsize|Multipurpose Vehicle|Van|3W Rural", vehicle_type), "Large Cars", veh)]
+  demandpkm[, veh := ifelse(grepl("Large|SUV|Midsize|Multipurpose Vehicle|Van|3W Rural", vehicle_type), "Large Cars", NA)]
   demandpkm[, veh := ifelse(grepl("Subcompact|Compact|Mini|Three-Wheeler_tmp_vehicletype", vehicle_type), "Small Cars", veh)]
   demandpkm[, veh := ifelse(grepl("Motorcycle|Moped|Scooter", vehicle_type), "Motorbikes", veh)]
   demandpkm[, veh := ifelse(grepl("bus|Bus", vehicle_type), "Bus", veh)]
+  demandpkm[, veh := ifelse(grepl("Truck", vehicle_type) & vehicle_type != "Light Truck and SUV", "Truck", veh)]
   demandpkm[, veh := ifelse(subsector_L3 == "Domestic Aviation", "Domestic Aviation", veh)]
   demandpkm[, veh := ifelse(subsector_L3 == "International Aviation", "International Aviation", veh)]
   demandpkm[, veh := ifelse(grepl("Freight Rail", vehicle_type), "Freight Rail", veh)]
   demandpkm[, veh := ifelse(grepl("Passenger Rail|HSR", vehicle_type), "Passenger Rail", veh)]
   demandpkm[, veh := ifelse(grepl("Ship", vehicle_type), "Shipping", veh)]
   demandpkm[, veh := ifelse(grepl("Cycle|Walk", subsector_L3), "Non motorized", veh)]
+  demandpkm[, veh := ifelse(is.na(veh), vehicle_type, veh)]
+
   demandpkm = demandpkm[,.(demand_F = sum(demand_F)), by = c("iso", "year", "aggr_mode", "veh")]
   setnames(demandpkm, old = "veh", new = "vehicle_type")
 
-
-  demandpkm[, vehicle_type_plot := factor(vehicle_type, levels = c("LDV","Freight Rail", "Trucks (<3.5t)", "Trucks (3.5t-16)",  "Trucks (>16)", "Trucks",
+  demandpkm[, vehicle_type_plot := factor(vehicle_type, levels = c("LDV","Truck",
+                                                                "Freight Rail",
                                                                 "Motorbikes", "Small Cars", "Large Cars", "Van",
                                                                 "Domestic Aviation", "International Aviation","Bus", "Passenger Rail",
                                                                 "Freight", "Non motorized", "Shipping"))]
 
 
-  demandpkm[, mode := ifelse(vehicle_type %in% c("Freight", "Freight Rail", "Trucks", "Trucks (3.5t-16)",  "Trucks (>16)", "Shipping"),"freight", "pass")]
+  demandpkm[, mode := ifelse(vehicle_type %in% c("Freight", "Freight Rail", "Truck", "Shipping"),"freight", "pass")]
 
   demandpkm = merge(demandpkm, REMIND2ISO_MAPPING, by = "iso")
   demandpkm = demandpkm[, .(demand_F = sum(demand_F)), by = c("region", "year", "vehicle_type_plot", "aggr_mode", "mode")]
@@ -416,6 +442,7 @@ CO2km_intensity_newsalesplotf = function(annual_sales, mj_km_data, sharesVS1, sh
   totalemi = totalemi[,.(gCO2_km_ave = sum(gCO2_km_ave)), by = c("year", "iso", "vehicle_type")]
 
   totalemi = merge(totalemi, sharesVS1, all.x = TRUE, by = c("iso", "year", "vehicle_type"))
+  totalemi = totalemi[year %in% unique(sharesVS1$year)]
 
   totalemi = totalemi[,.(gCO2_km_ave = sum(gCO2_km_ave*share)), by = c("iso", "year", "subsector_L1")]
 
-- 
GitLab