Pretty sure this can be done with clustering. Using AgglomerativeClustering
from sklearn
seems to do the trick. With it you can set a distance threshold to be your desired radius for grouping. Then set linkage='single'
to cluster via the minimum distances between all potential group members in a "greedy" fashion.
from sklearn.cluster import AgglomerativeClustering
def cluster_shapes_by_distance(geodf, distance):
"""
Make groups for all shapes within a defined distance. For a shape to be
excluded from a group, it must be greater than the defined distance
from *all* shapes in the group.
Distances are calculated using shape centroids.
Parameters
----------
geodf : data.frame
A geopandas data.frame of polygons. Should be a projected CRS where the
unit is in meters.
distance : float
Maximum distance between elements. In meters.
Returns
-------
np.array
Array of numeric labels assigned to each row in geodf.
"""
assert geodf.crs.is_projected, 'geodf should be a projected crs with meters as the unit'
centers = [p.centroid for p in geodf.geometry]
centers_xy = [[c.x, c.y] for c in centers]
cluster = AgglomerativeClustering(n_clusters=None,
linkage='single',
affinity='euclidean',
distance_threshold=distance)
cluster.fit(centers_xy)
return cluster.labels_
Here's an example using a dataset from the fiona packages
import geopandas as gpd
import pooch
shapefile_url = 'https://github.com/Toblerity/Fiona/raw/master/tests/data/coutwildrnp.zip'
file_paths = pooch.retrieve(
shapefile_url,
known_hash=None,
processor=pooch.Unzip()
)
shp_file = [f for f in file_paths if f.split('.')[-1] == 'shp'][0]
# Convert from 4326 to a local UTM CRS so units are in meters
shapes = gpd.read_file(shp_file).to_crs(epsg=32613)
shapes['group'] = cluster_shapes_by_distance(shapes, distance=125000) # 125km
# back to 4326 for plotting
shapes = shapes.to_crs(epsg=4326)
# Make the group a string so it's treated as categorical for plotting
shapes['group'] = shapes.group.astype(str)
shapes.plot(column='group', legend=True)
