Integrate build_bus_regions into base_network

Fixes a problem with the `build_bus_regions` writing to base.nc
without declaring base.nc as an output.

rules/build_electricity.smk

@@ -86,7 +86,9 @@ rule base_network:
         offshore_shapes=resources("offshore_shapes.geojson"),
         europe_shape=resources("europe_shape.geojson"),
     output:
-        resources("networks/base.nc"),
+        base_network=resources("networks/base.nc"),
+        regions_onshore=resources("regions_onshore.geojson"),
+        regions_offshore=resources("regions_offshore.geojson"),
     log:
         logs("base_network.log"),
     benchmark:
@@ -127,27 +129,6 @@ rule build_shapes:
         "../scripts/build_shapes.py"
 
 
-rule build_bus_regions:
-    params:
-        countries=config_provider("countries"),
-    input:
-        country_shapes=resources("country_shapes.geojson"),
-        offshore_shapes=resources("offshore_shapes.geojson"),
-        base_network=resources("networks/base.nc"),
-    output:
-        regions_onshore=resources("regions_onshore.geojson"),
-        regions_offshore=resources("regions_offshore.geojson"),
-    log:
-        logs("build_bus_regions.log"),
-    threads: 1
-    resources:
-        mem_mb=1000,
-    conda:
-        "../envs/environment.yaml"
-    script:
-        "../scripts/build_bus_regions.py"
-
-
 if config["enable"].get("build_cutout", False):
 
     rule build_cutout:

scripts/base_network.py

@@ -75,11 +75,11 @@ import shapely
 import shapely.prepared
 import shapely.wkt
 import yaml
-from _helpers import configure_logging, get_snapshots, set_scenario_config
+from _helpers import REGION_COLS, configure_logging, get_snapshots, set_scenario_config
 from packaging.version import Version, parse
 from scipy import spatial
 from scipy.sparse import csgraph
-from shapely.geometry import LineString, Point
+from shapely.geometry import LineString, Point, Polygon
 
 PD_GE_2_2 = parse(pd.__version__) >= Version("2.2")
 
