Cluster Points with DBSCAN and K-Means in PyQGIS

A nearest neighbour test can tell you that points are clustered; it cannot tell you where the clusters are or which points belong to each. Clustering algorithms answer those questions by assigning every point a group id. QGIS ships two of them in Processing, and they embody opposite ideas: DBSCAN finds dense areas of any shape and leaves sparse points unassigned, while k-means divides all points into a fixed number of compact groups. Choosing between them is the most important decision in the whole exercise.

This recipe belongs to Spatial Statistics & Pattern Analysis. It runs both algorithms, explains how to choose their parameters from the data, turns the resulting ids into cluster outlines and centres, and sets out which method suits which problem.

Density versus partitionLeft: DBSCAN on the same points finds two dense groups of irregular shape and labels scattered points as noise, which belong to no cluster. Right: k-means with k equals three divides every point into one of three compact groups around centres, including the scattered points, and cannot follow the curved shape of a group.DBSCAN finds density, k-means divides spacetwo groups + noise (hollow)k = 3, every point assigned

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series.
  • A point layer in a projected CRS; DBSCAN's distance parameter is in layer units.
  • Duplicates removed. Repeated points make dense groups out of nothing.

Run DBSCAN

DBSCAN needs two parameters: a search distance (epsilon) and a minimum number of points. A point with at least the minimum number of neighbours within epsilon is a core point; groups grow from core points through their neighbours; points that are neither core nor reachable from a core point are noise.

import processing
from qgis.core import QgsProject

incidents = QgsProject.instance().mapLayersByName("incidents")[0]
db = processing.run("native:dbscanclustering", {
    "INPUT": incidents,
    "EPS": 150,                # metres
    "MIN_SIZE": 8,
    "DBSCAN*": False,
    "FIELD_NAME": "CLUSTER_ID",
    "SIZE_FIELD_NAME": "CLUSTER_SIZE",
    "OUTPUT": "memory:incidents_dbscan",
})
out = db["OUTPUT"]
print(db["NUM_CLUSTERS"], "groups found")
noise = sum(1 for f in out.getFeatures() if f["CLUSTER_ID"] is None
            or (hasattr(f["CLUSTER_ID"], "isNull") and f["CLUSTER_ID"].isNull()))
print(noise, "noise points of", out.featureCount())
QgsProject.instance().addMapLayer(out)

Breakdown: Every input point appears in the output with two new fields: the group id and the size of its group. Noise points have a NULL id — not 0 or −1 — so test for NULL explicitly, as covered in handling NULL values. NUM_CLUSTERS gives the count without iterating. The DBSCAN* option changes how border points are treated: with it, only core points are assigned and border points become noise, which produces tighter groups and more noise. Styling the output with a categorized renderer on CLUSTER_ID shows the result immediately.

Choose epsilon from the data

DBSCAN's results depend heavily on epsilon. Too small and almost everything is noise; too large and separate groups merge into one. The k-distance plot is the standard way to choose it: sort every point's distance to its kth nearest neighbour, where k is the minimum size, and look for the elbow where distances start rising steeply.

The k-distance elbowDistances from each point to its eighth nearest neighbour, sorted ascending. The curve stays low and flat for points inside dense groups, then bends sharply upward for points in sparse areas. The distance at the bend, here about 150 metres, is a good epsilon: points left of it fall in dense groups and points right of it become noise.Pick epsilon at the bendelbow ≈ 150 m → EPSpoints sorted by distance to 8th neighbourm

from qgis.core import QgsSpatialIndex

K = 8
index = QgsSpatialIndex(incidents.getFeatures(),
                        flags=QgsSpatialIndex.FlagStoreFeatureGeometries)
pts = {f.id(): f.geometry() for f in incidents.getFeatures()}
kdist = []
for fid, g in pts.items():
    neighbours = [n for n in index.nearestNeighbor(g.asPoint(), K + 1) if n != fid][:K]
    kdist.append(g.distance(pts[neighbours[-1]]))
kdist.sort()
for q in (0.5, 0.75, 0.9, 0.95):
    print(f"{int(q * 100)}th percentile: {kdist[int(len(kdist) * q)]:.0f} m")

