Spatial Clustering Algorithms in Python

Clustering finds dense groups in point data without being told how many to expect — the core unsupervised technique within Spatial Analysis & Advanced Query Techniques, the parent guide this page sits under. It shares the metric-CRS discipline of its siblings Nearest Neighbor & KD-Tree Search and Proximity & Buffer Analysis: every distance-based algorithm here is only meaningful when coordinates live in a projected, metre-based system. This page covers the algorithm surface, the projection pipeline they demand, and the boundary and export steps that turn raw labels into shippable layers; for a head-to-head on the two density methods, see DBSCAN vs HDBSCAN for Spatial Clustering.

From raw points to labelled groups Undifferentiated points on the left pass through a density-based algorithm and emerge on the right partitioned into two coloured clusters plus scattered noise points. Density in, groups out raw points DBSCAN / HDBSCAN metric CRS only clusters + noise
Density-based clustering partitions points into groups and noise — but only meaningfully when coordinates are in a metric CRS.

Architecture & Data Structures

Every clustering algorithm on this page shares one output contract: a flat integer label array, one entry per input point, aligned by position to the rows of the source GeoDataFrame. Positive integers are cluster ids, and -1 is the reserved sentinel for noise — points too isolated to belong anywhere. This position-aligned array is the seam between scikit-learn's numeric world and the spatial world: you compute labels on a plain NumPy coordinate array, then attach them straight back as a column.

How DBSCAN labels each point: core, border, or noise A core point sits at the centre of an epsilon-radius circle that contains at least min_samples neighbours; a border point lies inside that radius but is not itself dense; an isolated point beyond every neighbourhood is labelled minus one for noise. Core, border, noise — the density test eps (m) -1 core — ≥ min_samples within eps border — inside a core's eps noise — labelled -1
DBSCAN classifies each point by its eps neighbourhood: dense cores seed clusters, borders attach to them, and everything else becomes noise (-1).

Density-based methods (DBSCAN, HDBSCAN) build a spatial index internally — a BallTree or KDTree — to answer the neighbourhood queries that define density. That index is the same data structure explored in Nearest Neighbor & KD-Tree Search; clustering is, at its core, a repeated radius query over a tree. The minimal initialization exposes the two parameters that govern everything downstream:

import geopandas as gpd
import numpy as np
from sklearn.cluster import DBSCAN

# sensor_readings: air-quality point stations in a projected metric CRS (metres)
sensor_readings = gpd.read_file("sensor_readings.gpkg").to_crs(epsg=32633)
coords_m = np.column_stack([sensor_readings.geometry.x, sensor_readings.geometry.y])

model = DBSCAN(eps=500, min_samples=8)          # eps in METRES here, not degrees
labels = model.fit_predict(coords_m)            # shape (n,), aligned to sensor_readings

sensor_readings["cluster_id"] = labels          # -1 marks noise
print(sensor_readings["cluster_id"].value_counts().head())

Because the label array is purely positional, never sort or filter the coordinate array between fit_predict and the column assignment — a single reindex silently misaligns every label. Treat the projected coordinate matrix as immutable once it is built.

Density is not the only family worth knowing, and picking the wrong family is a more expensive mistake than mis-tuning a parameter inside the right one. Four families cover almost every spatial question:

OPTICS deserves a specific note because it sits between the first two. It produces a reachability ordering rather than a partition, from which clusters can be extracted at any density afterwards (cluster_method="xi") or at a fixed one (cluster_method="dbscan", which reproduces a DBSCAN run without refitting). That makes it a genuinely useful exploration tool — one fit, many candidate eps values — but its memory grows with the square of the point count when max_eps is left unbounded, so always set max_eps to a plausible ceiling before running it on more than a few tens of thousands of points.

from sklearn.cluster import OPTICS

# One fit, then extract at several densities. max_eps bounds the memory.
optics = OPTICS(min_samples=8, max_eps=2000, cluster_method="xi", xi=0.05)
sensor_readings["optics"] = optics.fit_predict(coords_m)

# The reachability profile IS the diagnostic: valleys are clusters, peaks separate them
reachability = optics.reachability_[optics.ordering_]
print("finite reachability values:", int((reachability < np.inf).sum()))

The choice between the families is a question about the data-generating process, not about accuracy. If features have to be exhaustively assigned, density-based methods are out. If contiguity is a hard requirement, only the constrained hierarchical family delivers it. If the count of groups is fixed by a budget rather than discovered from the data, that is a partitioning problem wearing a clustering costume.

