Fixing a DuckDB Spatial Join That Never Finishes #

Problem statement #

The join has been running for forty minutes. The same join in GeoPandas took two, or crashed, and neither outcome tells you what is wrong.

A DuckDB spatial join that will not finish almost always has one of five causes, and they need different fixes:

• the predicate cannot be planned as a spatial join and has degraded to a cross product
• the geometries are enormous — a few multipolygons with hundreds of thousands of vertices
• the coordinate systems differ, so the engine is doing full work to produce nothing
• the join is genuinely many-to-many and the output is far larger than either input
• it is spilling because the build side does not fit

The first step is not to optimise. It is to find out which of the five it is, which takes about a minute.

Quick answer #

# 1. is it a spatial join at all, or a cross product?
plan = con.execute("explain " + sql).fetchall()[0][1]
print("SPATIAL_JOIN" in plan.upper(), "CROSS_PRODUCT" in plan.upper())

# 2. how big are the geometries?
con.execute("""
    select count(*)             as polygons,
           max(st_npoints(geom)) as max_vertices,
           avg(st_npoints(geom)) as mean_vertices
    from areas
""").df()

# 3. do the coordinate systems match?
con.execute("select typeof(geom) from areas limit 1").fetchone()
con.execute("select typeof(geom) from points limit 1").fetchone()

# 4. how big will the output be? test on a sample first
con.execute("""
    select count(*) from (select * from points using sample 10000 rows) p
    join areas a on st_intersects(p.geom, a.geom)
""").fetchone()

For reference, a healthy large join: 13,464,017 points against 4,596 polygons completed in 61.8 s using 305 MB. If yours is far slower than that shape suggests, something in the list above applies.

Step-by-step solution #

1. Look at the plan before anything else #

print(con.execute("explain " + sql).fetchall()[0][1])

Look for SPATIAL_JOIN. If you see CROSS_PRODUCT followed by a filter, DuckDB could not recognise the predicate as spatial and is testing every pair — which for a million rows against ten thousand is ten billion evaluations.

Predicates that plan properly are simple ones: ST_Intersects(a.geom, b.geom), ST_Within, ST_Contains, ST_DWithin. Predicates that often do not: anything wrapping a geometry in another function, ST_Distance(a, b) < 1000 written as an inequality, or an OR of two spatial conditions.

Rewrite ST_Distance(a.geom, b.geom) < 1000 as ST_DWithin(a.geom, b.geom, 1000) — same meaning, and the second one can be planned.

2. Check the coordinate systems #

A CRS mismatch does not raise. The engine does the whole join and finds nothing:

select typeof(geom) from points limit 1;   -- GEOMETRY  (no CRS)
select typeof(geom) from areas limit 1;    -- GEOMETRY('EPSG:4326')

If one side has been through ST_Transform, its type is a bare GEOMETRY and the annotation cannot help you. Transform both sides explicitly, with always_xy := true, and record what you did.

3. Measure the geometry complexity #

A spatial join's cost is driven by vertex count as much as by row count. One country polygon with 400,000 vertices costs more to test than ten thousand simple ones:

select count(*)                as n,
       sum(st_npoints(geom))   as total_vertices,
       max(st_npoints(geom))   as worst,
       avg(st_npoints(geom))   as mean
from areas;

If the maximum is in the hundreds of thousands, simplify for the join. A tolerance appropriate to your analysis usually removes most of the vertices without changing which points fall inside:

create table areas_simple as
select name, st_simplify(geom, 50) as geom from areas;

Be explicit that this changes results at the boundary. For counting points in countries it is harmless; for a legal boundary it is not.

4. Filter before joining, not after #

A WHERE clause on the joined result still requires the join to happen. Push the filter into the inputs:

-- slow: joins everything, then discards
select ... from points p join areas a on st_intersects(p.geom, a.geom)
where a.country = 'GB' and p.year = 2025;

-- fast: joins only what survives
with p as (select * from points where year = 2025),
     a as (select * from areas where country = 'GB')
select ... from p join a on st_intersects(p.geom, a.geom);

DuckDB's optimiser pushes many filters down automatically, but not all of them past a spatial join. Writing it explicitly costs nothing and removes the doubt.

5. Test the output size on a sample #

