complete structure for plotting electricity production

Snakefile

@@ -107,6 +107,6 @@ rule sync:
     shell:
         """
         rsync -uvarh --no-g --ignore-missing-args --files-from=.sync-send . {params.cluster}
-        rsync -uvarh --no-g --ignore-missing-args {params.cluster}/results results
+        rsync -uvarh --no-g {params.cluster}/results results
-        rsync -uvarh --no-g --ignore-missing-args {params.cluster}/logs logs
+        rsync -uvarh --no-g {params.cluster}/logs logs
         """

config/config.default.yaml

@@ -9,6 +9,10 @@ logging:
   level: INFO
   format: '%(levelname)s:%(name)s:%(message)s'
 
+private:
+  keys:
+    entsoe_api:
+
 remote:
   ssh: ""
   path: ""

matplotlibrc

@@ -4,3 +4,4 @@
 font.family: sans-serif
 font.sans-serif: Ubuntu, DejaVu Sans
 image.cmap: viridis
+figure.autolayout : True

rules/build_electricity.smk

@@ -331,7 +331,7 @@ rule add_electricity:
         BENCHMARKS + "add_electricity"
     threads: 1
     resources:
-        mem_mb=5000,
+        mem_mb=10000,
     conda:
         "../envs/environment.yaml"
     script:
@@ -365,7 +365,7 @@ rule simplify_network:
         BENCHMARKS + "simplify_network/elec_s{simpl}"
     threads: 1
     resources:
-        mem_mb=4000,
+        mem_mb=10000,
     conda:
         "../envs/environment.yaml"
     script:
@@ -406,7 +406,7 @@ rule cluster_network:
         BENCHMARKS + "cluster_network/elec_s{simpl}_{clusters}"
     threads: 1
     resources:
-        mem_mb=6000,
+        mem_mb=10000,
     conda:
         "../envs/environment.yaml"
     script:
@@ -429,7 +429,7 @@ rule add_extra_components:
         BENCHMARKS + "add_extra_components/elec_s{simpl}_{clusters}_ec"
     threads: 1
     resources:
-        mem_mb=3000,
+        mem_mb=4000,
     conda:
         "../envs/environment.yaml"
     script:

rules/collect.smk

@@ -79,6 +79,11 @@ rule validate_elec_networks:
     input:
         expand(
             RESULTS
-            + "figures/validate_electricity_production_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
+            + "figures/.statistics_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
+            **config["scenario"]
+        ),
+        expand(
+            RESULTS
+            + "figures/.validation_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
             **config["scenario"]
         ),

rules/common.smk

@@ -16,7 +16,7 @@ def memory(w):
             factor *= int(m.group(1)) / 8760
             break
     if w.clusters.endswith("m") or w.clusters.endswith("c"):
-        return int(factor * (18000 + 180 * int(w.clusters[:-1])))
+        return int(factor * (35000 + 180 * int(w.clusters[:-1])))
     elif w.clusters == "all":
         return int(factor * (18000 + 180 * 4000))
     else:

rules/postprocess.smk

@@ -148,12 +148,33 @@ rule plot_summary:
         "../scripts/plot_summary.py"
 
 
+STATISTICS_BARPLOTS = [
+    "capacity_factor",
+    "installed_capacity",
+    "optimal_capacity",
+    "capital_expenditure",
+    "operational_expenditure",
+    "curtailment",
+    "supply",
+    "withdrawal",
+    "market_value",
+]
+
+
 rule plot_statistics:
+    params:
+        plotting=config["plotting"],
+        barplots=STATISTICS_BARPLOTS,
     input:
         overrides="data/override_component_attrs",
         network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
     output:
-        bar=RESULTS
+        **{
-        + "figures/statistics_bar_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.pdf",
+            f"{plot}_bar": RESULTS
+            + f"figures/statistics_{plot}_bar_elec_s{{simpl}}_{{clusters}}_ec_l{{ll}}_{{opts}}.pdf"
+            for plot in STATISTICS_BARPLOTS
+        },
+        barplots_touch=RESULTS
+        + "figures/.statistics_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
     script:
         "../scripts/plot_statistics.py"

