Maps (Static)¶

PyPSA has a built-in method to map parameters and results in a static plot using n.plot(). For interactive maps, please see Interactive Maps

Input data¶

In [2]:

import geopandas as gpd

import pypsa

n = pypsa.examples.scigrid_de()

import geopandas as gpd import pypsa n = pypsa.examples.scigrid_de()

INFO:pypsa.network.io:Retrieving network data from https://github.com/PyPSA/PyPSA/raw/v1.0.5/examples/networks/scigrid-de/scigrid-de.nc.

INFO:pypsa.network.io:Imported network 'SciGrid-DE' has buses, carriers, generators, lines, loads, storage_units, transformers

Preparation¶

For illustrative purposes and a simpler postprocessing workflow, we cluster the network based on federal states. For a more detailed guide on clustering, please go to Network Clustering. Within the scope of this guide, you can ignore the following blocks.

In [3]:

n.calculate_dependent_values()
n.lines = n.lines.reindex(columns=n.components["Line"]["defaults"].index[1:])
n.lines["type"] = "Al/St 240/40 2-bundle 220.0"
n.buses = n.buses.reindex(columns=n.components["Bus"]["defaults"].index[1:])
n.buses["frequency"] = 50

url = "https://media.githubusercontent.com/media/wmgeolab/geoBoundaries/9469f09592ced973a3448cf66b6100b741b64c0d/releaseData/gbOpen/DEU/ADM1/geoBoundaries-DEU-ADM1-all.zip"
states = gpd.read_file(url, layer="geoBoundaries-DEU-ADM1_simplified")
states["shapeName"] = states["shapeName"].apply(
    lambda x: x.encode("latin1").decode("utf-8")
)  # fix encoding issue

bus_coords = gpd.GeoDataFrame(
    geometry=gpd.points_from_xy(n.buses.x, n.buses.y, crs=4326), index=n.buses.index
)
busmap = bus_coords.to_crs(3035).sjoin_nearest(states.to_crs(3035), how="left").shapeISO
nc = n.cluster.cluster_by_busmap(busmap)

n.calculate_dependent_values() n.lines = n.lines.reindex(columns=n.components["Line"]["defaults"].index[1:]) n.lines["type"] = "Al/St 240/40 2-bundle 220.0" n.buses = n.buses.reindex(columns=n.components["Bus"]["defaults"].index[1:]) n.buses["frequency"] = 50 url = "https://media.githubusercontent.com/media/wmgeolab/geoBoundaries/9469f09592ced973a3448cf66b6100b741b64c0d/releaseData/gbOpen/DEU/ADM1/geoBoundaries-DEU-ADM1-all.zip" states = gpd.read_file(url, layer="geoBoundaries-DEU-ADM1_simplified") states["shapeName"] = states["shapeName"].apply( lambda x: x.encode("latin1").decode("utf-8") ) # fix encoding issue bus_coords = gpd.GeoDataFrame( geometry=gpd.points_from_xy(n.buses.x, n.buses.y, crs=4326), index=n.buses.index ) busmap = bus_coords.to_crs(3035).sjoin_nearest(states.to_crs(3035), how="left").shapeISO nc = n.cluster.cluster_by_busmap(busmap)

By default, calling n.plot() will render all network components based on the x and y coordinates defined in n.buses. This allows us to get a first visual overview on the two networks, before and after clustering.

In [4]:

n.plot()

n.plot()

Out[4]:

{'nodes': {'Bus': <matplotlib.collections.PatchCollection at 0x792cf3f72900>},
 'branches': {'Line': <matplotlib.collections.LineCollection at 0x792cf3f72a50>,
  'Transformer': <matplotlib.collections.LineCollection at 0x792d10b80f50>},
 'flows': {}}

In [5]:

nc.plot()

nc.plot()

Out[5]:

{'nodes': {'Bus': <matplotlib.collections.PatchCollection at 0x792cf3fdefd0>},
 'branches': {'Line': <matplotlib.collections.LineCollection at 0x792cf3fdf110>},
 'flows': {}}

Retrieving Results Data¶

To map result to parameters of n.plot(), we first solve the network and then use n.statistics() to calculate relevant metrics.

In [6]:

# We reduce logging output for clarity
import logging

logging.getLogger("pypsa").setLevel(logging.ERROR)
logging.getLogger("linopy").setLevel(logging.ERROR)

nc.optimize()

# We reduce logging output for clarity import logging logging.getLogger("pypsa").setLevel(logging.ERROR) logging.getLogger("linopy").setLevel(logging.ERROR) nc.optimize()

Writing constraints.:   0%|          | 0/14 [00:00<?, ?it/s]

Writing constraints.:  14%|█▍        | 2/14 [00:00<00:00, 12.55it/s]

Writing constraints.:  93%|█████████▎| 13/14 [00:00<00:00, 57.85it/s]

Writing constraints.: 100%|██████████| 14/14 [00:00<00:00, 51.38it/s]

