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

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 = sites.to_crs("EPSG:27700")                          # metres; one row per location
xy = np.column_stack([pts.geometry.x, pts.geometry.y])

labels = DBSCAN(eps=700, min_samples=3).fit_predict(xy)   # eps in metres
pts["cluster"] = labels
print(f"{len(set(labels) - {-1})} clusters, {(labels == -1).sum()} noise points")

24 clusters, 61 noise points

Test for clustering, know the assumptions, then eps

1. Ask whether the points are clustered before clustering them

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

minx, miny, maxx, maxy = pts.total_bounds                 # the same extent and count, no structure
random_points = np.random.default_rng(0).uniform((minx, miny), (maxx, maxy), (len(pts), 2))
labels = DBSCAN(eps=700, min_samples=3).fit_predict(random_points)
print(f"{len(pts)} uniform random points -> {len(set(labels) - {-1})} clusters, {(labels == -1).sum()} noise")

226 uniform random points -> 16 clusters, 168 noise

Sixteen "clusters" in data with no structure whatsoever — the same number of points over the same extent as the Leeds shoplifting locations, placed at random — because uniformly random points are not evenly spaced; they clump by chance. The real locations gave 24 clusters and 61 noise points against the random 16 and 168: more structure, but the random run shows how much "clustering" appears from nothing. 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

Sweep eps on the Leeds shoplifting locations and watch the answer move:

for eps in (150, 250, 400, 600, 800, 1200, 2000):
    labels = DBSCAN(eps=eps, min_samples=3).fit_predict(xy)
    sizes = np.bincount(labels[labels >= 0]) if (labels >= 0).any() else np.array([0])
    print(f"  eps {eps:5}  clusters {len(set(labels) - {-1}):3}  noise {(labels == -1).sum():4}  largest {sizes.max():4}")

  eps   150  clusters   9  noise  180  largest   19
  eps   250  clusters  14  noise  146  largest   35
  eps   400  clusters  19  noise  124  largest   42
  eps   600  clusters  22  noise   89  largest   64
  eps   800  clusters  21  noise   51  largest   80
  eps  1200  clusters  11  noise   17  largest  160
  eps  2000  clusters   8  noise   10  largest  179

At eps=150 most locations are noise and the clusters are fragments. From 1,200 m the city centre swallows its neighbours — the largest cluster jumps from 80 locations to 160. In between, 22 and 21 clusters at 600 and 800 m is the nearest thing to a stable answer, and it is not a clean one.

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: 940 m
  p80: 1410 m
  p90: 2425 m
  p95: 3053 m

The knee sits between p80 and p90 — 1,410 to 2,425 m — which is above the 600–800 m range where the count held. On sparse, uneven data the curve and the sweep can disagree, which is exactly why both are worth running. 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 \times \text{dimensions}$, so 4 or 5 for 2D data. Larger values produce fewer, denser clusters and more noise.

DBSCAN sweeps and cluster polygons for Leeds shoplifting

The data used on this page

Every example runs on real data: the 488 Lower Layer Super Output Areas (LSOAs) of Leeds from the ONS, their Census 2021 population from Nomis, and the street-level crimes West Yorkshire Police published for June 2026 on data.police.uk. The clustering examples use the shoplifting offences, collapsed to one row per distinct location in sites — police.uk moves every crime to an anonymised point, so many offences share a coordinate, and clustering the raw records would find "clusters" of one spot.

import io
import time

import geopandas as gpd
import numpy as np
import pandas as pd
import requests

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query")
features, offset = [], 0
while True:                                   # the service pages at 1,000 features
    page = requests.get(ONS, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "LSOA21CD,LSOA21NM",
                                     "outSR": 27700, "f": "geojson", "resultOffset": offset}, timeout=120).json()
    features += page["features"]
    if not page.get("properties", {}).get("exceededTransferLimit"):
        break
    offset += len(page["features"])
zones = gpd.GeoDataFrame.from_features(features, crs=27700).set_index("LSOA21CD").sort_index()

population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
                         f"&geography={','.join(zones.index)}&c2021_restype_3=0&measures=20100"
                         "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
zones["population"] = population

west, south, east, north = zones.to_crs(4326).total_bounds
crimes = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):   # ≤10,000 a request
        for attempt in range(5):                        # the API answers bursts with 429
            reply = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                                  data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if reply.status_code != 429:
                break
            time.sleep(2 ** attempt)
        reply.raise_for_status()
        crimes += [{"id": c["id"], "category": c["category"], "lon": float(c["location"]["longitude"]),
                    "lat": float(c["location"]["latitude"])} for c in reply.json()]
crimes = pd.DataFrame(crimes).drop_duplicates("id")
incidents = gpd.GeoDataFrame(crimes, geometry=gpd.points_from_xy(crimes.lon, crimes.lat), crs=4326).to_crs(27700)
incidents = gpd.sjoin(incidents, zones[["geometry"]], predicate="within").rename(columns={"index_right": "lsoa", "LSOA21CD": "lsoa"})
zones["crimes"] = incidents["lsoa"].value_counts().reindex(zones.index, fill_value=0)
zones["rate"] = 1000 * zones["crimes"] / zones["population"]
print(f"{len(zones)} LSOAs, {zones['population'].sum():,} residents, {len(incidents):,} crimes; "
      f"rate median {zones['rate'].median():.1f}, max {zones['rate'].max():.1f} per 1,000")