@@ -779,9 +779,147 @@ def base_network(
     return n
 
 
+def voronoi_partition_pts(points, outline):
+    """
+    Compute the polygons of a voronoi partition of `points` within the polygon
+    `outline`. Taken from
+    https://github.com/FRESNA/vresutils/blob/master/vresutils/graph.py.
+
+    Attributes
+    ----------
+    points : Nx2 - ndarray[dtype=float]
+    outline : Polygon
+    Returns
+    -------
+    polygons : N - ndarray[dtype=Polygon|MultiPolygon]
+    """
+    points = np.asarray(points)
+
+    if len(points) == 1:
+        polygons = [outline]
+    else:
+        xmin, ymin = np.amin(points, axis=0)
+        xmax, ymax = np.amax(points, axis=0)
+        xspan = xmax - xmin
+        yspan = ymax - ymin
+
+        # to avoid any network positions outside all Voronoi cells, append
+        # the corners of a rectangle framing these points
+        vor = spatial.Voronoi(
+            np.vstack(
+                (
+                    points,
+                    [
+                        [xmin - 3.0 * xspan, ymin - 3.0 * yspan],
+                        [xmin - 3.0 * xspan, ymax + 3.0 * yspan],
+                        [xmax + 3.0 * xspan, ymin - 3.0 * yspan],
+                        [xmax + 3.0 * xspan, ymax + 3.0 * yspan],
+                    ],
+                )
+            )
+        )
+
+        polygons = []
+        for i in range(len(points)):
+            poly = Polygon(vor.vertices[vor.regions[vor.point_region[i]]])
+
+            if not poly.is_valid:
+                poly = poly.buffer(0)
+
+            with np.errstate(invalid="ignore"):
+                poly = poly.intersection(outline)
+
+            polygons.append(poly)
+
+    return polygons
+
+
+def build_bus_shapes(n, country_shapes, offshore_shapes, countries):
+    country_shapes = gpd.read_file(country_shapes).set_index("name")["geometry"]
+    offshore_shapes = gpd.read_file(offshore_shapes)
+    offshore_shapes = offshore_shapes.reindex(columns=REGION_COLS).set_index("name")[
+        "geometry"
+    ]
+
+    onshore_regions = []
+    offshore_regions = []
+
+    for country in countries:
+        c_b = n.buses.country == country
+
+        onshore_shape = country_shapes[country]
+        onshore_locs = (
+            n.buses.loc[c_b & n.buses.onshore_bus]
+            .sort_values(
+                by="substation_lv", ascending=False
+            )  # preference for substations
+            .drop_duplicates(subset=["x", "y"], keep="first")[["x", "y"]]
+        )
+        onshore_regions.append(
+            gpd.GeoDataFrame(
+                {
+                    "name": onshore_locs.index,
+                    "x": onshore_locs["x"],
+                    "y": onshore_locs["y"],
+                    "geometry": voronoi_partition_pts(
+                        onshore_locs.values, onshore_shape
+                    ),
+                    "country": country,
+                }
+            )
+        )
+
+        if country not in offshore_shapes.index:
+            continue
+        offshore_shape = offshore_shapes[country]
+        offshore_locs = n.buses.loc[c_b & n.buses.substation_off, ["x", "y"]]
+        offshore_regions_c = gpd.GeoDataFrame(
+            {
+                "name": offshore_locs.index,
+                "x": offshore_locs["x"],
+                "y": offshore_locs["y"],
+                "geometry": voronoi_partition_pts(offshore_locs.values, offshore_shape),
+                "country": country,
+            }
+        )
+        offshore_regions_c = offshore_regions_c.loc[offshore_regions_c.area > 1e-2]
+        offshore_regions.append(offshore_regions_c)
+
+    shapes = pd.concat(onshore_regions, ignore_index=True)
+
+    return onshore_regions, offshore_regions, shapes
+
+
+def append_bus_shapes(n, shapes, type):
+    """
+    Append shapes to the network. If shapes with the same component and type
+    already exist, they will be removed.
+
+    Parameters:
+        n (pypsa.Network): The network to which the shapes will be appended.
+        shapes (geopandas.GeoDataFrame): The shapes to be appended.
+        **kwargs: Additional keyword arguments used in `n.madd`.
+
+    Returns:
+        None
+    """
+    remove = n.shapes.query("component == 'Bus' and type == @type").index
+    n.mremove("Shape", remove)
+
+    offset = n.shapes.index.astype(int).max() + 1 if not n.shapes.empty else 0
+    shapes = shapes.rename(lambda x: int(x) + offset)
+    n.madd(
+        "Shape",
+        shapes.index,
+        geometry=shapes.geometry,
+        idx=shapes.name,
+        component="Bus",
+        type=type,
+    )
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
-
         from _helpers import mock_snakemake
 
         snakemake = mock_snakemake("base_network")
@@ -803,5 +941,22 @@ if __name__ == "__main__":
         snakemake.config,
     )
 
+    onshore_regions, offshore_regions, shapes = build_bus_shapes(
+        n,
+        snakemake.input.country_shapes,
+        snakemake.input.offshore_shapes,
+        snakemake.params.countries,
+    )
+
+    shapes.to_file(snakemake.output.regions_onshore)
+    append_bus_shapes(n, shapes, "onshore")
+
+    if offshore_regions:
+        shapes = pd.concat(offshore_regions, ignore_index=True)
+        shapes.to_file(snakemake.output.regions_offshore)
+        append_bus_shapes(n, shapes, "offshore")
+    else:
+        offshore_shapes.to_frame().to_file(snakemake.output.regions_offshore)
+
     n.meta = snakemake.config
-    n.export_to_netcdf(snakemake.output[0])
+    n.export_to_netcdf(snakemake.output.base_network)

