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: "" + 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 821a651..175bffe 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,9 +12,13 @@ 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**](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. 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). @@ -24,8 +28,8 @@ 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. +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 @@ -33,14 +37,12 @@ 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 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). @@ -52,10 +54,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. @@ -106,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 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] diff --git a/config/README.md b/config/README.md index e5166c9..07fd48f 100644 --- a/config/README.md +++ b/config/README.md @@ -1,49 +1,45 @@ -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`). Retrieven data are transferred to the cleaning phase and after that are used to build a topologically-clean network model. + +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.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 + +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: + files: + - data/custom/BE_lines_way.json +``` + +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 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 ... ``` -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. diff --git a/config/config.schema.json b/config/config.schema.json index 6b20258..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.", @@ -105,6 +111,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.", @@ -123,12 +135,59 @@ } } }, + "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" + }, + "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", + "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", "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", @@ -209,6 +268,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..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: "" @@ -20,15 +21,27 @@ retrieve: network: include_relations: true minimum_voltage_kv: 220.0 + minimum_voltage_dc_kv: 150.0 frequency_hz: AC: 50.0 DC: 0.0 + accepted_ac_frequencies_hz: + - 50.0 + - 60.0 + frequency_tolerance_hz: 0.1 + 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 regions: {} +custom_data: + files: [] + interactive_map: coordinate_decimals: 5 simplify_geometries: 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" 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] diff --git a/workflow/internal/tag_corrections.yaml b/workflow/internal/tag_corrections.yaml index ec83725..3f9aed0 100644 --- a/workflow/internal/tag_corrections.yaml +++ b/workflow/internal/tag_corrections.yaml @@ -1,29 +1,33 @@ # 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 +# 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. # -# 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. +# 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 @@ -40,6 +44,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: [":", ";"] @@ -55,11 +64,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 @@ -83,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 diff --git a/workflow/rules/network.smk b/workflow/rules/network.smk index 3d5476b..0386707 100644 --- a/workflow/rules/network.smk +++ b/workflow/rules/network.smk @@ -2,34 +2,42 @@ # # SPDX-License-Identifier: MIT +from scripts._helpers import custom_files + rule clean: input: lines_way=expand( "/retrieve/{country}_lines_way.json", country=config["countries"], - ), + ) + + custom_files(config["custom_data"]["files"], "lines_way"), cables_way=expand( "/retrieve/{country}_cables_way.json", country=config["countries"], - ), + ) + + custom_files(config["custom_data"]["files"], "cables_way"), substations_way=expand( "/retrieve/{country}_substations_way.json", country=config["countries"], - ), + ) + + custom_files(config["custom_data"]["files"], "substations_way"), substations_node=expand( "/retrieve/{country}_substations_node.json", country=config["countries"], - ), + ) + + custom_files(config["custom_data"]["files"], "substations_node"), substations_relation=expand( "/retrieve/{country}_substations_relation.json", country=config["countries"], - ), + ) + + custom_files(config["custom_data"]["files"], "substations_relation"), routes_relation=( expand( "/retrieve/{country}_routes_relation.json", country=config["countries"], ) + + custom_files(config["custom_data"]["files"], "routes_relation") if config["network"]["include_relations"] else [] ), @@ -64,9 +72,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: @@ -76,6 +86,9 @@ 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"], crs=config["crs"].model_dump(mode="json"), 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: 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/_helpers.py b/workflow/scripts/_helpers.py index 0ca446b..b30ffdd 100644 --- a/workflow/scripts/_helpers.py +++ b/workflow/scripts/_helpers.py @@ -26,9 +26,6 @@ logger = logging.getLogger(__name__) -GEO_CRS = "EPSG:4326" -BUS_TOL = 500 # metres; default station merge tolerance - def configure_logging(log_path: str) -> None: """Send rule and dependency logging to the Snakemake log file.""" @@ -59,6 +56,29 @@ 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) + + +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"] +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 4877b15..e931120 100644 --- a/workflow/scripts/_schema.py +++ b/workflow/scripts/_schema.py @@ -12,12 +12,29 @@ 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 _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. @@ -92,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) @@ -151,15 +176,81 @@ 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", ) + 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, + ) + 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", 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" ) @@ -169,6 +260,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.""" @@ -242,6 +370,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 " @@ -281,6 +410,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", diff --git a/workflow/scripts/build_interactive_map.py b/workflow/scripts/build_interactive_map.py index 4fb5e97..29f90c3 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, @@ -107,17 +106,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"), @@ -135,7 +137,6 @@ def polygon_layer( frame: gpd.GeoDataFrame, name: str, color: list[int], - *, geo_crs: str, distance_crs: str, coord_decimals: int, @@ -150,16 +151,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: @@ -188,7 +188,7 @@ def build_map( transformers: gpd.GeoDataFrame, stations: gpd.GeoDataFrame, bus_polygons: gpd.GeoDataFrame, - *, + converters: gpd.GeoDataFrame | None, geo_crs: str, distance_crs: str, stations_simplify_m: float | None, @@ -196,7 +196,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 an interactive network map.""" lines = lines.copy() if not lines.empty: lines["color"] = line_colors(lines["voltage_kv"]) @@ -240,6 +240,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 ] @@ -1086,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 e4803af..ba0f3ff 100644 --- a/workflow/scripts/build_network.py +++ b/workflow/scripts/build_network.py @@ -36,6 +36,65 @@ 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", + "dc", + "country", + "under_construction", + "start_date", + "osm_ids", + "geometry", +] + +LINE_COLUMNS = [ + "line_id", + "bus0", + "bus1", + "voltage_kv", + "dc", + "p_nom_mw", + "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", +] + +# 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"] + def _empty_geodataframe(columns: list[str], crs: str) -> gpd.GeoDataFrame: """Build an empty GeoDataFrame with the given non-geometry ``columns``.""" @@ -44,6 +103,24 @@ 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.""" + return [column for column in columns if column != "geometry"] + + +def _merge_country_codes(values: Any) -> str: + """Clean-up country codes for multy-country entries. + + This 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: @@ -66,11 +143,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] @@ -117,8 +194,21 @@ 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"]] + """Create deterministic virtual buses at each unique (voltage, endpoint) pair. + + A virtual bus inherits its attributes from the lines that meet there. + """ + line_data = lines[ + [ + "voltage", + "dc", + "geometry", + "line_id", + "country", + "under_construction", + "start_date", + ] + ] line_geoms = line_data["geometry"].apply(_remove_loops_from_multiline) endpoints0 = line_data.assign( @@ -144,10 +234,20 @@ 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"]), + "under_construction": bool(group["under_construction"].all()), + "start_date": group["start_date"].min(), + } + ) + # 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() ) @@ -157,7 +257,18 @@ 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", + "dc", + "geometry", + "contains", + "country", + "under_construction", + "start_date", + ] + ] def _split_linestring_by_point( @@ -182,14 +293,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 @@ -204,6 +315,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] @@ -266,15 +380,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 = ( @@ -304,7 +419,19 @@ 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", + "dc", + "p_nom_mw", + "circuits", + "length", + "geometry", + "underground", + "country", + "under_construction", + "start_date", + ], ] lines_to_merge_dict = [ (node, row.to_dict()) @@ -325,6 +452,13 @@ 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. + 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()] ) @@ -340,13 +474,29 @@ 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() + ) + 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}", "circuits": circuits, "voltage": voltage, + "dc": dc, + "p_nom_mw": p_nom_mw, "geometry": geometry, "underground": underground, + "country": country, + "under_construction": under_construction, + "start_date": start_date, "contains_lines": contains_lines, "contains_buses": contains_buses, } @@ -356,8 +506,13 @@ def _create_merge_mapping( "line_id", "circuits", "voltage", + "dc", + "p_nom_mw", "geometry", "underground", + "country", + "under_construction", + "start_date", "contains_lines", "contains_buses", ] @@ -388,7 +543,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] @@ -475,7 +629,11 @@ 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, + 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) @@ -483,39 +641,49 @@ 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): - poi_x_offset = poi_x + offset * np.sin( - np.pi / 4 + 2 * np.pi * idx / len(voltages) + for idx, (voltage, dc) in enumerate(levels): + 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( - np.pi / 4 + 2 * np.pi * idx / len(voltages) + 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)) - 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"] @@ -559,14 +727,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"] @@ -667,8 +840,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 @@ -708,6 +886,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, @@ -716,25 +992,47 @@ def build_network( remove_after: str | None, geo_crs: str, distance_crs: str, + 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, gpd.GeoDataFrame, 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 ``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"]) + # 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( - buses, remove_under_construction, remove_after - ).drop(columns=["start_date"]) + substations, remove_under_construction, remove_after + ) buses_polygon = substations_polygon[ substations_polygon["bus_id"].isin(buses["bus_id"]) @@ -745,25 +1043,29 @@ 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_converters = _empty_geodataframe( + _non_geometry(CONVERTER_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, empty_lines, empty_transformers, + empty_converters, empty_stations_polygon, 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( @@ -776,7 +1078,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] @@ -798,8 +1102,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, + offset_m=station_bus_offset_m, ) buses["geometry"] = gpd.points_from_xy( @@ -835,6 +1154,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 @@ -851,11 +1178,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() @@ -863,21 +1187,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() @@ -889,40 +1199,35 @@ 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 + ) + + 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_id", "geometry"]].copy() + stations_polygon_out = stations[STATION_POLYGON_COLUMNS].copy() buses_polygon_out = buses_polygon.copy() return ( buses_out, lines_out, transformers_out, + converters_out, stations_polygon_out, buses_polygon_out, ) @@ -945,21 +1250,28 @@ def _write_components( snakemake = mock_snakemake("build_network") configure_logging(snakemake.log[0]) - buses, lines, transformers, 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"], - snakemake.params.station_merge_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, and %d transformers.", + "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) @@ -968,5 +1280,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") diff --git a/workflow/scripts/clean.py b/workflow/scripts/clean.py index 32f5fcd..3950a23 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 @@ -41,6 +42,45 @@ "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() + +# 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.""" @@ -83,28 +123,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: @@ -116,7 +184,117 @@ 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 +) -> 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 _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: @@ -128,6 +306,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: @@ -155,12 +353,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) @@ -192,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: @@ -224,7 +422,11 @@ 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. @@ -257,17 +459,17 @@ 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, + dc_hz: str, + accepted_ac_hz: list[float], + tolerance_hz: float, ) -> pd.DataFrame: """Clean lines/cables heuristically, deriving circuits from whatever tags exist. @@ -285,9 +487,22 @@ 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 + # 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 ) @@ -299,6 +514,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) @@ -384,7 +611,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"], ), ) ), @@ -404,6 +632,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 @@ -426,20 +663,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 @@ -474,54 +710,34 @@ 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", } ) .reset_index() ) - return df_lines[ - [ - "line_id", - "circuits", - "voltage", - "underground", - "under_construction", - "start_date", - "geometry", - "contains", - ] - ] + return df_lines[LINE_COLUMNS] def _finalise_lines(df_lines: pd.DataFrame) -> pd.DataFrame: @@ -529,18 +745,9 @@ 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", - "underground", - "under_construction", - "start_date", - "geometry", - "contains", - ] - ] + # 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) return df_lines @@ -716,6 +923,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", @@ -742,6 +989,7 @@ def _create_line(row: pd.Series) -> tuple[Any, list[str]]: "cables", "frequency", "voltage", + "rating", "construction", "construction:power", "start_date", @@ -923,13 +1171,65 @@ 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: - """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.""" + 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: @@ -945,11 +1245,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( @@ -967,8 +1288,15 @@ def clean( ``include_relations`` is off) is simply skipped. """ crs = geo_crs - min_voltage_ac = network["minimum_voltage_kv"] * 1000 # V + # 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.") @@ -978,8 +1306,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"] = ( @@ -989,22 +1322,28 @@ 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=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( 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- @@ -1018,7 +1357,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"]) @@ -1057,18 +1398,32 @@ 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"]) - 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", @@ -1078,13 +1433,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]) @@ -1106,21 +1461,10 @@ 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 - lines_frames.append( - df_relation[ - [ - "line_id", - "circuits", - "voltage", - "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.") @@ -1129,6 +1473,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( @@ -1144,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=min_voltage_ac - ) + 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"]) @@ -1156,21 +1507,9 @@ 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), - ) - - 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 = _clean_lines(df_lines, list_voltages, dc_hz, **frequency_options) + 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) @@ -1183,10 +1522,14 @@ 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 = 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"]) 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, )