Breakdown: Asking for K + 1 neighbours and dropping the point itself gives the true k nearest. Printing percentiles is a quick numerical stand-in for the plot; a jump between the 75th and 90th percentile is the elbow. Plotting kdist with matplotlib, as in plotting layer data with matplotlib, makes the bend obvious. Choose the minimum size first from the problem — "a hot spot needs at least eight incidents" — then epsilon from the curve.

Run k-means

K-means needs only the number of groups. It places k centres, assigns each point to the nearest, moves each centre to the mean of its points, and repeats until assignments stop changing.

km = processing.run("native:kmeansclustering", {
    "INPUT": incidents,
    "CLUSTERS": 6,
    "FIELD_NAME": "CLUSTER_ID",
    "SIZE_FIELD_NAME": "CLUSTER_SIZE",
    "OUTPUT": "memory:incidents_kmeans",
})["OUTPUT"]

from collections import Counter
sizes = Counter(f["CLUSTER_ID"] for f in km.getFeatures())
print(sorted(sizes.values(), reverse=True))

Breakdown: Every point is assigned; there is no noise. Group sizes are often very uneven, because k-means minimises distances to centres rather than balancing counts. The algorithm starts from initial centres chosen by the k-means++ method, so results are usually stable, but different runs can differ slightly on ambiguous data. K-means suits partitioning tasks — dividing customers among six depots, splitting inspection sites into six daily routes — rather than discovering natural groups.

Choose k

There is no single right k, but the within-group spread falls as k rises, and the point where adding groups stops helping much is a reasonable choice.

import math

def within_spread(layer):
    groups = {}
    for f in layer.getFeatures():
        groups.setdefault(f["CLUSTER_ID"], []).append(f.geometry().asPoint())
    total = 0.0
    for pts in groups.values():
        cx = sum(p.x() for p in pts) / len(pts)
        cy = sum(p.y() for p in pts) / len(pts)
        total += sum((p.x() - cx) ** 2 + (p.y() - cy) ** 2 for p in pts)
    return math.sqrt(total / layer.featureCount())

for k in range(2, 11):
    res = processing.run("native:kmeansclustering", {
        "INPUT": incidents, "CLUSTERS": k, "OUTPUT": "memory:"})["OUTPUT"]
    print(k, round(within_spread(res)))

Breakdown: The root-mean-square distance of points to their group centre falls steeply for the first few values of k and then flattens; the bend is the "elbow" for k-means, analogous to the k-distance plot for DBSCAN. Often the problem fixes k anyway — the number of depots, vans or districts — and the elbow simply shows whether that number fits the data's natural structure.

Summarise each group

Point ids are hard to communicate. Convex hulls and centres, with counts, turn groups into polygons and points you can map, label and report.

From labelled points to outlines and centresPoints carrying a group id are aggregated per group: a convex hull outlines each group, its area and point count are calculated, and a centre point marks its mean location. Noise points are excluded. The resulting polygon and point layers can be styled, labelled and exported for reports.Points with ids → hulls and centrespointsCLUSTER_IDnoise = NULLnative:convexhullper CLUSTER_IDnative:meancoordinatesper CLUSTER_IDmap-readyoutline · countarea · centre

dense = out.materialize(QgsFeatureRequest().setFilterExpression('"CLUSTER_ID" IS NOT NULL'))
hulls = processing.run("native:minimumboundinggeometry", {
    "INPUT": dense, "FIELD": "CLUSTER_ID", "TYPE": 3,     # 3 = convex hull
    "OUTPUT": "memory:dbscan_hulls"})["OUTPUT"]
centres = processing.run("native:meancoordinates", {
    "INPUT": dense, "UID": "CLUSTER_ID", "WEIGHT": None,
    "OUTPUT": "memory:dbscan_centres"})["OUTPUT"]
for f in hulls.getFeatures():
    print(f["CLUSTER_ID"], f["count"], round(f["area"]))
QgsProject.instance().addMapLayers([hulls, centres])

