Fix raw OSM build and release process

doc/release_notes.rst

@@ -9,19 +9,23 @@ Release Notes
 Upcoming Release
 ================
 
+* Fix virtual bus naming when building the transmission network from raw OSM data to use invariant names.
+
+* Fix column selection when preparing OSM pre-built releases.
+
 * Fix compatibility of rules `build_gas_input_locations` and `build_gas_network` with pyogrio >=0.12.0 (https://github.com/PyPSA/pypsa-eur/pull/1955).
 
 * Added interactive (html) balance maps `results/maps/interactive/` (https://github.com/PyPSA/pypsa-eur/pull/1935) based on https://docs.pypsa.org/latest/user-guide/plotting/explore/. Settings for interactive maps can be found in `plotting.default.yaml` under `plotting["balance_map_interactive"]`.
 
-* Relocated and modified static (pdf) balance maps to `results/maps/static/` (https://github.com/PyPSA/pypsa-eur/pull/1935) for better organization. 
+* Relocated and modified static (pdf) balance maps to `results/maps/static/` (https://github.com/PyPSA/pypsa-eur/pull/1935) for better organization.
 
-* With https://github.com/PyPSA/pypsa-eur/pull/1935, note that bus carriers for balance maps containing spaces need to be specified with underscores `_` in the configuration file, e.g., `co2_stored` instead of `co2 stored`. This is to ensure compatibility with queue managers like slurm. 
+* With https://github.com/PyPSA/pypsa-eur/pull/1935, note that bus carriers for balance maps containing spaces need to be specified with underscores `_` in the configuration file, e.g., `co2_stored` instead of `co2 stored`. This is to ensure compatibility with queue managers like slurm.
 
 * Fix building osm network using overpass API (https://github.com/PyPSA/pypsa-eur/pull/1940).
 
 * Added configuration option to set overpass API URL, maximum retries, timeout and user agent information (https://github.com/PyPSA/pypsa-eur/pull/1940 and https://pypsa-eur.readthedocs.io/en/latest/configuration.html#overpass_api). For a list of public overpass APIs see `here <https://wiki.openstreetmap.org/wiki/Overpass_API#Public_Overpass_API_instances>`_.
 
-* Refactored `solve_network.py` and `solve_operations_network.py` to separate optimization problem preparation from solving, enabling inspection of optimization problems before solve execution. 
+* Refactored `solve_network.py` and `solve_operations_network.py` to separate optimization problem preparation from solving, enabling inspection of optimization problems before solve execution.
 
 * Added example configurations for rolling horizon and iterative optimization modes in `config/examples/`.
 

scripts/build_osm_network.py

@@ -12,7 +12,7 @@
 import pandas as pd
 import pypsa
 from pyproj import Transformer
-from shapely import prepare
+from shapely import get_point, prepare
 from shapely.algorithms.polylabel import polylabel
 from shapely.geometry import LineString, MultiLineString, Point
 from shapely.ops import linemerge, split
@@ -148,18 +148,23 @@ def _add_line_endings(buses, lines, add=0, name="line-end"):
     -------
         - pd.DataFrame: DataFrame containing the virtual bus endpoints with columns 'bus_id', 'voltage', 'geometry', and 'contains'.
     """
-    endpoints0 = lines[["voltage", "geometry"]].copy()
-    endpoints0["geometry"] = endpoints0["geometry"].apply(lambda x: x.boundary.geoms[0])
-
-    endpoints1 = lines[["voltage", "geometry"]].copy()
-    endpoints1["geometry"] = endpoints1["geometry"].apply(lambda x: x.boundary.geoms[1])
-
+    line_data = lines[["voltage", "geometry", "line_id"]]
+    line_geoms = line_data["geometry"].apply(_remove_loops_from_multiline)
+    endpoints0 = line_data.assign(
+        geometry=get_point(line_geoms.geometry, 0), endpoint=0
+    )
+    endpoints1 = line_data.assign(
+        geometry=get_point(line_geoms.geometry, -1), endpoint=1
+    )
     endpoints = pd.concat([endpoints0, endpoints1], ignore_index=True)
     endpoints.drop_duplicates(subset=["geometry", "voltage"], inplace=True)
     endpoints.reset_index(drop=True, inplace=True)
 
-    endpoints["bus_id"] = endpoints.index + add + 1
-    endpoints["bus_id"] = "virtual" + "-" + endpoints["bus_id"].astype(str)
+    # Create deterministic ID from line_id and endpoint
+    endpoints["bus_id"] = (
+        "virtual_" + endpoints["line_id"] + "_" + endpoints["endpoint"].astype(str)
+    )
+    endpoints.drop(columns=["line_id", "endpoint"], inplace=True)
 
     endpoints["contains"] = name
 
@@ -651,10 +656,9 @@ def _create_station_seeds(
     buses_to_rename = buses_to_rename.sort_values(
         by=["country", "lat", "lon"], ascending=[True, False, True]
     )
-    buses_to_rename["bus_id"] = buses_to_rename.groupby("country").cumcount() + 1
-    buses_to_rename["bus_id"] = buses_to_rename["country"] + buses_to_rename[
-        "bus_id"
-    ].astype(str)
+    buses_to_rename["bus_id"] = buses_to_rename[
+        "country"
+    ] + buses_to_rename.index.str.replace("virtual", "")
 
     # Dict to rename virtual buses
     dict_rename = buses_to_rename["bus_id"].to_dict()
@@ -729,7 +733,7 @@ def _merge_buses_to_stations(
         voltages = sorted(
             g_value["voltage"].unique(), reverse=True
         )  # Sort voltags in descending order
-
+        not_virtual = ~g_value.bus_id.str.startswith("virtual_")
         if len(voltages) > 1:
             poi_x, poi_y = geo_to_dist.transform(
                 g_value["poi"].values[0].x, g_value["poi"].values[0].y
@@ -745,10 +749,11 @@ def _merge_buses_to_stations(
 
                 poi_offset = Point(dist_to_geo.transform(poi_x_offset, poi_y_offset))
 
-                # Update bus_name
-                g_value.loc[g_value["voltage"] == v, "bus_id"] = (
-                    g_name + "-" + str(int(v / 1000))
-                )
+                # Update bus_name if not virtual (in which case the voltage suffix is already present)
+                g_value.loc[
+                    (g_value["voltage"] == v) & not_virtual,
+                    "bus_id",
+                ] = g_name + "-" + str(int(v / 1000))
 
                 # Update geometry
                 g_value.loc[g_value["voltage"] == v, "geometry"] = poi_offset
@@ -757,7 +762,9 @@ def _merge_buses_to_stations(
             buses_all.loc[g_value.index, "geometry"] = g_value["geometry"]
         else:
             v = voltages[0]
-            buses_all.loc[g_value.index, "bus_id"] = g_name + "-" + str(int(v / 1000))
+            buses_all.loc[g_value.loc[not_virtual].index, "bus_id"] = (
+                g_name + "-" + str(int(v / 1000))
+            )
             buses_all.loc[g_value.index, "geometry"] = g_value["poi"]
 
     return buses_all
@@ -889,8 +896,9 @@ def _map_endpoints_to_buses(
     for coord in range(2):
         # Obtain endpoints
         endpoints = lines_all[["voltage", "geometry"]].copy()
-        endpoints["geometry"] = endpoints["geometry"].apply(
-            lambda x: x.boundary.geoms[coord]
+        # -1 * coord returns 0 for coord=0 and -1 for coord=1
+        endpoints["geometry"] = get_point(
+            endpoints.geometry.apply(_remove_loops_from_multiline), -1 * coord
         )
         if sjoin == "intersects":
             endpoints = gpd.sjoin(

scripts/prepare_osm_network_release.py

@@ -7,6 +7,7 @@
 import folium
 import geopandas as gpd
 import numpy as np
+import pandas as pd
 import pypsa
 from shapely.wkt import loads
 
@@ -78,34 +79,29 @@
 ]
 
 
-def export_clean_csv(df, columns, output_file):
+def export_clean_csv(
+    df: pd.DataFrame, columns: list[str], output_file: str, rename_idx: str
+):
     """
     Export a cleaned DataFrame to a CSV file.
 
     Args:
         df (pandas.DataFrame): The DataFrame to be exported.
         columns (list): A list of column names to include in the exported CSV file.
         output_file (str): The path to the output CSV file.
+        rename_idx (str): The name to use for renaming the index column.
 
     Returns:
         None
     """
-    columns = [col for col in columns if col in df.columns]
     rename_dict = {
-        "Bus": "bus_id",
-        "Line": "line_id",
-        "Link": "link_id",
-        "Transformer": "transformer_id",
         "v_nom": "voltage",
         "num_parallel": "circuits",
     }
 
-    if "converter_id" in columns:
-        rename_dict["Link"] = "converter_id"
-
-    df.reset_index().rename(columns=rename_dict).loc[:, columns].replace(
-        {True: "t", False: "f"}
-    ).to_csv(output_file, index=False, quotechar="'")
+    df.rename_axis(index=rename_idx).reset_index().rename(columns=rename_dict).loc[
+        :, columns
+    ].replace({True: "t", False: "f"}).to_csv(output_file, index=False, quotechar="'")
 
     return None
 
@@ -131,22 +127,20 @@ def create_geometries(network, is_converter, crs=GEO_CRS):
     """
 
     buses_cols = [
-        "Bus",
         "v_nom",
         "dc",
         "symbol",
         "under_construction",
         "tags",
         "geometry",
     ]
-    buses = network.buses.reset_index()[
+    buses = network.buses[
         [c for c in buses_cols if c in network.buses.columns]
-    ]
+    ].reset_index()
     buses["geometry"] = buses.geometry.apply(lambda x: loads(x))
     buses = gpd.GeoDataFrame(buses, geometry="geometry", crs=crs)
 
     lines_cols = [
-        "Line",
         "bus0",
         "bus1",
         "v_nom",
@@ -163,15 +157,14 @@ def create_geometries(network, is_converter, crs=GEO_CRS):
         "tags",
         "geometry",
     ]
-    lines = network.lines.reset_index()[
+    lines = network.lines[
         [c for c in lines_cols if c in network.lines.columns]
-    ]
+    ].reset_index()
     # Create shapely linestring from geometry column
     lines["geometry"] = lines.geometry.apply(lambda x: loads(x))
     lines = gpd.GeoDataFrame(lines, geometry="geometry", crs=crs)
 
     links_cols = [
-        "Link",
         "bus0",
         "bus1",
         "v_nom",
@@ -184,16 +177,15 @@ def create_geometries(network, is_converter, crs=GEO_CRS):
     ]
     links = (
         network.links[~is_converter]
-        .reset_index()
         .rename(columns={"voltage": "v_nom"})[
             [c for c in links_cols if c in network.links.columns]
         ]
+        .reset_index()
     )
     links["geometry"] = links.geometry.apply(lambda x: loads(x))
     links = gpd.GeoDataFrame(links, geometry="geometry", crs=crs)
 
     converters_cols = [
-        "Link",
         "bus0",
         "bus1",
         "v_nom",
@@ -202,26 +194,25 @@ def create_geometries(network, is_converter, crs=GEO_CRS):
     ]
     converters = (
         network.links[is_converter]
-        .reset_index()
         .rename(columns={"voltage": "v_nom"})[
             [c for c in converters_cols if c in network.links.columns]
         ]
+        .reset_index()
     )
     converters["geometry"] = converters.geometry.apply(lambda x: loads(x))
     converters = gpd.GeoDataFrame(converters, geometry="geometry", crs=crs)
 
     transformers_cols = [
-        "Transformer",
         "bus0",
         "bus1",
         "voltage_bus0",
         "voltage_bus1",
         "s_nom",
         "geometry",
     ]
-    transformers = network.transformers.reset_index()[
+    transformers = network.transformers[
         [c for c in transformers_cols if c in network.transformers.columns]
-    ]
+    ].reset_index()
     transformers["geometry"] = transformers.geometry.apply(lambda x: loads(x))
     transformers = gpd.GeoDataFrame(transformers, geometry="geometry", crs=crs)
 
@@ -287,19 +278,21 @@ def create_geometries(network, is_converter, crs=GEO_CRS):
 
     # Export to clean csv for release
     logger.info(f"Exporting {len(network.buses)} buses to %s", snakemake.output.buses)
-    export_clean_csv(network.buses, BUSES_COLUMNS, snakemake.output.buses)
+    export_clean_csv(network.buses, BUSES_COLUMNS, snakemake.output.buses, "bus_id")
 
     logger.info(
         f"Exporting {len(network.transformers)} transformers to %s",
         snakemake.output.transformers,
     )
     export_clean_csv(
-        network.transformers, TRANSFORMERS_COLUMNS, snakemake.output.transformers
+        network.transformers,
+        TRANSFORMERS_COLUMNS,
+        snakemake.output.transformers,
+        "transformer_id",
     )
 
     logger.info(f"Exporting {len(network.lines)} lines to %s", snakemake.output.lines)
-    export_clean_csv(network.lines, LINES_COLUMNS, snakemake.output.lines)
-
+    export_clean_csv(network.lines, LINES_COLUMNS, snakemake.output.lines, "line_id")
     # Boolean that specifies if link element is a converter
     is_converter = network.links.index.str.startswith("conv") == True
 
@@ -308,15 +301,18 @@ def create_geometries(network, is_converter, crs=GEO_CRS):
         snakemake.output.links,
     )
     export_clean_csv(
-        network.links[~is_converter], LINKS_COLUMNS, snakemake.output.links
+        network.links[~is_converter], LINKS_COLUMNS, snakemake.output.links, "link_id"
     )
 
     logger.info(
         f"Exporting {len(network.links[is_converter])} converters to %s",
         snakemake.output.converters,
     )
     export_clean_csv(
-        network.links[is_converter], CONVERTERS_COLUMNS, snakemake.output.converters
+        network.links[is_converter],
+        CONVERTERS_COLUMNS,
+        snakemake.output.converters,
+        "converter_id",
     )
 
     ## Create interactive map