Spatial Clustering and Grid Binning in DuckDB

Reducing millions of raw points to a manageable set of groups is the first move in almost every density, hot-spot, or map-tiling workload, and DuckDB Spatial offers four structurally different ways to do it — density-based clustering with ST_ClusterDBSCAN, regular square-grid binning from integer cell keys, equal-area hexagon tessellation through the H3 extension, and prefix bucketing with geohash or quadkey strings. This guide, part of Modern Spatial SQL Query Patterns, treats those four as one family of “collapse many points into groups” operators, shows the canonical vectorized SQL for each, and gives the plan-validation and regression tooling that keeps the chosen method fast and correct as the input grows. The decision that actually matters is not how to write each query — every one is a few lines — but which grouping model matches the analytical question, the coordinate reference system, and the accuracy budget you are working under.

Runtime Configuration & Memory Guardrails #

Grouping workloads split cleanly into two cost profiles. Grid, hexagon, and prefix binning derive a scalar key per row and feed a single HASH_GROUP_BY, so they are memory-cheap and stream well. ST_ClusterDBSCAN is different: it is a window aggregate that buffers each partition’s geometry to compute reachability, so its peak working set scales with the largest partition, not the vector size. Configure the session for the heavier of the two before running either.

INSTALL spatial;          -- one-time download of the spatial extension into the local cache
LOAD spatial;             -- per-session load; ST_ClusterDBSCAN and ST_X/ST_Y resolve only after this

-- Match physical cores. DBSCAN's neighbour search and the hash-build for grid keys are both
-- memory-bandwidth bound, so hyperthread oversubscription raises contention without adding throughput.
SET threads = 8;

-- Ceiling for the largest partition DBSCAN buffers, plus the grid hash table. Too low → the
-- window operator spills mid-partition; too high on a shared host → the OS OOM-kills the process.
SET memory_limit = '12GB';

-- Release row-order guarantees so the hash group-by and window partitions parallelize. Grid and
-- hex keys are order-independent; re-impose order with an explicit ORDER BY on the final result.
SET preserve_insertion_order = false;

-- Spill target for an over-budget DBSCAN partition. Point at local NVMe — a network mount turns
-- a spilling window buffer into an I/O cliff that dwarfs the clustering work itself.
SET temp_directory = '/var/lib/duckdb/spill';

Before running a clustering or binning pipeline, confirm these preconditions hold:

• The coordinate reference system matches the method: metric CRS for square grids and DBSCAN eps in metres, WGS84 lat/lng for H3 and geohash.
• Geometry columns are native GEOMETRY (WKB-backed), not text re-parsed per query — materialize at ingestion as covered in the GeoParquet parsing reference.
• memory_limit exceeds the uncompressed footprint of the largest PARTITION BY group that ST_ClusterDBSCAN will hold in flight.
• The h3 community extension is installed separately if hexagon binning is in scope — it is not bundled with spatial (see the H3 hexagon binning walkthrough).

Trade-off: a generous memory_limit keeps a dense DBSCAN partition resident instead of spilling, but on a shared host it enlarges the blast radius — one skewed partition can starve co-tenants. Pair the limit with SET max_temp_directory_size rather than trusting a single global ceiling, and coarsen the partition key if one bucket dominates.

Choosing a Grouping Model #

Every method answers a slightly different question. Density clustering finds arbitrary-shaped concentrations and labels the sparse remainder as noise. Square-grid binning imposes a uniform lattice for the cheapest possible aggregation. Hexagon binning trades a little arithmetic for equal-area cells with uniform neighbour distance. Prefix bucketing groups by a truncated geohash or quadkey string, which is ideal when the bucket boundaries must align with an existing tiling scheme. Pick the model first, then the SQL is mechanical.

Density clustering with ST_ClusterDBSCAN #

When the interesting groups have no fixed shape — retail catchments, disease hot spots, vessel loitering zones — density-based spatial clustering is the right tool. ST_ClusterDBSCAN is a window function that assigns each geometry a numeric cluster id based on two parameters: eps, the neighbourhood radius in the geometry’s own units, and min_points, the count of neighbours needed to seed a dense region. Points that never reach the density threshold receive NULL — they are noise, not a group of one.