rules/retrieve.smk

@@ -141,11 +141,13 @@ if config["sector"]["gas_network"] or config["sector"]["H2_retrofit"]:
 
 
 rule retrieve_electricity_demand:
-    params:
-        version="2019-06-05" if config["snapshots"]["end"] < "2019" else "latest",
     input:
         HTTP.remote(
-            "data.open-power-system-data.org/time_series/{params.version}/time_series_60min_singleindex.csv",
+            "data.open-power-system-data.org/time_series/{version}/time_series_60min_singleindex.csv".format(
+            version="2019-06-05"
+                if config["snapshots"]["end"] < "2019"
+                else "latest"
+            ),
             keep_local=True,
             static=True,
         ),

rules/solve_electricity.smk

@@ -55,7 +55,7 @@ rule solve_operations_network:
         )
     threads: 4
     resources:
-        mem_mb=(lambda w: 5000 + 372 * int(w.clusters)),
+        mem_mb=(lambda w: 10000 + 372 * int(w.clusters)),
     shadow:
         "minimal"
     conda:

rules/validate.smk

@@ -18,15 +18,23 @@ rule build_electricity_production:
     resources:
         mem_mb=5000,
     script:
-        "../scripts/retrieve_electricity_production.py"
+        "../scripts/build_electricity_production.py"
+
+
+PLOTS = ["production_bar", "production_deviation_bar", "seasonal_operation_area"]
 
 
 rule plot_electricity_production:
     input:
         network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
-        electricity_production="data/historical_electricity_production.csv",
+        electricity_production=RESOURCES + "historical_electricity_production.csv",
     output:
-        electricity_producion=RESULTS
+        **{
-        + "figures/validate_electricity_production_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.pdf",
+            plot: RESULTS
+            + f"figures/validation_{plot}_elec_s{{simpl}}_{{clusters}}_ec_l{{ll}}_{{opts}}.pdf"
+            for plot in PLOTS
+        },
+        plots_touch=RESULTS
+        + "figures/.validation_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
     script:
-        "scripts/plot_electricity_production.py"
+        "../scripts/plot_electricity_production.py"

scripts/build_electricity_production.py

@@ -3,16 +3,11 @@
 # SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
-"""
-Created on Mon Jul  3 11:19:54 2023.
-
-@author: fabian
-"""
-
 
 import logging
 
 import pandas as pd
+from _helpers import configure_logging
 from entsoe import EntsoePandasClient
 from entsoe.exceptions import NoMatchingDataError
 
@@ -31,6 +26,10 @@ carrier_grouper = {
     "Fossil Brown coal/Lignite": "Lignite",
     "Fossil Peat": "Lignite",
     "Fossil Hard coal": "Coal",
+    "Wind Onshore": "Onshore Wind",
+    "Wind Offshore": "Offshore Wind",
+    "Other renewable": "Other",
+    "Marine": "Other",
 }
 
 
@@ -38,43 +37,37 @@ if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
 
-        snakemake = mock_snakemake("retrieve_historical_electricity_generation")
+        snakemake = mock_snakemake("build_electricity_production")
+    configure_logging(snakemake)
 
+    api_key = snakemake.config["private"]["keys"]["entsoe_api"]
+    client = EntsoePandasClient(api_key=api_key)
 
-api_key = "aeff3346-a240-40df-bd12-692772b845d0"
+    start = pd.Timestamp(snakemake.params.snapshots["start"], tz="Europe/Brussels")
-client = EntsoePandasClient(api_key=api_key)
+    end = pd.Timestamp(snakemake.params.snapshots["end"], tz="Europe/Brussels")
 
