Merge pull request #89 from martavp/distribute_CO2budget

Distribute CO2 budget among the planning horizons for the myopic option

Snakefile

@@ -8,7 +8,7 @@ wildcard_constraints:
     clusters="[0-9]+m?",
     sectors="[+a-zA-Z0-9]+",
     opts="[-+a-zA-Z0-9]*",
-    sector_opts="[-+a-zA-Z0-9]*"
+    sector_opts="[-+a-zA-Z0-9\.]*"
 
 
 
@@ -292,6 +292,7 @@ rule build_retro_cost:
     output:
         retro_cost="resources/retro_cost_{network}_s{simpl}_{clusters}.csv",
         floor_area="resources/floor_area_{network}_s{simpl}_{clusters}.csv"
+    resources: mem_mb=1000
     script: "scripts/build_retro_cost.py"
 
 

config.default.yaml

@@ -24,11 +24,17 @@ scenario:
   # B for biomass supply, I for industry, shipping and aviation
   # solarx or onwindx changes the available installable potential by factor x
   # dist{n} includes distribution grids with investment cost of n times cost in data/costs.csv
+  # for myopic/perfect foresight cb states the carbon budget in GtCO2 (cumulative 
+  # emissions throughout the transition path in the timeframe determined by the 
+  # planning_horizons), be:beta decay; ex:exponential decay
+  # cb40ex0 distributes a carbon budget of 40 GtCO2 following an exponential 
+  # decay with initial growth rate 0
   planning_horizons : [2030] # investment years for myopic and perfect; or costs year for overnight
   # for example, set to [2020, 2030, 2040, 2050] for myopic foresight
 
 # CO2 budget as a fraction of 1990 emissions
 # this is over-ridden if CO2Lx is set in sector_opts
+# this is also over-ridden if cb is set in sector_opts
 co2_budget:
   2020: 0.7011648746
   2025: 0.5241935484

doc/installation.rst

@@ -16,7 +16,7 @@ its dependencies. Clone the repository:
 
 .. code:: bash
 
-    projects % git clone git@github.com:PyPSA/pypsa-eur.git
+    projects % git clone https://github.com/PyPSA/pypsa-eur.git
 
 then download and unpack all the PyPSA-Eur data files by running the following snakemake rule:
 
@@ -32,7 +32,7 @@ Next install the technology assumptions database `technology-data <https://githu
 
 .. code:: bash
 
-    projects % git clone git@github.com:PyPSA/technology-data.git
+    projects % git clone https://github.com/PyPSA/technology-data.git
 
 
 Clone PyPSA-Eur-Sec repository
@@ -42,7 +42,7 @@ Create a parallel directory for `PyPSA-Eur-Sec <https://github.com/PyPSA/pypsa-e
 
 .. code:: bash
 
-    projects % git clone git@github.com:PyPSA/pypsa-eur-sec.git
+    projects % git clone https://github.com/PyPSA/pypsa-eur-sec.git
 
 Environment/package requirements
 ================================
@@ -54,6 +54,13 @@ The requirements are the same as `PyPSA-Eur <https://github.com/PyPSA/pypsa-eur>
 xarray version >= 0.15.1, you will need the latest master branch of
 atlite version 0.0.2.
 
+You can create an enviroment using the environment.yaml file in pypsa-eur/envs:
+
+.../pypsa-eur % conda env create -f envs/environment.yaml
+
+.../pypsa-eur % conda activate pypsa-eur
+
+See details in `PyPSA-Eur Installation <https://pypsa-eur.readthedocs.io/en/latest/installation.html>`_
 
 Data requirements
 =================

doc/myopic.rst

@@ -6,7 +6,7 @@ Myopic transition path
 
 The myopic code can be used to investigate progressive changes in a network, for instance, those taking place throughout a transition path. The capacities installed in a certain time step are maintained in the network until their operational lifetime expires.
 
