Fixing a DuckDB Spatial Join That Never Finishes

On this page

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 nested loop over every pair
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 nested loop over every pair?
plan = con.execute("explain " + sql).fetchall()[0][1]
print("SPATIAL_JOIN" in plan.upper(), "NL_JOIN" 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,472,195 GeoNames points against Natural Earth's 4,596 provinces, counted per province, completed in 66.2 s with a peak of 304 MB on 24 threads, with the machine at a load average of 22 (Example 4). If yours is far slower than that shape suggests, something in the list above applies.

A nested loop tests every pair: 62 billion tests for 13.5 million points against 4,596 provinces.

Plan, CRS, complexity, filters and output size

1. Look at the plan before anything else

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

Look for SPATIAL_JOIN. If you see BLOCKWISE_NL_JOIN instead, DuckDB could not recognise the predicate as spatial and is testing every pair — for 13.5 million points against 4,596 provinces, 62 billion tests. On a 2,000-point sample the nested loop took 6.9 s where the spatial join took 0.4 s; at that rate the full join would run for about thirteen hours instead of a minute.

Measured with EXPLAIN in DuckDB 1.5.5, ST_Intersects, ST_Within, ST_Contains and ST_DWithin plan as SPATIAL_JOIN — and so does a predicate with a function around one geometry, such as ST_Intersects(ST_Buffer(p.geom, 0.01), a.geom). Two forms do not: a distance written as a comparison, ST_Distance(a, b) < d, and 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. Split an OR into two joins and UNION them.

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;

For a single polygon the cost of the exact test grows with its vertex count: a million points against one polygon took 0.68 s at 1,001 vertices, 5.47 s at 10,001 and 25.0 s at 50,001 (Example 4). If one or a few polygons have tens of thousands of vertices and most of the points fall in them, simplify for the join:

create table areas_simple as
select name, st_simplify(geom, 0.01) as geom from areas;   -- CRS units: degrees here

The tolerance is in the units of the coordinates — degrees for lon/lat data, so a tolerance of 50 meant as metres would reduce every province to a triangle or nothing. Check before assuming it will help: across all 4,596 Natural Earth provinces, where the largest has 38,527 vertices, 0.01° removed 71% of the vertices and cut the join from 66.4 s to 59.6 s — and moved 16,191 points to a different province or none (Example 2). Be explicit that this changes results at the boundary. For counting points in countries it is usually harmless; for a legal boundary it is not.

Vertex count matters when a few huge polygons hold most of the points — not across a table of ordinary ones.

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:

-- looks slow: joins everything, then discards
select count(*) from points p join areas a on st_intersects(p.geom, a.geom)
where a.admin = 'United Kingdom' and p.feature_class = 'P';

-- explicit: joins only what survives
with p as (select * from points where feature_class = 'P'),
     a as (select * from areas where admin = 'United Kingdom')
select count(*) from p join a on st_intersects(p.geom, a.geom);

For these two filters it makes no difference: the EXPLAIN of the first query already shows feature_class='P' and admin='United Kingdom' applied in the scans below the SPATIAL_JOIN. The optimiser pushes simple column filters down on its own. Writing them into the inputs costs nothing and is still the habit to keep for anything the optimiser cannot see through — a filter on a computed column, or one inside a subquery with a LIMIT.

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. An in-memory connection already spills to .tmp in the working directory by default; point temp_directory at fast local disk, and never set it to an empty string, which disables spilling and turns it into an out-of-memory error. Reducing threads lowers the total working set when several threads are each holding buffers.

Same meaning, and only one of each pair can be planned.

Diagnosing a join of GeoNames against provinces

The data used on this page

The points are every GeoNames feature (allCountries, about 420 MB zipped); the areas are Natural Earth's 1:10m states and provinces, 4,596 polygons. points is a view that builds each point from the Parquet file as it is read, which is how a large join is usually written.

import os
import zipfile

import duckdb
import requests


def fetch(url, path):
    if not os.path.exists(path):
        with requests.get(url, stream=True, timeout=600) as reply:
            reply.raise_for_status()
            with open(path, "wb") as out:
                for chunk in reply.iter_content(1 << 20):
                    out.write(chunk)


if not os.path.exists("places.parquet"):
    fetch("https://download.geonames.org/export/dump/allCountries.zip", "allCountries.zip")
    zipfile.ZipFile("allCountries.zip").extract("allCountries.txt")
    duckdb.sql("""
        copy (select column01 as name, column04 as lat, column05 as lon, column06 as feature_class,
                     column08 as country
              from read_csv('allCountries.txt', delim = '\t', header = false, quote = '',
                            sample_size = -1))
        to 'places.parquet' (format parquet)""")
if not os.path.exists("provinces.shp"):
    fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip",
          "provinces.zip")
    with zipfile.ZipFile("provinces.zip") as z:
        for member in z.namelist():
            ext = os.path.splitext(member)[1]
            if ext in (".shp", ".shx", ".dbf", ".prj", ".cpg"):
                open("provinces" + ext, "wb").write(z.read(member))