Environment Configuration & Dependency Resolution

The stack chains geopandas for spatial I/O, scikit-learn for the algorithms, and numpy for the coordinate matrix. The one version detail that trips people up: HDBSCAN shipped inside scikit-learn as sklearn.cluster.HDBSCAN only from version 1.3. On older environments you must install the standalone hdbscan package (a separate import with a slightly different API), so pin explicitly rather than assume:

conda install -c conda-forge \
  "geopandas=0.14.*" \
  "scikit-learn=1.4.*" \
  "numpy=1.26.*"
# Confirm the in-tree HDBSCAN is available before relying on it
import sklearn
from importlib.util import find_spec

assert sklearn.__version__ >= "1.3", "Upgrade scikit-learn or install the standalone hdbscan"
print("sklearn", sklearn.__version__, "| standalone hdbscan present:", find_spec("hdbscan") is not None)

Conda-forge is strongly preferred over pip for this stack: geopandas pulls GEOS, GDAL, and PROJ as C-level libraries, and pip wheels can mix incompatible builds. If you also compute spatial weights for validation, add libpysal from the same channel so it shares the identical GEOS build.

Two optional additions are worth installing deliberately rather than discovering you need them halfway through. libpysal supplies the contiguity weights (Queen, Rook, KNN) that a constrained hierarchical clustering needs as its connectivity matrix, and spopt builds on it with purpose-built regionalization solvers. h3-py provides the hexagonal index used for grid binning at scale; it works in lon/lat and returns a stable string cell id per point, which sidesteps the projection question entirely for pure aggregation.

# Only if you need contiguity-constrained regions or hex binning
conda install -c conda-forge "libpysal=4.9.*" "spopt=0.6.*" "h3-py=3.7.*"

One environment detail affects results rather than convenience: scikit-learn's tree-based neighbour searches parallelise through OpenMP, and if n_jobs=-1 is combined with an OpenMP thread pool that is already using every core, the two layers oversubscribe and the run gets slower while memory climbs. Set OMP_NUM_THREADS explicitly in any environment where you also pass n_jobs, and set it before Python starts — the value is read when the native libraries load.

Vectorized Operations & Core Workflow

The canonical workflow is: validate and project, build the coordinate matrix, cluster, attach labels. Projection is not optional decoration — it is the step that makes Euclidean eps mean "metres on the ground". Wrap it once so no downstream call can forget it. Coordinate reference system alignment here follows the same rules as Coordinate Systems with PyProj:

import geopandas as gpd
import numpy as np


def validate_and_prepare_crs(
    gdf: gpd.GeoDataFrame, target_epsg: int = 32633
) -> gpd.GeoDataFrame:
    """Validate CRS, remove invalid geometries, and project to a local metric system."""
    if gdf.crs is None:
        raise ValueError(
            "Input GeoDataFrame lacks a defined CRS. Assign EPSG:4326 if lat/lon."
        )

    # Clean invalid geometries (self-intersections, empty shapes)
    gdf = gdf[~gdf.geometry.is_empty & gdf.geometry.is_valid].copy()
    gdf = gdf.drop_duplicates(subset=["geometry"])

    # Project to local metric CRS for accurate Euclidean distance
    return gdf.to_crs(epsg=target_epsg)


# Usage: df = validate_and_prepare_crs(raw_points, target_epsg=32618)

With clean projected points, the end-to-end run is three lines. DBSCAN remains the industry default for a single, known density scale; HDBSCAN adapts across densities and drops the eps parameter entirely, needing only a minimum cluster size:

from sklearn.cluster import DBSCAN, HDBSCAN

wildlife_sightings = validate_and_prepare_crs(
    gpd.read_file("wildlife_sightings.gpkg"), target_epsg=32633
)
coords_m = np.column_stack(
    [wildlife_sightings.geometry.x, wildlife_sightings.geometry.y]
)

# Fixed-scale density: eps is a hard radius in metres
wildlife_sightings["dbscan"] = DBSCAN(eps=750, min_samples=6).fit_predict(coords_m)

# Adaptive density: no eps, only a minimum viable group size
wildlife_sightings["hdbscan"] = HDBSCAN(min_cluster_size=12).fit_predict(coords_m)

n_groups = wildlife_sightings["hdbscan"].max() + 1
print(f"HDBSCAN found {n_groups} groups, "
      f"{(wildlife_sightings['hdbscan'] == -1).sum()} noise points")