-- Density clustering over a metric CRS: eps is in metres because the geometry is projected.
SELECT
    id,
    geom,
    ST_ClusterDBSCAN(geom, eps := 150.0, min_points := 5)
        OVER ()                       -- empty window: one global clustering pass over all rows
        AS cluster_id
FROM incidents;

The OVER () clause runs one clustering pass across the whole relation; partitioning the window (OVER (PARTITION BY city_id)) clusters each group independently and bounds the buffer per partition. Because eps is expressed in coordinate units, the input must be in a metric CRS for a metres-based radius to mean anything — the same unit discipline that governs grid binning. The parameter-tuning detail, noise handling, and partition-buffer behaviour are covered in depth in the ST_ClusterDBSCAN spatial grouping walkthrough.

The two parameters interact: eps sets how close points must be to be reachable, and min_points sets how many mutually reachable points constitute a dense core. Raising min_points makes the algorithm more conservative — more points fall to noise, and only genuinely dense concentrations survive as clusters. Lowering eps fragments a single loose group into several tight ones. Neither has a universal default: they encode a domain judgement about what “close” and “dense” mean for the phenomenon being clustered, so calibrate them against a labelled sample rather than copying a value from an example. Unlike a grid or hexagon, DBSCAN output is not reproducible across CRS changes — reprojecting the input silently rescales eps, so pin the projection alongside the parameters.

Square-grid binning with integer cell keys #

When you want a uniform lattice and the cheapest possible aggregation, snap each point to a fixed-size cell with integer arithmetic on its coordinates. DuckDB ships no built-in grid generator, but floor(ST_X(geom) / size) and floor(ST_Y(geom) / size) produce a fully vectorizable scalar key pair:

WITH binned AS (
    -- Snap each point to a 500-metre grid cell. Requires a metric CRS: floor(x / 500)
    -- only means "500 m" when ST_X returns metres, not degrees.
    SELECT
        floor(ST_X(geom) / 500) AS cell_x,
        floor(ST_Y(geom) / 500) AS cell_y,
        magnitude,
        geom
    FROM events
)
SELECT
    cell_x,
    cell_y,
    COUNT(*)          AS event_count,
    AVG(magnitude)    AS mean_magnitude,
    ST_Collect(geom)  AS members        -- ST_Collect bundles members; it is an aggregate, not a scalar
FROM binned
GROUP BY cell_x, cell_y;

Because the keys are plain integers, this collapses the work from a pairwise $O(N \times M)$ containment comparison to a single $O(N)$ hash aggregation. The reconstruction — turning a cell key back into a drawable square — is a projection over the group, never work inside the GROUP BY. The broader family of key-then-aggregate patterns is catalogued in the vectorized aggregations guide, and the metric-CRS requirement is spelled out in CRS mapping and transformations.

Hexagon binning with H3 #

Square grids have a defect that matters for density surfaces: a cell’s diagonal neighbours are $\sqrt{2}$ times farther than its edge neighbours, so distance-weighted statistics are biased by direction. Hexagons fix this — every neighbour sits at the same centre-to-centre distance, and the tessellation has no gaps or overlaps. Uber’s H3 grid also gives equal-area cells at each resolution. In DuckDB this comes from the separate h3 community extension:

INSTALL h3 FROM community;   -- separate community extension; NOT bundled with spatial
LOAD h3;

-- Convert WGS84 points to H3 cells at resolution 8 (~0.7 km² per cell), then aggregate by cell.
SELECT
    h3_latlng_to_cell(ST_Y(geom), ST_X(geom), 8) AS h3_cell,   -- note: lat first, then lng
    COUNT(*)       AS point_count,
    AVG(magnitude) AS mean_magnitude
FROM events
GROUP BY h3_cell;

The one gotcha worth flagging early is axis order: h3_latlng_to_cell expects latitude then longitude, whereas ST_X returns longitude and ST_Y returns latitude — swapping them silently mislocates every point. The resolution table, boundary reconstruction with h3_cell_to_boundary_wkt, and the full WGS84 requirement live in the H3 hexagon binning reference.

Prefix bucketing with geohash #

