DuckDB Spatial Analytical SQL for GIS

Replacing Python Loops with Set-Based SQL

A loop over query results is the most expensive thing a spatial pipeline can contain, and it is almost always standing in for a construct SQL already has — this walkthrough, part of the Python spatial UDFs and predicate pushdown reference, maps the four loop shapes that recur in geospatial Python onto their SQL equivalents, shows why the nested one degrades catastrophically rather than gradually, and gives the incremental workflow that makes the rewrite debuggable.

Root-Cause Analysis: why the loop was written and why it stops working

Loops are not written out of ignorance. They are written because five things are true at the moment of writing, and the fifth stops being true later.

  • The logic was easier to express procedurally. A condition with several branches reads naturally as an if and awkwardly as a CASE, at least until the CASE has been written once.
  • The data was small. On a thousand rows the loop finishes instantly, so nothing signals that it is the wrong shape.
  • The result had to go somewhere Python-shaped. A list, a dictionary, a call into a library — all of which felt like they required iteration to build.
  • A debugger was available. Stepping through a loop is a familiar workflow; inspecting an intermediate result set is not, until it is.
  • Nobody expected the input to grow. This is the one that changes. A nested loop is quadratic, so a tenfold growth in each input is a hundredfold growth in work, which turns thirty seconds into an hour without any code changing.

The distinguishing question is whether the loop is nested. A single loop over rows is linear and merely wasteful; a nested loop over two collections is quadratic and will eventually stop finishing.

Four loop shapes and the SQL each one is asking for A filtering loop is a WHERE; an accumulating loop is a GROUP BY; a nested loop is a join; a look-back loop is a window function. THE LOOP YOU WROTE WHAT IT IS ASKING FOR WHAT SQL ADDS filter each row WHERE index pruning before evaluation accumulate per key GROUP BY a parallel hash aggregation compare every pair a spatial join envelope pruning of the pairs look at the previous row a window function one sort instead of N lookups

In each row the SQL form is shorter as well as faster, which is the part that survives review.

Deterministic Configuration 

import duckdb

con = duckdb.connect("gis.duckdb")
con.execute("INSTALL spatial; LOAD spatial")

# The rewrite moves work into the engine, so the engine now needs the budget
# the Python process used to consume.
con.execute("SET memory_limit = '8GB'")
con.execute("SET threads = 8")
con.execute("SET temp_directory = '/var/tmp/duckdb_rewrite'")

Optimized Execution Pattern

The nested loop is the case worth working through, because it is the one that fails rather than merely wasting time.

# ANTI-PATTERN: every point against every polygon. The pruning an index would
# have done for free does not happen, so this is quadratic by construction.
points   = con.execute("SELECT id, ST_AsWKB(geom) FROM incidents").fetchall()
polygons = con.execute("SELECT zone_id, ST_AsWKB(geom) FROM zones").fetchall()

matches = []
for pid, pwkb in points:
    p = from_wkb(pwkb)
    for zid, zwkb in polygons:              # ← the quadratic
        if from_wkb(zwkb).contains(p):
            matches.append((pid, zid))
            break
-- PATTERN: the same question as a join. The bounding-box operator in ON lets
-- the R-tree eliminate the overwhelming majority of pairs before any vertex
-- is compared, which is exactly the work the loop was doing by hand.
CREATE OR REPLACE TABLE incident_zone AS
SELECT i.id AS incident_id, z.zone_id
FROM incidents i
JOIN zones z
  ON  z.geom && i.geom
WHERE ST_Contains(z.geom, i.geom);
A nested loop and a spatial join, as the inputs grow A Python nested loop grows with the product of the inputs while an indexed join grows with the number of genuine matches, so the gap widens from imperceptible to an hour. INPUTS PYTHON NESTED LOOP INDEXED SPATIAL JOIN 1,000 × 100 100,000 comparisons · <1 s <0.1 s — no visible difference 10,000 × 1,000 10,000,000 · ~30 s <0.5 s 100,000 × 10,000 1,000,000,000 · ~1 hour <1 s The loop is not slow because Python is slow. It is slow because it compares pairs that an envelope test would have eliminated for free.

The first row is why the loop was written. The third is why it has to be replaced.

The three rows of that comparison explain the whole life cycle of such a loop: it was written when the difference was invisible, it survived when the difference was tolerable, and it has to be replaced now that the difference is an hour. Nothing about the code changed in between.

Building the rewrite so it can be debugged

The genuine cost of moving from a loop to a query is the loss of a debugger, and the workflow that replaces it is to build the query as a chain of named CTEs, each of which can be selected from on its own. That gives the same incremental inspection a debugger provides, with the advantage that each stage is inspectable over the whole dataset rather than one row at a time.

