Merge pull request #200 from MariannaR/inconv_comp

GDPcap-pkmcap plots for EDGE-T, new cache number

config/default.cfg

@@ -22,7 +22,7 @@ cfg$title <- "default"
 cfg$regionmapping <- "config/regionmappingH12.csv"
 
 #### Current input data revision (<mainrevision>.<subrevision>) ####
-cfg$revision <- 5.945
+cfg$revision <- 5.946
 
 #### Force the model to download new input data ####
 cfg$force_download <- FALSE

modules/30_biomass/magpie_40/bounds.gms

@@ -22,7 +22,7 @@ $ifthen.edge_esm_transport "%transport%" == "edge_esm"
 *** Slightly relaxed extraction bounds for biofuels.
 vm_fuExtr.up(t,regi,"pebios","5")$(t.val ge 2045)  = 1.4*p30_datapebio(regi,"pebios","5","maxprod",t);
 vm_fuExtr.up(t,regi,"pebios","5")$(t.val ge 2055)  = p30_datapebio(regi,"pebios","5","maxprod",t);
-vm_fuExtr.up(t,regi,"pebioil","5")$(t.val ge 2030) = 1.5*p30_datapebio(regi,"pebioil","5","maxprod",t);
+vm_fuExtr.up(t,regi,"pebioil","5")$(t.val ge 2030) = 2*p30_datapebio(regi,"pebioil","5","maxprod",t);
 vm_fuExtr.up(t,regi,"pebioil","5")$(t.val ge 2050) = p30_datapebio(regi,"pebioil","5","maxprod",t);
 $endif.edge_esm_transport
 

scripts/output/comparison/EDGEcomparison.R

@@ -48,6 +48,7 @@ EJfuelsFrgt_all = NULL
 emipSource_all = NULL
 costs_all = NULL
 pref_FV_all = NULL
+demgdpcap_all = NULL
 
 scenNames <- getScenNames(outputdirs)
 EDGEdata_path  <- path(outputdirs, paste("EDGE-T/"))
@@ -488,6 +489,25 @@ costscompFun = function(newcomp, sharesVS1,  EF_shares, pref_FV, capcost4Wall, c
 }
 
 
