How to Cluster Points by Location with DBSCAN in Python

On this page

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

coords = np.column_stack([incidents.geometry.x, incidents.geometry.y])      # metres, EPSG:27700
labels = DBSCAN().fit_predict(coords)
print(len(set(labels)))

Fifty-three labels — 52 clusters and the noise label — from Leeds' 903 shoplifting offences in June 2026 (see The data used on this page). eps defaults to 0.5, which in metres is half a metre, so this should have labelled everything noise. It found clusters anyway because 903 offences sit on only 226 distinct points: every location with five or more offences became a "cluster" of identical coordinates, with an area of zero.

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

Quick answer

Project, collapse identical coordinates to one row per location, then cluster the locations:

import numpy as np
from sklearn.cluster import DBSCAN

pts = incidents.to_crs("EPSG:27700")                        # metres — mandatory
sites = (pts.groupby([pts.geometry.x, pts.geometry.y]).size()
            .rename("offences").rename_axis(["x", "y"]).reset_index())   # one row per location
xy = sites[["x", "y"]].to_numpy()

sites["cluster"] = DBSCAN(eps=700, min_samples=3).fit_predict(xy)

inside = sites[sites["cluster"] != -1]
summary = inside.groupby("cluster")["offences"].agg(locations="size", offences="sum")
print(f"{len(summary)} clusters of {len(inside)} locations holding {inside['offences'].sum()} "
      f"of {len(pts)} offences; {(sites['cluster'] == -1).sum()} locations are noise")
print(summary.sort_values("offences", ascending=False).head(6).to_string())

24 clusters of 165 locations holding 688 of 903 offences; 61 locations are noise
         locations  offences
cluster                     
10              71       363
7                7        38
16               5        34
2                4        25
6                7        21
13               4        21

Three parameters, and two of them are hard:

Parameter	Meaning	How to set it
one row per location	what a "point" is	collapse duplicates first, carry the count as a column
`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.

Three roles. Border points join a cluster without being dense themselves, which is why DBSCAN finds elongated shapes.

Project, collapse duplicates, eps from the k-distance curve

1. Project first, and assert it

points = incidents
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. Collapse identical coordinates

counts = pts.geometry.value_counts()
print(f"{len(pts)} records at {len(counts)} distinct locations; "
      f"{(counts == 1).sum()} locations have one record, the busiest has {counts.max()}")

903 records at 226 distinct locations; 102 locations have one record, the busiest has 38

Real point data is rarely one point per event. police.uk moves every crime to an anonymised "on or near" point — a street segment, a car park, a shop — so a supermarket's offences share one coordinate. Geocoded addresses do the same at building or postcode centroids. DBSCAN counts points, so five records at one spot satisfy min_samples=5 with a "cluster" of zero area: the default run in the problem statement found its 52 clusters that way. Cluster the distinct locations and carry the count along as a column, so "how many offences" stays an answer rather than a parameter.

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

Sorted 5th-nearest-neighbour distances rise slowly through the clustered locations then sharply for the isolated ones. The knee — where the rise steepens — is a reasonable eps. On the raw records this curve starts at 0 m, because most records have a duplicate; on the distinct locations it measures real spacing.

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

4. Sweep and find the plateau

for eps in (150, 250, 400, 600, 800, 1200, 2000):
    labels = DBSCAN(eps=eps, min_samples=3).fit_predict(xy)
    held = sites["offences"][labels != -1].sum() / sites["offences"].sum()
    print(f"eps={eps:5}  clusters {len(set(labels) - {-1}):3}  "
          f"noise locations {(labels == -1).sum():4}  offences in clusters {held:4.0%}")

eps=  150  clusters   9  noise locations  180  offences in clusters  28%
eps=  250  clusters  14  noise locations  146  offences in clusters  44%
eps=  400  clusters  19  noise locations  124  offences in clusters  52%
eps=  600  clusters  22  noise locations   89  offences in clusters  67%
eps=  800  clusters  21  noise locations   51  offences in clusters  81%
eps= 1200  clusters  11  noise locations   17  offences in clusters  94%
eps= 2000  clusters   8  noise locations   10  offences in clusters  98%

Pick from inside the plateau and report its width. At the top of the range everything merges into a few large groups and almost nothing is noise — which is what "no structure found" looks like.

5. Set `min_samples` deliberately

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

min_samples=  2  clusters  35  noise locations 17%
min_samples=  3  clusters  24  noise locations 27%
min_samples=  5  clusters   5  noise locations 61%
min_samples= 10  clusters   1  noise locations 77%

On distinct locations min_samples means "this many different places within eps", which is a question you can answer from the domain: three shops with offences within 700 m of each other is a shopping area; two is a coincidence.

6. Join the labels back to the records and check the totals

pts = pts.merge(sites[["x", "y", "cluster"]], left_on=[pts.geometry.x, pts.geometry.y],
                right_on=["x", "y"], how="left").drop(columns=["x", "y"])
clustered = pts[pts["cluster"] != -1]
print(f"{len(clustered)} of {len(pts)} offences in {clustered['cluster'].nunique()} clusters "
      f"({len(clustered) / len(pts):.0%})")
assert len(pts) == len(incidents), "row alignment broken"

688 of 903 offences in 24 clusters (76%)

The merge on coordinates restores one row per record. Merging on the rounded coordinate values rather than relying on row order means nothing can shift: fit_predict returns labels in input order, which is only safe while nobody reorders or filters the frame in between.

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

A DBSCAN function that reports its parameters

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. rate is crimes per 1,000 residents. The clustering examples use the shoplifting offences, which concentrate where the shops are.

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")

incidents = incidents[incidents["category"] == "shoplifting"].reset_index(drop=True)
print(f"{len(incidents):,} shoplifting offences at {incidents.geometry.nunique():,} 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 — 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=3, crs=None):
    """DBSCAN over distinct point locations; returns one row per location with a count and a cluster."""
    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]
    if len(valid) < len(working):
        print(f"dropped {len(working) - len(valid)} null or empty geometries")

    sites = (valid.groupby([valid.geometry.x, valid.geometry.y]).size()
                  .rename("records").rename_axis(["x", "y"]).reset_index())
    xy = sites[["x", "y"]].to_numpy()
    sites["cluster"] = DBSCAN(eps=eps, min_samples=min_samples).fit_predict(xy)

    sizes = sites[sites["cluster"] != -1].groupby("cluster")["records"].sum()
    print(f"eps={eps} min_samples={min_samples}: {len(valid)} records at {len(sites)} locations, "
          f"{len(sizes)} clusters holding {sizes.sum()} records ({sizes.sum() / len(valid):.0%}), "
          f"cluster sizes {sizes.min() if len(sizes) else 0}–{sizes.max() if len(sizes) else 0}")
    return gpd.GeoDataFrame(sites, geometry=gpd.points_from_xy(sites.x, sites.y), crs=working.crs)


def eps_candidates(points, k=5):
    xy = np.unique(np.column_stack([points.geometry.x, points.geometry.y]), axis=0)   # distinct locations
    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=700, min_samples=3)

{'p50': 940.2, 'p75': 1240.8, 'p80': 1410.1, 'p90': 2425.2, 'p95': 3052.7}
eps=700 min_samples=3: 903 records at 226 locations, 24 clusters holding 688 records (76%), cluster sizes 3–363

Collapsing duplicates before clustering is what makes min_samples mean places rather than records. Dropping null geometries before extracting coordinates keeps everything aligned.

Example 2 — describing each cluster

from shapely.geometry import MultiPoint


def describe_clusters(clustered):
    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
        records.append({
            "cluster": int(cid),
            "locations": len(group),
            "records": int(group["records"].sum()),
            "area_km2": round(area_km2, 2),
            "x": round(centre.x), "y": round(centre.y),
            "geometry": hull,
        })
    hulls = gpd.GeoDataFrame(records, geometry="geometry", crs=clustered.crs)
    return hulls.sort_values("records", ascending=False)


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

 cluster  locations  records  area_km2      x      y
      10         71      363      9.22 429992 433687
       7          7       38      0.38 426246 427373
      16          5       34      0.23 430307 431011
       2          4       25      0.08 421662 434444
       6          7       21      0.29 426021 435415
      13          4       21      0.08 428920 437034
       4          3       19      0.13 423438 437743
       1          4       18      0.17 420688 441205

The convex hull is a reasonable summary for compact clusters and a poor one for linear ones. For locations 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=3):
    wgs = points.to_crs("EPSG:4326")
    lonlat = np.unique(np.column_stack([wgs.geometry.x, wgs.geometry.y]), axis=0)   # distinct locations
    radians = np.radians(lonlat[:, ::-1])                                         # (lat, lon)
    return DBSCAN(
        eps=eps_metres / EARTH_RADIUS_M,
        min_samples=min_samples,
        metric="haversine",
        algorithm="ball_tree",           # haversine requires ball_tree
    ).fit_predict(radians)


for metres in (400, 700, 1200):
    labels = cluster_geographic(incidents, metres)
    print(f"eps={metres:5} m -> {len(set(labels) - {-1})} clusters, {(labels == -1).sum()} noise locations")

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

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. Within one city the answers match the projected run closely; reserve this for genuinely global data.

Defaults, min_samples and border points

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` matters so much on sparse data

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 on dense, well-separated clusters raising min_samples mostly trims the periphery and the clusters survive.

The Leeds shoplifting locations are not like that. Outside the city centre, shopping areas have three or four locations with offences, so the core threshold is the difference between a cluster and noise: 35 clusters at min_samples=2, 24 at 3, 5 at 5, and only the city centre at 10. Here min_samples is not a tuning knob but a definition — "how many shops make a shopping area" — and it has to be chosen and stated, not swept for stability.

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.

Given a long cluster, k-means will cut it in two rather than follow it.

Why to keep the noise

The instinct is to tune until nothing is noise. Resist it — at eps=2000 only 10 locations are noise, and the city centre has swallowed its suburbs, which conveys less than the 24 clusters at 700 m.

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

"76% of shoplifting offences fell within 24 concentrations of three or more locations; the remaining 24% were dispersed."

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

Noise before groupby, null geometry and HDBSCAN

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.