con = duckdb.connect()
con.install_extension("spatial")
con.load_extension("spatial")
con.execute("create or replace table areas as select name, admin, geom from st_read('provinces.shp')")
con.execute("""create or replace view points as
               select name, country, feature_class, st_point(lon, lat) as geom
               from read_parquet('places.parquet') where lat is not null""")
print(con.execute("select count(*) from points").fetchone()[0], "points,",
      con.execute("select count(*) from areas").fetchone()[0], "areas")

13472195 points, 4596 areas

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 any(op in plan for op in ("BLOCKWISE_NL_JOIN", "NESTED_LOOP_JOIN", "CROSS_PRODUCT")):
        print("   ! NESTED LOOP — the predicate was not recognised as spatial")
        print("     rewrite comparisons as ST_DWithin and split OR into a UNION")
    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}")


diagnose_join(con, "select count(*) from points p join areas a "
                   "on st_intersects(p.geom, a.geom)", "points", "areas")

print("--- how other predicates plan")
for predicate in ("st_within(p.geom, a.geom)", "st_dwithin(p.geom, a.geom, 0.01)",
                  "st_intersects(st_buffer(p.geom, 0.01), a.geom)",
                  "st_intersects(p.geom, st_convexhull(a.geom))",
                  "st_distance(p.geom, a.geom) < 0.01",
                  "st_intersects(p.geom, a.geom) or st_touches(p.geom, a.geom)"):
    plan = con.execute(f"explain select count(*) from points p join areas a on {predicate}"
                       ).fetchall()[0][1]
    operator = next(op for op in ("SPATIAL_JOIN", "BLOCKWISE_NL_JOIN", "NESTED_LOOP_JOIN",
                                  "CROSS_PRODUCT") if op in plan)
    print(f"   {predicate:60} {operator}")

plan = con.execute("""explain select count(*) from points p join areas a
                      on st_intersects(p.geom, a.geom)
                      where a.admin = 'United Kingdom' and p.feature_class = 'P'""").fetchall()[0][1]
print("--- filters below the join:",
      plan.index("SPATIAL_JOIN") < plan.index("feature_class='P'")
      and plan.index("SPATIAL_JOIN") < plan.index("admin='United Kingdom'"))

1. plan
   SPATIAL_JOIN — good
2. coordinate systems
   points               GEOMETRY
   areas                GEOMETRY('EPSG:4326')
3. geometry complexity
   points                 13,472,195 rows      13,472,195 vertices  worst         1
   areas                       4,596 rows       1,295,319 vertices  worst    38,527
4. expected output size
   9,625 pairs from 10,000 sampled rows → about 12,966,987 pairs in full
5. settings
   memory_limit     48.7 GiB
   temp_directory   .tmp
   threads          24