-The myopic approach was initially developed and used in the paper `Early decarbonisation of the European Energy system pays off (2020) <https://arxiv.org/abs/2004.11009>`__ but the current implementation complies with the pypsa-eur-sec standard working flow and is compatible with using the higher resolution electricity transmission model `PyPSA-Eur <https://github.com/PyPSA/pypsa-eur>`__ rather than a one-node-per-country model.
+The myopic approach was initially developed and used in the paper `Early decarbonisation of the European Energy system pays off (2020) <https://www.nature.com/articles/s41467-020-20015-4>`__ but the current implementation complies with the pypsa-eur-sec standard working flow and is compatible with using the higher resolution electricity transmission model `PyPSA-Eur <https://github.com/PyPSA/pypsa-eur>`__ rather than a one-node-per-country model.
 
 The current code applies the myopic approach to generators, storage technologies and links in the power sector and the space and water heating sector.
 
@@ -61,12 +61,15 @@ Wildcards
 The {planning_horizons} wildcard indicates the timesteps in which the network is optimized, e.g. planning_horizons: [2020, 2030, 2040, 2050]
 
 
+Options
+=============
+The total carbon budget for the entire transition path can be indicated in the ``scenario.sector_opts`` in ``config.yaml``.
+The carbon budget can be split among the ``planning_horizons`` following an exponential or beta decay. 
+E.g. ``'cb40ex0'`` splits the a carbon budget equal to 40 GtCO_2 following an exponential decay whose initial linear growth rate $r$ is zero
 
-**{co2_budget_name} wildcard**
+$e(t) = e_0 (1+ (r+m)t) e^(-mt)$
 
-The {co2_budget_name} wildcard indicates the name of the co2 budget.
+See details in Supplementary Note 1 of the paper `Early decarbonisation of the European Energy system pays off (2020) <https://www.nature.com/articles/s41467-020-20015-4>`__
-
-A csv file is used as input including the planning_horizons as index, the name of co2_budget as column name, and the maximum co2 emissions (relative to 1990) as values.
 
 Rules overview
 =================
@@ -74,17 +77,17 @@ Rules overview
 General myopic code structure
 ===============================
 
-The myopic code solves the network for the time steps included in planning_horizons in a recursive loop, so that:
+The myopic code solves the network for the time steps included in ``planning_horizons`` in a recursive loop, so that:
 
 1.The existing capacities (those installed before the base year are added as fixed capacities with p_nom=value, p_nom_extendable=False). E.g. for baseyear=2020, capacities installed before 2020 are added. In addition, the network comprises additional generator, storage, and link capacities with p_nom_extendable=True. The non-solved network is saved in ``results/run_name/networks/prenetworks-brownfield``.
 
-The base year is the first element in planning_horizons. Step 1 is implemented with the rule add_baseyear for the base year and with the rule add_brownfield for the remaining planning_horizons.
+The base year is the first element in ``planning_horizons``. Step 1 is implemented with the rule add_baseyear for the base year and with the rule add_brownfield for the remaining planning_horizons.
 
-2.The 2020 network is optimized. The solved network is saved in ‘results/run_name/networks/postnetworks’
+2.The 2020 network is optimized. The solved network is saved in ``results/run_name/networks/postnetworks``
 
 3.For the next planning horizon, e.g. 2030, the capacities from a previous time step are added if they are still in operation (i.e., if they fulfil planning horizon <= commissioned year + lifetime). In addition, the network comprises additional generator, storage, and link capacities with p_nom_extendable=True. The non-solved network is saved in ``results/run_name/networks/prenetworks-brownfield``.
 
-Steps 2 and 3 are solved recursively for all the planning_horizons included in the configuration file.
+Steps 2 and 3 are solved recursively for all the planning_horizons included in ``config.yaml``.
 
 
 rule add_existing baseyear
@@ -110,8 +113,8 @@ Then, the resulting network is saved in ``results/run_name/networks/prenetworks-
 rule add_brownfield
 ===================
 
-The rule add_brownfield loads the network in ‘results/run_name/networks/prenetworks’ and performs the following operation:
+The rule add_brownfield loads the network in ``results/run_name/networks/prenetworks`` and performs the following operation:
 
-1.Read the capacities optimized in the previous time step and add them to the network if they are still in operation (i.e., if they fulfil planning horizon < commissioned year + lifetime)
+1.Read the capacities optimized in the previous time step and add them to the network if they are still in operation (i.e., if they fulfill planning horizon < commissioned year + lifetime)
 
 Then, the resulting network is saved in ``results/run_name/networks/prenetworks_brownfield``.

doc/release_notes.rst

@@ -2,6 +2,9 @@
 Release Notes
 ##########################################
 
+Future release
+===================
+*For the myopic option, a carbon budget and a type of decay (exponential or beta) can be selected in the config file to distribute the budget across the planning_horizons.
 
 PyPSA-Eur-Sec 0.4.0 (11th December 2020)
 =========================================

doc/supply_demand.rst

@@ -106,7 +106,7 @@ Thermal energy storage using hot water tanks
 
 Small for decentral applications.
 
-Big pit storage for district heating.
+Big water pit storage for district heating.
 
 
 Hydrogen demand
@@ -122,7 +122,7 @@ Industry (ammonia, precursor to hydrocarbons for chemicals and iron/steel).
 Hydrogen supply
 =================
 
-SMR, SMR+CCS, electrolysers.
+Steam Methane Reforming (SMR), SMR+CCS, electrolysers.
 
 
 Methane demand

scripts/build_energy_totals.py

@@ -390,12 +390,12 @@ def build_energy_totals():
     return clean_df
 
 
-def build_eea_co2():
+def build_eea_co2(year=1990):
     # see ../notebooks/compute_1990_Europe_emissions_for_targets.ipynb
 
-    #https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-14
+    #https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-16
-    #downloaded 190222 (modified by EEA last on 181130)
+    #downloaded 201228 (modified by EEA last on 201221)
-    fn = "data/eea/UNFCCC_v21.csv"
+    fn = "data/eea/UNFCCC_v23.csv"
     df = pd.read_csv(fn, encoding="latin-1")
     df.loc[df["Year"] == "1985-1987","Year"] = 1986
     df["Year"] = df["Year"].astype(int)
@@ -418,16 +418,14 @@ def build_eea_co2():
     e['waste management'] = '5 - Waste management'
     e['other'] = '6 - Other Sector'
     e['indirect'] = 'ind_CO2 - Indirect CO2'
-    e["total wL"] = "Total (with LULUCF, with indirect CO2)"
+    e["total wL"] = "Total (with LULUCF)"
-    e["total woL"] = "Total (without LULUCF, with indirect CO2)"
+    e["total woL"] = "Total (without LULUCF)"
 
 
     pol = "CO2" #["All greenhouse gases - (CO2 equivalent)","CO2"]
 
     cts = ["CH","EUA","NO"] + eu28_eea
 
-    year = 1990
-
     emissions = df.loc[idx[cts,pol,year,e.values],"emissions"].unstack("Sector_name").rename(columns=pd.Series(e.index,e.values)).rename(index={"All greenhouse gases - (CO2 equivalent)" : "GHG"},level=1)
 
     #only take level 0, since level 1 (pol) and level 2 (year) are trivial
@@ -467,7 +465,7 @@ def build_eurostat_co2(year=1990):
     return eurostat_co2
 
 
-def build_co2_totals(year=1990):
+def build_co2_totals(eea_co2, eurostat_co2, year=1990):
 
     co2 = eea_co2.reindex(["EU28","NO","CH","BA","RS","AL","ME","MK"] + eu28)
 
@@ -486,10 +484,6 @@ def build_co2_totals(year=1990):
         #doesn't include non-energy emissions
         co2.loc[ct,'agriculture'] = eurostat_co2[ct,"+","+","Agriculture / Forestry"].sum()
 
-
-
-    co2.to_csv(snakemake.output.co2_name)
-
     return co2
 
 
@@ -547,7 +541,7 @@ if __name__ == "__main__":
         snakemake.output['transport_name'] = "data/transport_data.csv"
 
         snakemake.input = Dict()
-        snakemake.input['nuts3_shapes'] = 'resources/nuts3_shapes.geojson'
+        snakemake.input['nuts3_shapes'] = '../pypsa-eur/resources/nuts3_shapes.geojson'
 
     nuts3 = gpd.read_file(snakemake.input.nuts3_shapes).set_index('index')
     population = nuts3['pop'].groupby(nuts3.country).sum()