The full parameter trade-off — when a single eps fails and adaptive density wins — is worked through in DBSCAN vs HDBSCAN for Spatial Clustering.

Geometry / Data Processing Details

Two processing details separate a toy example from a production result: how distance is measured when you cannot project, and how point labels become usable polygons.

Haversine on the sphere. When projecting is impractical — a global dataset spanning many UTM zones — cluster directly on the sphere with metric="haversine". The metric expects coordinates in radians, shape (n, 2), in (lat, lon) order (not the (lon, lat)/x, y convention Shapely uses), and eps is then an angular distance you derive by dividing your radius by the Earth radius:

from sklearn.cluster import DBSCAN
import numpy as np


def run_geodesic_dbscan(
    coords_rad: np.ndarray, eps_km: float, min_samples: int = 5
) -> np.ndarray:
    """
    Execute DBSCAN using Haversine metric on radian coordinates.

    Args:
        coords_rad: Array of shape (n, 2) in (lat, lon) radians.
        eps_km: Neighbourhood radius in kilometres.
        min_samples: Minimum points to form a dense region.

    Returns:
        Cluster label array (-1 = noise).
    """
    # Convert km to radians for Haversine (Earth radius ~6371 km)
    eps_rad = eps_km / 6371.0

    model = DBSCAN(
        eps=eps_rad,
        min_samples=min_samples,
        metric="haversine",
        algorithm="ball_tree",
    )
    return model.fit_predict(coords_rad)


# Convert lat/lon degrees to radians before calling:
# coords_rad = np.radians(df[["lat", "lon"]].values)
# labels = run_geodesic_dbscan(coords_rad, eps_km=2.0, min_samples=10)

From labels to boundaries. Downstream consumers want polygons, not scattered points. Group by cluster_id, skip the -1 noise bucket, and wrap each group in a hull. A convex hull is fast and robust; alpha shapes (via the vectorized predicates in Shapely Geometry Operations) hug concave outlines more tightly when a detected cluster is not blob-shaped:

import geopandas as gpd
import numpy as np
from shapely.geometry import MultiPoint


def pipeline_generate_cluster_boundaries(
    gdf: gpd.GeoDataFrame, labels: np.ndarray
) -> gpd.GeoDataFrame:
    """
    Generate cluster geometries from point labels.

    Args:
        gdf: GeoDataFrame of point features (metric CRS).
        labels: Cluster label array from DBSCAN/HDBSCAN (-1 = noise).

    Returns:
        GeoDataFrame with one row per cluster, convex hull geometry.
    """
    gdf = gdf.copy()
    gdf["cluster_id"] = labels
    clusters = []

    for cid, group in gdf.groupby("cluster_id"):
        if cid == -1:  # Noise points — skip
            continue

        points = [(p.x, p.y) for p in group.geometry]
        # convex_hull yields Point/LineString for <3 points, Polygon otherwise
        geom = MultiPoint(points).convex_hull

        clusters.append({
            "cluster_id": int(cid),
            "geometry": geom,
            "count": len(group),
        })

    return gpd.GeoDataFrame(clusters, crs=gdf.crs)


# Performance: build R-tree spatial index for neighbourhood queries
# _ = gdf.sindex  # before clustering to accelerate proximity lookups

Note the two-point degenerate case in the comment: a labelled cluster of fewer than three points yields a Point or LineString, not a Polygon, so a downstream .area call returns 0. Filter or buffer those before treating every row as an area. Collinear points are the sneakier version of the same problem — three or more sightings along a road produce a valid MultiPoint whose convex hull is still a LineString, so the guard has to test the output geometry type rather than count the input points.

Adding time as a third axis. Most spatial event data is really space-time data, and clustering it in space alone merges things that never co-occurred: two vessel positions in the same harbour six weeks apart are not a group. There is no space-time DBSCAN in scikit-learn, but the useful approximation is a coordinate trick. Append the timestamp as a third column, scaled by a conversion factor that states — explicitly — how many metres one second of separation is worth. The algorithm then runs unchanged on a three-dimensional Euclidean space.

import numpy as np
import geopandas as gpd
from sklearn.cluster import DBSCAN

# vessel_pings: AIS positions with a UTC timestamp, projected to metres
vessel_pings = gpd.read_parquet("vessel_pings.parquet").to_crs(epsg=32633)
coords_m = np.column_stack([vessel_pings.geometry.x, vessel_pings.geometry.y])

