bjoern/pypsa-eur
1
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases Wiki Activity

plot_network: split into separate scripts for power, hydrogen, gas

Browse Source
This commit is contained in:
Fabian Neumann 2024-01-25 20:34:59 +01:00
parent ce4d18c861
commit ffd4e1f1af
9 changed files with 1086 additions and 895 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
doc/plotting.rst
View File

@ -22,7 +22,22 @@ Rule ``plot_summary``
.. _map_plot:
Rule ``plot_network``
========================
Rule ``plot_power_network``
===========================
.. automodule:: plot_network
.. automodule:: plot_power_network
Rule ``plot_power_network_perfect``
===================================
.. automodule:: plot_power_network_perfect
Rule ``plot_hydrogen_network``
==============================
.. automodule:: plot_hydrogen_network
Rule ``plot_gas_network``
=========================
.. automodule:: plot_gas_network

6
doc/release_notes.rst
View File

@ -37,6 +37,12 @@ Upcoming Release
* Add the option to customise map projection in plotting config.
* The rule ``plot_network`` has been split into separate rules for plotting
electricity, hydrogen and gas networks.
* Added new collection rule ``plot_all`` which should be used instead of
``plot_summary``. This allows running the rule :mod:`make_summary` and
:mod:`plot_summary` even if the network plotting rules fail.
PyPSA-Eur 0.9.0 (5th January 2024)
==================================

10
rules/collect.smk
View File