Writing continuous variables.:   0%|          | 0/5 [00:00<?, ?it/s]

Writing continuous variables.: 100%|██████████| 5/5 [00:00<00:00, 82.96it/s]

Running HiGHS 1.11.0 (git hash: 364c83a): Copyright (c) 2025 HiGHS under MIT licence terms
LP   linopy-problem-513ptrq3 has 76872 rows; 37608 cols; 117490 nonzeros
Coefficient ranges:
  Matrix [4e-01, 9e+01]
  Cost   [3e+00, 1e+02]
  Bound  [0e+00, 0e+00]
  RHS    [2e-11, 2e+04]
Presolving model
1656 rows, 27772 cols, 32438 nonzeros  0s
1512 rows, 5899 cols, 10325 nonzeros  0s
Dependent equations search running on 1488 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
1488 rows, 5875 cols, 10459 nonzeros  0s
Presolve : Reductions: rows 1488(-75384); columns 5875(-31733); elements 10459(-107031)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0    -1.1463614410e+00 Pr: 576(1.25327e+07) 0s
       2284     4.6934521029e+06 Pr: 0(0); Du: 0(1.02802e-12) 0s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-513ptrq3
Model status        : Optimal
Simplex   iterations: 2284
Objective value     :  4.6934521029e+06
P-D objective error :  3.6709580533e-15
HiGHS run time      :          0.14

Writing the solution to /tmp/linopy-solve-dxfk0doz.sol

Out[6]:

('ok', 'optimal')

From above we learned that bus_size accepts parameters of type float, dict, and pd.Series. When passing a multi-index pd.Series, its values will be mapped to pie chart slices.

In [7]:

eb = (
    nc.statistics.energy_balance(
        groupby=["bus", "carrier"],
        components=["Generator", "Load", "StorageUnit"],
    )
    .groupby(["bus", "carrier"])
    .sum()
)

eb = ( nc.statistics.energy_balance( groupby=["bus", "carrier"], components=["Generator", "Load", "StorageUnit"], ) .groupby(["bus", "carrier"]) .sum() )

We also extract branch results, e.g., line and link flows in this example.

In [8]:

line_flow = nc.lines_t.p0.sum(axis=0)
link_flow = nc.links_t.p0.sum(axis=0)

line_flow = nc.lines_t.p0.sum(axis=0) link_flow = nc.links_t.p0.sum(axis=0)

Note that for the pie slices to be plotted and colored correctly, passing a multi-index pd.Series requires all carrier colors to exist. Colors can be specified by their hex code representation or from the list of matplotlib names. In n.statistics.energy_balance() load is also included, so we also need to include a color for the load carrier.

In [9]:

colors = {
    "Multiple": "pink",
    "AC": "black",
    "Brown Coal": "saddlebrown",
    "Gas": "darkorange",
    "Geothermal": "firebrick",
    "Hard Coal": "darkslategray",
    "Nuclear": "mediumorchid",
    "Oil": "peru",
    "Other": "dimgray",
    "Pumped Hydro": "cornflowerblue",
    "Run of River": "royalblue",
    "Solar": "gold",
    "Storage Hydro": "navy",
    "Waste": "olive",
    "Wind Offshore": "teal",
    "Wind Onshore": "turquoise",
}

nc.carriers.color = nc.carriers.index.map(colors)

colors = { "Multiple": "pink", "AC": "black", "Brown Coal": "saddlebrown", "Gas": "darkorange", "Geothermal": "firebrick", "Hard Coal": "darkslategray", "Nuclear": "mediumorchid", "Oil": "peru", "Other": "dimgray", "Pumped Hydro": "cornflowerblue", "Run of River": "royalblue", "Solar": "gold", "Storage Hydro": "navy", "Waste": "olive", "Wind Offshore": "teal", "Wind Onshore": "turquoise", } nc.carriers.color = nc.carriers.index.map(colors)

As the carriers for loads are missing, we need to add them, manually.

In [10]:

nc.carriers.loc["", "color"] = "darkred"
nc.carriers.loc["-", "color"] = "darkred"

nc.carriers.loc["", "color"] = "darkred" nc.carriers.loc["-", "color"] = "darkred"

Balances and flow¶

We first choose a suitable projection for the plot by importing cartopy and matplotlib. Commonly used projections include ccrs.Mercator() or ccrs.EqualEarth(). We can also pass geomap_color=True to get default colorings for land and water bodies. Note that this requires cartopy to be installed.

In [11]:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 7), subplot_kw={"projection": ccrs.EqualEarth()})
nc.plot(ax=ax, geomap_color=True)

import cartopy.crs as ccrs import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(8, 7), subplot_kw={"projection": ccrs.EqualEarth()}) nc.plot(ax=ax, geomap_color=True)

Out[11]:

{'nodes': {'Bus': <matplotlib.collections.PatchCollection at 0x792cec458190>},
 'branches': {'Line': <matplotlib.collections.LineCollection at 0x792cec4582d0>},
 'flows': {}}

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/v1.0.5/lib/python3.13/site-packages/cartopy/io/__init__.py:242: DownloadWarning:

Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_land.zip

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/v1.0.5/lib/python3.13/site-packages/cartopy/io/__init__.py:242: DownloadWarning:

Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_ocean.zip

Next, let's pass the results to n.explore(). Setting bus_split_circle=True maps negative values to the bottom half and positive values to the positive half. If set to False, bottom half circles are not used and negative values will automatically be omitted. Note that the area of bus_size scales proportionally to the values passed. To get a useful map, we need to scale the values, according to personal preference.

In [12]:

bus_size_factor = 3e6
branch_width_factor = 2e4
branch_flow_factor = 5e4

nc.plot(
    ax=ax,
    bus_size=eb / bus_size_factor,
    bus_split_circle=True,
    line_width=line_flow / branch_width_factor,
    link_width=link_flow / branch_width_factor,
    line_flow=line_flow / branch_flow_factor,
)

fig

bus_size_factor = 3e6 branch_width_factor = 2e4 branch_flow_factor = 5e4 nc.plot( ax=ax, bus_size=eb / bus_size_factor, bus_split_circle=True, line_width=line_flow / branch_width_factor, link_width=link_flow / branch_width_factor, line_flow=line_flow / branch_flow_factor, ) fig

Out[12]:

Legends¶

We import additional functions from pypsa.plot to add legends to our figure. It makes sense to select values that are close to values represented in the figure, e.g. the maximum, minimum values and something in between.

In [13]:

from pypsa.plot import add_legend_lines, add_legend_patches, add_legend_semicircles

print(f"Max flow: {line_flow.abs().max()}")
print(f"Min flow: {line_flow.abs().min()}")

from pypsa.plot import add_legend_lines, add_legend_patches, add_legend_semicircles print(f"Max flow: {line_flow.abs().max()}") print(f"Min flow: {line_flow.abs().min()}")

Max flow: 181341.97489722894
Min flow: 62.45143026589699

Based on the values, we choose the following. Not that the values need to be scaled with the same factors determined before.

In [14]:

add_legend_lines(
    ax,
    sizes=[branch / branch_width_factor for branch in [150000, 100000, 10000]],
    labels=["150", "100", "10"],
    legend_kw={"loc": "lower right", "frameon": False, "title": "Line flow (GWh)"},
)
fig

add_legend_lines( ax, sizes=[branch / branch_width_factor for branch in [150000, 100000, 10000]], labels=["150", "100", "10"], legend_kw={"loc": "lower right", "frameon": False, "title": "Line flow (GWh)"}, ) fig

Out[14]:

Now we apply the same process for bus sizes.

In [15]:

print(f"Max gen./load: {eb.groupby('bus').sum().abs().max()}")
print(f"Max gen./load: {eb.groupby('bus').sum().abs().min()}")

print(f"Max gen./load: {eb.groupby('bus').sum().abs().max()}") print(f"Max gen./load: {eb.groupby('bus').sum().abs().min()}")

Max gen./load: 123128.39888000001
Max gen./load: 62.45143000000007

In [16]:

add_legend_semicircles(
    ax,
    sizes=[bus / bus_size_factor for bus in [120000, -120000]],
    labels=["+120 GWh", "-120 GWh"],
    legend_kw={
        "loc": "upper left",
        "frameon": False,
        "bbox_to_anchor": (0.02, 0.98),
    },
)
fig

add_legend_semicircles( ax, sizes=[bus / bus_size_factor for bus in [120000, -120000]], labels=["+120 GWh", "-120 GWh"], legend_kw={ "loc": "upper left", "frameon": False, "bbox_to_anchor": (0.02, 0.98), }, ) fig

/tmp/ipykernel_6199/932776199.py:1: UserWarning:

When combining n.plot() with other plots on a geographical axis, ensure n.plot() is called first or the final axis extent is set initially (ax.set_extent(boundaries, crs=crs)) for consistent legend semicircle sizes.

Out[16]:

To add legend entries for each bus carrier, we use add_legend_patches.

In [17]:

add_legend_patches(
    ax,
    colors=list(nc.carriers.color),  # colors
    labels=list(nc.carriers.index),  # labels
    legend_kw={
        "loc": "lower center",
        "bbox_to_anchor": (0.5, -0.25),  # For offsetting
        "ncol": 4,
        "frameon": False,
    },
)

fig

add_legend_patches( ax, colors=list(nc.carriers.color), # colors labels=list(nc.carriers.index), # labels legend_kw={ "loc": "lower center", "bbox_to_anchor": (0.5, -0.25), # For offsetting "ncol": 4, "frameon": False, }, ) fig

Out[17]:

Export¶

We can export the figure to any desired format, i.e., .png, .jpg, .pdf etc.

In [18]:

fig.savefig("static-plot.jpg", dpi=150)

fig.savefig("static-plot.jpg", dpi=150)