seconds = vessel_pings["timestamp"].astype("int64") / 1e9

# The conversion IS the modelling decision: one hour of separation costs
# the same as 500 m of separation.
M_PER_SECOND = 500 / 3600.0
space_time = np.column_stack([coords_m, seconds * M_PER_SECOND])

vessel_pings["st_cluster"] = DBSCAN(eps=500, min_samples=6).fit_predict(space_time)
print(vessel_pings["st_cluster"].nunique(), "space-time groups")

Be honest about what that buys and costs. The scaling makes space and time substitutable: a point 500 m away right now and a point in the same berth an hour later are equally close, so the neighbourhood is a ball in a mixed unit space rather than a cylinder. Published ST-DBSCAN variants instead take two independent thresholds and require both to hold, which is stricter and usually what a domain expert means. Use the scaled-axis version when a trade-off between proximity and recency is genuinely acceptable; when it is not, cluster in space and then split every resulting group by time gaps, which is a one-line groupby and needs no new algorithm. Anchor the timestamp column in UTC before scaling — a local-time column crossing a daylight-saving boundary produces a one-hour discontinuity that reads as a genuine gap.

Clustering more points than fit in memory. Both density algorithms hold the full coordinate array plus an index, so the ceiling arrives sooner than the file size suggests. The scaling pattern is to tile the extent, cluster each tile independently with an overlap band at least as wide as eps, and then merge the clusters that share a point between two tiles. The overlap is what makes this correct: any pair of points close enough to be neighbours is seen together in at least one tile, so no cluster can be silently cut in half at a tile edge.

import geopandas as gpd
import numpy as np
from shapely.geometry import box
from sklearn.cluster import DBSCAN

EPS_M, TILE_M, MIN_SAMPLES = 400, 20_000, 6


def cluster_by_tile(points: gpd.GeoDataFrame) -> np.ndarray:
    """DBSCAN per tile with an eps-wide halo, merging clusters that share a point."""
    assert points.crs.is_projected, "tile sizes are in metres"
    sindex = points.sindex
    minx, miny, maxx, maxy = points.total_bounds
    parent: dict[int, int] = {}      # union-find over globally unique cluster ids
    owner: dict[int, int] = {}       # row position -> the id it was first given
    next_id = 0

    def find(a: int) -> int:
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        return a

    for x0 in np.arange(minx, maxx, TILE_M):
        for y0 in np.arange(miny, maxy, TILE_M):
            halo = box(x0, y0, x0 + TILE_M, y0 + TILE_M).buffer(EPS_M)
            rows = sindex.query(halo)
            if len(rows) < MIN_SAMPLES:
                continue
            sub = points.iloc[rows]
            local = DBSCAN(eps=EPS_M, min_samples=MIN_SAMPLES).fit_predict(
                np.column_stack([sub.geometry.x, sub.geometry.y])
            )
            local_to_global: dict[int, int] = {}
            for row, label in zip(rows, local):
                if label == -1:
                    continue
                if label not in local_to_global:
                    local_to_global[label] = next_id
                    parent[next_id] = next_id
                    next_id += 1
                gid = local_to_global[label]
                if row in owner:                       # also seen in a neighbour's halo
                    a, b = find(owner[row]), find(gid)
                    parent[max(a, b)] = min(a, b)      # merge the two tile-local groups
                else:
                    owner[row] = gid

    labels = np.full(len(points), -1, dtype=np.int64)
    for row, gid in owner.items():
        labels[row] = find(gid)
    return labels

The tile size is a memory-versus-duplication trade: the halo re-processes a band of width eps around every tile, which at 400 m and 20 km tiles duplicates roughly 8% of the points, and halving the tile size doubles that overhead. The result is not bit-identical to a single global run — a point that would be core only by borrowing neighbours from three different tiles can end up labelled noise — but the differences are confined to the sparsest fringe, which is the part of the output you were least prepared to defend anyway. If the point set is large enough that even tiling is awkward, push the work to an engine built for it: Scaling with Dask-GeoPandas parallelises exactly this tile-and-merge shape, and DuckDB Spatial Analytics will do the grid-binning alternative in SQL without materialising anything.

CRS Alignment & Projection Pipeline

The single most common failure in spatial clustering is clustering in degrees. A DBSCAN(eps=0.01) on raw WGS84 (EPSG:4326) coordinates is not "roughly a kilometre" — one degree of longitude is ~111 km at the equator and collapses toward zero at the poles, so the same eps describes wildly different ground distances across a dataset. Always project first, and pick the projection deliberately. Related transformation gotchas are catalogued in Coordinate Reference System Transformations.