@ -11,7 +11,6 @@ localrules:
prepare_sector_networks,
solve_elec_networks,
solve_sector_networks,
plot_networks,
rule cluster_networks:
@ -69,15 +68,6 @@ rule solve_sector_networks_perfect:
),
rule plot_networks:
input:
expand(
RESULTS
+ "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-costs-all_{planning_horizons}.pdf",
**config["scenario"]
),
rule validate_elec_networks:
input:
expand(

87
rules/postprocess.smk
View File

@ -1,4 +1,4 @@
# SPDX-FileCopyrightText: : 2023 The PyPSA-Eur Authors
# SPDX-FileCopyrightText: : 2023-2024 The PyPSA-Eur Authors
#
# SPDX-License-Identifier: MIT
@ -9,9 +9,8 @@ localrules:
if config["foresight"] != "perfect":
rule plot_network:
rule plot_power_network:
params:
foresight=config["foresight"],
plotting=config["plotting"],
input:
network=RESULTS
@ -26,19 +25,66 @@ if config["foresight"] != "perfect":
benchmark:
(
BENCHMARKS
+ "plot_network/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}"
+ "plot_power_network/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}"
)
conda:
"../envs/environment.yaml"
script:
"../scripts/plot_network.py"
"../scripts/plot_power_network.py"
rule plot_hydrogen_network:
params:
plotting=config["plotting"],
input:
network=RESULTS
+ "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
regions=RESOURCES + "regions_onshore_elec_s{simpl}_{clusters}.geojson",
output:
map=RESULTS
+ "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-h2_network_{planning_horizons}.pdf",
threads: 2
resources:
mem_mb=10000,
benchmark:
(
BENCHMARKS
+ "plot_hydrogen_network/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}"
)
conda:
"../envs/environment.yaml"
script:
"../scripts/plot_hydrogen_network.py"
rule plot_gas_network:
params:
plotting=config["plotting"],
input:
network=RESULTS
+ "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
regions=RESOURCES + "regions_onshore_elec_s{simpl}_{clusters}.geojson",
output:
map=RESULTS
+ "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-ch4_network_{planning_horizons}.pdf",
threads: 2
resources:
mem_mb=10000,
benchmark:
(
BENCHMARKS
+ "plot_gas_network/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}"
)
conda:
"../envs/environment.yaml"
script:
"../scripts/plot_gas_network.py"
if config["foresight"] == "perfect":
rule plot_network:
rule plot_power_network_perfect:
params:
foresight=config["foresight"],
plotting=config["plotting"],
input:
network=RESULTS
@ -60,7 +106,7 @@ if config["foresight"] == "perfect":
conda:
"../envs/environment.yaml"
script:
"../scripts/plot_network.py"
"../scripts/plot_power_network_perfect.py"
rule copy_config:
@ -95,11 +141,6 @@ rule make_summary:
costs="data/costs_{}.csv".format(config["costs"]["year"])
if config["foresight"] == "overnight"
else "data/costs_{}.csv".format(config["scenario"]["planning_horizons"][0]),
plots=expand(
RESULTS
+ "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-costs-all_{planning_horizons}.pdf",
**config["scenario"]
),
output:
nodal_costs=RESULTS + "csvs/nodal_costs.csv",
nodal_capacities=RESULTS + "csvs/nodal_capacities.csv",
@ -161,6 +202,26 @@ rule plot_summary:
"../scripts/plot_summary.py"
rule plot_all:
input:
RESULTS + "graphs/costs.pdf",
expand(
RESULTS
+ "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-costs-all_{planning_horizons}.pdf",
**config["scenario"]
),
expand(
RESULTS
+ "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-h2_network_{planning_horizons}.pdf",
**config["scenario"]
),
expand(
RESULTS
+ "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-ch4_network_{planning_horizons}.pdf",
**config["scenario"]
),
STATISTICS_BARPLOTS = [
"capacity_factor",
"installed_capacity",

251
scripts/plot_gas_network.py Normal file
View File

@ -0,0 +1,251 @@
# -*- coding: utf-8 -*-
# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
#
# SPDX-License-Identifier: MIT
"""
Creates map of optimised gas 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
from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches
from plot_power_network import assign_location, load_projection
logger = logging.getLogger(__name__)
def plot_ch4_map(n):
# if "gas pipeline" not in n.links.carrier.unique():
# return
assign_location(n)
bus_size_factor = 8e7
linewidth_factor = 1e4
# MW below which not drawn
line_lower_threshold = 1e3
# Drop non-electric buses so they don't clutter the plot
n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True)
fossil_gas_i = n.generators[n.generators.carrier == "gas"].index
fossil_gas = (
n.generators_t.p.loc[:, fossil_gas_i]
.mul(n.snapshot_weightings.generators, axis=0)
.sum()
.groupby(n.generators.loc[fossil_gas_i, "bus"])
.sum()
/ bus_size_factor
)
fossil_gas.rename(index=lambda x: x.replace(" gas", ""), inplace=True)
fossil_gas = fossil_gas.reindex(n.buses.index).fillna(0)
# make a fake MultiIndex so that area is correct for legend
fossil_gas.index = pd.MultiIndex.from_product([fossil_gas.index, ["fossil gas"]])
methanation_i = n.links.query("carrier == 'Sabatier'").index
methanation = (
abs(
n.links_t.p1.loc[:, methanation_i].mul(
n.snapshot_weightings.generators, axis=0
)
)
.sum()
.groupby(n.links.loc[methanation_i, "bus1"])
.sum()
/ bus_size_factor
)
methanation = (
methanation.groupby(methanation.index)
.sum()
.rename(index=lambda x: x.replace(" gas", ""))
)
# make a fake MultiIndex so that area is correct for legend
methanation.index = pd.MultiIndex.from_product([methanation.index, ["methanation"]])
biogas_i = n.stores[n.stores.carrier == "biogas"].index
biogas = (
n.stores_t.p.loc[:, biogas_i]
.mul(n.snapshot_weightings.generators, axis=0)
.sum()
.groupby(n.stores.loc[biogas_i, "bus"])
.sum()
/ bus_size_factor
)
biogas = (
biogas.groupby(biogas.index)
.sum()
.rename(index=lambda x: x.replace(" biogas", ""))
)
# make a fake MultiIndex so that area is correct for legend
biogas.index = pd.MultiIndex.from_product([biogas.index, ["biogas"]])
bus_sizes = pd.concat([fossil_gas, methanation, biogas])
bus_sizes.sort_index(inplace=True)
to_remove = n.links.index[~n.links.carrier.str.contains("gas pipeline")]
n.links.drop(to_remove, inplace=True)
link_widths_rem = n.links.p_nom_opt / linewidth_factor
link_widths_rem[n.links.p_nom_opt < line_lower_threshold] = 0.0
link_widths_orig = n.links.p_nom / linewidth_factor
link_widths_orig[n.links.p_nom < line_lower_threshold] = 0.0
max_usage = n.links_t.p0.abs().max(axis=0)
link_widths_used = max_usage / linewidth_factor
link_widths_used[max_usage < line_lower_threshold] = 0.0
tech_colors = snakemake.params.plotting["tech_colors"]
pipe_colors = {
"gas pipeline": "#f08080",
"gas pipeline new": "#c46868",
"gas pipeline (in 2020)": "lightgrey",
"gas pipeline (available)": "#e8d1d1",
}
link_color_used = n.links.carrier.map(pipe_colors)
n.links.bus0 = n.links.bus0.str.replace(" gas", "")
n.links.bus1 = n.links.bus1.str.replace(" gas", "")
bus_colors = {
"fossil gas": tech_colors["fossil gas"],
"methanation": tech_colors["methanation"],
"biogas": "seagreen",
}
fig, ax = plt.subplots(figsize=(7, 6), subplot_kw={"projection": proj})
n.plot(
bus_sizes=bus_sizes,
bus_colors=bus_colors,
link_colors=pipe_colors["gas pipeline (in 2020)"],
link_widths=link_widths_orig,
branch_components=["Link"],
ax=ax,
**map_opts,
)
n.plot(
ax=ax,
bus_sizes=0.0,
link_colors=pipe_colors["gas pipeline (available)"],
link_widths=link_widths_rem,
branch_components=["Link"],
color_geomap=False,
boundaries=map_opts["boundaries"],
)
n.plot(
ax=ax,
bus_sizes=0.0,
link_colors=link_color_used,
link_widths=link_widths_used,
branch_components=["Link"],
color_geomap=False,
boundaries=map_opts["boundaries"],
)
sizes = [100, 10]
labels = [f"{s} TWh" for s in sizes]
sizes = [s / bus_size_factor * 1e6 for s in sizes]
legend_kw = dict(
loc="upper left",
bbox_to_anchor=(0, 1.03),
labelspacing=0.8,
frameon=False,
handletextpad=1,
title="gas sources",
)
add_legend_circles(
ax,
sizes,
labels,
srid=n.srid,
patch_kw=dict(facecolor="lightgrey"),
legend_kw=legend_kw,
)
sizes = [50, 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.25, 1.03),
frameon=False,
labelspacing=0.8,
handletextpad=1,
title="gas pipeline",
)
add_legend_lines(
ax,
sizes,
labels,
patch_kw=dict(color="lightgrey"),
legend_kw=legend_kw,
)
colors = list(pipe_colors.values()) + list(bus_colors.values())
labels = list(pipe_colors.keys()) + list(bus_colors.keys())
# legend on the side
# legend_kw = dict(
# bbox_to_anchor=(1.47, 1.04),
# frameon=False,
# )
legend_kw = dict(
loc="upper left",
bbox_to_anchor=(0, 1.24),
ncol=2,
frameon=False,
)
add_legend_patches(
ax,
colors,
labels,
legend_kw=legend_kw,
)
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_gas_network",
simpl="",
opts="",
clusters="37",
ll="v1.0",
sector_opts="4380H-T-H-B-I-A-dist1",
)
configure_logging(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_ch4_map(n)

267
scripts/plot_hydrogen_network.py Normal file
View File

@ -0,0 +1,267 @@
# -*- 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
from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches
from plot_power_network import assign_location, load_projection
logger = logging.getLogger(__name__)
def group_pipes(df, drop_direction=False):
"""
Group pipes which connect same buses and return overall capacity.
"""
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 = 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"})
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:
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
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)
n.links = n.links.groupby(level=0).first()
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",
simpl="",
opts="",
clusters="37",
ll="v1.0",
sector_opts="4380H-T-H-B-I-A-dist1",
)
configure_logging(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)

869
scripts/plot_network.py
View File

@ -1,869 +0,0 @@
# -*- coding: utf-8 -*-
# SPDX-FileCopyrightText: : 2020-2023 The PyPSA-Eur Authors
#
# SPDX-License-Identifier: MIT
"""
Creates plots for optimised network topologies, including electricity, gas and
hydrogen networks, and regional generation, storage and conversion capacities
built.
This rule plots a map of the network with technology capacities at the
nodes.
"""
import logging
import cartopy.crs as ccrs
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import pypsa
from make_summary import assign_carriers
from plot_summary import preferred_order, rename_techs
from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches
logger = logging.getLogger(__name__)
plt.style.use(["ggplot"])
def rename_techs_tyndp(tech):
tech = rename_techs(tech)
if "heat pump" in tech or "resistive heater" in tech:
return "power-to-heat"
elif tech in ["H2 Electrolysis", "methanation", "H2 liquefaction"]:
return "power-to-gas"
elif tech == "H2":
return "H2 storage"
elif tech in ["NH3", "Haber-Bosch", "ammonia cracker", "ammonia store"]:
return "ammonia"
elif tech in ["OCGT", "CHP", "gas boiler", "H2 Fuel Cell"]:
return "gas-to-power/heat"
# elif "solar" in tech:
# return "solar"
elif tech in ["Fischer-Tropsch", "methanolisation"]:
return "power-to-liquid"
elif "offshore wind" in tech:
return "offshore wind"
elif "CC" in tech or "sequestration" in tech:
return "CCS"
else:
return tech
def assign_location(n):
for c in n.iterate_components(n.one_port_components | n.branch_components):
ifind = pd.Series(c.df.index.str.find(" ", start=4), c.df.index)
for i in ifind.value_counts().index:
# these have already been assigned defaults
if i == -1:
continue
names = ifind.index[ifind == i]
c.df.loc[names, "location"] = names.str[:i]
def plot_map(
network,
components=["links", "stores", "storage_units", "generators"],
bus_size_factor=1.7e10,
transmission=False,
with_legend=True,
):
tech_colors = snakemake.params.plotting["tech_colors"]
n = network.copy()
assign_location(n)
# Drop non-electric buses so they don't clutter the plot
n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True)
costs = pd.DataFrame(index=n.buses.index)
for comp in components:
df_c = getattr(n, comp)
if df_c.empty:
continue
df_c["nice_group"] = df_c.carrier.map(rename_techs_tyndp)
attr = "e_nom_opt" if comp == "stores" else "p_nom_opt"
costs_c = (
(df_c.capital_cost * df_c[attr])
.groupby([df_c.location, df_c.nice_group])
.sum()
.unstack()
.fillna(0.0)
)
costs = pd.concat([costs, costs_c], axis=1)
logger.debug(f"{comp}, {costs}")
costs = costs.groupby(costs.columns, axis=1).sum()
costs.drop(list(costs.columns[(costs == 0.0).all()]), axis=1, inplace=True)
new_columns = preferred_order.intersection(costs.columns).append(
costs.columns.difference(preferred_order)
)
costs = costs[new_columns]
for item in new_columns:
if item not in tech_colors:
logger.warning(f"{item} not in config/plotting/tech_colors")
costs = costs.stack() # .sort_index()
# hack because impossible to drop buses...
eu_location = snakemake.params.plotting.get("eu_node_location", dict(x=-5.5, y=46))
n.buses.loc["EU gas", "x"] = eu_location["x"]
n.buses.loc["EU gas", "y"] = eu_location["y"]
n.links.drop(
n.links.index[(n.links.carrier != "DC") & (n.links.carrier != "B2B")],
inplace=True,
)
# drop non-bus
to_drop = costs.index.levels[0].symmetric_difference(n.buses.index)
if len(to_drop) != 0:
logger.info(f"Dropping non-buses {to_drop.tolist()}")
costs.drop(to_drop, level=0, inplace=True, axis=0, errors="ignore")
# make sure they are removed from index
costs.index = pd.MultiIndex.from_tuples(costs.index.values)
threshold = 100e6 # 100 mEUR/a
carriers = costs.groupby(level=1).sum()
carriers = carriers.where(carriers > threshold).dropna()
carriers = list(carriers.index)
# PDF has minimum width, so set these to zero
line_lower_threshold = 500.0
line_upper_threshold = 1e4
linewidth_factor = 4e3
ac_color = "rosybrown"
dc_color = "darkseagreen"
title = "added grid"
if snakemake.wildcards["ll"] == "v1.0":
# should be zero
line_widths = n.lines.s_nom_opt - n.lines.s_nom
link_widths = n.links.p_nom_opt - n.links.p_nom
if transmission:
line_widths = n.lines.s_nom_opt
link_widths = n.links.p_nom_opt
linewidth_factor = 2e3
line_lower_threshold = 0.0
title = "current grid"
else:
line_widths = n.lines.s_nom_opt - n.lines.s_nom_min
link_widths = n.links.p_nom_opt - n.links.p_nom_min
if transmission:
line_widths = n.lines.s_nom_opt
link_widths = n.links.p_nom_opt
title = "total grid"
line_widths = line_widths.clip(line_lower_threshold, line_upper_threshold)
link_widths = link_widths.clip(line_lower_threshold, line_upper_threshold)
line_widths = line_widths.replace(line_lower_threshold, 0)
link_widths = link_widths.replace(line_lower_threshold, 0)
fig, ax = plt.subplots(subplot_kw={"projection": proj})
fig.set_size_inches(7, 6)
n.plot(
bus_sizes=costs / bus_size_factor,
bus_colors=tech_colors,
line_colors=ac_color,
link_colors=dc_color,
line_widths=line_widths / linewidth_factor,
link_widths=link_widths / linewidth_factor,
ax=ax,
**map_opts,
)
sizes = [20, 10, 5]
labels = [f"{s} bEUR/a" for s in sizes]
sizes = [s / bus_size_factor * 1e9 for s in sizes]
legend_kw = dict(
loc="upper left",
bbox_to_anchor=(0.01, 1.06),
labelspacing=0.8,
frameon=False,
handletextpad=0,
title="system cost",
)
add_legend_circles(
ax,
sizes,
labels,
srid=n.srid,
patch_kw=dict(facecolor="lightgrey"),
legend_kw=legend_kw,
)
sizes = [10, 5]
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.27, 1.06),
frameon=False,
labelspacing=0.8,
handletextpad=1,
title=title,
)
add_legend_lines(
ax, sizes, labels, patch_kw=dict(color="lightgrey"), legend_kw=legend_kw
)
legend_kw = dict(
bbox_to_anchor=(1.52, 1.04),
frameon=False,
)
if with_legend:
colors = [tech_colors[c] for c in carriers] + [ac_color, dc_color]
labels = carriers + ["HVAC line", "HVDC link"]
add_legend_patches(
ax,
colors,
labels,
legend_kw=legend_kw,
)
fig.savefig(snakemake.output.map, transparent=True, bbox_inches="tight")
def group_pipes(df, drop_direction=False):
"""
Group pipes which connect same buses and return overall capacity.
"""
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 = 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"})
def plot_h2_map(network, regions):
n = network.copy()
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:
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
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)
n.links = n.links.groupby(level=0).first()
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.replace("-costs-all", "-h2_network"), bbox_inches="tight"
)
def plot_ch4_map(network):
n = network.copy()
if "gas pipeline" not in n.links.carrier.unique():
return
assign_location(n)
bus_size_factor = 8e7
linewidth_factor = 1e4
# MW below which not drawn
line_lower_threshold = 1e3
# Drop non-electric buses so they don't clutter the plot
n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True)
fossil_gas_i = n.generators[n.generators.carrier == "gas"].index
fossil_gas = (
n.generators_t.p.loc[:, fossil_gas_i]
.mul(n.snapshot_weightings.generators, axis=0)
.sum()
.groupby(n.generators.loc[fossil_gas_i, "bus"])
.sum()
/ bus_size_factor
)
fossil_gas.rename(index=lambda x: x.replace(" gas", ""), inplace=True)
fossil_gas = fossil_gas.reindex(n.buses.index).fillna(0)
# make a fake MultiIndex so that area is correct for legend
fossil_gas.index = pd.MultiIndex.from_product([fossil_gas.index, ["fossil gas"]])
methanation_i = n.links.query("carrier == 'Sabatier'").index
methanation = (
abs(
n.links_t.p1.loc[:, methanation_i].mul(
n.snapshot_weightings.generators, axis=0
)
)
.sum()
.groupby(n.links.loc[methanation_i, "bus1"])
.sum()
/ bus_size_factor
)
methanation = (
methanation.groupby(methanation.index)
.sum()
.rename(index=lambda x: x.replace(" gas", ""))
)
# make a fake MultiIndex so that area is correct for legend
methanation.index = pd.MultiIndex.from_product([methanation.index, ["methanation"]])
biogas_i = n.stores[n.stores.carrier == "biogas"].index
biogas = (
n.stores_t.p.loc[:, biogas_i]
.mul(n.snapshot_weightings.generators, axis=0)
.sum()
.groupby(n.stores.loc[biogas_i, "bus"])
.sum()
/ bus_size_factor
)
biogas = (
biogas.groupby(biogas.index)
.sum()
.rename(index=lambda x: x.replace(" biogas", ""))
)
# make a fake MultiIndex so that area is correct for legend
biogas.index = pd.MultiIndex.from_product([biogas.index, ["biogas"]])
bus_sizes = pd.concat([fossil_gas, methanation, biogas])
bus_sizes.sort_index(inplace=True)
to_remove = n.links.index[~n.links.carrier.str.contains("gas pipeline")]
n.links.drop(to_remove, inplace=True)
link_widths_rem = n.links.p_nom_opt / linewidth_factor
link_widths_rem[n.links.p_nom_opt < line_lower_threshold] = 0.0
link_widths_orig = n.links.p_nom / linewidth_factor
link_widths_orig[n.links.p_nom < line_lower_threshold] = 0.0
max_usage = n.links_t.p0.abs().max(axis=0)
link_widths_used = max_usage / linewidth_factor
link_widths_used[max_usage < line_lower_threshold] = 0.0
tech_colors = snakemake.params.plotting["tech_colors"]
pipe_colors = {
"gas pipeline": "#f08080",
"gas pipeline new": "#c46868",
"gas pipeline (in 2020)": "lightgrey",
"gas pipeline (available)": "#e8d1d1",
}
link_color_used = n.links.carrier.map(pipe_colors)
n.links.bus0 = n.links.bus0.str.replace(" gas", "")
n.links.bus1 = n.links.bus1.str.replace(" gas", "")
bus_colors = {
"fossil gas": tech_colors["fossil gas"],
"methanation": tech_colors["methanation"],
"biogas": "seagreen",
}
fig, ax = plt.subplots(figsize=(7, 6), subplot_kw={"projection": proj})
n.plot(
bus_sizes=bus_sizes,
bus_colors=bus_colors,
link_colors=pipe_colors["gas pipeline (in 2020)"],
link_widths=link_widths_orig,
branch_components=["Link"],
ax=ax,
**map_opts,
)
n.plot(
ax=ax,
bus_sizes=0.0,
link_colors=pipe_colors["gas pipeline (available)"],
link_widths=link_widths_rem,
branch_components=["Link"],
color_geomap=False,
boundaries=map_opts["boundaries"],
)
n.plot(
ax=ax,
bus_sizes=0.0,
link_colors=link_color_used,
link_widths=link_widths_used,
branch_components=["Link"],
color_geomap=False,
boundaries=map_opts["boundaries"],
)
sizes = [100, 10]
labels = [f"{s} TWh" for s in sizes]
sizes = [s / bus_size_factor * 1e6 for s in sizes]
legend_kw = dict(
loc="upper left",
bbox_to_anchor=(0, 1.03),
labelspacing=0.8,
frameon=False,
handletextpad=1,
title="gas sources",
)
add_legend_circles(
ax,
sizes,
labels,
srid=n.srid,
patch_kw=dict(facecolor="lightgrey"),
legend_kw=legend_kw,
)
sizes = [50, 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.25, 1.03),
frameon=False,
labelspacing=0.8,
handletextpad=1,
title="gas pipeline",
)
add_legend_lines(
ax,
sizes,
labels,
patch_kw=dict(color="lightgrey"),
legend_kw=legend_kw,
)
colors = list(pipe_colors.values()) + list(bus_colors.values())
labels = list(pipe_colors.keys()) + list(bus_colors.keys())
# legend on the side
# legend_kw = dict(
# bbox_to_anchor=(1.47, 1.04),
# frameon=False,
# )
legend_kw = dict(
loc="upper left",
bbox_to_anchor=(0, 1.24),
ncol=2,
frameon=False,
)
add_legend_patches(
ax,
colors,
labels,
legend_kw=legend_kw,
)
fig.savefig(
snakemake.output.map.replace("-costs-all", "-ch4_network"), bbox_inches="tight"
)
def plot_map_perfect(
network,
components=["Link", "Store", "StorageUnit", "Generator"],
bus_size_factor=1.7e10,
):
n = network.copy()
assign_location(n)
# Drop non-electric buses so they don't clutter the plot
n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True)
# investment periods
investments = n.snapshots.levels[0]
costs = {}
for comp in components:
df_c = n.df(comp)
if df_c.empty:
continue
df_c["nice_group"] = df_c.carrier.map(rename_techs_tyndp)
attr = "e_nom_opt" if comp == "Store" else "p_nom_opt"
active = pd.concat(
[n.get_active_assets(comp, inv_p).rename(inv_p) for inv_p in investments],
axis=1,
).astype(int)
capital_cost = n.df(comp)[attr] * n.df(comp).capital_cost
capital_cost_t = (
(active.mul(capital_cost, axis=0))
.groupby([n.df(comp).location, n.df(comp).nice_group])
.sum()
)
capital_cost_t.drop("load", level=1, inplace=True, errors="ignore")
costs[comp] = capital_cost_t
costs = pd.concat(costs).groupby(level=[1, 2]).sum()
costs.drop(costs[costs.sum(axis=1) == 0].index, inplace=True)
new_columns = preferred_order.intersection(costs.index.levels[1]).append(
costs.index.levels[1].difference(preferred_order)
)
costs = costs.reindex(new_columns, level=1)
for item in new_columns:
if item not in snakemake.config["plotting"]["tech_colors"]:
print(
"Warning!",
item,
"not in config/plotting/tech_colors, assign random color",
)
snakemake.config["plotting"]["tech_colors"] = "pink"
n.links.drop(
n.links.index[(n.links.carrier != "DC") & (n.links.carrier != "B2B")],
inplace=True,
)
# drop non-bus
to_drop = costs.index.levels[0].symmetric_difference(n.buses.index)
if len(to_drop) != 0:
print("dropping non-buses", to_drop)
costs.drop(to_drop, level=0, inplace=True, axis=0, errors="ignore")
# make sure they are removed from index
costs.index = pd.MultiIndex.from_tuples(costs.index.values)
# PDF has minimum width, so set these to zero
line_lower_threshold = 500.0
line_upper_threshold = 1e4
linewidth_factor = 2e3
ac_color = "gray"
dc_color = "m"
line_widths = n.lines.s_nom_opt
link_widths = n.links.p_nom_opt
linewidth_factor = 2e3
line_lower_threshold = 0.0
title = "Today's transmission"
line_widths[line_widths < line_lower_threshold] = 0.0
link_widths[link_widths < line_lower_threshold] = 0.0
line_widths[line_widths > line_upper_threshold] = line_upper_threshold
link_widths[link_widths > line_upper_threshold] = line_upper_threshold
for year in costs.columns:
fig, ax = plt.subplots(subplot_kw={"projection": proj})
fig.set_size_inches(7, 6)
fig.suptitle(year)
n.plot(
bus_sizes=costs[year] / bus_size_factor,
bus_colors=snakemake.config["plotting"]["tech_colors"],
line_colors=ac_color,
link_colors=dc_color,
line_widths=line_widths / linewidth_factor,
link_widths=link_widths / linewidth_factor,
ax=ax,
**map_opts,
)
sizes = [20, 10, 5]
labels = [f"{s} bEUR/a" for s in sizes]
sizes = [s / bus_size_factor * 1e9 for s in sizes]
legend_kw = dict(
loc="upper left",
bbox_to_anchor=(0.01, 1.06),
labelspacing=0.8,
frameon=False,
handletextpad=0,
title="system cost",
)
add_legend_circles(
ax,
sizes,
labels,
srid=n.srid,
patch_kw=dict(facecolor="lightgrey"),
legend_kw=legend_kw,
)
sizes = [10, 5]
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.27, 1.06),
frameon=False,
labelspacing=0.8,
handletextpad=1,
title=title,
)
add_legend_lines(
ax, sizes, labels, patch_kw=dict(color="lightgrey"), legend_kw=legend_kw
)
legend_kw = dict(
bbox_to_anchor=(1.52, 1.04),
frameon=False,
)
fig.savefig(
snakemake.output[f"map_{year}"], transparent=True, bbox_inches="tight"
)
if __name__ == "__main__":
if "snakemake" not in globals():
from _helpers import mock_snakemake
snakemake = mock_snakemake(
"plot_network",
simpl="",
opts="",
clusters="37",
ll="v1.0",
sector_opts="4380H-T-H-B-I-A-dist1",
)
logging.basicConfig(level=snakemake.config["logging"]["level"])
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_kwargs = snakemake.params.plotting.get("projection", dict(name="EqualEarth"))
proj_func = getattr(ccrs, proj_kwargs.pop("name"))
proj = proj_func(**proj_kwargs)
if snakemake.params["foresight"] == "perfect":
plot_map_perfect(
n,
components=["Link", "Store", "StorageUnit", "Generator"],
bus_size_factor=2e10,
)
else:
plot_map(
n,
components=["generators", "links", "stores", "storage_units"],
bus_size_factor=2e10,
transmission=False,
)
plot_h2_map(n, regions)
plot_ch4_map(n)

271
scripts/plot_power_network.py Normal file
View File

@ -0,0 +1,271 @@
# -*- coding: utf-8 -*-
# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
#
# SPDX-License-Identifier: MIT
"""
Creates plots for optimised power network topologies and regional generation,
storage and conversion capacities built.
"""
import logging
import cartopy.crs as ccrs
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import pypsa
from plot_summary import preferred_order, rename_techs
from _helpers import configure_logging
from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches
logger = logging.getLogger(__name__)
def rename_techs_tyndp(tech):
tech = rename_techs(tech)
if "heat pump" in tech or "resistive heater" in tech:
return "power-to-heat"
elif tech in ["H2 Electrolysis", "methanation", "H2 liquefaction"]:
return "power-to-gas"
elif tech == "H2":
return "H2 storage"
elif tech in ["NH3", "Haber-Bosch", "ammonia cracker", "ammonia store"]:
return "ammonia"
elif tech in ["OCGT", "CHP", "gas boiler", "H2 Fuel Cell"]:
return "gas-to-power/heat"
# elif "solar" in tech:
# return "solar"
elif tech in ["Fischer-Tropsch", "methanolisation"]:
return "power-to-liquid"
elif "offshore wind" in tech:
return "offshore wind"
elif "CC" in tech or "sequestration" in tech:
return "CCS"
else:
return tech
def assign_location(n):
for c in n.iterate_components(n.one_port_components | n.branch_components):
ifind = pd.Series(c.df.index.str.find(" ", start=4), c.df.index)
for i in ifind.value_counts().index:
# these have already been assigned defaults
if i == -1:
continue
names = ifind.index[ifind == i]
c.df.loc[names, "location"] = names.str[:i]
def load_projection(plotting_params):
proj_kwargs = plotting_params.get("projection", dict(name="EqualEarth"))
proj_func = getattr(ccrs, proj_kwargs.pop("name"))
return proj_func(**proj_kwargs)
def plot_map(
n,
components=["links", "stores", "storage_units", "generators"],
bus_size_factor=2e10,
transmission=False,
with_legend=True,
):
tech_colors = snakemake.params.plotting["tech_colors"]
assign_location(n)
# Drop non-electric buses so they don't clutter the plot
n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True)
costs = pd.DataFrame(index=n.buses.index)
for comp in components:
df_c = getattr(n, comp)
if df_c.empty:
continue
df_c["nice_group"] = df_c.carrier.map(rename_techs_tyndp)
attr = "e_nom_opt" if comp == "stores" else "p_nom_opt"
costs_c = (
(df_c.capital_cost * df_c[attr])
.groupby([df_c.location, df_c.nice_group])
.sum()
.unstack()
.fillna(0.0)
)
costs = pd.concat([costs, costs_c], axis=1)
logger.debug(f"{comp}, {costs}")
costs = costs.groupby(costs.columns, axis=1).sum()
costs.drop(list(costs.columns[(costs == 0.0).all()]), axis=1, inplace=True)
new_columns = preferred_order.intersection(costs.columns).append(
costs.columns.difference(preferred_order)
)
costs = costs[new_columns]
for item in new_columns:
if item not in tech_colors:
logger.warning(f"{item} not in config/plotting/tech_colors")
costs = costs.stack() # .sort_index()
# hack because impossible to drop buses...
eu_location = snakemake.params.plotting.get("eu_node_location", dict(x=-5.5, y=46))
n.buses.loc["EU gas", "x"] = eu_location["x"]
n.buses.loc["EU gas", "y"] = eu_location["y"]
n.links.drop(
n.links.index[(n.links.carrier != "DC") & (n.links.carrier != "B2B")],
inplace=True,
)
# drop non-bus
to_drop = costs.index.levels[0].symmetric_difference(n.buses.index)
if len(to_drop) != 0:
logger.info(f"Dropping non-buses {to_drop.tolist()}")
costs.drop(to_drop, level=0, inplace=True, axis=0, errors="ignore")
# make sure they are removed from index
costs.index = pd.MultiIndex.from_tuples(costs.index.values)
threshold = 100e6 # 100 mEUR/a
carriers = costs.groupby(level=1).sum()
carriers = carriers.where(carriers > threshold).dropna()
carriers = list(carriers.index)
# PDF has minimum width, so set these to zero
line_lower_threshold = 500.0
line_upper_threshold = 1e4
linewidth_factor = 4e3
ac_color = "rosybrown"
dc_color = "darkseagreen"
title = "added grid"
if snakemake.wildcards["ll"] == "v1.0":
# should be zero
line_widths = n.lines.s_nom_opt - n.lines.s_nom
link_widths = n.links.p_nom_opt - n.links.p_nom
if transmission:
line_widths = n.lines.s_nom_opt
link_widths = n.links.p_nom_opt
linewidth_factor = 2e3
line_lower_threshold = 0.0
title = "current grid"
else:
line_widths = n.lines.s_nom_opt - n.lines.s_nom_min
link_widths = n.links.p_nom_opt - n.links.p_nom_min
if transmission:
line_widths = n.lines.s_nom_opt
link_widths = n.links.p_nom_opt
title = "total grid"
line_widths = line_widths.clip(line_lower_threshold, line_upper_threshold)
link_widths = link_widths.clip(line_lower_threshold, line_upper_threshold)
line_widths = line_widths.replace(line_lower_threshold, 0)
link_widths = link_widths.replace(line_lower_threshold, 0)
fig, ax = plt.subplots(subplot_kw={"projection": proj})
fig.set_size_inches(7, 6)
n.plot(
bus_sizes=costs / bus_size_factor,
bus_colors=tech_colors,
line_colors=ac_color,
link_colors=dc_color,
line_widths=line_widths / linewidth_factor,
link_widths=link_widths / linewidth_factor,
ax=ax,
**map_opts,
)
sizes = [20, 10, 5]
labels = [f"{s} bEUR/a" for s in sizes]
sizes = [s / bus_size_factor * 1e9 for s in sizes]
legend_kw = dict(
loc="upper left",
bbox_to_anchor=(0.01, 1.06),
labelspacing=0.8,
frameon=False,
handletextpad=0,
title="system cost",
)
add_legend_circles(
ax,
sizes,
labels,
srid=n.srid,
patch_kw=dict(facecolor="lightgrey"),
legend_kw=legend_kw,
)
sizes = [10, 5]
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.27, 1.06),
frameon=False,
labelspacing=0.8,
handletextpad=1,
title=title,
)
add_legend_lines(
ax, sizes, labels, patch_kw=dict(color="lightgrey"), legend_kw=legend_kw
)
legend_kw = dict(
bbox_to_anchor=(1.52, 1.04),
frameon=False,
)
if with_legend:
colors = [tech_colors[c] for c in carriers] + [ac_color, dc_color]
labels = carriers + ["HVAC line", "HVDC link"]
add_legend_patches(
ax,
colors,
labels,
legend_kw=legend_kw,
)
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_power_network",
simpl="",
opts="",
clusters="37",
ll="v1.0",
sector_opts="4380H-T-H-B-I-A-dist1",
)
configure_logging(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_map(n)

199
scripts/plot_power_network_perfect.py Normal file
View File

@ -0,0 +1,199 @@
# -*- coding: utf-8 -*-
# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
#
# SPDX-License-Identifier: MIT
"""
Creates plots for optimised power network topologies and regional generation,
storage and conversion capacities built for the perfect foresight scenario.
"""
import logging
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import pypsa
from _helpers import configure_logging
from pypsa.plot import add_legend_circles, add_legend_lines
from plot_power_network import assign_location, rename_techs_tyndp, load_projection
from plot_summary import preferred_order
logger = logging.getLogger(__name__)
def plot_map_perfect(
n,
components=["Link", "Store", "StorageUnit", "Generator"],
bus_size_factor=2e10,
):
assign_location(n)
# Drop non-electric buses so they don't clutter the plot
n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True)
# investment periods
investments = n.snapshots.levels[0]
costs = {}
for comp in components:
df_c = n.df(comp)
if df_c.empty:
continue
df_c["nice_group"] = df_c.carrier.map(rename_techs_tyndp)
attr = "e_nom_opt" if comp == "Store" else "p_nom_opt"
active = pd.concat(
[n.get_active_assets(comp, inv_p).rename(inv_p) for inv_p in investments],
axis=1,
).astype(int)
capital_cost = n.df(comp)[attr] * n.df(comp).capital_cost
capital_cost_t = (
(active.mul(capital_cost, axis=0))
.groupby([n.df(comp).location, n.df(comp).nice_group])
.sum()
)
capital_cost_t.drop("load", level=1, inplace=True, errors="ignore")
costs[comp] = capital_cost_t
costs = pd.concat(costs).groupby(level=[1, 2]).sum()
costs.drop(costs[costs.sum(axis=1) == 0].index, inplace=True)
new_columns = preferred_order.intersection(costs.index.levels[1]).append(
costs.index.levels[1].difference(preferred_order)
)
costs = costs.reindex(new_columns, level=1)
for item in new_columns:
if item not in snakemake.config["plotting"]["tech_colors"]:
print(
"Warning!",
item,
"not in config/plotting/tech_colors, assign random color",
)
snakemake.config["plotting"]["tech_colors"] = "pink"
n.links.drop(
n.links.index[(n.links.carrier != "DC") & (n.links.carrier != "B2B")],
inplace=True,
)
# drop non-bus
to_drop = costs.index.levels[0].symmetric_difference(n.buses.index)
if len(to_drop) != 0:
print("dropping non-buses", to_drop)
costs.drop(to_drop, level=0, inplace=True, axis=0, errors="ignore")
# make sure they are removed from index
costs.index = pd.MultiIndex.from_tuples(costs.index.values)
# PDF has minimum width, so set these to zero
line_lower_threshold = 500.0
line_upper_threshold = 1e4
linewidth_factor = 2e3
ac_color = "gray"
dc_color = "m"
line_widths = n.lines.s_nom_opt
link_widths = n.links.p_nom_opt
linewidth_factor = 2e3
line_lower_threshold = 0.0
title = "Today's transmission"
line_widths[line_widths < line_lower_threshold] = 0.0
link_widths[link_widths < line_lower_threshold] = 0.0
line_widths[line_widths > line_upper_threshold] = line_upper_threshold
link_widths[link_widths > line_upper_threshold] = line_upper_threshold
for year in costs.columns:
fig, ax = plt.subplots(subplot_kw={"projection": proj})
fig.set_size_inches(7, 6)
fig.suptitle(year)
n.plot(
bus_sizes=costs[year] / bus_size_factor,
bus_colors=snakemake.config["plotting"]["tech_colors"],
line_colors=ac_color,
link_colors=dc_color,
line_widths=line_widths / linewidth_factor,
link_widths=link_widths / linewidth_factor,
ax=ax,
**map_opts,
)
sizes = [20, 10, 5]
labels = [f"{s} bEUR/a" for s in sizes]
sizes = [s / bus_size_factor * 1e9 for s in sizes]
legend_kw = dict(
loc="upper left",
bbox_to_anchor=(0.01, 1.06),
labelspacing=0.8,
frameon=False,
handletextpad=0,
title="system cost",
)
add_legend_circles(
ax,
sizes,
labels,
srid=n.srid,
patch_kw=dict(facecolor="lightgrey"),
legend_kw=legend_kw,
)
sizes = [10, 5]
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.27, 1.06),
frameon=False,
labelspacing=0.8,
handletextpad=1,
title=title,
)
add_legend_lines(
ax, sizes, labels, patch_kw=dict(color="lightgrey"), legend_kw=legend_kw
)
legend_kw = dict(
bbox_to_anchor=(1.52, 1.04),
frameon=False,
)
fig.savefig(snakemake.output[f"map_{year}"], bbox_inches="tight")
if __name__ == "__main__":
if "snakemake" not in globals():
from _helpers import mock_snakemake
snakemake = mock_snakemake(
"plot_power_network_perfect",
simpl="",
opts="",
clusters="37",
ll="v1.0",
sector_opts="4380H-T-H-B-I-A-dist1",
)
configure_logging(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_map_perfect(n)