A join is not a lookup. If each point falls in several polygons — overlapping zones, buffers, a coverage where boundaries are shared — the output can be many times the input.

select count(*) from (select * from points using sample 10000 rows) p
join areas a on st_intersects(p.geom, a.geom);

Multiply by the sampling ratio. If ten thousand points produce eighty thousand rows, ten million will produce eighty million, and the join is not slow — it is large.

6. Give it a memory limit and somewhere to spill #

con.execute("set memory_limit = '4GB'")
con.execute("set temp_directory = '/fast/scratch'")
con.execute("set threads = 8")

A join whose build side does not fit will spill, which is slow but finite. Without a temp directory it fails instead. Reducing threads lowers the total working set when several threads are each holding buffers.

Code examples #

Example 1 — a diagnostic that checks all five causes #

import duckdb


def diagnose_join(con, sql, left, right, geom="geom", sample=10_000):
    print("1. plan")
    plan = con.execute("explain " + sql).fetchall()[0][1].upper()
    if "SPATIAL_JOIN" in plan:
        print("   SPATIAL_JOIN — good")
    elif "CROSS_PRODUCT" in plan or "NESTED_LOOP" in plan:
        print("   ! CROSS PRODUCT — the predicate was not recognised as spatial")
        print("     rewrite it as ST_Intersects / ST_Within / ST_DWithin "
              "with bare geometry arguments")
    else:
        print("   neither operator found; read the plan by hand")

    print("2. coordinate systems")
    for table in (left, right):
        t = con.execute(f"select typeof({geom}) from {table} limit 1").fetchone()
        print(f"   {table:20} {t[0] if t else 'empty'}")

    print("3. geometry complexity")
    for table in (left, right):
        row = con.execute(f"""
            select count(*), sum(st_npoints({geom})), max(st_npoints({geom}))
            from {table}""").fetchone()
        note = "  <- simplify for the join" if (row[2] or 0) > 100_000 else ""
        print(f"   {table:20} {row[0]:>12,} rows  "
              f"{row[1] or 0:>14,} vertices  worst {row[2] or 0:>9,}{note}")

    print("4. expected output size")
    est = con.execute(f"""
        select count(*) from (select * from {left} using sample {sample} rows) l
        join {right} r on st_intersects(l.{geom}, r.{geom})
    """).fetchone()[0]
    total = con.execute(f"select count(*) from {left}").fetchone()[0]
    print(f"   {est:,} pairs from {sample:,} sampled rows "
          f"→ about {int(est * total / sample):,} pairs in full")

    print("5. settings")
    for setting in ("memory_limit", "temp_directory", "threads"):
        value = con.execute(
            f"select current_setting('{setting}')").fetchone()[0]
        print(f"   {setting:16} {value}")

Example 2 — simplifying only for the join #

def join_with_simplified(con, points, areas, tolerance_m, name_column="name"):
    """Simplify for the join; keep the original geometry for anything measured."""
    before = con.execute(
        f"select sum(st_npoints(geom)) from {areas}").fetchone()[0]
    con.execute(f"""
        create or replace table _areas_join as
        select {name_column}, st_simplify(geom, {tolerance_m}) as geom
        from {areas}
    """)
    after = con.execute(
        "select sum(st_npoints(geom)) from _areas_join").fetchone()[0]
    print(f"vertices {before:,} → {after:,} ({100 * after / before:.0f}%) "
          f"at {tolerance_m} m tolerance")
    print("  boundary results will differ slightly — check that is acceptable")

    return con.execute(f"""
        select a.{name_column}, count(*) as n
        from {points} p join _areas_join a
          on st_intersects(p.geom, a.geom)
        group by 1 order by n desc
    """).df()

Example 3 — running the join in chunks when nothing else helps #

def join_in_chunks(con, points, areas, chunk_rows=2_000_000, out_table="joined"):
    """A last resort: bound the working set by processing the probe side in slices."""
    total = con.execute(f"select count(*) from {points}").fetchone()[0]
    con.execute(f"drop table if exists {out_table}")

    first = True
    for offset in range(0, total, chunk_rows):
        sql = f"""
            select a.name, count(*) as n
            from (select * from {points} limit {chunk_rows} offset {offset}) p
            join {areas} a on st_intersects(p.geom, a.geom)
            group by 1
        """
        if first:
            con.execute(f"create table {out_table} as {sql}")
            first = False
        else:
            con.execute(f"insert into {out_table} {sql}")
        print(f"  {min(offset + chunk_rows, total):,}/{total:,}")

    return con.execute(f"""
        select name, sum(n) as n from {out_table} group by 1 order by n desc
    """).df()