--- how other predicates plan
   st_within(p.geom, a.geom)                                    SPATIAL_JOIN
   st_dwithin(p.geom, a.geom, 0.01)                             SPATIAL_JOIN
   st_intersects(st_buffer(p.geom, 0.01), a.geom)               SPATIAL_JOIN
   st_intersects(p.geom, st_convexhull(a.geom))                 SPATIAL_JOIN
   st_distance(p.geom, a.geom) < 0.01                           BLOCKWISE_NL_JOIN
   st_intersects(p.geom, a.geom) or st_touches(p.geom, a.geom)  BLOCKWISE_NL_JOIN
--- filters below the join: True

Example 2 — simplifying only for the join

def join_with_simplified(con, points, areas, tolerance, name_column="name"):
    """Simplify for the join; keep the original geometry for anything measured.

    tolerance is in the coordinates' units: degrees for lon/lat data.
    """
    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}) 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 a tolerance of {tolerance} (in the CRS units)")
    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()


import time

start = time.time()
original = con.execute("""select a.name, count(*) as n from points p join areas a
                          on st_intersects(p.geom, a.geom) group by 1""").df()
print(f"original:   {original['n'].sum():,} points matched in {time.time() - start:.1f} s")
start = time.time()
simple = join_with_simplified(con, "points", "areas", 0.01)
print(f"simplified: {simple['n'].sum():,} points matched in {time.time() - start:.1f} s")
print(simple.head(3).to_string(index=False))

original:   12,950,244 points matched in 66.4 s
vertices 1,295,319 → 378,854 (29%) at a tolerance of 0.01 (in the CRS units)
  boundary results will differ slightly — check that is acceptable
simplified: 12,934,053 points matched in 59.6 s
         name      n
Uttar Pradesh 172470
   California 125909
        Texas 113796

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

def join_in_chunks(con, points, areas, column="lon",
                   edges=(-180, -108, -36, 36, 108, 180), out_table="joined"):
    """A last resort: bound the working set by joining one slice of a column at a time.

    Slice on a column range, not LIMIT/OFFSET: a LIMIT feeds the join through a
    single-threaded STREAMING_LIMIT, and every OFFSET re-reads the rows before it.
    """
    con.execute(f"drop table if exists {out_table}")
    for i, (low, high) in enumerate(zip(edges[:-1], edges[1:])):
        upper = "<=" if i == len(edges) - 2 else "<"          # the last band keeps its edge
        sql = f"""
            select a.name, count(*) as n
            from (select * from {points} where {column} >= {low} and {column} {upper} {high}) p
            join {areas} a on st_intersects(p.geom, a.geom)
            group by 1
        """
        con.execute(f"create table {out_table} as {sql}" if i == 0
                    else f"insert into {out_table} {sql}")
        print(f"  {column} {low} to {high}: done")

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


con.execute("""create or replace view points_lon as
               select name, lon, st_point(lon, lat) as geom
               from read_parquet('places.parquet') where lat is not null""")
start = time.time()
chunked = join_in_chunks(con, "points_lon", "areas")
print(f"{int(chunked['n'].sum()):,} points matched in {time.time() - start:.1f} s")

  lon -180 to -108: done
  lon -108 to -36: done
  lon -36 to 36: done
  lon 36 to 108: done
  lon 108 to 180: done
12,950,244 points matched in 82.1 s

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. How you slice matters. An earlier version of this example sliced with LIMIT 5000000 OFFSET …: it matched the same 12,950,244 points but took 215 s against 62 s unchunked in the same run, because EXPLAIN shows a STREAMING_LIMIT under the join, which reads its input on one thread. Slices on a column keep every thread busy — and a half-open last band would have lost the ten points that sit exactly on 180°.

Example 4 — the timings on this page

Each case runs in a fresh process so its peak memory is its own. The nested-loop case is timed on 2,000 points, because on all 13.5 million it would not finish in a working day.

%%writefile join_bench.py
import sys
import time

import duckdb


def peak_mb():
    """VmHWM resets at exec; ru_maxrss would report the parent's peak."""
    for line in open("/proc/self/status"):
        if line.startswith("VmHWM:"):
            return int(line.split()[1]) / 1024

con = duckdb.connect()
con.load_extension("spatial")
con.execute("create table areas as select name, geom from st_read('provinces.shp')")
points = "(select st_point(lon, lat) as geom from read_parquet('places.parquet') where lat is not null)"
sample = f"(select * from {points} using sample 2000 rows (reservoir, 42))"
case = sys.argv[1]
if case == "crs mismatch":
    con.execute("create or replace table areas as select name, "
                "st_transform(geom, 'EPSG:4326', 'EPSG:3857', always_xy := true) as geom from areas")
sql = {
    "healthy join": f"select count(*) from (select a.name, count(*) from {points} p "
                    f"join areas a on st_intersects(p.geom, a.geom) group by 1)",
    "crs mismatch": f"select count(*) from {points} p join areas a on st_intersects(p.geom, a.geom)",
    "sample, spatial join": f"select count(*) from {sample} p join areas a on st_intersects(p.geom, a.geom)",
    "sample, nested loop": f"select count(*) from {sample} p join areas a "
                           f"on st_distance(p.geom, a.geom) < 0.0",
}[case]
start = time.time()
rows = con.execute(sql).fetchone()[0]
print(f"{case:21} {rows:>10,} rows  {time.time() - start:6.1f} s  "
      f"peak {peak_mb():,.0f} MB")

import subprocess
import sys

print(f"load average {os.getloadavg()[0]:.1f} on {os.cpu_count()} threads")
for case in ("healthy join", "crs mismatch", "sample, spatial join", "sample, nested loop"):
    print(subprocess.run([sys.executable, "join_bench.py", case],
                         capture_output=True, text=True).stdout.strip())

# vertex count, isolated: one polygon, a million random points
con.execute("create or replace table pts as "
            "select st_point(random() * 2 - 1, random() * 2 - 1) as geom from range(1000000)")
for segments in (250, 2500, 12500):
    con.execute(f"create or replace table poly as "
                f"select st_buffer(st_point(0, 0), 0.9, {segments}) as geom")
    vertices = con.execute("select st_npoints(geom) from poly").fetchone()[0]
    start = time.time()
    con.execute("select count(*) from pts p join poly a on st_intersects(p.geom, a.geom)").fetchone()
    print(f"one polygon of {vertices:>6,} vertices: {time.time() - start:.2f} s")

load average 21.9 on 24 threads
healthy join               4,419 rows    66.2 s  peak 304 MB
crs mismatch                   0 rows    42.6 s  peak 310 MB
sample, spatial join       1,921 rows     0.4 s  peak 167 MB
sample, nested loop            0 rows     6.9 s  peak 166 MB
one polygon of  1,001 vertices: 0.68 s
one polygon of 10,001 vertices: 5.47 s
one polygon of 50,001 vertices: 24.97 s

A mismatched CRS cost more than half the healthy join to return nothing. The nested loop spent about 3.5 ms on each sampled point against the whole province table; the spatial join handled all 13.5 million at under 5 µs each — a difference of more than 600 times per point.

At 6.9 s per 2,000 points, the nested loop would need about thirteen hours for the full table.

Nested loops, vertices and silent mismatches

Why a nested loop 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 nested-loop join tests every pair. For 13.5 million points against 4,596 provinces that is 62 billion exact geometric predicate evaluations, each one non-trivial. Measured per point, the difference was more than 600 times, which is why the job appears to hang rather than to be slow.

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

When vertex count matters as much as row count

An exact ST_Intersects between a point and a polygon walks the polygon's rings. Measured on one polygon, fifty times the vertices cost thirty-seven times the time — close to linear.

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.

For the large case, move the aggregation into the join so pairs are discarded as counted.

Plannable predicates and wrapped geometries

explain first, always. A nested loop is the only cause with a penalty of hundreds of times.
ST_DWithin(a, b, d) plans; ST_Distance(a, b) < d may not.
A function around one geometry — ST_Buffer, ST_ConvexHull — still plans as SPATIAL_JOIN in DuckDB 1.5.5; precompute it into a column if it is expensive or used by more than one query.
Invalid geometry makes predicate results unreliable; run 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 BLOCKWISE_NL_JOIN 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.