shoplifting = incidents[incidents["category"] == "shoplifting"]
xy_all = np.column_stack([shoplifting.geometry.x, shoplifting.geometry.y])
locations, offences = np.unique(xy_all, axis=0, return_counts=True)
sites = gpd.GeoDataFrame({"offences": offences}, crs=incidents.crs,
                         geometry=gpd.points_from_xy(locations[:, 0], locations[:, 1]))
print(f"{len(shoplifting):,} shoplifting offences at {len(sites):,} distinct locations")

488 LSOAs, 811,957 residents, 9,229 crimes; rate median 8.4, max 174.4 per 1,000
903 shoplifting offences at 226 distinct locations

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))
    if len(plateau["eps"]) < 2:
        print("\nno plateau: no two neighbouring eps values give the same cluster count")
    else:
        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(sites))
dbscan_sweep(sites, [150, 250, 400, 600, 800, 1200, 2000], min_samples=3)

{'p50': 940.2, 'p80': 1410.1, 'p90': 2425.2, 'p95': 3052.7}
 eps  clusters  noise noise_pct  largest
 150         9    180       80%       19
 250        14    146       65%       35
 400        19    124       55%       42
 600        22     89       39%       64
 800        21     51       23%       80
1200        11     17        8%      160
2000         8     10        4%      179

no plateau: no two neighbouring eps values give the same cluster count

plateau: 4 clusters for eps 150–250

Report the plateau, not one run — and when the function says there is none, say that. "About 21 clusters for eps between 600 and 800 m, with 23–39% of locations classed as noise; no count repeats exactly" is a defensible sentence. "24 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 (400, 700, 1200):
    labels = dbscan_geographic(sites, metres, min_samples=3)
    print(f"eps={metres:4} m -> {len(set(labels) - {-1})} clusters, "
          f"{(labels == -1).sum()} noise")

eps= 400 m -> 19 clusters, 124 noise
eps= 700 m -> 24 clusters, 61 noise
eps=1200 m -> 11 clusters, 17 noise

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

wgs = sites.to_crs("EPSG:4326")
degrees = np.column_stack([wgs.geometry.y, wgs.geometry.x])
labels = DBSCAN(eps=700 / 111_000, min_samples=3).fit_predict(degrees)   # 700 m of latitude, in degrees
print(f"euclidean on degrees: {len(set(labels) - {-1})} clusters, {(labels == -1).sum()} noise")

euclidean on degrees: 24 clusters, 90 noise

The same cluster count here, but not the same clusters: 90 locations are noise against 61 with the haversine metric. At Leeds' 53.8°N a degree of longitude is 66 km against 111 km for latitude, so a "circular" eps in degrees is 1.69× 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=700, min_samples=3).fit_predict(
    np.column_stack([sites.geometry.x, sites.geometry.y])
)
hulls = cluster_hulls(sites, labels, min_size=3)
print(hulls[["cluster", "n", "area_km2", "per_km2"]].to_string(index=False))

 cluster  n  area_km2  per_km2
      10 71     9.218      8.0
       6  7     0.290     24.0
       7  7     0.377     19.0
       0  5     0.112     45.0
      16  5     0.234     21.0
      23  5     0.422     12.0
      12  5     0.450     11.0
      20  5     0.087     57.0
      22  5     0.273     18.0
       1  4     0.173     23.0
       2  4     0.078     51.0
       5  4     0.070     57.0
      13  4     0.076     53.0
      21  4     0.031    129.0
       4  3     0.134     22.0
       3  3     0.080     38.0
      15  3     0.015    204.0
      14  3     0.121     25.0
       9  3     0.002   1695.0
      11  3     0.009    339.0
       8  3     0.026    114.0
      19  3     0.015    195.0
      17  3     0.009    318.0
      18  3     0.012    244.0

The per_km2 column shows the varying-density problem plainly: from 8 locations per km² in the city-centre cluster — 71 locations spread over 9.2 km² — to 1,695 in a cluster of three locations on one block, yet a single eps had to accept all of them. A convex hull also overstates a linear cluster badly — for clusters strung along a road, use a buffered union or a concave hull instead.

K-means for locations, density and noise

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.

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: "three quarters of shoplifting offences fall within 24 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.8°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.

Haversine, unstable labels and convex hulls

• 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.

Internal links

• How to cluster points by location with DBSCAN in Python — the practical implementation
• DBSCAN returns one giant cluster or labels everything noise — diagnosing the parameter failures
• How to find hotspots with Getis-Ord Gi* in Python — a significance test, not a description
• Kernel density explained — the continuous view of the same question
• The modifiable areal unit problem explained — why binning first changes the answer
• Choose a projected CRS for your area — so eps means metres
• How to create buffers in GeoPandas for spatial analysis — a better outline for linear clusters
• Nearest-neighbour joins explained — the distance machinery underneath
• How to build management zones with clustering in Python — clustering a field into manageable units
• How to classify urban fabric types — clustering buildings by form
• Forest fragmentation explained — patches, edges and core area

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 \times \text{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 found 16 "clusters" in 226 uniformly random points spread over Leeds. Use a statistical test if you need to claim significance.