diff --git a/README.md b/README.md index 54db8cd5e4d2385fc9b2004c87e9fbbe0cf3f133..87a46ece1de31ecd7899751504852c57321b31cd 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@ The model documentation for version 2.1 can be found at XXX. A most recent version of the documentation can also be extracted from the model source code via the R package goxygen (https://github.com/pik-piam/goxygen). To extract the documentation, install the -package and run the main function `(goxygen(unitPattern = c("\\[","\\]"), includeCore=T))` +package and run the main function `(goxygen(unitPattern = c("\\[","\\]"), includeCore=T, use_advanced_interfacePlot_function=T))` in the main folder of the model. The resulting documentation can be found in the folder "doc". diff --git a/config/default.cfg b/config/default.cfg index a579bd5cda0aa9842a1aa7d9fbd17bc0f27ff99f..65c57a9f0cabe997d3b03ca82c1c310ab9a51c19 100755 --- a/config/default.cfg +++ b/config/default.cfg @@ -22,7 +22,7 @@ cfg$title <- "default" cfg$regionmapping <- "config/regionmappingH12.csv" #### Current input data revision (<mainrevision>.<subrevision>) #### -cfg$revision <- 5.938 +cfg$revision <- 5.9385 #### Force the model to download new input data #### cfg$force_download <- FALSE diff --git a/core/sets.gms b/core/sets.gms index 84e1cbcc6e73cdb4ad2a48bd714ead8120f78aab..732842a198b5f84198ab66cec473fff0a53f2887 100755 --- a/core/sets.gms +++ b/core/sets.gms @@ -782,6 +782,7 @@ $ENDIF.RegScenCapt ***######################## R SECTION START (MODULES) ############################### *** THIS CODE IS CREATED AUTOMATICALLY, DO NOT MODIFY THESE LINES DIRECTLY *** ANY DIRECT MODIFICATION WILL BE LOST AFTER NEXT MODEL START +*** CHANGES CAN BE DONE USING THE RESPECTIVE LINES IN scripts/start_functions.R sets @@ -822,7 +823,7 @@ sets codePerformance / - module2realisation(modules,*) "mapping of modules and active realisations" / +module2realisation(modules,*) "mapping of modules and active realisations" / macro . %macro% welfare . %welfare% PE_FE_parameters . %PE_FE_parameters% diff --git a/main.gms b/main.gms index fe30cace39e2c03b81c0eefa41bf7d54e9fbe5cd..5f3f6647b0343b43f2e9e3a380637a3c49e3b315 100644 --- a/main.gms +++ b/main.gms @@ -7,64 +7,74 @@ *** SOF ./main.gms *' @title REMIND - REgional Model of INvestments and Development *' -*' @description REMIND is a global multi-regional model incorporating the economy, the climate system and a detailed representation of the energy sector. -*' It solves for an inter-temporal Pareto optimum in economic and energy investments in the model regions, fully accounting for interregional trade in goods, -*' energy carriers and emissions allowances. REMIND allows for the analysis of technology options and policy proposals for climate mitigation. +*' @description REMIND is a global multi-regional model incorporating the economy, the climate system +*' and a detailed representation of the energy sector. It solves for an intertemporal Pareto optimum +*' in economic and energy investments in the model regions, fully accounting for interregional trade in +*' goods, energy carriers and emissions allowances. REMIND enables analyses of technology options and +*' policy proposals for climate change mitigation. *' -*' The macro-economic core of REMIND is a Ramsey-type optimal growth model in which intertemporal global welfare is optimized subject to equilibrium constraints ([02_welfare]). -*' Intertemporal optimization ([80_optimization]) with perfect foresight is subject to market clearing. The model explicitly represents trade in final goods, primary energy carriers, -*' and in the case of climate policy, emissions allowances. Macro-economic production factors are capital, labor, and final energy. -*' The production function with constant elasticity of substitution (nested CES production function) determines the final energy demand ([29_CES_parameters]). -*' REMIND uses economic output for investments in the macro-economic capital stock as well as consumption, trade, and energy system expenditures. +*' The macro-economic core of REMIND is a Ramsey-type optimal growth model in which intertemporal global +*' welfare is optimized subject to equilibrium constraints ([02_welfare]). Intertemporal optimization +*' ([80_optimization]) with perfect foresight is subject to market clearing. The model explicitly represents +*' trade in final goods, primary energy carriers, and when certain climate policies are enabled, emissions +*' allowances ([24_trade]). The macro-economic production factors are capital, labor, and final energy. +*' A nested production function with constant elasticity of substitution determines the final energy demand +*' ([01_macro], [29_CES_parameters]). REMIND uses economic output for investments in the macro-economic +*' capital stock as well as for consumption, trade, and energy system expenditures. *' -*' The macro-economic core and the energy system part are hard-linked via the final energy demand and costs incurred by the energy system. -*' Economic activity results in demand for final energy in different sectors such as transport energy ([35_transport]), electricity ([32_power]), -*' and non-electric energy for stationary end uses ([38_stationary]) splitted into an industy ([37_industry]) and buildings ([36_buildings]) sector. - -*' The primary energy carriers in REMIND include both exhaustible and renewable resources. Exhaustible resources comprise uranium as well as three fossil resources ([31_fossil]), namely coal, oil, and gas. -*' Renewable resources include hydro, wind, solar, geothermal, and biomass ([30_biomass]). -*' More than 50 technologies are available for the conversion of primary energy into secondary energy carriers as well as for the distribution of secondary energy carriers into final energy. +*' The macro-economic core and the energy system part are hard-linked via the final energy demand and the +*' costs incurred by the energy system. Economic activity results in demand for final energy in different +*' sectors (transport ([35_transport]), industry ([37_industry]), buildings ([36_buildings])...) and of +*' different type (electric ([32_power]) and non-electric). The primary energy carriers in REMIND include +*' both exhaustible and renewable resources. Exhaustible resources comprise uranium as well as three fossil +*' resources ([31_fossil]), namely coal, oil, and gas. Renewable resources include hydro, wind, solar, +*' geothermal, and biomass ([30_biomass]). More than 50 technologies are available for the conversion of +*' primary energy into secondary energy carriers as well as for the distribution of secondary energy carriers +*' into final energy. *' -*' The model accounts for the full range of anthropogenic greenhouse gas (GHG) emissions, most of which are represented by source. -*' REMIND simulates emissions from long-lived GHGs (CO2, CH4, N2O), short-lived GHGs (CO, NOx, VOC) and aerosols (SO2, BC, OC). -*' It accounts for these emissions with different levels of detail depending on the types and sources of emissions. -*' It calculates CO2 emissions from fuel combustion, CH4 emissions from fossil fuel extraction and residential energy use and N2O emissions from energy supply based on sources. +*' The model accounts for the full range of anthropogenic greenhouse gas (GHG) emissions, most of which are +*' represented by source. REMIND simulates emissions from long-lived GHGs (CO2, CH4, N2O), short-lived GHGs +*' (CO, NOx, VOC) and aerosols (SO2, BC, OC). It accounts for these emissions with different levels of detail +*' depending on the types and sources of emissions. It calculates CO2 emissions from fuel combustion, CH4 +*' emissions from fossil fuel extraction and residential energy use, and N2O emissions from energy supply +*' based on sources. *' -*' The numerical code is structured in a modular way. The technical structure looky as follows: At the top level you find the folders config, core, modules and scripts. -*' The overall structure is build in the file main.gms. All settings and configuration information is given in the config folder. -*' The core folder contains all files that are part of the core of the REMIND model. For each module there exists a sub-folder in the modules folder. -*' Helpful scripts for e.g. starting a run or analysing results you find in the scripts folder. +*' The code is structured in a modular way, with code belonging either to the model's core, or to one of the +*' modules. The folder structure is as follows: at the top level are the folders config, core, modules and +*' scripts. The config folder contains the REMIND settings and configuration information. The core folder +*' contains all the files that are part of the core. The modules folder holds all the files that belong to +*' the modules, with numbered sub-folders for every module. The scripts folder contains helpful scripts for +*' starting a model run and analysing results. *' -*' In the main.gms the technical structure of REMIND can be found. -*' First, the *.gms files from the core folder are included and afterward the *.gms files from the activated module realization, -*' beginning with the one with the smallest module-number. The technical structure of REMIND looks as follows: +*' REMIND is run by executing the main.gms file, which loads the configuration information and builds the model, +*' by concatenating all necessary files from the core and modules folders into a single file called full.gms. +*' The concatenation process starts with files from the core and continues with files from activated modules, +*' in increasing order of module-number. It observes the following structure: *' *' { width=100% } *' -*' In general, the .gms-files in each module realization can be the same as in the core. -*' For each module it has to be clearly defined what kind of interfaces it has with the core part of the model. *' -*' The REMIND GAMS code folllows a coding etiquette including the following prefixes: +*' The GAMS code follows a naming etiquette based on the following prefixes: *' -*' * q_ eQuations -*' * v_ Variables -*' * s_ Scalars -*' * f_ File parameters - these parameters contain data as it was read from file -*' * o_ Output parameters - only being influenced by optimization but without effect on the optimization -*' * c_ switches from the Config.gms - parameters, that are switches to choose different scenarios +*' * "q_" to designate equations, +*' * "v_" to designate variables, +*' * "s_" to designate scalars, +*' * "f_" to designate file parameters (parameters that contain unaltered data read in from input files), +*' * "o_" to designate output parameters (parameters that do not affect the optimization, but are affected by it), +*' * "c_" to designate switches (parameters that enable different configuration choices), +*' * "s_FIRSTUNIT_2_SECONDUNIT" to designate a scalar used to convert from the FIRSTUNIT to the SECONDUNIT +*' through multiplication, e.g. s_GWh_2_EJ. *' -*' The prefixes are extended in some cases by a second letter: +*' These prefixes are extended in some cases by a second letter: *' -*' * ?m_ module-relevant object - This object is used by at least one module and the core code. Changes related to this object have to be performed carefully. -*' * ?00_ (a 2-digit number) module-only object - This 2-digit number defines the module the object belongs to. The number is used here to make sure that different modules cannot have the same object. +*' * "?m_" to designate objects used in the core and in at least one module. +*' * "?00_" to designate objects used in a single module, exclusively, with the 2-digit number corresponding +*' to the module number. *' -*' Sets are treated slightly different: Instead of adding a prefix sets should get a 2-digit number suffix giving the number of the module in which the set is exclusively used. -*' If the set is used in more than one module no suffix should be given. For specific sets also prefixes exist: -*' -*' * s_FIRSTUNIT_2_SECONDUNIT unit conversion scalar - a scalar that is used to convert from FIRSTUNIT to SECONDUNIT by multiplying - example: s_GWh_2_EJ. -*' * c_@ - configuration switch, must be defined and assigned in the config.gms file. It's practically the former switches we had (emiscen, climscen, etc). +*' Sets are treated differently: instead of a prefix, sets exclusively used within a module get that module's +*' number added as a suffix. If the set is used in more than one module no suffix is given. *' -*' The units (e.g., TWa, EJ, GtC, GtCO2, ...) of varialbles and parameters are documented at the location of the variable and parameter declaration in \[ \]. +*' The units (e.g., TWa, EJ, GtC, GtCO2, ...) of variables and parameters are documented in the declaration files. @@ -72,9 +82,9 @@ * * Regionscode: 690d3718e151be1b450b394c1064b1c5 * -* Input data revision: 5.938 +* Input data revision: 5.9385 * -* Last modification (input data): Fri Feb 14 10:03:25 2020 +* Last modification (input data): Mon Feb 24 15:24:12 2020 * *###################### R SECTION END (VERSION INFO) ########################### diff --git a/modules/24_trade/standard/bounds.gms b/modules/24_trade/standard/bounds.gms index 2f31953cac48ccd2a51712643253d33cd8048fa5..91f04b1313f9da1152e662de3b801282563ecd13 100644 --- a/modules/24_trade/standard/bounds.gms +++ b/modules/24_trade/standard/bounds.gms @@ -92,6 +92,30 @@ vm_Xport.up("2010",regi,"peoil") = 1.05 * pm_IO_trade("2010",regi,"peoil","Xport vm_Xport.lo("2015",regi,"peoil") = 0.95 * pm_IO_trade("2015",regi,"peoil","Xport"); vm_Xport.up("2015",regi,"peoil") = 1.05 * pm_IO_trade("2015",regi,"peoil","Xport"); +*** bounds on gas exports in 2010 and 2015 +vm_Xport.lo("2010",regi,"pegas") = 0.95 * pm_IO_trade("2010",regi,"pegas","Xport"); +vm_Xport.up("2010",regi,"pegas") = 1.05 * pm_IO_trade("2010",regi,"pegas","Xport"); +vm_Xport.lo("2015",regi,"pegas") = 0.95 * pm_IO_trade("2015",regi,"pegas","Xport"); +vm_Xport.up("2015",regi,"pegas") = 1.05 * pm_IO_trade("2015",regi,"pegas","Xport"); + +*** bounds on gas imports in 2010 and 2015 +vm_Mport.lo("2010",regi,"pegas") = 0.95 * pm_IO_trade("2010",regi,"pegas","Mport"); +vm_Mport.up("2010",regi,"pegas") = 1.05 * pm_IO_trade("2010",regi,"pegas","Mport"); +vm_Mport.lo("2015",regi,"pegas") = 0.95 * pm_IO_trade("2015",regi,"pegas","Mport"); +vm_Mport.up("2015",regi,"pegas") = 1.05 * pm_IO_trade("2015",regi,"pegas","Mport"); + +*** bounds on coal exports in 2010 and 2015 +vm_Xport.lo("2010",regi,"pecoal") = 0.95 * pm_IO_trade("2010",regi,"pecoal","Xport"); +vm_Xport.up("2010",regi,"pecoal") = 1.05 * pm_IO_trade("2010",regi,"pecoal","Xport"); +vm_Xport.lo("2015",regi,"pecoal") = 0.95 * pm_IO_trade("2015",regi,"pecoal","Xport"); +vm_Xport.up("2015",regi,"pecoal") = 1.05 * pm_IO_trade("2015",regi,"pecoal","Xport"); + +*** bounds on coal imports in 2010 and 2015 +vm_Mport.lo("2010",regi,"pecoal") = 0.95 * pm_IO_trade("2010",regi,"pecoal","Mport"); +vm_Mport.up("2010",regi,"pecoal") = 1.05 * pm_IO_trade("2010",regi,"pecoal","Mport"); +vm_Mport.lo("2015",regi,"pecoal") = 0.95 * pm_IO_trade("2015",regi,"pecoal","Mport"); +vm_Mport.up("2015",regi,"pecoal") = 1.05 * pm_IO_trade("2015",regi,"pecoal","Mport"); + *** upper bounds ( 1% yearly growth rate) on all big oil exporters (more than 15EJ in 2010) in 2020, 2025 and 2030 loop(regi, if( (pm_IO_trade("2010",regi,"peoil","Xport") ge (15*sm_EJ_2_TWa)), @@ -102,4 +126,4 @@ loop(regi, ); -*** EOF ./modules/24_trade/standard/bounds.gms \ No newline at end of file +*** EOF ./modules/24_trade/standard/bounds.gms diff --git a/modules/35_transport/edge_esm/output.gms b/modules/35_transport/edge_esm/output.gms new file mode 100644 index 0000000000000000000000000000000000000000..49c34856e554c3fdaabef32622f91e31189bb5c4 --- /dev/null +++ b/modules/35_transport/edge_esm/output.gms @@ -0,0 +1,10 @@ +*** | (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 +*** SOF ./modules/35_transport/edge_esm/output.gms + +Execute "Rscript EDGE_transport.R --reporting"; +*** EOF ./modules/35_transport/edge_esm/output.gms diff --git a/modules/35_transport/edge_esm/realization.gms b/modules/35_transport/edge_esm/realization.gms index 90ede08e76b4840a88170c17195d245e7210f4ad..6ada3a91d37df1332df5cb49b29919f86613e43a 100644 --- a/modules/35_transport/edge_esm/realization.gms +++ b/modules/35_transport/edge_esm/realization.gms @@ -25,5 +25,6 @@ $Ifi "%phase%" == "preloop" $include "./modules/35_transport/edge_esm/preloop.gm $Ifi "%phase%" == "bounds" $include "./modules/35_transport/edge_esm/bounds.gms" $Ifi "%phase%" == "presolve" $include "./modules/35_transport/edge_esm/presolve.gms" $Ifi "%phase%" == "postsolve" $include "./modules/35_transport/edge_esm/postsolve.gms" +$Ifi "%phase%" == "output" $include "./modules/35_transport/edge_esm/output.gms" *######################## R SECTION END (PHASES) ############################### *** EOF ./modules/35_transport/edge_esm.gms diff --git a/modules/45_carbonprice/diffPhaseIn2LinFlex/datainput.gms b/modules/45_carbonprice/diffPhaseIn2LinFlex/datainput.gms deleted file mode 100644 index f5a89931e0007db0ca88091aa9dc19bded6f5af0..0000000000000000000000000000000000000000 --- a/modules/45_carbonprice/diffPhaseIn2LinFlex/datainput.gms +++ /dev/null @@ -1,101 +0,0 @@ -*** | (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 -*** SOF ./modules/45_carbonprice/diffPhaseIn2LinFlex/datainput.gms -***------------------------------------------------------------------------------------------------------------------------ -*** *BS* 20190930 linear convergence with starting points differentiated by GDP/capita, global price from 2040 -***----------------------------------------------------------------------------------------------------------------------- - -*** can make this flexible later -s45_stagestart = 2020; - -*** price from stageend onwards (value set here is for first iteration only, will be adjusted afterwards) -s45_constantCO2price = 500 * sm_DptCO2_2_TDpGtC; - -*** convergence to global CO2 price depends on GDP per capita (in 1e3 $ PPP 2005). -p45_gdppcap2015_PPP(regi) = pm_gdp("2015",regi)/pm_shPPPMER(regi) / pm_pop("2015",regi); -display p45_gdppcap2015_PPP; -*** suggestion by Robert: differentiate "zero-crossing" of linear convergence path -*** earlier zero-crossing --> higher starting price during convergence period -*** for now GDP/cap differentiation hardcoded -*** BS: modified limits to have SSA in first category -p45_phasein_zeroyear(regi)$(p45_gdppcap2015_PPP(regi) le 3) = 2024; -p45_phasein_zeroyear(regi)$(p45_gdppcap2015_PPP(regi) gt 3 and p45_gdppcap2015_PPP(regi) le 5) = 2023; -p45_phasein_zeroyear(regi)$(p45_gdppcap2015_PPP(regi) gt 5 and p45_gdppcap2015_PPP(regi) le 8) = 2022; -p45_phasein_zeroyear(regi)$(p45_gdppcap2015_PPP(regi) gt 8 and p45_gdppcap2015_PPP(regi) le 11) = 2021; -p45_phasein_zeroyear(regi)$(p45_gdppcap2015_PPP(regi) gt 11 and p45_gdppcap2015_PPP(regi) le 14) = 2020; -p45_phasein_zeroyear(regi)$(p45_gdppcap2015_PPP(regi) gt 14 and p45_gdppcap2015_PPP(regi) le 19) = 2018; -p45_phasein_zeroyear(regi)$(p45_gdppcap2015_PPP(regi) gt 19 and p45_gdppcap2015_PPP(regi) le 24) = 2016; -p45_phasein_zeroyear(regi)$(p45_gdppcap2015_PPP(regi) gt 24 and p45_gdppcap2015_PPP(regi) le 30) = 2013; -p45_phasein_zeroyear(regi)$(p45_gdppcap2015_PPP(regi) gt 30) = 2010; -display p45_phasein_zeroyear; - -p45_phasein_2025ratio(regi)$(p45_gdppcap2015_PPP(regi) le 3) = 0.1; -p45_phasein_2025ratio(regi)$(p45_gdppcap2015_PPP(regi) gt 3 and p45_gdppcap2015_PPP(regi) le 5) = 0.2; -p45_phasein_2025ratio(regi)$(p45_gdppcap2015_PPP(regi) gt 5 and p45_gdppcap2015_PPP(regi) le 8) = 0.3; -p45_phasein_2025ratio(regi)$(p45_gdppcap2015_PPP(regi) gt 8 and p45_gdppcap2015_PPP(regi) le 11) = 0.5; -p45_phasein_2025ratio(regi)$(p45_gdppcap2015_PPP(regi) gt 11 and p45_gdppcap2015_PPP(regi) le 14) = 0.65; -p45_phasein_2025ratio(regi)$(p45_gdppcap2015_PPP(regi) gt 14 and p45_gdppcap2015_PPP(regi) le 19) = 0.8; -p45_phasein_2025ratio(regi)$(p45_gdppcap2015_PPP(regi) gt 19 and p45_gdppcap2015_PPP(regi) le 24) = 0.9; -p45_phasein_2025ratio(regi)$(p45_gdppcap2015_PPP(regi) gt 24) = 1; -display p45_phasein_2025ratio; - - - - - -*** get CO2 price before transition stage from reference (NDC) run -***Execute_Loadpoint 'input_ref' p45_tauCO2_ref = pm_taxCO2eq; -***pm_taxCO2eq(ttot,regi)$(ttot.val le s45_stagestart) = p45_tauCO2_ref(ttot,regi); -***display p45_tauCO2_ref; - -*** for the current implementation, use the following trajectory for rich countries: -*** global price is linear from 2010 until the pkBudgYr, then increases with cm_taxCO2inc_after_peakBudgYr -if(cm_co2_tax_2020 lt 0, - abort "please choose a valid cm_co2_tax_2020" -elseif cm_co2_tax_2020 ge 0, -*** convert tax value from $/t CO2eq to T$/GtC - p45_CO2priceTrajDeveloped("2040")= 3 * cm_co2_tax_2020 * sm_DptCO2_2_TDpGtC; !! shifted to 2040 to make sure that even in delay scenarios the fixpoint of the linear price path is inside the "t" range, otherwise the CO2 prices from reference run may be overwritten -*** The factor 3 comes from shifting the 2020 value 20 years into the future at linear increase of 10% of 2020 value per year. -); - - - -p45_CO2priceTrajDeveloped(ttot)$((ttot.val gt 2005) AND (ttot.val ge cm_startyear)) = p45_CO2priceTrajDeveloped("2040")*( 1 + 0.1/3 * (ttot.val-2040)); !! no CO2 price in 2005 and only change CO2 prices after ; -*** annual increase by (10/3)% of the 2040 value is the same as a 10% increase of the 2020 value is the same as a linear increase from 0 in 2010 to the 2020/2040 value - -loop(t2$(t2.val eq cm_peakBudgYr), - p45_CO2priceTrajDeveloped(t)$(t.val gt cm_peakBudgYr) = p45_CO2priceTrajDeveloped(t2) + (t.val - t2.val) * cm_taxCO2inc_after_peakBudgYr * sm_DptCO2_2_TDpGtC; !! increase by cm_taxCO2inc_after_peakBudgYr per year -); - -*** Then create regional phase-in: -*** loop(ttot$((ttot.val ge cm_startyear) AND (ttot.val le cm_CO2priceRegConvEndYr) ), -*** p45_regCO2priceFactor(ttot,regi) = max(0, p45_phasein_2025ratio(regi) + (1-p45_phasein_2025ratio(regi)) * (ttot.val - 2025) / (cm_CO2priceRegConvEndYr - 2025) ); -*** ); -*** p45_regCO2priceFactor(ttot,regi)$(ttot.val ge cm_CO2priceRegConvEndYr) = 1; - -loop(ttot$((ttot.val ge cm_startyear) AND (ttot.val le cm_CO2priceRegConvEndYr) ), - p45_regCO2priceFactor(ttot,regi) = - min(1, - max(0, - p45_phasein_2025ratio(regi) + (1 - p45_phasein_2025ratio(regi)) - * ( ( (ttot.val - 2025) + (cm_CO2priceRegConvEndYr - 2025) * 0.1 ) - / ( (cm_CO2priceRegConvEndYr - 2025) * 1.1 ) - ) ** 2 - ) - ); -); -p45_regCO2priceFactor(ttot,regi)$(ttot.val ge cm_CO2priceRegConvEndYr) = 1; - - -*** linear transition to global price - starting point depends on GDP/cap -pm_taxCO2eq(t,regi) = p45_regCO2priceFactor(t,regi) * p45_CO2priceTrajDeveloped(t); - - - -display p45_regCO2priceFactor, p45_CO2priceTrajDeveloped, pm_taxCO2eq; - -*** EOF ./modules/45_carbonprice/diffPhaseIn2LinFlex/datainput.gms diff --git a/modules/45_carbonprice/diffPhaseIn2LinFlex/declarations.gms b/modules/45_carbonprice/diffPhaseIn2LinFlex/declarations.gms deleted file mode 100644 index d541ed76fd5777ba9f1a31aed8b475effc139d29..0000000000000000000000000000000000000000 --- a/modules/45_carbonprice/diffPhaseIn2LinFlex/declarations.gms +++ /dev/null @@ -1,30 +0,0 @@ -*** | (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 -*** SOF ./modules/45_carbonprice/diffPhaseIn2LinFlex/declarations.gms -***------------------------------------------------------------------------------------------------------------------------ -*** *BS* 20190930 linear convergence with starting points differentiated by GDP/capita, global price from 2040 -***----------------------------------------------------------------------------------------------------------------------- - -parameters -p45_tauCO2_ref(ttot, all_regi) "CO2 tax path of reference policy (NDC)" -p45_gdppcap2015_PPP(all_regi) "2015 GDP per capita (k $ PPP 2005)" -p45_phasein_zeroyear(all_regi) "year when CO2 price convergence line crosses zero" -p45_phasein_2025ratio(all_regi) "ratio of CO2 price to that of developed region in 2025" - -p45_regCO2priceFactor(ttot,all_regi) "regional multiplicative factor to the CO2 price of the developed countries" -p45_CO2priceTrajDeveloped(ttot) "CO2 price trajectory for developed/rich countries" -; - -scalars -s45_stagestart "first time-step fixed to ref. / beginning of staged accession period" - -s45_constantCO2price "initial value for constant global CO2 price" -s45_convergenceCO2price "price to which the regional values converge" - - -; -*** EOF ./modules/45_carbonprice/diffPhaseIn2LinFlex/declarations.gms diff --git a/modules/45_carbonprice/diffPhaseIn2LinFlex/not_used.txt b/modules/45_carbonprice/diffPhaseIn2LinFlex/not_used.txt deleted file mode 100644 index c0d73713a745a09828651bf8fb218affd32af91b..0000000000000000000000000000000000000000 --- a/modules/45_carbonprice/diffPhaseIn2LinFlex/not_used.txt +++ /dev/null @@ -1,23 +0,0 @@ -name,type,reason -sm_c_2_co2,input,questionnaire -vm_co2eq,input,questionnaire -vm_emiFgas,input,questionnaire -vm_cesIO,input,questionnaire -vm_prodSe,input,questionnaire -vm_prodFe,input,questionnaire -pm_pvp,input,questionnaire -pm_globalMeanTemperature,input,questionnaire -pm_temperatureImpulseResponseCO2,input,questionnaire -pm_consPC,input,questionnaire -pm_GDPGross,input,questionnaire -pm_ttot_val,input,questionnaire -pm_ts,input,questionnaire -pm_ttot_2_tall,input,questionnaire -pm_prtp,input,questionnaire -cm_emiscen,input,questionnaire -cm_co2_tax_growth,input,questionnaire -cm_iterative_target_adj,input,questionnaire -cm_expoLinear_yearStart,input,questionnaire -cm_carbonprice_temperatureLimit,input,questionnaire -pm_prodCouple,input,questionnaire -pm_bunker_share_in_nonldv_fe,input,questionnaire diff --git a/modules/45_carbonprice/diffPhaseIn2LinFlex/postsolve.gms b/modules/45_carbonprice/diffPhaseIn2LinFlex/postsolve.gms deleted file mode 100644 index c47766461ad7c31146470fecc38cd7a154bf72ab..0000000000000000000000000000000000000000 --- a/modules/45_carbonprice/diffPhaseIn2LinFlex/postsolve.gms +++ /dev/null @@ -1,32 +0,0 @@ -*** | (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 -*** SOF ./modules/45_carbonprice/diffPhaseIn2LinFlex/postsolve.gms -***------------------------------------------------------------------------------------------------------------------------ -*** *BS* 20190930 linear convergence with starting points differentiated by GDP/capita, global price from 2040 -***----------------------------------------------------------------------------------------------------------------------- - -*** *** updated constant global price as scalar (regional prices are the same anyway) -*** s45_constantCO2price = sum((ttot,regi)$(ttot.val eq s45_stageend), pm_taxCO2eq(ttot,regi))/card(regi) ; -*** *** entire path has been shifted in update, so have to set these again -*** pm_taxCO2eq(ttot,regi)$(ttot.val le s45_stagestart) = p45_tauCO2_ref(ttot, regi); -*** pm_taxCO2eq(ttot,regi)$(ttot.val gt s45_stagestart and ttot.val lt s45_stageend) -*** = s45_constantCO2price * (ttot.val - p45_phasein_zeroyear(regi))/(s45_stageend - p45_phasein_zeroyear(regi)); -*** *** price trajectory should be constant anyway but let's be explicit here -*** pm_taxCO2eq(ttot,regi)$(ttot.val ge s45_stageend) = s45_constantCO2price; - - -*** re-create the regional differentation, use path from developed countries as the basis. -*** This doesn't need to be a loop, but it will be correct for any cycle of the loop, so also for the last cycle. -loop(regi$(p45_gdppcap2015_PPP(regi) gt 30), - p45_CO2priceTrajDeveloped(t) = pm_taxCO2eq(t,regi); -); - -*** linear transition to global price - starting point depends on GDP/cap -pm_taxCO2eq(t,regi) = p45_regCO2priceFactor(t,regi) * p45_CO2priceTrajDeveloped(t); - -display pm_taxCO2eq; -*** EOF ./modules/45_carbonprice/diffPhaseIn2LinFlex/postsolve.gms diff --git a/modules/45_carbonprice/diffPhaseIn2LinFlex/realization.gms b/modules/45_carbonprice/diffPhaseIn2LinFlex/realization.gms deleted file mode 100644 index 27beb4e0bd5a678e32fcca6b48221ebb40104598..0000000000000000000000000000000000000000 --- a/modules/45_carbonprice/diffPhaseIn2LinFlex/realization.gms +++ /dev/null @@ -1,12 +0,0 @@ -*** | (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 - -*####################### R SECTION START (PHASES) ############################## -$Ifi "%phase%" == "declarations" $include "./modules/45_carbonprice/diffPhaseIn2LinFlex/declarations.gms" -$Ifi "%phase%" == "datainput" $include "./modules/45_carbonprice/diffPhaseIn2LinFlex/datainput.gms" -$Ifi "%phase%" == "postsolve" $include "./modules/45_carbonprice/diffPhaseIn2LinFlex/postsolve.gms" -*######################## R SECTION END (PHASES) ############################### diff --git a/modules/45_carbonprice/module.gms b/modules/45_carbonprice/module.gms index 5440ea0b84eb07a30354fa6a67718265a3c0be0c..268947d965c1d5f2725ee1bade6e6bebdec0e459 100644 --- a/modules/45_carbonprice/module.gms +++ b/modules/45_carbonprice/module.gms @@ -23,7 +23,6 @@ $Ifi "%carbonprice%" == "NPi2018" $include "./modules/45_carbonprice/NPi2018/rea $Ifi "%carbonprice%" == "diffCurvPhaseIn2Lin" $include "./modules/45_carbonprice/diffCurvPhaseIn2Lin/realization.gms" $Ifi "%carbonprice%" == "diffPhaseIn2Constant" $include "./modules/45_carbonprice/diffPhaseIn2Constant/realization.gms" $Ifi "%carbonprice%" == "diffPhaseIn2Lin" $include "./modules/45_carbonprice/diffPhaseIn2Lin/realization.gms" -$Ifi "%carbonprice%" == "diffPhaseIn2LinFlex" $include "./modules/45_carbonprice/diffPhaseIn2LinFlex/realization.gms" $Ifi "%carbonprice%" == "diffPhaseInLin2LinFlex" $include "./modules/45_carbonprice/diffPhaseInLin2LinFlex/realization.gms" $Ifi "%carbonprice%" == "diffPriceSameCost" $include "./modules/45_carbonprice/diffPriceSameCost/realization.gms" $Ifi "%carbonprice%" == "exogenous" $include "./modules/45_carbonprice/exogenous/realization.gms" diff --git a/scripts/iterative/EDGE_transport.R b/scripts/iterative/EDGE_transport.R index 887e99369871eb4fd539c8b466b6035f29c7428f..1a7f4d340a6ee1f0b9db5062ef46a1ea0f8fd334 100644 --- a/scripts/iterative/EDGE_transport.R +++ b/scripts/iterative/EDGE_transport.R @@ -1,19 +1,31 @@ -require(data.table) -require(gdx) -require(gdxdt) -require(edgeTrpLib) -require(rmndt) +library(optparse) + +opt_parser = OptionParser( + description = "Coupled version of EDGE-T, to be run within a REMIND output folder.", + option_list = list( + make_option( + "--reporting", action="store_true", default=FALSE, + help="Store output files in subfolder EDGE-T"))); +opt = parse_args(opt_parser); + +library(data.table) +library(gdx) +library(gdxdt) +library(edgeTrpLib) +library(rmndt) +library(moinput) + ## use cached input data for speed purpose -require(moinput) setConfig(forcecache=T) +data_folder <- "EDGE-T" mapspath <- function(fname){ file.path("../../modules/35_transport/edge_esm/input", fname) } datapath <- function(fname){ - file.path("input_EDGE", fname) + file.path(data_folder, fname) } REMINDpath <- function(fname){ @@ -35,12 +47,9 @@ EDGE_scenario <- cfg$gms$cm_EDGEtr_scen EDGEscenarios <- fread("../../modules/35_transport/edge_esm/input/EDGEscenario_description.csv")[scenario_name == EDGE_scenario] -merge_traccs <<- EDGEscenarios[options == "merge_traccs", switch] -addvintages <<- EDGEscenarios[options == "addvintages", switch] -inconvenience <<- EDGEscenarios[options == "inconvenience", switch] -selfmarket_taxes <<- EDGEscenarios[options == "selfmarket_taxes", switch] -selfmarket_policypush <<- EDGEscenarios[options == "selfmarket_policypush", switch] -selfmarket_acceptancy <<- EDGEscenarios[options == "selfmarket_acceptancy", switch] +inconvenience <- EDGEscenarios[options == "inconvenience", switch] +selfmarket_policypush <- EDGEscenarios[options == "selfmarket_policypush", switch] +selfmarket_acceptancy <- EDGEscenarios[options == "selfmarket_acceptancy", switch] if (EDGE_scenario == "Conservative_liquids") { techswitch <<- "Liquids" @@ -53,22 +62,22 @@ if (EDGE_scenario == "Conservative_liquids") { exit() } -endogeff <<- EDGEscenarios[options== "endogeff", switch] -enhancedtech <<- EDGEscenarios[options== "enhancedtech", switch] -rebates_febates <<- EDGEscenarios[options== "rebates_febates", switch] ##NB THEY ARE ONLY IN PSI! ONLY WORKING IN EUROPE -savetmpinput <<- FALSE -smartlifestyle <<- EDGEscenarios[options== "smartlifestyle", switch] - - REMIND2ISO_MAPPING <- fread(REMINDpath(cfg$regionmapping))[, .(iso = CountryCode, region = RegionCode)] EDGE2teESmap <- fread(mapspath("mapping_EDGE_REMIND_transport_categories.csv")) ## input data loading -input_path = paste0("../../modules/35_transport/edge_esm/input/") +input_folder = paste0("../../modules/35_transport/edge_esm/input/") + +if (length(list.files(path = data_folder, pattern = "RDS")) < 7) { + createRDS(input_folder, data_folder, + SSP_scenario = scenario, + EDGE_scenario = EDGE_scenario) +} +inputdata <- loadInputData(data_folder) + -inputdata = createRDS(input_path, SSP_scenario = scenario, EDGE_scenario = EDGE_scenario) vot_data = inputdata$vot_data sw_data = inputdata$sw_data inco_data = inputdata$inco_data @@ -77,27 +86,21 @@ int_dat = inputdata$int_dat nonfuel_costs = inputdata$nonfuel_costs price_nonmot = inputdata$price_nonmot -## add learning optional -setlearning = TRUE -## add optional vintages -addvintages = TRUE ## optional average of prices average_prices = FALSE -## inconvenience costs instead of preference factors -inconvenience = TRUE -if (setlearning | addvintages){ - ES_demand = readREMINDdemand(gdx, REMIND2ISO_MAPPING, EDGE2teESmap, REMINDyears) - ## select from total demand only the passenger sm - ES_demand = ES_demand[sector == "trn_pass",] -} + +ES_demand_all = readREMINDdemand(gdx, REMIND2ISO_MAPPING, EDGE2teESmap, REMINDyears) +## select from total demand only the passenger sm +ES_demand = ES_demand_all[sector == "trn_pass",] + -if (setlearning & file.exists("demand_previousiter.RDS")) { +if (file.exists(datapath("demand_previousiter.RDS"))) { ## load previous iteration number of cars - demand_BEVtmp = readRDS("demand_BEV.RDS") + demand_BEVtmp = readRDS(datapath("demand_BEV.RDS")) ## load previous iteration demand - ES_demandpr = readRDS("demand_previousiter.RDS") + ES_demandpr = readRDS(datapath("demand_previousiter.RDS")) ## calculate non fuel costs and nonfuel_costs = applylearning(gdx,REMINDmapping,EDGE2teESmap, demand_BEVtmp, ES_demandpr) saveRDS(nonfuel_costs, "nonfuel_costs_learning.RDS") @@ -111,6 +114,7 @@ REMIND_prices <- merge_prices( intensity_data = int_dat, nonfuel_costs = nonfuel_costs) + ## save prices ## read last iteration count keys <- c("iso", "year", "technology", "vehicle_type") @@ -122,13 +126,13 @@ iter <- as.vector(gdxrrw::rgdx(gdx, list(name="o_iterationNumber"))$val) REMIND_prices[, iternum := iter] ## save REMIND prices (before dampening) -saveRDS(REMIND_prices, paste0("REMINDprices", iter, ".RDS")) +saveRDS(REMIND_prices, datapath(paste0("REMINDprices", iter, ".RDS"))) if(average_prices){ if(max(unique(REMIND_prices$iternum)) >= 20 & max(unique(REMIND_prices$iternum)) <= 30){ - old_prices <- readRDS(pfile) + old_prices <- readRDS(datapath(pfile)) all_prices <- rbind(old_prices, REMIND_prices) setkeyv(all_prices, keys) ## apply moving avg @@ -138,10 +142,10 @@ if(average_prices){ }else{ all_prices <- REMIND_prices } - saveRDS(all_prices, pfile) + saveRDS(all_prices, datapath(pfile)) ## save REMIND prices (after dampening) - saveRDS(REMIND_prices,paste0("REMINDpricesDampened", iter, ".RDS")) + saveRDS(REMIND_prices, datapath(paste0("REMINDpricesDampened", iter, ".RDS"))) } @@ -157,7 +161,9 @@ if (inconvenience) { inco_data = inco_data, logit_params = logit_params, intensity_data = int_dat, - price_nonmot = price_nonmot) + price_nonmot = price_nonmot, + selfmarket_policypush = selfmarket_policypush, + selfmarket_acceptancy = selfmarket_acceptancy) } else{ @@ -177,47 +183,65 @@ shares <- logit_data[["share_list"]] ## shares of alternatives for each level of mj_km_data <- logit_data[["mj_km_data"]] ## energy intensity at a technology level prices <- logit_data[["prices_list"]] ## prices at each level of the logit function, 1990USD/pkm -if(addvintages){ - ## calculate vintages (new shares, prices, intensity) - vintages = calcVint(shares = shares, - totdem_regr = ES_demand, - prices = prices, - mj_km_data = mj_km_data, - years = REMINDyears) - - shares$FV_shares = vintages[["shares"]]$FV_shares - prices = vintages[["prices"]] - mj_km_data = vintages[["mj_km_data"]] -} + +## calculate vintages (new shares, prices, intensity) +vintages = calcVint(shares = shares, + totdem_regr = ES_demand, + prices = prices, + mj_km_data = mj_km_data, + years = REMINDyears) + +shares$FV_shares = vintages[["shares"]]$FV_shares +prices = vintages[["prices"]] +mj_km_data = vintages[["mj_km_data"]] + ## use logit to calculate shares and intensities (on tech level) EDGE2CESmap <- fread(mapspath("mapping_CESnodes_EDGE.csv")) -shares_intensity_demand <- shares_intensity_and_demand( - logit_shares=shares, - MJ_km_base=mj_km_data, - EDGE2CESmap=EDGE2CESmap, - REMINDyears=REMINDyears, - scenario=scenario, - REMIND2ISO_MAPPING=REMIND2ISO_MAPPING) -demByTech <- shares_intensity_demand[["demand"]] ##in [-] -intensity <- shares_intensity_demand[["demandI"]] ##in million pkm/EJ -norm_demand <- shares_intensity_demand$demandF_plot_pkm ## total demand is 1, required for costs - -if (setlearning) { - demand_BEV=calc_num_vehicles( norm_dem_BEV = norm_demand[technology == "BEV" & ## battery vehicles - subsector_L1 == "trn_pass_road_LDV_4W", ## only 4wheelers - c("iso", "year", "sector", "vehicle_type", "demand_F") ], - ES_demand = ES_demand) - - ## save number of vehicles for next iteration - saveRDS(demand_BEV, "demand_BEV.RDS") - ## save the demand for next iteration renaming the column - setnames(ES_demand, old ="demand", new = "demandpr") - saveRDS(ES_demand, "demand_previousiter.RDS") +shares_int_dem <- shares_intensity_and_demand( + logit_shares=shares, + MJ_km_base=mj_km_data, + EDGE2CESmap=EDGE2CESmap, + REMINDyears=REMINDyears, + scenario=scenario, + REMIND2ISO_MAPPING=REMIND2ISO_MAPPING, + demand_input = if (opt$reporting) ES_demand_all) + +demByTech <- shares_int_dem[["demand"]] ##in [-] +intensity <- shares_int_dem[["demandI"]] ##in million pkm/EJ +norm_demand <- shares_int_dem[["demandF_plot_pkm"]] ## total demand is 1, required for costs + + +if (opt$reporting) { + saveRDS(vintages[["vintcomp"]], file = datapath("vintcomp.RDS")) + saveRDS(vintages[["newcomp"]], file = datapath("newcomp.RDS")) + saveRDS(shares, file = datapath("shares.RDS")) + saveRDS(logit_data$EF_shares, file = datapath("EF_shares.RDS")) + saveRDS(logit_data$mj_km_data, file = datapath("mj_km_data.RDS")) + saveRDS(logit_data$inconv_cost, file=datapath("inco_costs.RDS")) + saveRDS(shares_int_dem$demandF_plot_EJ, + file=datapath("demandF_plot_EJ.RDS")) + saveRDS(shares_int_dem$demandF_plot_pkm, + datapath("demandF_plot_pkm.RDS")) + saveRDS(logit_data$annual_sales, file = datapath("annual_sales.RDS")) + quit() } +demand_BEV=calc_num_vehicles( + norm_dem_BEV = norm_demand[ + technology == "BEV" & ## battery vehicles + subsector_L1 == "trn_pass_road_LDV_4W", ## only 4wheelers + c("iso", "year", "sector", "vehicle_type", "demand_F") ], + ES_demand = ES_demand) + +## save number of vehicles for next iteration +saveRDS(demand_BEV, datapath("demand_BEV.RDS")) +## save the demand for next iteration renaming the column +setnames(ES_demand, old ="demand", new = "demandpr") +saveRDS(ES_demand, datapath("demand_previousiter.RDS")) + ## use logit to calculate costs budget <- calculate_capCosts( diff --git a/scripts/output/single/EDGETransportValidation.R b/scripts/output/single/EDGETransportValidation.R index 927aadbfc68fab3a0d7738efa721e70521bf3df4..2fdfab5fda67f8d52cd162448da9c8530501d319 100644 --- a/scripts/output/single/EDGETransportValidation.R +++ b/scripts/output/single/EDGETransportValidation.R @@ -1,4 +1,19 @@ +# | (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) + +if(!exists("source_include")) { + ## Define arguments that can be read from command line + readArgs("outputdir","gdx_name","gdx_ref_name") +} + +load(file.path(outputdir, "config.Rdata")) ## run EDGE transport validation output if required if(cfg$gms$transport == "edge_esm"){ diff --git a/scripts/output/single/notebook_templates/EDGETransportReport.Rmd b/scripts/output/single/notebook_templates/EDGETransportReport.Rmd index 20e7b74c78da0e3abc926c4181f50af9a0a0feee..4286f64d6baf73942208786c49d5af4eabf8df2c 100644 --- a/scripts/output/single/notebook_templates/EDGETransportReport.Rmd +++ b/scripts/output/single/notebook_templates/EDGETransportReport.Rmd @@ -14,645 +14,426 @@ require(remind) require(gdxdt) require(gdx) require(rmndt) -require(data.table) - +require(magclass) +require(quitte) +require(ggpubr) +require(gridExtra) require(edgeTrpLib) ``` ```{r, echo=FALSE, warning=FALSE} +iso_plot = "DEU" +output_folder = "EDGE-T/" -dem_shares <- list() -intensity <- list() -demand_km <- list() -demand_ej <- list() -sw_tech <- list() -prices_FV <- list() +cols <- c("NG" = "#d11141", + "Liquids" = "#8c8c8c", + "Hybrid Liquids" = "#ffc425", + "Hybrid Electric" = "#f37735", + "BEV" = "#00b159", + "FCEV" = "#00aedb") datapath <- function(fname){ - file.path("./input_EDGE/", fname) + file.path("./EDGE-T/", fname) } mapspath <- function(fname, scenariopath=""){ file.path("../../modules/35_transport/edge_esm/input", fname) } + ## Load mappings -EDGE2CESmap <- fread(mapspath("mapping_CESnodes_EDGE.csv")) +REMIND2ISO_MAPPING <- fread("../../config/regionmappingH12.csv")[, .(iso = CountryCode, region = RegionCode)] +EDGE2teESmap <- fread(mapspath("mapping_EDGE_REMIND_transport_categories.csv")) -REMIND2ISO_MAPPING <- fread("../../config/regionmappingH12.csv")[, .(iso = CountryCode, - region = RegionCode)] +## load input data from last EDGE run +demand_km <- readRDS(datapath(fname = "demandF_plot_pkm.RDS")) ## detailed energy services demand, million km +demand_ej <- readRDS(datapath(fname = "demandF_plot_EJ.RDS")) ## detailed final energy demand, EJ +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")) +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")) + +name_mif = list.files(pattern = "REMIND_generic", full.names = F) +name_mif = name_mif[!grepl("withoutPlu", name_mif)] +miffile <- as.data.table(read.quitte(name_mif)) +``` -EDGE2teESmap <- fread(mapspath("mapping_EDGE_REMIND_transport_categories.csv")) +# LDVs vintages -years <- c(1990, - seq(2005, 2060, by = 5), - seq(2070, 2110, by = 10), - 2130, 2150) - -REMINDyears <- c(1990, - seq(2005, 2060, by = 5), - seq(2070, 2110, by = 10), - 2130, 2150) - -## include the paths to the scenarios you want to compare -output_folders <- "./" - -for(output_folder in output_folders){ - ## load gdx for fuel prices and demand - gdx = file.path(output_folder, "fulldata.gdx") - ## load policy scenario - load(file.path(output_folder, "config.Rdata")) - REMIND_scenario <- cfg$gms$cm_GDPscen - EDGE_scenario <- cfg$gms$cm_EDGEtr_scen - policy_scenario <- cfg$gms$c_expname - - scen <- paste0(REMIND_scenario, "-", EDGE_scenario, "-", policy_scenario) - - ## load demand - ES_demand = readREMINDdemand(gdx, REMIND2ISO_MAPPING, EDGE2teESmap, years) - - ## load input data - int_dat <- readRDS(datapath("harmonized_intensities.RDS")) - nonfuel_costs <- readRDS(datapath("UCD_NEC_iso.RDS")) - sw_data <- readRDS(datapath("SW.RDS")) - vot_data <- readRDS(datapath("VOT_iso.RDS")) - logit_params <- readRDS(datapath("logit_exp.RDS")) - price_nonmot <- readRDS(datapath("price_nonmot.RDS")) - - ## FIXME: hotfix to make the (empty) vot_data$value_time_VS1 with the right column types. Probably there is another way to do that, did not look for it. - vot_data$value_time_VS1$iso = as.character(vot_data$value_time_VS1$iso) - vot_data$value_time_VS1$subsector_L1 = as.character(vot_data$value_time_VS1$subsector_L1) - vot_data$value_time_VS1$vehicle_type = as.character(vot_data$value_time_VS1$vehicle_type) - vot_data$value_time_VS1$year = as.numeric(vot_data$value_time_VS1$year) - vot_data$value_time_VS1$time_price = as.numeric(vot_data$value_time_VS1$time_price) - - ## calculate prices - REMIND_prices <- merge_prices( - gdx = gdx, - REMINDmapping = REMIND2ISO_MAPPING, - REMINDyears = REMINDyears, - intensity_data = int_dat, - nonfuel_costs = nonfuel_costs) - - ## calculates logit - logit_data <- calculate_logit( - REMIND_prices[tot_price > 0], - REMIND2ISO_MAPPING, - vot_data = vot_data, - sw_data = sw_data, - logit_params = logit_params, - intensity_data = int_dat, - price_nonmot = price_nonmot) - - - shares <- logit_data[["share_list"]] ## shares of alternatives for each level of the logit function - mj_km_data <- logit_data[["mj_km_data"]] ## energy intensity at a technology level - prices_FV[[scen]] <- REMIND_prices[, EDGE_scenario := scen] ## prices at each level of the logit function, 1990USD/pkm - - ## calculate energy intensity and FE demand at a REMIND-region level for the desired level of aggregation - res <- shares_intensity_and_demand( - logit_shares=shares, - MJ_km_base=mj_km_data, - REMIND2ISO_MAPPING=REMIND2ISO_MAPPING, - EDGE2CESmap=EDGE2CESmap, - REMINDyears=REMINDyears, - demand_input=ES_demand) - - dem_shares[[scen]] <- res$demand[, EDGE_scenario := scen] - intensity[[scen]] <- res$demandI[, EDGE_scenario := scen] - demand_km[[scen]] <- res$demandF_plot_pkm[, EDGE_scenario := scen] - demand_ej[[scen]] <- res$demandF_plot_EJ[, EDGE_scenario := scen] - sw_tech[[scen]] <- sw_data$FV_final_SW[, EDGE_scenario := scen] -} +```{r, echo=FALSE, warning=FALSE} -dem_shares <- rbindlist(dem_shares) -intensity <- rbindlist(intensity) -demand_km <- rbindlist(demand_km) -demand_ej <- rbindlist(demand_ej) -sw_tech <- rbindlist(sw_tech) -prices_FV <- rbindlist(prices_FV) -``` +plotVint = 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 = 2 + 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)] + + + p = ggplot()+ + geom_bar(data = allfleet[year %in% c(2015,2030,2050)], + aes(x=as.character(year),y=value, group=interaction(variable, technology), + fill = technology), alpha = 0.5, position="stack", stat = "identity", width = 0.5)+ + geom_bar(data = allfleet[year %in% c(2015,2030,2050)], + aes(x=as.character(year),y=value, group=interaction(variable, technology), + fill = technology, alpha = factor(alphaval)), position="stack", stat = "identity", width = 0.5, color = "black")+ + guides(fill = guide_legend(reverse=TRUE))+ + theme_minimal()+ + facet_wrap(~region, nrow = 4)+ + 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))+ + scale_x_discrete(breaks = c(2015,2030,2050))+ + 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 = "LDV fleet [million Veh]", x="") + return(p) + } -```{r, echo=FALSE} -## plot settings -years_plot = c(2010,2015,2020,2025,2030,2040,2050) ## in bar charts, these are the time steps that are represented -year_single = 2050 -region_plot = "NEU" ## in case is a region specific plot, this region is represented -sector_plot ="trn_pass" ## in case is a sector specific plot, this sector is represented +p = plotVint(vintcomp, newcomp, shares$VS1_shares) +p -##conversion rate 2005->1990 USD -CONV_2005USD_1990USD=0.67 +``` -# print(paste0("Scenario: ", REMIND_scenario)) -print(paste0("Regional plots are about ",region_plot)) -print(paste0("Sectoral plots are about ",sector_plot)) +# Inconvenience cost trend -## maps -cesmap <- data.table(CES_parent=c("_p_sm", "_p_lo", "_f_lo", "_f_sm"), - CES_label=c("Passenger, Short-to-Medium Distances", - "Passenger, Long Distances", - "Freight, Long Distances", - "Freight, Short-to-Medium Distances")) +```{r, echo=FALSE, warning=FALSE} -EDGE_sectormap <- data.table(sector=c("trn_pass", "trn_freight", "trn_aviation_intl", "trn_shipping_intl"), - CES_label=c("Passenger, Short-to-Medium Distances", - "Passenger, Long Distances", - "Freight, Long Distances", - "Freight, Short-to-Medium Distances")) +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")+ + facet_grid(~technology)+ + theme_minimal()+ + scale_fill_manual(values = cols)+ + expand_limits(y = c(0,0.8))+ + 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 ``` +# Endogenous intensity for Liquids -```{r, echo=FALSE} -## aggregate demands to REMIND regions -demandF_plot_EJ <- demand_ej[,c("EDGE_scenario", "sector","subsector_L3","subsector_L2", - "subsector_L1","vehicle_type","technology", "iso","year","demand_EJ")] -demandF_plot_pkm <- demand_km[,c("EDGE_scenario", "sector","subsector_L3","subsector_L2", - "subsector_L1","vehicle_type","technology","iso","year","demand_F")] - -demandF_plot_EJ=aggregate_dt(demandF_plot_EJ,REMIND2ISO_MAPPING, - datacols = c("EDGE_scenario", "sector", "subsector_L3", "subsector_L2", "subsector_L1", - "vehicle_type", "technology"), - valuecol = "demand_EJ") - -demandF_plot_pkm=aggregate_dt(demandF_plot_pkm,REMIND2ISO_MAPPING, - datacol = c("EDGE_scenario", "sector","subsector_L3","subsector_L2", - "subsector_L1","vehicle_type","technology"), - valuecol = "demand_F") -## add GLO -glo <- demandF_plot_pkm[,.(region="GLO", demand_F=sum(demand_F)), - by=eval(names(demandF_plot_pkm)[2:9])] -demandF_plot_pkm <- rbind(demandF_plot_pkm, glo) -glo <- demandF_plot_EJ[,.(region="GLO", demand_EJ=sum(demand_EJ)), - by=eval(names(demandF_plot_EJ)[2:9])] -demandF_plot_EJ <- rbind(demandF_plot_EJ, glo) +```{r, echo=FALSE, message=FALSE, warning=FALSE} +## Choice of the energy intensity (of the new sales) +intcompplotf = function(EF_shares, FV_shares, VS1_shares){ + EF_shares = EF_shares[,c("iso", "year", "technology", "vehicle_type", "subsector_L1", "subsector_L2", "subsector_L3", "sector", "share","type")] + setnames(EF_shares, old="share", new = "shareINT") + FV_shares = FV_shares[iso == iso_plot & subsector_L1 == "trn_pass_road_LDV_4W" & technology == "Liquids"] + setnames(FV_shares, old="share", new = "shareF") + VS1_shares = VS1_shares[iso == iso_plot & subsector_L1 == "trn_pass_road_LDV_4W"] + shares_LDV = merge(FV_shares, EF_shares, all = FALSE, by = c("iso", "year", "technology", "vehicle_type", "subsector_L1")) + + shares_LDV[, shareIF := shareF*shareINT] + shares_LDV <- shares_LDV[,.(shareIF=sum(shareIF)),by=c("iso","technology","type","vehicle_type","subsector_L1", "year")] + + shares_LDV = merge(shares_LDV, VS1_shares, all = TRUE, by = c("iso", "year", "vehicle_type", "subsector_L1")) + shares_LDV[, shareIS1 := shareIF*share] + shares_LDV <- shares_LDV[,.(shareIS1=sum(shareIS1)),by=c("iso","type", "technology","subsector_L1","year")] + p = ggplot()+ + geom_bar(data = shares_LDV[year<=2100 & year>=2025], aes(x=year,y=shareIS1, group = technology, fill = technology), alpha = 0.5, position = position_fill(), stat = "identity")+ + geom_bar(data = shares_LDV[year<=2100 & year>=2025], aes(x=year,y=shareIS1, group = technology, fill = technology, alpha = type), position = position_fill(), stat = "identity")+ + facet_wrap("technology")+ + theme_minimal()+ + expand_limits(y = c(0,1))+ + scale_fill_manual("technology", values = cols)+ + scale_alpha_discrete("Type")+ + labs(y = "Share [%]", title = paste0("Energy intensity new sales of Liquids, example for ", iso_plot)) -``` + return(p) -## ES - - -```{r, echo=FALSE} -##chunk of code that plots the ES -ES_modes_bar1=function(demandpkm){ - #group by subsector_L3 and summarise the demand - df=demandpkm[, .(demand_F=sum(demand_F) - ), by = c("EDGE_scenario", "region", "year","sector","subsector_L1")] - df[,demand_F:=demand_F ## in millionkm - *1e-6 ## in trillion km - ] - df=df[order(year)] - # #filter only 2020, 2050 and 2100 - df=df[year %in% years_plot,] - #separate into passenger and freight categories - pass=c("trn_pass","trn_aviation_intl") - freight=c("trn_freight","trn_shipping_intl") - ## give proper names to the categories - df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1") - #plot - plot_p=ggplot()+ - geom_bar(data=df%>%filter(sector %in% pass, year %in% years_plot, region == region_plot), - aes(x=year,y=demand_F,group=mode,fill=mode),position=position_stack(),stat="identity")+ - facet_wrap(~EDGE_scenario)+ - ggtitle("Energy Services Demand - Passenger Transport Modes")+ - theme_light()+ - theme(axis.text.x = element_text(angle = 90, hjust = 1))+ - scale_x_continuous(breaks=years_plot)+ - xlab("Year")+ - ylab("Energy Services Demand (trillion pkm)")+ - guides(fill=guide_legend(title="Transport mode"))+ - theme(axis.text.x = element_text(angle = 90, hjust = 1), - axis.text = element_text(size=13), - title = element_text(size=13), - legend.text = element_text(size=13)) - - plot_f=ggplot()+ - geom_bar(data=df%>%filter(sector %in% freight,year %in% years_plot, region == region_plot), - aes(x=year,y=demand_F,group=mode,fill=mode),position=position_stack(),stat="identity")+ - facet_wrap(~EDGE_scenario)+ - ggtitle("Energy Services Demand - Freight Transport Modes")+ - theme_light()+ - theme(axis.text.x = element_text(angle = 90, hjust = 1))+ - scale_x_continuous(breaks=years_plot)+ - xlab("Year")+ - ylab("Energy Services Demand (trillion tkm)")+ - guides(fill=guide_legend(title="Transport mode"))+ - theme(axis.text.x = element_text(angle = 90, hjust = 1), - axis.text = element_text(size=13), - title = element_text(size=13), - legend.text = element_text(size=13)) - - plot=list(plot_p,plot_f) - return(plot) } -p=ES_modes_bar1(demandpkm=demandF_plot_pkm) -p[[1]] -p[[2]] - +intcompplotf(EF_shares, shares$FV_shares, shares$VS1_shares) ``` -```{r, echo=FALSE} -##chunk of code that plots the ES -ES_modes_bar=function(demandpkm){ - demandpkm[technology == "LA-BEV", technology := "BEV"] - ## use proper non-fuel mode names - demandpkm[technology %in% c("Cycle_tmp_technology", "Walk_tmp_technology"), technology := "Human Powered"] - #group by subsector_L3 and summarise the demand - df=demandpkm[, .(demand_F=sum(demand_F)), - by = c("EDGE_scenario", "region", "year","sector","subsector_L1", "technology")] - df[,demand_F:=demand_F ## in millionkm - *1e-6 ## in trillion km - ] - df=df[order(year)] - #separate into passenger and freight categories - pass=c("trn_pass","trn_aviation_intl") - freight=c("trn_freight","trn_shipping_intl") - ## give proper names to the categories - df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1") - #plot - plot_p=ggplot()+ - geom_bar(data=df%>%filter(sector %in% pass,year %in% years_plot,region==region_plot), - aes(x=year,y=demand_F,group=technology,fill=technology), - position=position_stack(),stat="identity")+ - ggtitle("Energy Services Demand by Technology, Passenger Transport, EUR")+ - theme_light()+ - facet_wrap(~EDGE_scenario)+ - scale_x_continuous(breaks=years_plot)+ - xlab("Year")+ - ylab("Energy Services Demand (trillion pkm)")+ - guides(fill=guide_legend(title="Technology"))+ - theme(axis.text.x = element_text(angle = 90, hjust = 1), - axis.text = element_text(size=13), - title = element_text(size=13), - legend.text = element_text(size=13)) - - plot_f=ggplot()+ - geom_bar(data=df%>%filter(sector %in% freight,year %in% years_plot, region==region_plot), - aes(x=year,y=demand_F, group=technology, fill=technology), - position=position_stack(),stat="identity")+ - ggtitle(paste0("Energy Services Demand by Technology, Freight Transport, ", region_plot))+ - facet_wrap(~EDGE_scenario) + - theme_light()+ - scale_x_continuous(breaks=years_plot)+ - xlab("Year")+ - ylab("Energy Services Demand (trllion tkm)")+ - guides(fill=guide_legend(title="Technology"))+ - theme(axis.text.x = element_text(angle = 90, hjust = 1), - axis.text = element_text(size=13), - title = element_text(size=13), - legend.text = element_text(size=13)) - - plot=list(plot_p,plot_f) - return(plot) +# 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<=2050 & year>=2015 & 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))+ + 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) } -p=ES_modes_bar(demandpkm=demandF_plot_pkm) -p[[1]] -p[[2]] +salesplot(annual_sales) ``` -```{r, echo=FALSE} -## plot ES for LDVs only divided by fuel -ES_modes_LDV_bar=function(demandpkm){ - demandpkm[technology == "LA-BEV", technology := "BEV"] - ## use proper non-fuel mode names - demandpkm[technology %in% c("Cycle_tmp_technology", "Walk_tmp_technology"), technology := "Human Powered"] - #group by subsector_L3 and summarise the demand - df=demandpkm[, .(demand_F=sum(demand_F)), - by = c("EDGE_scenario", "region", "year","sector","subsector_L1", "technology")] - df[,demand_F:=demand_F ## in millionkm - *1e-6 ## in trillion km - ] - df=df[order(year)] - ## give proper names to the categories - df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1") - - - ## select order of facets - df$technology = factor(df$technology, levels=c("Liquids","Hybrid Liquids","NG","BEV","FCEV")) - - #plot - plot_LDV=ggplot()+ - geom_bar(data=df%>%filter(year %in% years_plot,region==region_plot, mode %in% c("4W","2W")), - aes(x=year,y=demand_F,group=technology,fill=technology), alpha = 0.9, - position=position_stack(),stat="identity")+ - ggtitle(paste0("Energy Services Demand by Technology, LDVs", region_plot))+ - theme_light()+ - facet_wrap(~EDGE_scenario)+ - theme(axis.text.x = element_text(angle = 90, hjust = 1))+ - xlab("Year")+ - ylab("Energy Services Demand (trillion pkm)")+ - guides(fill=guide_legend(title="Technology"))+ - theme(axis.text.x = element_text(angle = 90, hjust = 1), - axis.text = element_text(size=13), - title = element_text(size=13), - legend.text = element_text(size=13), - strip.text.x = element_text(size = 13, color = "black"), - strip.background=element_rect(fill="white"))+ - scale_x_continuous(breaks=years_plot)+ - scale_fill_brewer(palette = "Set1") - - - return(plot_LDV) -} +# Final energy demand -p=ES_modes_LDV_bar(demandpkm=demandF_plot_pkm) +```{r, echo=FALSE, warning=FALSE} +demandEJplotf = function(demandEJ){ + ## EDGE results + demandEJ <- demandEJ[, c("sector", "subsector_L3", "subsector_L2", "subsector_L1", "vehicle_type", "technology", "iso", "year", "demand_EJ")] + + ## attribute aggregated mode and vehicle names for plotting purposes, and aggregate + 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(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("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("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", "Trucks (<3.5t)", "Trucks (3.5t-16)", "Truck (>12t)", "Trucks (>16)", "Trucks","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", + "Motorbikes", "Small Cars", "Large Cars", "Van", + "International Aviation", "Domestic Aviation","Bus", "Passenger Rail", + "Freight", "LDV", "Pass non LDV", "Freight (Inland)", "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")] -p + p=ggplot()+ + geom_area(data = demandEJ[year > 2010], aes(x=year, y=demand_EJ, group = interaction(vehicle_type_plot,aggr_mode), fill = vehicle_type_plot), color = "black", position= position_stack())+ + facet_wrap(~region, nrow = 4) + labs(x = "", y = "Final Energy demand [EJ]")+ + theme_minimal()+ + # scale_fill_manual("Vehicle Type",values = cols, breaks=legend_ord)+ + theme(axis.text.x = element_text(size = 8), + axis.text.y = element_text(size=8), + axis.title = element_text(size = 9), + title = element_text(size = 9), + legend.text = element_text(size = 9), + legend.title = element_text(size =9), + strip.text = element_text(size=9)) -``` -## FE - -```{r, echo=FALSE} -FE_modes_bar=function(demandEJ){ - #group by subsector_L1 and summarise the demand - df=demandEJ[, .(demand_EJ=sum(demand_EJ)), - by = c("EDGE_scenario", "region", "year","subsector_L1","subsector_L3","sector")] - df=df[order(year)] - df=df[year %in% years_plot,] - #separate into passenger and freight categories - pass=c("trn_pass","trn_aviation_intl") - freight=c("trn_freight","trn_shipping_intl") - ## give proper names to the categories - df <- merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1") - #plot - plot_p=ggplot()+ - geom_bar(data=df%>%filter(sector %in% pass, region==region_plot), - aes(x=year,y=demand_EJ,group=mode,fill=mode),position=position_stack(),stat="identity",color="black")+ - facet_wrap(~EDGE_scenario)+ - ylab("Energy (EJ)") + - ggtitle(paste0("Final Energy Demand by Mode, Passenger, ", region_plot))+ - theme(axis.text.x = element_text(angle = 90))+ - scale_x_continuous(breaks=years_plot) - - plot_f=ggplot()+ - geom_bar(data=df%>%filter(sector %in% freight,region==region_plot), - aes(x=year,y=demand_EJ,group=mode,fill=mode),position=position_stack(),stat="identity",color="black")+ - facet_wrap(~EDGE_scenario) + - ggtitle(paste0("Final Energy Demand by Mode, Freight, ", region_plot))+ - ylab("Energy (EJ)") + - theme(axis.text.x = element_text(angle = 90))+ - scale_x_continuous(breaks=years_plot) - - plot=list(plot_p,plot_f) - return(plot) + return(p) } -p=FE_modes_bar(demandEJ = demandF_plot_EJ) -p[[1]] -p[[2]] +## Final Energy demand +demandEJplotf(demand_ej) ``` -```{r, echo=FALSE} -## function that calculates FE split and splits out the liquids by source -FE_modes_bar_oilcomponent=function(demandEJ, msect="trn_pass", region_plot){ - if(msect == "trn_pass") - sector_display = "Passenger" - if(msect == "trn_freight") - sector_display = "Freight" - - #group by subsector_L1 and summarise the demand - df=demandEJ[, .(demand_EJ=sum(demand_EJ)), - by = c("EDGE_scenario", "region", "year","subsector_L1","subsector_L3","sector","technology")] - df=df[order(year)] - df=df[year %in% years_plot,] - ## give proper names to the categories - df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1") - - df[,technology := ifelse(technology == "LA-BEV", "BEV", technology)] - df[,technology := ifelse(technology == "Electric", "El. Trains", technology)] - ## select order of facerts - df$technology = factor(df$technology, levels=c("BEV","FCEV","Hybrid Liquids", "El. Trains", "NG","Liquids", "Coal")) - - #plot - plot_psm = ggplot()+ - geom_bar(data=df%>%filter(sector == msect, region == region_plot), - aes(x=year,y=demand_EJ,group=technology,fill=technology), - position=position_stack(),stat="identity", alpha = 0.9)+ - theme_light()+ - ggtitle(paste0("Final Energy Demand by Tech, ", sector_display, ", ", region_plot))+ - theme(axis.text.x = element_text(angle = 90), - strip.text.x = element_text(size = 13, color = "black"))+ - scale_x_continuous(breaks=years_plot)+ - scale_fill_brewer(palette = "Set1")+ - xlab("Year")+ - ylab("Final energy demand [EJ]")+ - facet_wrap(~EDGE_scenario) + - guides(fill=guide_legend(title="Technology")) - - plot_psm_LDV = ggplot()+ - geom_bar(data=df%>%filter(sector == msect, region == region_plot, mode == "4W"), - aes(x=year,y=demand_EJ,group=technology,fill=technology), - position=position_stack(),stat="identity", alpha = 0.9)+ - theme_light()+ - ggtitle(paste0("Final Energy Demand by Tech, LDVs, ", region_plot))+ - theme(axis.text.x = element_text(angle = 90), - strip.text.x = element_text(size = 13, color = "black"), - strip.background=element_rect(fill="white"), - axis.text = element_text(size=13), - title = element_text(size=13), - legend.text = element_text(size=13))+ - scale_x_continuous(breaks=years_plot)+ - scale_fill_brewer(palette = "Set1")+ - xlab("Year")+ - ylab("Final energy demand [EJ]")+ - facet_wrap(~EDGE_scenario) + - guides(fill=guide_legend(title="Technology")) - - plot_list = list(plot_psm, plot_psm_LDV) - return(plot_list) -} +# LDVs final energy demand +```{r, echo=FALSE, warning=FALSE} +## demand EJ for LDV, divided by fuel type + +demandEJLDVplotf <- 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")] + + p = ggplot()+ + geom_area(data = demandEJ[year > 2010], aes(x=year, y=demand_EJ, group = technology, fill = technology), color="black",position= position_stack())+ + labs(x = "", y = "Final energy demand for LDVs [EJ]")+ + facet_wrap(~region, nrow = 4) + theme_minimal()+ + # scale_fill_manual("Vehicle Type",values = cols, breaks=legend_ord)+ + theme(axis.text.x = element_text(size = 7), + axis.text.y = element_text(size=7), + 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)) + + return(p) +} -FE_modes_bar_oilcomponent(demandEJ = demandF_plot_EJ, msect="trn_pass", region=region_plot) -FE_modes_bar_oilcomponent(demandEJ = demandF_plot_EJ, msect="trn_pass", region="GLO") +demandEJLDVplotf(demand_ej) ``` +# Energy services demand + +```{r, echo=FALSE, warning=FALSE} +demandpkmplotf = function(demandpkm){ + ## REMIND-EDGE results + demandpkm <- demandpkm[,c("sector","subsector_L3","subsector_L2", + "subsector_L1","vehicle_type","technology", "iso","year","demand_F")] + demandpkm[,demand_F:=demand_F ## in millionkm + *1e-6 ## in trillion km + ] + ## attribute aggregated mode and vehicle names for plotting purposes, and aggregate + demandpkm[, aggr_mode := ifelse(subsector_L1 %in% c("Three-Wheeler", "trn_pass_road_LDV_4W"), "LDV", NA)] + 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("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(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 = 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", + "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 = merge(demandpkm, REMIND2ISO_MAPPING, by = "iso") + demandpkm = demandpkm[, .(demand_F = sum(demand_F)), by = c("region", "year", "vehicle_type_plot", "aggr_mode", "mode")] + + demandpkm = demandpkm[order(aggr_mode)] + + p = ggplot()+ + geom_area(data = demandpkm[mode =="pass"& year > 2010], aes(x=year, y=demand_F, group = interaction(vehicle_type_plot,aggr_mode), fill = vehicle_type_plot), color="black",position= position_stack())+ + labs(x = "", y = "Energy Services demand [trillion pkm]")+ + facet_wrap(~region, nrow = 4) + theme_minimal()+ + # scale_fill_manual("Vehicle Type",values = cols, breaks=legend_ord)+ + theme(axis.text.x = element_text(size = 7), + axis.text.y = element_text(size=7), + 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)) + + + return(p) -```{r, echo=FALSE} -FE_modes_bar1=function(demandEJ, msect="trn_pass", region_plot){ - if(msect == "trn_pass") - sector_display = "Passenger" - if(msect == "trn_freight") - sector_display = "Freight" - - ## group by subsector_L1 and summarise the demand - df=demandEJ[, .(demand_EJ=sum(demand_EJ)), - by = c("EDGE_scenario", "region", "year","subsector_L1","sector")] - df=df[order(year) & year %in% years_plot] - ## give proper names to the categories - df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1") - #plot - plot_p=ggplot()+ - geom_bar(data=df%>%filter(sector == msect, region == region_plot), - aes(x=year,y=demand_EJ,group=mode,fill=mode), - position=position_stack(),stat="identity",color="black")+ - facet_wrap(~EDGE_scenario) + - ylab("Energy (EJ)") + - ggtitle(paste0("Final Energy Demand by Mode, ", sector_display, ", ", region_plot))+ - theme(axis.text.x = element_text(angle = 90))+ - scale_x_continuous(breaks=years_plot) - - return(plot_p) } -FE_modes_bar1(demandEJ = demandF_plot_EJ, region_plot = region_plot) +## energy services demand +demandpkmplotf(demand_km) ``` +# CO2 intensity of new sales +```{r, echo=FALSE, warning=FALSE} -```{r, echo=FALSE} -FE_modes_bar_oilVSelec=function(demandEJ, region_plot){ - sector_display = "Total transport" - ## group by subsector_L1 and summarise the demand - df=demandEJ[, .(demand_EJ=sum(demand_EJ)), - by = c("EDGE_scenario", "region", "year","subsector_L1","sector","technology")] - df[, tech_plot := ifelse(technology %in% c("BEV","Electric"), "Electriciy", NA)] - df[, tech_plot := ifelse(technology %in% c("Liquids", "Hybrid Liquids"), "Liquids", tech_plot)] - df=df[!is.na(tech_plot),] ## only liquids and electric driven entries interesting - df=df[order(year) & year %in% years_plot] - ## give proper names to the categories - df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1") - df[,short_names:=ifelse(mode %in% c("Buses","Rail Passenger","High Speed Rail"),"Other Passenger",NA)] - df[,short_names:=ifelse(mode %in% c("2W","4W","Three Wheelers"),"LDV",short_names)] - df[,short_names:=ifelse(mode %in% c("Domestic Aviation","International Aviation"),"Aviation",short_names)] - df[,short_names:=ifelse(mode %in% c("International Shipping","Domestic Shipping"),"Shipping",short_names)] - df[,short_names:=ifelse(mode %in% c("Road Freight","Rail Freight"),"Road and Rail Freight",short_names)] - - #plot - plot_p=ggplot()+ - geom_bar(data=df%>%filter(region == region_plot, year ==2050), - aes(x=tech_plot,y=demand_EJ,group=short_names,fill=short_names), - position=position_stack(),stat="identity",alpha=0.95)+ - facet_wrap(~EDGE_scenario) + - ylab("Energy (EJ)") + - ggtitle(paste0("Final Energy Demand by Mode in 2050, total transport ", region_plot))+ - theme_light()+ - theme(axis.text.x = element_text(angle = 90), - axis.title.x = element_blank(), - strip.text.x = element_text(size = 13, color = "black"), - strip.background=element_rect(fill="white"), - axis.text = element_text(size=13), - title = element_text(size=13), - legend.text = element_text(size=13))+ - scale_fill_brewer(palette = "Set2")+ - guides(fill=guide_legend(title="Transport mode")) - - return(plot_p) -} +CO2km_intensity_newsalesplotf = function(annual_sales, mj_km_data, sharesVS1, shares_source_liquids){ + shares_source_liquids[, technology := ifelse(variable %in% c("FE|Transport|Liquids|Oil", "FE|Transport|Liquids|Coal"), "Oil", "Biodiesel")] + shares_source_liquids = shares_source_liquids[,.(value = sum(value)), by = c("model","scenario","region", "period", "unit","technology")] + shares_source_liquids = shares_source_liquids[region != "World"] + shares_source_liquids[, region:=as.character(region)] + shares_source_liquids[, year := period] + shares_source_liquids[, period:=NULL] + gdp <- getRMNDGDP(scenario = "SSP2", usecache = T) + shares_source_liquids <- disaggregate_dt(shares_source_liquids, REMIND2ISO_MAPPING, + valuecol="value", + datacols=c("model","scenario", "unit","technology"), + weights=gdp) -FE_modes_bar_oilVSelec(demandEJ = demandF_plot_EJ, region_plot = region_plot) -``` + shares_source_liquids[, shareliq := value/sum(value),by=c("iso", "year")] + # ## CO2 content + # CO2_petrol = 3.1 ## gCO2/gFUEL + # CO2_biodiesel = 2.7 ## TODO this number is made up! + # CO2_cng = 2.7 ## gCO2/gFUEL -## FE composition - -```{r, echo=FALSE} - -FE_modeshares_area=function(demandEJ){ - #group by subsector_L3 and summarise the demand - df=demandEJ[, .(demand_EJ=sum(demand_EJ)), - by = c("EDGE_scenario", "region", "year","subsector_L3")] - #order by year - df=df[order(year)] - df=df[year>=2005,] - #plot - plot=ggplot()+ - geom_area(data=df %>% filter(year <= max(years_plot), region == region_plot), - aes(x=year,y=demand_EJ,group=subsector_L3,fill=subsector_L3), - color="black")+ - facet_wrap(~EDGE_scenario)+ - ylab("Energy (EJ)") + - ggtitle("Final Energy Demand, Mode composition")+ - theme(axis.text.x = element_text(angle = 90, hjust = 1)) - - return(plot) -} + ## TODO of CO2 content of biodiesel is made up! gCO2/gFUEL + emi_fuel = data.table(technology = c("Oil", "Biodiesel", "NG"), ei_gF_MJ = c(20, 20, 20), emi_cGO2_gF = c(3.1, 3.1, 2.7)) -p=FE_modeshares_area(demandEJ = demandF_plot_EJ) -p -## ggsave("FE_modeshares.png") -``` + emi_liquids = merge(shares_source_liquids, emi_fuel, all.x = TRUE, by = "technology") + emi_liquids = emi_liquids[, .(ei_gF_MJ = sum(shareliq*ei_gF_MJ), emi_cGO2_gF = sum(shareliq*emi_cGO2_gF)), by = c("iso", "year")][, technology := "Liquids"] + emi_NG = cbind(emi_fuel[technology == "NG"], unique(shares_source_liquids[,c("year", "iso")])) + emi_fuel = rbind(emi_NG, emi_liquids) + emi_fuel[, gCO2_MJ := ei_gF_MJ*emi_cGO2_gF] + emi_fuel = merge(mj_km_data[subsector_L1 == "trn_pass_road_LDV_4W"], emi_fuel, all.x = TRUE, by = c("iso", "year", "technology")) + emi_fuel[is.na(gCO2_MJ) & !technology %in% c("Liquids", "NG"), gCO2_MJ := 0] + emi_fuel[, gCO2_km := MJ_km * gCO2_MJ] -```{r, echo=FALSE} -## fuel use by sector -fuel_shares_area=function(demandEJ, msect="trn_pass", region_plot){ - if(msect == "trn_pass") - sector_display = "Passenger" - if(msect == "trn_freight") - sector_display = "Freight" - ##group by sector and technology and summarise demand - df=demandEJ[, .(demand_EJ=sum(demand_EJ)), - by = c("EDGE_scenario", "region", "year","technology","sector")] - - df=df[order(year) & year>=2005,] - #plot - plot1=ggplot()+ - geom_area(data=df%>%filter(sector == msect, year <= max(years_plot), region == region_plot), - aes(x=year,y=demand_EJ,group=technology,fill=technology),position="fill")+ - facet_wrap(~EDGE_scenario)+ - ggtitle(paste0("Final Energy Demand, ", sector_display, ", Fuel Composition"))+ - ylab("Share") - theme(axis.text.x = element_text(angle = 90, hjust = 1)) - - return(plot1) -} + totalemi = merge(emi_fuel, annual_sales, all.y = TRUE, by = c("iso", "year", "technology", "vehicle_type", "subsector_L1"), all.x = TRUE) + totalemi = totalemi[!is.na(share) & !is.na(gCO2_km)] + totalemi[, gCO2_km_ave := gCO2_km*share] -p=fuel_shares_area(demandEJ = demandF_plot_EJ, region_plot = region_plot) -p -``` + ##totalemi = merge(totalemi, demand_ej_plot) -```{r, echo=FALSE} + totalemi = totalemi[,.(gCO2_km_ave = sum(gCO2_km_ave)), by = c("year", "iso", "vehicle_type")] -SW_trend_plot = function(FV_SW,sector_plot){ - if (sector_plot == "trn_pass") { - FV_SW=FV_SW[iso=="DEU" & vehicle_type =="Large Car and SUV",] - } else if (sector_plot =="trn_freight"){ - FV_SW=FV_SW[iso=="DEU" & vehicle_type =="Truck (16-32t)",] - } else if (sector_plot =="trn_aviation_intl"){ - FV_SW=FV_SW[iso=="DEU" & subsector_L3 =="International Aviation",] - } else if (sector_plot =="trn_shipping_intl"){ - FV_SW=FV_SW[iso=="DEU" & subsector_L3 =="International Ship",] - } - FV_SW[,type:=ifelse(technology=="Liquids", "Conventional ICE (Liquid fuels)",NA)] - FV_SW[,type:=ifelse(technology=="NG", "Natural Gas ICE",type)] - FV_SW[,type:=ifelse(technology=="BEV", "Alternative fuels: BEV",type)] - FV_SW[,type:=ifelse(technology=="FCEV", "Alternative fuels: FCEV",type)] - FV_SW[,type:=ifelse(technology=="Hybrid Liquids", "Unconventional ICE (Hybrid)",type)] - - p=ggplot()+ - geom_line(data=FV_SW%>%filter(year>= min(years_plot), year<=max(years_plot)),aes(x=year,y=sw,group=type,color=type),alpha = 0.8,size=1.5)+ - theme_light()+ - facet_wrap(~EDGE_scenario)+ - theme(axis.text.x = element_text(angle = 90, hjust = 1), - axis.text = element_text(size=13), - title = element_text(size=14), - legend.text = element_text(size=13))+ - scale_x_continuous(breaks=years_plot)+ - xlab("Year")+ - ylab ("Preference factors tech. types [-]")+ - ggtitle(paste0("Preference factors trend for tech. types for ", sector_plot, " [-]"))+ - theme(strip.text.x = element_text(size=13,color="black"), - strip.background = element_rect(fill="white",color = "black"))+ - scale_color_discrete(name="Technology type") - return(p) - } + totalemi = merge(totalemi, sharesVS1, all.x = TRUE, by = c("iso", "year", "vehicle_type")) -p=SW_trend_plot(FV_SW=sw_tech,sector_plot) -p -``` + totalemi = totalemi[,.(gCO2_km_ave = sum(gCO2_km_ave*share)), by = c("iso", "year", "subsector_L1")] + + + totalemi = merge(totalemi, REMIND2ISO_MAPPING, by="iso") + totalemi = merge(totalemi, gdp, all.x=TRUE, by = c("iso", "year")) + totalemi[, share := weight/sum(weight), by = c("year", "region")] + totalemi = totalemi[,.(gCO2_km_ave = sum(gCO2_km_ave*share)), by = c("year", "region")] + p = ggplot()+ + geom_line(data = totalemi, aes(x = year, y = gCO2_km_ave))+ + labs(title = "gCO2/km average", y = "Average gCO2/km LDVs new additions")+ + facet_wrap(~region, nrow = 4)+ + theme_minimal() + + return(p) +} +shares_source_liquids = miffile[variable %in% c("FE|Transport|Liquids|Biomass", "FE|Transport|Liquids|Coal", "FE|Transport|Liquids|Oil"),] +CO2km_intensity_newsalesplotf(annual_sales, mj_km_data, sharesVS1 = shares$VS1_shares, shares_source_liquids) +``` \ No newline at end of file