DuckDB Spatial Analytical SQL for GIS

Tuning ST_DWithin Radius and Join Order

Two decisions dominate the cost of a proximity join and neither is visible in the SQL: how large the radius is relative to the local density, and which relation becomes the hash build side. This walkthrough, part of the spatial joins and proximity filters reference, covers deriving the radius from density rather than convention, measuring which side should build, and the plan signals that say the optimizer chose differently from you.

Root-Cause Analysis: why a proximity join costs what it does

  • The radius sets the candidate count, quadratically. Doubling the radius roughly quadruples the area searched and therefore the candidate pairs. A radius chosen with a safety margin is not conservative; it is expensive by the square of the margin.
  • Density varies across the layer. One radius that yields forty candidates downtown yields four in a suburb and none in the countryside, so a single global value is wrong nearly everywhere.
  • The build side is chosen by rows, not vertices. The optimizer builds its hash table from the relation it believes is smaller, and for spatial data the useful measure of size is total vertex count. The two disagree routinely.
  • The unit is whatever the coordinates are. A radius of 500 means 500 degrees on a geographic layer, which is most of the planet, and 500 metres on a projected one. The same query changes meaning entirely with the frame.
  • The exact predicate runs on survivors, not on candidates. ST_DWithin can answer from envelopes when they are further apart than the radius, so its cost is dominated by the pairs that survive the envelope test rather than by all pairs.

The distinguishing question is how many candidates the envelope stage produces per genuine match. That ratio is the tuning target, and both decisions on this page move it.

Candidate pairs scale with the square of the radius On two million points: 50 m yields 800k pairs in a second, 500 m yields 80M pairs in nearly two minutes. RADIUS CANDIDATE PAIRS RUNTIME 50 m ~800,000 ~1 s 100 m ~3,200,000 ~4 s 250 m ~20,000,000 ~26 s 500 m ~80,000,000 ~110 s A safety margin of two costs a factor of four; a margin of ten costs a factor of a hundred. Margins on a radius are not cheap.

Nothing about the data changed across these four rows. Only the number in the predicate did.

Deterministic Configuration 

INSTALL spatial; LOAD spatial;

SET memory_limit = '8GB';
SET threads = 8;
SET temp_directory = '/var/tmp/duckdb_join';

-- Build the index on the side that will be probed, after the data is loaded.
CREATE INDEX idx_shops_geom ON shops USING RTREE (geom);

Optimized Execution Pattern

The pattern is a two-stage predicate with a density-derived radius, and an explicit check on which side is heavier before assuming the optimizer chose well.

-- ANTI-PATTERN: a round-number radius applied globally, and no envelope stage
-- for the index to serve. Both problems compound.
SELECT i.incident_id, s.shop_id
FROM incidents i JOIN shops s
  ON ST_DWithin(i.geom, s.geom, 500);        -- 500 what? and no && to prune with
-- PATTERN: envelope stage first, then the distance guard, with the radius
-- taken per cell from a density table rather than fixed globally.
SELECT i.incident_id, s.shop_id, ST_Distance(i.geom, s.geom) AS metres
FROM incidents i
JOIN cell_radius r
  ON  r.cx = floor(ST_X(i.geom) / 1000)::INT
 AND  r.cy = floor(ST_Y(i.geom) / 1000)::INT
JOIN shops s
  ON  s.geom && ST_Expand(i.geom, r.r0)      -- index-servable square window
WHERE ST_DWithin(i.geom, s.geom, r.r0);      -- circular truth, on survivors

ST_Expand produces a square window that is a superset of the circle, which is exactly what makes it safe as a pre-filter: it lets through some corner pairs that the distance test then removes, and it never excludes a pair the distance test would have kept.

Choosing the build side by weight, not by rows Row count is the default signal and the wrong one; total vertex count predicts the memory the hash table occupies. SIGNAL WHAT IT PREDICTS USE IT? row count nothing, for geometry no — the default, and wrong total vertex count build-side memory yes average vertices per row which layer is dense as a diagnostic rows × avg vertices the number to compare yes — this is the answer

Two million points weigh 30 MB; fifty thousand boundaries weigh 1.1 GB.

Deriving the radius from density

A radius is only meaningful relative to how far apart the features are, and that varies by orders of magnitude within a single national layer. The derivation is straightforward: count features per grid cell once, invert to get a mean spacing, and scale it so the expected candidate count lands somewhere useful — a small multiple of however many neighbours the question wants.

Stored as a per-cell column, that estimate costs one join and removes the whole problem of a global value being wrong in most places. It also makes the query self-documenting: the radius is visibly a function of local density rather than a number somebody once chose.

