# -*- coding: utf-8 -*- # SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors # # SPDX-License-Identifier: MIT """ Creates map of optimised hydrogen network, storage and selected other infrastructure. """ import logging import geopandas as gpd import matplotlib.pyplot as plt import pandas as pd import pypsa from _helpers import configure_logging, set_scenario_config from plot_power_network import assign_location, load_projection from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches logger = logging.getLogger(__name__) def group_pipes(df, drop_direction=False): """ Group pipes which connect same buses and return overall capacity. """ df = df.copy() if drop_direction: positive_order = df.bus0 < df.bus1 df_p = df[positive_order] swap_buses = {"bus0": "bus1", "bus1": "bus0"} df_n = df[~positive_order].rename(columns=swap_buses) df = pd.concat([df_p, df_n]) # there are pipes for each investment period rename to AC buses name for plotting df["index_orig"] = df.index df.index = df.apply( lambda x: f"H2 pipeline {x.bus0.replace(' H2', '')} -> {x.bus1.replace(' H2', '')}", axis=1, ) return df.groupby(level=0).agg( {"p_nom_opt": "sum", "bus0": "first", "bus1": "first", "index_orig": "first"} ) def plot_h2_map(n, regions): # if "H2 pipeline" not in n.links.carrier.unique(): # return assign_location(n) h2_storage = n.stores.query("carrier == 'H2'") regions["H2"] = ( h2_storage.rename(index=h2_storage.bus.map(n.buses.location)) .e_nom_opt.groupby(level=0) .sum() .div(1e6) ) # TWh regions["H2"] = regions["H2"].where(regions["H2"] > 0.1) bus_size_factor = 1e5 linewidth_factor = 7e3 # MW below which not drawn line_lower_threshold = 750 # Drop non-electric buses so they don't clutter the plot n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True) carriers = ["H2 Electrolysis", "H2 Fuel Cell"] elec = n.links[n.links.carrier.isin(carriers)].index bus_sizes = ( n.links.loc[elec, "p_nom_opt"].groupby([n.links["bus0"], n.links.carrier]).sum() / bus_size_factor ) # make a fake MultiIndex so that area is correct for legend bus_sizes.rename(index=lambda x: x.replace(" H2", ""), level=0, inplace=True) # drop all links which are not H2 pipelines n.links.drop( n.links.index[~n.links.carrier.str.contains("H2 pipeline")], inplace=True ) h2_new = n.links[n.links.carrier == "H2 pipeline"] h2_retro = n.links[n.links.carrier == "H2 pipeline retrofitted"] if snakemake.params.foresight == "myopic": # sum capacitiy for pipelines from different investment periods h2_new = group_pipes(h2_new) if not h2_retro.empty: h2_retro = ( group_pipes(h2_retro, drop_direction=True) .reindex(h2_new.index) .fillna(0) ) if not h2_retro.empty: if snakemake.params.foresight != "myopic": positive_order = h2_retro.bus0 < h2_retro.bus1 h2_retro_p = h2_retro[positive_order] swap_buses = {"bus0": "bus1", "bus1": "bus0"} h2_retro_n = h2_retro[~positive_order].rename(columns=swap_buses) h2_retro = pd.concat([h2_retro_p, h2_retro_n]) h2_retro["index_orig"] = h2_retro.index h2_retro.index = h2_retro.apply( lambda x: f"H2 pipeline {x.bus0.replace(' H2', '')} -> {x.bus1.replace(' H2', '')}", axis=1, ) retro_w_new_i = h2_retro.index.intersection(h2_new.index) h2_retro_w_new = h2_retro.loc[retro_w_new_i] retro_wo_new_i = h2_retro.index.difference(h2_new.index) h2_retro_wo_new = h2_retro.loc[retro_wo_new_i] h2_retro_wo_new.index = h2_retro_wo_new.index_orig.apply( lambda x: x.split("-2")[0] ) to_concat = [h2_new, h2_retro_w_new, h2_retro_wo_new] h2_total = pd.concat(to_concat).p_nom_opt.groupby(level=0).sum() else: h2_total = h2_new.p_nom_opt link_widths_total = h2_total / linewidth_factor n.links.rename(index=lambda x: x.split("-2")[0], inplace=True) # group links by summing up p_nom values and taking the first value of the rest of the columns other_cols = dict.fromkeys(n.links.columns.drop(["p_nom_opt", "p_nom"]), "first") n.links = n.links.groupby(level=0).agg( {"p_nom_opt": "sum", "p_nom": "sum", **other_cols} ) link_widths_total = link_widths_total.reindex(n.links.index).fillna(0.0) link_widths_total[n.links.p_nom_opt < line_lower_threshold] = 0.0 retro = n.links.p_nom_opt.where( n.links.carrier == "H2 pipeline retrofitted", other=0.0 ) link_widths_retro = retro / linewidth_factor link_widths_retro[n.links.p_nom_opt < line_lower_threshold] = 0.0 n.links.bus0 = n.links.bus0.str.replace(" H2", "") n.links.bus1 = n.links.bus1.str.replace(" H2", "") regions = regions.to_crs(proj.proj4_init) fig, ax = plt.subplots(figsize=(7, 6), subplot_kw={"projection": proj}) color_h2_pipe = "#b3f3f4" color_retrofit = "#499a9c" bus_colors = {"H2 Electrolysis": "#ff29d9", "H2 Fuel Cell": "#805394"} n.plot( geomap=True, bus_sizes=bus_sizes, bus_colors=bus_colors, link_colors=color_h2_pipe, link_widths=link_widths_total, branch_components=["Link"], ax=ax, **map_opts, ) n.plot( geomap=True, bus_sizes=0, link_colors=color_retrofit, link_widths=link_widths_retro, branch_components=["Link"], ax=ax, **map_opts, ) regions.plot( ax=ax, column="H2", cmap="Blues", linewidths=0, legend=True, vmax=6, vmin=0, legend_kwds={ "label": "Hydrogen Storage [TWh]", "shrink": 0.7, "extend": "max", }, ) sizes = [50, 10] labels = [f"{s} GW" for s in sizes] sizes = [s / bus_size_factor * 1e3 for s in sizes] legend_kw = dict( loc="upper left", bbox_to_anchor=(0, 1), labelspacing=0.8, handletextpad=0, frameon=False, ) add_legend_circles( ax, sizes, labels, srid=n.srid, patch_kw=dict(facecolor="lightgrey"), legend_kw=legend_kw, ) sizes = [30, 10] labels = [f"{s} GW" for s in sizes] scale = 1e3 / linewidth_factor sizes = [s * scale for s in sizes] legend_kw = dict( loc="upper left", bbox_to_anchor=(0.23, 1), frameon=False, labelspacing=0.8, handletextpad=1, ) add_legend_lines( ax, sizes, labels, patch_kw=dict(color="lightgrey"), legend_kw=legend_kw, ) colors = [bus_colors[c] for c in carriers] + [color_h2_pipe, color_retrofit] labels = carriers + ["H2 pipeline (total)", "H2 pipeline (repurposed)"] legend_kw = dict( loc="upper left", bbox_to_anchor=(0, 1.13), ncol=2, frameon=False, ) add_legend_patches(ax, colors, labels, legend_kw=legend_kw) ax.set_facecolor("white") fig.savefig(snakemake.output.map, bbox_inches="tight") if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake snakemake = mock_snakemake( "plot_hydrogen_network", opts="", clusters="37", ll="v1.0", sector_opts="4380H-T-H-B-I-A-dist1", ) configure_logging(snakemake) set_scenario_config(snakemake) n = pypsa.Network(snakemake.input.network) regions = gpd.read_file(snakemake.input.regions).set_index("name") map_opts = snakemake.params.plotting["map"] if map_opts["boundaries"] is None: map_opts["boundaries"] = regions.total_bounds[[0, 2, 1, 3]] + [-1, 1, -1, 1] proj = load_projection(snakemake.params.plotting) plot_h2_map(n, regions)