Skip to content
Snippets Groups Projects
Commit 4415c63d authored by Marianna Rottoli's avatar Marianna Rottoli
Browse files

Reporting improvemente: dashboard new plots, loadFactor is loaded, emissions.

parent 029338a3
No related branches found
No related tags found
1 merge request!194Updates to EDGE-Transport
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 787 additions and 173 deletions
scripts/output/comparison/EDGEcomparison.R
+ 167
− 29
View file @ 4415c63d
......@@ -12,6 +12,8 @@ require(data.table)
require(rmndt)
require(moinput)
require(edgeTrpLib)
require(gdx)
require(gdxdt)
setConfig(forcecache = TRUE)
if(!exists("source_include")) {
......@@ -43,7 +45,8 @@ CO2km_int_newsales_all = NULL
emidem_all = NULL
EJfuelsPass_all = NULL
EJfuelsFrgt_all = NULL
emidemPass_all = NULL
emipdem_all = NULL
emipUp_all = NULL
scenNames <- getScenNames(outputdirs)
EDGEdata_path <- path(outputdirs, paste("EDGE-T/"))
......@@ -58,7 +61,7 @@ REMIND2ISO_MAPPING <- fread("config/regionmappingH12.csv")[, .(iso = CountryCode
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)
## I need the total demand for each region to get the average composition in the region (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]
......@@ -90,7 +93,7 @@ SalesFun = function(shares_LDV, newcomp, sharesVS1){
}
fleetFun = function(vintcomp, newcomp, sharesVS1){
fleetFun = function(vintcomp, newcomp, sharesVS1, loadFactor){
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)]
......@@ -102,9 +105,9 @@ fleetFun = function(vintcomp, newcomp, sharesVS1){
"year"), measure.vars = c("vintdem", "newdem"))
allfleet[,alpha:=ifelse(variable == "vintdem", 0, 1)]
load_factor = 1.5
allfleet = merge(allfleet, loadFactor, all.x = TRUE, by = c("iso", "vehicle_type", "year"))
annual_mileage = 15000
allfleet = allfleet[,.(value = sum(value/load_factor/annual_mileage)), by = c("iso", "technology", "variable", "year")]
allfleet = allfleet[,.(value = sum(value/loadFactor/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")]
......@@ -165,8 +168,8 @@ EJmodeFun = function(demandEJ){
ESmodeFun = function(demandkm, POP){
## REMIND-EDGE results
demandkm <- demandkm[,c("sector","subsector_L3","subsector_L2",
"subsector_L1","vehicle_type","technology", "iso","year","demand_F")]
"subsector_L1","vehicle_type","technology", "iso","year","demand_F")]
## attribute aggregated mode and vehicle names for plotting purposes, and aggregate
demandkm[, aggr_mode := ifelse(subsector_L1 %in% c("Three-Wheeler", "trn_pass_road_LDV_4W"), "LDV", NA)]
demandkm[, aggr_mode := ifelse(sector %in% c("trn_freight", "trn_shipping_intl"), "Freight", aggr_mode)]
......@@ -191,9 +194,9 @@ ESmodeFun = function(demandkm, POP){
demandkm[, vehicle_type_plot := factor(vehicle_type, levels = c("LDV","Freight Rail", "Truck", "Domestic Ship", "International Ship",
"Motorbikes", "Small Cars", "Large Cars", "Van",
"Domestic Aviation", "International Aviation","Bus", "Passenger Rail",
"Freight", "Non motorized", "Shipping"))]
"Motorbikes", "Small Cars", "Large Cars", "Van",
"Domestic Aviation", "International Aviation","Bus", "Passenger Rail",
"Freight", "Non motorized", "Shipping"))]
## attribute aggregate mode (passenger, freight)
demandkm[, mode := ifelse(vehicle_type %in% c("Freight", "Freight Rail", "Truck", "Shipping") ,"freight", "pass")]
......@@ -224,7 +227,7 @@ ESmodeFun = function(demandkm, POP){
}
FEliq_sourceFun = function(FEliq_source, gdp){
## Attribute oil and biodiesel (TODO Coal2Liquids is accounted for as Oil!
## Attribute oil and biodiesel (TODO Coal2Liquids is accounted for as Oil!
FEliq_source[, technology := ifelse(variable %in% c("FE|Transport|Liquids|Oil", "FE|Transport|Liquids|Coal"), "Oil", NA)]
FEliq_source[, technology := ifelse(variable %in% c("FE|Transport|Liquids|Biomass"), "Biodiesel", technology)]
FEliq_source[, technology := ifelse(variable %in% c("FE|Transport|Liquids|Hydrogen"), "Synfuel", technology)]
......@@ -233,15 +236,16 @@ FEliq_sourceFun = function(FEliq_source, gdp){
FEliq_sourceR = FEliq_source[][, shareliq := value/sum(value),by=c("region", "year")]
## to ISO level
FEliq_sourceISO <- disaggregate_dt(FEliq_source, REMIND2ISO_MAPPING,
valuecol="value",
datacols=c("model","scenario", "unit","technology"),
weights=gdp)
valuecol="value",
datacols=c("model","scenario", "unit","technology"),
weights=gdp)
## calculate share
FEliq_sourceISO[, shareliq := value/sum(value),by=c("iso", "year")]
return(list(FEliq_sourceISO = FEliq_sourceISO, FEliq_sourceR = FEliq_sourceR))
}
CO2km_int_newsales_Fun = function(shares_LDV, mj_km_data, sharesVS1, FEliq_source, gdp){
## energy intensity https://en.wikipedia.org/wiki/Energy_density
# emi_petrol = 45 ## MJ/gFUEL
......@@ -331,19 +335,139 @@ emidemFun = function(emidem){
emidemPassFun = function(emidemPass){
emidemPass = emidemPass[region!="World" & year >= 2015 & year <= 2100]
emidemPass[, variable := as.character(variable)]
emidemPass[, c("model", "scenario", "variable", "unit") := NULL]
setnames(emidemPass, old = "value", new = "emi_sum")
emidemPass[, type := "Demand"]
return(emidemPass)
}
emipUpFun = function(prodSe, prodFe, vemi, pemi, pass_sm_el_emi){
prodSe <- readgdx(gdx, "vm_prodSe")[
tall >= minyr & tall <= maxyr]
prodFe <- readgdx(gdx, "vm_prodFe")[
ttot >= minyr & ttot <= maxyr]
prodFefos <- prodFe[all_te %in% c("tdfosdie", "tdfospet")]
TWa_2_EJ <- 31.536
GtC_2_MtCO2 <- 44 / 12 * 1000
prodSe[, value := value*TWa_2_EJ]
prodFe[, value := value*TWa_2_EJ]
setnames(prodFe, c("year", "region", "se", "fe", "te", "value"))
setnames(prodFefos, c("year", "region", "se", "fe", "te", "value"))
setnames(prodSe, c("year", "region", "pe", "se", "te", "value"))
## Fe
prodFe <- prodFe[fe %in% c("fepet", "fedie")]
prodFe[, all_liq := sum(value), by=.(year, region)]
prodFe[, liq_val := sum(value), by=.(year, region, se)]
prodFe[, c("fe", "te", "value") := NULL]
prodFe <- unique(prodFe)
prodSyn <- prodSe[te == "MeOH"]
prodSyn[, se := "synliq"]
prodSyn[, c("pe", "te") := NULL]
setnames(prodSyn, "value", "syn_val")
prodSyn <- rbind(prodSyn,
CJ(year=c(2020), region = unique(prodFe$region), se="synliq", syn_val=0))
allLiq <- merge(prodSyn, prodFe, all=T)
allLiq[is.na(syn_val), syn_val := 0]
allLiq[, all_liq := max(all_liq, na.rm=T), by=.(year, region)]
allLiq[, syn_val := max(syn_val, na.rm=T), by=.(year, region)]
## separate synthetic fuels
allLiq[se == "synliq", liq_val := syn_val]
allLiq[se == "seliqfos", liq_val := liq_val - syn_val]
allLiq <- allLiq[
data.table(se=c("seliqbio", "seliqfos", "synliq"),
lname=c("biofuels", "fossil fuels", "synfuels")),
on="se"]
## Electricity production - not really important but interesting
prodEl <- prodSe[se == "seel"]
## Electricity consumption for H2 production
prodEl[, allEl := sum(value), by=.(year, region)][
, c("pe", "se", "te", "value") := NULL
]
prodEl <- unique(prodEl)
prodElH2 <- prodSe[pe == "seel" & se == "seh2"]
prodElH2[, elh2 := sum(value), by=.(year, region)]
prodElH2[, c("pe", "se", "te", "value") := NULL]
prodElH2 <- unique(prodElH2)
prodShare <- merge(prodEl, prodElH2, all=T)
prodShare[, netProd := allEl - elh2]
prodVert <- melt(prodShare, id.vars = c("year", "region"))[
data.table(variable=c("elh2", "netProd"),
label=c("Electricity for H2", "Other Demands")),
on="variable"
]
setnames(vemi, c("year", "region", "pe", "se", "te", "emi", "val"))
emico2 <- vemi[se == "seel" & emi == "co2"][
, co2val := sum(val) * GtC_2_MtCO2, by=.(year, region)][
, c("pe", "se", "te", "emi", "val") := NULL
]
emico2 <- unique(emico2)
emico2 <- merge(emico2, prodShare, by=c("year", "region"))
emico2[, int := co2val/allEl] # MtCO2/EJ -> tCO2/MJ
## emissions from electricity for hydrogen production
emisyn <- emico2[, .(se="seh2", co2 = int * elh2), by=.(year, region)] # MtCO2
## synfuel share
syn_fos <- allLiq[se != "seliqbio"][
, .(syn_share = syn_val/sum(liq_val)), by=.(year, region)]
syn_fos <- unique(syn_fos)
## final energy from prodFefos
setnames(prodFefos, c("year", "region", "se", "fe", "te", "prodFE"))
prodFefos <- prodFefos[syn_fos, on=.(year, region)]
## substract synfuels
prodFefos[, no_syn := prodFE * (1-syn_share)]
setnames(pemi, c("fe", "emi_fac"))
prodFefos <- pemi[prodFefos, on="fe"]
prodFefos[, co2 := emi_fac * no_syn * TWa_2_EJ]
## add synthetic fuel/ upstream emissions
emi_all <- merge(emisyn, prodFefos[, .(year, region, se=fe, co2)], by=c("year", "region", "se", "co2"), all=T)
## add emissions from electricity
pass_sm_el_emi[, c("model", "scenario", "variable", "unit") := NULL]
pass_sm_el_emi[, se := "electricity"]
pass_sm_el_emi = pass_sm_el_emi[year %in% unique(emi_all$year)]
setnames(pass_sm_el_emi, c("region", "co2", "year", "se"))
emi_all = rbind(emi_all, pass_sm_el_emi)
emi_all = emi_all[, .(emi_upstr = sum(co2)), by=.(year, region)]
emi_all[, type := "Upstream"]
return(emi_all)
}
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)))
miffile[, region:=as.character(region)]
miffile[, region := as.character(region)]
miffile[, year := period]
miffile[, period:=NULL]
miffile[, period := NULL]
miffile = miffile[region != "World"]
## load gdx file
gdx = paste0(outputdir, "/fulldata.gdx")
## load RDS files
sharesVS1 = readRDS(paste0(outputdir, "/EDGE-T/", "shares.RDS"))[["VS1_shares"]]
newcomp = readRDS(paste0(outputdir, "/EDGE-T/", "newcomp.RDS"))
......@@ -352,25 +476,34 @@ for (outputdir in outputdirs) {
demandEJ = readRDS(paste0(outputdir, "/EDGE-T/", "demandF_plot_EJ.RDS"))
demandkm = readRDS(paste0(outputdir, "/EDGE-T/", "demandF_plot_pkm.RDS"))
mj_km_data = readRDS(paste0(outputdir, "/EDGE-T/", "mj_km_data.RDS"))
loadFactor = readRDS(paste0(outputdir, "/EDGE-T/", "loadFactor.RDS"))
## load data from gdx file
minyr <- 2020
maxyr <- 2100
prodSe <- readgdx(gdx, "vm_prodSe")[tall >= minyr & tall <= maxyr]
prodFe <- readgdx(gdx, "vm_prodFe")[ttot >= minyr & ttot <= maxyr]
vemi <- readgdx(gdx, "vm_emiTeDetail")[ttot >= minyr & ttot <= maxyr]
pemi <- readgdx(gdx, "p_ef_dem") ## emission factors
## load population and GDP
POP_country=calcOutput("Population", aggregate = F)[,, "pop_SSP2"]
POP <- magpie2dt(POP_country, regioncol = "iso",
yearcol = "year", datacols = "POP")
yearcol = "year", datacols = "POP")
gdp <- getRMNDGDP(scenario = "gdp_SSP2", usecache = T)
## select useful entries from mif file
FEliq_source = miffile[variable %in% c("FE|Transport|Liquids|Biomass", "FE|Transport|Liquids|Hydrogen", "FE|Transport|Liquids|Coal", "FE|Transport|Liquids|Oil"),]
emidem = miffile[variable %in% c("Emi|CO2|Transport|Demand"),]
emidemPass = miffile[variable %in% c("Emi|CO2|Transport|Pass|Short-Medium Distance|Demand"),]
pass_sm_el_emi = miffile[variable %in% c("Emi|CO2|Transport|Pass|Short-Medium Distance|Electricity")]
## modify mif file entries to be used in the functions
FEliq_source = FEliq_sourceFun(FEliq_source, gdp)
## calculate sales
salescomp = SalesFun(shares_LDV, newcomp, sharesVS1)
## calculate fleet compositons
fleet = fleetFun(vintcomp, newcomp, sharesVS1)
fleet = fleetFun(vintcomp, newcomp, sharesVS1, loadFactor)
## calculate EJ from LDVs by technology
EJroad = EJroadFun(demandEJ)
## calculate FE demand by mode
......@@ -388,7 +521,9 @@ for (outputdir in outputdirs) {
## calculate demand emissions
emidem = emidemFun(emidem)
## calculate emissions from passenger demand SM
emidemPass = emidemPassFun(emidemPass)
emipdem = emidemPassFun(emidemPass)
## calculate upstream emissions (accounting for synfuels, hydrogen and electricity production)
emipUp = emipUpFun(prodSe, prodFe, vemi, pemi, pass_sm_el_emi)
## add scenario dimension to the results
fleet[, scenario := as.character(unique(miffile$scenario))]
salescomp[, scenario := unique(miffile$scenario)]
......@@ -400,7 +535,8 @@ for (outputdir in outputdirs) {
emidem[, scenario := as.character(unique(miffile$scenario))]
EJfuelsPass[, scenario := as.character(unique(miffile$scenario))]
EJfuelsFrgt[, scenario := as.character(unique(miffile$scenario))]
emidemPass[, scenario := as.character(unique(miffile$scenario))]
emipdem[, scenario := as.character(unique(miffile$scenario))]
emipUp[, scenario := as.character(unique(miffile$scenario))]
## rbind scenarios
salescomp_all = rbind(salescomp_all, salescomp)
fleet_all = rbind(fleet_all, fleet)
......@@ -412,7 +548,8 @@ for (outputdir in outputdirs) {
emidem_all = rbind(emidem_all, emidem)
EJfuelsPass_all = rbind(EJfuelsPass_all, EJfuelsPass)
EJfuelsFrgt_all = rbind(EJfuelsFrgt_all, EJfuelsFrgt)
emidemPass_all = rbind(emidemPass_all, emidemPass)
emipdem_all = rbind(emipdem_all, emipdem)
emipUp_all = rbind(emipUp_all, emipUp)
}
outdir = paste0("output/comparerunEDGE", gsub(" | ^([[:alpha:]]*).*","", Sys.time()))
......@@ -430,7 +567,8 @@ saveRDS(CO2km_int_newsales_all, paste0(outdir, "/CO2km_int_newsales_all.RDS"))
saveRDS(emidem_all, paste0(outdir, "/emidem_all.RDS"))
saveRDS(EJfuelsPass_all, paste0(outdir, "/EJfuelsPass_all.RDS"))
saveRDS(EJfuelsFrgt_all, paste0(outdir, "/EJfuelsFrgt_all.RDS"))
saveRDS(emidemPass_all, paste0(outdir, "/emidemPass_all.RDS"))
saveRDS(emipdem_all, paste0(outdir, "/emipdem_all.RDS"))
saveRDS(emipUp_all, paste0(outdir, "/emipUp_all.RDS"))
file.copy(file.path("./scripts/output/comparison/notebook_templates", md_template), outdir)
rmarkdown::render(path(outdir, md_template), output_format="pdf_document")
......@@ -441,8 +579,8 @@ if (length(outputdirs) == 5 &
isTRUE(any(grepl("Budg1100_ElecEra", outputdirs))) &
isTRUE(any(grepl("Budg1100_HydrHype", outputdirs))) &
isTRUE(any(grepl("Base_ConvCase", outputdirs)))){
file.copy(file.path("./scripts/output/comparison/notebook_templates/helper_dashboard.R"), outdir)
file.copy(file.path("./scripts/output/comparison/notebook_templates", dash_template), outdir)
rmarkdown::render(path(outdir, dash_template))
file.copy(file.path("./scripts/output/comparison/notebook_templates/helper_dashboard.R"), outdir)
file.copy(file.path("./scripts/output/comparison/notebook_templates", dash_template), outdir)
rmarkdown::render(path(outdir, dash_template))
}
scripts/output/comparison/notebook_templates/EDGETransportComparison.Rmd
+ 230
− 5
View file @ 4415c63d
......@@ -117,7 +117,7 @@ legend_ord = c(legend_ord_modes, legend_ord_fuels, legend_ord_costs)
regionplot = "EUR"
```
```{r, echo=FALSE, message=FALSE, warning=FALSE}
......@@ -155,11 +155,10 @@ vintcomparisonpf = function(dt){
}
vintcomparisonpf(fleet_all)
## Sales composition
```
## Sales composition
```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width=14, fig.height=12}
salescompf = function(dt){
......@@ -230,7 +229,7 @@ EJroadpf = function(dt){
strip.background = element_rect(color = "grey"))
plotTrucks = ggplot()+
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()+
......@@ -248,7 +247,7 @@ EJroadpf = function(dt){
strip.text = element_text(size = 14),
strip.background = element_rect(color = "grey"))
return(plotlist = list(plotLDV = plotLDV, plotBus = plotBus, plotTrucks = plotTrucks))
return(plotlist = list(plotLDV = plotLDV, plotBus = plotBus, plotTruck = plotTruck))
}
EJroadpf(EJroad_all)
......@@ -452,3 +451,229 @@ emidem_pf = function(dt){
emidem_pf(emidem_all)
```
## Focus on slected region
## vintages
```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width=14, fig.height=12}
vintcomparison_regi_pf = function(dt, rp){
dt = dt[year %in% c(2015, 2030, 2050, 2100) & region == rp]
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_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"))+
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)
}
p = vintcomparison_regi_pf(fleet_all, rp = regionplot)
p
aspect_ratio <- 1.5
height <- 6
ggsave("pvint.png", p, dpi=500, height = height , width = height * aspect_ratio)
```
## Sales composition
```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width=14, fig.height=12}
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)
```
```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width=14, fig.height=12}
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)
```
```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width=14, fig.height=12}
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)
```
scripts/output/single/notebook_templates/EDGETransportReport.Rmd
+ 5
− 5
View file @ 4415c63d
......@@ -123,7 +123,7 @@ 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"))
loadFactor <- readRDS(datapath(fname = "loadFactor.RDS"))
## Load population to calculate per capita values
POP_country=calcOutput("Population", aggregate = F)[,, "pop_SSP2"]
POP <- magpie2dt(POP_country, regioncol = "iso",
......@@ -138,7 +138,7 @@ miffile <- as.data.table(read.quitte(name_mif))
```{r, echo=FALSE, warning=FALSE}
plotVint = function(vintcomp, newcomp, sharesVS1){
plotVint = function(vintcomp, newcomp, sharesVS1, loadFactor){
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)]
......@@ -150,9 +150,9 @@ plotVint = function(vintcomp, newcomp, sharesVS1){
"year"), measure.vars = c("vintdem", "newdem"))
allfleet[,alpha:=ifelse(variable == "vintdem", 0, 1)]
load_factor = 2
allfleet = merge(allfleet, loadFactor, all.x = TRUE, by = c("iso", "vehicle_type", "year"))
annual_mileage = 15000
allfleet = allfleet[,.(value = sum(value/load_factor/annual_mileage)), by = c("iso", "technology", "variable", "year")]
allfleet = allfleet[,.(value = sum(value/loadFactor/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")]
......@@ -184,7 +184,7 @@ plotVint = function(vintcomp, newcomp, sharesVS1){
}
p = plotVint(vintcomp, newcomp, shares$VS1_shares)
p = plotVint(vintcomp, newcomp, shares$VS1_shares, loadFactor)
p
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment