From 1cf85b0be889683935ec3898b5eac494aff5dd1b Mon Sep 17 00:00:00 2001
From: Bobby Xiong <36541459+bobbyxng@users.noreply.github.com>
Date: Wed, 14 Aug 2024 14:22:03 +0200
Subject: [PATCH] Fix for Corsica in simplify_network: Include local substation
 (#1215)

* Fixed simplify_network.py to handle more complex topologies, i.e. if two links are connected to nearby AC buses separated by an AC line.

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Added fix for Corsica node: If region containing Corsica is modelled, include substation on Corsica as supernode (end point).

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
---
 scripts/simplify_network.py | 45 +++++++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py
index 2b759a732..c54e3418b 100644
--- a/scripts/simplify_network.py
+++ b/scripts/simplify_network.py
@@ -306,6 +306,14 @@ def split_links(nodes):
 
         seen = set()
 
+        # Corsica substation
+        node_corsica = find_closest_bus(
+            n,
+            x=9.44802,
+            y=42.52842,
+            tol=2000,  # Tolerance needed to only return the bus if the region is actually modelled
+        )
+
         # Supernodes are endpoints of links, identified by having lass then two neighbours or being an AC Bus
         # An example for the latter is if two different links are connected to the same AC bus.
         supernodes = {
@@ -314,6 +322,7 @@ def split_links(nodes):
             if (
                 (len(G.adj[m]) < 2 or (set(G.adj[m]) - nodes))
                 or (n.buses.loc[m, "carrier"] == "AC")
+                or (m == node_corsica)
             )
         }
 
@@ -530,6 +539,42 @@ def cluster(
     return clustering.network, clustering.busmap
 
 
+def find_closest_bus(n, x, y, tol=2000):
+    """
+    Find the index of the closest bus to the given coordinates within a specified tolerance.
+    Parameters:
+        n (pypsa.Network): The network object.
+        x (float): The x-coordinate (longitude) of the target location.
+        y (float): The y-coordinate (latitude) of the target location.
+        tol (float): The distance tolerance in meters. Default is 2000 meters.
+
+    Returns:
+        int: The index of the closest bus to the target location within the tolerance.
+             Returns None if no bus is within the tolerance.
+    """
+    # Conversion factors
+    meters_per_degree_lat = 111139  # Meters per degree of latitude
+    meters_per_degree_lon = 111139 * np.cos(
+        np.radians(y)
+    )  # Meters per degree of longitude at the given latitude
+
+    x0 = np.array(n.buses.x)
+    y0 = np.array(n.buses.y)
+
+    # Calculate distances in meters
+    dist = np.sqrt(
+        ((x - x0) * meters_per_degree_lon) ** 2
+        + ((y - y0) * meters_per_degree_lat) ** 2
+    )
+
+    # Find the closest bus within the tolerance
+    min_dist = dist.min()
+    if min_dist <= tol:
+        return n.buses.index[dist.argmin()]
+    else:
+        return None
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake