From a05881479ccd363cf729e74e9947dcc45fbede33 Mon Sep 17 00:00:00 2001 From: Fabian Neumann Date: Mon, 28 Mar 2022 15:17:55 +0200 Subject: [PATCH 1/5] 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 From d5db3b8d8060dfd7d4a33c5749a4c2c70ef64864 Mon Sep 17 00:00:00 2001 From: Fabian Neumann Date: Tue, 7 Jun 2022 10:57:01 +0200 Subject: [PATCH 2/5] Update scripts/build_bus_regions.py Co-authored-by: Fabian Hofmann --- scripts/build_bus_regions.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/scripts/build_bus_regions.py b/scripts/build_bus_regions.py index 4f2369b6..382a32e8 100644 --- a/scripts/build_bus_regions.py +++ b/scripts/build_bus_regions.py @@ -117,9 +117,7 @@ def voronoi_partition_pts(points, outline, no_multipolygons=False): return poly polygons = [demultipolygon(poly) for poly in polygons] - polygons_arr = np.empty((len(polygons),), 'object') - polygons_arr[:] = polygons - return polygons_arr + return np.array(polygons, dtype=object) if __name__ == "__main__": From aa867cb70489041cd7fdbed61ebec85a52e38ed6 Mon Sep 17 00:00:00 2001 From: Fabian Hofmann Date: Tue, 7 Jun 2022 15:00:57 +0200 Subject: [PATCH 3/5] Update scripts/build_bus_regions.py Co-authored-by: Fabian Neumann --- scripts/build_bus_regions.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/scripts/build_bus_regions.py b/scripts/build_bus_regions.py index 382a32e8..89765ed9 100644 --- a/scripts/build_bus_regions.py +++ b/scripts/build_bus_regions.py @@ -107,15 +107,6 @@ def voronoi_partition_pts(points, outline, no_multipolygons=False): 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] return np.array(polygons, dtype=object) From 97fbf77ff8ddb7da8bb75602b73853e319854cb5 Mon Sep 17 00:00:00 2001 From: Fabian Hofmann Date: Tue, 7 Jun 2022 15:01:18 +0200 Subject: [PATCH 4/5] Update scripts/build_bus_regions.py --- scripts/build_bus_regions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/build_bus_regions.py b/scripts/build_bus_regions.py index 89765ed9..b6d2d129 100644 --- a/scripts/build_bus_regions.py +++ b/scripts/build_bus_regions.py @@ -62,7 +62,7 @@ def save_to_geojson(s, fn): s.to_file(fn, driver='GeoJSON', schema=schema) -def voronoi_partition_pts(points, outline, no_multipolygons=False): +def voronoi_partition_pts(points, outline): """ Compute the polygons of a voronoi partition of `points` within the polygon `outline`. Taken from From bdd094d796b3f896b42ba8d6fb27c0093f467a67 Mon Sep 17 00:00:00 2001 From: Fabian Hofmann Date: Tue, 7 Jun 2022 15:01:40 +0200 Subject: [PATCH 5/5] Update scripts/build_bus_regions.py --- scripts/build_bus_regions.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/scripts/build_bus_regions.py b/scripts/build_bus_regions.py index b6d2d129..8003d370 100644 --- a/scripts/build_bus_regions.py +++ b/scripts/build_bus_regions.py @@ -71,8 +71,6 @@ def voronoi_partition_pts(points, outline): ---------- 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]