Additional fixes on plots

Snakefile

@@ -109,6 +109,39 @@ rule plot_electricity_for_heats:
         ),
 
 
+rule plot_heat_saving:
+    params:
+        clusters=config["plotting"]["clusters"],
+        planning_horizon=config["plotting"]["planning_horizon"],
+    output:
+        figure=RESULTS+"plot_heat_savings_{clusters}_{planning_horizon}.png",
+        figure_full=RESULTS+"plot_heat_savings_full_year_{clusters}_{planning_horizon}.png",
+        figure_flexibility=RESULTS+"plot_heat_flexibility_{clusters}_{planning_horizon}.png",
+    resources:
+        mem_mb=10000,
+    script:
+        "plots/plot_heat_savings.py"
+
+
+rule plot_heat_savings:
+    input:
+        expand(
+            RESULTS
+            + "plot_heat_savings_{clusters}_{planning_horizon}.png",
+            **config["plotting"],
+        ),
+        expand(
+            RESULTS
+            + "plot_heat_savings_full_year_{clusters}_{planning_horizon}.png",
+            **config["plotting"],
+        ),
+        expand(
+            RESULTS
+            + "plot_heat_flexibility_{clusters}_{planning_horizon}.png",
+            **config["plotting"],
+        ),
+
+
 rule plot_electricity_generation:
     params:
         clusters=config["plotting"]["clusters"],
@@ -366,6 +399,21 @@ rule plot_all:
             + "plot_electricity_for_heat_{clusters}_{planning_horizon}.png",
             **config["plotting"],
         ),
+        expand(
+            RESULTS
+            + "plot_heat_savings_{clusters}_{planning_horizon}.png",
+            **config["plotting"],
+        ),
+        expand(
+            RESULTS
+            + "plot_heat_savings_full_year_{clusters}_{planning_horizon}.png",
+            **config["plotting"],
+        ),
+        expand(
+            RESULTS
+            + "plot_heat_flexibility_{clusters}_{planning_horizon}.png",
+            **config["plotting"],
+        ),
         expand(
             RESULTS
             + "table_heat_pumps_{clusters}.csv",

configs/EEE_study/config.retro_tes-industry_2050.yaml

@@ -121,3 +121,4 @@ solving:
       threads: 8
       BarConvTol: 1.e-4
       BarHomogeneous: 1
+      NumericFocus: 3

plots/_helpers.py

@@ -171,6 +171,13 @@
      ]
 )
 
+# day-ahead prices for 2021 URL:https://www.ffe.de/en/publications/european-day-ahead-electricity-prices-in-2022/
+HISTORIC_PRICES = {"AT": 106.9, "BE": 104.1, "CH": 114.9, "CZ": 100.7, "DE": 96.8,
+                   "DK": 88.1, "EE": 86.7, "ES": 111.9, "FI": 72.3, "FR": 109.2,
+                   "GR": 116.4, "HU": 113.9, "IT": 125.2, "LT": 90.5, "LV": 88.8,
+                   "NL": 103.0, "NO": 74.7, "PL": 87.0, "PT": 112.0, "RO": 110.9,
+                   "RS": 114.0, "SE": 66.0, "SI": 115.0, "SK": 102.8}
+
 
 def rename_techs(label):
 

plots/plot_capacity_expansion.py

@@ -39,19 +39,29 @@ def plot_capacities(caps_df, clusters, planning_horizon, plot_width=7):
     df = df[df.index != 'solid biomass transport']
 
     # convert to GW
-    df = df / 1e3
+    df = df / 1e6
     df = df.groupby(df.index.map(rename_techs)).sum()
 
-    caps_threshold = 10
+    caps_threshold = 10e-3
     to_drop = df.index[df.max(axis=1) < caps_threshold]  #df <
 
     logger.info(
-        f"Dropping technology with capacity below {caps_threshold} GW"
+        f"Dropping technology with capacity below {caps_threshold} TW"
     )
     logger.debug(df.loc[to_drop])
 
     df = df.drop(to_drop)
 
+    # remove extra technologies that are not expanded
+    drop_techs = ["gas boiler", "process emissions",
+                  "solid biomass", "solid biomass CHP",
+                  "gas storage",
+                  "kerosene for aviation", "land transport oil", "naphtha for industry",
+                  "shipping methanol", "shipping oil",
+                  "coal", "lignite", "oil",
+                  "coal for industry", "gas for industry"]
+    df = df.drop(set(df.index).intersection(set(drop_techs)))
+
     logger.info(f"Total optimal capacity is {round(df.sum())} GW")
 
     new_index = PREFERRED_ORDER.intersection(df.index).append(
@@ -75,11 +85,11 @@ def plot_capacities(caps_df, clusters, planning_horizon, plot_width=7):
 
     plt.xticks(rotation=0, fontsize=10)
 
-    ax.set_ylabel("Capacity expansion [GW]")
+    ax.set_ylabel("Capacity expansion [TW]")
 
     ax.set_xlabel("")
-    ax.set_ylim([0,16000])
-    ax.set_yticks(np.arange(0, 17000, 2000))
+    ax.set_ylim([0,16])
+    ax.set_yticks(np.arange(0, 17, 2))
 
     x_ticks = list(df.columns)
     if planning_horizon in ["2040", "2050"] and "Limited \nRenovation &\nCost-Optimal Heating" in x_ticks:

plots/plot_co2_level.py

@@ -110,18 +110,18 @@
         "solid biomass for industry",
         "gas for industry",
         "coal for industry",
-        "methanol",
-        "oil",
-        "lignite",
-        "coal",
         "shipping oil",
         "shipping methanol",
         "naphtha for industry",
         "land transport oil",
         "kerosene for aviation",
-
-        "transmission lines",
-        "distribution lines",
+        "process emissions",
+        "SMR",
+        "coal",
+        "lignite",
+        "methanol",
+        "oil",
+
         "gas pipeline",
         "H2 pipeline",
 
@@ -132,26 +132,17 @@
         "methanation",
         "BEV charger",
         "V2G",
-        "SMR",
         "methanolisation",
 
         "battery storage",
         "gas storage",
         "H2 storage",
         "TES",
-
-        "hydroelectricity",
         "OCGT",
         "CCGT",
-        "onshore wind",
-        "offshore wind",
-        "solar PV",
-        "solar thermal",
-        "solar rooftop",
 
         "co2",
         "CO2 sequestration",
-        "process emissions",
 
         "gas CHP",
         "resistive heater",
@@ -195,12 +186,14 @@ def get_co2_balance(n, nice_name):
 
 
 def plot_co2_balance(co2_df, clusters, planning_horizon, plot_width=7):
+    # convert to Billion tCO2_eq
+    co2_df = co2_df / 1e9
     # filter out technologies with very small emission
     max_emissions = co2_df.abs().sum().max() / 2
     co2_threshold = max_emissions / 100 # 1% of max
     to_drop = co2_df.index[co2_df.abs().max(axis=1) < co2_threshold]
     logger.info(
-        f"Dropping technology with co2 balance below {co2_threshold} ton CO2_eq per year"
+        f"Dropping technology with co2 balance below {co2_threshold} Bton CO2_eq per year"
     )
     logger.debug(co2_df.loc[to_drop])
     co2_df = co2_df.drop(to_drop)
@@ -229,10 +222,11 @@ def plot_co2_balance(co2_df, clusters, planning_horizon, plot_width=7):
     labels = labels[:-neg_co2_length] + labels[-neg_co2_length:][::-1]
 
     plt.xticks(rotation=0, fontsize=10)
-    ax.set_ylabel("CO$_2$ emissions [tCO$_{2-eq}$]")
+    ax.set_ylabel("CO$_2$ emissions [BtCO$_{2-eq}$]")
     ax.set_xlabel("")
-    ax.set_ylim([-4e9,4e9])
-    ax.set_yticks(np.arange(-4.0e9, 4.0e9, 5e8))
+    ax.set_yticks(np.arange(-4.0, 4.0, 0.5))
+    ax.set_ylim([-3.6,3.6])
+
     x_ticks = list(co2_df.columns)
     if planning_horizon in ["2040", "2050"] and "Limited \nRenovation &\nOptimal Heating" in x_ticks:
         # replace name for Limited Renovation scenario for 2030 to be LROH
@@ -372,7 +366,7 @@ def fill_table_df(df, planning_horizon, scenarios, values):
             co2_savings_df = co2_emission_BAU - table_emissions_df
             co2_savings_df.index = ["CO2 emission savings [tCO2_eq]"]
             # saving in cubic meter of gas (1 m3 natural gas = 1,9 kg CO2) Source: https://www.eeagrants.gov.pt/media/2776/conversion-guidelines.pdf
-            co2_savings_df.loc["Equivalent gas savings [m3]"] = co2_savings_df.loc["CO2 emission savings [tCO2_eq]", :] * 1000 / 1.9
+            co2_savings_df.loc["Equivalent gas savings [trillion m3]"] = co2_savings_df.loc["CO2 emission savings [tCO2_eq]", :] * 1000 / 1.9 / 1e12
 
 
     # move to base directory

plots/plot_electricity_bills.py

@@ -14,7 +14,7 @@
 warnings.filterwarnings("ignore")
 from _helpers import mock_snakemake, update_config_from_wildcards, load_network, \
                      change_path_to_pypsa_eur, change_path_to_base, \
-                     CO2L_LIMITS, LINE_LIMITS, BAU_HORIZON
+                     CO2L_LIMITS, LINE_LIMITS, BAU_HORIZON, HISTORIC_PRICES
 
 
 def get_households():
@@ -123,6 +123,10 @@ def electricity_bills(network, households):
     total_cost_country += total_gas_cost_country
     # electricity bill per household [EUR/household] (households given in thousands)
     elec_bills_household = total_cost_country / (households*1e3)
+
+    # add EU average bills
+    average_bills_household = total_cost_country.sum() / (households.sum() * 1e3)
+    elec_bills_household["EUR"] = average_bills_household
 
     return elec_bills_household
 
@@ -141,7 +145,14 @@ def plot_electricity_cost(df_prices, name):
                    "Optimal Renovation and Electric Heating":"limegreen", 
                    "Limited Renovation and Cost-Optimal Heating":"royalblue", 
                    "No Renovation and Electric Heating":"#f4b609",
-                   "BAU": "grey"}
+                   "BAU": "grey",
+                   "Historic data": '#900C3F'}
+
+    # add historic electricity price for BAU plot
+    if "BAU" in df_prices.index and name == "prices":
+        historic_prices = pd.Series(HISTORIC_PRICES).to_frame().T
+        historic_prices.index = ["Historic data"]
+        sorted_df_prices = pd.concat([sorted_df_prices, historic_prices], axis=0)
 
     # plot as bar plot
     fig, ax = plt.subplots(figsize=(7,3))
@@ -243,9 +254,9 @@ def electricity_prices(network, households):
     # electricity bill per MWh [EUR/MWh]
     energy_price_MWh = prices / total_load_country
 
-    # drop EU
-    if "EU" in energy_price_MWh.index:
-        energy_price_MWh.drop("EU", axis=0, inplace=True)
+    # add EU average 
+    average_price_MWh = total_cost.sum() / total_load_country.sum()
+    energy_price_MWh["EUR"] = average_price_MWh
 
     return energy_price_MWh
 
@@ -256,7 +267,7 @@ def electricity_prices(network, households):
         snakemake = mock_snakemake(
             "plot_electricity_bill", 
             clusters="48",
-            planning_horizon="2030",
+            planning_horizon="2020",
         )
     # update config based on wildcards
     config = update_config_from_wildcards(snakemake.config, snakemake.wildcards)

plots/plot_electricity_for_heat.py

@@ -66,10 +66,10 @@ def get_elec_consumption_for_heat(n):
     return elec_demand_f_heating
 
 
-def plot_elec_consumption_for_heat(dict_elec):
+def plot_elec_consumption_for_heat(dict_elec, horizon):
     # set heights for each subplots
     if "BAU" in dict_elec.keys():
-        heights = [1.4]
+        heights = [1.5]
     else:
         heights = [1.4] * 4
     fig = plt.figure(figsize=(6.4, sum(heights)))
@@ -85,7 +85,7 @@ def plot_elec_consumption_for_heat(dict_elec):
         ax.set_facecolor("whitesmoke")
 
         print("Hard coded coordinates on x-axis, selected for 3H granularity")
-        where = [359, 527]
+        where = [359, 695]
 
         elec_demand_f_heating = elec_demand_f_heating.reset_index()[["ground heat pump", "air heat pump", "resistive heater", "gas boiler"]]
         colors = ["#2fb537", "#48f74f", "#d8f9b8", "#db6a25"]
@@ -95,20 +95,20 @@ def plot_elec_consumption_for_heat(dict_elec):
             linewidth=0, ax=ax
         )
 
+        cumulative = (elec_demand_f_heating / 1e3).cumsum(axis=1)
+        ax.plot(cumulative.iloc[where[0]:where[1]], color='black', linewidth=0.5)
+
         ax.set_xlabel("", fontsize=12)
-        ax.set_ylim([1,2000])
-        ax.set_yscale("log")  # Set the y-axis to logarithmic scale
+        # get in y limits for all horizons
+        y_limits = {'2020':1600, '2030':1100, '2040':900, '2050': 900}
 
-        # Custom y-ticks for logarithmic scale
-        ax.yaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=10))
-        ax.yaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs='auto'))
-        ax.yaxis.set_minor_formatter(ticker.NullFormatter())
+        ax.set_ylim([1, y_limits[horizon]])
 
         if i < 3:
             ax.set_xticks([])
             ax.set_xlabel("")
         if i == 3 or "BAU" in dict_elec.keys():
-            ticks = [i for i in range(where[0], where[1], 24)]
+            ticks = [i for i in range(where[0], where[1], 48)]
             ax.set_xticks(ticks)  # Set the tick positions
             ticks = [
                 str(n.snapshots[i]).split(" ")[0].split("-")[2]+"."+str(n.snapshots[i]).split(" ")[0].split("-")[1]
@@ -136,10 +136,10 @@ def plot_elec_consumption_for_heat(dict_elec):
         loc=[1.02, -.2], fontsize=10
     )
     if len(dict_elec.keys()) == 1:
-        ylabel = "Electricity and Gas \nConsumption to Cover \nHeating Demands [GW]"
+        ylabel = "Electricity and Gas \nConsumption to Cover \nHeating Demands [$GW_{el}$]"
         axes[0].set_ylabel(ylabel, fontsize=10)
     else:
-        ylabel = "Electricity and Gas Consumption to \n Cover Heating Demands [GW]"
+        ylabel = "Electricity and Gas Consumption to \n Cover Heating Demands [$GW_{el}$]"
         axes[2].set_ylabel(ylabel, fontsize=10)
         axes[2].yaxis.set_label_coords(-0.1, 1)
     plt.subplots_adjust(hspace=0.3)
@@ -203,5 +203,5 @@ def plot_elec_consumption_for_heat(dict_elec):
         total_elec_consumption_for_heat[name] = elec_consumption_for_heat
 
     # plot and store electricity prices
-    plot_elec_consumption_for_heat(total_elec_consumption_for_heat)
+    plot_elec_consumption_for_heat(total_elec_consumption_for_heat, planning_horizon)