-- One pass over the probed layer gives a per-cell radius that lands a
-- comparable candidate count everywhere, instead of forty downtown and none
-- in the countryside.
CREATE OR REPLACE TABLE cell_radius AS
SELECT floor(ST_X(geom) / 1000)::INT AS cx,
       floor(ST_Y(geom) / 1000)::INT AS cy,
       2.0 * sqrt(1000.0 * 1000.0 / greatest(count(*), 1)) AS r0
FROM shops
GROUP BY cx, cy;

Diagnostic Queries & Plan Validation

The plan says which side built and how many candidates survived, and both numbers are more useful than the runtime.

-- Look for the hash build side and the candidate count out of the join.
-- A build that spills before any output row appears means the heavy side
-- was chosen to build.
EXPLAIN ANALYZE
SELECT count(*)
FROM incidents i JOIN shops s ON s.geom && ST_Expand(i.geom, 250)
WHERE ST_DWithin(i.geom, s.geom, 250);

Compare candidates against matches. A ratio in the low tens is healthy for a proximity join; a ratio in the thousands means the radius is far larger than the density warrants, and shrinking it is worth more than any other change.

Four proximity-join symptoms No results means degrees; spilling before output means the wrong build side; a huge candidate ratio means the radius; regional variation means density. SYMPTOM POINTS AT FIX returns nothing at all the radius is in degrees reproject to a metric frame spills before any output the heavy side is building reduce or swap the build side candidates ›› matches the radius is too large derive it from density fast in one region, slow in another density variation a per-cell radius

Only the first is a correctness bug. The other three are the same bug about size.

Reducing the heavy side before the join

When the vertex-heavy relation genuinely has to participate, the productive move is to make it lighter before the join rather than to fight the optimizer over which side builds. Two reductions apply. The first is an attribute filter: most proximity questions are asked about a subset — a category of shop, a period of incidents — and applying that before the join reduces both sides. The second is simplification: a boundary layer used only to answer “is this within 250 m” does not need centimetre vertices, and ST_SimplifyPreserveTopology at a tolerance well below the radius changes no answers while removing most of the weight.

Both are worth trying before reaching for a hint, partly because DuckDB offers no join hints and partly because they are usually larger wins than the join order would have been.

Geometry Validation & Fallback Routing

Where the join still will not fit, partition it on a key both sides share and run the pieces.

-- Partition-aligned proximity join. Correct whenever no genuine match can
-- cross a region boundary; where one can, widen the region rather than
-- distribute the join.
SELECT i.incident_id, s.shop_id
FROM incidents i
JOIN shops s
  ON  s.region = i.region                      -- restrict the pair space first
 AND  s.geom && ST_Expand(i.geom, 250)
WHERE ST_DWithin(i.geom, s.geom, 250);

Frequently Asked Questions

Why does my ST_DWithin return no rows?

Almost always a unit mismatch. On a geographic layer the radius is in degrees, so 500 means five hundred degrees and 0.0005 means roughly fifty metres. Check the coordinate range before adjusting the number: values inside plus or minus 180 mean the layer is still in degrees and the fix is a reprojection.

How should I choose the radius?

From local density rather than as a round number. Count features per cell once, invert to a mean spacing, and scale so the expected candidate count is a small multiple of the neighbours you want. A global radius is wrong nearly everywhere, and its cost grows with the square of how wrong it is.

Do I still need the bounding-box operator with ST_DWithin?

For the index to be used, yes. ST_DWithin can answer from envelopes internally, but the optimizer needs an index-eligible expression to route through the R-tree, and ST_Expand plus && is that expression. The distance test then removes the corner pairs the square window let through.

Which side should build the hash table?

The one with fewer total vertices, which is frequently not the one with fewer rows. Two million points weigh about thirty megabytes as a build side; fifty thousand municipal boundaries weigh over a gigabyte. Measure sum(ST_NPoints(geom)) before assuming the row counts explain anything.

Can I force the join order?

There are no join hints, so the productive route is to change what the optimizer sees: filter the heavy relation first, simplify its geometry to a tolerance well below the radius, or materialise a reduced version of it. Those are usually larger wins than the join order would have been anyway.

Is ST_Expand plus && the same as ST_DWithin?

No, and the difference is what makes the pair work. ST_Expand gives a square window; ST_DWithin describes a circle. The square is a superset, so it is safe as a pre-filter and unsafe as a replacement — used alone it returns corner features further away than the radius allows.

Up: Spatial Joins & Proximity Filters