correct the electricity trade calculaton

workflow/Snakefile

@@ -570,6 +570,8 @@ rule plot_ariadne_variables:
         transmission_investment_csv=RESULTS + "ariadne/transmission_investment.csv",
         trassenlaenge_csv=RESULTS + "ariadne/trassenlaenge.csv",
         Kernnetz_Investment_plot=RESULTS + "ariadne/Kernnetz_Investment_plot.png",
+        elec_trade=RESULTS + "ariadne/elec-trade-DE.png",
+        h2_trade=RESULTS + "ariadne/h2-trade-DE.png",
     log:
         RESULTS + "logs/plot_ariadne_variables.log",
     script:

workflow/scripts/export_ariadne_variables.py

@@ -4330,72 +4330,92 @@ def get_trade(n, region):
     var = pd.Series()
 
     def get_export_import_links(n, region, carriers):
-        exporting = n.links.index[
+        # note: links can also used bidirectional if efficiency=1 (e.g. "H2 pipeline retrofitted")
+        outgoing = n.links.index[
             (n.links.carrier.isin(carriers))
             & (n.links.bus0.str[:2] == region)
             & (n.links.bus1.str[:2] != region)
         ]
 
-        importing = n.links.index[
+        incoming = n.links.index[
             (n.links.carrier.isin(carriers))
             & (n.links.bus0.str[:2] != region)
             & (n.links.bus1.str[:2] == region)
         ]
 
         exporting_p = (
-            n.links_t.p0.loc[:, exporting]
+            # if p0 > 0 (=clip(lower=0)) system is withdrawing from bus0 (DE) and feeding into bus1 (non-DE) -> export
+            n.links_t.p0.loc[:, outgoing]
             .clip(lower=0)
             .multiply(n.snapshot_weightings.generators, axis=0)
             .values.sum()
-            - n.links_t.p0.loc[:, importing]
-            .clip(upper=0)
+            +
+            # if p1 > 0 system is withdrawing from bus1 (DE) and feeding into bus0 (non-DE) -> export
+            n.links_t.p1.loc[:, incoming]
+            .clip(lower=0)
             .multiply(n.snapshot_weightings.generators, axis=0)
             .values.sum()
         )
+
         importing_p = (
-            n.links_t.p0.loc[:, importing]
-            .clip(lower=0)
+            # if p1 < 0 (=clip(upper=0)) system is feeding into bus1 (DE) and withdrawing from bus0 (non-DE) -> import (with negative sign here)
+            n.links_t.p1.loc[:, incoming]
+            .clip(upper=0)
             .multiply(n.snapshot_weightings.generators, axis=0)
             .values.sum()
-            - n.links_t.p0.loc[:, exporting]
+            * -1
+            +
+            # if p0 < 0 (=clip(upper=0)) system is feeding into bus0 (DE) and withdrawing from bus1 (non-DE) -> import (with negative sign here)
+            n.links_t.p0.loc[:, outgoing]
             .clip(upper=0)
             .multiply(n.snapshot_weightings.generators, axis=0)
             .values.sum()
+            * -1
         )
 
         return exporting_p, importing_p
 
     # Trade|Secondary Energy|Electricity|Volume
-    exporting_ac = n.lines.index[
+    outgoing_ac = n.lines.index[
         (n.lines.carrier == "AC")
         & (n.lines.bus0.str[:2] == region)
         & (n.lines.bus1.str[:2] != region)
     ]
 
-    importing_ac = n.lines.index[
+    incoming_ac = n.lines.index[
         (n.lines.carrier == "AC")
         & (n.lines.bus0.str[:2] != region)
         & (n.lines.bus1.str[:2] == region)
     ]
+
     exporting_p_ac = (
-        n.lines_t.p0.loc[:, exporting_ac]
+        # if p0 > 0 (=clip(lower=0)) system is withdrawing from bus0 (DE) and feeding into bus1 (non-DE) -> export
+        n.lines_t.p0.loc[:, outgoing_ac]
         .clip(lower=0)
         .multiply(n.snapshot_weightings.generators, axis=0)
         .values.sum()
-        - n.lines_t.p0.loc[:, importing_ac]
-        .clip(upper=0)
+        +
+        # if p1 > 0 system is withdrawing from bus1 (DE) and feeding into bus0 (non-DE) -> export
+        n.lines_t.p1.loc[:, incoming_ac]
+        .clip(lower=0)
         .multiply(n.snapshot_weightings.generators, axis=0)
         .values.sum()
     )
+
     importing_p_ac = (
-        n.lines_t.p0.loc[:, importing_ac]
-        .clip(lower=0)
+        # if p1 < 0 (=clip(upper=0)) system is feeding into bus1 (DE) and withdrawing from bus0 (non-DE) -> import (with negative sign here)
+        n.lines_t.p1.loc[:, incoming_ac]
+        .clip(upper=0)
         .multiply(n.snapshot_weightings.generators, axis=0)
         .values.sum()
-        - n.lines_t.p0.loc[:, exporting_ac]
+        * -1
+        +
+        # if p0 < 0 (=clip(upper=0)) system is feeding into bus0 (DE) and withdrawing from bus1 (non-DE) -> import (with negative sign here)
+        n.lines_t.p0.loc[:, outgoing_ac]
         .clip(upper=0)
         .multiply(n.snapshot_weightings.generators, axis=0)
         .values.sum()
+        * -1
     )
 
     exports_dc, imports_dc = get_export_import_links(n, region, ["DC"])

workflow/scripts/plot_ariadne_report.py

@@ -2020,128 +2020,6 @@ def plot_elec_map_de(
     plt.close()
 
 
-def plot_elec_trade(
-    networks,
-    planning_horizons,
-    tech_colors,
-    savepath,
-):
-    incoming_elec = []
-    outgoing_elec = []
-    for year in planning_horizons:
-        n = networks[planning_horizons.index(year)]
-        incoming_lines = n.lines[
-            (n.lines.bus0.str[:2] != "DE") & (n.lines.bus1.str[:2] == "DE")
-        ].index
-        outgoing_lines = n.lines[
-            (n.lines.bus0.str[:2] == "DE") & (n.lines.bus1.str[:2] != "DE")
-        ].index
-        incoming_links = n.links[
-            (n.links.carrier == "DC")
-            & (n.links.bus0.str[:2] != "DE")
-            & (n.links.bus1.str[:2] == "DE")
-        ].index
-        outgoing_links = n.links[
-            (n.links.carrier == "DC")
-            & (n.links.bus0.str[:2] == "DE")
-            & (n.links.bus1.str[:2] != "DE")
-        ].index
-        # positive when withdrawing power from bus0/bus1
-        incoming = n.lines_t.p1[incoming_lines].sum(axis=1).mul(
-            n.snapshot_weightings.generators, axis=0
-        ) + n.links_t.p1[incoming_links].sum(axis=1).mul(
-            n.snapshot_weightings.generators, axis=0
-        )
-        outgoing = n.lines_t.p0[outgoing_lines].sum(axis=1).mul(
-            n.snapshot_weightings.generators, axis=0
-        ) + n.links_t.p0[outgoing_links].sum(axis=1).mul(
-            n.snapshot_weightings.generators, axis=0
-        )
-        incoming_elec.append(
-            incoming[incoming < 0].abs().sum() + outgoing[outgoing > 0].sum()
-        )
-        outgoing_elec.append(
-            incoming[incoming > 0].sum() + outgoing[outgoing < 0].abs().sum()
-        )
-    elec_import = np.array(incoming_elec) / 1e6
-    elec_export = -np.array(outgoing_elec) / 1e6
-    net = elec_import + elec_export
-    x = np.arange(len(planning_horizons))
-
-    fig, ax = plt.subplots(figsize=(6, 4))
-    ax.bar(x, elec_import, color=tech_colors["AC"], label="Import")
-    ax.bar(x, elec_export, color="#3f630f", label="Export")
-    ax.scatter(x, net, color="black", marker="x", label="Netto")
-
-    ax.set_xticks(x)
-    ax.set_xticklabels(planning_horizons)
-    ax.axhline(0, color="black", linewidth=0.5)
-    ax.set_ylabel("Strom [TWh]")
-    ax.set_title("Strom Import/Export Deutschland")
-    ax.legend()
-
-    plt.tight_layout()
-    fig.savefig(savepath, bbox_inches="tight")
-
-
-def plot_h2_trade(
-    networks,
-    planning_horizons,
-    tech_colors,
-    savepath,
-):
-    incoming_h2 = []
-    outgoing_h2 = []
-    for year in planning_horizons:
-        n = networks[planning_horizons.index(year)]
-        incoming_links = n.links[
-            (n.links.carrier.str.contains("H2 pipeline"))
-            & (n.links.bus0.str[:2] != "DE")
-            & (n.links.bus1.str[:2] == "DE")
-        ].index
-        outgoing_links = n.links[
-            (n.links.carrier.str.contains("H2 pipeline"))
-            & (n.links.bus0.str[:2] == "DE")
-            & (n.links.bus1.str[:2] != "DE")
-        ].index
-        # positive when withdrawing power from bus0/bus1
-        incoming = (
-            n.links_t.p1[incoming_links]
-            .sum(axis=1)
-            .mul(n.snapshot_weightings.generators, axis=0)
-        )
-        outgoing = (
-            n.links_t.p0[outgoing_links]
-            .sum(axis=1)
-            .mul(n.snapshot_weightings.generators, axis=0)
-        )
-        incoming_h2.append(
-            incoming[incoming < 0].abs().sum() + outgoing[outgoing > 0].sum()
-        )
-        outgoing_h2.append(
-            incoming[incoming > 0].sum() + outgoing[outgoing < 0].abs().sum()
-        )
-    h2_import = np.array(incoming_h2) / 1e6
-    h2_export = -np.array(outgoing_h2) / 1e6
-    net = h2_import + h2_export
-    x = np.arange(len(planning_horizons))
-
-    fig, ax = plt.subplots(figsize=(6, 4))
-    ax.bar(x, h2_import, color=tech_colors["H2 pipeline"], label="Import")
-    ax.bar(x, h2_export, color=tech_colors["H2 pipeline (Kernnetz)"], label="Export")
-    ax.scatter(x, net, color="black", marker="x", label="Netto")
-
-    ax.set_xticks(x)
-    ax.set_xticklabels(planning_horizons)
-    ax.axhline(0, color="black", linewidth=0.5)
-    ax.set_ylabel("H2 [TWh]")
-    ax.set_title("Wasserstoff Import/Export Deutschland")
-    ax.legend()
-
-    plt.tight_layout()
-    fig.savefig(savepath, bbox_inches="tight")
-
-
 # electricity capacity map
 def plot_cap_map_de(
     network,
@@ -2623,13 +2501,6 @@ def plot_cap_map_de(
             )
             del network
 
-    plot_h2_trade(
-        networks,
-        planning_horizons,
-        tech_colors,
-        savepath=f"{snakemake.output.h2_transmission}/h2-trade-DE.png",
-    )
-
     ## electricity transmission
     logger.info("Plotting electricity transmission")
     for year in planning_horizons:
@@ -2650,13 +2521,6 @@ def plot_cap_map_de(
             savepath=f"{snakemake.output.elec_transmission}/elec-cap-DE-{year}.png",
         )
 
-    plot_elec_trade(
-        networks,
-        planning_horizons,
-        tech_colors,
-        savepath=f"{snakemake.output.elec_transmission}/elec-trade-DE.png",
-    )
-
     ## nodal balances general (might not be very robust)
     logger.info("Plotting nodal balances")
     plt.style.use(["bmh", snakemake.input.rc])

workflow/scripts/plot_ariadne_variables.py

@@ -5,6 +5,80 @@
 import numpy as np
 import pandas as pd
 
+TWh2PJ = 3.6
+
+
+def plot_h2_trade(
+    df,
+    savepath,
+):
+    # load data and convert to TWh
+    h2_balance = (
+        df.loc["Trade|Secondary Energy|Hydrogen|Volume"] / TWh2PJ * -1
+    )  # exports-imports
+    h2_import = df.loc["Trade|Secondary Energy|Hydrogen|Gross Import|Volume"] / TWh2PJ
+    h2_export = h2_balance - h2_import
+
+    fig, ax = plt.subplots(figsize=(6, 4))
+    ax.bar(h2_import.columns, h2_import.loc["PJ/yr"], color="#f081dc", label="Import")
+    ax.bar(h2_export.columns, h2_export.loc["PJ/yr"], color="#6b3161", label="Export")
+    ax.scatter(
+        h2_balance.columns,
+        h2_balance.loc["PJ/yr"],
+        color="black",
+        marker="x",
+        label="Netto",
+    )
+
+    ax.set_xticks(h2_balance.columns)
+    ax.set_xticklabels(h2_balance.columns)
+    ax.axhline(0, color="black", linewidth=0.5)
+    ax.set_ylabel("TWh")
+    ax.set_title("Wasserstoff Import/Export Deutschland")
+    ax.legend()
+
+    plt.tight_layout()
+    fig.savefig(savepath, bbox_inches="tight")
+
+
+def plot_elec_trade(
+    df,
+    savepath,
+):
+    # load data and convert to TWh
+    elec_balance = (
+        df.loc["Trade|Secondary Energy|Electricity|Volume"] / TWh2PJ * -1
+    )  # exports-imports
+    elec_import = (
+        df.loc["Trade|Secondary Energy|Electricity|Gross Import|Volume"] / TWh2PJ
+    )
+    elec_export = elec_balance - elec_import
+
+    fig, ax = plt.subplots(figsize=(6, 4))
+    ax.bar(
+        elec_import.columns, elec_import.loc["PJ/yr"], color="#70af1d", label="Import"
+    )
+    ax.bar(
+        elec_export.columns, elec_export.loc["PJ/yr"], color="#3f630f", label="Export"
+    )
+    ax.scatter(
+        elec_balance.columns,
+        elec_balance.loc["PJ/yr"],
+        color="black",
+        marker="x",
+        label="Netto",
+    )
+
+    ax.set_xticks(elec_balance.columns)
+    ax.set_xticklabels(elec_balance.columns)
+    ax.axhline(0, color="black", linewidth=0.5)
+    ax.set_ylabel("Twh")
+    ax.set_title("Strom Import/Export Deutschland")
+    ax.legend(loc="upper left")
+
+    plt.tight_layout()
+    fig.savefig(savepath, bbox_inches="tight")
+
 
 def plot_Kernnetz(df, savepath=None, currency_year=2020):
     key = "Investment|Energy Supply|Hydrogen|Transmission and Distribution|"
@@ -678,7 +752,7 @@ def elec_val_plot(df, savepath):
         snakemake = mock_snakemake(
             "plot_ariadne_variables",
             simpl="",
-            clusters=22,
+            clusters=49,
             opts="",
             ll="v1.2",
             sector_opts="None",
@@ -931,3 +1005,13 @@ def elec_val_plot(df, savepath):
     plot_Kernnetz(
         df, savepath=snakemake.output.Kernnetz_Investment_plot, currency_year=2020
     )
+
+    plot_elec_trade(
+        df,
+        savepath=snakemake.output.elec_trade,
+    )
+
+    plot_h2_trade(
+        df,
+        savepath=snakemake.output.h2_trade,
+    )