-start = pd.Timestamp(snakemake.params.snapshots["start"], tz="Europe/Brussels")
+    countries = snakemake.params.countries
-end = pd.Timestamp(snakemake.params.snapshots["end"], tz="Europe/Brussels")
 
-countries = snakemake.params.countries
+    generation = []
+    unavailable_countries = []
 
+    for country in countries:
+        country_code = country
 
-generation = []
+        try:
-unavailable_countries = []
+            gen = client.query_generation(country, start=start, end=end, nett=True)
+            gen = gen.tz_localize(None).resample("1h").mean()
+            gen = gen.loc[start.tz_localize(None) : end.tz_localize(None)]
+            gen = gen.rename(columns=carrier_grouper).groupby(level=0, axis=1).sum()
+            generation.append(gen)
+        except NoMatchingDataError:
+            unavailable_countries.append(country)
 
-for country in countries:
+    if unavailable_countries:
-    country_code = country
+        logger.warning(
+            f"Historical electricity production for countries {', '.join(unavailable_countries)} not available."
+        )
 
-    try:
+    keys = [c for c in countries if c not in unavailable_countries]
-        gen = client.query_generation(country, start=start, end=end, nett=True)
+    generation = pd.concat(generation, keys=keys, axis=1)
-        gen = gen.tz_localize(None).resample("1h").mean()
+    generation.to_csv(snakemake.output[0])
-        gen = gen.rename(columns=carrier_grouper).groupby(level=0, axis=1).sum()
-        generation.append(gen)
-    except NoMatchingDataError:
-        unavailable_countries.append(country)
-
-
-if unavailable_countries:
-    logger.warning(
-        f"Historical electricity production for countries {', '.join(unavailable_countries)} not available."
-    )
-
-keys = [c for c in countries if c not in unavailable_countries]
-generation = pd.concat(generation, keys=keys, axis=1)
-generation = generation.loc[start.tz_localize(None) : end.tz_localize(None)]
-# -*- coding: utf-8 -*-
-# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors
-#
-# SPDX-License-Identifier: MIT
-generation.to_csv(snakemake.output[0])

scripts/plot_electricity_production.py

@@ -0,0 +1,146 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+import matplotlib.pyplot as plt
+import pandas as pd
+import pypsa
+import seaborn as sns
+from _helpers import configure_logging
+from pypsa.statistics import get_bus_and_carrier
+
+sns.set_theme("paper", style="whitegrid")
+
+carrier_groups = {
+    "Offshore Wind (AC)": "Offshore Wind",
+    "Offshore Wind (DC)": "Offshore Wind",
+    "Open-Cycle Gas": "Gas",
+    "Combined-Cycle Gas": "Gas",
+    "Reservoir & Dam": "Hydro",
+    "Pumped Hydro Storage": "Hydro",
+}
+
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake(
+            "plot_electricity_production",
+            simpl="",
+            opts="Ept-12h",
+            clusters="37",
+            ll="v1.0",
+        )
+    configure_logging(snakemake)
+
+    n = pypsa.Network(snakemake.input.network)
+    n.loads.carrier = "load"
+
+    historic = pd.read_csv(
+        snakemake.input.electricity_production,
+        index_col=0,
+        header=[0, 1],
+        parse_dates=True,
+    )
+    historic = historic.drop("Other renewable", axis=1, level=1)
+    historic = historic.drop("Marine", axis=1, level=1)
+
+    colors = n.carriers.set_index("nice_name").color.where(
+        lambda s: s != "", "lightgrey"
+    )
+    colors["Offshore Wind"] = colors["Offshore Wind (AC)"]
+    colors["Gas"] = colors["Combined-Cycle Gas"]
+    colors["Hydro"] = colors["Reservoir & Dam"]
+    colors["Other"] = "lightgray"
+
+    if len(historic.index) > len(n.snapshots):
+        historic = historic.resample(n.snapshots.inferred_freq).mean().loc[n.snapshots]
+
+    optimized = n.statistics.dispatch(
+        groupby=get_bus_and_carrier, aggregate_time=False
+    ).T
+    optimized = optimized[["Generator", "StorageUnit"]].droplevel(0, axis=1)
+    optimized = optimized.rename(columns=n.buses.country, level=0)
+    optimized = optimized.rename(columns=carrier_groups, level=1)
+    optimized = optimized.groupby(axis=1, level=[0, 1]).sum()
+
+    data = pd.concat([historic, optimized], keys=["Historic", "Optimized"], axis=1)
+    data.columns.names = ["Kind", "Country", "Carrier"]
+    data = data.mul(n.snapshot_weightings.generators, axis=0)
+
+    # %% total production per carrier
+    fig, ax = plt.subplots(figsize=(6, 6))
+
+    df = data.groupby(level=["Kind", "Carrier"], axis=1).sum().sum().unstack().T
+    df = df / 1e6  # TWh
+    df.plot.barh(ax=ax, xlabel="Electricity Production [TWh]", ylabel="")
+    ax.grid(axis="y")
+    fig.savefig(snakemake.output.production_bar, bbox_inches="tight")
+
+    # %% highest diffs
+
+    fig, ax = plt.subplots(figsize=(6, 10))
+
+    df = data.sum() / 1e6  # TWh
+    df = df["Optimized"] - df["Historic"]
+    df = df.dropna().sort_values()
+    df = pd.concat([df.iloc[:5], df.iloc[-5:]])
+    c = colors[df.index.get_level_values(1)]
+    df.plot.barh(
+        xlabel="Optimized Production - Historic Production [TWh]", ax=ax, color=c.values
+    )
+    ax.set_title("Strongest Deviations")
+    ax.grid(axis="y")
+    fig.savefig(snakemake.output.production_deviation_bar, bbox_inches="tight")
+
+    # %% seasonal operation
+
+    fig, axes = plt.subplots(3, 1, figsize=(9, 9))
+
+    df = (
+        data.groupby(level=["Kind", "Carrier"], axis=1)
+        .sum()
+        .resample("1W")
+        .mean()
+        .clip(lower=0)
+    )
+    df = df / 1e3
+
+    order = (
+        (df["Historic"].diff().abs().sum() / df["Historic"].sum()).sort_values().index
+    )
+    c = colors[order]
+    optimized = df["Optimized"].reindex(order, axis=1, level=1)
+    historical = df["Historic"].reindex(order, axis=1, level=1)
+
+    kwargs = dict(color=c, legend=False, ylabel="Production [GW]", xlabel="")
+
+    optimized.plot.area(ax=axes[0], **kwargs, title="Optimized")
+    historical.plot.area(ax=axes[1], **kwargs, title="Historic")
+
+    diff = historical - optimized
+    diff.clip(lower=0).plot.area(
+        ax=axes[2], **kwargs, title="$\Delta$ (Optimized - Historic)"
+    )
+    lim = axes[2].get_ylim()[1]
+    diff.clip(upper=0).plot.area(ax=axes[2], **kwargs)
+    axes[2].set_ylim(bottom=-lim, top=lim)
+
+    h, l = axes[0].get_legend_handles_labels()
+    fig.legend(
+        h[::-1],
+        l[::-1],
+        loc="center left",
+        bbox_to_anchor=(1, 0.5),
+        ncol=1,
+        frameon=False,
+        labelspacing=1,
+    )
+    fig.savefig(snakemake.output.seasonal_operation_area, bbox_inches="tight")
+
+    # touch file
+    with open(snakemake.output.plots_touch, "a"):
+        pass

scripts/plot_historic_comparison.py

@@ -1,47 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors
-#
-# SPDX-License-Identifier: MIT
-"""
-Created on Mon Jul  3 12:50:26 2023.
-
-@author: fabian
-"""
-
-import pandas as pd
-import pypsa
-from pypsa.statistics import get_bus_and_carrier
-
-carrier_groups = {
-    "Offshore Wind (AC)": "Offshore Wind",
-    "Offshore Wind (DC)": "Offshore Wind",
-    "Open-Cycle Gas": "Gas",
-    "Combined-Cycle Gas": "Gas",
-    "Reservoir & Dam": "Hydro",
-    "Pumped Hydro Storage": "Hydro",
-}
-
-if __name__ == "__main__":
-    if "snakemake" not in globals():
-        from _helpers import mock_snakemake
-
-        snakemake = mock_snakemake(
-            "plot_statistics",
-            simpl="",
-            opts="Co2L-3h",
-            clusters="37c",
-            ll="v1.0",
-        )
-
-n = pypsa.Network(snakemake.input.network)
-historic = pd.read_csv(
-    snakemake.input.historic_electricity_generation, index_col=0, header=[0, 1]
-)
-
-
-simulated = n.statistics.dispatch(groupby=get_bus_and_carrier, aggregate_time=False).T
-simulated = simulated[["Generator", "StorageUnit"]].droplevel(0, axis=1)
-simulated = simulated.rename(columns=n.buses.country, level=0)
-simulated = simulated.rename(carrier_groups, level=1)
-simulated = simulated.groupby(axis=1, level=[0, 1]).sum()

scripts/plot_statistics.py

@@ -3,15 +3,11 @@
 # SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
-"""
-Created on Fri Jun 30 10:50:53 2023.
-
-@author: fabian
-"""
 
 import matplotlib.pyplot as plt
 import pypsa
 import seaborn as sns
+from _helpers import configure_logging
 
 sns.set_theme("paper", style="whitegrid")
 
@@ -23,80 +19,94 @@ if __name__ == "__main__":
         snakemake = mock_snakemake(
             "plot_statistics",
             simpl="",
-            opts="Co2L-3h",
+            opts="Ept-12h",
-            clusters="37c",
+            clusters="37",
             ll="v1.0",
         )
+    configure_logging(snakemake)
 
+    n = pypsa.Network(snakemake.input.network)
 
-n = pypsa.Network(snakemake.network)
+    n.loads.carrier = "load"
+    n.carriers.loc["load", ["nice_name", "color"]] = "Load", "darkred"
+    colors = n.carriers.set_index("nice_name").color.where(
+        lambda s: s != "", "lightgrey"
+    )
 
+    # %%
 
-n.loads.carrier = "load"
+    def rename_index(ds):
-n.carriers.loc["load", ["nice_name", "color"]] = "Load", "darkred"
+        return ds.set_axis(ds.index.map(lambda x: f"{x[1]}\n({x[0].lower()})"))
-colors = n.carriers.set_index("nice_name").color.where(lambda s: s != "", "lightgrey")
 
+    def plot_static_per_carrier(ds, ax, drop_zero=True):
+        if drop_zero:
+            ds = ds[ds != 0]
+        ds = ds.dropna()
+        c = colors[ds.index.get_level_values("carrier")]
+        ds = ds.pipe(rename_index)
+        label = f"{ds.attrs['name']} [{ds.attrs['unit']}]"
+        ds.plot.barh(color=c.values, xlabel=label, ax=ax)
+        ax.grid(axis="y")
 
-def rename_index(ds):
+    fig, ax = plt.subplots()
-    return ds.set_axis(ds.index.map(lambda x: f"{x[1]}\n({x[0].lower()})"))
+    ds = n.statistics.capacity_factor().dropna()
+    plot_static_per_carrier(ds, ax)
+    fig.savefig(snakemake.output.capacity_factor_bar)
 
+    fig, ax = plt.subplots()
+    ds = n.statistics.installed_capacity().dropna()
+    ds = ds.drop("Line")
+    ds = ds.drop(("Generator", "Load"))
+    ds = ds / 1e3
+    ds.attrs["unit"] = "GW"
+    plot_static_per_carrier(ds, ax)
+    fig.savefig(snakemake.output.installed_capacity_bar)
 
