Removing Duplicate Spatial Points with Tolerance Thresholds in Python

pip install geopandas>=0.12 shapely>=2.0 scipy>=1.9 numpy>=1.24

import shapely
print(shapely.geos_version_string)  # expect 3.11.x or later

import logging
import geopandas as gpd
import numpy as np
from scipy.spatial import cKDTree
from shapely.geometry import Point

logger = logging.getLogger(__name__)


def deduplicate_points_tolerance(
    gdf: gpd.GeoDataFrame,
    tolerance_meters: float,
    method: str = "first",
) -> gpd.GeoDataFrame:
    """
    Remove near-duplicate point features within a tolerance threshold.

    Parameters
    ----------
    gdf : GeoDataFrame
        Input point layer; must contain only Point geometries.
    tolerance_meters : float
        Maximum centre-to-centre distance (m) to treat two points as duplicates.
    method : {'first', 'centroid'}
        'first'    — keep the row with the lowest original index per cluster.
        'centroid' — replace geometry with the arithmetic mean of cluster coords
                     and copy non-spatial attributes from the first row.

    Returns
    -------
    GeoDataFrame
        Deduplicated point layer in the same CRS as the input.

    Raises
    ------
    ValueError
        If the GeoDataFrame contains non-Point geometries or an unknown method.
    """
    if gdf.empty:
        logger.info("deduplicate_points_tolerance: empty GeoDataFrame, returning as-is")
        return gdf.copy()

    if not all(gdf.geometry.geom_type == "Point"):
        raise ValueError(
            "deduplicate_points_tolerance requires a Point-only GeoDataFrame. "
            "Found mixed or non-Point geometry types."
        )

    original_crs = gdf.crs

    # Project to metric CRS for accurate Euclidean distance in metres.
    # estimate_utm_crs() picks the locally best UTM zone (GeoPandas ≥ 0.12).
    if gdf.crs is None or gdf.crs.is_geographic:
        metric_crs = gdf.estimate_utm_crs() if gdf.crs is not None else "EPSG:3857"
        gdf = gdf.to_crs(metric_crs)
        logger.debug("Projected to %s for distance calculations", gdf.crs.to_string())

    # Shapely 2.0 vectorized coordinate extraction — no .apply()
    coords = np.column_stack([gdf.geometry.x.to_numpy(), gdf.geometry.y.to_numpy()])

    # O(n log n) spatial index
    tree = cKDTree(coords)

    # Each entry: list of indices within tolerance_meters of point i
    neighbor_lists = tree.query_ball_point(coords, r=tolerance_meters)

    # ------------------------------------------------------------------ #
    # Union-Find with path compression
    # Neighbour relations are transitive: if A~B and B~C, then {A,B,C}
    # form one duplicate cluster even if dist(A,C) > tolerance.
    # ------------------------------------------------------------------ #
    parent = np.arange(len(gdf), dtype=np.int64)

    def find(i: int) -> int:
        """Iterative find with path halving."""
        while parent[i] != i:
            parent[i] = parent[parent[i]]  # path halving
            i = parent[i]
        return i

    def union(i: int, j: int) -> None:
        ri, rj = find(i), find(j)
        if ri != rj:
            parent[ri] = rj

    for i, nbrs in enumerate(neighbor_lists):
        for j in nbrs:
            if i < j:
                union(i, j)

    # Group original indices by cluster root
    clusters: dict[int, list[int]] = {}
    for i in range(len(gdf)):
        clusters.setdefault(find(i), []).append(i)

    # Build one representative row per cluster
    keep_rows = []
    for indices in clusters.values():
        if method == "first":
            keep_rows.append(gdf.iloc[indices[0]])
        elif method == "centroid":
            sub = gdf.iloc[indices]
            mean_x = float(sub.geometry.x.mean())
            mean_y = float(sub.geometry.y.mean())
            row = sub.iloc[0].copy()
            row.geometry = Point(mean_x, mean_y)
            keep_rows.append(row)
        else:
            raise ValueError(f"Unknown method '{method}'. Choose 'first' or 'centroid'.")

    result = gpd.GeoDataFrame(keep_rows, crs=gdf.crs).reset_index(drop=True)

    # Reproject back to original CRS if we changed it
    if result.crs != original_crs and original_crs is not None:
        result = result.to_crs(original_crs)

    removed = len(gdf) - len(result)
    logger.info(
        "deduplicate_points_tolerance: removed %d duplicate(s) from %d input points "
        "(tolerance=%.3f m, method=%s)",
        removed, len(gdf), tolerance_meters, method,
    )
    return result

import numpy as np
from scipy.spatial import cKDTree

# Quick threshold diagnostic — run on a sample before processing the full dataset
sample_coords = coords[:5000]
sample_tree = cKDTree(sample_coords)
neighbor_lists_sample = sample_tree.query_ball_point(sample_coords, r=tolerance_meters)

cluster_radii = []
for i, nbrs in enumerate(neighbor_lists_sample):
    if len(nbrs) > 1:
        dists = np.linalg.norm(sample_coords[nbrs] - sample_coords[i], axis=1)
        cluster_radii.append(float(dists.max()))

if cluster_radii:
    print(f"Max intra-cluster radius: {max(cluster_radii):.3f} m")
    print(f"Median intra-cluster radius: {np.median(cluster_radii):.3f} m")
    print(f"Clusters with >1 member: {len(cluster_radii)}")