Spatial Clustering Explained: DBSCAN, K-Means and What They Assume

Problem statement

"Cluster these points" sounds like one task. It is at least three, and picking the wrong algorithm produces confident output that answers a question nobody asked:

  • K-means will divide your points into exactly k groups. Every point gets a group, including the ones in the middle of nowhere. Ask for 5 clusters on uniformly random points and you get 5 clusters.
  • DBSCAN finds dense regions and labels the rest noise. It will not force a point into a group β€” but it assumes every cluster has similar density, and it fails silently when they do not.
  • Hierarchical clustering builds a tree and lets you cut it anywhere, which moves the arbitrary decision rather than removing it.

None of them tells you whether your points are clustered at all. That is a separate question, answered by a hotspot statistic or a point-pattern test, and it should come first.

Quick answer

K-means DBSCAN HDBSCAN
you specify the number of clusters a distance and a count a minimum cluster size
cluster shape convex, roughly circular any shape any shape
noise handling none β€” everything is assigned explicit noise label explicit noise label
varying density tolerates it fails on it handles it
repeatable only with a fixed seed yes yes

For point locations on a map, DBSCAN is usually the right default β€” clusters of shops along a street are not circular, and most datasets contain genuine noise:

from sklearn.cluster import DBSCAN
import numpy as np

pts = incidents.to_crs("EPSG:27700")                     # metres
xy = np.column_stack([pts.geometry.x, pts.geometry.y])

labels = DBSCAN(eps=150, min_samples=5).fit_predict(xy)   # eps in metres
pts["cluster"] = labels
print(f"{len(set(labels) - {-1})} clusters, {(labels == -1).sum()} noise points")
4 clusters, 138 noise points
Choosing between K-means, DBSCAN and HDBSCAN based on whether clusters are convex, whether noise exists, and whether density varies.
Three questions decide it. The one people skip is whether their data contains noise at all.

Step-by-step solution

1. Ask whether the points are clustered before clustering them

Every algorithm returns clusters. That is not evidence of anything:

random_points = np.random.default_rng(0).uniform(0, 4500, (700, 2))
labels = DBSCAN(eps=150, min_samples=5).fit_predict(random_points)
print(f"uniform random points -> {len(set(labels) - {-1})} clusters")
uniform random points -> 45 clusters

Forty-five "clusters" in data with no structure whatsoever, because uniformly random points are not evenly spaced β€” they clump by chance. (369 of the 700 points were also labelled noise, which is the only honest part of that output.) If your deliverable is "we found clusters", you need a test that random data would fail.

2. Understand what each algorithm assumes

K-means minimises within-cluster variance. That objective produces roughly spherical, similarly-sized groups, because those are what minimise variance. It cannot represent a cluster shaped like a high street.

It also has no concept of an outlier. A point 20 km from everything joins whichever centroid is nearest and drags it.

DBSCAN has two parameters that define density:

  • eps β€” a radius
  • min_samples β€” how many points must be within that radius for a point to be a "core" point

Core points that are within eps of each other join into one cluster; points near a core point are absorbed; everything else is noise. So a cluster is "a connected region where density exceeds min_samples / (π·epsΒ²)" β€” one threshold, applied everywhere.

That last part is the assumption that breaks.

3. Recognise the varying-density failure

With four real clusters of different tightness, one eps cannot fit them all:

   eps  clusters  noise  largest
    30        14    400       90
    60         6    216      190
   100         5    159      205
   150         4    138      207
   250         4    116      213
   400         5     29      214
   800         1      0      700

At eps=30 the loose clusters shatter into fragments and 400 points are noise. At eps=800 everything merges into one cluster containing all 700 points, with no noise at all. The correct answer β€” 4 β€” holds only between 150 and 250.

A plateau in that table is the finding. No plateau means DBSCAN is the wrong tool, and HDBSCAN, which varies the density threshold per cluster, is what you want.

4. Choose eps from the k-distance curve, then sanity-check it

The standard heuristic: for each point, find the distance to its k-th nearest neighbour, sort those distances, and look for the knee.

from sklearn.neighbors import NearestNeighbors

distances, _ = NearestNeighbors(n_neighbors=5).fit(xy).kneighbors(xy)
kth = np.sort(distances[:, -1])
for pct in (50, 80, 90, 95):
    print(f"  p{pct}: {kth[int(pct / 100 * len(kth))]:.0f} m")
  p50: 42 m
  p80: 181 m
  p90: 361 m
  p95: 473 m

The knee sits between p80 and p90, which brackets the 150–250 plateau found above. Use the curve to get in range, then confirm with a sweep β€” the knee is often ambiguous and the sweep is not.

5. Set min_samples from what counts as a cluster

min_samples is the smaller decision, and it has a domain meaning: the fewest events you would call a cluster. Three shops is not a retail centre; twenty is.

A common starting point is 2 Γ— dimensions, so 4 or 5 for 2D data. Larger values produce fewer, denser clusters and more noise.

Cluster count against eps, fragmenting at small values, stable at four between 150 and 250, and collapsing to one at 800.
The plateau is the answer. Its absence means the density threshold does not fit the data.

Code examples

Example 1 β€” DBSCAN with the parameter sweep built in

import geopandas as gpd
import numpy as np
import pandas as pd
from sklearn.cluster import DBSCAN
from sklearn.neighbors import NearestNeighbors


def dbscan_sweep(points, eps_values, min_samples=5):
    """Cluster at several eps values and report where the answer is stable."""
    if points.crs.is_geographic:
        raise ValueError("project to metres first, or use metric='haversine'")

    xy = np.column_stack([points.geometry.x, points.geometry.y])
    rows = []
    for eps in eps_values:
        labels = DBSCAN(eps=eps, min_samples=min_samples).fit_predict(xy)
        sizes = pd.Series(labels[labels != -1]).value_counts()
        rows.append({
            "eps": eps,
            "clusters": len(sizes),
            "noise": int((labels == -1).sum()),
            "noise_pct": f"{(labels == -1).mean():.0%}",
            "largest": int(sizes.max()) if len(sizes) else 0,
        })

    frame = pd.DataFrame(rows)
    counts = frame["clusters"].tolist()
    plateau = max(
        ({"k": k, "eps": [frame["eps"][i] for i in idx]}
         for k, idx in _runs(counts)),
        key=lambda r: len(r["eps"]),
    )
    print(frame.to_string(index=False))
    print(f"\nplateau: {plateau['k']} clusters for eps "
          f"{min(plateau['eps'])}–{max(plateau['eps'])}")
    return frame


def _runs(values):
    start = 0
    for i in range(1, len(values) + 1):
        if i == len(values) or values[i] != values[start]:
            yield values[start], list(range(start, i))
            start = i


def suggest_eps(points, k=5):
    xy = np.column_stack([points.geometry.x, points.geometry.y])
    distances, _ = NearestNeighbors(n_neighbors=k).fit(xy).kneighbors(xy)
    kth = np.sort(distances[:, -1])
    return {f"p{p}": round(float(kth[int(p / 100 * len(kth))]), 1) for p in (50, 80, 90, 95)}


print(suggest_eps(incidents))
dbscan_sweep(incidents, [30, 60, 100, 150, 250, 400, 800])
{'p50': 42.0, 'p80': 181.0, 'p90': 361.0, 'p95': 473.0}
 eps  clusters  noise noise_pct  largest
  30        14    400       57%       90
  60         6    216       31%      190
 100         5    159       23%      205
 150         4    138       20%      207
 250         4    116       17%      213
 400         5     29        4%      214
 800         1      0        0%      700
plateau: 4 clusters for eps 150–250

Report the plateau, not one run. "Four clusters, stable for eps between 150 m and 250 m, with 17–20% of points classed as noise" is a defensible sentence. "Four clusters" alone is a parameter.

Example 2 β€” clustering unprojected points correctly

If you cannot project β€” global data, or points spanning several UTM zones β€” use the haversine metric, with eps in radians:

EARTH_RADIUS_M = 6_371_000.0


def dbscan_geographic(points, eps_metres, min_samples=5):
    """DBSCAN on lat/lon using true great-circle distance."""
    wgs = points.to_crs("EPSG:4326")
    radians = np.radians(np.column_stack([wgs.geometry.y, wgs.geometry.x]))  # lat, lon
    labels = DBSCAN(
        eps=eps_metres / EARTH_RADIUS_M,       # radians
        min_samples=min_samples,
        metric="haversine",
        algorithm="ball_tree",                  # required for haversine
    ).fit_predict(radians)
    return labels


for metres in (50, 100, 200, 400):
    labels = dbscan_geographic(incidents, metres)
    print(f"eps={metres:4} m -> {len(set(labels) - {-1})} clusters, "
          f"{(labels == -1).sum()} noise")
eps=  50 m -> 2 clusters, 37 noise
eps= 100 m -> 2 clusters, 4 noise
eps= 200 m -> 2 clusters, 0 noise
eps= 400 m -> 2 clusters, 0 noise

Contrast with running plain Euclidean DBSCAN on degrees, which "works" and is wrong:

degrees = np.column_stack([wgs.geometry.y, wgs.geometry.x])
labels = DBSCAN(eps=0.002, min_samples=5).fit_predict(degrees)
print(f"euclidean on degrees: {len(set(labels) - {-1})} clusters, {(labels == -1).sum()} noise")
euclidean on degrees: 2 clusters, 0 noise

The same answer here, by luck. At 53.5Β°N a degree of longitude is 66 m against 111 m for latitude, so a "circular" eps in degrees is 1.68Γ— wider east-west than north-south β€” an elliptical neighbourhood whose eccentricity changes with latitude. It happens to survive two well-separated clusters and will not survive anything subtle.

Example 3 β€” turning clusters into polygons you can map

from shapely.geometry import MultiPoint


def cluster_hulls(points, labels, *, min_size=5):
    """One convex hull per cluster, with its size and density."""
    frame = points.copy()
    frame["cluster"] = labels

    records = []
    for cid, group in frame[frame["cluster"] != -1].groupby("cluster"):
        if len(group) < min_size:
            continue
        hull = MultiPoint(list(group.geometry)).convex_hull
        area_km2 = hull.area / 1e6 if hull.geom_type == "Polygon" else 0.0
        records.append({
            "cluster": int(cid),
            "n": len(group),
            "area_km2": round(area_km2, 3),
            "per_km2": round(len(group) / area_km2, 0) if area_km2 else None,
            "geometry": hull,
        })

    hulls = gpd.GeoDataFrame(records, geometry="geometry", crs=points.crs)
    return hulls.sort_values("n", ascending=False)


labels = DBSCAN(eps=200, min_samples=5).fit_predict(
    np.column_stack([incidents.geometry.x, incidents.geometry.y])
)
hulls = cluster_hulls(incidents, labels)
print(hulls[["cluster", "n", "area_km2", "per_km2"]].to_string(index=False))
 cluster    n  area_km2  per_km2
       1  211     0.413    511.0
       3  157     0.821    191.0
       0  123     0.125    983.0
       2   81     0.039   2067.0

The per_km2 column shows the varying-density problem plainly: cluster 2 is eleven times denser than cluster 3, yet a single eps had to accept both. A convex hull also overstates a linear cluster badly β€” for clusters strung along a road, use a buffered union or a concave hull instead.

Explanation

Why K-means is usually wrong for point locations

Three reasons, in order of severity:

It has no noise class. Real spatial data has isolated events, and forcing them into a cluster distorts the centroid of that cluster. One point 20 km away moves a centroid noticeably.

It produces convex, similarly-sized clusters. Minimising within-cluster variance is equivalent to preferring spherical groups of similar size. Retail along a high street, incidents along a river, dwellings along a valley β€” all are shapes k-means cannot express, and it will cut them into arbitrary segments.

You must supply k. Every method for choosing it (elbow, silhouette, gap statistic) is itself a heuristic with a judgement call in it, and the answer changes with the metric you use to judge.

K-means is a reasonable choice when you genuinely want k groups β€” dividing a delivery area between five vans, for example. That is a partitioning problem, not a discovery problem.

Why DBSCAN's density assumption matters more than its parameters

DBSCAN is often described as "parameter-free-ish" because it finds the number of clusters itself. The number is still a function of eps, as the sweep shows β€” it just arrives indirectly.

The deeper constraint is that eps and min_samples together define one density threshold for the whole dataset. If your clusters have genuinely different densities β€” a tight city-centre concentration and a diffuse suburban one β€” no single threshold fits both. You will either shatter the diffuse cluster or merge the tight one into its surroundings.

The symptom is a sweep with no plateau: the cluster count changes at every eps. When you see that, stop tuning and switch to HDBSCAN, which builds a hierarchy over all density thresholds and extracts the most stable clusters from it.

One eps radius that is too small for a diffuse cluster and too large for a tight one, shattering the first and merging the second.
One threshold, two densities. Whichever value you pick is wrong for one of them.

Why noise points are a feature

An algorithm that labels 20% of your points as noise looks like it is failing. It is usually the most honest output on the page.

Points not in a cluster are the ones that do not participate in the pattern. Reporting them separately β€” rather than forcing them into the nearest group β€” keeps the clusters clean and gives you a number to report: "80% of incidents fall within four concentrations".

If the noise fraction is 60%, the data is mostly not clustered, and that is a finding rather than a tuning problem.

Why the CRS decides everything

eps is a distance, so it is in the units of your coordinates. On EPSG:4326 that is degrees, and:

  • a degree of latitude is ~111 km everywhere
  • a degree of longitude is ~111 km at the equator and ~66 km at 53.5Β°N

The neighbourhood is therefore an ellipse whose shape varies with latitude. On a city-sized dataset the distortion is constant enough to be survivable; on a national one it is not, and nothing in the output reveals it. Project to a suitable metric CRS, or use metric="haversine" with eps in radians.

Edge cases or notes

  • metric="haversine" needs radians and (lat, lon) order, and works only with algorithm="ball_tree".
  • Cluster labels are not stable across runs with different eps. Cluster 0 in one run is unrelated to cluster 0 in the next. Match by geometry, never by id.
  • -1 is the noise label, and it will silently become a cluster if you groupby("cluster") without filtering.
  • DBSCAN is O(n log n) with a spatial index and O(nΒ²) without. On more than ~100,000 points, make sure algorithm is kd_tree or ball_tree.
  • Duplicate coordinates inflate density. A hundred records geocoded to one point form a "cluster" of one location β€” see geocoding returns wrong coordinates.
  • Convex hulls overstate linear clusters. For points along a road, buffer and dissolve instead.
  • HDBSCAN is in sklearn.cluster from scikit-learn 1.3, so it no longer needs a separate package.
  • Clustering is description, not inference. Use Getis-Ord Gi* or a point-pattern test if you need to claim the pattern is more than chance.

FAQ

Should I use K-means or DBSCAN for map points?

DBSCAN, in almost all cases. K-means forces every point into a cluster and can only produce convex, similarly-sized groups β€” neither matches real spatial patterns.

How do I choose eps?

Start from the k-distance curve (the knee sits between the 80th and 90th percentile of 5th-nearest-neighbour distances), then sweep a range and report the plateau where the cluster count is stable.

What does min_samples mean?

The minimum number of points within eps for a point to be a cluster core. Set it from what you would call a cluster in the domain; 2 Γ— dimensions is a common starting point.

Why does DBSCAN label so many points as noise?

Because they are not in dense regions. That is usually correct and worth reporting. A very high noise fraction means the data is mostly unclustered.

Why do my clusters change every time I adjust eps?

If there is no stable plateau, your clusters have genuinely different densities and no single threshold fits. Use HDBSCAN.

Can I run DBSCAN on latitude and longitude?

Only with metric="haversine", radians, and (lat, lon) order. Euclidean distance on degrees gives a latitude-dependent elliptical neighbourhood.

Does finding clusters mean the points are clustered?

No. DBSCAN finds nine "clusters" in 700 uniformly random points. Use a statistical test if you need to claim significance.