How to Cluster Points by Location with DBSCAN in Python

Problem statement

You have several thousand locations — incidents, sightings, deliveries, sensor readings — and you need to identify the concentrations. Not a heatmap: actual groups, with membership you can join back to the original records, count, and hand to someone.

The two-line version runs and gives you nothing useful:

from sklearn.cluster import DBSCAN

labels = DBSCAN().fit_predict(coords)
print(len(set(labels)))

One cluster, because eps defaults to 0.5 and your coordinates are in metres — so every point is its own island, or everything merges, depending which way the mismatch runs.

DBSCAN is the right algorithm for this job. It needs three things done correctly first: a projected CRS, a defensible eps, and a check that the answer is stable.

Quick answer

import numpy as np
from sklearn.cluster import DBSCAN

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

pts["cluster"] = DBSCAN(eps=200, min_samples=5).fit_predict(xy)

summary = pts[pts["cluster"] != -1].groupby("cluster").size().sort_values(ascending=False)
print(f"{len(summary)} clusters, {(pts['cluster'] == -1).sum()} noise points")
print(summary.to_string())

4 clusters, 128 noise points
cluster
1    211
3    157
0    123
2     81

Two parameters, and only one of them is hard:

Parameter	Meaning	How to set it
`eps`	neighbourhood radius, in coordinate units	k-distance curve, then a stability sweep
`min_samples`	points within `eps` to be a cluster core	the smallest group you would call a cluster

-1 is the noise label. It is not a cluster and it must be filtered before any groupby.

Core points with enough neighbours inside eps, border points adjacent to a core, and noise points with neither. — Three roles. Border points join a cluster without being dense themselves, which is why DBSCAN finds elongated shapes.

Step-by-step solution

1. Project first, and assert it

if points.crs is None or points.crs.is_geographic:
    raise ValueError(f"eps is in coordinate units; {points.crs} gives degrees")
pts = points.to_crs("EPSG:27700")

This single check prevents the commonest DBSCAN failure. eps=200 on degrees is 200 degrees — more than half the planet — so everything becomes one cluster. eps=0.002 on metres is 2 mm, so everything becomes noise. Both run without complaint.

2. Get `eps` in the right range from the k-distance curve

from sklearn.neighbors import NearestNeighbors

k = 5
distances, _ = NearestNeighbors(n_neighbors=k).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

Sorted 5th-nearest-neighbour distances rise slowly through the clustered points then sharply for the isolated ones. The knee — where the rise steepens — is a reasonable eps. Here it is somewhere between 181 and 361 m.

That is a range, not a value, and the knee is frequently ambiguous. Use it to bracket the sweep.

3. Sweep and find the plateau

for eps in (30, 60, 100, 150, 250, 400, 800):
    labels = DBSCAN(eps=eps, min_samples=5).fit_predict(xy)
    print(f"eps={eps:4}  clusters {len(set(labels) - {-1}):3}  "
          f"noise {(labels == -1).sum():4} ({(labels == -1).mean():4.0%})")

eps=  30  clusters  14  noise  400 ( 57%)
eps=  60  clusters   6  noise  216 ( 31%)
eps= 100  clusters   5  noise  159 ( 23%)
eps= 150  clusters   4  noise  138 ( 20%)
eps= 250  clusters   4  noise  116 ( 17%)
eps= 400  clusters   5  noise   29 (  4%)
eps= 800  clusters   1  noise    0 (  0%)

Four clusters from 150 to 250. Below that the loose clusters fragment; above 400 they start merging, and at 800 everything is one group with no noise at all — which is what "no structure found" looks like.

Pick from inside the plateau and report its width.

4. Set `min_samples` deliberately

for ms in (3, 5, 10, 20):
    labels = DBSCAN(eps=200, min_samples=ms).fit_predict(xy)
    print(f"min_samples={ms:3}  clusters {len(set(labels) - {-1})}  "
          f"noise {(labels == -1).mean():.0%}")

