From a05881479ccd363cf729e74e9947dcc45fbede33 Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Mon, 28 Mar 2022 15:17:55 +0200
Subject: [PATCH] build_bus_regions: move voronoi partition from vresutils to
 script

---
 scripts/build_bus_regions.py | 65 ++++++++++++++++++++++++++++++++++--
 1 file changed, 63 insertions(+), 2 deletions(-)

diff --git a/scripts/build_bus_regions.py b/scripts/build_bus_regions.py
index d91d0575..4f2369b6 100644
--- a/scripts/build_bus_regions.py
+++ b/scripts/build_bus_regions.py
@@ -47,9 +47,10 @@ from _helpers import configure_logging
 import pypsa
 import os
 import pandas as pd
+import numpy as np
 import geopandas as gpd
-
-from vresutils.graph import voronoi_partition_pts
+from shapely.geometry import Polygon
+from scipy.spatial import Voronoi
 
 logger = logging.getLogger(__name__)
 
@@ -61,6 +62,66 @@ def save_to_geojson(s, fn):
     s.to_file(fn, driver='GeoJSON', schema=schema)
 
 
+def voronoi_partition_pts(points, outline, no_multipolygons=False):
+    """
+    Compute the polygons of a voronoi partition of `points` within the
+    polygon `outline`. Taken from
+    https://github.com/FRESNA/vresutils/blob/master/vresutils/graph.py
+    Attributes
+    ----------
+    points : Nx2 - ndarray[dtype=float]
+    outline : Polygon
+    no_multipolygons : bool (default: False)
+        If true, replace each MultiPolygon by its largest component
+    Returns
+    -------
+    polygons : N - ndarray[dtype=Polygon|MultiPolygon]
+    """
+
+    points = np.asarray(points)
+
+    if len(points) == 1:
+        polygons = [outline]
+    else:
+        xmin, ymin = np.amin(points, axis=0)
+        xmax, ymax = np.amax(points, axis=0)
+        xspan = xmax - xmin
+        yspan = ymax - ymin
+
+        # to avoid any network positions outside all Voronoi cells, append
+        # the corners of a rectangle framing these points
+        vor = Voronoi(np.vstack((points,
+                                 [[xmin-3.*xspan, ymin-3.*yspan],
+                                  [xmin-3.*xspan, ymax+3.*yspan],
+                                  [xmax+3.*xspan, ymin-3.*yspan],
+                                  [xmax+3.*xspan, ymax+3.*yspan]])))
+
+        polygons = []
+        for i in range(len(points)):
+            poly = Polygon(vor.vertices[vor.regions[vor.point_region[i]]])
+
+            if not poly.is_valid:
+                poly = poly.buffer(0)
+
+            poly = poly.intersection(outline)
+
+            polygons.append(poly)
+
+    if no_multipolygons:
+        def demultipolygon(poly):
+            try:
+                # for a MultiPolygon pick the part with the largest area
+                poly = max(poly.geoms, key=lambda pg: pg.area)
+            except:
+                pass
+            return poly
+        polygons = [demultipolygon(poly) for poly in polygons]
+
+    polygons_arr = np.empty((len(polygons),), 'object')
+    polygons_arr[:] = polygons
+    return polygons_arr
+
+
 if __name__ == "__main__":
     if 'snakemake' not in globals():
         from _helpers import mock_snakemake