Fix bug myopic co2

Snakefile

@@ -442,14 +442,14 @@ rule build_population_weighted_energy_totals:
 
 
 rule build_transport_demand:
-    input: 
+    input:
         clustered_pop_layout="resources/pop_layout_elec_s{simpl}_{clusters}.csv",
         pop_weighted_energy_totals="resources/pop_weighted_energy_totals_s{simpl}_{clusters}.csv",
         transport_data='resources/transport_data.csv',
         traffic_data_KFZ="data/emobility/KFZ__count",
         traffic_data_Pkw="data/emobility/Pkw__count",
         temp_air_total="resources/temp_air_total_elec_s{simpl}_{clusters}.nc",
-    output: 
+    output:
         transport_demand="resources/transport_demand_s{simpl}_{clusters}.csv",
         transport_data="resources/transport_data_s{simpl}_{clusters}.csv",
         avail_profile="resources/avail_profile_s{simpl}_{clusters}.csv",
@@ -464,12 +464,14 @@ rule prepare_sector_network:
         overrides="data/override_component_attrs",
         network=pypsaeur('networks/elec_s{simpl}_{clusters}_ec_lv{lv}_{opts}.nc'),
         energy_totals_name='resources/energy_totals.csv',
+        eurostat=input_eurostat,
         pop_weighted_energy_totals="resources/pop_weighted_energy_totals_s{simpl}_{clusters}.csv",
         transport_demand="resources/transport_demand_s{simpl}_{clusters}.csv",
         transport_data="resources/transport_data_s{simpl}_{clusters}.csv",
         avail_profile="resources/avail_profile_s{simpl}_{clusters}.csv",
         dsm_profile="resources/dsm_profile_s{simpl}_{clusters}.csv",
         co2_totals_name='resources/co2_totals.csv',
+        co2="data/eea/UNFCCC_v23.csv",
         biomass_potentials='resources/biomass_potentials_s{simpl}_{clusters}.csv',
         heat_profile="data/heat_load_profile_BDEW.csv",
         costs=CDIR + "costs_{planning_horizons}.csv",
@@ -568,7 +570,9 @@ rule plot_summary:
     input:
         costs=SDIR + '/csvs/costs.csv',
         energy=SDIR + '/csvs/energy.csv',
-        balances=SDIR + '/csvs/supply_energy.csv'
+        balances=SDIR + '/csvs/supply_energy.csv',
+        eurostat=input_eurostat,
+        country_codes='data/Country_codes.csv',
     output:
         costs=SDIR + '/graphs/costs.pdf',
         energy=SDIR + '/graphs/energy.pdf',

scripts/add_existing_baseyear.py

@@ -193,24 +193,33 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
 
     for grouping_year, generator in df.index:
 
+
         # capacity is the capacity in MW at each node for this
         capacity = df.loc[grouping_year, generator]
         capacity = capacity[~capacity.isna()]
         capacity = capacity[capacity > snakemake.config['existing_capacities']['threshold_capacity']]
-
+        suffix = '-ac' if generator == 'offwind' else ''
+        name_suffix = f' {generator}{suffix}-{baseyear}'
+        asset_i = capacity.index + name_suffix
         if generator in ['solar', 'onwind', 'offwind']:
 
-            suffix = '-ac' if generator == 'offwind' else ''
-            name_suffix = f' {generator}{suffix}-{baseyear}'
-
             # to consider electricity grid connection costs or a split between
             # solar utility and rooftop as well, rather take cost assumptions
             # from existing network than from the cost database
             capital_cost = n.generators.loc[n.generators.carrier==generator+suffix, "capital_cost"].mean()
 
+            # check if assets are already in network (e.g. for 2020)
+            already_build = n.generators.index.intersection(asset_i)
+            new_build = asset_i.difference(n.generators.index)
+
+            # this is for the year 2020
+            if not already_build.empty:
+                n.generators.loc[already_build, "p_nom_min"] = capacity.loc[already_build.str.replace(name_suffix, "")].values
+            new_capacity = capacity.loc[new_build.str.replace(name_suffix, "")]
+
             if 'm' in snakemake.wildcards.clusters:
 
-                for ind in capacity.index:
+                for ind in new_capacity.index:
 
                     # existing capacities are split evenly among regions in every country
                     inv_ind = [i for i in inv_busmap[ind]]
@@ -225,7 +234,7 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
                         [i + name_suffix for i in inv_ind],
                         bus=ind,
                         carrier=generator,
-                        p_nom=capacity[ind] / len(inv_ind), # split among regions in a country
+                        p_nom=new_capacity[ind] / len(inv_ind), # split among regions in a country
                         marginal_cost=costs.at[generator,'VOM'],
                         capital_cost=capital_cost,
                         efficiency=costs.at[generator, 'efficiency'],
@@ -238,43 +247,54 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
 
                 p_max_pu = n.generators_t.p_max_pu[capacity.index + name_suffix]
 
-                n.madd("Generator",
-                    capacity.index,
-                    suffix=' ' + generator +"-"+ str(grouping_year),
-                    bus=capacity.index,
-                    carrier=generator,
-                    p_nom=capacity,
-                    marginal_cost=costs.at[generator, 'VOM'],
-                    capital_cost=capital_cost,
-                    efficiency=costs.at[generator, 'efficiency'],
-                    p_max_pu=p_max_pu.rename(columns=n.generators.bus),
-                    build_year=grouping_year,
-                    lifetime=costs.at[generator, 'lifetime']
-                )
+                if not new_build.empty:
+                    n.madd("Generator",
+                        new_capacity.index,
+                        suffix=' ' + generator + name_suffix,
+                        bus=new_capacity.index,
+                        carrier=generator,
+                        p_nom=new_capacity,
+                        marginal_cost=costs.at[generator, 'VOM'],
+                        capital_cost=capital_cost,
+                        efficiency=costs.at[generator, 'efficiency'],
+                        p_max_pu=p_max_pu.rename(columns=n.generators.bus),
+                        build_year=grouping_year,
+                        lifetime=costs.at[generator, 'lifetime']
+                    )
 
         else:
             bus0 = vars(spatial)[carrier[generator]].nodes
             if "EU" not in vars(spatial)[carrier[generator]].locations:
                 bus0 = bus0.intersection(capacity.index + " gas")
 
-            n.madd("Link",
-                capacity.index,
-                suffix= " " + generator +"-" + str(grouping_year),
-                bus0=bus0,
-                bus1=capacity.index,
-                bus2="co2 atmosphere",
-                carrier=generator,
-                marginal_cost=costs.at[generator, 'efficiency'] * costs.at[generator, 'VOM'], #NB: VOM is per MWel
-                capital_cost=costs.at[generator, 'efficiency'] * costs.at[generator, 'fixed'], #NB: fixed cost is per MWel
-                p_nom=capacity / costs.at[generator, 'efficiency'],
-                efficiency=costs.at[generator, 'efficiency'],
-                efficiency2=costs.at[carrier[generator], 'CO2 intensity'],
-                build_year=grouping_year,
-                lifetime=costs.at[generator, 'lifetime']
-            )
+            already_build = n.links.index.intersection(asset_i)
+            new_build = asset_i.difference(n.links.index)
+
+            # this is for the year 2020
+            if not already_build.empty:
+                n.links.loc[already_build, "p_nom_min"] = capacity.loc[already_build.str.replace(name_suffix, "")].values
+
+            if not new_build.empty:
+                new_capacity = capacity.loc[new_build.str.replace(name_suffix, "")]
+
+                n.madd("Link",
+                    new_capacity.index,
+                    suffix= name_suffix,
+                    bus0=bus0,
+                    bus1=new_capacity.index,
+                    bus2="co2 atmosphere",
+                    carrier=generator,
+                    marginal_cost=costs.at[generator, 'efficiency'] * costs.at[generator, 'VOM'], #NB: VOM is per MWel
+                    capital_cost=costs.at[generator, 'efficiency'] * costs.at[generator, 'fixed'], #NB: fixed cost is per MWel
+                    p_nom=capacity / costs.at[generator, 'efficiency'],
+                    efficiency=costs.at[generator, 'efficiency'],
+                    efficiency2=costs.at[carrier[generator], 'CO2 intensity'],
+                    build_year=grouping_year,
+                    lifetime=costs.at[generator, 'lifetime']
+                )
 
 
-def add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years, ashp_cop, gshp_cop, time_dep_hp_cop, costs, default_lifetime):
+def add_heating_capacities_installed_before_baseyear(n, baseyear, ashp_cop, gshp_cop, time_dep_hp_cop, costs, default_lifetime):
     """
     Parameters
     ----------
@@ -294,7 +314,7 @@ def add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years
     # https://ec.europa.eu/energy/studies/mapping-and-analyses-current-and-future-2020-2030-heatingcooling-fuel-deployment_en?redir=1
     # file: "WP2_DataAnnex_1_BuildingTechs_ForPublication_201603.xls" -> "existing_heating_raw.csv".
     # TODO start from original file
-
+    grouping_years = [1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2019]
     # retrieve existing heating capacities
     techs = [
         'gas boiler',
@@ -454,10 +474,10 @@ def add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years
         snakemake = mock_snakemake(
             'add_existing_baseyear',
             simpl='',
-            clusters="37",
+            clusters="45",
             lv=1.0,
             opts='',
-            sector_opts='168H-T-H-B-I-solar+p3-dist1',
+            sector_opts='cb40ex0-365H-T-H-B-I-A-solar+p3-dist1',
             planning_horizons=2020,
         )
 
@@ -493,6 +513,6 @@ def add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years
         ashp_cop = xr.open_dataarray(snakemake.input.cop_air_total).to_pandas().reindex(index=n.snapshots)
         gshp_cop = xr.open_dataarray(snakemake.input.cop_soil_total).to_pandas().reindex(index=n.snapshots)
         default_lifetime = snakemake.config['costs']['lifetime']
-        add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years, ashp_cop, gshp_cop, time_dep_hp_cop, costs, default_lifetime)
+        add_heating_capacities_installed_before_baseyear(n, baseyear, ashp_cop, gshp_cop, time_dep_hp_cop, costs, default_lifetime)
 
     n.export_to_netcdf(snakemake.output[0])

scripts/build_energy_totals.py

@@ -127,17 +127,16 @@ def reverse(dictionary):
 }
 
 
-def build_eurostat(countries, year):
+def build_eurostat(input_eurostat, countries, report_year,  year):
     """Return multi-index for all countries' energy data in TWh/a."""
 