Rather than hard-code a UTM zone, let GeoPandas choose the correct metric zone from the data's own extent — this avoids the classic error of applying a neighbouring zone to a dataset that straddles a boundary:

import geopandas as gpd

# Incidents arrive as lon/lat WGS84 from an API
traffic_incidents = gpd.read_file("traffic_incidents.geojson")
assert traffic_incidents.crs.to_epsg() == 4326

# Let the data pick its own UTM zone (metric, ~metre accuracy locally)
metric_crs = traffic_incidents.estimate_utm_crs()
print("Selected", metric_crs.to_epsg(), metric_crs.name)

traffic_incidents_m = traffic_incidents.to_crs(metric_crs)
assert traffic_incidents_m.crs.is_projected, "Must be projected before clustering"

Three CRS rules govern this pipeline. First, never use Web Mercator (EPSG:3857) for metric clustering — its distances distort by the secant of latitude, so an eps calibrated in London is wrong in Oslo. Use a local UTM zone or a national grid instead. Second, when you fall back to Haversine, remember the axis-order flip: PyProj-era transforms honour authority axis order, and the Haversine metric wants (lat, lon) radians, the reverse of Shapely's (x, y). Third, assert crs.is_projected immediately before the fit_predict call — a cheap guard that catches the degrees-clustering bug at the exact line it would otherwise corrupt.

The ground distance an eps of 0.01 degrees describes, by latitude Two panels of horizontal bars compare the same neighbourhood radius in two coordinate systems. On the left, an eps of 0.01 degrees in EPSG:4326 shrinks steadily with latitude: about 1113 metres at the equator, 964 metres at 30 degrees north, 787 metres at 45, 557 metres at 60 and 381 metres at 70. On the right, an eps of 1000 metres in a local UTM zone gives an identical bar at every latitude. The same parameter therefore means five different things in degrees and one consistent thing in a projected metric system. One eps, five different ground distances DBSCAN(eps=0.01) on EPSG:4326 Euclidean distance on degrees 0° equator 1113 m 30° N 964 m 45° N 787 m 60° N 557 m 70° N 381 m DBSCAN(eps=1000) on a UTM zone Euclidean distance on metres 0° equator 1000 m 30° N 1000 m 45° N 1000 m 60° N 1000 m 70° N 1000 m A radius in degrees is a different neighbourhood in every row; a radius in metres is the same one everywhere.
An eps expressed in degrees silently rescales with latitude — one reason the projection step must come before fit_predict, never after.

Production Export & Integration

Clustered outputs feed dashboards and tiles, so export in formats built for streaming rather than legacy Shapefiles. FlatGeobuf (.fgb) offers fast spatial-range reads; GeoParquet is the analytics-friendly columnar option — both are covered in Cloud-Native Geospatial Formats. Precompute centroids server-side so the browser never recomputes label positions:

import geopandas as gpd


def prepare_web_export(
    cluster_gdf: gpd.GeoDataFrame, output_stem: str
) -> gpd.GeoDataFrame:
    """Export cluster boundaries to web-optimized formats with precomputed centroids."""
    cluster_gdf = cluster_gdf.copy()

    # Precompute centroids for tile generation & label placement
    centroids = cluster_gdf.geometry.centroid
    cluster_gdf["centroid_lat"] = centroids.y
    cluster_gdf["centroid_lon"] = centroids.x

    # Export to FlatGeobuf (fast spatial queries) and Parquet (analytics)
    fgb_path = f"{output_stem}.fgb"
    parquet_path = f"{output_stem}.parquet"

    cluster_gdf.to_file(fgb_path, driver="FlatGeobuf")
    cluster_gdf.to_parquet(parquet_path)

    print(f"Exported {len(cluster_gdf)} clusters to {fgb_path} and {parquet_path}")
    return cluster_gdf


# Datashader hexbin prep (server-side — requires datashader installed separately)
# import datashader as ds
# cvs = ds.Canvas(plot_width=1000, plot_height=1000)
# agg = cvs.points(df, "lon", "lat", ds.count())

To enrich clusters with administrative or demographic attributes, attribute-join them against reference layers using Geometric Intersections & Overlays. When the point set outgrows a single machine, push both the clustering and the overlay into a database or a partitioned engine: PostGIS Integration with Python for indexed server-side queries, DuckDB Spatial Analytics for in-process columnar out-of-core reads, or Scaling with Dask-GeoPandas for parallel partition-wise processing.