Chunking is rarely necessary — the engine streams — but it converts an unbounded run into a series of bounded ones with visible progress, which is sometimes worth more than raw speed.

Explanation #

Why a cross product is the difference between minutes and never #

SPATIAL_JOIN builds a spatial structure over one side and probes it with the other, so each row on the probe side tests against a small candidate set.

A cross product tests every pair. For a million rows against ten thousand polygons that is ten billion exact geometric predicate evaluations, each one non-trivial. The difference is not a factor of two; it is a factor of thousands, which is why the job appears to hang rather than to be slow.

Checking the plan first is therefore always the right first move.

Why vertex count matters as much as row count #

An exact ST_Intersects between a point and a polygon walks the polygon's rings. A polygon with 400,000 vertices costs roughly four hundred times a polygon with a thousand.

The bounding-box filter that precedes the exact test removes most candidates, so the cost concentrates on the pairs that actually overlap — which for a coverage means every point pays the full cost of the polygon it falls in. Simplifying the polygons for the join attacks exactly that cost.

Why a CRS mismatch looks like slowness #

The engine does not know the coordinate systems differ. It builds the structure, probes it, tests candidates, and finds nothing — the full cost of a join, for an empty result.

On large inputs that can take minutes before returning zero rows, which reads as "slow" until you notice the result. Checking typeof(geom) on both sides takes a second and rules it out.

Why the output size is often the real answer #

Not every slow join is inefficient. A join that legitimately produces eighty million rows from ten million is doing eight times the work of one that produces ten million, and there is no optimisation for that beyond aggregating earlier.

The fix, when it applies, is to move the aggregation into the join: count(*) grouped by area rather than a materialised list of pairs. The engine can then discard each matched pair as soon as it has counted it.

Edge cases or notes #

• explain first, always. A cross product is the only cause with a thousand-fold penalty.
• ST_DWithin(a, b, d) plans; ST_Distance(a, b) < d may not.
• Wrapping a geometry in a function — ST_Buffer, ST_Centroid — inside the predicate can defeat the planner. Precompute it into a column.
• Invalid geometry can slow the exact test dramatically; ST_IsValid on untrusted input.
• A CRS mismatch costs the full join for zero rows.
• using sample n rows is the cheapest way to estimate output size.
• Aggregate inside the join rather than materialising pairs when the output is large.
• Simplify for the join only, and say so — boundary results change.

Internal links #

• How to run a spatial join in DuckDB — the operation done properly
• How DuckDB uses (and does not use) a spatial index — why an index will not save you
• CRS in DuckDB: why ST_Transform moves your data to the wrong place — the zero-row cause
• DuckDB's memory model: streaming, spilling and limits — the settings that stop a spill from failing
• Fixing DuckDB out of memory on a large spatial query — when it fails instead of dragging
• How to simplify geometry in GeoPandas — the simplification step
• Spatial join cardinality explained — why the output grows
• Fixing a slow PostGIS query — the same diagnosis in the other engine

FAQ #

Why is my DuckDB spatial join so slow? #

Check the plan first. If it shows a cross product rather than SPATIAL_JOIN, the predicate was not recognised as spatial and every pair is being tested.

Would a spatial index help? #

No. DuckDB's table-to-table spatial join uses a dedicated operator, not the R-tree. The index accelerates filters against a fixed geometry.

Why does my join return zero rows after running for minutes? #

A CRS mismatch. DuckDB does not raise on one, so it performs the entire join and finds nothing. Check typeof(geom) on both sides.

Should I simplify my polygons? #

For the join, often yes — cost scales with vertex count, and one polygon with 400,000 vertices dominates. Keep the originals for anything you measure, and accept that boundary results change.

How do I know if the join is slow or just large? #

Run it on a sample: using sample 10000 rows. Multiply the pair count by the sampling ratio. If the output is eight times the input, the join is large, not slow.

What settings should I check? #

memory_limit, temp_directory and threads. Without a temp directory a spilling join fails instead of finishing; reducing threads lowers the total working set.