+demgdpcap_Fun = function(demkm, REMIND2ISO_MAPPING) {
+GDP_POP = getRMNDGDPcap()
+GDP_POP = merge(GDP_POP, REMIND2ISO_MAPPING, by = "iso")
+GDP_POP = merge(GDP_POP, REMIND2ISO_MAPPING, by = "iso")
+
+demcap_gdp = merge(demkm, GDP_POP, by = c("iso", "year"))
+demcap_gdp = merge(demcap_gdp, REMIND2ISO_MAPPING, by = "iso")
+demcap_gdp = demcap_gdp[,.(dem = sum(demand_F), gdp = sum(weight), pop = sum(POP_val)), by = .(region, sector, year)]
+demcap_gdp[, GDP_cap := gdp/pop]
+demcap_gdp[, demcap := dem*    ## in trillion km
+                  1e+6/   ## in million km
+                  pop]    ## in million km/million people=pkm/person
+
+return(demcap_gdp)
+}
+
+
+
+
 for (outputdir in outputdirs) {
   ## load mif file
   name_mif = list.files(path = outputdir, pattern = "REMIND_generic", full.names = F)
@@ -562,6 +582,8 @@ for (outputdir in outputdirs) {
   emipSource =  emipSourceFun(miffile)
   ## calculate costs by component
   costs = costscompFun(newcomp = newcomp, sharesVS1 = sharesVS1, EF_shares = EF_shares, pref_FV = pref_FV, capcost4Wall = capcost4Wall, capcost4W_BEVFCEV = capcost4W_BEVFCEV, nonf = nonf, totp = totp, REMIND2ISO_MAPPING)
+  ## per capita demand-gdp per capita
+  demgdpcap = demgdpcap_Fun(demkm = demandkm, REMIND2ISO_MAPPING)
   ## add scenario dimension to the results
   fleet[, scenario := as.character(unique(miffile$scenario))]
   salescomp[, scenario := unique(miffile$scenario)]
@@ -576,6 +598,7 @@ for (outputdir in outputdirs) {
   emipSource[, scenario := as.character(unique(miffile$scenario))]
   costs[, scenario := as.character(unique(miffile$scenario))]
   pref_FV[, scenario := as.character(unique(miffile$scenario))]
+  demgdpcap[,  scenario := as.character(unique(miffile$scenario))]
   ## rbind scenarios
   salescomp_all = rbind(salescomp_all, salescomp)
   fleet_all = rbind(fleet_all, fleet)
@@ -590,6 +613,7 @@ for (outputdir in outputdirs) {
   emipSource_all = rbind(emipSource_all, emipSource)
   costs_all = rbind(costs_all, costs)
   pref_FV_all = rbind(pref_FV_all, pref_FV)
+  demgdpcap_all = rbind(demgdpcap_all, demgdpcap)
 }
 
 ## create string with date and time
@@ -614,6 +638,7 @@ saveRDS(EJfuelsFrgt_all, paste0(outdir, "/EJfuelsFrgt_all.RDS"))
 saveRDS(emipSource_all, paste0(outdir, "/emipSource_all.RDS"))
 saveRDS(costs_all, paste0(outdir, "/costs_all.RDS"))
 saveRDS(pref_FV_all, paste0(outdir, "/pref_FV_all.RDS"))
+saveRDS(demgdpcap_all, paste0(outdir, "/demgdpcap_all.RDS"))
 file.copy(file.path("./scripts/output/comparison/notebook_templates", md_template), outdir)
 rmarkdown::render(path(outdir, md_template), output_format="pdf_document")
 

scripts/output/comparison/notebook_templates/EDGETransportComparison.Rmd

@@ -35,7 +35,7 @@ EJfrgt_all = readRDS("EJfuelsFrgt_all.RDS")
 emidem_all = readRDS("emidem_all.RDS")
 costs_all = readRDS("costs_all.RDS")
 pref_FV_all = readRDS("pref_FV_all.RDS")
-
+dempkm_cap_all = readRDS("demgdpcap_all.RDS")
 setConfig(forcecache=T)
 
 cols <- c("NG" = "#d11141",
@@ -463,7 +463,64 @@ emidem_pf = function(dt){
 emidem_pf(emidem_all)
 ```
 
+## demand per capita VS gdp per capita
+```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width=14, fig.height=12}
+demcapgdpcap_pf = function(dt) {
+
+dt = dt[year >= 2005]
+dt = dt[, year := as.character(year)]
+ppass_sm = ggplot()+
+  geom_line(data = dt[sector == "trn_pass"], aes(x = GDP_cap, y = demcap, color = region, group = region))+
+  geom_point(data = dt[sector == "trn_pass" & year %in% c("2020", "2030", "2050", "2070", "2100")], aes(x = GDP_cap, y = demcap, shape = year, group = region, color = region))+
+  theme_minimal()+
+  facet_grid(scenario~sector)+
+  theme(axis.text.x = element_text(angle = 90, hjust = 1),
+        axis.text = element_text(size=7),
+        title = element_text(size=8),
+        legend.text = element_text(size=8))+
+  labs(y = "Per capita demand [km/capita]", x="GDP per capita [$/person]", title = "Passenger short-medium")
+
+pfreigt_sm = ggplot()+
+  geom_line(data = dt[sector == "trn_freight"], aes(x = GDP_cap, y = demcap, color = region, group = region))+
+  geom_point(data = dt[sector == "trn_freight" & year %in% c("2020", "2030", "2050", "2070", "2100")], aes(x = GDP_cap, y = demcap, shape = year, group = region, color = region))+
+  theme_minimal()+
+  facet_grid(scenario~sector)+
+  theme(axis.text.x = element_text(angle = 90, hjust = 1),
+        axis.text = element_text(size=7),
+        title = element_text(size=8),
+        legend.text = element_text(size=8))+
+  labs(y = "Per capita demand [km/capita]", x="GDP per capita [$/person]", title = "Freight short-medium")
+
+
+ppass_lo = ggplot()+
+  geom_line(data = dt[sector == "trn_aviation_intl"], aes(x = GDP_cap, y = demcap, color = region, group = region))+
+  geom_point(data = dt[sector == "trn_aviation_intl" & year %in% c("2020", "2030", "2050", "2070", "2100")], aes(x = GDP_cap, y = demcap, shape = year, group = region, color = region))+
+  theme_minimal()+
+  facet_grid(scenario~sector)+
+  theme(axis.text.x = element_text(angle = 90, hjust = 1),
+        axis.text = element_text(size=7),
+        title = element_text(size=8),
+        legend.text = element_text(size=8))+
+  labs(y = "Per capita demand [km/capita]", x="GDP per capita [$/person]", title = "International aviation")
+
+pfreigt_lo = ggplot()+
+  geom_line(data = dt[sector == "trn_shipping_intl"], aes(x = GDP_cap, y = demcap, color = region, group = region))+
+  geom_point(data = dt[sector == "trn_shipping_intl" & year %in% c("2020", "2030", "2050", "2070", "2100")], aes(x = GDP_cap, y = demcap, shape = year, group = region, color = region))+
+  theme_minimal()+
+  facet_grid(scenario~sector)+
+  theme(axis.text.x = element_text(angle = 90, hjust = 1),
+        axis.text = element_text(size=7),
+        title = element_text(size=8),
+        legend.text = element_text(size=8))+
+  labs(y = "Per capita demand [km/capita]", x="GDP per capita [$/person]", title = "International shipping")
+
+p = list(pfreigt_lo = pfreigt_lo, ppass_lo = ppass_lo, pfreigt_sm = pfreigt_sm, ppass_sm = ppass_sm)
+
+return(p)
+}
 
+demcapgdpcap_pf(dempkm_cap_all)
+```
 
 ## Focus on slected region
 
@@ -544,9 +601,9 @@ costspf = function(dt, rp){
   dt_inc = dt[type == "inc"]
   dt_inc = dt_inc[,.(cost = sum(cost)), by = c("region", "year", "type", "technology", "scenario")]
   dt_tot = rbind(dt_inc[, c("name", "logit_type") := list("Inconvenience cost", "pinc")], dt[type == "real"]) 
-  dt[, technology := factor(technology, levels = legend_ord)]
-  dt_tot[, technology := factor(technology, levels = legend_ord)]
-
+  dt[, name := factor(name, levels = legend_ord)]
+  dt_tot[, name := factor(name, levels = legend_ord)]
+  dt_tot[, alph := ifelse(logit_type=="pinc", 0.6, 1)]
   plot1 = ggplot()+
          geom_bar(data = dt[type == "inc"], aes(x = year, y = cost, group = name, fill = name), position = "stack", stat = "identity")+
 theme_minimal()+
@@ -564,7 +621,7 @@ facet_grid(scenario~technology)+
     labs(y = "Costs [$/pkm]", x="")
 
   plot2 = ggplot()+
-         geom_bar(data = dt_tot, aes(x = year, y = cost, group = name, fill = name), position = "stack", stat = "identity")+
+         geom_bar(data = dt_tot, aes(x = year, y = cost, group = name, fill = name, alpha = alph), position = "stack", stat = "identity")+
 theme_minimal()+
 facet_grid(scenario~technology)+
     theme(axis.text.x = element_text(angle = 90, size=14, vjust=0.5, hjust=1),
@@ -575,9 +632,11 @@ facet_grid(scenario~technology)+
           legend.text = element_text(size=14),
           strip.text = element_text(size=12),
           strip.background = element_rect(color = "grey"))+
-    guides(fill=guide_legend(title="Cost component"))+
     scale_fill_manual(values = cols)+
-    labs(y = "Costs [$/pkm]", x="")
+    labs(y = "Costs [$/pkm]", x="")+
+      scale_alpha(range=c(0.4,1)) + 
+  guides(alpha=FALSE, linetype=FALSE,
+         fill=guide_legend(reverse=TRUE, title="Cost component"))
 
 plot = list(plot1 = plot1, plot2 = plot2)
 
@@ -782,12 +841,12 @@ ggsave("pTruck.png", pTruck, dpi=500, height = height , width = height * aspect_
 ```{r, echo=FALSE, warning=FALSE}
 
 prefBusTrucksplotf = function(pref){
-pref = pref[iso=="DEU" & technology %in% c("Electric", "FCEV") & vehicle_type %in% c("Bus_tmp_vehicletype")]
-pref = pref[year >= 2020 & year <= 2100 & vehicle_type == "Bus_tmp_vehicletype"]
-pref[, vehicle_type := "Buses and Trucks"]
+pref = pref[iso=="USA" & technology %in% c("Electric", "FCEV") & vehicle_type %in% c("Bus_tmp_vehicletype", "Truck (0-2.7t)")]
+pref = pref[year >= 2020 & year <= 2050]
+pref[, vehicle_type := ifelse(vehicle_type == "Bus_tmp_vehicletype", "Large Trucks and Buses", "Small Trucks")]
 
 p = ggplot()+
-  geom_line(data = pref, aes(x = year, y = value, group = technology, color = technology, linetype = technology))+
+  geom_line(data = pref, aes(x = year, y = value, group = interaction(technology, vehicle_type), color = technology, linetype = vehicle_type))+
   facet_wrap(~scenario, ncol = 1)+
   theme_minimal()+
   theme(axis.text.x = element_text(angle = 90, hjust = 1),
@@ -796,9 +855,9 @@ p = ggplot()+
         legend.text = element_text(size=7),
         strip.text = element_text(size=7))+
   labs(y = "Preference factor [-]", x="")+
-  scale_color_manual("Technology", values = cols,labels = c("Electric", "FCEV")) +
-  scale_linetype_manual("Technology", values = c(1,2),
-                        labels = c("Electric", "FCEV"))
+  scale_color_manual("Technology", values = cols,labels = c("Electric", "FCEV")) 
+  scale_linetype_manual("Vehicle", values = c(1,2),
+                        labels = c("Large Trucks and Buses", "Small Trucks"))
 return(p)
 }
 p = prefBusTrucksplotf(pref_FV_all)