@@ -566,6 +560,7 @@ if __name__ == "__main__":
 
     eurostat_co2 = build_eurostat_co2()
 
-    build_co2_totals()
+    co2=build_co2_totals(eea_co2, eurostat_co2, year)
+    co2.to_csv(snakemake.output.co2_name)
     
     build_transport_data()

scripts/make_summary.py

@@ -196,6 +196,24 @@ def calculate_costs(n,label,costs):
 
     return costs
 
+def calculate_cumulative_cost():        
+    planning_horizons = snakemake.config['scenario']['planning_horizons']
+
+    cumulative_cost = pd.DataFrame(index = df["costs"].sum().index,
+                                  columns=pd.Series(data=np.arange(0,0.1, 0.01), name='social discount rate'))
+    
+    #discount cost and express them in money value of planning_horizons[0]
+    for r in cumulative_cost.columns:
+        cumulative_cost[r]=[df["costs"].sum()[index]/((1+r)**(index[-1]-planning_horizons[0])) for index in cumulative_cost.index]
+    
+    #integrate cost throughout the transition path
+    for r in cumulative_cost.columns:       
+        for cluster in cumulative_cost.index.get_level_values(level=0).unique():
+            for lv in cumulative_cost.index.get_level_values(level=1).unique():
+                for sector_opts in cumulative_cost.index.get_level_values(level=2).unique():
+                    cumulative_cost.loc[(cluster, lv, sector_opts,'cumulative cost'),r] = np.trapz(cumulative_cost.loc[idx[cluster, lv, sector_opts,planning_horizons],r].values, x=planning_horizons)
+
+    return cumulative_cost 
     
 def calculate_nodal_capacities(n,label,nodal_capacities):
     #Beware this also has extraneous locations for country (e.g. biomass) or continent-wide (e.g. fossil gas/oil) stuff
@@ -564,16 +582,17 @@ if __name__ == "__main__":
             snakemake.config = yaml.safe_load(f)
 
         #overwrite some options
-        snakemake.config["run"] = "test"
+        snakemake.config["run"] = "version-8"
         snakemake.config["scenario"]["lv"] = [1.0]
-        snakemake.config["scenario"]["sector_opts"] = ["Co2L0-168H-T-H-B-I-solar3-dist1"]
+        snakemake.config["scenario"]["sector_opts"] = ["3H-T-H-B-I-solar3-dist1"]
         snakemake.config["planning_horizons"] = ['2020', '2030', '2040', '2050']
         snakemake.input = Dict()
         snakemake.input['heat_demand_name'] = 'data/heating/daily_heat_demand.h5'
+        snakemake.input['costs'] = snakemake.config['costs_dir'] + "costs_{}.csv".format(snakemake.config['scenario']['planning_horizons'][0])
         snakemake.output = Dict()
         for item in outputs:
             snakemake.output[item] = snakemake.config['summary_dir'] + '/{name}/csvs/{item}.csv'.format(name=snakemake.config['run'],item=item)