When the buckets must align with an existing tiling scheme — a CDN tile pyramid, a partitioned object store, a geohash-keyed cache — group by a truncated geohash. DuckDB Spatial produces a geohash string with ST_GeoHash, and truncating the prefix coarsens the bucket: a 5-character prefix is roughly a 5 km cell, a 6-character prefix roughly 1 km.

-- Bucket by a 6-character geohash prefix (~1.2 km cells). Prefix length is the resolution knob.
SELECT
    substr(ST_GeoHash(geom), 1, 6) AS gh6,
    COUNT(*)      AS point_count
FROM events
GROUP BY gh6;

Prefix bucketing shares the scalar-key advantage of grid binning — the group key is a short string, the aggregation is a single hash pass — but its cells are neither equal-area nor square; they narrow toward the poles like any lat/lng-derived scheme. It earns its place only when downstream systems already speak geohash. The quadkey scheme used by slippy-map tile pyramids behaves the same way: a quadkey is a base-4 path down a quadtree, and truncating it to a shorter length selects a coarser zoom level, so grouping by a fixed-length quadkey prefix buckets points to map tiles at that zoom. Both schemes are hierarchical by string length, which is exactly what makes them cheap to re-aggregate: a coarser bucket is a prefix of a finer one, so a roll-up needs no re-scan of the source geometry.

Execution Plan Validation #

Reading the plan is the only reliable way to confirm a grouping query stayed vectorized. The three key-based methods — grid, hexagon, prefix — should all resolve to a HASH_GROUP_BY over scalar keys feeding from a plain TABLE_SCAN. ST_ClusterDBSCAN instead shows a WINDOW operator; that is expected, not a regression.

EXPLAIN
SELECT cell_x, cell_y, COUNT(*)
FROM binned
GROUP BY cell_x, cell_y;

What to assert in the output:

• Operator type. For key-based binning, a HASH_GROUP_BY over the scalar key is correct. A NESTED_LOOP_JOIN means a stray cross join crept in; a sort-based group means the key is not scalar — verify you grouped on the integer cell coordinates or the H3/geohash key, not on a geometry.
• No row-at-a-time fallback. A ROW_EXECUTION node or a SCALAR_FUNCTION_EVAL of an ST_ function inside the aggregation phase means the query dropped out of SIMD kernels — usually from a scalar Python UDF or an implicit GEOMETRY → VARCHAR cast in the key expression.
• DBSCAN window shape. For ST_ClusterDBSCAN, expect a WINDOW node whose partition matches your PARTITION BY. A single unpartitioned window over a huge relation is the classic DBSCAN memory trap; watch for a non-empty spill indicator on that node.

-- Measured plan with per-operator timing and peak memory.
EXPLAIN ANALYZE
SELECT h3_latlng_to_cell(ST_Y(geom), ST_X(geom), 8) AS h3_cell, COUNT(*)
FROM events
GROUP BY h3_cell;

A sudden jump in operator_timing on the HASH_GROUP_BY, or a spill indicator appearing on the DBSCAN WINDOW, localizes a regression to exactly one node.

Performance Trade-offs #

The four methods trade accuracy, cost, and CRS requirements against each other. Quantify against your own data, but these are the characteristic ranges:

The single highest-leverage rule across all four is keeping the group key scalar. An integer cell pair, an H3 UBIGINT, or a short geohash string all hash cheaply and stay cache-resident; a geometry-valued GROUP BY hashes serialized WKB byte-for-byte and can defeat the hash aggregate entirely, forcing a sort-based group. Project the cheap key first, aggregate, reconstruct geometry last.

Accuracy and cost pull in opposite directions and the right balance is workload-specific. A square grid is the cheapest and the least faithful: it is a single division per axis, but its directional neighbour bias distorts any spread or nearest-neighbour statistic. H3 costs a per-row encode on top of the hash pass but repays it with equal-area cells and a uniform neighbour ring, which is why it is the usual choice when the binned surface will feed adjacency analysis or be joined across data sources on a shared cell key. DBSCAN is the most expensive and the most expressive — it is the only member of the family that discovers group shape from the data rather than imposing a fixed lattice, and the only one that distinguishes signal from noise. When two methods would answer the question equally well, take the cheaper one; reach for DBSCAN only when the shape of the grouping is itself the finding.