min_samples=  3  clusters 14  noise 12%
min_samples=  5  clusters  4  noise 18%
min_samples= 10  clusters  4  noise 19%
min_samples= 20  clusters  4  noise 20%

From 5 upwards the cluster count is rock stable and the noise fraction barely moves. At 3 it explodes to fourteen — three points within 200 m of each other is a common accident, so the algorithm starts finding groups in the background scatter.

Choose it from the domain: if three events is not a cluster to you, min_samples=3 will produce ten groups you have to explain away.

5. Join the labels back and check the totals

pts["cluster"] = DBSCAN(eps=200, min_samples=5).fit_predict(xy)

clustered = pts[pts["cluster"] != -1]
print(f"{len(clustered)} of {len(pts)} points in {clustered['cluster'].nunique()} clusters "
      f"({len(clustered) / len(pts):.0%})")
assert len(pts) == len(incidents), "row alignment broken"

572 of 700 points in 4 clusters (82%)

fit_predict returns labels in input order, so assigning them straight to a column is safe provided you have not reordered or filtered the frame between building xy and assigning. If you dropped rows with null geometry, do it before extracting coordinates.

Report the plateau. A single run reports the parameter, not the data.

Code examples

Example 1 — a clustering function that refuses to hide its parameters

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


def cluster_points(points, eps, *, min_samples=5, crs=None):
    """DBSCAN over point geometry, returning the frame with a cluster column."""
    working = points.to_crs(crs) if crs else points
    if working.crs is None or working.crs.is_geographic:
        raise ValueError(f"eps is in coordinate units; {working.crs} gives degrees — "
                         f"pass crs= a projected CRS")

    valid = working[working.geometry.notna() & ~working.geometry.is_empty].copy()
    if len(valid) < len(working):
        print(f"dropped {len(working) - len(valid)} null or empty geometries")

    xy = np.column_stack([valid.geometry.x, valid.geometry.y])
    labels = DBSCAN(eps=eps, min_samples=min_samples).fit_predict(xy)
    valid["cluster"] = labels

    sizes = pd.Series(labels[labels != -1]).value_counts()
    print(f"eps={eps} min_samples={min_samples}: {len(sizes)} clusters, "
          f"{(labels == -1).sum()} noise ({(labels == -1).mean():.0%}), "
          f"sizes {sizes.min() if len(sizes) else 0}–{sizes.max() if len(sizes) else 0}")
    return valid


def eps_candidates(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, 75, 80, 90, 95)}


pts = incidents.to_crs("EPSG:27700")
print(eps_candidates(pts))
clustered = cluster_points(pts, eps=200, min_samples=5)

{'p50': 42.0, 'p75': 137.0, 'p80': 181.0, 'p90': 361.0, 'p95': 473.0}
eps=200 min_samples=5: 4 clusters, 128 noise (18%), sizes 81–211

Dropping null geometries before extracting coordinates is what keeps the label array aligned with the rows. Doing it afterwards is a silent off-by-n that assigns every label to the wrong record.

Example 2 — describing each cluster

from shapely.geometry import MultiPoint


def describe_clusters(clustered, value=None):
    records = []
    for cid, group in clustered[clustered["cluster"] != -1].groupby("cluster"):
        hull = MultiPoint(list(group.geometry)).convex_hull
        area_km2 = hull.area / 1e6 if hull.geom_type == "Polygon" else 0.0
        centre = group.geometry.union_all().centroid
        record = {
            "cluster": int(cid),
            "n": len(group),
            "area_km2": round(area_km2, 3),
            "per_km2": round(len(group) / area_km2) if area_km2 else None,
            "x": round(centre.x), "y": round(centre.y),
            "geometry": hull,
        }
        if value:
            record[f"mean_{value}"] = round(group[value].mean(), 1)
        records.append(record)

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


