3

I have polygons as shapefile sprawl on very large area in Brazil, with many "blank" area in the middle:
enter image description here

I would like to split this shapefile to smaller groups of polygons based on their distance from each other.

For example, if I set the distance to be 70 km, than plots that are located inside this range of 70 meters will be grouped into one shapefile/one class so I know they are within 70km radius from each other.

This is somehow resembling a little the near tool by ArcMap. However, I want to do that in Python, without ArcMap.

I'm looking here for idea/tools/ solutions that people used if they had the same challenge.

4
  • Should all be withing 7 km of eachother or can two on each side be further away? Which python library do you want to use? Try something
    – BERA
    Dec 19, 2021 at 13:21
  • @BERA I want to clusterize the polygons based on that range. Given plots x,y,z if plot y is 50 km away from plot x, and plot z is 50km away from th eother side (100km between z and z), they will still be clustered in the same cluster as they are less than 70km away from plot inside the cluster ((x). Is it more clear?
    – ReutKeller
    Dec 19, 2021 at 13:51
  • @BERA I haven't tried as I'm looking for idea. any python library is welcome.
    – ReutKeller
    Dec 19, 2021 at 13:51

1 Answer 1

4

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)

enter image description here

Not the answer you're looking for? Browse other questions tagged or ask your own question.