Edge Cases & Anti-Patterns #

Square grid over a geographic CRS (silent wrong answers). floor(ST_X(geom) / 500) over EPSG:4326 divides degrees, producing cells hundreds of kilometres wide near the equator and collapsing toward the poles. Confirm the SRID before binning by metres:

SELECT DISTINCT ST_SRID(geom) AS srid FROM events;  -- expect a projected CRS, not 4326

DBSCAN eps in the wrong unit. Passing eps := 150 against degree-space geometry treats 150 degrees as the radius — the entire dataset collapses into one cluster. Project to a metric CRS first, or express eps in the degree magnitude that corresponds to your target distance (roughly 0.00135 degrees per 150 m at the equator, and less as latitude rises).

Treating DBSCAN noise as a real group. ST_ClusterDBSCAN returns NULL for points that never reach min_points; a GROUP BY cluster_id that does not filter or bucket the NULLs silently folds all noise into one phantom group. Handle it explicitly:

-- Separate real clusters from noise instead of letting NULL become a group.
SELECT
    coalesce(CAST(cluster_id AS VARCHAR), 'noise') AS grp,
    COUNT(*) AS n
FROM clustered
GROUP BY grp;

H3 axis order reversed. Calling h3_latlng_to_cell(ST_X(geom), ST_Y(geom), 8) passes longitude as latitude; every point lands in the wrong cell, and because valid H3 indexes still come back the error is silent. Latitude (ST_Y) must come first.

Geometry in the GROUP BY key. Grouping directly on geom — instead of a derived cell key — hashes raw WKB and defeats vectorization. Derive the scalar key in a projection, then group on it, exactly as the vectorized aggregations guide details.

Query Regression Analysis #

A production grouping pipeline needs a baseline plan and an automated diff so a degradation surfaces in CI rather than on a pager. Capture the plan as JSON, walk the tree, and assert the expected operator is present and no row-at-a-time or nested-loop node appears. This harness slots into the orchestration covered in Python and DuckDB integration workflows.

import duckdb
import json

con = duckdb.connect(":memory:")
con.execute("INSTALL spatial; LOAD spatial;")
con.execute("SET threads = 8; SET memory_limit = '12GB';")

QUERY = """
SELECT floor(ST_X(geom) / 500) AS cx, floor(ST_Y(geom) / 500) AS cy, COUNT(*) AS n
FROM events
GROUP BY cx, cy
"""

# FORMAT JSON yields a machine-readable plan tree for diffing across builds.
plan = json.loads(
    con.execute(f"EXPLAIN (FORMAT JSON) {QUERY}").fetchone()[1]
)

BANNED = {"NESTED_LOOP_JOIN", "ROW_EXECUTION"}
REQUIRED = "HASH_GROUP_BY"   # for DBSCAN pipelines, assert "WINDOW" instead

def walk(node):
    yield node.get("name", "")
    for child in node.get("children", []):
        yield from walk(child)

found = set(walk(plan[0])) if isinstance(plan, list) else set(walk(plan))
assert REQUIRED in found, f"Expected {REQUIRED}; grouping key may not be scalar. Saw: {found}"
assert not (found & BANNED), f"Plan regression: {found & BANNED} appeared in grouping plan"

Three fields are worth tracking across builds: operator_name (a shift away from HASH_GROUP_BY, or a NESTED_LOOP_JOIN, is a hard regression), operator_timing (per-node deltas localize a slowdown), and peak_memory plus any spill indicator — the last is the early warning for a DBSCAN partition that has outgrown memory_limit. To run this capture without blocking an event loop, see async execution patterns.

Related #

See also

• H3 hexagon binning in DuckDB — the equal-area hexagon method, resolution selection, and boundary reconstruction in full.
• ST_ClusterDBSCAN for spatial grouping — parameter tuning, noise handling, and window-partition sizing for density clustering.
• Vectorized aggregations — the key-then-aggregate mechanics every binning method shares.
• CRS mapping and transformations — the projection rules that decide whether a grid cell or eps radius means what you think.

Up: Modern Spatial SQL Query Patterns