hulls = describe_clusters(clustered)
print(hulls.drop(columns="geometry").to_string(index=False))
hulls.to_file("clusters.gpkg", layer="hulls", driver="GPKG")

 cluster   n  area_km2  per_km2      x      y
       1 211     0.413      511 386998 397199
       3 157     0.821      191 387496 399201
       0 123     0.125      983 385003 395002
       2  81     0.039     2067 385501 398000

Four clusters spanning 191 to 2,067 points per km² — an elevenfold density range that one eps had to accommodate. If that spread were wider, the sweep would show no plateau and HDBSCAN would be the right tool.

The convex hull is a reasonable summary for compact clusters and a poor one for linear ones. For points along a road, buffer and dissolve instead:

outline = group.buffer(eps).union_all()          # follows the shape of the cluster

Example 3 — clustering points that cannot be projected

For data spanning several UTM zones, use great-circle distance directly:

EARTH_RADIUS_M = 6_371_000.0


def cluster_geographic(points, eps_metres, min_samples=5):
    wgs = points.to_crs("EPSG:4326")
    valid = wgs[wgs.geometry.notna()].copy()
    radians = np.radians(np.column_stack([valid.geometry.y, valid.geometry.x]))  # lat, lon

    valid["cluster"] = DBSCAN(
        eps=eps_metres / EARTH_RADIUS_M,
        min_samples=min_samples,
        metric="haversine",
        algorithm="ball_tree",           # haversine requires ball_tree
    ).fit_predict(radians)
    return valid


