Add functionality to merge each isolated bus to the backbone network

scripts/simplify_network.py

@@ -89,6 +89,7 @@
 import sys
 from functools import reduce
 
+import geopandas as gpd
 import numpy as np
 import pandas as pd
 import pypsa
@@ -104,6 +105,8 @@
 )
 from pypsa.io import import_components_from_dataframe, import_series_from_dataframe
 from scipy.sparse.csgraph import connected_components, dijkstra
+from scipy.spatial import cKDTree
+from shapely.geometry import Point
 
 sys.settrace
 
@@ -738,6 +741,118 @@ def drop_isolated_nodes(n, threshold):
     return n
 
 
+def merge_into_network(n, aggregation_strategies=dict()):
+    """
+    Find isolated nodes in the network and merge those of them which have load
+    value below than a specified threshold into a single isolated node which
+    represents all the remote generation.
+
+    Parameters
+    ----------
+    n : PyPSA.Network
+        Original network
+    threshold : float
+        Load power used as a threshold to merge isolated nodes
+
+    Returns
+    -------
+    modified network
+    """
+    # keep original values of the overall load and generation in the network
+    # to track changes due to drop of buses
+    generators_mean_origin = n.generators.p_nom.mean()
+    load_mean_origin = n.loads_t.p_set.mean().mean()
+
+    n.determine_network_topology()
+
+    # duplicated sub-networks mean that there is at least one interconnection between buses
+    i_islands = n.buses[
+        (~n.buses.duplicated(subset=["sub_network"], keep=False))
+        & (n.buses.carrier == "AC")
+    ].index
+
+    # TODO filtering may be applied to decide if the isolated buses should be fetched
+    # # isolated buses with load below than a specified threshold should be merged
+    # i_load_islands = n.loads_t.p_set.columns.intersection(i_islands)
+    # i_suffic_load = i_load_islands[
+    #     n.loads_t.p_set[i_load_islands].mean(axis=0) <= threshold
+    # ]
+
+    # return the original network if no isolated nodes are detected
+    if len(i_islands) == 0:
+        return n, n.buses.index.to_series()
+
+    i_connected = n.buses.loc[n.buses.carrier == "AC"].index.difference(i_islands)
+
+    points_buses = np.array(
+        list(zip(n.buses.loc[i_connected].x, n.buses.loc[i_connected].y))
+    )
+    islands_points = np.array(
+        list(zip(n.buses.loc[i_islands].x, n.buses.loc[i_islands].y))
+    )
+
+    gds_buses = gpd.GeoSeries(map(Point, points_buses))
+    gds_islands = gpd.GeoSeries(map(Point, islands_points))
+
+    gdf_buses = gpd.GeoDataFrame(geometry=gds_buses)
+    gdf_islands = gpd.GeoDataFrame(geometry=gds_islands)
+
+    gdf_map = gpd.sjoin_nearest(gdf_islands, gdf_buses, how="left")
+
+    nearest_bus_list = [
+        n.buses.loc[(n.buses.x == x) & (n.buses.y == y)]
+        for x, y in zip(gdf_map["geometry"].x, gdf_map["geometry"].y)
+    ]
+    nearest_bus_df = pd.concat(nearest_bus_list)
+
+    # each isolated node should be mapped into the closes non-isolated node
+    map_isolated_node_by_country = (
+        n.buses.assign(bus_id=n.buses.index)
+        .loc[nearest_bus_df.index]
+        .groupby("country")["bus_id"]
+        .first()
+        .to_dict()
+    )
+    isolated_buses_mapping = n.buses.loc[i_islands, "country"].replace(
+        map_isolated_node_by_country
+    )
+    busmap = (
+        n.buses.index.to_series()
+        .replace(isolated_buses_mapping)
+        .astype(str)
+        .rename("busmap")
+    )
+
+    # return the original network if no changes are detected
+    if (busmap.index == busmap).all():
+        return n, n.buses.index.to_series()
+
+    bus_strategies, generator_strategies = get_aggregation_strategies(
+        aggregation_strategies
+    )
+
+    clustering = get_clustering_from_busmap(
+        n,
+        busmap,
+        bus_strategies=bus_strategies,
+        aggregate_generators_weighted=True,
+        aggregate_generators_carriers=None,
+        aggregate_one_ports=["Load", "StorageUnit"],
+        line_length_factor=1.0,
+        generator_strategies=generator_strategies,
+        scale_link_capital_costs=False,
+    )
+
+    load_mean_final = n.loads_t.p_set.mean().mean()
+    generators_mean_final = n.generators.p_nom.mean()
+
+    logger.info(
+        f"Fetched {len(islands_points)} isolated buses into the network. Load attached to a single bus with discrepancies of {(100 * ((load_mean_final - load_mean_origin)/load_mean_origin)):2.1E}% and {(100 * ((generators_mean_final - generators_mean_origin)/generators_mean_origin)):2.1E}% for load and generation capacity, respectively"
+    )
+
+    return clustering.network, busmap
+
+
 def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
     """
     Find isolated nodes in the network and merge those of them which have load
@@ -981,6 +1096,13 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
         )
         busmaps.append(merged_nodes_map)
 
+    # TODO Add a configuration option
+    n, fetched_nodes_map = merge_into_network(
+        n,
+        aggregation_strategies=aggregation_strategies,
+    )
+    busmaps.append(fetched_nodes_map)
+
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
     n.export_to_netcdf(snakemake.output.network)