scripts/build_bus_regions.py

@@ -1,218 +0,0 @@
-# -*- coding: utf-8 -*-
-# SPDX-FileCopyrightText: : 2017-2024 The PyPSA-Eur Authors
-#
-# SPDX-License-Identifier: MIT
-"""
-Creates Voronoi shapes for each bus representing both onshore and offshore
-regions.
-
-Relevant Settings
------------------
-
-.. code:: yaml
-
-    countries:
-
-.. seealso::
-    Documentation of the configuration file ``config/config.yaml`` at
-    :ref:`toplevel_cf`
-
-Inputs
-------
-
-- ``resources/country_shapes.geojson``: confer :ref:`shapes`
-- ``resources/offshore_shapes.geojson``: confer :ref:`shapes`
-- ``networks/base.nc``: confer :ref:`base`
-
-Outputs
--------
-
-- ``resources/regions_onshore.geojson``:
-
-    .. image:: img/regions_onshore.png
-        :scale: 33 %
-
-- ``resources/regions_offshore.geojson``:
-
-    .. image:: img/regions_offshore.png
-        :scale: 33 %
-
-Description
------------
-"""
-
-import logging
-
-import geopandas as gpd
-import numpy as np
-import pandas as pd
-import pypsa
-from _helpers import REGION_COLS, configure_logging, set_scenario_config
-from scipy.spatial import Voronoi
-from shapely.geometry import Polygon
-
-logger = logging.getLogger(__name__)
-
-
-def voronoi_partition_pts(points, outline):
-    """
-    Compute the polygons of a voronoi partition of `points` within the polygon
-    `outline`. Taken from
-    https://github.com/FRESNA/vresutils/blob/master/vresutils/graph.py.
-
-    Attributes
-    ----------
-    points : Nx2 - ndarray[dtype=float]
-    outline : Polygon
-    Returns
-    -------
-    polygons : N - ndarray[dtype=Polygon|MultiPolygon]
-    """
-    points = np.asarray(points)
-
-    if len(points) == 1:
-        polygons = [outline]
-    else:
-        xmin, ymin = np.amin(points, axis=0)
-        xmax, ymax = np.amax(points, axis=0)
-        xspan = xmax - xmin
-        yspan = ymax - ymin
-
-        # to avoid any network positions outside all Voronoi cells, append
-        # the corners of a rectangle framing these points
-        vor = Voronoi(
-            np.vstack(
-                (
-                    points,
-                    [
-                        [xmin - 3.0 * xspan, ymin - 3.0 * yspan],
-                        [xmin - 3.0 * xspan, ymax + 3.0 * yspan],
-                        [xmax + 3.0 * xspan, ymin - 3.0 * yspan],
-                        [xmax + 3.0 * xspan, ymax + 3.0 * yspan],
-                    ],
-                )
-            )
-        )
-
-        polygons = []
-        for i in range(len(points)):
-            poly = Polygon(vor.vertices[vor.regions[vor.point_region[i]]])
-
-            if not poly.is_valid:
-                poly = poly.buffer(0)
-
-            with np.errstate(invalid="ignore"):
-                poly = poly.intersection(outline)
-
-            polygons.append(poly)
-
-    return polygons
-
-
-def append_bus_shapes(n, shapes, type):
-    """
-    Append shapes to the network. If shapes with the same component and type
-    already exist, they will be removed.
-
-    Parameters:
-        n (pypsa.Network): The network to which the shapes will be appended.
-        shapes (geopandas.GeoDataFrame): The shapes to be appended.
-        **kwargs: Additional keyword arguments used in `n.madd`.
-
-    Returns:
-        None
-    """
-    remove = n.shapes.query("component == 'Bus' and type == @type").index
-    n.mremove("Shape", remove)
-
-    offset = n.shapes.index.astype(int).max() + 1 if not n.shapes.empty else 0
-    shapes = shapes.rename(lambda x: int(x) + offset)
-    n.madd(
-        "Shape",
-        shapes.index,
-        geometry=shapes.geometry,
-        idx=shapes.name,
-        component="Bus",
-        type=type,
-    )
-
-
-if __name__ == "__main__":
-    if "snakemake" not in globals():
-        from _helpers import mock_snakemake
-
-        snakemake = mock_snakemake("build_bus_regions")
-    configure_logging(snakemake)
-    set_scenario_config(snakemake)
-
-    countries = snakemake.params.countries
-
-    base_network = snakemake.input.base_network
-    n = pypsa.Network(base_network)
-
-    country_shapes = gpd.read_file(snakemake.input.country_shapes).set_index("name")[
-        "geometry"
-    ]
-    offshore_shapes = gpd.read_file(snakemake.input.offshore_shapes)
-    offshore_shapes = offshore_shapes.reindex(columns=REGION_COLS).set_index("name")[
-        "geometry"
-    ]
-
-    onshore_regions = []
-    offshore_regions = []
-
-    for country in countries:
-        c_b = n.buses.country == country
-
-        onshore_shape = country_shapes[country]
-        onshore_locs = (
-            n.buses.loc[c_b & n.buses.onshore_bus]
-            .sort_values(
-                by="substation_lv", ascending=False
-            )  # preference for substations
-            .drop_duplicates(subset=["x", "y"], keep="first")[["x", "y"]]
-        )
-        onshore_regions.append(
-            gpd.GeoDataFrame(
-                {
-                    "name": onshore_locs.index,
-                    "x": onshore_locs["x"],
-                    "y": onshore_locs["y"],
-                    "geometry": voronoi_partition_pts(
-                        onshore_locs.values, onshore_shape
-                    ),
-                    "country": country,
-                }
-            )
-        )
-
-        if country not in offshore_shapes.index:
-            continue
-        offshore_shape = offshore_shapes[country]
-        offshore_locs = n.buses.loc[c_b & n.buses.substation_off, ["x", "y"]]
-        offshore_regions_c = gpd.GeoDataFrame(
-            {
-                "name": offshore_locs.index,
-                "x": offshore_locs["x"],
-                "y": offshore_locs["y"],
-                "geometry": voronoi_partition_pts(offshore_locs.values, offshore_shape),
-                "country": country,
-            }
-        )
-        offshore_regions_c = offshore_regions_c.loc[offshore_regions_c.area > 1e-2]
-        offshore_regions.append(offshore_regions_c)
-
-    shapes = pd.concat(onshore_regions, ignore_index=True)
-    shapes.to_file(snakemake.output.regions_onshore)
-    append_bus_shapes(n, shapes, "onshore")
-
-    if offshore_regions:
-        shapes = pd.concat(offshore_regions, ignore_index=True)
-        shapes.to_file(snakemake.output.regions_offshore)
-        append_bus_shapes(n, shapes, "offshore")
-
-    else:
-        offshore_shapes.to_frame().to_file(snakemake.output.regions_offshore)
-
-    # save network with shapes
-    n.export_to_netcdf(base_network)

scripts/cluster_network.py

@@ -135,7 +135,7 @@ import pypsa
 import seaborn as sns
 from _helpers import configure_logging, set_scenario_config, update_p_nom_max
 from add_electricity import load_costs
-from build_bus_regions import append_bus_shapes
+from base_network import append_bus_shapes
 from packaging.version import Version, parse
 from pypsa.clustering.spatial import (
     busmap_by_greedy_modularity,

scripts/simplify_network.py

@@ -95,7 +95,7 @@ import pypsa
 import scipy as sp
 from _helpers import configure_logging, set_scenario_config, update_p_nom_max
 from add_electricity import load_costs
-from build_bus_regions import append_bus_shapes
+from base_network import append_bus_shapes
 from cluster_network import cluster_regions, clustering_for_n_clusters
 from pypsa.clustering.spatial import (
     aggregateoneport,