-    report_year = snakemake.config["energy"]["eurostat_report_year"]
     filenames = {
         2016: f"/{year}-Energy-Balances-June2016edition.xlsx",
         2017: f"/{year}-ENERGY-BALANCES-June2017edition.xlsx"
     }
 
     dfs = pd.read_excel(
-        snakemake.input.eurostat + filenames[report_year],
+        input_eurostat + filenames[report_year],
         sheet_name=None,
         skiprows=1,
         index_col=list(range(4)),
@@ -563,18 +562,18 @@ def build_energy_totals(countries, eurostat, swiss, idees):
     return df
 
 
-def build_eea_co2(year=1990):
+def build_eea_co2(input_co2, year=1990, emissions_scope="CO2"):
 
     # 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)
-    df = pd.read_csv(snakemake.input.co2, encoding="latin-1")
+    df = pd.read_csv(input_co2, encoding="latin-1")
 
     df.replace(dict(Year="1985-1987"), 1986, inplace=True)
     df.Year = df.Year.astype(int)
     index_col = ["Country_code", "Pollutant_name", "Year", "Sector_name"]
     df = df.set_index(index_col).sort_index()
 
-    emissions_scope = snakemake.config["energy"]["emissions"]
+    emissions_scope = emissions_scope
 
     cts = ["CH", "EUA", "NO"] + eu28_eea
 
@@ -611,9 +610,9 @@ def build_eea_co2(year=1990):
     return emissions / 1e3
 
 
-def build_eurostat_co2(countries, year=1990):
+def build_eurostat_co2(input_eurostat, countries, report_year, year=1990):
 
-    eurostat = build_eurostat(countries, year)
+    eurostat = build_eurostat(input_eurostat, countries, report_year, year)
 
     specific_emissions = pd.Series(index=eurostat.columns, dtype=float)
 
@@ -702,16 +701,19 @@ def build_transport_data(countries, population, idees):
     idees_countries = countries.intersection(eu28)
 
     data_year = config["energy_totals_year"]
-    eurostat = build_eurostat(countries, data_year)
+    report_year = snakemake.config["energy"]["eurostat_report_year"]
+    input_eurostat = snakemake.input.eurostat
+    eurostat = build_eurostat(input_eurostat, countries, report_year, data_year)
     swiss = build_swiss(data_year)
     idees = build_idees(idees_countries, data_year)
 
     energy = build_energy_totals(countries, eurostat, swiss, idees)
     energy.to_csv(snakemake.output.energy_name)
 
     base_year_emissions = config["base_emissions_year"]
-    eea_co2 = build_eea_co2(base_year_emissions)
-    eurostat_co2 = build_eurostat_co2(countries, base_year_emissions)
+    emissions_scope = snakemake.config["energy"]["emissions"]
+    eea_co2 = build_eea_co2(snakemake.input.co2, base_year_emissions, emissions_scope)
+    eurostat_co2 = build_eurostat_co2(input_eurostat, countries, report_year, base_year_emissions)
 
     co2 = build_co2_totals(countries, eea_co2, eurostat_co2)
     co2.to_csv(snakemake.output.co2_name)

scripts/plot_summary.py

@@ -202,7 +202,7 @@ def plot_energy():
     new_index = preferred_order.intersection(df.index).append(df.index.difference(preferred_order))
 
     new_columns = df.columns.sort_values()
-    
+
     fig, ax = plt.subplots(figsize=(12,8))
 
     print(df.loc[new_index, new_columns])
@@ -363,7 +363,7 @@ def historical_emissions(cts):
 
 
 
-def plot_carbon_budget_distribution():
+def plot_carbon_budget_distribution(input_eurostat):
     """
     Plot historical carbon emissions in the EU and decarbonization path
     """
@@ -385,9 +385,10 @@ def plot_carbon_budget_distribution():
     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)
+    pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
+    countries=pd.read_csv(snakemake.input.country_codes, index_col=1)
     cts=countries.index.to_list()
-    e_1990 = co2_emissions_year(cts, opts, year=1990)
+    e_1990 = co2_emissions_year(cts, input_eurostat, opts, year=1990)
     CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv',
                         index_col=0)
 
@@ -438,8 +439,7 @@ def plot_carbon_budget_distribution():
     if 'snakemake' not in globals():
         from helper import mock_snakemake
         snakemake = mock_snakemake('plot_summary')
-
-    update_config_with_sector_opts(snakemake.config, snakemake.wildcards.sector_opts)
+
 
     n_header = 4
 
@@ -453,4 +453,4 @@ def plot_carbon_budget_distribution():
         opts=sector_opts.split('-')
         for o in opts:
             if "cb" in o:
-                plot_carbon_budget_distribution()
+                plot_carbon_budget_distribution(snakemake.input.eurostat)