Judging the result before you ship it. A clustering has no ground truth, so the question is not "is it correct?" but "would it come back the same on next week's data?" Two checks answer that cheaply. The first is a stability test: refit on repeated random subsamples and compare each result to the full-data labels with the adjusted Rand index, which measures agreement on which points share a group and is immune to the label integers being renumbered.

import numpy as np
from sklearn.cluster import DBSCAN
from sklearn.metrics import adjusted_rand_score

rng = np.random.default_rng(42)
base = DBSCAN(eps=400, min_samples=6).fit_predict(coords_m)

scores = []
for _ in range(10):
    keep = rng.choice(len(coords_m), size=int(0.8 * len(coords_m)), replace=False)
    resampled = DBSCAN(eps=400, min_samples=6).fit_predict(coords_m[keep])
    scores.append(adjusted_rand_score(base[keep], resampled))

print(f"stability {np.mean(scores):.2f} ± {np.std(scores):.2f}")
# stability 0.91 ± 0.03

Read roughly 0.8 and above as a partition that is a property of the data; below about 0.6 you are reporting an artefact of the parameters, and the honest response is to widen the parameter or say the structure is not there. The second check is choosing the right internal metric. The silhouette score is the reflex, and it is the wrong instrument here: it rewards compact, roughly spherical groups, so it will consistently prefer a KMeans partition over a correct density-based one that found a crescent-shaped hotspot along a river. Use a density-based validity index instead — the standalone hdbscan package ships one — or, more simply, judge by how persistent the groupings stay across parameter values. Whatever you use, publish the noise fraction and how many groups were found alongside the map: a clustering reported without them cannot be reproduced or argued with.

Production performance checklist:

Windows / Platform Edge Cases & Debugging

Clustering pipelines fail in a small number of recognisable ways — most are CRS or environment issues, not algorithm bugs.

Frequently Asked Questions

How do I choose between density-based clustering and hex binning? Ask whether you need object identity. Density clustering gives you things — a hotspot with a boundary, a member count, and a life across time that you can track — at the cost of parameters you must defend. Hex binning gives you a surface: deterministic, parameter-free apart from resolution, trivially joinable, and immune to the argument about whether the groups are real. For dashboards and heat maps, bin. For anything that will be counted, named, or monitored over time, cluster.

Should the clustering happen in Python or in the database? Move it to the data when the data is large and static, keep it in Python when the parameters are still moving. PostGIS ships ST_ClusterDBSCAN and ST_ClusterKMeans as window functions, which cluster millions of rows without ever leaving the server — see PostGIS Integration with Python. The trade is expressiveness: there is no HDBSCAN in SQL, no probability output, and no easy way to iterate on parameters over a remote table.

Do I have to remove duplicate points before clustering? Yes, or at least know how many there are. Duplicate coordinates — several records geocoded to the same building centroid — count as neighbours of each other, so a handful of duplicates can push a location over min_samples and manufacture a spurious grouping from a single address. De-duplicate on rounded coordinates, and if the duplicates are meaningful (three incidents at one address really is denser), weight them explicitly with sample_weight rather than letting the row count decide silently.

What resolution should I round coordinates to before clustering? Match the accuracy of the source, not the precision of the file. Consumer GPS is good to a few metres, a geocoded street address to the length of a building, an administrative centroid to the size of a parcel. Rounding to that accuracy before the fit removes phantom structure at scales the data never measured, and it makes eps values below the rounding step obviously meaningless rather than subtly wrong.

How do I track clusters over time — the same hotspot week to week? Not through the labels, which are renumbered on every fit. Cluster each period independently, build the hull or centroid for each group, then match periods with a spatial join between consecutive hull sets and carry a stable id across the matches. The overlay mechanics are in Geometric Intersections & Overlays; the important discipline is that continuity is something you assert with geometry, never something the algorithm gives you.

Can I cluster polygons instead of points? Yes, by clustering their representative points — geometry.representative_point() rather than centroid, since a centroid can fall outside a concave or ring-shaped polygon. That reduces each feature to one location and discards size and shape, which is fine for "where are the concentrations of parcels" and wrong for "which parcels touch each other". The second question is contiguity, not density: build a Queen weights matrix and take connected components, or dissolve on a shared attribute as in dissolving and aggregating features by attribute.

Related