Breakdown: native:minimumboundinggeometry with a group field builds one geometry per group and adds count, area and perimeter fields; type 3 is the convex hull, other values give envelopes, oriented rectangles and minimum circles. Excluding noise first prevents a meaningless hull around scattered points. Mean coordinates per group are the centres; adding a weight field makes them weighted centres, as covered in mean centre and standard distance. Add from qgis.core import QgsFeatureRequest when running this alone.

Go further with scikit-learn

The two Processing algorithms cover the common cases. When densities vary across the map, or when you want clustering that respects attributes as well as location, the scikit-learn library — installable into QGIS's Python as described in installing Python packages into QGIS — offers HDBSCAN, which needs no epsilon at all.

import numpy as np
from sklearn.cluster import HDBSCAN
from qgis.core import edit, QgsField
from qgis.PyQt.QtCore import QVariant

feats = list(incidents.getFeatures())
xy = np.array([[f.geometry().asPoint().x(), f.geometry().asPoint().y()] for f in feats])
labels = HDBSCAN(min_cluster_size=8).fit_predict(xy)     # -1 = noise

layer = incidents.materialize(QgsFeatureRequest())
layer.dataProvider().addAttributes([QgsField("hdb_id", QVariant.Int)])
layer.updateFields()
idx = layer.fields().indexOf("hdb_id")
changes = {f.id(): {idx: (int(l) if l >= 0 else None)}
           for f, l in zip(layer.getFeatures(), labels)}
layer.dataProvider().changeAttributeValues(changes)
QgsProject.instance().addMapLayer(layer)
print(len(set(labels)) - (1 if -1 in labels else 0), "groups from HDBSCAN")

Breakdown: HDBSCAN builds a hierarchy of DBSCAN results over all epsilon values and keeps the most stable groups, so dense cores and sparse outskirts can both be found in one run; min_cluster_size plays the role of the minimum size. scikit-learn marks noise as −1, which the code converts to NULL to match the convention of the QGIS algorithms. Materialising the layer gives a memory copy with the same feature order, so labels can be written back by position. The same pattern — coordinates out as a NumPy array, labels back as a field — works for any scikit-learn estimator, which the Python data stack guide covers in general.

Which method for which question

DBSCAN suits discovery: where are the dense areas, how many are there, which points are outliers? It finds groups of any shape — along a street, around a station — and it is honest about points that belong nowhere. Its weakness is varying density: one epsilon cannot fit a dense city centre and sparse suburbs at once, which is where HDBSCAN from the hdbscan or scikit-learn Python packages helps.

K-means suits allocation: divide these points into exactly k compact groups. It is fast and always assigns everything, but it assumes roughly round groups of similar spread and it will happily split a natural group in two or merge two into one to reach the requested k.

Neither method tests significance. To ask whether a dense area is denser than chance would produce, combine clustering with a hot spot analysis.

QGIS version compatibility

native:dbscanclustering and native:kmeansclustering have been available since QGIS 3.8 and are unchanged on 3.34 LTR, 3.40 LTR and QGIS 4. The NUM_CLUSTERS output exists for both. native:minimumboundinggeometry and native:meancoordinates are available on all current releases.

Troubleshooting

  • Everything is noise. Epsilon is too small or the layer is in degrees; check the CRS and the k-distance percentiles.
  • One giant group. Epsilon is too large; lower it towards the elbow.
  • K-means groups split obvious clusters. k is larger than the natural number of groups, or groups are elongated; try DBSCAN.
  • Group ids change between runs. Ids are arbitrary labels; compare groups by membership, not by number.

Conclusion

Use DBSCAN to discover dense groups and outliers, with minimum size chosen from the problem and epsilon from the k-distance elbow; use k-means to divide points into a fixed number of compact groups, checking the spread curve for a sensible k; summarise groups as hulls and centres; and test significance separately.

Frequently Asked Questions

Can I cluster polygons? Cluster their centroids or points on surface; the algorithms only accept points.

Does DBSCAN consider attributes? No, only location. For attribute-aware clustering, use scikit-learn with location and attributes as features.

How do I cluster by travel distance instead of straight-line distance? Compute a network distance matrix as in building an origin–destination matrix and run DBSCAN on it in scikit-learn with a precomputed metric.

Is the point cluster renderer the same thing? No. The point cluster renderer groups overlapping symbols for display at the current scale; it does not analyse anything.