-- Each CTE is a checkpoint. Comment out everything after any one of them and
-- SELECT * FROM it to see exactly what that stage produced.
WITH valid_points AS (
    SELECT id, geom FROM incidents WHERE ST_IsValid(geom)
), candidates AS (
    SELECT p.id, z.zone_id, z.geom AS zgeom, p.geom AS pgeom
    FROM valid_points p JOIN zones z ON z.geom && p.geom          -- envelope stage
), confirmed AS (
    SELECT id, zone_id FROM candidates WHERE ST_Contains(zgeom, pgeom)
)
SELECT zone_id, count(*) AS incidents FROM confirmed GROUP BY zone_id;

The count at each stage is also the diagnostic: if candidates is close to the product of the two inputs, the envelope stage is not pruning and the index is not being used, which is a different problem from the query being wrong.

Three objections to the rewrite, weighed Complexity is usually unfamiliarity; the debugging objection is real and addressable with incremental CTEs; the "small enough" objection is about today rather than about the code. OBJECTION WHAT IT IS WORTH THE ANSWER “too complex for SQL” usually unfamiliarity CTEs and window functions cover it “I need to step through it” a real cost build it as a chain of CTEs “the data is small” true today datasets grow; the code stays

Only the middle objection survives, and it has a workflow rather than a rebuttal.

Diagnostic Queries & Plan Validation

The rewrite has to be verified against the loop it replaces, and row count alone is not enough — a join that fans out produces more rows and a predicate written backwards produces fewer, and both are plausible.

-- Compare against the loop's output on the same input. The checksum catches a
-- reversed predicate that the row count would not: ST_Contains and ST_Within
-- with swapped arguments can return the same number of rows and different rows.
SELECT count(*)                       AS rows_out,
       sum(hash(incident_id, zone_id)) AS checksum
FROM incident_zone;

Diagnostic — did the envelope stage prune? Compare the candidate count against the product of the input sizes. Anything close to the product means the index was not consulted, which usually means the && is missing from the ON clause or the geometry column is a raw blob.

SELECT (SELECT count(*) FROM incidents) * (SELECT count(*) FROM zones) AS pairs_possible,
       (SELECT count(*) FROM incidents i JOIN zones z ON z.geom && i.geom) AS pairs_after_bbox;

Geometry Validation & Fallback Routing

Where a genuine gap remains — a library call with no SQL equivalent — the rewrite does not disappear, it shrinks. The set-based query does everything it can, and Python receives the reduced result rather than the raw inputs.

# The hybrid that survives: SQL reduces, Python handles only what SQL cannot,
# and the reduced result is small enough that a loop over it is harmless.
reduced = con.execute("""
    SELECT i.id, ST_AsWKB(i.geom) AS wkb
    FROM incidents i JOIN zones z ON z.geom && i.geom
    WHERE ST_Contains(z.geom, i.geom) AND z.zone_id = 'W12'
""").fetch_arrow_table()

# A few thousand rows rather than a few million — and even here, the
# vectorised call beats a loop.
values = _library_call(from_wkb(reduced.column("wkb").to_numpy(zero_copy_only=False)))

Frequently Asked Questions

Is every loop worth rewriting?

No — a loop over a result set of a few hundred rows costs nothing and may well be clearer. What is always worth rewriting is a nested loop, because its cost grows with the product of the inputs, and a dataset that doubles turns thirty seconds into two minutes and then into an hour without anything in the code changing.

How do I debug a query the way I debugged a loop?

Build it as a chain of named CTEs and select from each one in turn. That gives the same incremental inspection, over the whole dataset rather than one row at a time, and the row count at each stage doubles as a diagnostic — a candidate count close to the product of the inputs means the pruning stage is not working.

The logic has several branches. Can SQL express that?

Yes, with CASE, and usually more compactly than the if/elif chain it replaces. The genuine limits are recursion and stateful accumulation that cannot be expressed as a window, and both are rarer in spatial work than they feel — most apparent state turns out to be “the previous row”, which is lag.

How do I know the rewrite is correct?

Compare a row count and a checksum against the loop’s output on the same input. The count alone misses a reversed predicate — ST_Contains and ST_Within with swapped arguments can return the same number of rows and a different set of them — and that is exactly the mistake a mechanical translation makes.

What if part of the work genuinely cannot move?

Then it shrinks rather than disappears. Let SQL do every filter, join and aggregation it can, and hand Python the reduced result. A loop over a few thousand rows that a query has already narrowed is a different proposition from a loop over the raw input, and it is usually the honest end state rather than a compromise.

Does the rewrite always use less memory?

Usually much less, because the intermediate results stay inside the engine where they can spill, rather than accumulating as Python objects that cannot. The exception is a query that materialises a large intermediate the loop was consuming incrementally — which is worth knowing about, and is what a bounded batch stream is for.

Up: Python Spatial UDFs and Predicate Pushdown