-def plot_static_per_carrier(ds, ax, drop_zero=True):
+    fig, ax = plt.subplots()
-    if drop_zero:
+    ds = n.statistics.optimal_capacity()
-        ds = ds[ds != 0]
+    ds = ds.drop("Line")
-    ds = ds.dropna()
+    ds = ds.drop(("Generator", "Load"))
-    c = colors[ds.index.get_level_values("carrier")]
+    ds = ds / 1e3
-    ds = ds.pipe(rename_index)
+    ds.attrs["unit"] = "GW"
-    label = f"{ds.attrs['name']} [{ds.attrs['unit']}]"
+    plot_static_per_carrier(ds, ax)
-    ds.plot.barh(color=c.values, xlabel=label, ax=ax)
+    fig.savefig(snakemake.output.optimal_capacity_bar)
 
+    fig, ax = plt.subplots()
+    ds = n.statistics.capex()
+    plot_static_per_carrier(ds, ax)
+    fig.savefig(snakemake.output.capital_expenditure_bar)
 
-fig, ax = plt.subplots()
+    fig, ax = plt.subplots()
-ds = n.statistics.capacity_factor().dropna()
+    ds = n.statistics.opex()
-plot_static_per_carrier(ds, ax)
+    plot_static_per_carrier(ds, ax)
-# fig.savefig("")
+    fig.savefig(snakemake.output.operational_expenditure_bar)
 
-fig, ax = plt.subplots()
+    fig, ax = plt.subplots()
-ds = n.statistics.installed_capacity().dropna()
+    ds = n.statistics.curtailment()
-ds = ds.drop("Line")
+    plot_static_per_carrier(ds, ax)
-ds = ds / 1e3
+    fig.savefig(snakemake.output.curtailment_bar)
-ds.attrs["unit"] = "GW"
-plot_static_per_carrier(ds, ax)
-# fig.savefig("")
 
+    fig, ax = plt.subplots()
+    ds = n.statistics.supply()
+    ds = ds.drop("Line")
+    ds = ds / 1e6
+    ds.attrs["unit"] = "TWh"
+    plot_static_per_carrier(ds, ax)
+    fig.savefig(snakemake.output.supply_bar)
 
-fig, ax = plt.subplots()
+    fig, ax = plt.subplots()
-ds = n.statistics.optimal_capacity()
+    ds = n.statistics.withdrawal()
-ds = ds.drop("Line")
+    ds = ds.drop("Line")
-ds = ds / 1e3
+    ds = ds / -1e6
-ds.attrs["unit"] = "GW"
+    ds.attrs["unit"] = "TWh"
-plot_static_per_carrier(ds, ax)
+    plot_static_per_carrier(ds, ax)
-# fig.savefig("")
+    fig.savefig(snakemake.output.withdrawal_bar)
 
+    fig, ax = plt.subplots()
+    ds = n.statistics.market_value()
+    plot_static_per_carrier(ds, ax)
+    fig.savefig(snakemake.output.market_value_bar)
 
-fig, ax = plt.subplots()
+    # touch file
-ds = n.statistics.capex()
+    with open(snakemake.output.barplots_touch, "a"):
-plot_static_per_carrier(ds, ax)
+        pass
-# fig.savefig("")
-
-fig, ax = plt.subplots()
-ds = n.statistics.curtailment()
-plot_static_per_carrier(ds, ax)
-# fig.savefig("")
-
-fig, ax = plt.subplots()
-ds = n.statistics.supply()
-ds = ds.drop("Line")
-ds = ds / 1e6
-ds.attrs["unit"] = "TWh"
-plot_static_per_carrier(ds, ax)
-# fig.savefig("")
-
-
-fig, ax = plt.subplots()
-ds = n.statistics.withdrawal()
-ds = ds.drop("Line")
-ds = ds / -1e6
-ds.attrs["unit"] = "TWh"
-plot_static_per_carrier(ds, ax)
-# fig.savefig("")