From 2579c251ca2a2ae13f9ee1bbb940c3df5ca569d6 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 13:11:14 +0200 Subject: [PATCH 01/54] Restore country field --- README.md | 2 ++ workflow/scripts/build_network.py | 39 +++++++++++++++++++++++++++---- workflow/scripts/clean.py | 3 +++ 3 files changed, 40 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 821a651..172dfc2 100644 --- a/README.md +++ b/README.md @@ -16,6 +16,8 @@ A modular Snakemake workflow for retrieving OpenStreetMap power infrastructure. The workflow retains AC substations, overhead lines, and cables at configured voltage levels, then creates generic buses, connected line segments, and voltage-pair transformers. The outputs preserve OSM provenance and geometry but contain no PyPSA-specific line types, capacities, or electrical-component assumptions. +Buses and lines carry the country they belong to, along with their construction status and planned start date. Country information is essential to resolve assign a correct line types which is strongly regional-specific. + This module follows the Modelblocks conventions (https://www.modelblocks.org). For more information, consult the [integration example](./tests/integration/Snakefile) and the `snakemake` [modularisation documentation](https://snakemake.readthedocs.io/en/stable/snakefiles/modularization.html). ## Overview diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index e4803af..50efaed 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -118,7 +118,14 @@ def _remove_loops_from_multiline(multiline: Any) -> Any: def _add_line_endings(lines: gpd.GeoDataFrame) -> pd.DataFrame: """Create deterministic virtual buses at each unique (voltage, endpoint) combination.""" - line_data = lines[["voltage", "geometry", "line_id"]] + line_data = lines[ + [ + "voltage", + "geometry", + "line_id", + "country", + ] + ] line_geoms = line_data["geometry"].apply(_remove_loops_from_multiline) endpoints0 = line_data.assign( @@ -144,7 +151,13 @@ def create_bus_data(group: pd.DataFrame) -> pd.Series: candidates = endpoint_names[numeric_parts == min_numeric] bus_id = candidates.sort_values().iloc[0] osm_ids = list(set(group["osm_id"].tolist())) - return pd.Series({"bus_id": bus_id, "contains": osm_ids}) + return pd.Series( + { + "bus_id": bus_id, + "contains": osm_ids, + "country": _merge_country_codes(group["country"]), + } + ) endpoints = ( endpoints.groupby(["voltage", "geometry"]) @@ -157,7 +170,15 @@ def create_bus_data(group: pd.DataFrame) -> pd.Series: + "-" + (endpoints["voltage"] / 1000).astype(int).astype(str) ) - return endpoints[["bus_id", "voltage", "geometry", "contains"]] + return endpoints[ + [ + "bus_id", + "voltage", + "geometry", + "contains", + "country", + ] + ] def _split_linestring_by_point( @@ -304,7 +325,15 @@ def _create_merge_mapping( unique_lines = pd.Series(itertools.chain(*buses_to_remove["line_id"])).unique() lines_to_merge = lines.loc[ lines["line_id"].isin(unique_lines), - ["line_id", "voltage", "circuits", "length", "geometry", "underground"], + [ + "line_id", + "voltage", + "circuits", + "length", + "geometry", + "underground", + "country", + ], ] lines_to_merge_dict = [ (node, row.to_dict()) @@ -347,6 +376,7 @@ def _create_merge_mapping( "voltage": voltage, "geometry": geometry, "underground": underground, + "country": country, "contains_lines": contains_lines, "contains_buses": contains_buses, } @@ -358,6 +388,7 @@ def _create_merge_mapping( "voltage", "geometry", "underground", + "country", "contains_lines", "contains_buses", ] diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 32f5fcd..dfb5e67 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -515,6 +515,7 @@ def _aggregate_lines(df_lines: pd.DataFrame) -> pd.DataFrame: "line_id", "circuits", "voltage", + "country", "underground", "under_construction", "start_date", @@ -534,6 +535,7 @@ def _finalise_lines(df_lines: pd.DataFrame) -> pd.DataFrame: "line_id", "circuits", "voltage", + "country", "underground", "under_construction", "start_date", @@ -1113,6 +1115,7 @@ def clean( "line_id", "circuits", "voltage", + "country", "underground", "under_construction", "start_date", From cdfe34eea8b1b0f016c7a00ab95d99dce10eb8b4 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 13:25:25 +0200 Subject: [PATCH 02/54] Polish LLMish wording --- workflow/internal/tag_corrections.yaml | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/workflow/internal/tag_corrections.yaml b/workflow/internal/tag_corrections.yaml index ec83725..395ced8 100644 --- a/workflow/internal/tag_corrections.yaml +++ b/workflow/internal/tag_corrections.yaml @@ -1,12 +1,9 @@ # Module data that users cannot modify: ordered, literal string corrections # for known-bad OSM tag values, applied by clean.py's `_clean_*` -# functions before their general syntax cleanup. Order matters — a later -# entry can depend on an earlier one having already run — and each list is +# functions before their general syntax cleanup. Order matters as a later +# entry can depend on an earlier one having already run. Each list is # applied verbatim, so add new entries at the point in the sequence where -# they belong, not just at the end. `lower: true` marks where the original -# PyPSA-Eur code lowercases the column; for most columns that happens before -# any corrections, but `circuits` lowercases only after its first two -# (case-sensitive) corrections, so entries must stay in that exact order. +# they belong, not just at the end. # # Ported verbatim from PyPSA-Eur's clean_osm_data.py (fix-osm branch), # including one inert entry: voltage's "400/220/110 kV'" pattern contains an From 07e58b7720a59c06de6597905d7df779b6958f22 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 13:28:28 +0200 Subject: [PATCH 03/54] Removing excessive details which stem from agentic testing --- workflow/internal/tag_corrections.yaml | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/workflow/internal/tag_corrections.yaml b/workflow/internal/tag_corrections.yaml index 395ced8..9a6e12e 100644 --- a/workflow/internal/tag_corrections.yaml +++ b/workflow/internal/tag_corrections.yaml @@ -11,16 +11,6 @@ # it can never match. Preserved as-is rather than "fixed" to keep behaviour # identical to the reference. # -# One deliberate deviation: PyPSA-Eur applies "low"/"minor"/"medium"/"med"/ -# "m"/"high" as substring replacements (like everything else in this list), -# but they're meant to catch a tag whose *entire* value is that word (some -# OSM contributors classify voltage qualitatively instead of numerically). -# As a substring match, "m" in particular corrupts any unrelated freeform -# text containing the letter m - seen for real on a Philippines line tagged -# voltage="...Max Amperes 100 Amps...", where every stray "m" expanded into -# "33000" and produced a 30+ digit value that overflowed int64 downstream. -# `exact` (whole-value match) below fixes that without changing behaviour -# for any tag that's genuinely just "medium"/"m"/etc. voltage: - lower: true From 918107d87d3003d001cab2c04fe46e273a49609b Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 13:30:53 +0200 Subject: [PATCH 04/54] Bringing in global mappings --- workflow/internal/tag_corrections.yaml | 50 ++++++++++++++++++++++++-- 1 file changed, 48 insertions(+), 2 deletions(-) diff --git a/workflow/internal/tag_corrections.yaml b/workflow/internal/tag_corrections.yaml index 9a6e12e..35ef0cf 100644 --- a/workflow/internal/tag_corrections.yaml +++ b/workflow/internal/tag_corrections.yaml @@ -5,8 +5,17 @@ # applied verbatim, so add new entries at the point in the sequence where # they belong, not just at the end. # -# Ported verbatim from PyPSA-Eur's clean_osm_data.py (fix-osm branch), -# including one inert entry: voltage's "400/220/110 kV'" pattern contains an +# The code merges the mappings from clean_osm_data: +# - global ones from PyPSA-Earth's (repl_circuits, repl_cables, repl_voltage) +# - regional from PyPSA-Eur's regional values (fix-osm branch) +# Those cover values that encode arithmetic ("2x3", "3+3") or a range ("2-1") +# which would otherwise turn into a concatenation, e.g. transformimg "2x3" +# into 23 cables instead of 6. + +# Entries are `exact` (whole-value) rather than `replace` (substring) so they +# cannot corrupt a longer freeform tag that happens to contain them. +# +# The port includes one inert entry: voltage's "400/220/110 kV'" pattern contains an # uppercase V, but by the time it runs the column is already lowercased, so # it can never match. Preserved as-is rather than "fixed" to keep behaviour # identical to the reference. @@ -27,6 +36,11 @@ voltage: - exact: ["med", "33000"] - exact: ["m", "33000"] - exact: ["high", "150000"] + # earth: repl_voltage. Must precede the "kv" -> "000" step below, which + # would otherwise read "19.1 kv" as 191000 instead of 19100. + - exact: ["19.1 kv", "19100"] + - exact: ["2*220000", "220000;220000"] + - exact: ["kv30", "30000"] - replace: ["23000-109000", "109000"] - replace: ["380000>220000", "380000;220000"] - replace: [":", ";"] @@ -42,11 +56,43 @@ circuits: - lower: true - replace: ["1,5", "3"] - replace: ["1/3", "1"] + # earth: repl_circuits. "2/3" and "2-1" would strip to 23 and 21. + - exact: ["2/3", "2"] + - exact: ["2-1", "2"] + - exact: ["1;1 disused", "1;0"] + - exact: ["single", "1"] + - exact: ["1.", "1"] + # A typo still means at least one circuit, per PyPSA-Earth's comment. + - exact: ["`", "1"] + - exact: ["^1", "1"] + - exact: ["e", "1"] + - exact: ["d", "1"] cables: - lower: true - replace: ["1/3", "1"] - replace: ["3x2;2", "3"] + # earth: repl_cables. Every arithmetic form below would otherwise strip to + # a concatenation: "3+3" to 33 cables, "6+1" to 61, "2x3" to 23. + - exact: ["3+3", "6"] + - exact: ["6+1", "6"] + - exact: ["2x3", "6"] + - exact: ["3x2", "6"] + - exact: ["2x2", "4"] + - exact: ["2-1", "3"] + - exact: ["1 (looped - haul & return) + 1 power wire", "1"] + # Zero live conductors: the way is a ground wire or a retired cable. + - exact: ["1 disused", "0"] + - exact: ["ground", "0"] + - exact: ["3;3 disused", "3;0"] + - exact: ["single", "1"] + - exact: ["triple", "3"] + - exact: ["line", "1"] + - exact: ["partial", "1"] + - exact: ["`", "1"] + - exact: ["^1", "1"] + - exact: ["e", "1"] + - exact: ["d", "1"] wires: - lower: true From e4816faad87f41764f7690fc6605dc8b0eae58bd Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 13:32:05 +0200 Subject: [PATCH 05/54] Add human TODOs to OSM corrections --- workflow/internal/tag_corrections.yaml | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/workflow/internal/tag_corrections.yaml b/workflow/internal/tag_corrections.yaml index 35ef0cf..fa7c729 100644 --- a/workflow/internal/tag_corrections.yaml +++ b/workflow/internal/tag_corrections.yaml @@ -20,6 +20,14 @@ # it can never match. Preserved as-is rather than "fixed" to keep behaviour # identical to the reference. # +# NB PyPSA-Earth's float artefacts ("50.0" -> "50") are deliberately not ported, +# as they should be covered by earth_osm's GeoJSON round-trip usaga of raw OSM +# strings (TODO may need some additional checking). +# +# NB2 PyPSA-Eur applies "low"/"minor"/"medium"/"med"/"m"/"high" as substring +# replacements,but they're meant to catch a tag whose *entire* value is that +# word as some OSM contributors classify voltage qualitatively instead of +# numerically (TODO Needs verification with MapYourGrid colleagues). voltage: - lower: true From fd923dac2bd507bc01cd693dc72901fdc0f0debd Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 13:37:34 +0200 Subject: [PATCH 06/54] Add OSM features to config validator --- workflow/rules/retrieve.smk | 14 ++++++-------- workflow/scripts/_schema.py | 10 ++++++++++ 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/workflow/rules/retrieve.smk b/workflow/rules/retrieve.smk index cc6144f..4211abf 100644 --- a/workflow/rules/retrieve.smk +++ b/workflow/rules/retrieve.smk @@ -4,6 +4,8 @@ from pathlib import Path +from scripts._schema import OSM_FEATURES + # Both retrieve_osm_pbf and retrieve_osm_overpass produce this same fixed # set of six files per country (routes_relation included even when # network.include_relations is off, just empty) — see either script's @@ -13,14 +15,10 @@ from pathlib import Path # retrieve.source picks exactly one implementation for the same output # paths — defining both unconditionally would make Snakemake's DAG # ambiguous about which one produces a given {country}_{feature}.json. -_OSM_FEATURES = [ - "lines_way", - "cables_way", - "substations_way", - "substations_node", - "substations_relation", - "routes_relation", -] +# The feature list itself lives in scripts/_schema.py, so the custom-data +# validator there checks names against the same list these paths are built +# from; the Snakefile imports it before including this file. +_OSM_FEATURES = list(OSM_FEATURES) _OSM_OUTPUTS = { feature: f"/retrieve/{{country}}_{feature}.json" for feature in _OSM_FEATURES diff --git a/workflow/scripts/_schema.py b/workflow/scripts/_schema.py index 4877b15..72697b6 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -18,6 +18,16 @@ _VALID_REGIONS: frozenset[str] = frozenset(get_all_valid_codes()) +#: The fixed set of raw OSM features to enable compartibility with custom data. +OSM_FEATURES: tuple[str, ...] = ( + "lines_way", + "cables_way", + "substations_way", + "substations_node", + "substations_relation", + "routes_relation", +) + def _validate_countries(countries: list[str]) -> None: """Raise a clear error for any country not recognised by earth-osm's region list. From 37bf53f1eb00c0e91fa27d436fc8c36231758c44 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 13:49:10 +0200 Subject: [PATCH 07/54] Add treatment of country codes for cross-border elements --- workflow/scripts/build_network.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 50efaed..4e1c7aa 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -44,6 +44,18 @@ def _empty_geodataframe(columns: list[str], crs: str) -> gpd.GeoDataFrame: ) +def _merge_country_codes(values: Any) -> str: + """Clean-up country codes for multy-country entries which is essential + for cross-border elements. + """ + co_codes: set[str] = set() + for value in values: + if value is None or pd.isna(value): + continue + co_codes.update(str(value).split(";")) + return ";".join(sorted(code for code in co_codes if code)) + + def _treat_under_construction( df: pd.DataFrame, remove_under_construction: bool, remove_after: str | None ) -> pd.DataFrame: @@ -369,6 +381,9 @@ def _create_merge_mapping( if not isinstance(geometry, LineString) or geometry.is_closed: continue + country = _merge_country_codes( + graph.nodes[node].get("country") for node in subgraph.nodes() + ) subgraph_data.append( { "line_id": f"merged_{node_longest}+{len(contains_lines) - 1}", From d683649e617936cd19460992f69f777152c8ba83 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 13:49:31 +0200 Subject: [PATCH 08/54] Clarify docstring --- workflow/scripts/build_network.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 4e1c7aa..0c785ca 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -129,7 +129,11 @@ def _remove_loops_from_multiline(multiline: Any) -> Any: def _add_line_endings(lines: gpd.GeoDataFrame) -> pd.DataFrame: - """Create deterministic virtual buses at each unique (voltage, endpoint) combination.""" + """Create deterministic virtual buses at each unique (voltage, endpoint) + combination. + + A virtual bus inherits its attributes from the lines that meet there. + """ line_data = lines[ [ "voltage", From b997f1bfa81d7993b6ce2f22212662b33780dedb Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 13:52:51 +0200 Subject: [PATCH 09/54] Restore under_construction values --- workflow/scripts/build_network.py | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 0c785ca..1b861ba 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -140,6 +140,8 @@ def _add_line_endings(lines: gpd.GeoDataFrame) -> pd.DataFrame: "geometry", "line_id", "country", + "under_construction", + "start_date", ] ] line_geoms = line_data["geometry"].apply(_remove_loops_from_multiline) @@ -172,6 +174,8 @@ def create_bus_data(group: pd.DataFrame) -> pd.Series: "bus_id": bus_id, "contains": osm_ids, "country": _merge_country_codes(group["country"]), + "under_construction": bool(group["under_construction"].all()), + "start_date": group["start_date"].min(), } ) @@ -193,6 +197,8 @@ def create_bus_data(group: pd.DataFrame) -> pd.Series: "geometry", "contains", "country", + "under_construction", + "start_date", ] ] @@ -349,6 +355,8 @@ def _create_merge_mapping( "geometry", "underground", "country", + "under_construction", + "start_date", ], ] lines_to_merge_dict = [ @@ -388,6 +396,14 @@ def _create_merge_mapping( country = _merge_country_codes( graph.nodes[node].get("country") for node in subgraph.nodes() ) + under_construction = any( + bool(graph.nodes[node].get("under_construction")) + for node in subgraph.nodes() + ) + start_date = pd.Series( + [graph.nodes[node].get("start_date") for node in subgraph.nodes()] + ).max() + subgraph_data.append( { "line_id": f"merged_{node_longest}+{len(contains_lines) - 1}", @@ -396,6 +412,8 @@ def _create_merge_mapping( "geometry": geometry, "underground": underground, "country": country, + "under_construction": under_construction, + "start_date": start_date, "contains_lines": contains_lines, "contains_buses": contains_buses, } @@ -408,6 +426,8 @@ def _create_merge_mapping( "geometry", "underground", "country", + "under_construction", + "start_date", "contains_lines", "contains_buses", ] @@ -438,7 +458,6 @@ def _merge_lines_over_virtual_buses( buses_merged = buses_merged[~buses_merged["bus_id"].isin(buses_to_remove)] lines_to_add = merged_lines_map.copy().reset_index(drop=True) - lines_to_add["under_construction"] = False lines_to_add["length"] = lines_to_add["geometry"].to_crs(distance_crs).length lines_to_add["contains"] = lines_to_add["contains_lines"] lines_to_add = lines_to_add[lines_merged.columns] @@ -781,10 +800,9 @@ def build_network( ``station_id``) and ``buses_polygon`` (the substation polygons scoped to the buses that made it into the output, keyed by ``bus_id``). """ - buses = substations.drop(columns=["country"]) buses = _treat_under_construction( - buses, remove_under_construction, remove_after - ).drop(columns=["start_date"]) + substations.copy(), remove_under_construction, remove_after + ) buses_polygon = substations_polygon[ substations_polygon["bus_id"].isin(buses["bus_id"]) @@ -811,9 +829,7 @@ def build_network( buses_polygon, ) - lines = _treat_under_construction( - lines, remove_under_construction, remove_after - ).drop(columns=["start_date"]) + lines = _treat_under_construction(lines, remove_under_construction, remove_after) lines = _merge_identical_lines(lines) buses["voltage"] = (np.floor(buses["voltage"] / 1000) * 1000).astype( From 69f9d488c3c3d36af9fe0181bc4384fe205322b3 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 14:03:15 +0200 Subject: [PATCH 10/54] Maintain output structure for empty entries --- workflow/scripts/build_network.py | 97 ++++++++++++++++++------------- 1 file changed, 55 insertions(+), 42 deletions(-) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 1b861ba..88a9a01 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -36,6 +36,46 @@ COORD_PRECISION = 8 +# Output schemas shared by the populated and the empty-input paths so both +# always produce the same columns. +BUS_COLUMNS = [ + "bus_id", + "station_id", + "voltage_kv", + "country", + "under_construction", + "start_date", + "osm_ids", + "geometry", +] + +LINE_COLUMNS = [ + "line_id", + "bus0", + "bus1", + "voltage_kv", + "circuits", + "length_m", + "underground", + "country", + "under_construction", + "start_date", + "osm_ids", + "geometry", +] + +TRANSFORMER_COLUMNS = [ + "transformer_id", + "station_id", + "bus0", + "bus1", + "voltage_bus0_kv", + "voltage_bus1_kv", + "geometry", +] + +STATION_POLYGON_COLUMNS = ["station_id", "geometry"] + def _empty_geodataframe(columns: list[str], crs: str) -> gpd.GeoDataFrame: """Build an empty GeoDataFrame with the given non-geometry ``columns``.""" @@ -44,6 +84,11 @@ def _empty_geodataframe(columns: list[str], crs: str) -> gpd.GeoDataFrame: ) +def _non_geometry(columns: list[str]) -> list[str]: + """Drop the geometry column, which ``_empty_geodataframe`` supplies itself.""" + return [column for column in columns if column != "geometry"] + + def _merge_country_codes(values: Any) -> str: """Clean-up country codes for multy-country entries which is essential for cross-border elements. @@ -813,13 +858,15 @@ def build_network( buses_polygon = buses_polygon.drop(columns=["voltage"]) if lines.empty: - empty_buses = _empty_geodataframe(["bus_id", "geometry"], crs=geo_crs) - empty_lines = _empty_geodataframe(["line_id", "geometry"], crs=geo_crs) + # Same columns as the populated path, so a consumer reading an + # empty result doesn't hit a different schema. + empty_buses = _empty_geodataframe(_non_geometry(BUS_COLUMNS), crs=geo_crs) + empty_lines = _empty_geodataframe(_non_geometry(LINE_COLUMNS), crs=geo_crs) empty_transformers = _empty_geodataframe( - ["transformer_id", "geometry"], crs=geo_crs + _non_geometry(TRANSFORMER_COLUMNS), crs=geo_crs ) empty_stations_polygon = _empty_geodataframe( - ["station_id", "geometry"], crs=geo_crs + _non_geometry(STATION_POLYGON_COLUMNS), crs=geo_crs ) return ( empty_buses, @@ -929,21 +976,7 @@ def _contains_to_osm_ids(value: Any) -> str: lines_out["osm_ids"] = lines_out["contains_lines"].apply(_contains_to_osm_ids) lines_out["length_m"] = lines_out["length"].round(2) lines_out = gpd.GeoDataFrame( - lines_out[ - [ - "line_id", - "bus0", - "bus1", - "voltage_kv", - "circuits", - "length_m", - "underground", - "osm_ids", - "geometry", - ] - ], - geometry="geometry", - crs=geo_crs, + lines_out[LINE_COLUMNS], geometry="geometry", crs=geo_crs ) transformers_out = transformers.copy() @@ -955,34 +988,14 @@ def _contains_to_osm_ids(value: Any) -> str: transformers_out["voltage_bus1"] / 1000 ).astype(int) transformers_out = gpd.GeoDataFrame( - transformers_out[ - [ - "transformer_id", - "station_id", - "bus0", - "bus1", - "voltage_bus0_kv", - "voltage_bus1_kv", - "geometry", - ] - ], - geometry="geometry", - crs=geo_crs, + transformers_out[TRANSFORMER_COLUMNS], geometry="geometry", crs=geo_crs ) else: transformers_out = _empty_geodataframe( - [ - "transformer_id", - "station_id", - "bus0", - "bus1", - "voltage_bus0_kv", - "voltage_bus1_kv", - ], - crs=geo_crs, + _non_geometry(TRANSFORMER_COLUMNS), crs=geo_crs ) - stations_polygon_out = stations[["station_id", "geometry"]].copy() + stations_polygon_out = stations[STATION_POLYGON_COLUMNS].copy() buses_polygon_out = buses_polygon.copy() return ( From 9cf45ee93009278100ceff192388c75aedbfc041 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 15:02:48 +0200 Subject: [PATCH 11/54] Correct citation --- CITATION.cff | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/CITATION.cff b/CITATION.cff index be879f0..36a6672 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -10,3 +10,16 @@ authors: - given-names: "Bobby" family-names: "Xiong" orcid: "https://orcid.org/0000-0003-2854-0730" + - given-names: "Daniele" + family-names: "Lerede" + orcid: "" + - given-names: "Davide" + family-names: "Fioriti" + orcid: "" + - given-names: "Ekaterina" + family-names: "Fedotova" + orcid: "" + - given-names: "Emmanuel" + family-names: "Bolarinwa" + orcid: "" + From 932fad7bf4d900979a7a7bedba58811ddf4f4dc7 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 15:08:24 +0200 Subject: [PATCH 12/54] Sketch custom data functionality --- config/config.schema.json | 13 ++++++++++++ config/config.yaml | 3 +++ workflow/rules/network.smk | 32 ++++++++++++++++++++++++----- workflow/scripts/_schema.py | 41 +++++++++++++++++++++++++++++++++++++ 4 files changed, 84 insertions(+), 5 deletions(-) diff --git a/config/config.schema.json b/config/config.schema.json index 6b20258..9a55950 100644 --- a/config/config.schema.json +++ b/config/config.schema.json @@ -209,6 +209,19 @@ "description": "Country-specific network overrides loaded from config/regions", "type": "object" }, + "custom_data": { + "additionalProperties": false, + "description": "Non-OSM elements to clean alongside the retrieved ones.\n\nOSM's high-voltage coverage is uneven, and today the only way to correct\na missing or mistagged asset is to edit OSM upstream and wait for the\nnext extract, which also makes a study hard to reproduce. Files listed\nhere are read by ``clean`` exactly like retrieved ones, because both\nretrieval backends already write the same raw shape.", + "properties": { + "files": { + "description": "Paths to extra raw files in the same shape retrieval writes, named ``{country}_{feature}.json`` so ``clean`` picks up the country and feature from the filename. Their elements are added to the retrieved ones", + "items": { + "type": "string" + }, + "type": "array" + } + } + }, "interactive_map": { "additionalProperties": false, "description": "Geometry simplification and coordinate rounding for build_interactive_map.py.", diff --git a/config/config.yaml b/config/config.yaml index d508946..0d043d5 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -29,6 +29,9 @@ network: regions: {} +custom_data: + files: [] + interactive_map: coordinate_decimals: 5 simplify_geometries: diff --git a/workflow/rules/network.smk b/workflow/rules/network.smk index 3d5476b..957512a 100644 --- a/workflow/rules/network.smk +++ b/workflow/rules/network.smk @@ -2,34 +2,56 @@ # # SPDX-License-Identifier: MIT +from pathlib import Path + + +def _custom_files(feature): + """Custom raw files for one feature, cleaned alongside the retrieved ones. + + Selected by filename, the same ``{country}_{feature}.json`` convention + retrieval writes and ``clean`` reads the country back out of, so a + custom file needs no special handling downstream. + """ + return [ + path + for path in config["custom_data"]["files"] + if Path(path).stem.endswith(f"_{feature}") + ] + rule clean: input: lines_way=expand( "/retrieve/{country}_lines_way.json", country=config["countries"], - ), + ) + + _custom_files("lines_way"), cables_way=expand( "/retrieve/{country}_cables_way.json", country=config["countries"], - ), + ) + + _custom_files("cables_way"), substations_way=expand( "/retrieve/{country}_substations_way.json", country=config["countries"], - ), + ) + + _custom_files("substations_way"), substations_node=expand( "/retrieve/{country}_substations_node.json", country=config["countries"], - ), + ) + + _custom_files("substations_node"), substations_relation=expand( "/retrieve/{country}_substations_relation.json", country=config["countries"], - ), + ) + + _custom_files("substations_relation"), routes_relation=( expand( "/retrieve/{country}_routes_relation.json", country=config["countries"], ) + + _custom_files("routes_relation") if config["network"]["include_relations"] else [] ), diff --git a/workflow/scripts/_schema.py b/workflow/scripts/_schema.py index 72697b6..16c2c95 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -179,6 +179,43 @@ class NetworkConfig(ConfigModel): ) +class CustomDataConfig(ConfigModel): + """Non-OSM elements to clean alongside the retrieved ones. + + OSM's high-voltage coverage is uneven, and today the only way to correct + a missing or mistagged asset is to edit OSM upstream and wait for the + next extract, which also makes a study hard to reproduce. Files listed + here are read by ``clean`` exactly like retrieved ones, because both + retrieval backends already write the same raw shape. + """ + + model_config = ConfigDict(extra="forbid") + + files: list[str] = Field( + default_factory=list, + description=( + "Paths to extra raw files in the same shape retrieval writes, " + "named ``{country}_{feature}.json`` so ``clean`` picks up the " + "country and feature from the filename. Their elements are added " + "to the retrieved ones" + ), + ) + + @field_validator("files") + @classmethod + def validate_file_names(cls, v: list[str]) -> list[str]: + """Reject names ``clean`` could not map back to a country and feature.""" + for path in v: + stem = Path(path).stem + if not any(stem.endswith(f"_{feature}") for feature in OSM_FEATURES): + raise ValueError( + f"Custom data file {path!r} must be named " + f"'{{country}}_{{feature}}.json', where feature is one of " + f"{', '.join(OSM_FEATURES)}." + ) + return v + + class SimplifyGeometriesConfig(ConfigModel): """Douglas-Peucker simplification tolerances for the interactive map.""" @@ -291,6 +328,10 @@ def validate_country_identifiers(cls, v: list[str]) -> list[str]: default_factory=dict, description="Country-specific network overrides loaded from config/regions", ) + custom_data: CustomDataConfig = Field( + default_factory=CustomDataConfig, + description="Non-OSM raw files cleaned alongside the retrieved ones", + ) interactive_map: InteractiveMapConfig = Field( default_factory=InteractiveMapConfig, description="Settings for build_interactive_map.py", From e2a2db94d4bb83f2b52429766cd1b2a73ee284ff Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 15:16:05 +0200 Subject: [PATCH 13/54] Fix LMM artifacts in README --- README.md | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 172dfc2..54832ce 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # grid-builder -A modular Snakemake workflow for retrieving OpenStreetMap power infrastructure. +A modular Snakemake workflow for building a model of transmission power grid for any country of the world. OpenStreetMap is used an a major source of original data and can be supplemented by custom inputs.

@@ -12,7 +12,7 @@ A modular Snakemake workflow for retrieving OpenStreetMap power infrastructure. ## About -`grid-builder` is a modular `snakemake` workflow that retrieves OpenStreetMap power infrastructure and builds a generic high-voltage network. It can be imported into another `snakemake` workflow. +`grid-builder` is a modular `snakemake` workflow that retrieves OpenStreetMap power infrastructure and builds a generic high-voltage network. Custom data can be injected by providing input files of unified structure. The `grid-builder` workflow can be imported into another `snakemake`based project. The workflow retains AC substations, overhead lines, and cables at configured voltage levels, then creates generic buses, connected line segments, and voltage-pair transformers. The outputs preserve OSM provenance and geometry but contain no PyPSA-specific line types, capacities, or electrical-component assumptions. @@ -27,7 +27,7 @@ Currently implemented: 1. Retrieve OSM substations, lines, cables, and (optionally) circuit relations by country, either from a cached local Geofabrik PBF extract or the live Overpass API. 2. Clean the raw retrieval output, filtering voltage, frequency, construction status, and future assets, and grouping relation member ways into one line per real-world circuit. 3. Merge nearby stations and line endpoints into generic buses, AC lines, and transformers. -4. Build a self-contained interactive map of the resulting network (`map.html`), with layer toggles, voltage/text filtering, and click-through OSM links — this is the workflow's default target. +4. Build an interactive map of the resulting network (`map.html`). ## Configuration @@ -35,8 +35,6 @@ Configuration lives in [`config/config.yaml`](./config/config.yaml), validated a ## Input / output structure -Please consult the [interface file](./INTERFACE.yaml) for more information. - Raw retrieval outputs use `/retrieve/{country}_{feature}.json`, one file per country and feature (`lines_way`, `cables_way`, `substations_way`, `substations_node`, `substations_relation`, `routes_relation`). Both retrieval @@ -54,10 +52,9 @@ integration example sets these roots to `resources/grid-builder` and `logs/grid-builder`. Downloaded PBF files (used for `retrieve.source: geofabrik`) are cached in `data/earth-osm` in this checkout. -DC assets (links, converters, switching stations) are out of scope: this workflow -builds a generic AC topology only, with no PyPSA-specific line types or capacities. +Please consult the [interface file](./INTERFACE.yaml) for more information. -## Development +## Dependency management We use [`pixi`](https://pixi.sh/) as our package manager for development. Once installed, run the following to clone this repository and install all dependencies. From f42224bacb162e44fcfca40690bd04fdbd855255 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 15:20:54 +0200 Subject: [PATCH 14/54] Fix LLM artifacts in config/README --- config/README.md | 66 ++++++++++++++++++++---------------------------- 1 file changed, 28 insertions(+), 38 deletions(-) diff --git a/config/README.md b/config/README.md index e5166c9..4f967e5 100644 --- a/config/README.md +++ b/config/README.md @@ -1,49 +1,39 @@ -Set `countries` to the ISO country codes that define the retrieval scope. The -workflow first loads the default `config/config.yaml`, then an optional -`config/regions/config..yaml` file for every selected country. A `regions` -mapping in the calling configuration overrides values from those regional files. - -`retrieve.source` picks the retrieval backend: `geofabrik` reads a cached local -PBF extract (`retrieve_osm_pbf.py`), `overpass` queries the live Overpass API -(`retrieve_osm_overpass.py`). Both produce the same output shape, so -`clean` doesn't need to know which one ran. `network.include_relations` -decides whether the network should consider `route=power`/`power=circuit` -relations, grouping their member ways into one line per real-world circuit; -retrieval respects this too, so relations aren't fetched at all when it's off. -`network` also controls the minimum retained AC voltage, station merge buffer -radius, construction filtering, and planned-asset cutoff date. - -The [BE+NL example](./examples/config.BE-NL.yaml) is a small European development -scope. Country files under `config/regions` are intentionally small defaults for -now; community-maintained local corrections belong there rather than in workflow -code. - -`interactive_map` controls the size of `map.html`: `coordinate_decimals` rounds -embedded coordinates, and `simplify_geometries` sets per-geometry-type -Douglas-Peucker tolerances (in metres) for station polygons, bus polygons, and -lines, or disables simplification entirely via `simplify_geometries.enable`. +Set `countries` to the ISO country codes that define the retrieval scope. The workflow first loads the default `config/config.yaml`, then an optional `config/regions/config..yaml` file for every selected country. A `regions` mapping in the calling configuration overrides values from those regional files. + +`retrieve.source` picks the retrieval backend: `geofabrik` reads a cached local PBF extract (`retrieve_osm_pbf.py`), `overpass` queries the live Overpass API (`retrieve_osm_overpass.py`). Both produce the same output shape, so `clean` doesn't need to know which one ran. `network.include_relations` decides whether the network should consider `route=power`/`power=circuit` relations, grouping their member ways into one line per real-world circuit; retrieval respects this too, so relations aren't fetched at all when it's off. `network` also controls the minimum retained AC voltage, station merge buffer radius, construction filtering, and planned-asset cutoff date. + +`network.station_merge_radius_m` is a buffer radius, not a merge distance: both sides of a pair are buffered by it before the buffers are dissolved, so two elements merge once they are within *twice* the value of each other. The default of 500 therefore merges substations up to a kilometre apart. + +`network.max_station_voltage_ratio` guards against the over-merging that radius can cause. Buffer-and-union clustering is transitive, so a dense corridor of buses can chain into one oversized station, and an implausibly wide voltage span is the usual symptom. Setting it splits any station whose highest and lowest bus voltage differ by more than that factor into one station per voltage. It is null by default, which disables the guard. + +Choose the value above the largest step a real substation in your grid steps down, and verify it against your own countries before enabling it. Too low a value splits legitimate multi-voltage yards, and because transformers are only created between buses of the same station, every wrongly split yard silently loses its transformer. Colombia is a worked example: its backbone pairs 500 kV with 230 kV, a ratio of 2.17, so `2.0` splits 14 perfectly normal stations and drops the country from 20 transformers to 4, while `3.0` changes nothing. + +The [BE+NL example](./examples/config.BE-NL.yaml) is a small European development scope. Country files under `config/regions` are intentionally small defaults for now; community-maintained local corrections belong there rather than in workflow code. + +### Adding custom data + +OSM's high-voltage coverage is uneven, and editing OSM upstream and waiting for the next extract is a slow way to correct a missing or mistagged asset. It also makes a study hard to reproduce, since upstream keeps changing. List extra raw files under `custom_data.files` to have `clean` read them alongside the retrieved ones: + +```yaml +custom_data: + files: + - data/custom/BE_lines_way.json +``` + +Each file holds raw elements in the same `{"elements": [...]}` shape both retrieval backends write, and must be named `{country}_{feature}.json`, because `clean` reads the country and feature straight back out of the filename. The feature has to be one the retrieval step produces: `lines_way`, `cables_way`, `substations_way`, `substations_node`, `substations_relation`, or `routes_relation`. Their elements are added to the retrieved ones, and the files are declared inputs of the `clean` rule, so a missing one fails the run instead of being skipped silently. + +`interactive_map` controls the size of `map.html`: `coordinate_decimals` rounds embedded coordinates, and `simplify_geometries` sets per-geometry-type Douglas-Peucker tolerances (in metres) for station polygons, bus polygons, and lines, or disables simplification entirely via `simplify_geometries.enable`. ### Personal settings and Overpass fair use -Keep `config/config.yaml` as pure defaults — a test enforces that it matches the -schema, and it's tracked in git, so it's not the place for anything -environment- or person-specific. For local overrides (a custom Overpass -endpoint, contact details, a smaller `countries` scope for development), create -an untracked `config/config.local.yaml` and pass it alongside the default: +Keep `config/config.yaml` as pure defaults — a test enforces that it matches the schema, and it's tracked in git, so it's not the place for anything environment- or person-specific. For local overrides (a custom Overpass endpoint, contact details, a smaller `countries` scope for development), create an untracked `config/config.local.yaml` and pass it alongside the default: ```shell snakemake --configfile config/config.local.yaml ... ``` -Snakemake deep-merges it on top of `config/config.yaml`, so you only need to -list the keys you're overriding. +Snakemake deep-merges it on top of `config/config.yaml`, so you only need to list the keys you're overriding. -If you use `retrieve.source: overpass`, set `retrieve.overpass_api.user_agent` -to your own project name, contact email, and website. The [Overpass API fair -use policy](https://wiki.openstreetmap.org/wiki/Overpass_API#Fair_use_policy) -expects automated queries to be identifiable and reachable; a generic or -missing user agent risks being rate-limited or blocked. `retrieve.overpass_api.url` -also lets you point at your own or a faster mirror instance instead of the -shared public endpoint, without touching the checked-in default. +If you use `retrieve.source: overpass`, set `retrieve.overpass_api.user_agent` to your own project name, contact email, and website. The [Overpass API fair use policy](https://wiki.openstreetmap.org/wiki/Overpass_API#Fair_use_policy) expects automated queries to be identifiable and reachable; a generic or missing user agent risks being rate-limited or blocked. `retrieve.overpass_api.url` also lets you point at your own or a faster mirror instance instead of the shared public endpoint, without touching the checked-in default. The generated [schema](./config.schema.json) describes every option. From a7edddd4084fb9da96e17ebcf920d44f5ae8c37e Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 15:21:51 +0200 Subject: [PATCH 15/54] Add TODO for CRS hardcoding in schema --- workflow/scripts/_schema.py | 1 + 1 file changed, 1 insertion(+) diff --git a/workflow/scripts/_schema.py b/workflow/scripts/_schema.py index 16c2c95..275a391 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -289,6 +289,7 @@ class CrsConfig(ConfigModel): "EPSG:4326", description="Geographic CRS used to store and exchange coordinates" ) distance: str = Field( + # TODO Mind European-centric hardcoding "EPSG:3035", description=( "Equal-area/equal-distance CRS used for buffering and length " From 544cfdad9beeabca40ec68c34098de55a6183312 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 16:17:57 +0200 Subject: [PATCH 16/54] Increase quality of README --- config/README.md | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/config/README.md b/config/README.md index 4f967e5..4acce0e 100644 --- a/config/README.md +++ b/config/README.md @@ -1,18 +1,16 @@ Set `countries` to the ISO country codes that define the retrieval scope. The workflow first loads the default `config/config.yaml`, then an optional `config/regions/config..yaml` file for every selected country. A `regions` mapping in the calling configuration overrides values from those regional files. -`retrieve.source` picks the retrieval backend: `geofabrik` reads a cached local PBF extract (`retrieve_osm_pbf.py`), `overpass` queries the live Overpass API (`retrieve_osm_overpass.py`). Both produce the same output shape, so `clean` doesn't need to know which one ran. `network.include_relations` decides whether the network should consider `route=power`/`power=circuit` relations, grouping their member ways into one line per real-world circuit; retrieval respects this too, so relations aren't fetched at all when it's off. `network` also controls the minimum retained AC voltage, station merge buffer radius, construction filtering, and planned-asset cutoff date. +`retrieve.source` picks the retrieval backend: `geofabrik` reads a cached local PBF extract (`retrieve_osm_pbf.py`), `overpass` queries the live Overpass API (`retrieve_osm_overpass.py`). Retrieven data are transferred to the cleaning phase and after that are used to build a topologically-clean network model. -`network.station_merge_radius_m` is a buffer radius, not a merge distance: both sides of a pair are buffered by it before the buffers are dissolved, so two elements merge once they are within *twice* the value of each other. The default of 500 therefore merges substations up to a kilometre apart. +A parameter `network.include_relations` defines whether the network should consider OSM relations `route=power`/`power=circuit` , grouping their member ways into one line per real-world circuit. In the network-building phase, `network` scripts accepts custom values the minimum retained AC voltage, station merge buffer radius and a construction status. -`network.max_station_voltage_ratio` guards against the over-merging that radius can cause. Buffer-and-union clustering is transitive, so a dense corridor of buses can chain into one oversized station, and an implausibly wide voltage span is the usual symptom. Setting it splits any station whose highest and lowest bus voltage differ by more than that factor into one station per voltage. It is null by default, which disables the guard. - -Choose the value above the largest step a real substation in your grid steps down, and verify it against your own countries before enabling it. Too low a value splits legitimate multi-voltage yards, and because transformers are only created between buses of the same station, every wrongly split yard silently loses its transformer. Colombia is a worked example: its backbone pairs 500 kV with 230 kV, a ratio of 2.17, so `2.0` splits 14 perfectly normal stations and drops the country from 20 transformers to 4, while `3.0` changes nothing. +`network.station_merge_radius_m` is a buffer radius with the merge distance being *twice* as high.E.g.the default value of 500 m merges substations up to one kilometre apart. The [BE+NL example](./examples/config.BE-NL.yaml) is a small European development scope. Country files under `config/regions` are intentionally small defaults for now; community-maintained local corrections belong there rather than in workflow code. ### Adding custom data -OSM's high-voltage coverage is uneven, and editing OSM upstream and waiting for the next extract is a slow way to correct a missing or mistagged asset. It also makes a study hard to reproduce, since upstream keeps changing. List extra raw files under `custom_data.files` to have `clean` read them alongside the retrieved ones: +The workflow provides an option to add custom data which can be handy to deal with inputs which are out of scope for OpenStreetMap, such as planned lines. To inject custom files into the worklow, a filed `custom_data.files` can be used: ```yaml custom_data: @@ -20,13 +18,13 @@ custom_data: - data/custom/BE_lines_way.json ``` -Each file holds raw elements in the same `{"elements": [...]}` shape both retrieval backends write, and must be named `{country}_{feature}.json`, because `clean` reads the country and feature straight back out of the filename. The feature has to be one the retrieval step produces: `lines_way`, `cables_way`, `substations_way`, `substations_node`, `substations_relation`, or `routes_relation`. Their elements are added to the retrieved ones, and the files are declared inputs of the `clean` rule, so a missing one fails the run instead of being skipped silently. +Enabling this functionality makes `clean` read the custom files alongside the retrieved ones. The expected format correspond to raw elements and must be named `{country}_{feature}.json`. The feature has to be one the retrieval step produces: `lines_way`, `cables_way`, `substations_way`, `substations_node`, `substations_relation`, or `routes_relation`. -`interactive_map` controls the size of `map.html`: `coordinate_decimals` rounds embedded coordinates, and `simplify_geometries` sets per-geometry-type Douglas-Peucker tolerances (in metres) for station polygons, bus polygons, and lines, or disables simplification entirely via `simplify_geometries.enable`. +`interactive_map` controls the parameters of `map.html` with `coordinate_decimals` for displayed precision of coordinates, and `simplify_geometries` in meteres applied to station polygons, bus polygons, and lines. ### Personal settings and Overpass fair use -Keep `config/config.yaml` as pure defaults — a test enforces that it matches the schema, and it's tracked in git, so it's not the place for anything environment- or person-specific. For local overrides (a custom Overpass endpoint, contact details, a smaller `countries` scope for development), create an untracked `config/config.local.yaml` and pass it alongside the default: +Keep a git-tracked `config/config.yaml` as pure defaults. A test enforces that this file matches the schema. For local overrides, such as a custom Overpass endpoint, contact details, a smaller `countries` scope, please create an untracked `config/config.local.yaml` and pass it alongside the default: ```shell snakemake --configfile config/config.local.yaml ... From 5e2c23fb63bd1df3a242625dd133165c835323d2 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 16:18:26 +0200 Subject: [PATCH 17/54] Fix refuse --- workflow/scripts/build_network.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 88a9a01..78b5cfb 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -964,11 +964,8 @@ def _contains_to_osm_ids(value: Any) -> str: buses_out = buses.copy() buses_out["voltage_kv"] = (buses_out["voltage"] / 1000).astype(int) buses_out["osm_ids"] = buses_out["contains"].apply(_contains_to_osm_ids) - buses_out = buses_out.rename(columns={"station_id": "station_id"}) buses_out = gpd.GeoDataFrame( - buses_out[["bus_id", "station_id", "voltage_kv", "osm_ids", "geometry"]], - geometry="geometry", - crs=geo_crs, + buses_out[BUS_COLUMNS], geometry="geometry", crs=geo_crs ) lines_out = lines.copy() From 7f55d299bd45da9160b1d5c49e3df425007f30a8 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 16:27:52 +0200 Subject: [PATCH 18/54] Quick fix of interactive frontend --- workflow/scripts/build_interactive_map.py | 34 ++++++++++++----------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/workflow/scripts/build_interactive_map.py b/workflow/scripts/build_interactive_map.py index 4fb5e97..9a72622 100644 --- a/workflow/scripts/build_interactive_map.py +++ b/workflow/scripts/build_interactive_map.py @@ -107,17 +107,20 @@ def path_layer( frame["geometry"] = ( frame.geometry.to_crs(distance_crs).simplify(simplify_m).to_crs(geo_crs) ) + # A PathLayer datum is one flat list of positions, so a MultiLineString has + # to become one row per part. Handing deck.gl a nested list instead throws + # inside its tesselator and takes the whole canvas down with it, not just + # the offending row, so a single multi-part line blanks the entire map. + if (frame.geom_type == "MultiLineString").any(): + frame = frame.explode(index_parts=False).reset_index(drop=True) data = tooltip(frame) data["path"] = data.geometry.map( - lambda line: ( - [ - [_coord(point, coord_decimals) for point in item.coords] - for item in line.geoms - ] - if line.geom_type == "MultiLineString" - else [_coord(point, coord_decimals) for point in line.coords] - ) + lambda line: [_coord(point, coord_decimals) for point in line.coords] ) + # Zero-length parts carry no picture and give deck.gl degenerate normals. + data = data[data["path"].map(lambda path: len({tuple(p) for p in path}) > 1)] + if data.empty: + return None return pdk.Layer( "PathLayer", data=data.drop(columns="geometry"), @@ -150,16 +153,15 @@ def polygon_layer( frame["geometry"] = ( frame.geometry.to_crs(distance_crs).simplify(simplify_m).to_crs(geo_crs) ) + # deck.gl reads a nested list as one polygon with holes, so the parts of a + # MultiPolygon would be punched out of the first part instead of drawn. + if (frame.geom_type == "MultiPolygon").any(): + frame = frame.explode(index_parts=False).reset_index(drop=True) data = tooltip(frame) data["polygon"] = data.geometry.map( - lambda polygon: ( - [ - [_coord(point, coord_decimals) for point in item.exterior.coords] - for item in polygon.geoms - ] - if polygon.geom_type == "MultiPolygon" - else [_coord(point, coord_decimals) for point in polygon.exterior.coords] - ) + lambda polygon: [ + _coord(point, coord_decimals) for point in polygon.exterior.coords + ] ) extrusion_kwargs: dict[str, Any] = {} if extruded: From 35a4d35811dfcfbbe41e3eb551849c2700f6ad61 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 16:34:00 +0200 Subject: [PATCH 19/54] Port frequency correction mappings --- workflow/internal/tag_corrections.yaml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/workflow/internal/tag_corrections.yaml b/workflow/internal/tag_corrections.yaml index fa7c729..3f9aed0 100644 --- a/workflow/internal/tag_corrections.yaml +++ b/workflow/internal/tag_corrections.yaml @@ -124,6 +124,9 @@ frequency: - replace: ["?", ""] - replace: ["hz", ""] - replace: [" ", ""] + # earth: repl_freq. Two traction entries run together; left alone this + # parses as no number at all instead of two 16.7 Hz circuits. + - exact: ["50;50;16.716.7", "50;50;16.7;16.7"] date: - lower: true From b1c1dab555951094c71500bf44a0f6258c73e0d4 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 16:42:47 +0200 Subject: [PATCH 20/54] Generalize treatment of frequency --- workflow/scripts/_schema.py | 20 +++++ workflow/scripts/clean.py | 165 +++++++++++++++++++++++++++++++----- 2 files changed, 165 insertions(+), 20 deletions(-) diff --git a/workflow/scripts/_schema.py b/workflow/scripts/_schema.py index 275a391..f55baaf 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -165,6 +165,26 @@ class NetworkConfig(ConfigModel): default_factory=FrequencyConfig, description="AC/DC frequency in Hz; override per country in config/regions for e.g. 60 Hz grids", ) + accepted_ac_frequencies_hz: list[PositiveFloat] = Field( + [50.0, 60.0], + description=( + "Frequency tag values treated as public-grid AC. A tagged value " + "within frequency_tolerance_hz of one of these is kept and " + "normalised to the region's AC frequency, since 50 Hz in a 60 Hz " + "country is far likelier a tagging slip than a separate grid. Any " + "other non-DC value is dropped: it marks a separate system such as " + "16.7 Hz railway traction, which must not be read as mains AC" + ), + min_length=1, + ) + frequency_tolerance_hz: float = Field( + 0.1, + description=( + "Tolerance for matching a frequency tag, so that 50.0 matches 50 " + "and 0.0 matches the DC marker" + ), + ge=0, + ) station_merge_radius_m: float = Field( 500.0, description="Buffer radius used to merge nearby substations and line endpoints, in metres", diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index dfb5e67..67e4f8b 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -41,6 +41,9 @@ "tag_corrections.yaml" ) +# Country values whose AC frequency overrides conflict, tracked so the +# warning in _region_ac_hz is logged once per border rather than once per row. +_AC_FREQUENCY_CONFLICTS: set[str] = set() def _create_linestring(row: pd.Series) -> LineString: """Build a LineString from a raw OSM geometry list of ``{lon, lat}`` points.""" @@ -118,6 +121,91 @@ def _clean_frequency(column: pd.Series) -> pd.Series: column = _apply_corrections(_to_str(column), _TAG_CORRECTIONS["frequency"]) return column.str.replace(r"[^0-9;.]", "", regex=True) +def _frequency_kind( + value: str, dc_hz: float, accepted_ac_hz: list[float], tolerance_hz: float +) -> str: + """Classify one cleaned frequency value as ``"ac"``, ``"dc"`` or ``"other"``. + + Comparison is numeric with a tolerance, so "50.0" reads as AC and "0.0" as + DC where a string comparison would miss both. An untagged value counts as + AC, as in PyPSA-Earth, which fills a missing frequency with the mains + default. A value that does not parse is ``"other"``, again as in + PyPSA-Earth, which coerces it to NaN and drops it. + """ + if value == "": + return "ac" + try: + hz = float(value) + except ValueError: + return "other" + if abs(hz - dc_hz) <= tolerance_hz: + return "dc" + if any(abs(hz - ac) <= tolerance_hz for ac in accepted_ac_hz): + return "ac" + return "other" + + +def _normalise_frequency( + df: pd.DataFrame, + dc_hz: str, + accepted_ac_hz: list[float], + tolerance_hz: float, + label: str, +) -> pd.DataFrame: + """Drop rows on a non-grid frequency and normalise the rest to AC or DC. + + Each row must already hold a single frequency value and a per-row + ``_ac_hz``. AC rows are rewritten to their region's AC frequency and DC + rows to the DC marker, so later code can compare plain strings. Rows on + any other frequency belong to a separate system, typically 16.7 Hz or + 25 Hz railway traction; rewriting them to mains AC, as a blanket + "invalid means AC" rule would, splices a railway into the public grid. + """ + if df.empty: + return df + dc_value = float(dc_hz) + kind = df.apply( + lambda row: _frequency_kind( + row["frequency"], + dc_value, + [*accepted_ac_hz, float(row["_ac_hz"])], + tolerance_hz, + ), + axis=1, + ) + other = kind == "other" + if other.any(): + logger.info( + "Dropping %d %s on a non-grid frequency: %s", + int(other.sum()), + label, + ", ".join(sorted(df.loc[other, "frequency"].unique())[:10]), + ) + df = df[~other].copy() + kind = kind[~other] + df.loc[kind == "ac", "frequency"] = df.loc[kind == "ac", "_ac_hz"] + df.loc[kind == "dc", "frequency"] = dc_hz + return df + + +def _frequency_for_split(row: pd.Series) -> str: + """Pick the frequency belonging to this row's voltage split. + + OSM lists one frequency per circuit in the same order as voltage, e.g. + voltage=380000;110000 with frequency=50;16.7 for a mains circuit sharing + towers with a railway traction circuit. Only voltage is split into rows, + so without this each split would carry the whole list. A shorter list is + padded with its last value, PyPSA-Earth's rule; a longer one is + truncated to the voltages present. + """ + values = row["frequency"].split(";") + if len(values) == 1: + return values[0] + position = ( + int(row["id"].rsplit("-", 1)[1]) - 1 if row["split_elements"] > 1 else 0 + ) + return values[min(position, len(values) - 1)] + def _clean_date(column: pd.Series) -> pd.Series: """Parse a raw ``start_date`` tag column to datetimes, coercing invalid values to NaT.""" @@ -225,6 +313,7 @@ def _filter_by_voltage( def _clean_substations( df_substations: pd.DataFrame, list_voltages: Any, dc_hz: str + accepted_ac_hz: list[float], ) -> pd.DataFrame: """Split multi-voltage substation rows and normalise each split's frequency. @@ -257,17 +346,16 @@ def _clean_substations( ) df_substations = _split_cells(df_substations, cols=["frequency"]) - bool_invalid_frequency = df_substations.apply( - lambda row: row["frequency"] not in (row["_ac_hz"], dc_hz), axis=1 + return _normalise_frequency( + df_substations, dc_hz, accepted_ac_hz, tolerance_hz, "substation elements" ) - df_substations.loc[bool_invalid_frequency, "frequency"] = df_substations.loc[ - bool_invalid_frequency, "_ac_hz" - ] - return df_substations def _clean_lines( - df_lines: pd.DataFrame, list_voltages: Any, dc_hz: str + df_lines: pd.DataFrame, + list_voltages: Any, + accepted_ac_hz: list[float], + tolerance_hz: float, ) -> pd.DataFrame: """Clean lines/cables heuristically, deriving circuits from whatever tags exist. @@ -285,6 +373,15 @@ def _clean_lines( _check_voltage, list_voltages=list_voltages ) df_lines = df_lines[bool_voltages] + if df_lines.empty: + return df_lines + + df_lines["frequency"] = df_lines.apply(_frequency_for_split, axis=1) + df_lines = _normalise_frequency( + df_lines, dc_hz, accepted_ac_hz, tolerance_hz, "line elements" + ) + if df_lines.empty: + return df_lines bool_ac = df_lines["frequency"] != dc_hz bool_dc = ~bool_ac @@ -925,6 +1022,21 @@ def _import_substations( # --------------------------------------------------------------------------- +def _split_country_codes(country: Any) -> list[str]: + """Split a country value into individual codes. + + ``_drop_duplicate_lines`` joins the codes of a cross-border element it + saw in more than one country's retrieval, so a single value can read + ``"BE;NL"``. Looking that up as one key silently misses every regional + override, which is why both resolvers below go through here first. + """ + # pd.isna covers None, float nan and pd.NA alike; a bare float check + # would let pd.NA through and stringify it into a bogus "" code. + if country is None or pd.isna(country): + return [] + return [code for code in str(country).split(";") if code] + + def _region_min_voltage( country: str, network: dict[str, Any], regions: dict[str, Any] ) -> float: @@ -947,11 +1059,32 @@ def _region_ac_hz( DC has no equivalent per-country override: OSM always tags DC as ``frequency=0`` worldwide, so unlike AC (50 Hz vs 60 Hz by continent) it isn't something a region should need to change. + + A cross-border element names more than one country. An AC line only ever + runs at one frequency, so disagreeing neighbours mean either a config + error or a tie that isn't a plain AC line; the lowest value is taken so + the result is deterministic, and the conflict is logged once. """ - frequency_override = regions.get(country, {}).get("frequency_hz") or {} - override = frequency_override.get("AC") - hz = override if override is not None else network["frequency_hz"]["AC"] - return _format_hz(hz) + values = [] + for code in _split_country_codes(country): + override = (regions.get(code, {}).get("frequency_hz") or {}).get("AC") + values.append( + override if override is not None else network["frequency_hz"]["AC"] + ) + if not values: + values = [network["frequency_hz"]["AC"]] + + unique = sorted(set(values)) + if len(unique) > 1 and country not in _AC_FREQUENCY_CONFLICTS: + _AC_FREQUENCY_CONFLICTS.add(country) + logger.warning( + "Countries %s disagree on AC frequency (%s Hz); using %s Hz for the " + "elements they share.", + _split_country_codes(country), + unique, + unique[0], + ) + return _format_hz(unique[0]) def clean( @@ -1159,15 +1292,7 @@ def clean( df_lines["_ac_hz"] = df_lines["country"].map( lambda c: _region_ac_hz(c, network, regions) ) - df_lines = _clean_lines(df_lines, list_voltages, dc_hz) - - len_before = len(df_lines) - df_lines = df_lines[df_lines["frequency"] != dc_hz] - logger.info( - "Dropped %d DC lines. Keeping %d AC lines.", - len_before - len(df_lines), - len(df_lines), - ) + df_lines = _clean_lines(df_lines, list_voltages, dc_hz, **frequency_options) if not df_lines.empty: row_min = df_lines["country"].map( From 0d4f431da3b9b8072ce002343115020984e0e233 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:28:59 +0200 Subject: [PATCH 21/54] Enhance cleanup --- workflow/scripts/clean.py | 39 ++++++++++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 5 deletions(-) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 67e4f8b..d46ff69 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -86,28 +86,56 @@ def _apply_corrections(column: pd.Series, steps: list[dict[str, Any]]) -> pd.Ser return column +def _strip_to(column: pd.Series, tag: str, allowed: str) -> pd.Series: + """Drop every character outside ``allowed``, reporting what was unmapped first. + + The strip is a catch-all: whatever the corrections above did not handle, + it reduces to bare digits. That is right for stray units and punctuation + but wrong for a value encoding arithmetic, where "2x3" becomes 23 rather + than 6. PyPSA-Earth's replacement tables grew precisely because it logs + the values it could not map, so log them here too instead of silently + turning a tag into a plausible but wrong number. + """ + # Only characters sitting *between* digits are reported: those are the ones + # whose removal splices two numbers into one ("2x3" -> 23). A stray unit or + # bracket at either end ("220000 V") strips away harmlessly and would + # otherwise drown the real signal in noise. + leftover = column[column.str.contains(f"[0-9][^{allowed}]+[0-9]", regex=True)] + if not leftover.empty: + counts = leftover.value_counts() + logger.warning( + "%d %s value(s) were not matched by any correction and will be reduced " + "to their digits, which may be wrong. Add them to " + "tag_corrections.yaml if so: %s", + int(counts.sum()), + tag, + ", ".join(f"{value!r} (x{n})" for value, n in counts.head(10).items()), + ) + return column.str.replace(f"[^{allowed}]", "", regex=True) + + def _clean_voltage(column: pd.Series) -> pd.Series: """Normalise a raw ``voltage`` tag column to semicolon-separated volts.""" column = _apply_corrections(_to_str(column), _TAG_CORRECTIONS["voltage"]) - return column.str.replace(r"[^0-9;]", "", regex=True) + return _strip_to(column, "voltage", "0-9;") def _clean_circuits(column: pd.Series) -> pd.Series: """Normalise a raw ``circuits`` tag column to semicolon-separated integers.""" column = _apply_corrections(_to_str(column), _TAG_CORRECTIONS["circuits"]) - return column.str.replace(r"[^0-9;]", "", regex=True) + return _strip_to(column, "circuits", "0-9;") def _clean_cables(column: pd.Series) -> pd.Series: """Normalise a raw ``cables`` tag column to semicolon-separated integers.""" column = _apply_corrections(_to_str(column), _TAG_CORRECTIONS["cables"]) - return column.str.replace(r"[^0-9;]", "", regex=True) + return _strip_to(column, "cables", "0-9;") def _clean_wires(column: pd.Series) -> pd.Series: """Normalise a raw ``wires`` tag column to semicolon-separated integers.""" column = _apply_corrections(_to_str(column), _TAG_CORRECTIONS["wires"]) - return column.str.replace(r"[^0-9;]", "", regex=True) + return _strip_to(column, "wires", "0-9;") def _check_voltage(voltage: str, list_voltages: Any) -> bool: @@ -119,7 +147,8 @@ def _check_voltage(voltage: str, list_voltages: Any) -> bool: def _clean_frequency(column: pd.Series) -> pd.Series: """Normalise a raw ``frequency`` tag column to semicolon-separated Hz values.""" column = _apply_corrections(_to_str(column), _TAG_CORRECTIONS["frequency"]) - return column.str.replace(r"[^0-9;.]", "", regex=True) + return _strip_to(column, "frequency", "0-9;.") + def _frequency_kind( value: str, dc_hz: float, accepted_ac_hz: list[float], tolerance_hz: float From a26a6c7ad2629edfee73dafa6aa3b714556782e2 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:32:09 +0200 Subject: [PATCH 22/54] Improve treatment of cables --- workflow/scripts/clean.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index d46ff69..52b4f13 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -272,12 +272,12 @@ def generate_new_id(row: pd.Series) -> str: def _distribute_to_circuits(row: pd.Series) -> str: - """Split a row's circuits (or cables/3) evenly across its ``split_elements``.""" + """Split a row's circuits (or cables per circuit) evenly across its ``split_elements``.""" circuits: float if row["circuits"] != "": circuits = int(row["circuits"]) else: - circuits = int(row["cables"]) / 3 + circuits = int(row["cables"]) / row["_cables_per_circuit"] single_circuit = int(max(1, np.floor_divide(circuits, row["split_elements"]))) return str(single_circuit) @@ -510,7 +510,8 @@ def _clean_lines( max( 1, np.floor_divide( - int(row["cables"].split(";")[int(row["id"].split("-")[-1]) - 1]), 3 + int(row["cables"].split(";")[int(row["id"].split("-")[-1]) - 1]), + row["_cables_per_circuit"], ), ) ), @@ -530,6 +531,15 @@ def _clean_lines( df_lines.loc[bool_leftover & bool_dc, "frequency"] = dc_hz df_lines.loc[bool_leftover, "cleaned"] = True + df_lines = df_lines.drop(columns=["_cables_per_circuit"]) + no_conductors = df_lines["circuits"] == "0" + if no_conductors.any(): + logger.info( + "Dropping %d line(s) tagged with no live conductor.", + int(no_conductors.sum()), + ) + df_lines = df_lines[~no_conductors] + return df_lines From d17ec6531eaad50139d6238a138733c167bedea3 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:34:24 +0200 Subject: [PATCH 23/54] Refactor output for clean script --- workflow/scripts/clean.py | 66 +++++++++++++++++++++++---------------- 1 file changed, 39 insertions(+), 27 deletions(-) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 52b4f13..5a18394 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -45,6 +45,42 @@ # warning in _region_ac_hz is logged once per border rather than once per row. _AC_FREQUENCY_CONFLICTS: set[str] = set() +# Output schemas, shared by the populated and the empty-input paths so a +# country with no substations or no lines still yields the same columns. +SUBSTATION_COLUMNS = [ + "bus_id", + "voltage", + "dc", + "converter", + "country", + "under_construction", + "start_date", + "geometry", + "polygon", + "contains", +] +LINE_COLUMNS = [ + "line_id", + "circuits", + "voltage", + "dc", + "p_nom_mw", + "country", + "underground", + "under_construction", + "start_date", + "geometry", + "contains", +] + + +def _empty_frame(columns: list[str], crs: str) -> gpd.GeoDataFrame: + """An empty GeoDataFrame with the given non-geometry ``columns``.""" + return gpd.GeoDataFrame( + {column: [] for column in columns}, geometry=gpd.GeoSeries([], crs=crs), crs=crs + ) + + def _create_linestring(row: pd.Series) -> LineString: """Build a LineString from a raw OSM geometry list of ``{lon, lat}`` points.""" coords = [(coord["lon"], coord["lat"]) for coord in row["geometry"]] @@ -646,19 +682,7 @@ def _aggregate_lines(df_lines: pd.DataFrame) -> pd.DataFrame: ) .reset_index() ) - return df_lines[ - [ - "line_id", - "circuits", - "voltage", - "country", - "underground", - "under_construction", - "start_date", - "geometry", - "contains", - ] - ] + return df_lines[LINE_COLUMNS] def _finalise_lines(df_lines: pd.DataFrame) -> pd.DataFrame: @@ -666,19 +690,7 @@ def _finalise_lines(df_lines: pd.DataFrame) -> pd.DataFrame: df_lines = df_lines.rename(columns={"id": "line_id", "power": "tag_type"}) df_lines["underground"] = df_lines["tag_type"] == "cable" df_lines["contains"] = df_lines["line_id"].apply(lambda x: [x.split("-")[0]]) - df_lines = df_lines[ - [ - "line_id", - "circuits", - "voltage", - "country", - "underground", - "under_construction", - "start_date", - "geometry", - "contains", - ] - ] + df_lines = df_lines[LINE_COLUMNS] df_lines["circuits"] = df_lines["circuits"].astype(int) df_lines["voltage"] = df_lines["voltage"].astype(int) return df_lines @@ -1353,7 +1365,7 @@ def clean( clean_lines = _extend_lines_to_substations(clean_lines, substation_polygons) clean_lines = gpd.GeoDataFrame(clean_lines, geometry="geometry", crs=crs) else: - clean_lines = gpd.GeoDataFrame(geometry=gpd.GeoSeries([], crs=crs), crs=crs) + clean_lines = _empty_frame([c for c in LINE_COLUMNS if c != "geometry"], crs) if "polygon" in substation_polygons.columns: substation_polygons = substation_polygons.drop(columns=["geometry"]) From b6fffa64570caa706be0f341ab6e5eca514af0d0 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:38:00 +0200 Subject: [PATCH 24/54] Add frequency-related parameters to config and schema --- config/config.schema.json | 19 +++++++++++++++++++ config/config.yaml | 4 ++++ 2 files changed, 23 insertions(+) diff --git a/config/config.schema.json b/config/config.schema.json index 9a55950..ff490ea 100644 --- a/config/config.schema.json +++ b/config/config.schema.json @@ -123,6 +123,25 @@ } } }, + "accepted_ac_frequencies_hz": { + "default": [ + 50.0, + 60.0 + ], + "description": "Frequency tag values treated as public-grid AC. A tagged value within frequency_tolerance_hz of one of these is kept and normalised to the region's AC frequency, since 50 Hz in a 60 Hz country is far likelier a tagging slip than a separate grid. Any other non-DC value is dropped: it marks a separate system such as 16.7 Hz railway traction, which must not be read as mains AC", + "items": { + "exclusiveMinimum": 0, + "type": "number" + }, + "minItems": 1, + "type": "array" + }, + "frequency_tolerance_hz": { + "default": 0.1, + "description": "Tolerance for matching a frequency tag, so that 50.0 matches 50 and 0.0 matches the DC marker", + "minimum": 0, + "type": "number" + }, "station_merge_radius_m": { "default": 500.0, "description": "Buffer radius used to merge nearby substations and line endpoints, in metres", diff --git a/config/config.yaml b/config/config.yaml index 0d043d5..6e9f263 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -23,6 +23,10 @@ network: frequency_hz: AC: 50.0 DC: 0.0 + accepted_ac_frequencies_hz: + - 50.0 + - 60.0 + frequency_tolerance_hz: 0.1 station_merge_radius_m: 500.0 remove_under_construction: true remove_after: 2026-12-31 From d1d1628a12b82f8abc7d673d33c0afde4a665b17 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:43:46 +0200 Subject: [PATCH 25/54] Bug-fix to enable customization of voltage filtering threshold --- config/config.schema.json | 6 ++++++ config/config.yaml | 1 + workflow/scripts/clean.py | 20 ++++++++++++++++---- 3 files changed, 23 insertions(+), 4 deletions(-) diff --git a/config/config.schema.json b/config/config.schema.json index ff490ea..661f0e7 100644 --- a/config/config.schema.json +++ b/config/config.schema.json @@ -105,6 +105,12 @@ "exclusiveMinimum": 0, "type": "number" }, + "minimum_voltage_dc_kv": { + "default": 150.0, + "description": "Minimum nominal DC voltage retained from OSM, in kV. Separate from the AC floor because HVDC links commonly run below 220 kV, e.g. the 150 kV Estlink 1 and Gotland links; PyPSA-Eur uses 150 kV", + "exclusiveMinimum": 0, + "type": "number" + }, "frequency_hz": { "additionalProperties": false, "description": "AC/DC frequency in Hz, used to classify and normalise OSM ``frequency`` tags.", diff --git a/config/config.yaml b/config/config.yaml index 6e9f263..e1cc31c 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -20,6 +20,7 @@ retrieve: network: include_relations: true minimum_voltage_kv: 220.0 + minimum_voltage_dc_kv: 150.0 frequency_hz: AC: 50.0 DC: 0.0 diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 5a18394..d0f0311 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -1096,6 +1096,18 @@ def _region_min_voltage( kv = override if override is not None else network["minimum_voltage_kv"] return float(kv) * 1000 # kV -> V +def _lowest_voltage_floor(network: dict[str, Any], regions: dict[str, Any]) -> float: + """The lowest voltage [V] any row could be kept at: global AC, any regional AC, or DC.""" + regional = [ + region.get("minimum_voltage_kv") + for region in regions.values() + if region.get("minimum_voltage_kv") + ] + return ( + min([network["minimum_voltage_kv"], network["minimum_voltage_dc_kv"], *regional]) + * 1000 + ) + def _format_hz(value: float) -> str: """Render a Hz value the way OSM tags whole numbers: "50", not "50.0".""" @@ -1153,7 +1165,7 @@ def clean( ``include_relations`` is off) is simply skipped. """ crs = geo_crs - min_voltage_ac = network["minimum_voltage_kv"] * 1000 # V + lowest_floor = _lowest_voltage_floor(network, regions) # V dc_hz = _format_hz(network["frequency_hz"]["DC"]) # --- Substations ------------------------------------------------- @@ -1176,7 +1188,7 @@ def clean( df_substations["start_date"] = _clean_date(df_substations["start_date"]) df_substations, list_voltages = _filter_by_voltage( - df_substations, min_voltage=min_voltage_ac + df_substations, min_voltage=lowest_floor ) df_substations["frequency"] = _clean_frequency(df_substations["frequency"]) df_substations["_ac_hz"] = df_substations["country"].map( @@ -1243,7 +1255,7 @@ def clean( df_relation["start_date"] = _clean_date(df_relation["start_date"]) df_relation["voltage"] = _clean_voltage(df_relation["voltage"]) df_relation, list_voltages = _filter_by_voltage( - df_relation, min_voltage=min_voltage_ac + df_relation, min_voltage=lowest_floor ) if not df_relation.empty: df_relation["frequency"] = _clean_frequency(df_relation["frequency"]) @@ -1332,7 +1344,7 @@ def clean( ) df_lines["start_date"] = _clean_date(df_lines["start_date"]) df_lines, list_voltages = _filter_by_voltage( - df_lines, min_voltage=min_voltage_ac + df_lines, min_voltage=lowest_floor ) if not df_lines.empty: From a403ac27d6928c89bb2c4816fbdc180afb26edca Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:46:20 +0200 Subject: [PATCH 26/54] Integrate PyPSA-Earth approach for processing DC lines --- workflow/scripts/clean.py | 46 +++++++++++++++++++++++++++++++++------ 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index d0f0311..75f3e79 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -377,8 +377,12 @@ def _filter_by_voltage( def _clean_substations( - df_substations: pd.DataFrame, list_voltages: Any, dc_hz: str + df_substations: pd.DataFrame, + list_voltages: Any, + dc_hz: str, + *, accepted_ac_hz: list[float], + tolerance_hz: float, ) -> pd.DataFrame: """Split multi-voltage substation rows and normalise each split's frequency. @@ -419,6 +423,8 @@ def _clean_substations( def _clean_lines( df_lines: pd.DataFrame, list_voltages: Any, + dc_hz: str, + *, accepted_ac_hz: list[float], tolerance_hz: float, ) -> pd.DataFrame: @@ -450,6 +456,10 @@ def _clean_lines( bool_ac = df_lines["frequency"] != dc_hz bool_dc = ~bool_ac + # Three conductors per AC circuit, two per DC circuit, as in PyPSA-Earth's + # fill_circuits. The single-cable branches below already split on this; + # the multi-voltage branches read it from here. + df_lines["_cables_per_circuit"] = np.where(bool_dc, 2, 3) bool_invalid_frequency = df_lines.apply( lambda row: row["frequency"] not in (row["_ac_hz"], dc_hz), axis=1 ) @@ -461,6 +471,18 @@ def _clean_lines( ] df_lines.loc[bool_noinfo, "cleaned"] = True + # An explicit zero means the way carries no live conductor: a ground wire + # ("ground") or a retired one ("1 disused"). Every branch below floors the + # count at one circuit, so without this the way would enter the network as + # a live single-circuit line. PyPSA-Earth reaches the same end by letting + # the zero through and dropping it in filter_circuits. + bool_no_conductors = (~df_lines["cleaned"]) & ( + ((df_lines["cables"] == "0") & (df_lines["circuits"] == "")) + | (df_lines["circuits"] == "0") + ) + df_lines.loc[bool_no_conductors, "circuits"] = "0" + df_lines.loc[bool_no_conductors, "cleaned"] = True + bool_cables_ac = ( (df_lines["cables"] != "") & (df_lines["split_elements"] == 1) @@ -1165,8 +1187,15 @@ def clean( ``include_relations`` is off) is simply skipped. """ crs = geo_crs + # Admit everything down to the lowest floor in play; each row's own + # floor (regional AC, or DC) is applied once its polarity is known. lowest_floor = _lowest_voltage_floor(network, regions) # V dc_hz = _format_hz(network["frequency_hz"]["DC"]) + dc_lines = network["dc_lines"] + frequency_options = { + "accepted_ac_hz": network["accepted_ac_frequencies_hz"], + "tolerance_hz": network["frequency_tolerance_hz"], + } # --- Substations ------------------------------------------------- logger.info("Importing substations.") @@ -1194,15 +1223,18 @@ def clean( df_substations["_ac_hz"] = df_substations["country"].map( lambda c: _region_ac_hz(c, network, regions) ) - df_substations = _clean_substations(df_substations, list_voltages, dc_hz) - - # Regional per-country minimum voltage override. - row_min = df_substations["country"].map( - lambda c: _region_min_voltage(c, network, regions) + df_substations = _clean_substations( + df_substations, list_voltages, dc_hz, **frequency_options ) + # A DC substation element only means something while DC lines are + # kept; otherwise it stays an ordinary bus, as it was before DC + # support existed. df_substations = df_substations[ - df_substations["voltage"].astype(int) >= row_min + _above_voltage_floor(df_substations, network, regions, dc_hz) ] + df_substations["dc"] = (df_substations["frequency"] == dc_hz) & ( + dc_lines == "keep" + ) # Nodes (earth-osm/Geofabrik-only, or our own overpass addition) skip # the polygon-only pipeline (touching-polygon merge, PoI, line- From 53e2c996f47328ff286ea4159eb7d1f215607d27 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:47:45 +0200 Subject: [PATCH 27/54] Integrate DC field into filtering --- workflow/scripts/clean.py | 44 +++++++++++++-------------------------- 1 file changed, 15 insertions(+), 29 deletions(-) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 75f3e79..4bed212 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -620,20 +620,19 @@ def _create_substations_poi( def _aggregate_substations(df_substations: pd.DataFrame) -> pd.DataFrame: - """One row per (original id, voltage, country), even after voltage-splitting.""" + """One row per (original id, voltage, polarity, country), even after splitting. + + Polarity is part of the key so that a converter station tagged e.g. + voltage=320000;320000 frequency=50;0 keeps both its AC and DC side. + """ df_substations = df_substations.copy() df_substations["id"] = df_substations["id"].apply( lambda x: x.split("-")[0] if "-" in x else x ) + keys = ["id", "voltage", "dc", "country"] df_substations = ( - df_substations.groupby(["id", "voltage", "country"]) - .agg( - { - col: "first" - for col in df_substations.columns - if col not in ["id", "voltage", "country"] - } - ) + df_substations.groupby(keys) + .agg({col: "first" for col in df_substations.columns if col not in keys}) .reset_index() ) return df_substations @@ -668,36 +667,27 @@ def _finalise_substations(df_substations: pd.DataFrame) -> gpd.GeoDataFrame: df_substations["contains"] = df_substations["bus_id"].apply( lambda x: x.split("-")[0] ) - columns = [ - "bus_id", - "voltage", - "country", - "under_construction", - "start_date", - "geometry", - "polygon", - "contains", - ] - df_substations = df_substations[columns] + df_substations = df_substations[SUBSTATION_COLUMNS] if not df_substations.empty: df_substations["voltage"] = df_substations["voltage"].astype(int) return df_substations def _aggregate_lines(df_lines: pd.DataFrame) -> pd.DataFrame: - """One row per (original line_id, voltage), summing circuits across splits.""" + """One row per (original line_id, voltage, polarity), summing circuits across splits.""" df_lines = df_lines.copy() df_lines["line_id"] = df_lines["line_id"].apply( lambda x: x.split("-")[0] if "-" in x else x ) + keys = ["line_id", "voltage", "dc"] df_lines = ( - df_lines.groupby(["line_id", "voltage"]) + df_lines.groupby(keys) .agg( { **{ col: "first" for col in df_lines.columns - if col not in ["line_id", "voltage", "circuits"] + if col not in [*keys, "circuits"] }, "circuits": "sum", } @@ -1388,12 +1378,8 @@ def clean( lambda c: _region_ac_hz(c, network, regions) ) df_lines = _clean_lines(df_lines, list_voltages, dc_hz, **frequency_options) - - if not df_lines.empty: - row_min = df_lines["country"].map( - lambda c: _region_min_voltage(c, network, regions) - ) - df_lines = df_lines[df_lines["voltage"].astype(int) >= row_min] + df_lines = df_lines[_above_voltage_floor(df_lines, network, regions, dc_hz)] + df_lines = _apply_dc_lines_mode(df_lines, dc_lines, dc_hz, "lines/cables") if not df_lines.empty: df_lines = _create_lines_geometry(df_lines) From 1a9c19fe85bfeac3c70da12001f3270bddf76972 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:50:25 +0200 Subject: [PATCH 28/54] Implement relationships into DC treatment --- workflow/scripts/clean.py | 110 +++++++++++++++++++++++++++++--------- 1 file changed, 86 insertions(+), 24 deletions(-) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 4bed212..bc37e7a 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -281,6 +281,26 @@ def _clean_date(column: pd.Series) -> pd.Series: return pd.to_datetime(column, errors="coerce", format="mixed") +def _clean_rating(column: pd.Series) -> pd.Series: + """Parse a ``rating`` tag into MW, summing ``;``-separated entries, as PyPSA-Eur does. + + HVDC relations carry their transfer capacity here, e.g. "1000 MW". A + value in GW is scaled to MW. Anything that yields no number is NaN. + """ + + def parse(value: str) -> float: + value = value.strip().lower() + scale = 1000.0 if "gw" in value else 1.0 + numbers = [ + float(part) + for part in re.sub(r"[^0-9.;]", "", value).split(";") + if re.fullmatch(r"\d+(\.\d+)?", part) + ] + return sum(numbers) * scale if numbers else np.nan + + return _to_str(column).map(parse).astype(float) + + def _split_cells(df: pd.DataFrame, cols: list[str] | None = None) -> pd.DataFrame: """Split semicolon-separated cells into new, identically-tagged rows.""" if cols is None: @@ -878,6 +898,46 @@ def _create_line(row: pd.Series) -> tuple[Any, list[str]]: return line, members +# Relation member roles that carry the conductor itself, as opposed to the +# terminal substations and electrodes an HVDC relation also lists. +_LINK_MEMBER_ROLES = ("line", "cable", "section") + + +def _create_single_link(row: pd.Series) -> tuple[Any, list[str]]: + """One LineString for an HVDC relation, following PyPSA-Eur's ``_create_single_link``. + + A bipole relation lists each pole as its own member, so merging all + members, as ``_create_line`` does for AC, yields parallel parts that + never form a single line, and the relation would be discarded. Instead, + keep only conductor members, keep one member per pair of endpoints so + parallel poles collapse to one, merge, and keep the longest connected + part. Every conductor member is still reported, so all of the poles' + ways are replaced by the link rather than surviving beside it. + """ + df = pd.json_normalize(row["members"]) + if "geometry" not in df.columns or "role" not in df.columns: + return linemerge([]), [] + df = df[df["role"].isin(_LINK_MEMBER_ROLES)].dropna(subset=["geometry"]) + if df.empty: + return linemerge([]), [] + df["ref"] = df["ref"].astype(str) + members = ("way/" + df["ref"]).tolist() + df["geometry"] = df.apply(_create_linestring, axis=1) + df["length"] = df["geometry"].apply(lambda line: line.length) + df["endpoints"] = df["geometry"].apply( + lambda line: tuple( + round(value, 3) + for point in sorted([line.coords[0], line.coords[-1]]) + for value in point + ) + ) + shortest = df.loc[df.groupby("endpoints")["length"].idxmin()] + link = linemerge(shortest["geometry"].tolist()) + if isinstance(link, MultiLineString): + link = max(link.geoms, key=lambda part: part.length) + return link, members + + _LINE_TAG_COLUMNS = [ "power", "cables", @@ -904,6 +964,7 @@ def _create_line(row: pd.Series) -> tuple[Any, list[str]]: "cables", "frequency", "voltage", + "rating", "construction", "construction:power", "start_date", @@ -1281,14 +1342,28 @@ def clean( ) if not df_relation.empty: df_relation["frequency"] = _clean_frequency(df_relation["frequency"]) - df_relation = df_relation[df_relation["frequency"] != dc_hz] df_relation["_ac_hz"] = df_relation["country"].map( lambda c: _region_ac_hz(c, network, regions) ) - df_relation["frequency"] = df_relation["_ac_hz"] df_relation["circuits"] = _clean_circuits(df_relation["circuits"]) df_relation["cables"] = _clean_cables(df_relation["cables"]) - df_relation = _clean_lines(df_relation, list_voltages, dc_hz) + df_relation = _clean_lines( + df_relation, list_voltages, dc_hz, **frequency_options + ) + df_relation = df_relation[ + _above_voltage_floor(df_relation, network, regions, dc_hz) + ] + # PyPSA-Eur's relationship concept: an HVDC relation is one + # link, whatever its poles and sections, rated by its tag. + # Decided before dc_lines is applied, since force_ac relabels + # the frequency but the relation is still physically a bipole. + df_relation["_hvdc"] = df_relation["frequency"] == dc_hz + df_relation["p_nom_mw"] = _clean_rating(df_relation["rating"]).where( + df_relation["_hvdc"] + ) + df_relation = _apply_dc_lines_mode( + df_relation, dc_lines, dc_hz, "route relations" + ) df_relation = df_relation.drop( columns=[ "voltage_original", @@ -1298,13 +1373,13 @@ def clean( ] ) - row_min = df_relation["country"].map( - lambda c: _region_min_voltage(c, network, regions) - ) - df_relation = df_relation[df_relation["voltage"].astype(int) >= row_min] - if not df_relation.empty: - components = df_relation.apply(_create_line, axis=1) + components = df_relation.apply( + lambda row: ( + _create_single_link(row) if row["_hvdc"] else _create_line(row) + ), + axis=1, + ) df_relation["geometry"] = components.apply(lambda x: x[0]) df_relation["contains"] = components.apply(lambda x: x[1]) @@ -1327,21 +1402,8 @@ def clean( df_relation["circuits"] = df_relation["circuits"].astype(int) df_relation["voltage"] = df_relation["voltage"].astype(int) df_relation["underground"] = False - lines_frames.append( - df_relation[ - [ - "line_id", - "circuits", - "voltage", - "country", - "underground", - "under_construction", - "start_date", - "geometry", - "contains", - ] - ] - ) + relation_lines = df_relation[LINE_COLUMNS].copy() + lines_frames.append(relation_lines) # --- AC lines/cables via individual ways --------------------------- logger.info("Importing lines and cables.") From b3b9db5e7c84d619a53e3af6cfc6ca53d7d6980a Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:51:22 +0200 Subject: [PATCH 29/54] Amend power field for HVDC --- workflow/scripts/clean.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index bc37e7a..640cfd5 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -722,6 +722,8 @@ def _finalise_lines(df_lines: pd.DataFrame) -> pd.DataFrame: df_lines = df_lines.rename(columns={"id": "line_id", "power": "tag_type"}) df_lines["underground"] = df_lines["tag_type"] == "cable" df_lines["contains"] = df_lines["line_id"].apply(lambda x: [x.split("-")[0]]) + # Capacity is tagged on HVDC relations, not on their member ways. + df_lines["p_nom_mw"] = np.nan df_lines = df_lines[LINE_COLUMNS] df_lines["circuits"] = df_lines["circuits"].astype(int) df_lines["voltage"] = df_lines["voltage"].astype(int) From 2ba595d967aa95d15dc53aba5d6c83fdea3772c3 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:52:30 +0200 Subject: [PATCH 30/54] Define function for force_ac treatment --- workflow/scripts/clean.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 640cfd5..02da17e 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -272,6 +272,30 @@ def _frequency_for_split(row: pd.Series) -> str: return values[min(position, len(values) - 1)] +def _apply_dc_lines_mode( + df: pd.DataFrame, dc_lines: str, dc_hz: str, label: str +) -> pd.DataFrame: + """Set the ``dc`` flag from the normalised frequency and apply ``network.dc_lines``.""" + df = df.copy() + is_dc = df["frequency"] == dc_hz + if dc_lines == "drop": + logger.info("Dropped %d DC %s (network.dc_lines: drop).", int(is_dc.sum()), label) + df = df[~is_dc].copy() + df["dc"] = False + elif dc_lines == "force_ac": + logger.info( + "Relabelled %d DC %s as AC (network.dc_lines: force_ac).", + int(is_dc.sum()), + label, + ) + df.loc[is_dc, "frequency"] = df.loc[is_dc, "_ac_hz"] + df["dc"] = False + else: + logger.info("Keeping %d DC %s.", int(is_dc.sum()), label) + df["dc"] = is_dc + return df + + def _clean_date(column: pd.Series) -> pd.Series: """Parse a raw ``start_date`` tag column to datetimes, coercing invalid values to NaT.""" column = _apply_corrections(_to_str(column), _TAG_CORRECTIONS["date"]) From 23b278846895d9ddd936631dec0ee547aeaa687a Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:53:02 +0200 Subject: [PATCH 31/54] Improve voltage filtering --- workflow/scripts/clean.py | 34 ++++++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 02da17e..badb8ae 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -1190,10 +1190,36 @@ def _split_country_codes(country: Any) -> list[str]: def _region_min_voltage( country: str, network: dict[str, Any], regions: dict[str, Any] ) -> float: - """Minimum AC voltage [V] for ``country``: its regional override, else the network default.""" - override = regions.get(country, {}).get("minimum_voltage_kv") - kv = override if override is not None else network["minimum_voltage_kv"] - return float(kv) * 1000 # kV -> V + """Minimum AC voltage [V] for ``country``: its regional override, else the network default. + + For a cross-border element the most permissive threshold wins, so an + interconnector survives whenever either side would keep it. + """ + thresholds = [ + regions.get(code, {}).get("minimum_voltage_kv") or network["minimum_voltage_kv"] + for code in _split_country_codes(country) + ] or [network["minimum_voltage_kv"]] + return float(min(thresholds)) * 1000 # kV -> V + + +def _above_voltage_floor( + df: pd.DataFrame, + network: dict[str, Any], + regions: dict[str, Any], + dc_hz: str, +) -> pd.Series: + """True where a row meets its own floor: the DC floor for DC, the regional AC floor otherwise. + + Rows must already carry a normalised ``frequency``. This runs after an + initial filter at the lowest floor in play (see ``_lowest_voltage_floor``), + so that neither a lower regional AC floor nor the DC floor is cut off + by the global AC one first. + """ + ac_floor = df["country"].map(lambda c: _region_min_voltage(c, network, regions)) + dc_floor = network["minimum_voltage_dc_kv"] * 1000 + floor = ac_floor.where(df["frequency"] != dc_hz, dc_floor) + return df["voltage"].astype(int) >= floor + def _lowest_voltage_floor(network: dict[str, Any], regions: dict[str, Any]) -> float: """The lowest voltage [V] any row could be kept at: global AC, any regional AC, or DC.""" From 4055f763268b7d0929eba58584dc8a39950af157 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:54:00 +0200 Subject: [PATCH 32/54] Improve treatment of empty geometries --- workflow/scripts/clean.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index badb8ae..bee2dba 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -1308,8 +1308,13 @@ def clean( inputs.get("substations_relation", []), ) if df_substations.empty: - empty_buses = gpd.GeoDataFrame(geometry=gpd.GeoSeries([], crs=crs), crs=crs) - empty_polygons = gpd.GeoDataFrame(geometry=gpd.GeoSeries([], crs=crs), crs=crs) + # Same columns as the populated path: build_network selects polygons + # by bus_id before it checks for empty input, so a geometry-only + # frame would fail there for a country with no substations. + empty_buses = _empty_frame( + [c for c in SUBSTATION_COLUMNS if c not in ("geometry", "polygon")], crs + ) + empty_polygons = _empty_frame(["bus_id", "voltage"], crs) else: df_substations["voltage"] = _clean_voltage(df_substations["voltage"]) df_substations["under_construction"] = ( From 922f240e22163cc2acc0694caf98b5037d0ecc4e Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:57:45 +0200 Subject: [PATCH 33/54] Add cables-to-underground transformation --- workflow/scripts/clean.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index bee2dba..8e732f8 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -1469,6 +1469,15 @@ def clean( ) if not df_lines.empty: df_lines = _drop_duplicate_lines(df_lines) + # A relation is assumed to be underground when all its member ways + # are cables (originates from PyPSA-Eur) + # TODO: check applicability on the global scale + if lines_frames: + way_power = df_lines.set_index("id")["power"] + relation_lines["underground"] = relation_lines["contains"].apply( + lambda ways: set(way_power.reindex(ways).dropna()) == {"cable"} + ) + lines_frames[0] = relation_lines len_before = len(df_lines) df_lines = df_lines[~df_lines["id"].isin(ways_to_replace)] logger.info( From 035520edb2a11ddc2e5cc41a13a5a020e9201cd6 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 17:57:54 +0200 Subject: [PATCH 34/54] Minor polishes --- workflow/scripts/clean.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 8e732f8..a637c21 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -21,6 +21,7 @@ import itertools import json import logging +import re from pathlib import Path from typing import TYPE_CHECKING, Any @@ -1324,6 +1325,9 @@ def clean( ) df_substations["start_date"] = _clean_date(df_substations["start_date"]) + df_substations["converter"] = ( + _to_str(df_substations["substation"]).str.contains("converter") + ) df_substations, list_voltages = _filter_by_voltage( df_substations, min_voltage=lowest_floor ) @@ -1458,6 +1462,7 @@ def clean( df_relation = df_relation.rename(columns={"id": "line_id"}) df_relation["circuits"] = df_relation["circuits"].astype(int) df_relation["voltage"] = df_relation["voltage"].astype(int) + # Set from the member ways' own power tag once they are read. df_relation["underground"] = False relation_lines = df_relation[LINE_COLUMNS].copy() lines_frames.append(relation_lines) From 5dec475b1178b1d791ebd1597ee3c93ed7ce3530 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:02:41 +0200 Subject: [PATCH 35/54] Integrate DC treatment into network build --- config/config.schema.json | 10 +++ config/config.yaml | 1 + workflow/scripts/build_network.py | 121 +++++++++++++++++++++++------- 3 files changed, 105 insertions(+), 27 deletions(-) diff --git a/config/config.schema.json b/config/config.schema.json index 661f0e7..cbfc29e 100644 --- a/config/config.schema.json +++ b/config/config.schema.json @@ -148,6 +148,16 @@ "minimum": 0, "type": "number" }, + "dc_lines": { + "default": "keep", + "description": "Treatment of DC lines and cables. keep carries them as dc=true, gives them their own buses and adds converters where a station holds both AC and DC buses. drop removes them. force_ac keeps them but relabels them AC, for downstream models that cannot handle DC", + "enum": [ + "keep", + "drop", + "force_ac" + ], + "type": "string" + }, "station_merge_radius_m": { "default": 500.0, "description": "Buffer radius used to merge nearby substations and line endpoints, in metres", diff --git a/config/config.yaml b/config/config.yaml index e1cc31c..cde76e9 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -28,6 +28,7 @@ network: - 50.0 - 60.0 frequency_tolerance_hz: 0.1 + dc_lines: keep station_merge_radius_m: 500.0 remove_under_construction: true remove_after: 2026-12-31 diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 78b5cfb..06b8c7b 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -42,6 +42,7 @@ "bus_id", "station_id", "voltage_kv", + "dc", "country", "under_construction", "start_date", @@ -54,6 +55,7 @@ "bus0", "bus1", "voltage_kv", + "dc", "circuits", "length_m", "underground", @@ -74,6 +76,22 @@ "geometry", ] +# A converter joins a DC bus (bus0) to an AC bus (bus1). ``pairing`` records +# which rule found the AC side: "same_station" (PyPSA-Earth's get_converters) +# or "nearest_station" (PyPSA-Eur's converter-hall mapping). ``p_nom_mw`` is +# the summed rating of the HVDC links on the DC bus, as in PyPSA-Eur. +CONVERTER_COLUMNS = [ + "converter_id", + "station_id", + "bus0", + "bus1", + "voltage_bus0_kv", + "voltage_bus1_kv", + "p_nom_mw", + "pairing", + "geometry", +] + STATION_POLYGON_COLUMNS = ["station_id", "geometry"] @@ -123,11 +141,11 @@ def _treat_under_construction( def _merge_identical_lines(lines: gpd.GeoDataFrame) -> gpd.GeoDataFrame: - """Aggregate lines with identical geometry and voltage (e.g. duplicated across a border).""" + """Aggregate lines with identical geometry, voltage and polarity (e.g. duplicated across a border).""" lines_all = lines.copy() lines_to_drop = [] - for _, group in lines_all.groupby(["geometry", "voltage"]): + for _, group in lines_all.groupby(["geometry", "voltage", "dc"]): line_ids = list(group["line_id"]) if len(line_ids) > 1: lid_old = line_ids[0] @@ -182,6 +200,7 @@ def _add_line_endings(lines: gpd.GeoDataFrame) -> pd.DataFrame: line_data = lines[ [ "voltage", + "dc", "geometry", "line_id", "country", @@ -224,8 +243,10 @@ def create_bus_data(group: pd.DataFrame) -> pd.Series: } ) + # Polarity is part of the key: an AC and a DC line ending at the same + # point at the same voltage must not share a bus. endpoints = ( - endpoints.groupby(["voltage", "geometry"]) + endpoints.groupby(["voltage", "dc", "geometry"]) .apply(create_bus_data, include_groups=False) .reset_index() ) @@ -239,6 +260,7 @@ def create_bus_data(group: pd.DataFrame) -> pd.Series: [ "bus_id", "voltage", + "dc", "geometry", "contains", "country", @@ -292,6 +314,9 @@ def split_overpassing_lines( continue nearby_buses = buses_epsgmod.iloc[possible_matches] + # A line can only connect to buses of its own polarity, so splitting + # it at a bus of the other one would only strand a stub there. + nearby_buses = nearby_buses[nearby_buses["dc"] == lines.at[line_index, "dc"]] bus_in_tol = nearby_buses[nearby_buses.geometry.distance(line_geom) <= tol] endpoint0 = line_geom.boundary.geoms[0] @@ -354,15 +379,16 @@ def _create_merge_mapping( buses_virtual = gpd.sjoin( buses_virtual, - lines[["line_id", "geometry", "voltage", "circuits"]], + lines[["line_id", "geometry", "voltage", "dc", "circuits"]], how="left", predicate="touches", ) buses_virtual = buses_virtual[ - buses_virtual["voltage_left"] == buses_virtual["voltage_right"] + (buses_virtual["voltage_left"] == buses_virtual["voltage_right"]) + & (buses_virtual["dc_left"] == buses_virtual["dc_right"]) ] - buses_virtual = buses_virtual.drop(columns=["voltage_right"]).rename( - columns={"voltage_left": "voltage"} + buses_virtual = buses_virtual.drop(columns=["voltage_right", "dc_right"]).rename( + columns={"voltage_left": "voltage", "dc_left": "dc"} ) counts = ( @@ -395,6 +421,7 @@ def _create_merge_mapping( [ "line_id", "voltage", + "dc", "circuits", "length", "geometry", @@ -423,6 +450,9 @@ def _create_merge_mapping( first_node = next(iter(component)) circuits = graph.nodes[first_node].get("circuits") voltage = graph.nodes[first_node].get("voltage") + dc = bool(graph.nodes[first_node].get("dc")) + # Segments of one rated link share its rating; take it from any + # segment that carries one. geometry = linemerge( [graph.nodes[node].get("geometry") for node in subgraph.nodes()] ) @@ -454,6 +484,7 @@ def _create_merge_mapping( "line_id": f"merged_{node_longest}+{len(contains_lines) - 1}", "circuits": circuits, "voltage": voltage, + "dc": dc, "geometry": geometry, "underground": underground, "country": country, @@ -468,6 +499,7 @@ def _create_merge_mapping( "line_id", "circuits", "voltage", + "dc", "geometry", "underground", "country", @@ -597,39 +629,50 @@ def _merge_buses_to_stations( stations_all["polygon"] = stations_all["geometry"].copy() buses_all = gpd.sjoin(buses_all, stations_all, how="left", predicate="within") - buses_all = buses_all.drop_duplicates(subset=["station_id", "voltage"]) + buses_all = buses_all.drop_duplicates(subset=["station_id", "voltage", "dc"]) offset = 15 # metres geo_to_dist = Transformer.from_crs(geo_crs, distance_crs, always_xy=True) dist_to_geo = Transformer.from_crs(distance_crs, geo_crs, always_xy=True) + def level_bus_id(station_id: str, voltage: float, dc: bool) -> str: + # AC keeps the existing "{station}-{kV}" form; DC gets a suffix so it + # cannot collide with an AC bus of the same voltage. + return f"{station_id}-{int(voltage / 1000)}" + ("-dc" if dc else "") + for station_id, group in buses_all.groupby("station_id"): - voltages = sorted(group["voltage"].unique(), reverse=True) + # One bus per (voltage, polarity) level, highest voltage first and + # AC before DC, which reproduces the previous order exactly when a + # station has no DC side. + levels = sorted( + set(zip(group["voltage"], group["dc"])), key=lambda lv: (-lv[0], lv[1]) + ) not_virtual = ~group.bus_id.str.startswith("virtual_") - if len(voltages) > 1: + if len(levels) > 1: poi_x, poi_y = geo_to_dist.transform( group["poi"].values[0].x, group["poi"].values[0].y ) - for idx, voltage in enumerate(voltages): + for idx, (voltage, dc) in enumerate(levels): poi_x_offset = poi_x + offset * np.sin( - np.pi / 4 + 2 * np.pi * idx / len(voltages) + np.pi / 4 + 2 * np.pi * idx / len(levels) ).round(4) poi_y_offset = poi_y + offset * np.cos( - np.pi / 4 + 2 * np.pi * idx / len(voltages) + np.pi / 4 + 2 * np.pi * idx / len(levels) ).round(4) poi_offset = Point(dist_to_geo.transform(poi_x_offset, poi_y_offset)) - group.loc[(group["voltage"] == voltage) & not_virtual, "bus_id"] = ( - station_id + "-" + str(int(voltage / 1000)) + at_level = (group["voltage"] == voltage) & (group["dc"] == dc) + group.loc[at_level & not_virtual, "bus_id"] = level_bus_id( + station_id, voltage, dc ) - group.loc[group["voltage"] == voltage, "geometry"] = poi_offset + group.loc[at_level, "geometry"] = poi_offset buses_all.loc[group.index, "bus_id"] = group["bus_id"] buses_all.loc[group.index, "geometry"] = group["geometry"] else: - voltage = voltages[0] - buses_all.loc[group.loc[not_virtual].index, "bus_id"] = ( - station_id + "-" + str(int(voltage / 1000)) + voltage, dc = levels[0] + buses_all.loc[group.loc[not_virtual].index, "bus_id"] = level_bus_id( + station_id, voltage, dc ) buses_all.loc[group.index, "geometry"] = group["poi"] @@ -673,14 +716,19 @@ def _map_endpoints_to_buses( lines_all = connection.copy().set_index(id_col) for coord in range(2): - endpoints = lines_all[["voltage", "geometry"]].copy() + endpoints = lines_all[["voltage", "dc", "geometry"]].copy() endpoints["geometry"] = get_point( endpoints.geometry.apply(_remove_loops_from_multiline), -1 * coord ) endpoints = gpd.sjoin(endpoints, buses_all, how="left", predicate="intersects") - endpoints = endpoints[endpoints["voltage_left"] == endpoints["voltage_right"]] - endpoints = endpoints.drop(columns=["voltage_right"]).rename( - columns={"voltage_left": "voltage"} + # A line attaches only to a bus of its own voltage and polarity, as + # in PyPSA-Earth's set_lines_ids, which groups by (voltage, dc). + endpoints = endpoints[ + (endpoints["voltage_left"] == endpoints["voltage_right"]) + & (endpoints["dc_left"] == endpoints["dc_right"]) + ] + endpoints = endpoints.drop(columns=["voltage_right", "dc_right"]).rename( + columns={"voltage_left": "voltage", "dc_left": "dc"} ) lines_all[f"poi_perimeter{coord}"] = endpoints["poi_perimeter"] @@ -781,8 +829,13 @@ def _extend_lines_to_buses( def _add_transformers(buses: gpd.GeoDataFrame, geo_crs: str) -> gpd.GeoDataFrame: - """All-pairs transformers between voltage-level buses of the same real station.""" - buses_all = buses.copy().set_index("bus_id") + """All-pairs transformers between the AC voltage-level buses of each station. + + DC buses are excluded: the AC and DC sides of a station are joined by + converters (see ``_add_converters``), never by a transformer, as in + PyPSA-Earth's get_transformers. + """ + buses_all = buses[~buses["dc"]].copy().set_index("bus_id") columns = ["bus0", "bus1", "voltage_bus0", "voltage_bus1", "station_id", "geometry"] all_transformers = gpd.GeoDataFrame( columns=[*columns, "transformer_id"], crs=geo_crs @@ -845,8 +898,22 @@ def build_network( ``station_id``) and ``buses_polygon`` (the substation polygons scoped to the buses that made it into the output, keyed by ``bus_id``). """ + # Inputs written before DC support existed carry no dc column; they + # were AC-only by construction, so default to AC rather than fail. + substations = substations.copy() + lines = lines.copy() + for frame in (substations, lines): + if "dc" not in frame.columns: + frame["dc"] = False + frame["dc"] = frame["dc"].fillna(False).astype(bool) + if "converter" not in substations.columns: + substations["converter"] = False + substations["converter"] = substations["converter"].fillna(False).astype(bool) + if "p_nom_mw" not in lines.columns: + lines["p_nom_mw"] = np.nan + buses = _treat_under_construction( - substations.copy(), remove_under_construction, remove_after + substations, remove_under_construction, remove_after ) buses_polygon = substations_polygon[ @@ -1032,7 +1099,7 @@ def _write_components( snakemake.params.station_merge_radius_m, ) logger.info( - "Built %d buses, %d lines, and %d transformers.", + "Built %d buses, %d lines (%d DC), %d transformers, and %d converters.", len(buses), len(lines), len(transformers), From 57aae7336a8ba88e5986ec656a9a15b89943c34b Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:06:12 +0200 Subject: [PATCH 36/54] Implement adding converters --- config/config.schema.json | 6 ++ config/config.yaml | 1 + workflow/rules/network.smk | 3 + workflow/scripts/build_network.py | 162 +++++++++++++++++++++++++++++- 4 files changed, 169 insertions(+), 3 deletions(-) diff --git a/config/config.schema.json b/config/config.schema.json index cbfc29e..69c230c 100644 --- a/config/config.schema.json +++ b/config/config.schema.json @@ -148,6 +148,12 @@ "minimum": 0, "type": "number" }, + "converter_search_radius_m": { + "default": 50000.0, + "description": "How far, in metres, a station tagged substation=converter looks for an AC station when it has DC buses but no AC bus of its own. Converter halls often sit apart from the AC substation they feed, further than station_merge_radius_m merges; PyPSA-Eur uses 50 km", + "exclusiveMinimum": 0, + "type": "number" + }, "dc_lines": { "default": "keep", "description": "Treatment of DC lines and cables. keep carries them as dc=true, gives them their own buses and adds converters where a station holds both AC and DC buses. drop removes them. force_ac keeps them but relabels them AC, for downstream models that cannot handle DC", diff --git a/config/config.yaml b/config/config.yaml index cde76e9..04592a9 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -28,6 +28,7 @@ network: - 50.0 - 60.0 frequency_tolerance_hz: 0.1 + converter_search_radius_m: 50000.0 dc_lines: keep station_merge_radius_m: 500.0 remove_under_construction: true diff --git a/workflow/rules/network.smk b/workflow/rules/network.smk index 957512a..9fa1441 100644 --- a/workflow/rules/network.smk +++ b/workflow/rules/network.smk @@ -86,9 +86,11 @@ rule build_network: buses="/build/csv/buses.csv", lines="/build/csv/lines.csv", transformers="/build/csv/transformers.csv", + converters="/build/csv/converters.csv", buses_geojson="/build/geojson/buses.geojson", lines_geojson="/build/geojson/lines.geojson", transformers_geojson="/build/geojson/transformers.geojson", + converters_geojson="/build/geojson/converters.geojson", stations_polygon="/build/geojson/stations_polygon.geojson", buses_polygon="/build/geojson/buses_polygon.geojson", log: @@ -98,6 +100,7 @@ rule build_network: threads: 1 params: station_merge_radius_m=config["network"]["station_merge_radius_m"], + converter_search_radius_m=config["network"]["converter_search_radius_m"], remove_under_construction=config["network"]["remove_under_construction"], remove_after=config["network"]["remove_after"], crs=config["crs"].model_dump(mode="json"), diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 06b8c7b..fcb66ca 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -875,6 +875,104 @@ def _add_transformers(buses: gpd.GeoDataFrame, geo_crs: str) -> gpd.GeoDataFrame return all_transformers[["transformer_id", *columns]] +def _add_converters( + buses: gpd.GeoDataFrame, + lines: gpd.GeoDataFrame, + converter_stations: set[str], + search_radius_m: float | None, + distance_crs: str, + geo_crs: str, +) -> gpd.GeoDataFrame: + """Join each DC bus to the AC grid, combining PyPSA-Earth's and PyPSA-Eur's rules. + + - PyPSA-Earth (``get_converters``): a station holding both AC and DC + buses is a converter station, and each DC bus pairs with the AC bus + of that station whose voltage is numerically closest. + - PyPSA-Eur (``_add_dc_buses``): a converter hall often stands apart + from the AC substation it feeds, further than station merging + reaches. A station holding an OSM ``substation=converter`` but no AC + bus therefore pairs its DC bus with the highest-voltage bus of the + nearest AC station within ``search_radius_m``. Requiring the tag + keeps a link that is merely cut at a border, or a cable-to-overhead + transition, from being wired to whatever AC happens to be near. + + Each converter's ``p_nom_mw`` is the summed rating of the HVDC links + on its DC bus, as in PyPSA-Eur, and NaN when none is rated. + """ + columns = [ + "converter_id", + "station_id", + "bus0", + "bus1", + "voltage_bus0", + "voltage_bus1", + "p_nom_mw", + "pairing", + ] + ac_all = buses[~buses["dc"]] + ac_points = ac_all.geometry.to_crs(distance_crs) + + def nearest_ac_bus(dc_bus: pd.Series) -> pd.Series | None: + if ac_all.empty or search_radius_m is None: + return None + origin = gpd.GeoSeries([dc_bus.geometry], crs=geo_crs).to_crs(distance_crs) + distances = ac_points.distance(origin.iloc[0]) + distances = distances[distances <= search_radius_m] + if distances.empty: + return None + station = ac_all.loc[distances.idxmin(), "station_id"] + at_station = ac_all[ac_all["station_id"] == station] + return at_station.loc[at_station["voltage"].idxmax()] + + def link_rating(bus_id: str) -> float: + attached = lines[(lines["bus0"] == bus_id) | (lines["bus1"] == bus_id)] + return attached["p_nom_mw"].sum(min_count=1) + + records = [] + for station_id, group in buses.groupby("station_id"): + dc_buses = group[group["dc"]] + if dc_buses.empty: + continue + ac_buses = group[~group["dc"]] + for _, dc_bus in dc_buses.sort_values("voltage", ascending=False).iterrows(): + if not ac_buses.empty: + ac_bus = ac_buses.loc[ + (ac_buses["voltage"] - dc_bus["voltage"]).abs().idxmin() + ] + pairing = "same_station" + elif station_id in converter_stations: + ac_bus = nearest_ac_bus(dc_bus) + if ac_bus is None: + logger.info( + "Converter station %s has no AC station within %s m.", + station_id, + search_radius_m, + ) + continue + pairing = "nearest_station" + else: + continue + records.append( + { + "converter_id": ( + f"{station_id}-{int(dc_bus['voltage'] / 1000)}dc" + f"-{int(ac_bus['voltage'] / 1000)}" + ), + "station_id": station_id, + "bus0": dc_bus["bus_id"], + "bus1": ac_bus["bus_id"], + "voltage_bus0": dc_bus["voltage"], + "voltage_bus1": ac_bus["voltage"], + "p_nom_mw": link_rating(dc_bus["bus_id"]), + "pairing": pairing, + "geometry": LineString([dc_bus.geometry, ac_bus.geometry]), + } + ) + if not records: + return _empty_geodataframe(columns, crs=geo_crs) + return gpd.GeoDataFrame(records, geometry="geometry", crs=geo_crs) + + def build_network( substations: gpd.GeoDataFrame, substations_polygon: gpd.GeoDataFrame, @@ -890,8 +988,14 @@ def build_network( gpd.GeoDataFrame, gpd.GeoDataFrame, gpd.GeoDataFrame, + gpd.GeoDataFrame, ]: - """Create buses, AC lines, and transformers from clean's output. + """Create buses, lines, transformers and converters from clean's output. + + Lines and buses carry a ``dc`` flag. DC gets its own buses at each + station, lines attach only to buses of their own polarity, transformers + join AC levels only, and converters join the DC side of a station to + its AC side. Also returns two polygon views for visualisation: ``stations_polygon`` (the clustered station shapes from station-seed buffering, keyed by @@ -978,8 +1082,23 @@ def build_network( geo_crs=geo_crs, tol=station_merge_radius_m, ) + # Stations holding an OSM converter hall, found before buses are merged + # to one per level, which would keep only one substation's tags. + converter_halls = buses[buses["converter"].fillna(False).astype(bool)] + converter_stations = set( + gpd.sjoin( + converter_halls[["geometry"]], + stations[["station_id", "geometry"]], + predicate="intersects", + )["station_id"] + ) + buses = _merge_buses_to_stations( - buses, stations, distance_crs=distance_crs, geo_crs=geo_crs + buses, + stations, + distance_crs=distance_crs, + geo_crs=geo_crs, + max_station_voltage_ratio=max_station_voltage_ratio, ) buses["geometry"] = gpd.points_from_xy( @@ -1015,6 +1134,14 @@ def build_network( buses = buses[~bool_not_connected].reset_index(drop=True) transformers = _add_transformers(buses, geo_crs=geo_crs) + converters = _add_converters( + buses, + lines, + converter_stations, + converter_search_radius_m, + distance_crs=distance_crs, + geo_crs=geo_crs, + ) lines["length"] = lines.to_crs(distance_crs).length @@ -1059,6 +1186,20 @@ def _contains_to_osm_ids(value: Any) -> str: _non_geometry(TRANSFORMER_COLUMNS), crs=geo_crs ) + if not converters.empty: + converters_out = converters.copy() + for side in ("0", "1"): + converters_out[f"voltage_bus{side}_kv"] = ( + converters_out[f"voltage_bus{side}"] / 1000 + ).astype(int) + converters_out = gpd.GeoDataFrame( + converters_out[CONVERTER_COLUMNS], geometry="geometry", crs=geo_crs + ) + else: + converters_out = _empty_geodataframe( + _non_geometry(CONVERTER_COLUMNS), crs=geo_crs + ) + stations_polygon_out = stations[STATION_POLYGON_COLUMNS].copy() buses_polygon_out = buses_polygon.copy() @@ -1066,6 +1207,7 @@ def _contains_to_osm_ids(value: Any) -> str: buses_out, lines_out, transformers_out, + converters_out, stations_polygon_out, buses_polygon_out, ) @@ -1088,7 +1230,14 @@ def _write_components( snakemake = mock_snakemake("build_network") configure_logging(snakemake.log[0]) - buses, lines, transformers, stations_polygon, buses_polygon = build_network( + ( + buses, + lines, + transformers, + converters, + stations_polygon, + buses_polygon, + ) = build_network( gpd.read_file(snakemake.input.substations), gpd.read_file(snakemake.input.substations_polygon), gpd.read_file(snakemake.input.lines), @@ -1097,12 +1246,16 @@ def _write_components( snakemake.params.crs["geo"], snakemake.params.crs["distance"], snakemake.params.station_merge_radius_m, + snakemake.params.max_station_voltage_ratio, + snakemake.params.converter_search_radius_m, ) logger.info( "Built %d buses, %d lines (%d DC), %d transformers, and %d converters.", len(buses), len(lines), + int(lines["dc"].sum()) if not lines.empty else 0, len(transformers), + len(converters), ) _write_components(buses, snakemake.output.buses, snakemake.output.buses_geojson) _write_components(lines, snakemake.output.lines, snakemake.output.lines_geojson) @@ -1111,5 +1264,8 @@ def _write_components( snakemake.output.transformers, snakemake.output.transformers_geojson, ) + _write_components( + converters, snakemake.output.converters, snakemake.output.converters_geojson + ) stations_polygon.to_file(snakemake.output.stations_polygon, driver="GeoJSON") buses_polygon.to_file(snakemake.output.buses_polygon, driver="GeoJSON") From 8b46d6c9b8c07cdee2376fbdbccecac7e14a4fd6 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:09:03 +0200 Subject: [PATCH 37/54] Add DC part into script validation schema --- workflow/scripts/_schema.py | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/workflow/scripts/_schema.py b/workflow/scripts/_schema.py index f55baaf..df1f72a 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -12,7 +12,14 @@ from typing import Any, Literal from earth_osm.regions import get_all_valid_codes, get_region_tuple -from pydantic import BaseModel, ConfigDict, Field, ValidationError, field_validator +from pydantic import ( + BaseModel, + ConfigDict, + Field, + PositiveFloat, + ValidationError, + field_validator, +) from ruamel.yaml import YAML from ruamel.yaml.comments import CommentedMap @@ -185,6 +192,26 @@ class NetworkConfig(ConfigModel): ), ge=0, ) + converter_search_radius_m: float = Field( + 50000.0, + description=( + "How far, in metres, a station tagged substation=converter looks " + "for an AC station when it has DC buses but no AC bus of its own. " + "Converter halls often sit apart from the AC substation they feed, " + "further than station_merge_radius_m merges; PyPSA-Eur uses 50 km" + ), + gt=0, + ) + dc_lines: Literal["keep", "drop", "force_ac"] = Field( + "keep", + description=( + "Treatment of DC lines and cables. keep carries them as dc=true, " + "gives them their own buses and adds converters where a station " + "holds both AC and DC buses. drop removes them. force_ac keeps " + "them but relabels them AC, for downstream models that cannot " + "handle DC" + ), + ) station_merge_radius_m: float = Field( 500.0, description="Buffer radius used to merge nearby substations and line endpoints, in metres", From 92dea8c0dcb706644d0b02055fc9a86a91dc9dcf Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:09:31 +0200 Subject: [PATCH 38/54] Add voltage threshold into script validation schema --- workflow/scripts/_schema.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/workflow/scripts/_schema.py b/workflow/scripts/_schema.py index df1f72a..27bb3b9 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -168,6 +168,15 @@ class NetworkConfig(ConfigModel): minimum_voltage_kv: float = Field( 220.0, description="Minimum nominal AC voltage retained from OSM, in kV", gt=0 ) + minimum_voltage_dc_kv: float = Field( + 150.0, + description=( + "Minimum nominal DC voltage retained from OSM, in kV. Separate from " + "the AC floor because HVDC links commonly run below 220 kV, e.g. " + "the 150 kV Estlink 1 and Gotland links; PyPSA-Eur uses 150 kV" + ), + gt=0, + ) frequency_hz: FrequencyConfig = Field( default_factory=FrequencyConfig, description="AC/DC frequency in Hz; override per country in config/regions for e.g. 60 Hz grids", From f04e049b81fbe837153ffee4b7b3d5f0bf9b948a Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:09:52 +0200 Subject: [PATCH 39/54] Add converters to plotting --- workflow/rules/plot.smk | 1 + 1 file changed, 1 insertion(+) diff --git a/workflow/rules/plot.smk b/workflow/rules/plot.smk index 3c06d71..83e8c6d 100644 --- a/workflow/rules/plot.smk +++ b/workflow/rules/plot.smk @@ -10,6 +10,7 @@ rule build_interactive_map: transformers=rules.build_network.output.transformers_geojson, stations_polygon=rules.build_network.output.stations_polygon, buses_polygon=rules.build_network.output.buses_polygon, + converters=rules.build_network.output.converters_geojson, output: map="/map.html", log: From 9d4bf911de9024a5f556a4b9417868e437634887 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:14:20 +0200 Subject: [PATCH 40/54] Document DC and frequency treatment --- INTERFACE.yaml | 2 +- README.md | 2 +- config/README.md | 8 ++++++++ 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/INTERFACE.yaml b/INTERFACE.yaml index 03794ca..1163b76 100644 --- a/INTERFACE.yaml +++ b/INTERFACE.yaml @@ -9,7 +9,7 @@ pathvars: description: >- Raw OSM files are written to retrieve/{country}_{feature}.json. Clean features are written below clean; generic network outputs - below build (buses/lines/transformers as CSV under build/csv + below build (buses/lines/transformers/converters as CSV under build/csv and GeoJSON, including substation polygons, under build/geojson). An interactive PyDeck map of the network, and the default target of this workflow, is written to map.html. diff --git a/README.md b/README.md index 54832ce..9f10ec9 100644 --- a/README.md +++ b/README.md @@ -40,7 +40,7 @@ file per country and feature (`lines_way`, `cables_way`, `substations_way`, `substations_node`, `substations_relation`, `routes_relation`). Both retrieval backends write the same raw-Overpass-JSON shape, so downstream cleaning doesn't need to know which one ran. Clean features use `/clean/*.geojson`; -generic network components use `/build/csv/{buses,lines,transformers}.csv` +generic network components use `/build/csv/{buses,lines,transformers,converters}.csv` and matching GeoJSON files under `/build/geojson/`, which also includes `stations_polygon.geojson` (clustered station shapes) and `buses_polygon.geojson` (substation polygons scoped to the buses in the output). diff --git a/config/README.md b/config/README.md index 4acce0e..07fd48f 100644 --- a/config/README.md +++ b/config/README.md @@ -6,6 +6,14 @@ A parameter `network.include_relations` defines whether the network should consi `network.station_merge_radius_m` is a buffer radius with the merge distance being *twice* as high.E.g.the default value of 500 m merges substations up to one kilometre apart. +### Frequencies and DC lines + +Frequency tags are matched numerically within `network.frequency_tolerance_hz`, so `50.0` counts as 50 Hz and `0.0` as DC. Values in `network.accepted_ac_frequencies_hz` are treated as frequencies of a public-grid alternated current (AC) and normalised to the region's nominal AC frequency. Any other value, such as 16.7 Hz corresponding railway traction, is dropped, since it belongs to a separate grid. Where a line lists several circuits, e.g. `voltage=380000;110000` with `frequency=50;16.7`, each frequency is paired with the voltage in the same position. + +`network.dc_lines` controls DC lines and cables. With `keep`, the default, they carry `dc: true` and get their own buses at each station. `drop` removes DC lines, and `force_ac` keeps them relabelled as AC. DC has its own voltage floor, `network.minimum_voltage_dc_kv`, since HVDC links often run below the AC floor. + +HVDC links combine two approaches. As in PyPSA-Eur, a DC `route=power` relation becomes a single link: parallel poles collapse into one line, the member ways are replaced by it, and its `rating` tag is kept as `p_nom_mw`. DC ways outside any relation are kept too, as in PyPSA-Earth. Converters are written to `converters.csv` by two rules. A station holding both AC and DC buses pairs each DC bus with its AC bus of the closest voltage, as in PyPSA-Earth. A station tagged `substation=converter` with no AC bus of its own pairs with the highest-voltage bus of the nearest AC station within `network.converter_search_radius_m`, as in PyPSA-Eur. The `pairing` column records which rule applied. + The [BE+NL example](./examples/config.BE-NL.yaml) is a small European development scope. Country files under `config/regions` are intentionally small defaults for now; community-maintained local corrections belong there rather than in workflow code. ### Adding custom data From a29aaa501c1afb861cf7f3e7c3cd271a9be51066 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:17:37 +0200 Subject: [PATCH 41/54] Fetch power field for DC lines --- workflow/scripts/build_network.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index fcb66ca..b7d5747 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -56,6 +56,7 @@ "bus1", "voltage_kv", "dc", + "p_nom_mw", "circuits", "length_m", "underground", @@ -422,6 +423,7 @@ def _create_merge_mapping( "line_id", "voltage", "dc", + "p_nom_mw", "circuits", "length", "geometry", @@ -453,6 +455,10 @@ def _create_merge_mapping( dc = bool(graph.nodes[first_node].get("dc")) # Segments of one rated link share its rating; take it from any # segment that carries one. + p_nom_mw = pd.Series( + [graph.nodes[node].get("p_nom_mw") for node in subgraph.nodes()], + dtype=float, + ).max() geometry = linemerge( [graph.nodes[node].get("geometry") for node in subgraph.nodes()] ) @@ -485,6 +491,7 @@ def _create_merge_mapping( "circuits": circuits, "voltage": voltage, "dc": dc, + "p_nom_mw": p_nom_mw, "geometry": geometry, "underground": underground, "country": country, @@ -500,6 +507,7 @@ def _create_merge_mapping( "circuits", "voltage", "dc", + "p_nom_mw", "geometry", "underground", "country", From ae2c1f01bd907ba6da8251fb70fc33820c3b13dd Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:20:57 +0200 Subject: [PATCH 42/54] Add processing of empty converters --- workflow/scripts/build_network.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index b7d5747..bf0031a 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -1044,6 +1044,9 @@ def build_network( empty_transformers = _empty_geodataframe( _non_geometry(TRANSFORMER_COLUMNS), crs=geo_crs ) + empty_converters = _empty_geodataframe( + _non_geometry(CONVERTER_COLUMNS), crs=geo_crs + ) empty_stations_polygon = _empty_geodataframe( _non_geometry(STATION_POLYGON_COLUMNS), crs=geo_crs ) @@ -1051,6 +1054,7 @@ def build_network( empty_buses, empty_lines, empty_transformers, + empty_converters, empty_stations_polygon, buses_polygon, ) From 1098f5061ac5c94b0a720869b654a8e320ece6fb Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:21:20 +0200 Subject: [PATCH 43/54] Improve clarity of docstring --- workflow/scripts/build_network.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index bf0031a..fda342d 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -104,7 +104,8 @@ def _empty_geodataframe(columns: list[str], crs: str) -> gpd.GeoDataFrame: def _non_geometry(columns: list[str]) -> list[str]: - """Drop the geometry column, which ``_empty_geodataframe`` supplies itself.""" + """Drop the geometry column from ``columns`` maintaining structure + otherwise""" return [column for column in columns if column != "geometry"] From 765b537f4d2c864628137149e12844e63d2cfcc7 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:21:57 +0200 Subject: [PATCH 44/54] Add converters to plotting --- workflow/scripts/build_interactive_map.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/workflow/scripts/build_interactive_map.py b/workflow/scripts/build_interactive_map.py index 9a72622..45b14c3 100644 --- a/workflow/scripts/build_interactive_map.py +++ b/workflow/scripts/build_interactive_map.py @@ -190,6 +190,7 @@ def build_map( transformers: gpd.GeoDataFrame, stations: gpd.GeoDataFrame, bus_polygons: gpd.GeoDataFrame, + converters: gpd.GeoDataFrame | None = None, *, geo_crs: str, distance_crs: str, @@ -198,7 +199,7 @@ def build_map( lines_simplify_m: float | None, coord_decimals: int, ) -> pdk.Deck: - """Create the generic AC map without requiring DC links or converters.""" + """Create the network map; ``converters`` is optional so AC-only callers can omit it.""" lines = lines.copy() if not lines.empty: lines["color"] = line_colors(lines["voltage_kv"]) @@ -242,6 +243,14 @@ def build_map( distance_crs=distance_crs, coord_decimals=coord_decimals, ), + path_layer( + converters if converters is not None else gpd.GeoDataFrame(), + "Converters", + [0, 255, 255, 220], + geo_crs=geo_crs, + distance_crs=distance_crs, + coord_decimals=coord_decimals, + ), ) if layer ] From 00b7bc94b1494896245fd70899772d00e20b7773 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:28:52 +0200 Subject: [PATCH 45/54] Minor clean-up --- workflow/scripts/build_network.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index fda342d..6d445ce 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -630,7 +630,10 @@ def _create_station_seeds( def _merge_buses_to_stations( - buses: gpd.GeoDataFrame, stations: gpd.GeoDataFrame, distance_crs: str, geo_crs: str + buses: gpd.GeoDataFrame, + stations: gpd.GeoDataFrame, + distance_crs: str, + geo_crs: str, ) -> gpd.GeoDataFrame: """Keep one bus per (station, voltage); offset multi-voltage stations for visual clarity.""" buses_all = buses.copy().reset_index(drop=True) @@ -991,6 +994,7 @@ def build_network( geo_crs: str, distance_crs: str, station_merge_radius_m: float = BUS_TOL, + converter_search_radius_m: float | None = None, ) -> tuple[ gpd.GeoDataFrame, gpd.GeoDataFrame, @@ -1111,7 +1115,6 @@ def build_network( stations, distance_crs=distance_crs, geo_crs=geo_crs, - max_station_voltage_ratio=max_station_voltage_ratio, ) buses["geometry"] = gpd.points_from_xy( @@ -1259,7 +1262,6 @@ def _write_components( snakemake.params.crs["geo"], snakemake.params.crs["distance"], snakemake.params.station_merge_radius_m, - snakemake.params.max_station_voltage_ratio, snakemake.params.converter_search_radius_m, ) logger.info( From d659cc1ca7416a3011fea3d73ed36a6faa0d687c Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 18:29:43 +0200 Subject: [PATCH 46/54] Add details to description of buffer radius --- workflow/scripts/_schema.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/workflow/scripts/_schema.py b/workflow/scripts/_schema.py index 27bb3b9..1c4bf3c 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -223,7 +223,13 @@ class NetworkConfig(ConfigModel): ) station_merge_radius_m: float = Field( 500.0, - description="Buffer radius used to merge nearby substations and line endpoints, in metres", + description=( + "Buffer radius used to merge nearby substations and line endpoints, " + "in metres. Both sides of a pair are buffered by this radius before " + "the buffers are dissolved, so two elements merge once they are " + "within twice this value of each other: the default of 500 merges " + "substations up to 1 km apart" + ), gt=0, ) remove_under_construction: bool = Field( From 6cdbb66a4b3e3de5419f0251f92e362cb3dcbb72 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 21:11:12 +0200 Subject: [PATCH 47/54] Fix hardcoding --- config/config.schema.json | 18 ++++++++++++++++ config/config.yaml | 3 +++ workflow/rules/network.smk | 2 ++ workflow/scripts/_helpers.py | 4 ++-- workflow/scripts/_schema.py | 25 +++++++++++++++++++++++ workflow/scripts/build_network.py | 24 ++++++++++++++-------- workflow/scripts/clean.py | 10 ++++++--- workflow/scripts/retrieve_osm_overpass.py | 10 ++++++--- 8 files changed, 80 insertions(+), 16 deletions(-) diff --git a/config/config.schema.json b/config/config.schema.json index 69c230c..0ec9cb2 100644 --- a/config/config.schema.json +++ b/config/config.schema.json @@ -65,6 +65,12 @@ "exclusiveMinimum": 0, "type": "integer" }, + "backoff_factor": { + "default": 2.0, + "description": "Exponential backoff between retries, in seconds: the n-th retry waits backoff_factor * 2**(n-1)", + "minimum": 0, + "type": "number" + }, "user_agent": { "additionalProperties": false, "description": "Identifies this tool to the Overpass API, per its fair-use policy.", @@ -170,6 +176,18 @@ "exclusiveMinimum": 0, "type": "number" }, + "station_bus_offset_m": { + "default": 15.0, + "description": "Distance, in metres, by which the buses of a multi-voltage station are spread around its centre so each level stays distinguishable on a map; 0 places them all at the centre", + "minimum": 0, + "type": "number" + }, + "overpassing_lines_tolerance_m": { + "default": 1.0, + "description": "A line passing within this distance, in metres, of a bus it does not end at is split there and connected to it", + "exclusiveMinimum": 0, + "type": "number" + }, "remove_under_construction": { "default": true, "description": "Whether assets tagged as under construction are dropped", diff --git a/config/config.yaml b/config/config.yaml index 04592a9..86d6c58 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -12,6 +12,7 @@ retrieve: url: "https://overpass-api.de/api/interpreter" max_tries: 5 timeout: 600 + backoff_factor: 2.0 user_agent: project_name: "grid-builder" email: "" @@ -31,6 +32,8 @@ network: converter_search_radius_m: 50000.0 dc_lines: keep station_merge_radius_m: 500.0 + station_bus_offset_m: 15.0 + overpassing_lines_tolerance_m: 1.0 remove_under_construction: true remove_after: 2026-12-31 diff --git a/workflow/rules/network.smk b/workflow/rules/network.smk index 9fa1441..7d37160 100644 --- a/workflow/rules/network.smk +++ b/workflow/rules/network.smk @@ -100,6 +100,8 @@ rule build_network: threads: 1 params: station_merge_radius_m=config["network"]["station_merge_radius_m"], + station_bus_offset_m=config["network"]["station_bus_offset_m"], + overpassing_lines_tolerance_m=config["network"]["overpassing_lines_tolerance_m"], converter_search_radius_m=config["network"]["converter_search_radius_m"], remove_under_construction=config["network"]["remove_under_construction"], remove_after=config["network"]["remove_after"], diff --git a/workflow/scripts/_helpers.py b/workflow/scripts/_helpers.py index 0ca446b..81cd932 100644 --- a/workflow/scripts/_helpers.py +++ b/workflow/scripts/_helpers.py @@ -26,8 +26,8 @@ logger = logging.getLogger(__name__) -GEO_CRS = "EPSG:4326" -BUS_TOL = 500 # metres; default station merge tolerance +GEO_CRS: str = _CONFIG_DEFAULTS["crs"]["geo"] +BUS_TOL: float = _CONFIG_DEFAULTS["network"]["station_merge_radius_m"] # metres def configure_logging(log_path: str) -> None: diff --git a/workflow/scripts/_schema.py b/workflow/scripts/_schema.py index 1c4bf3c..bc23fd6 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -109,6 +109,14 @@ class OverpassApiConfig(ConfigModel): 5, description="Maximum number of attempts per query before giving up", ge=1 ) timeout: int = Field(600, description="Per-request timeout in seconds", gt=0) + backoff_factor: float = Field( + 2.0, + description=( + "Exponential backoff between retries, in seconds: the n-th retry " + "waits backoff_factor * 2**(n-1)" + ), + ge=0, + ) user_agent: OverpassUserAgentConfig = Field(default_factory=OverpassUserAgentConfig) @@ -232,6 +240,23 @@ class NetworkConfig(ConfigModel): ), gt=0, ) + station_bus_offset_m: float = Field( + 15.0, + description=( + "Distance, in metres, by which the buses of a multi-voltage station " + "are spread around its centre so each level stays distinguishable " + "on a map; 0 places them all at the centre" + ), + ge=0, + ) + overpassing_lines_tolerance_m: float = Field( + 1.0, + description=( + "A line passing within this distance, in metres, of a bus it does " + "not end at is split there and connected to it" + ), + gt=0, + ) remove_under_construction: bool = Field( True, description="Whether assets tagged as under construction are dropped" ) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 6d445ce..14b5e20 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -294,14 +294,14 @@ def _alpha_suffix(i: int) -> str: def split_overpassing_lines( - lines: gpd.GeoDataFrame, buses: gpd.GeoDataFrame, distance_crs: str, tol: float = 1 + lines: gpd.GeoDataFrame, buses: gpd.GeoDataFrame, distance_crs: str, tol: float ) -> gpd.GeoDataFrame: """Split a line at any bus it geometrically overpasses without a shared OSM node.""" lines = lines.copy() lines_to_add = [] lines_to_split = [] - high_voltage_lines = lines.query("voltage >= 220000") + high_voltage_lines = lines if high_voltage_lines.empty: return lines @@ -634,6 +634,7 @@ def _merge_buses_to_stations( stations: gpd.GeoDataFrame, distance_crs: str, geo_crs: str, + offset_m: float, ) -> gpd.GeoDataFrame: """Keep one bus per (station, voltage); offset multi-voltage stations for visual clarity.""" buses_all = buses.copy().reset_index(drop=True) @@ -643,7 +644,6 @@ def _merge_buses_to_stations( buses_all = gpd.sjoin(buses_all, stations_all, how="left", predicate="within") buses_all = buses_all.drop_duplicates(subset=["station_id", "voltage", "dc"]) - offset = 15 # metres geo_to_dist = Transformer.from_crs(geo_crs, distance_crs, always_xy=True) dist_to_geo = Transformer.from_crs(distance_crs, geo_crs, always_xy=True) @@ -665,10 +665,10 @@ def level_bus_id(station_id: str, voltage: float, dc: bool) -> str: group["poi"].values[0].x, group["poi"].values[0].y ) for idx, (voltage, dc) in enumerate(levels): - poi_x_offset = poi_x + offset * np.sin( + poi_x_offset = poi_x + offset_m * np.sin( np.pi / 4 + 2 * np.pi * idx / len(levels) ).round(4) - poi_y_offset = poi_y + offset * np.cos( + poi_y_offset = poi_y + offset_m * np.cos( np.pi / 4 + 2 * np.pi * idx / len(levels) ).round(4) poi_offset = Point(dist_to_geo.transform(poi_x_offset, poi_y_offset)) @@ -993,7 +993,10 @@ def build_network( remove_after: str | None, geo_crs: str, distance_crs: str, + *, station_merge_radius_m: float = BUS_TOL, + station_bus_offset_m: float, + overpassing_lines_tolerance_m: float, converter_search_radius_m: float | None = None, ) -> tuple[ gpd.GeoDataFrame, @@ -1077,7 +1080,9 @@ def build_network( buses_line_endings = _add_line_endings(lines) buses = pd.concat([buses, buses_line_endings], ignore_index=True) - lines = split_overpassing_lines(lines, buses, distance_crs=distance_crs) + lines = split_overpassing_lines( + lines, buses, distance_crs=distance_crs, tol=overpassing_lines_tolerance_m + ) bool_virtual = buses["bus_id"].str.startswith("virtual") buses = buses[~bool_virtual] @@ -1115,6 +1120,7 @@ def build_network( stations, distance_crs=distance_crs, geo_crs=geo_crs, + offset_m=station_bus_offset_m, ) buses["geometry"] = gpd.points_from_xy( @@ -1261,8 +1267,10 @@ def _write_components( snakemake.params.remove_after, snakemake.params.crs["geo"], snakemake.params.crs["distance"], - snakemake.params.station_merge_radius_m, - snakemake.params.converter_search_radius_m, + station_merge_radius_m=snakemake.params.station_merge_radius_m, + station_bus_offset_m=snakemake.params.station_bus_offset_m, + overpassing_lines_tolerance_m=snakemake.params.overpassing_lines_tolerance_m, + converter_search_radius_m=snakemake.params.converter_search_radius_m, ) logger.info( "Built %d buses, %d lines (%d DC), %d transformers, and %d converters.", diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index a637c21..be939a5 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -390,7 +390,7 @@ def _drop_duplicate_lines(df_lines: pd.DataFrame) -> pd.DataFrame: def _filter_by_voltage( - df: pd.DataFrame, min_voltage: float = 220000 + df: pd.DataFrame, min_voltage: float ) -> tuple[pd.DataFrame, Any]: """Keep only rows at or above ``min_voltage`` [V]; return the surviving voltage set too.""" if df.empty: @@ -1360,7 +1360,9 @@ def clean( df_way = df_way.drop(columns=["is_node"]) df_way = _create_substations_geometry(df_way) df_way = _merge_touching_polygons(df_way, crs=crs) - df_way = _create_substations_poi(df_way) + df_way = _create_substations_poi( + df_way, tol=network["station_merge_radius_m"] / 2 + ) if not df_node.empty: df_node = df_node.drop(columns=["is_node"]) @@ -1525,7 +1527,9 @@ def clean( clean_lines = gpd.GeoDataFrame(df_lines_all, geometry="geometry", crs=crs) clean_lines = _remove_lines_within_substations(clean_lines, substation_polygons) if not clean_lines.empty and not substation_polygons.empty: - clean_lines = _extend_lines_to_substations(clean_lines, substation_polygons) + clean_lines = _extend_lines_to_substations( + clean_lines, substation_polygons, tol=network["station_merge_radius_m"] / 2 + ) clean_lines = gpd.GeoDataFrame(clean_lines, geometry="geometry", crs=crs) else: clean_lines = _empty_frame([c for c in LINE_COLUMNS if c != "geometry"], crs) diff --git a/workflow/scripts/retrieve_osm_overpass.py b/workflow/scripts/retrieve_osm_overpass.py index 6495dde..3d1acf7 100644 --- a/workflow/scripts/retrieve_osm_overpass.py +++ b/workflow/scripts/retrieve_osm_overpass.py @@ -112,13 +112,15 @@ def _normalise_node_geometry(payload: dict[str, Any]) -> dict[str, Any]: return {**payload, "elements": elements} -def _session(max_tries: int, user_agent: str) -> requests.Session: +def _session( + max_tries: int, backoff_factor: float, user_agent: str +) -> requests.Session: """Build a ``requests`` session that retries transient errors with backoff.""" session = requests.Session() session.headers.update({"User-Agent": user_agent}) retry = Retry( total=max_tries, - backoff_factor=2, + backoff_factor=backoff_factor, status_forcelist=(429, 500, 502, 503, 504), allowed_methods=False, # type: ignore[arg-type] # retry POST too; urllib3 stubs miss this value respect_retry_after_header=True, @@ -134,6 +136,7 @@ def retrieve_from_overpass( include_relations: bool, url: str, max_tries: int, + backoff_factor: float, timeout: int, user_agent: str, ) -> dict[str, dict[str, Any]]: @@ -149,7 +152,7 @@ def retrieve_from_overpass( if not include_relations: features.remove("routes_relation") - session = _session(max_tries, user_agent) + session = _session(max_tries, backoff_factor, user_agent) payloads: dict[str, dict[str, Any]] = {"routes_relation": {"elements": []}} for feature in features: logger.info("Querying Overpass for %s in %s", feature, iso_code) @@ -194,6 +197,7 @@ def retrieve_from_overpass( include_relations, url=overpass_api["url"], max_tries=overpass_api["max_tries"], + backoff_factor=overpass_api["backoff_factor"], timeout=overpass_api["timeout"], user_agent=user_agent, ) From 43950b9c524bb3ee4d480f3b689f1a37f0741fe7 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 21:30:07 +0200 Subject: [PATCH 48/54] Minor style fixes --- workflow/scripts/_helpers.py | 16 +++++++++++++--- workflow/scripts/_schema.py | 8 +------- workflow/scripts/build_interactive_map.py | 7 ++----- workflow/scripts/build_network.py | 3 +-- workflow/scripts/clean.py | 2 -- 5 files changed, 17 insertions(+), 19 deletions(-) diff --git a/workflow/scripts/_helpers.py b/workflow/scripts/_helpers.py index 81cd932..c438fa5 100644 --- a/workflow/scripts/_helpers.py +++ b/workflow/scripts/_helpers.py @@ -26,9 +26,6 @@ logger = logging.getLogger(__name__) -GEO_CRS: str = _CONFIG_DEFAULTS["crs"]["geo"] -BUS_TOL: float = _CONFIG_DEFAULTS["network"]["station_merge_radius_m"] # metres - def configure_logging(log_path: str) -> None: """Send rule and dependency logging to the Snakemake log file.""" @@ -59,6 +56,19 @@ def load_internal_yaml(filename: str) -> Any: return yaml.safe_load(handle) +def load_config_defaults() -> Any: + """Load config/config.yaml, the defaults generated from the schema in _schema.py.""" + path = Path(__file__).resolve().parents[2] / "config" / "config.yaml" + with open(path) as handle: + return yaml.safe_load(handle) + + +# Defaults for functions called outside Snakemake, whose rules pass their own +# resolved config values instead. +GEO_CRS: str = load_config_defaults()["crs"]["geo"] +BUS_TOL: float = load_config_defaults()["network"]["station_merge_radius_m"] # metres + + def mock_snakemake( rulename: str, root_dir: str | Path | None = None, diff --git a/workflow/scripts/_schema.py b/workflow/scripts/_schema.py index bc23fd6..e931120 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -231,13 +231,7 @@ class NetworkConfig(ConfigModel): ) station_merge_radius_m: float = Field( 500.0, - description=( - "Buffer radius used to merge nearby substations and line endpoints, " - "in metres. Both sides of a pair are buffered by this radius before " - "the buffers are dissolved, so two elements merge once they are " - "within twice this value of each other: the default of 500 merges " - "substations up to 1 km apart" - ), + description="Buffer radius used to merge nearby substations and line endpoints, in metres", gt=0, ) station_bus_offset_m: float = Field( diff --git a/workflow/scripts/build_interactive_map.py b/workflow/scripts/build_interactive_map.py index 45b14c3..20e28ef 100644 --- a/workflow/scripts/build_interactive_map.py +++ b/workflow/scripts/build_interactive_map.py @@ -92,7 +92,6 @@ def path_layer( frame: gpd.GeoDataFrame, name: str, color: list[int] | str, - *, geo_crs: str, distance_crs: str, coord_decimals: int, @@ -138,7 +137,6 @@ def polygon_layer( frame: gpd.GeoDataFrame, name: str, color: list[int], - *, geo_crs: str, distance_crs: str, coord_decimals: int, @@ -190,8 +188,7 @@ def build_map( transformers: gpd.GeoDataFrame, stations: gpd.GeoDataFrame, bus_polygons: gpd.GeoDataFrame, - converters: gpd.GeoDataFrame | None = None, - *, + converters: gpd.GeoDataFrame | None, geo_crs: str, distance_crs: str, stations_simplify_m: float | None, @@ -199,7 +196,7 @@ def build_map( lines_simplify_m: float | None, coord_decimals: int, ) -> pdk.Deck: - """Create the network map; ``converters`` is optional so AC-only callers can omit it.""" + """Create an interactive network map""" lines = lines.copy() if not lines.empty: lines["color"] = line_colors(lines["voltage_kv"]) diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 14b5e20..796f0e3 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -993,10 +993,9 @@ def build_network( remove_after: str | None, geo_crs: str, distance_crs: str, - *, - station_merge_radius_m: float = BUS_TOL, station_bus_offset_m: float, overpassing_lines_tolerance_m: float, + station_merge_radius_m: float = BUS_TOL, converter_search_radius_m: float | None = None, ) -> tuple[ gpd.GeoDataFrame, diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index be939a5..39e4b64 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -425,7 +425,6 @@ def _clean_substations( df_substations: pd.DataFrame, list_voltages: Any, dc_hz: str, - *, accepted_ac_hz: list[float], tolerance_hz: float, ) -> pd.DataFrame: @@ -469,7 +468,6 @@ def _clean_lines( df_lines: pd.DataFrame, list_voltages: Any, dc_hz: str, - *, accepted_ac_hz: list[float], tolerance_hz: float, ) -> pd.DataFrame: From 2fe33e866fd7c7b36cfb57c394860100a28e30be Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 21:33:57 +0200 Subject: [PATCH 49/54] Apply formatting fixes --- workflow/scripts/build_interactive_map.py | 12 +++++-- workflow/scripts/build_network.py | 44 ++++++++++------------- workflow/scripts/clean.py | 33 +++++++++-------- 3 files changed, 45 insertions(+), 44 deletions(-) diff --git a/workflow/scripts/build_interactive_map.py b/workflow/scripts/build_interactive_map.py index 20e28ef..29f90c3 100644 --- a/workflow/scripts/build_interactive_map.py +++ b/workflow/scripts/build_interactive_map.py @@ -196,7 +196,7 @@ def build_map( lines_simplify_m: float | None, coord_decimals: int, ) -> pdk.Deck: - """Create an interactive network map""" + """Create an interactive network map.""" lines = lines.copy() if not lines.empty: lines["color"] = line_colors(lines["voltage_kv"]) @@ -1094,8 +1094,16 @@ def compress_json_in_script(match: re.Match[str]) -> str: configure_logging(snakemake.log[0]) geo_crs = snakemake.params.crs["geo"] simplify = snakemake.params.interactive_map["simplify_geometries"] + buses, lines, transformers, stations, bus_polygons, converters = ( + gpd.read_file(path).to_crs(geo_crs) for path in snakemake.input + ) deck = build_map( - *(gpd.read_file(path).to_crs(geo_crs) for path in snakemake.input), + buses, + lines, + transformers, + stations, + bus_polygons, + converters, geo_crs=geo_crs, distance_crs=snakemake.params.crs["distance"], stations_simplify_m=simplify["stations_m"] if simplify["enable"] else None, diff --git a/workflow/scripts/build_network.py b/workflow/scripts/build_network.py index 796f0e3..ba0f3ff 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -104,14 +104,14 @@ def _empty_geodataframe(columns: list[str], crs: str) -> gpd.GeoDataFrame: def _non_geometry(columns: list[str]) -> list[str]: - """Drop the geometry column from ``columns`` maintaining structure - otherwise""" + """Drop the geometry column from ``columns`` maintaining structure otherwise.""" return [column for column in columns if column != "geometry"] def _merge_country_codes(values: Any) -> str: - """Clean-up country codes for multy-country entries which is essential - for cross-border elements. + """Clean-up country codes for multy-country entries. + + This is essential for cross-border elements. """ co_codes: set[str] = set() for value in values: @@ -194,8 +194,7 @@ def _remove_loops_from_multiline(multiline: Any) -> Any: def _add_line_endings(lines: gpd.GeoDataFrame) -> pd.DataFrame: - """Create deterministic virtual buses at each unique (voltage, endpoint) - combination. + """Create deterministic virtual buses at each unique (voltage, endpoint) pair. A virtual bus inherits its attributes from the lines that meet there. """ @@ -1251,25 +1250,20 @@ def _write_components( snakemake = mock_snakemake("build_network") configure_logging(snakemake.log[0]) - ( - buses, - lines, - transformers, - converters, - stations_polygon, - buses_polygon, - ) = build_network( - gpd.read_file(snakemake.input.substations), - gpd.read_file(snakemake.input.substations_polygon), - gpd.read_file(snakemake.input.lines), - snakemake.params.remove_under_construction, - snakemake.params.remove_after, - snakemake.params.crs["geo"], - snakemake.params.crs["distance"], - station_merge_radius_m=snakemake.params.station_merge_radius_m, - station_bus_offset_m=snakemake.params.station_bus_offset_m, - overpassing_lines_tolerance_m=snakemake.params.overpassing_lines_tolerance_m, - converter_search_radius_m=snakemake.params.converter_search_radius_m, + (buses, lines, transformers, converters, stations_polygon, buses_polygon) = ( + build_network( + gpd.read_file(snakemake.input.substations), + gpd.read_file(snakemake.input.substations_polygon), + gpd.read_file(snakemake.input.lines), + snakemake.params.remove_under_construction, + snakemake.params.remove_after, + snakemake.params.crs["geo"], + snakemake.params.crs["distance"], + station_merge_radius_m=snakemake.params.station_merge_radius_m, + station_bus_offset_m=snakemake.params.station_bus_offset_m, + overpassing_lines_tolerance_m=snakemake.params.overpassing_lines_tolerance_m, + converter_search_radius_m=snakemake.params.converter_search_radius_m, + ) ) logger.info( "Built %d buses, %d lines (%d DC), %d transformers, and %d converters.", diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 39e4b64..3950a23 100644 --- a/workflow/scripts/clean.py +++ b/workflow/scripts/clean.py @@ -267,9 +267,7 @@ def _frequency_for_split(row: pd.Series) -> str: values = row["frequency"].split(";") if len(values) == 1: return values[0] - position = ( - int(row["id"].rsplit("-", 1)[1]) - 1 if row["split_elements"] > 1 else 0 - ) + position = int(row["id"].rsplit("-", 1)[1]) - 1 if row["split_elements"] > 1 else 0 return values[min(position, len(values) - 1)] @@ -280,7 +278,9 @@ def _apply_dc_lines_mode( df = df.copy() is_dc = df["frequency"] == dc_hz if dc_lines == "drop": - logger.info("Dropped %d DC %s (network.dc_lines: drop).", int(is_dc.sum()), label) + logger.info( + "Dropped %d DC %s (network.dc_lines: drop).", int(is_dc.sum()), label + ) df = df[~is_dc].copy() df["dc"] = False elif dc_lines == "force_ac": @@ -1202,10 +1202,7 @@ def _region_min_voltage( def _above_voltage_floor( - df: pd.DataFrame, - network: dict[str, Any], - regions: dict[str, Any], - dc_hz: str, + df: pd.DataFrame, network: dict[str, Any], regions: dict[str, Any], dc_hz: str ) -> pd.Series: """True where a row meets its own floor: the DC floor for DC, the regional AC floor otherwise. @@ -1228,7 +1225,9 @@ def _lowest_voltage_floor(network: dict[str, Any], regions: dict[str, Any]) -> f if region.get("minimum_voltage_kv") ] return ( - min([network["minimum_voltage_kv"], network["minimum_voltage_dc_kv"], *regional]) + min( + [network["minimum_voltage_kv"], network["minimum_voltage_dc_kv"], *regional] + ) * 1000 ) @@ -1323,9 +1322,9 @@ def clean( ) df_substations["start_date"] = _clean_date(df_substations["start_date"]) - df_substations["converter"] = ( - _to_str(df_substations["substation"]).str.contains("converter") - ) + df_substations["converter"] = _to_str( + df_substations["substation"] + ).str.contains("converter") df_substations, list_voltages = _filter_by_voltage( df_substations, min_voltage=lowest_floor ) @@ -1474,7 +1473,7 @@ def clean( ) if not df_lines.empty: df_lines = _drop_duplicate_lines(df_lines) - # A relation is assumed to be underground when all its member ways + # A relation is assumed to be underground when all its member ways # are cables (originates from PyPSA-Eur) # TODO: check applicability on the global scale if lines_frames: @@ -1498,9 +1497,7 @@ def clean( | (df_lines["power"] == "construction") ) df_lines["start_date"] = _clean_date(df_lines["start_date"]) - df_lines, list_voltages = _filter_by_voltage( - df_lines, min_voltage=lowest_floor - ) + df_lines, list_voltages = _filter_by_voltage(df_lines, min_voltage=lowest_floor) if not df_lines.empty: df_lines["circuits"] = _clean_circuits(df_lines["circuits"]) @@ -1526,7 +1523,9 @@ def clean( clean_lines = _remove_lines_within_substations(clean_lines, substation_polygons) if not clean_lines.empty and not substation_polygons.empty: clean_lines = _extend_lines_to_substations( - clean_lines, substation_polygons, tol=network["station_merge_radius_m"] / 2 + clean_lines, + substation_polygons, + tol=network["station_merge_radius_m"] / 2, ) clean_lines = gpd.GeoDataFrame(clean_lines, geometry="geometry", crs=crs) else: From 14caf00a9a93e16d3577f99ed2618d1765f2d3fb Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 26 Sep 2026 21:51:46 +0200 Subject: [PATCH 50/54] Revise README --- README.md | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 9f10ec9..175bffe 100644 --- a/README.md +++ b/README.md @@ -12,9 +12,11 @@ A modular Snakemake workflow for building a model of transmission power grid for ## About -`grid-builder` is a modular `snakemake` workflow that retrieves OpenStreetMap power infrastructure and builds a generic high-voltage network. Custom data can be injected by providing input files of unified structure. The `grid-builder` workflow can be imported into another `snakemake`based project. +`grid-builder` is a modular `snakemake` workflow that retrieves [**OpenStreetMap**](https://osm.org) power infrastructure and builds a generic high-voltage network. Custom data can be injected by providing input files of unified structure. The `grid-builder` workflow can be imported into another `snakemake`based project. -The workflow retains AC substations, overhead lines, and cables at configured voltage levels, then creates generic buses, connected line segments, and voltage-pair transformers. The outputs preserve OSM provenance and geometry but contain no PyPSA-specific line types, capacities, or electrical-component assumptions. +You can help to impove quality of [**OpenStreetMap data**] by joining the [**MapYourGrid initiative**](https://mapyourgrid.org). Many resources as [video tutorials](https://www.youtube.com/channel/UC52jOcw_6_7iTMW-lXwLrQQ) or [starter-kit](https://mapyourgrid.org/starter-kit/) help to improve open data that is used by GridBuilder to build the grid topology. + +The GridBuilder workflow retains AC substations, overhead lines, and cables at configured voltage levels, then creates generic buses, connected line segments, and voltage-pair transformers. The outputs preserve OSM provenance and geometry but contain no PyPSA-specific line types, capacities, or electrical-component assumptions. Buses and lines carry the country they belong to, along with their construction status and planned start date. Country information is essential to resolve assign a correct line types which is strongly regional-specific. @@ -26,7 +28,7 @@ Currently implemented: 1. Retrieve OSM substations, lines, cables, and (optionally) circuit relations by country, either from a cached local Geofabrik PBF extract or the live Overpass API. 2. Clean the raw retrieval output, filtering voltage, frequency, construction status, and future assets, and grouping relation member ways into one line per real-world circuit. -3. Merge nearby stations and line endpoints into generic buses, AC lines, and transformers. +3. Merge nearby stations and line endpoints into generic buses, AC and DC lines, and transformers. 4. Build an interactive map of the resulting network (`map.html`). ## Configuration @@ -105,6 +107,17 @@ rebuild the installed environment with `pixi reinstall --locked`. ## References & related work +GridBuilder is built on top of other initiatives which created and improved open power infrastructure data, developed and ways to integrate those data into energy modelling workflows. + +The list bellow is contains a non-exaustive list of links and references, and can be absolutely expanded and improved. + +### Open data and open source projects + +* [OpenStreetMap](https://osmfoundation.org/) initiative which is the biggest crowd-sourced database of geospatial information +* [MapYourGrid](https://mapyourgrid.org/) initiative that empower individuals, communities and nations around the world to map the electrical grid + +### Academic publications + * Jonas Hörsch et al. 2018. PyPSA-Eur: An open optimisation model of the European transmission system, *Energy Strategy Reviews*, Volume 22. https://doi.org/10.1016/j.esr.2018.08.012 * Maximilian Parzen et al. 2023. PyPSA-Earth: A new global open energy system optimization model demonstrated in Africa, *Applied Energy*, Volume 341. https://doi.org/10.1016/j.apenergy.2023.121096 * Bobby Xiong et al. 2025. Modelling the high-voltage grid using open data for Europe and beyond. *Sci Data* 12, 277. https://doi.org/10.1038/s41597-025-04550-7 From 78fac480a08417166b76070aff7ecfd0579ea6d3 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 27 Sep 2026 10:58:35 +0200 Subject: [PATCH 51/54] Add minimal testing config for Colombia --- config.CO.yaml | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 config.CO.yaml diff --git a/config.CO.yaml b/config.CO.yaml new file mode 100644 index 0000000..53621f0 --- /dev/null +++ b/config.CO.yaml @@ -0,0 +1,5 @@ +# Colombia run used for the main vs integrate-composite-features comparison. +# Everything else is left at the shipped defaults, so the two branches differ +# only by code. Run with: +# snakemake --cores 2 --configfile config.CO.yaml +countries: [CO] From fdcc6cd88d59d6b17f84b952b43b2cc4ec4a80bf Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 27 Sep 2026 10:59:00 +0200 Subject: [PATCH 52/54] Enhance test config --- tests/test_config.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/test_config.py b/tests/test_config.py index 05be5e2..d9f10b1 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -32,6 +32,9 @@ {"network": {"remove_under_construction": "not-a-bool"}}, {"crs": {"typo": True}}, {"regions": {"BE": {"frequency_hz": {"typo": True}}}}, + {"custom_data": {"typo": True}}, + {"custom_data": {"files": ["data/lines.json"]}}, + {"custom_data": {"files": ["data/BE_not_a_feature.json"]}}, ], ) def test_invalid_config(config): @@ -40,6 +43,18 @@ def test_invalid_config(config): validate_config(config) +def test_custom_data_accepts_retrieval_style_filenames(): + """A custom file named like a retrieved one validates. + + ``clean`` reads the country back out of the filename, so that naming is + what lets a custom file need no special handling downstream. + """ + config = validate_config( + {"custom_data": {"files": ["data/custom/BE_lines_way.json"]}} + ) + assert config.custom_data.files == ["data/custom/BE_lines_way.json"] + + def test_generated_config_matches_repository(tmp_path): """Keep the shipped defaults and JSON schema in sync with validation.""" config_dir = Path(__file__).resolve().parents[1] / "config" From 07d7d24d6eb427d6bc3bba060754ccac178dcded Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 27 Sep 2026 10:59:59 +0200 Subject: [PATCH 53/54] Sketch integration of newly introduced features --- tests/test_network_processing.py | 670 ++++++++++++++++++++++++++++++- 1 file changed, 660 insertions(+), 10 deletions(-) diff --git a/tests/test_network_processing.py b/tests/test_network_processing.py index e896207..de59911 100644 --- a/tests/test_network_processing.py +++ b/tests/test_network_processing.py @@ -1,16 +1,25 @@ """Focused tests for generic OSM cleaning and topology construction.""" import json +import logging import geopandas as gpd import pandas as pd +import pytest from shapely.geometry import LineString from workflow.scripts.build_network import build_network from workflow.scripts.clean import ( _apply_corrections, + _clean_cables, + _clean_circuits, + _clean_lines, + _clean_rating, + _clean_voltage, + _create_single_link, _filter_by_voltage, _region_ac_hz, + _region_min_voltage, clean, ) @@ -29,7 +38,18 @@ def _ring(lon0, lat0, lon1, lat1): ] -_NETWORK = {"minimum_voltage_kv": 220, "frequency_hz": {"AC": 50.0, "DC": 0.0}} +_NETWORK = { + "minimum_voltage_kv": 220, + "minimum_voltage_dc_kv": 150, + "frequency_hz": {"AC": 50.0, "DC": 0.0}, + "accepted_ac_frequencies_hz": [50.0, 60.0], + "frequency_tolerance_hz": 0.1, + "dc_lines": "keep", + "station_merge_radius_m": 500.0, +} +_FREQUENCY = {"accepted_ac_hz": [50.0, 60.0], "tolerance_hz": 0.1} +# build_network's config-driven settings that have no default. +_BUILD = {"station_bus_offset_m": 15.0, "overpassing_lines_tolerance_m": 1.0} _GEO_CRS = "EPSG:4326" _DISTANCE_CRS = "EPSG:3035" @@ -197,6 +217,7 @@ def test_builder_merges_compatible_segments_through_virtual_bus(monkeypatch): "line_id": "way/1-220", "circuits": 1, "voltage": 220000, + "country": "BE", "underground": False, "under_construction": False, "start_date": pd.NaT, @@ -207,6 +228,7 @@ def test_builder_merges_compatible_segments_through_virtual_bus(monkeypatch): "line_id": "way/2-220", "circuits": 1, "voltage": 220000, + "country": "BE", "underground": False, "under_construction": False, "start_date": pd.NaT, @@ -228,15 +250,18 @@ def capture_station_merge_radius(*args, **kwargs): build_network.__globals__, "_create_station_seeds", capture_station_merge_radius ) - buses, built_lines, transformers, stations_polygon, buses_polygon = build_network( - substations, - substations_polygon, - lines, - False, - None, - _GEO_CRS, - _DISTANCE_CRS, - station_merge_radius_m=1, + (buses, built_lines, transformers, _converters, stations_polygon, buses_polygon) = ( + build_network( + substations, + substations_polygon, + lines, + False, + None, + _GEO_CRS, + _DISTANCE_CRS, + station_merge_radius_m=1, + **_BUILD, + ) ) assert len(buses) == 2 @@ -247,3 +272,628 @@ def capture_station_merge_radius(*args, **kwargs): assert len(stations_polygon) == 2 assert list(buses_polygon.columns) == ["bus_id", "geometry"] assert buses_polygon.empty + + +def test_clean_keeps_country_on_lines(): + """Lines carry their country out of cleaning, as substations already did. + + A consumer resolves line types per country, so dropping it here would + leave no way to recover it downstream. + """ + import tempfile + from pathlib import Path + + with tempfile.TemporaryDirectory() as tmp: + tmp_path = Path(tmp) + lines_path = tmp_path / "BE_lines_way.json" + _write( + lines_path, + [ + { + "type": "way", + "id": 10, + "tags": {"power": "line", "voltage": "380000", "circuits": "1"}, + "geometry": [{"lon": 4.0, "lat": 50.0}, {"lon": 4.1, "lat": 50.0}], + } + ], + ) + _, _, lines = clean({"lines_way": [str(lines_path)]}, _NETWORK, {}, _GEO_CRS) + + assert "country" in lines.columns + assert lines.iloc[0]["country"] == "BE" + + +def test_builder_carries_provenance_columns_to_output(): + """country/under_construction/start_date reach the built components. + + Without them a consumer cannot resolve line types, and + ``remove_under_construction: false`` produces output in which a planned + asset is indistinguishable from a commissioned one. + """ + substations = gpd.GeoDataFrame( + { + "bus_id": pd.Series(dtype=str), + "voltage": pd.Series(dtype=int), + "country": pd.Series(dtype=str), + "under_construction": pd.Series(dtype=bool), + "start_date": pd.Series(dtype="datetime64[ns]"), + "contains": pd.Series(dtype=str), + }, + geometry=gpd.GeoSeries([], crs=_GEO_CRS), + crs=_GEO_CRS, + ) + substations_polygon = gpd.GeoDataFrame( + {"bus_id": pd.Series(dtype=str)}, + geometry=gpd.GeoSeries([], crs=_GEO_CRS), + crs=_GEO_CRS, + ) + lines = gpd.GeoDataFrame( + [ + { + "line_id": "way/1-380", + "circuits": 1, + "voltage": 380000, + "country": "BE;NL", + "underground": False, + "under_construction": True, + "start_date": pd.Timestamp("2030-01-01"), + "contains": ["way/1"], + "geometry": LineString([(4.0, 50.0), (4.3, 50.0)]), + } + ], + crs=_GEO_CRS, + ) + + buses, built_lines, _, _, _, _ = build_network( + substations, + substations_polygon, + lines, + remove_under_construction=False, + remove_after=None, + geo_crs=_GEO_CRS, + distance_crs=_DISTANCE_CRS, + station_merge_radius_m=1, + **_BUILD, + ) + + for frame in (buses, built_lines): + for column in ("country", "under_construction", "start_date"): + assert column in frame.columns + + assert built_lines.iloc[0]["under_construction"] + assert built_lines.iloc[0]["country"] == "BE;NL" + # Virtual buses are created from the line, so they inherit its country + # rather than arriving with a null a consumer would later drop. + assert set(buses["country"]) == {"BE;NL"} + assert buses["under_construction"].all() + + +def test_region_lookups_resolve_cross_border_country_values(): + """A ";"-joined country still resolves its regional overrides. + + ``_drop_duplicate_lines`` merges the codes of an element seen by two + countries, and looking that value up as a single key used to miss every + override and fall back to the global defaults. For a Brazilian + cross-border line that silently meant a 220 kV threshold instead of 60, + deleting a real interconnector, and 50 Hz instead of 60. + """ + regions = { + "BR": {"minimum_voltage_kv": 60.0, "frequency_hz": {"AC": 60.0, "DC": None}}, + "PY": {"minimum_voltage_kv": 60.0, "frequency_hz": {"AC": 60.0, "DC": None}}, + } + + assert _region_min_voltage("BR;PY", _NETWORK, regions) == 60_000 + assert _region_ac_hz("BR;PY", _NETWORK, regions) == "60" + + # Single-country values behave exactly as before. + assert _region_min_voltage("BR", _NETWORK, regions) == 60_000 + assert _region_min_voltage("BE", _NETWORK, regions) == 220_000 + + +def test_region_min_voltage_takes_the_most_permissive_side(): + """An interconnector survives if either country would keep it.""" + regions = {"BR": {"minimum_voltage_kv": 60.0, "frequency_hz": None}} + assert _region_min_voltage("BR;BE", _NETWORK, regions) == 60_000 + + +def test_clean_merges_custom_elements_with_retrieved_ones(tmp_path): + """A custom raw file is cleaned exactly like a retrieved one. + + Both retrieval backends already write the same shape, and the loader + takes a list of paths per feature while reading the country back out of + each filename, so custom data needs no special handling here. + """ + retrieved = tmp_path / "BE_lines_way.json" + _write( + retrieved, + [ + { + "type": "way", + "id": 10, + "tags": {"power": "line", "voltage": "380000", "circuits": "1"}, + "geometry": [{"lon": 4.0, "lat": 50.0}, {"lon": 4.1, "lat": 50.0}], + } + ], + ) + custom = tmp_path / "custom" / "BE_lines_way.json" + custom.parent.mkdir() + _write( + custom, + [ + { + "type": "way", + "id": 11, + "tags": {"power": "line", "voltage": "380000", "circuits": "1"}, + "geometry": [{"lon": 5.0, "lat": 51.0}, {"lon": 5.1, "lat": 51.0}], + } + ], + ) + + _, _, lines = clean( + {"lines_way": [str(retrieved), str(custom)]}, _NETWORK, {}, _GEO_CRS + ) + + assert set(lines["line_id"]) == {"way/10", "way/11"} + assert set(lines["country"]) == {"BE"} + + +def _circuits_for(cables=None, circuits=None, voltage="220000"): + """Run one raw line through cleaning and return its total circuits.""" + frame = pd.DataFrame( + [ + { + "id": "way/1", + "voltage": voltage, + "cables": cables, + "circuits": circuits, + "frequency": "50", + "_ac_hz": "50", + } + ] + ) + frame["voltage"] = _clean_voltage(frame["voltage"]) + frame["cables"] = _clean_cables(frame["cables"]) + frame["circuits"] = _clean_circuits(frame["circuits"]) + cleaned = _clean_lines(frame, ["220000"], "0", **_FREQUENCY) + return sum(int(value) for value in cleaned["circuits"]) + + +@pytest.mark.parametrize( + ("cables", "expected"), + [("3+3", 2), ("6+1", 2), ("2x3", 2), ("3x2", 2), ("2x2", 1), ("2-1", 1)], +) +def test_arithmetic_cable_tags_are_evaluated_not_concatenated(cables, expected): + """Stripping non-digits would read "2x3" as 23 cables instead of 6. + + These values are whole-value corrections precisely because the digit + strip at the end of _clean_cables cannot know they encode arithmetic. + """ + assert _circuits_for(cables=cables) == expected + + +@pytest.mark.parametrize(("circuits", "expected"), [("2/3", 2), ("2-1", 2)]) +def test_range_circuit_tags_are_not_concatenated(circuits, expected): + """Stripping non-digits would read "2/3" as 23 circuits instead of 2.""" + assert _circuits_for(circuits=circuits) == expected + + +@pytest.mark.parametrize("cables", ["ground", "1 disused"]) +def test_lines_without_live_conductors_are_dropped(cables): + """A ground wire or a retired cable is not a live single-circuit line.""" + assert _circuits_for(cables=cables) == 0 + + +@pytest.mark.parametrize( + ("cables", "expected"), [("triple", 1), ("single", 1), ("3;3 disused", 1)] +) +def test_word_and_partly_disused_cable_tags_keep_their_meaning(cables, expected): + """Word forms and partly disused cables resolve to the live circuit count.""" + assert _circuits_for(cables=cables) == expected + + +def test_unmapped_corrupting_tag_values_are_reported(caplog): + """The tables only grow if the workflow says what it could not map.""" + with caplog.at_level(logging.WARNING, logger="workflow.scripts.clean"): + result = list(_clean_cables(pd.Series(["4x5", "4x5", "7+2"]))) + assert result == ["45", "45", "72"], "the strip still runs" + assert "4x5" in caplog.text + assert "7+2" in caplog.text + + +def test_harmless_and_mapped_tag_values_are_not_reported(caplog): + """A trailing unit strips away cleanly, so warning about it is noise.""" + with caplog.at_level(logging.WARNING, logger="workflow.scripts.clean"): + assert list(_clean_voltage(pd.Series(["220000 V"]))) == ["220000"] + assert list(_clean_cables(pd.Series(["2x3"]))) == ["6"] + assert caplog.text == "" + + +# --- Frequency and DC ------------------------------------------------------- + + +def _raw_line(way_id, voltage, lon_lats, frequency=None, power="line"): + tags = {"power": power, "voltage": voltage} + if frequency is not None: + tags["frequency"] = frequency + return { + "type": "way", + "id": way_id, + "tags": tags, + "geometry": [{"lon": lon, "lat": lat} for lon, lat in lon_lats], + } + + +def _converter_station_inputs(tmp_path): + """A converter station fed by 380 kV and 220 kV AC lines and a 320 kV DC line.""" + substations_path = tmp_path / "BE_substations_way.json" + _write( + substations_path, + [ + { + "type": "way", + "id": 1, + "tags": { + "power": "substation", + "voltage": "380000;220000;320000", + "frequency": "50;50;0", + }, + "geometry": _ring(3.999, 49.999, 4.001, 50.001), + } + ], + ) + lines_path = tmp_path / "BE_lines_way.json" + _write( + lines_path, + [ + _raw_line(10, "380000", [(3.9, 50.0), (4.0, 50.0)]), + _raw_line(11, "220000", [(4.0, 49.9), (4.0, 50.0)]), + _raw_line(12, "320000", [(4.1, 50.0), (4.0, 50.0)], frequency="0"), + ], + ) + return {"substations_way": [str(substations_path)], "lines_way": [str(lines_path)]} + + +def _build(buses, polygons, lines): + return build_network( + buses, polygons, lines, False, None, _GEO_CRS, _DISTANCE_CRS, **_BUILD + ) + + +def test_converter_station_gets_dc_bus_and_converter_not_transformer(tmp_path): + """PyPSA-Earth's model: DC has its own bus, joined to AC only by a converter.""" + buses, polygons, lines = clean( + _converter_station_inputs(tmp_path), _NETWORK, {}, _GEO_CRS + ) + assert lines.set_index("voltage")["dc"].to_dict() == { + 380000: False, + 220000: False, + 320000: True, + } + + built_buses, built_lines, transformers, converters, _, _ = _build( + buses, polygons, lines + ) + station = built_buses[built_buses["station_id"] == "way/1"] + assert set(station["bus_id"]) == {"way/1-380", "way/1-220", "way/1-320-dc"} + assert station.set_index("bus_id")["dc"].to_dict()["way/1-320-dc"] + + dc_line = built_lines[built_lines["dc"]] + assert len(dc_line) == 1 + assert "way/1-320-dc" in set(dc_line[["bus0", "bus1"]].iloc[0]) + + # Transformers join AC levels only. + assert list(transformers["transformer_id"]) == ["way/1-380-220"] + + # The DC side converts to the AC level closest in voltage: 380 kV, not 220. + assert len(converters) == 1 + converter = converters.iloc[0] + assert (converter["bus0"], converter["bus1"]) == ("way/1-320-dc", "way/1-380") + assert (converter["voltage_bus0_kv"], converter["voltage_bus1_kv"]) == (320, 380) + + +@pytest.mark.parametrize("mode", ["drop", "force_ac"]) +def test_dc_lines_can_be_dropped_or_relabelled_as_ac(tmp_path, mode): + """network.dc_lines drop removes DC lines; force_ac keeps them as AC.""" + network = {**_NETWORK, "dc_lines": mode} + buses, polygons, lines = clean( + _converter_station_inputs(tmp_path), network, {}, _GEO_CRS + ) + assert not lines["dc"].any() + assert not buses["dc"].any() + assert (320000 in set(lines["voltage"])) == (mode == "force_ac") + + _, _, _, converters, _, _ = _build(buses, polygons, lines) + assert converters.empty + + +def test_railway_traction_frequency_is_dropped_not_read_as_mains(tmp_path): + """A 16.7 Hz traction line is a separate grid, not a 50 Hz line.""" + lines_path = tmp_path / "DE_lines_way.json" + _write( + lines_path, + [ + _raw_line(1, "220000", [(10.0, 50.0), (10.1, 50.0)], frequency="16.7"), + _raw_line(2, "220000", [(10.0, 51.0), (10.1, 51.0)], frequency="50"), + ], + ) + _, _, lines = clean({"lines_way": [str(lines_path)]}, _NETWORK, {}, _GEO_CRS) + assert list(lines["line_id"]) == ["way/2"] + + +def test_frequency_list_is_paired_with_voltage_list(tmp_path): + """voltage=380000;220000 frequency=50;16.7 is a mains circuit plus a railway one.""" + lines_path = tmp_path / "DE_lines_way.json" + _write( + lines_path, + [ + _raw_line( + 1, "380000;220000", [(10.0, 50.0), (10.1, 50.0)], frequency="50;16.7" + ) + ], + ) + _, _, lines = clean({"lines_way": [str(lines_path)]}, _NETWORK, {}, _GEO_CRS) + assert list(lines["voltage"]) == [380000] + + +@pytest.mark.parametrize( + ("frequency", "expected_dc"), [("50.0", False), ("60", False), ("0.0", True)] +) +def test_frequency_matching_is_numeric_not_textual(tmp_path, frequency, expected_dc): + """A tag of 50.0 is mains and 0.0 is DC, though neither equals its marker.""" + lines_path = tmp_path / "BE_lines_way.json" + _write( + lines_path, + [_raw_line(1, "220000", [(4.0, 50.0), (4.1, 50.0)], frequency=frequency)], + ) + _, _, lines = clean({"lines_way": [str(lines_path)]}, _NETWORK, {}, _GEO_CRS) + assert list(lines["dc"]) == [expected_dc] + + +def test_empty_substation_input_keeps_the_output_schema(tmp_path): + """build_network selects polygons by bus_id, so the columns must exist.""" + lines_path = tmp_path / "BE_lines_way.json" + _write(lines_path, [_raw_line(1, "220000", [(4.0, 50.0), (4.1, 50.0)])]) + buses, polygons, lines = clean( + {"lines_way": [str(lines_path)]}, _NETWORK, {}, _GEO_CRS + ) + assert buses.empty + assert {"bus_id", "dc"} <= set(buses.columns) + assert "bus_id" in polygons.columns + built_buses, built_lines, *_ = _build(buses, polygons, lines) + assert len(built_lines) == 1 + assert not built_lines["dc"].any() + + +def test_dc_cables_count_two_per_circuit_on_split_lines(tmp_path): + """A DC circuit has two conductors, so 4;4 cables mean two circuits each, not one.""" + lines_path = tmp_path / "BE_lines_way.json" + _write( + lines_path, + [ + { + **_raw_line( + 1, "320000;320000", [(4.0, 50.0), (4.1, 50.0)], frequency="0" + ), + "tags": { + "power": "line", + "voltage": "320000;320000", + "frequency": "0", + "cables": "4;4", + }, + } + ], + ) + _, _, lines = clean({"lines_way": [str(lines_path)]}, _NETWORK, {}, _GEO_CRS) + assert list(lines["circuits"]) == [4] + assert list(lines["dc"]) == [True] + + +# --- PyPSA-Eur's relation-based HVDC links --------------------------------- + + +def _member(way_id, lon_lats, role): + return { + "type": "way", + "ref": way_id, + "role": role, + "geometry": [{"lon": lon, "lat": lat} for lon, lat in lon_lats], + } + + +def _hvdc_relation(relation_id, members, voltage="500000", rating="1000 MW"): + tags = {"type": "route", "route": "power", "voltage": voltage, "frequency": "0"} + if rating is not None: + tags["rating"] = rating + return {"type": "relation", "id": relation_id, "tags": tags, "members": members} + + +@pytest.mark.parametrize( + ("raw", "expected"), + [("1000 MW", 1000.0), ("1.2 GW", 1200.0), ("500;500", 1000.0), ("", None)], +) +def test_rating_parses_to_megawatts(raw, expected): + """HVDC ratings sum across entries and scale GW to MW.""" + [value] = _clean_rating(pd.Series([raw])).tolist() + assert (pd.isna(value) and expected is None) or value == expected + + +def test_single_link_collapses_parallel_poles_and_skips_terminals(): + """A bipole's two poles would never merge into one line, so AC merging drops it.""" + pole = [(20.0, -5.0), (21.0, -5.0)] + row = pd.Series( + { + "members": [ + _member(1, pole, "line"), + _member(2, [(20.0, -5.0), (20.5, -5.0001), (21.0, -5.0)], "line"), + _member( + 3, + [(19.99, -5.01), (20.01, -5.01), (20.01, -4.99), (19.99, -5.01)], + "substation", + ), + ] + } + ) + link, members = _create_single_link(row) + assert link.geom_type == "LineString" + assert set(members) == {"way/1", "way/2"}, "both poles are replaced by the link" + + +def _write_hvdc_case(tmp_path, converter_tag=True, gap_deg=0.03): + """A 500 kV bipole into a converter hall standing apart from its 220 kV substation.""" + hall_tags = {"power": "substation", "voltage": "500000", "frequency": "0"} + if converter_tag: + hall_tags["substation"] = "converter" + ac_lon = 21.0 + gap_deg + _write( + tmp_path / "CD_substations_way.json", + [ + { + "type": "way", + "id": 1, + "tags": hall_tags, + "geometry": _ring(20.999, -5.001, 21.001, -4.999), + }, + { + "type": "way", + "id": 2, + "tags": {"power": "substation", "voltage": "220000"}, + "geometry": _ring(ac_lon - 0.001, -5.001, ac_lon + 0.001, -4.999), + }, + ], + ) + poles = [ + _member(10, [(19.0, -5.0), (21.0, -5.0)], "line"), + _member(11, [(19.0, -5.0), (20.0, -5.0002), (21.0, -5.0)], "line"), + ] + _write(tmp_path / "CD_routes_relation.json", [_hvdc_relation(100, poles)]) + _write( + tmp_path / "CD_lines_way.json", + [ + _raw_line(10, "500000", [(19.0, -5.0), (21.0, -5.0)], frequency="0"), + _raw_line( + 11, + "500000", + [(19.0, -5.0), (20.0, -5.0002), (21.0, -5.0)], + frequency="0", + ), + _raw_line(20, "220000", [(ac_lon, -5.0), (ac_lon, -4.0)]), + ], + ) + return { + "substations_way": [str(tmp_path / "CD_substations_way.json")], + "routes_relation": [str(tmp_path / "CD_routes_relation.json")], + "lines_way": [str(tmp_path / "CD_lines_way.json")], + } + + +def test_hvdc_relation_becomes_one_rated_link_replacing_its_poles(tmp_path): + """PyPSA-Eur's relationship concept: one rated link per HVDC relation.""" + _, _, lines = clean(_write_hvdc_case(tmp_path), _NETWORK, {}, _GEO_CRS) + dc = lines[lines["dc"]] + assert list(dc["line_id"]) == ["relation/100"] + assert list(dc["p_nom_mw"]) == [1000.0] + assert lines["p_nom_mw"].isna().sum() == len(lines) - 1, "AC lines carry no rating" + + +def test_converter_hall_apart_from_its_substation_is_wired_to_nearest_ac(tmp_path): + """PyPSA-Eur's rule: a tagged converter hall 3 km from the AC yard still connects.""" + buses, polygons, lines = clean(_write_hvdc_case(tmp_path), _NETWORK, {}, _GEO_CRS) + *_, converters, _, _ = build_network( + buses, + polygons, + lines, + False, + None, + _GEO_CRS, + _DISTANCE_CRS, + converter_search_radius_m=50000, + **_BUILD, + ) + assert len(converters) == 1 + converter = converters.iloc[0] + assert converter["pairing"] == "nearest_station" + assert converter["bus1"] == "way/2-220" + assert converter["p_nom_mw"] == 1000.0 + + +def test_untagged_dc_terminal_is_not_wired_to_nearby_ac(tmp_path): + """Without the converter tag a DC end may be a border cut, so no guessing.""" + buses, polygons, lines = clean( + _write_hvdc_case(tmp_path, converter_tag=False), _NETWORK, {}, _GEO_CRS + ) + *_, converters, _, _ = build_network( + buses, + polygons, + lines, + False, + None, + _GEO_CRS, + _DISTANCE_CRS, + converter_search_radius_m=50000, + **_BUILD, + ) + assert converters.empty + + +def test_dc_voltage_floor_is_separate_from_ac(tmp_path): + """150 kV HVDC survives the 220 kV AC floor; a 150 kV AC line does not.""" + _write( + tmp_path / "EE_lines_way.json", + [ + _raw_line(1, "150000", [(24.0, 59.0), (24.1, 59.0)], frequency="0"), + _raw_line(2, "150000", [(24.0, 58.0), (24.1, 58.0)], frequency="50"), + ], + ) + _, _, lines = clean( + {"lines_way": [str(tmp_path / "EE_lines_way.json")]}, _NETWORK, {}, _GEO_CRS + ) + assert list(lines["line_id"]) == ["way/1"] + + +def test_lower_regional_ac_floor_takes_effect(tmp_path): + """A region's 60 kV floor used to be cut off by the global 220 kV filter first.""" + _write( + tmp_path / "BR_lines_way.json", + [_raw_line(1, "69000", [(-47.0, -15.0), (-47.1, -15.0)], frequency="60")], + ) + regions = { + "BR": {"minimum_voltage_kv": 60.0, "frequency_hz": {"AC": 60.0, "DC": None}} + } + _, _, lines = clean( + {"lines_way": [str(tmp_path / "BR_lines_way.json")]}, + _NETWORK, + regions, + _GEO_CRS, + ) + assert list(lines["voltage"]) == [69000] + + +def test_relation_of_cable_ways_is_underground(tmp_path): + """A relation whose member ways are all cables is underground.""" + _write( + tmp_path / "BE_routes_relation.json", + [ + _hvdc_relation( + 100, [_member(1, [(4.0, 50.0), (4.1, 50.0)], "cable")], voltage="320000" + ) + ], + ) + _write( + tmp_path / "BE_cables_way.json", + [ + _raw_line( + 1, "320000", [(4.0, 50.0), (4.1, 50.0)], frequency="0", power="cable" + ) + ], + ) + _, _, lines = clean( + { + "routes_relation": [str(tmp_path / "BE_routes_relation.json")], + "cables_way": [str(tmp_path / "BE_cables_way.json")], + }, + _NETWORK, + {}, + _GEO_CRS, + ) + assert list(lines["line_id"]) == ["relation/100"] + assert list(lines["underground"]) == [True] From 9745dfa35baa958835597937d35b9d7d7410d0ae Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 27 Sep 2026 11:01:14 +0200 Subject: [PATCH 54/54] Externalize custom files input functionality --- workflow/rules/network.smk | 28 +++++++--------------------- workflow/scripts/_helpers.py | 10 ++++++++++ 2 files changed, 17 insertions(+), 21 deletions(-) diff --git a/workflow/rules/network.smk b/workflow/rules/network.smk index 7d37160..0386707 100644 --- a/workflow/rules/network.smk +++ b/workflow/rules/network.smk @@ -2,21 +2,7 @@ # # SPDX-License-Identifier: MIT -from pathlib import Path - - -def _custom_files(feature): - """Custom raw files for one feature, cleaned alongside the retrieved ones. - - Selected by filename, the same ``{country}_{feature}.json`` convention - retrieval writes and ``clean`` reads the country back out of, so a - custom file needs no special handling downstream. - """ - return [ - path - for path in config["custom_data"]["files"] - if Path(path).stem.endswith(f"_{feature}") - ] +from scripts._helpers import custom_files rule clean: @@ -25,33 +11,33 @@ rule clean: "/retrieve/{country}_lines_way.json", country=config["countries"], ) - + _custom_files("lines_way"), + + custom_files(config["custom_data"]["files"], "lines_way"), cables_way=expand( "/retrieve/{country}_cables_way.json", country=config["countries"], ) - + _custom_files("cables_way"), + + custom_files(config["custom_data"]["files"], "cables_way"), substations_way=expand( "/retrieve/{country}_substations_way.json", country=config["countries"], ) - + _custom_files("substations_way"), + + custom_files(config["custom_data"]["files"], "substations_way"), substations_node=expand( "/retrieve/{country}_substations_node.json", country=config["countries"], ) - + _custom_files("substations_node"), + + custom_files(config["custom_data"]["files"], "substations_node"), substations_relation=expand( "/retrieve/{country}_substations_relation.json", country=config["countries"], ) - + _custom_files("substations_relation"), + + custom_files(config["custom_data"]["files"], "substations_relation"), routes_relation=( expand( "/retrieve/{country}_routes_relation.json", country=config["countries"], ) - + _custom_files("routes_relation") + + custom_files(config["custom_data"]["files"], "routes_relation") if config["network"]["include_relations"] else [] ), diff --git a/workflow/scripts/_helpers.py b/workflow/scripts/_helpers.py index c438fa5..b30ffdd 100644 --- a/workflow/scripts/_helpers.py +++ b/workflow/scripts/_helpers.py @@ -63,6 +63,16 @@ def load_config_defaults() -> Any: return yaml.safe_load(handle) +def custom_files(files: list[str], feature: str) -> list[str]: + """Custom raw files for one feature, cleaned alongside the retrieved ones. + + Selected by filename, the same ``{country}_{feature}.json`` convention + retrieval writes and ``clean`` reads the country back out of, so a + custom file needs no special handling downstream. + """ + return [path for path in files if Path(path).stem.endswith(f"_{feature}")] + + # Defaults for functions called outside Snakemake, whose rules pass their own # resolved config values instead. GEO_CRS: str = load_config_defaults()["crs"]["geo"]