-
+        snakemake.output['cumulative_cost'] = snakemake.config['summary_dir'] + '/{name}/csvs/cumulative_cost.csv'.format(name=snakemake.config['run'])
     networks_dict = {(cluster, lv, opt+sector_opt, planning_horizon) :
                      snakemake.config['results_dir'] + snakemake.config['run'] + '/postnetworks/elec_s{simpl}_{cluster}_lv{lv}_{opt}_{sector_opt}_{planning_horizon}.nc'\
                      .format(simpl=simpl,
@@ -592,6 +611,7 @@ if __name__ == "__main__":
     print(networks_dict)
 
     Nyears = 1
+    
     costs_db = prepare_costs(snakemake.input.costs,
                              snakemake.config['costs']['USD2013_to_EUR2013'],
                              snakemake.config['costs']['discountrate'],
@@ -603,3 +623,10 @@ if __name__ == "__main__":
     df["metrics"].loc["total costs"] =  df["costs"].sum()
 
     to_csv(df)
+    
+    if snakemake.config["foresight"]=='myopic':
+        cumulative_cost=calculate_cumulative_cost()
+        cumulative_cost.to_csv(snakemake.config['summary_dir'] + '/' + snakemake.config['run'] + '/csvs/cumulative_cost.csv')
+
+
+                                                     

scripts/plot_summary.py

@@ -1,6 +1,6 @@
 
 
-
+import numpy as np
 import pandas as pd
 
 #allow plotting without Xwindows
@@ -9,7 +9,7 @@ matplotlib.use('Agg')
 
 import matplotlib.pyplot as plt
 
-
+from prepare_sector_network import co2_emissions_year
 
 #consolidate and rename
 def rename_techs(label):
@@ -237,7 +237,137 @@ def plot_balances():
 
         fig.savefig(snakemake.output.balances[:-10] + k + ".pdf",transparent=True)
 
+def historical_emissions(cts):
+    """
+    read historical emissions to add them to the carbon budget plot
+    """
+    #https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-16
+    #downloaded 201228 (modified by EEA last on 201221)
+    fn = "data/eea/UNFCCC_v23.csv"
+    df = pd.read_csv(fn, encoding="latin-1")
+    df.loc[df["Year"] == "1985-1987","Year"] = 1986
+    df["Year"] = df["Year"].astype(int)
+    df = df.set_index(['Year', 'Sector_name', 'Country_code', 'Pollutant_name']).sort_index()
 
+    e = pd.Series()
+    e["electricity"] = '1.A.1.a - Public Electricity and Heat Production'
+    e['residential non-elec'] = '1.A.4.b - Residential'
+    e['services non-elec'] = '1.A.4.a - Commercial/Institutional'
+    e['rail non-elec'] = "1.A.3.c - Railways"
+    e["road non-elec"] = '1.A.3.b - Road Transportation'
+    e["domestic navigation"] = "1.A.3.d - Domestic Navigation"
+    e['international navigation'] = '1.D.1.b - International Navigation'
+    e["domestic aviation"] = '1.A.3.a - Domestic Aviation'
+    e["international aviation"] = '1.D.1.a - International Aviation'   
+    e['total energy'] = '1 - Energy'
+    e['industrial processes'] = '2 - Industrial Processes and Product Use'
+    e['agriculture'] = '3 - Agriculture'
+    e['LULUCF'] = '4 - Land Use, Land-Use Change and Forestry'
+    e['waste management'] = '5 - Waste management'
+    e['other'] = '6 - Other Sector'
+    e['indirect'] = 'ind_CO2 - Indirect CO2'
+    e["total wL"] = "Total (with LULUCF)"
+    e["total woL"] = "Total (without LULUCF)"
+       
+    pol = ["CO2"] # ["All greenhouse gases - (CO2 equivalent)"] 
+    cts
+    if "GB" in cts:
+        cts.remove("GB")
+        cts.append("UK")
+         
+    year = np.arange(1990,2018).tolist()
+
+    idx = pd.IndexSlice
+    co2_totals = df.loc[idx[year,e.values,cts,pol],"emissions"].unstack("Year").rename(index=pd.Series(e.index,e.values)) 
+    
+    co2_totals = (1/1e6)*co2_totals.groupby(level=0, axis=0).sum() #Gton CO2
+
+    co2_totals.loc['industrial non-elec'] = co2_totals.loc['total energy'] - co2_totals.loc[['electricity', 'services non-elec','residential non-elec', 'road non-elec',
+                                                                              'rail non-elec', 'domestic aviation', 'international aviation', 'domestic navigation',
+                                                                              'international navigation']].sum()
+
+    emissions = co2_totals.loc["electricity"]   
+    if "T" in opts:
+        emissions += co2_totals.loc[[i+ " non-elec" for i in ["rail","road"]]].sum()
+    if "H" in opts:
+        emissions += co2_totals.loc[[i+ " non-elec" for i in ["residential","services"]]].sum()
+    if "I" in opts:
+        emissions += co2_totals.loc[["industrial non-elec","industrial processes",
+                                          "domestic aviation","international aviation",
+                                          "domestic navigation","international navigation"]].sum()           
+    return emissions
+
+
+
+def plot_carbon_budget_distribution():
+    """
+    Plot historical carbon emissions in the EU and decarbonization path
+    """ 
+    
+    import matplotlib.gridspec as gridspec
+    import seaborn as sns; sns.set()
+    sns.set_style('ticks')
+    plt.style.use('seaborn-ticks')
+    plt.rcParams['xtick.direction'] = 'in'
+    plt.rcParams['ytick.direction'] = 'in'
+    plt.rcParams['xtick.labelsize'] = 20
+    plt.rcParams['ytick.labelsize'] = 20   
+
+    plt.figure(figsize=(10, 7))
+    gs1 = gridspec.GridSpec(1, 1)
+    ax1 = plt.subplot(gs1[0,0])
+    ax1.set_ylabel('CO$_2$ emissions (Gt per year)',fontsize=22)
+    ax1.set_ylim([0,5])
+    ax1.set_xlim([1990,snakemake.config['scenario']['planning_horizons'][-1]+1])
+    
+    path_cb = snakemake.config['results_dir'] + snakemake.config['run'] + '/csvs/'
+    countries=pd.read_csv(path_cb + 'countries.csv',  index_col=1) 
+    cts=countries.index.to_list()
+    e_1990 = co2_emissions_year(cts, opts, year=1990)     
+    CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv',  
+                        index_col=0) 
+    
+    
+    ax1.plot(e_1990*CO2_CAP[o],linewidth=3, 
+                         color='dodgerblue', label=None)
+            
+    emissions = historical_emissions(cts)
+
+    ax1.plot(emissions, color='black', linewidth=3, label=None) 
+    
+    #plot commited and uder-discussion targets  
+    #(notice that historical emissions include all countries in the
+    # network, but targets refer to EU)
+    ax1.plot([2020],[0.8*emissions[1990]],
+                     marker='*', markersize=12, markerfacecolor='black',
+                     markeredgecolor='black')    
+            
+    ax1.plot([2030],[0.45*emissions[1990]],
+                     marker='*', markersize=12, markerfacecolor='white',
+                     markeredgecolor='black')    
+            
+    ax1.plot([2030],[0.6*emissions[1990]],
+                     marker='*', markersize=12, markerfacecolor='black',
+                     markeredgecolor='black')
+            
+    ax1.plot([2050, 2050],[x*emissions[1990] for x in [0.2, 0.05]],
+                  color='gray', linewidth=2, marker='_', alpha=0.5) 
+            
+    ax1.plot([2050],[0.01*emissions[1990]],
+                     marker='*', markersize=12, markerfacecolor='white', 
+                     linewidth=0, markeredgecolor='black', 
+                     label='EU under-discussion target', zorder=10, 
+                     clip_on=False) 
+            
+    ax1.plot([2050],[0.125*emissions[1990]],'ro',
+                     marker='*', markersize=12, markerfacecolor='black',
+                     markeredgecolor='black', label='EU commited target')
+            
+    ax1.legend(fancybox=True, fontsize=18, loc=(0.01,0.01), 
+                       facecolor='white', frameon=True) 
+            
+    path_cb_plot = snakemake.config['results_dir'] + snakemake.config['run'] + '/graphs/'             
+    plt.savefig(path_cb_plot+'carbon_budget_plot.pdf', dpi=300) 
 
 if __name__ == "__main__":
     # Detect running outside of snakemake and mock snakemake for testing
@@ -249,12 +379,15 @@ if __name__ == "__main__":
             snakemake.config = yaml.safe_load(f)
         snakemake.input = Dict()
         snakemake.output = Dict()
+        snakemake.wildcards = Dict()
+        #snakemake.wildcards['sector_opts']='3H-T-H-B-I-solar3-dist1-cb48be3'        
         
         for item in ["costs", "energy"]:
             snakemake.input[item] = snakemake.config['summary_dir'] + '/{name}/csvs/{item}.csv'.format(name=snakemake.config['run'],item=item)
             snakemake.output[item] = snakemake.config['summary_dir'] + '/{name}/graphs/{item}.pdf'.format(name=snakemake.config['run'],item=item)
-        snakemake.input["balances"] = snakemake.config['summary_dir'] + '/test/csvs/supply_energy.csv'
+        snakemake.input["balances"] = snakemake.config['summary_dir'] + '/{name}/csvs/supply_energy.csv'.format(name=snakemake.config['run'],item=item)
-        snakemake.output["balances"] = snakemake.config['summary_dir'] + '/test/graphs/balances-energy.csv'
+        snakemake.output["balances"] = snakemake.config['summary_dir'] + '/{name}/graphs/balances-energy.csv'.format(name=snakemake.config['run'],item=item)
+        
         
     n_header = 4
 
@@ -263,3 +396,9 @@ if __name__ == "__main__":
     plot_energy()
 
     plot_balances()
+    
+    for sector_opts in snakemake.config['scenario']['sector_opts']:
+        opts=sector_opts.split('-')
+        for o in opts:
+            if "cb" in o:
+                plot_carbon_budget_distribution()

scripts/prepare_sector_network.py

@@ -20,6 +20,8 @@ import pytz
 
 from vresutils.costdata import annuity
 
+from scipy.stats import beta
+from build_energy_totals import build_eea_co2, build_eurostat_co2, build_co2_totals
 
 #First tell PyPSA that links can have multiple outputs by
 #overriding the component_attrs. This can be done for
@@ -45,6 +47,83 @@ override_component_attrs["Generator"].loc["lifetime"] = ["float","years",np.nan,
 override_component_attrs["Store"].loc["build_year"] = ["integer","year",np.nan,"build year","Input (optional)"]
 override_component_attrs["Store"].loc["lifetime"] = ["float","years",np.nan,"lifetime","Input (optional)"]
 
+
+def co2_emissions_year(cts, opts, year):
+    """
+    calculate co2 emissions in one specific year (e.g. 1990 or 2018).
+    """
+    eea_co2 = build_eea_co2(year)
+    
+    #TODO: read Eurostat data from year>2014, this only affects the estimation of 
+    # CO2 emissions for "BA","RS","AL","ME","MK"
+    if year > 2014:
+        eurostat_co2 = build_eurostat_co2(year=2014)
+    else:
+        eurostat_co2 = build_eurostat_co2(year)
+        
+    co2_totals=build_co2_totals(eea_co2, eurostat_co2, year)
+
+    co2_emissions = co2_totals.loc[cts, "electricity"].sum()
+
+    if "T" in opts:
+        co2_emissions += co2_totals.loc[cts, [i+ " non-elec" for i in ["rail","road"]]].sum().sum()
+    if "H" in opts:
+        co2_emissions += co2_totals.loc[cts, [i+ " non-elec" for i in ["residential","services"]]].sum().sum()
+    if "I" in opts:
+        co2_emissions += co2_totals.loc[cts, ["industrial non-elec","industrial processes",
+                                              "domestic aviation","international aviation",
+                                              "domestic navigation","international navigation"]].sum().sum()
+    co2_emissions *=0.001 #MtCO2 to GtCO2
+    return co2_emissions
+
+
+def build_carbon_budget(o):
+    #distribute carbon budget following beta or exponential transition path
+    if "be" in o:    
+        #beta decay
+        carbon_budget = float(o[o.find("cb")+2:o.find("be")])
+        be=float(o[o.find("be")+2:])        
+    if "ex" in o:   
+        #exponential decay
+        carbon_budget = float(o[o.find("cb")+2:o.find("ex")])
+        r=float(o[o.find("ex")+2:])
+    
+
+    pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
+    pop_layout["ct"] = pop_layout.index.str[:2]
+    cts = pop_layout.ct.value_counts().index          
+    
+    e_1990 = co2_emissions_year(cts, opts, year=1990) 
+            
+    #emissions at the beginning of the path (last year available 2018)
+    e_0 = co2_emissions_year(cts, opts, year=2018) 
+    #emissions in 2019 and 2020 assumed equal to 2018 and substracted
+    carbon_budget -= 2*e_0
+    planning_horizons = snakemake.config['scenario']['planning_horizons']
+    CO2_CAP = pd.DataFrame(index = pd.Series(data=planning_horizons, 
+                                             name='planning_horizon'),
+                                             columns=pd.Series(data=[],
+                                                              name='paths', 
+                                                              dtype='float'))
+    t_0 = planning_horizons[0]
+    if "be" in o: 
+        #beta decay                
+        t_f = t_0 + (2*carbon_budget/e_0).round(0) # final year in the path 
+        #emissions (relative to 1990)
+        CO2_CAP[o] = [(e_0/e_1990)*(1-beta.cdf((t-t_0)/(t_f-t_0), be, be)) for t in planning_horizons]
+            
+    if "ex" in o:   
+        #exponential decay without delay
+        T=carbon_budget/e_0
+        m=(1+np.sqrt(1+r*T))/T
+        CO2_CAP[o] = [(e_0/e_1990)*(1+(m+r)*(t-t_0))*np.exp(-m*(t-t_0)) for t in planning_horizons]
+                
+    CO2_CAP.to_csv(path_cb + 'carbon_budget_distribution.csv', sep=',', 
+                   line_terminator='\n',  float_format='%.3f') 
+    countries=pd.Series(data=cts)
+    countries.to_csv(path_cb + 'countries.csv', sep=',', 
+               line_terminator='\n',  float_format='%.3f')
+            
 def add_lifetime_wind_solar(n):
     """
     Add lifetime for solar and wind generators
@@ -1775,15 +1854,15 @@ def get_parameter(item):
         return item
 
 
-#%%
+
 if __name__ == "__main__":
     # Detect running outside of snakemake and mock snakemake for testing
     if 'snakemake' not in globals():
         from vresutils.snakemake import MockSnakemake
         snakemake = MockSnakemake(
             wildcards=dict(network='elec', simpl='', clusters='37', lv='1.0',
-                           opts='', planning_horizons='2030', co2_budget_name="go",
+                           opts='', planning_horizons='2020', 
-                           sector_opts='Co2L0-120H-T-H-B-I-solar3-dist1'),
+                           sector_opts='120H-T-H-B-I-solar3-dist1-cb48be3'),
             input=dict( network='../pypsa-eur/networks/{network}_s{simpl}_{clusters}_ec_lv{lv}_{opts}.nc',
                         energy_totals_name='resources/energy_totals.csv',
                         co2_totals_name='resources/co2_totals.csv',
@@ -1819,10 +1898,10 @@ if __name__ == "__main__":
                         solar_thermal_total="resources/solar_thermal_total_{network}_s{simpl}_{clusters}.nc",
                         solar_thermal_urban="resources/solar_thermal_urban_{network}_s{simpl}_{clusters}.nc",
                         solar_thermal_rural="resources/solar_thermal_rural_{network}_s{simpl}_{clusters}.nc",
-                	retro_cost_energy = "resources/retro_cost_{network}_s{simpl}_{clusters}.csv",
+                	    retro_cost_energy = "resources/retro_cost_{network}_s{simpl}_{clusters}.csv",
                         floor_area = "resources/floor_area_{network}_s{simpl}_{clusters}.csv"
             ),
-            output=['pypsa-eur-sec/results/test/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{co2_budget_name}_{planning_horizons}.nc']
+            output=['results/version-cb48be3/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc']
         )
         import yaml
         with open('config.yaml', encoding='utf8') as f:
@@ -1926,6 +2005,20 @@ if __name__ == "__main__":
     limit = get_parameter(snakemake.config["co2_budget"])
     print("CO2 limit set to",limit)
 
+    for o in opts:
+        if "cb" in o:    
+            path_cb = snakemake.config['results_dir'] + snakemake.config['run'] + '/csvs/'
+            if not os.path.exists(path_cb):
+                os.makedirs(path_cb)
+            try:
+                CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv', index_col=0)
+            except:    
+                build_carbon_budget(o)
+                CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv', index_col=0)
+            
+            limit=CO2_CAP.loc[investment_year]            
+            print("overriding CO2 limit with scenario limit",limit)
+            
     for o in opts:
         if "Co2L" in o:
             limit = o[o.find("Co2L")+4:]