Automated polygon splitting based on points is a common task in spatial analysis and geographic information systems (GIS). It is often used in various applications such as urban planning, resource allocation, and spatial clustering. The goal is to divide a larger geographic area (represented as a polygon) into smaller, more manageable areas based on specific criteria or points of interest.
In this tutorial, we will walk through the process of automating polygon splitting using Voronoi diagrams and clustering techniques in Python. Voronoi diagrams help to partition a space into regions based on proximity to a set of points, which aligns well with the task of splitting polygons based on points of interest.
Requirements
- Python (3.x preferred)
- Geopandas
- Shapely
- Numpy
- Scikit-learn (for KMeans clustering)
Understanding Polygon Splitting with K-Means Clustering and Voronoi Diagrams
Scenario:
- You have a large polygon representing a geographical area.
- You want to subdivide this area into smaller, more manageable regions based on the distribution of points of interest within it.
Approach:
- Generate Random Points:
- Clustering Points with K-Means:
- Separating Points into Clusters:
- Constructing Convex Hulls:
- Calculating Centroids:
- Creating Voronoi Diagrams:
- Aligning Voronoi Polygons with Main Polygon:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 |
import geopandas as gpd from shapely.geometry import Polygon, Point, MultiPoint, MultiPolygon, LineString import numpy as np import random from sklearn.cluster import KMeans from shapely import wkt import math from shapely.ops import voronoi_diagram # WKT for main polygon for which random points needs to be generated wkt_str = 'POLYGON ((-85.00037094509452 34.43022863177922, -84.66465765300998 34.59550286788238, -84.07328577695337 34.51028333989169, -83.74531971468618 34.3166025944583, -83.52323245992255 33.96797725267821, -83.59295752827857 33.56253889223765, -83.87960503151999 33.35078127723049, -84.17141735463963 33.23715523990957, -84.54844920574995 33.28105620887447, -84.86092080838247 33.41017670583006, -84.86092080838247 33.41017670583006, -85.09333770290255 33.62451673077634, -85.2353702495537 33.8569336252964, -85.2353702495537 34.13841630865959, -85.00037094509452 34.43022863177922))' polygon = wkt.loads(wkt_str) """ Create Random Points """ # Define the number of random points you want num_points = 10000 # Generate random points inside the polygon min_x, min_y, max_x, max_y = polygon.bounds points = [] while len(points) < num_points: random_x = np.random.uniform(min_x, max_x) random_y = np.random.uniform(min_y, max_y) point = Point(random_x, random_y) if polygon.contains(point): points.append(point) print(f"Generated {len(points)} random points inside the polygon.") # You can now use random_points list containing the generated points # Converting point geometry to plain coordinates points = [[point.x,point.y] for point in points] n_clusters=20 kmeans = KMeans(n_clusters=n_clusters) kmeans.fit(points) # Get cluster labels labels = kmeans.labels_ # Create empty lists to store points for each cluster cluster_points = [[] for _ in range(n_clusters)] # Separate points into clusters based on labels for point_idx, label in enumerate(labels): cluster_points[label].append(points[point_idx]) # Print number of points in each cluster for i, points in enumerate(cluster_points): print(f"Cluster {i+1} has {len(points)} points.") # Get convex_hulls of clusters convex_hulls = [MultiPoint(x).convex_hull for x in cluster_points] # Get Centroid of convex_hulls convex_hulls_centroids = [x.centroid for x in convex_hulls] # Get voronoi for given centroids voronoi_polygons = list(voronoi_diagram(MultiPoint([centroid for centroid in convex_hulls_centroids])).geoms) # Perform Intersection with main polygon for alignment aligned_polygon = MultiPolygon([x.buffer(0).intersection(polygon) for x in voronoi_polygons]) print(aligned_polygon) |
I hope this tutorial will create a good foundation for you. If you want tutorials on another GIS topic or you have any queries, please send an mail at contact@spatial-dev.guru.