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
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 radiusmin_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.
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.
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 withalgorithm="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. -1is the noise label, and it will silently become a cluster if yougroupby("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
algorithmiskd_treeorball_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.clusterfrom 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
epsmeans metres - How to create buffers in GeoPandas for spatial analysis β a better outline for linear clusters
- Nearest-neighbour joins explained β the distance machinery underneath
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.