split build_shapes rule into offshore, nuts3, country/europe rules

doc/preparation.rst

@@ -14,7 +14,9 @@ Instead we provide separate data bundles which can be obtained
 using the ``retrieve*`` rules (:ref:`data`).
 Having downloaded the necessary data, it can build a base PyPSA network with the following rules
 
-- :mod:`build_shapes` generates GeoJSON files with shapes of the countries, exclusive economic zones and `NUTS3 <https://en.wikipedia.org/wiki/Nomenclature_of_Territorial_Units_for_Statistics>`__ areas.
+- :mod:`build_offshore_shapes` generates GeoJSON files of offshore exclusive economic zones (EEZ).
+- :mod:`build_nuts3_shapes` generates GeoJSON files of `NUTS3 <https://en.wikipedia.org/wiki/Nomenclature_of_Territorial_Units_for_Statistics>`__ and OSM ADM1 areas enriched with GDP and population data.
+- :mod:`build_shapes` generates GeoJSON files of country boundaries and the Europe bounding shape.
 - :mod:`base_network` builds and stores the base network with all buses, HVAC lines and HVDC links, and determines `Voronoi cells <https://en.wikipedia.org/wiki/Voronoi_diagram>`__ for all substations.
 
 The network is then simplified by preparing **approximations** of the network model, for which it is computationally viable to co-optimize generation, storage and transmission capacities.
@@ -101,6 +103,16 @@ Rule ``build_bidding_zones``
 
 .. _shapes:
 
+Rule ``build_offshore_shapes``
+===============================
+
+.. automodule:: build_offshore_shapes
+
+Rule ``build_nuts3_shapes``
+============================
+
+.. automodule:: build_nuts3_shapes
+
 Rule ``build_shapes``
 =============================
 

doc/release_notes.rst

@@ -8,6 +8,8 @@ Release Notes
 
 .. Upcoming Release
 .. =================
+* refactor: Split :mod:`build_shapes` into three independent rules: :mod:`build_offshore_shapes` (EEZ only), :mod:`build_nuts3_shapes` (NUTS3/ADM1 regions with GDP/population data), and :mod:`build_shapes` (country and Europe boundary aggregation).
+
 * Fix: Re-introduce capital costs for non-bicharging discharge links in ``add_electricity.py``, e.g. fuel cells.
 
 * The lockfile update workflow now excludes packages published within the last 7 days to reduce the risk of pulling in broken or yanked releases (https://github.com/PyPSA/pypsa-eur/pull/2130).

rules/build_electricity.smk

@@ -153,16 +153,36 @@ rule build_bidding_zones:
         scripts("build_bidding_zones.py")
 
 
-rule build_shapes:
+rule build_offshore_shapes:
     input:
         eez=ancient(rules.retrieve_eez.output["gpkg"]),
+    output:
+        offshore_shapes=resources("offshore_shapes.geojson"),
+    log:
+        logs("build_offshore_shapes.log"),
+    benchmark:
+        benchmarks("build_offshore_shapes")
+    threads: 1
+    resources:
+        mem_mb=1500,
+    params:
+        countries=config_provider("countries"),
+    message:
+        "Building offshore shapes"
+    script:
+        scripts("build_offshore_shapes.py")
+
+
+rule build_nuts3_shapes:
+    input:
         nuts3_2021=rules.retrieve_eu_nuts_2021.output["shapes_level_3"],
         ba_adm1=f"data/osm_boundaries/build/{OSM_BOUNDARIES_DATASET['version']}/BA_adm1.geojson",
         md_adm1=f"data/osm_boundaries/build/{OSM_BOUNDARIES_DATASET['version']}/MD_adm1.geojson",
         ua_adm1=f"data/osm_boundaries/build/{OSM_BOUNDARIES_DATASET['version']}/UA_adm1.geojson",
         xk_adm1=f"data/osm_boundaries/build/{OSM_BOUNDARIES_DATASET['version']}/XK_adm1.geojson",
         nuts3_gdp=rules.retrieve_jrc_ardeco.output["ardeco_gdp"],
         nuts3_pop=rules.retrieve_jrc_ardeco.output["ardeco_pop"],
+        offshore_shapes=resources("offshore_shapes.geojson"),
         bidding_zones=lambda w: (
             resources("bidding_zones.geojson")
             if config_provider("clustering", "mode")(w) == "administrative"
@@ -171,10 +191,29 @@ rule build_shapes:
         other_gdp=rules.retrieve_gdp_per_capita.output["gdp"],
         other_pop=rules.retrieve_population_count.output["tif"],
     output:
-        country_shapes=resources("country_shapes.geojson"),
+        nuts3_shapes=resources("nuts3_shapes.geojson"),
+    log:
+        logs("build_nuts3_shapes.log"),
+    benchmark:
+        benchmarks("build_nuts3_shapes")
+    threads: 1
+    resources:
+        mem_mb=1500,
+    params:
+        countries=config_provider("countries"),
+    message:
+        "Building NUTS3 shapes"
+    script:
+        scripts("build_nuts3_shapes.py")
+
+
+rule build_shapes:
+    input:
+        nuts3_shapes=resources("nuts3_shapes.geojson"),
         offshore_shapes=resources("offshore_shapes.geojson"),
+    output:
+        country_shapes=resources("country_shapes.geojson"),
         europe_shape=resources("europe_shape.geojson"),
-        nuts3_shapes=resources("nuts3_shapes.geojson"),
     log:
         logs("build_shapes.log"),
     benchmark:
@@ -183,10 +222,9 @@ rule build_shapes:
     resources:
         mem_mb=1500,
     params:
-        config_provider("clustering", "mode"),
         countries=config_provider("countries"),
     message:
-        "Building geographical shapes"
+        "Building country and Europe shapes"
     script:
         scripts("build_shapes.py")
 

scripts/_helpers.py

@@ -11,6 +11,8 @@
 import time
 from collections.abc import Callable
 from functools import partial, wraps
+from itertools import takewhile
+from operator import attrgetter
 from pathlib import Path
 from tempfile import NamedTemporaryFile
 from typing import Literal
@@ -1088,3 +1090,31 @@ def load_costs(cost_file: str) -> pd.DataFrame:
     """
 
     return pd.read_csv(cost_file, index_col=0)
+
+
+def _simplify_polys(
+    polys, minarea=100 * 1e6, maxdistance=None, tolerance=None, filterremote=True
+):  # 100*1e6 = 100 km² if CRS is DISTANCE_CRS
+    from shapely.geometry import MultiPolygon
+
+    if isinstance(polys, MultiPolygon):
+        polys = sorted(polys.geoms, key=attrgetter("area"), reverse=True)
+        mainpoly = polys[0]
+        mainlength = mainpoly.area**0.5
+
+        if maxdistance is not None:
+            mainlength = maxdistance
+
+        if mainpoly.area > minarea:
+            polys = MultiPolygon(
+                [
+                    p
+                    for p in takewhile(lambda p: p.area > minarea, polys)
+                    if not filterremote or (mainpoly.distance(p) < mainlength)
+                ]
+            )
+        else:
+            polys = mainpoly
+    if tolerance is not None:
+        polys = polys.simplify(tolerance=tolerance)
+    return polys