for metres in (50, 100, 200, 400):
    result = cluster_geographic(incidents, metres)
    n_clusters = result.loc[result["cluster"] != -1, "cluster"].nunique()
    print(f"eps={metres:4} m -> {n_clusters} clusters, "
          f"{(result['cluster'] == -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

Three things must all be right here: radians, not degrees; (lat, lon) order, not (x, y); and algorithm="ball_tree", which is the only one supporting haversine. Get any of them wrong and it either raises or — worse — silently measures the wrong thing.

For a single city, projecting is simpler and faster. Reserve this for genuinely global data.

Explanation

Why the default parameters never work

DBSCAN() defaults to eps=0.5, min_samples=5. Half a coordinate unit.

On EPSG:27700 (metres) that is 50 cm: every point is isolated, everything is noise.
On EPSG:4326 (degrees) that is about 55 km: everything within a city merges into one cluster.

Neither raises. Both produce a labels array of the right length that you can happily join back and map. This is why the CRS assertion belongs at the top of the function rather than in a comment.

Why `min_samples` changes the noise fraction more than the cluster count

min_samples decides how many neighbours a point needs to be a core point. Non-core points adjacent to a core point still join the cluster as border points — so raising min_samples shrinks the core skeleton, but the clusters it finds are largely the same ones.

What changes is the periphery. In the sweep above the cluster count was 4 from min_samples=5 to 20, and the noise fraction moved only from 18% to 20%. The four clusters were dense enough that shrinking the core skeleton barely changed which points they reached.

The exception is at the bottom of the range. At min_samples=3 the count jumped to 14, because a threshold that low is met by chance triples in the background scatter. So min_samples is forgiving above a sensible floor and sharply wrong below it, while eps needs the sweep across its whole range.

Why border points matter

DBSCAN's three roles — core, border, noise — are what let it find non-convex clusters. A chain of core points, each within eps of the next, forms a single cluster of any shape at all: a line of shops, a river's floodplain, a road corridor.

K-means cannot represent those, because minimising within-cluster variance prefers round blobs. This is the main reason DBSCAN is the default for point locations on a map.

One consequence: a border point equidistant from two clusters is assigned to whichever core reached it first, which depends on the order of the input. It is a small non-determinism and it affects only points on the boundary between adjacent clusters.

A chain of overlapping eps neighbourhoods forming one elongated cluster that no circular method could represent. — Connectivity, not compactness. This is what DBSCAN can express and K-means cannot.

Why to keep the noise

The instinct is to tune until nothing is noise. Resist it — eps=800 in the sweep gives exactly that, and the result is one cluster containing every point, which conveys nothing.

Noise is the set of points not participating in any concentration, and its size is a result:

"82% of incidents fall within four concentrations; the remaining 18% are dispersed."

That sentence is more informative than any cluster count on its own, and it comes free from the -1 label.

Edge cases or notes

Filter -1 before every groupby. Noise is not cluster number −1; treating it as a cluster produces a spurious largest group.
Cluster ids are arbitrary and unstable. They change between runs with different parameters. Match clusters across runs by geometry, never by id.
Drop null geometries before extracting coordinates, or the label array misaligns with the rows.
Duplicate coordinates create fake clusters. A hundred records geocoded to one centroid satisfy any min_samples — see geocoding returns wrong coordinates.
eps should not exceed roughly a tenth of the study extent, or clusters merge simply because everything is within reach of everything.
algorithm="ball_tree" or "kd_tree" keeps DBSCAN near O(n log n). The brute-force fallback is O(n²) and will not finish on large datasets.
HDBSCAN lives in sklearn.cluster from scikit-learn 1.3 — use it when the sweep shows no plateau.
Convex hulls overstate linear clusters. Buffer by eps and dissolve for a shape that follows the points.

Internal links

Spatial clustering explained: DBSCAN, K-means and what they assume — why DBSCAN, and what it assumes
DBSCAN returns one giant cluster or labels everything noise — diagnosing a run that went wrong
Choose a projected CRS for your area — so eps means metres
Kernel density explained — the continuous view of the same pattern
How to find hotspots with Getis-Ord Gi* in Python — testing significance rather than describing
How to create buffers in GeoPandas for spatial analysis — better outlines for linear clusters
How to build a distance matrix between two layers in Python — the distance machinery underneath
Geocoding returns wrong or missing coordinates — where fake clusters come from

FAQ

Why did DBSCAN return one cluster containing everything?

eps is too large for your coordinate units — usually degrees being treated as metres, or the default eps=0.5. Project to metres and sweep a sensible range.

Why is everything labelled noise?

eps is too small, or min_samples too high. Check the k-distance curve: if your eps is below the 50th percentile of 5th-nearest-neighbour distances, almost nothing will be dense enough.

What value should `eps` be?

Whatever the k-distance knee suggests, confirmed by a sweep. Report the range over which the cluster count is stable rather than one number.

How do I attach cluster labels back to my GeoDataFrame?

fit_predict returns labels in input order, so assign directly — but drop null geometries before building the coordinate array, or the alignment breaks.

Can I cluster latitude and longitude directly?

Only with metric="haversine", eps in radians, (lat, lon) order and algorithm="ball_tree". Projecting is simpler for anything smaller than a continent.

What does the `-1` label mean?

Noise — points in no cluster. Filter it before aggregating, and report its share as a result rather than tuning it to zero.

My clusters have very different densities. What should I use?

HDBSCAN, in sklearn.cluster. DBSCAN applies one density threshold everywhere, so no single eps fits clusters that differ substantially.

How to Cluster Points by Location with DBSCAN in Python #

Problem statement #

Quick answer #

Step-by-step solution #

1. Project first, and assert it #

2. Get eps in the right range from the k-distance curve #

3. Sweep and find the plateau #

4. Set min_samples deliberately #

5. Join the labels back and check the totals #

Code examples #

Example 1 — a clustering function that refuses to hide its parameters #

Example 2 — describing each cluster #

Example 3 — clustering points that cannot be projected #

Explanation #

Why the default parameters never work #

Why min_samples changes the noise fraction more than the cluster count #

Why border points matter #

Why to keep the noise #

Edge cases or notes #

Internal links #

FAQ #

Why did DBSCAN return one cluster containing everything? #

Why is everything labelled noise? #

What value should eps be? #

How do I attach cluster labels back to my GeoDataFrame? #

Can I cluster latitude and longitude directly? #

What does the -1 label mean? #

My clusters have very different densities. What should I use? #