⚙️ How-To DuckDB & analytic SQL Joins & relationships DuckDB for Spatial Data · step 6 of 20 Updated

How to Run a Spatial Join in DuckDB

Problem statement

A spatial join is the operation that dominates most GIS workloads: which points fall in which polygons, which parcels touch which flood zone, which roads cross which boundary.

In GeoPandas it is one line and it holds both datasets in memory. That is fine until it is not — and the point at which it stops being fine is measurable. Joining 13,464,017 points to 4,596 polygons:

engine       time      peak memory     result
DuckDB      61.8 s        305 MB       12,942,217 pairs
GeoPandas  117.5 s      4,739 MB       12,942,217 pairs

Same answer, 1.9× the speed and 15.5× less memory. The memory figure is the one that decides whether the job runs at all on an ordinary machine.

But DuckDB is not automatically the right choice: on a small join — 7,342 points against the same 4,596 polygons — GeoPandas was faster, 0.104 s against 0.168 s.

Quick answer

A spatial join is a JOIN with a spatial predicate:

select p.name, a.name as area
from read_parquet('places.parquet') p
join st_read('areas.shp') a
  on st_intersects(st_point(p.lon, p.lat), a.geom);

DuckDB recognises the predicate and plans a dedicated SPATIAL_JOIN operator rather than a nested loop — which is what makes 62 billion candidate pairs finish in a minute.

import duckdb

con = duckdb.connect()
con.execute("install spatial; load spatial;")
con.execute("set enable_progress_bar = false")
result = con.execute(query).df()
Bar chart of spatial join times at three scales for DuckDB and GeoPandas.
Below the crossover the fixed costs of planning and structure building dominate.

Step-by-step solution

1. Choose the predicate deliberately

Predicate True when Typical use
ST_Intersects they share any point the default; points in polygons
ST_Within / ST_Contains one is entirely inside the other strict containment
ST_Touches they share a boundary but no interior adjacency
ST_Crosses interiors cross a road crossing a boundary
ST_DWithin within a distance "shops within 500 m"

ST_Intersects and ST_Within differ on the boundary: a point exactly on a polygon edge intersects it and is not within it. On a coverage of adjacent polygons, ST_Intersects can therefore produce two matches for one point, which is usually not what you want in a count.

2. Build point geometry inline when the source is coordinates

Most large point datasets arrive as latitude and longitude columns, not as geometry. Constructing the point in the join is cheaper than materialising a geometry column first:

join areas a on st_intersects(st_point(p.lon, p.lat), a.geom)

ST_Point(x, y) is longitude first, matching every GIS file convention.

3. Check the coordinate systems before you believe the answer

DuckDB does not raise on a CRS mismatch. It returns zero rows:

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

A join across mismatched systems produces a count of zero that looks exactly like "nothing matched". Check the types, or transform one side explicitly with always_xy := true.

4. Aggregate in the query

The reason to use a database engine is that the aggregate never leaves it. Counting points per polygon returns 51 rows rather than 2.2 million:

select a.name, count(*) as points, avg(p.population) as mean_population
from read_parquet('places.parquet') p
join st_read('states.shp') a
  on st_intersects(st_point(p.lon, p.lat), a.geom)
where p.country = 'US'
group by 1
order by points desc;

Measured: this exact shape ran in 4.27 s with 227 MB peak memory on 2,241,395 points against 51 polygons, against 15.5 s and 988 MB for the GeoPandas equivalent.

5. Count what did not match

An inner join silently drops unmatched rows, and the number of them is usually informative:

select count(*) filter (where a.name is null) as unmatched,
       count(*) filter (where a.name is not null) as matched
from read_parquet('places.parquet') p
left join st_read('states.shp') a
  on st_intersects(st_point(p.lon, p.lat), a.geom);

In the measured US join, 2,241,395 points produced 2,198,704 matches — 42,691 points, 1.9%, fell outside every state polygon. Offshore features, territories and generalised coastlines account for most of it, and knowing the number stops it becoming a surprise in a total.

6. Watch the row count on many-to-many joins

A join is not a lookup. A point on a shared boundary matches two polygons; a road crossing three counties produces three rows. The output can be larger than either input, which is correct and is the usual reason a "join" produces more rows than expected.

If you need one row per input feature, aggregate — count(*), string_agg, or any_value — rather than assuming the join is one-to-one.

Grid of five spatial predicates with when each is true and its typical use.
The boundary case is what makes a count differ between the first two.

Code examples

Example 1 — the join, with the checks that make it trustworthy 

import duckdb


def spatial_join_summary(points_source, polygons_source, polygon_name="name",
                         lon="lon", lat="lat", predicate="st_intersects"):
    """Join, aggregate, and report what did not match."""
    con = duckdb.connect()
    con.execute("install spatial; load spatial;")
    con.execute("set enable_progress_bar = false")

    con.execute(f"create table areas as select * from {polygons_source}")
    crs = con.execute("select typeof(geom) from areas limit 1").fetchone()[0]
    print(f"polygon geometry type: {crs}")
    if crs == "GEOMETRY":
        print("  ! no CRS on the polygons — a mismatch will return zero rows silently")

    totals = con.execute(f"""
        select count(*) as total,
               count(a.{polygon_name}) as matched
        from {points_source} p
        left join areas a
          on {predicate}(st_point(p.{lon}, p.{lat}), a.geom)
    """).fetchone()
    unmatched = totals[0] - totals[1]
    print(f"{totals[0]:,} point-rows, {totals[1]:,} matched, "
          f"{unmatched:,} unmatched ({100 * unmatched / totals[0]:.1f}%)")

    return con.execute(f"""
        select a.{polygon_name} as area, count(*) as n
        from {points_source} p
        join areas a on {predicate}(st_point(p.{lon}, p.{lat}), a.geom)
        group by 1 order by n desc
    """).df()

The left join count runs first deliberately. A zero-match result is far cheaper to diagnose before you have spent a minute on the aggregate.

Example 2 — a distance join with an explicit radius 

def within_distance(con, points_table, sites_table, metres, crs="EPSG:27700"):
    """Everything within `metres` of a site. Requires a projected CRS —
    a distance in degrees is not a distance."""
    con.execute(f"""
        create or replace table points_proj as
        select * exclude geom,
               st_transform(geom, 'EPSG:4326', '{crs}', always_xy := true) as geom
        from {points_table}
    """)
    con.execute(f"""
        create or replace table sites_proj as
        select * exclude geom,
               st_transform(geom, 'EPSG:4326', '{crs}', always_xy := true) as geom
        from {sites_table}
    """)
    return con.execute(f"""
        select s.name as site, count(*) as within_{metres}m,
               round(min(st_distance(p.geom, s.geom)), 1) as nearest_m
        from points_proj p
        join sites_proj s on st_dwithin(p.geom, s.geom, {metres})
        group by 1 order by 2 desc
    """).df()

Reprojecting both sides into a projected CRS before a distance join is not optional. ST_DWithin on lon/lat measures degrees, and a degree is 111 km at the equator and 56 km at 60° north.

Example 3 — verifying against GeoPandas on a sample 

def cross_check(con, points_source, polygons_path, sample=50_000):
    """Run both engines on a sample and assert they agree."""
    import geopandas as gpd

    duck = con.execute(f"""
        select count(*) from (select * from {points_source} limit {sample}) p
        join st_read('{polygons_path}') a
          on st_intersects(st_point(p.lon, p.lat), a.geom)
    """).fetchone()[0]

    frame = con.execute(f"select lon, lat from {points_source} limit {sample}").df()
    polygons = gpd.read_file(polygons_path)
    points = gpd.GeoDataFrame(frame,
                              geometry=gpd.points_from_xy(frame.lon, frame.lat),
                              crs=polygons.crs)
    pandas_count = len(gpd.sjoin(points, polygons, predicate="intersects"))

    print(f"duckdb {duck:,}   geopandas {pandas_count:,}   "
          f"{'agree' if duck == pandas_count else 'DISAGREE'}")
    assert duck == pandas_count, (
        "the two engines disagree — usually a CRS difference or a boundary "
        "predicate difference; check before trusting either")

Run this once when you move a pipeline from GeoPandas to DuckDB. A faster answer that differs from the one you trusted is not a speed-up.

Explanation

Why DuckDB's spatial join is not an index lookup

EXPLAIN on a spatial join shows a dedicated SPATIAL_JOIN operator, not an index scan. Even with an R-tree index present on the polygon table, the join plan does not use it:

┌───────────────────────────┐
│        SPATIAL_JOIN       │
│      Join Type: INNER     │
│ ST_Intersects(geom, geom) │
└───────────────────────────┘

The operator builds its own spatial structure over one side and probes it with the other, which is the right strategy for joining two whole tables. The R-tree index is for filtering a table by a fixed geometry — a bounding-box query — and DuckDB uses it there.

Knowing this stops a common wasted afternoon: creating an index and finding the join is no faster, because it was never going to use it.

Why the memory difference is larger than the speed difference

The global join was 1.9× faster and 15.5× lighter. GeoPandas materialises both GeoDataFrames and the join result in Python objects; DuckDB streams the probe side through the join and emits aggregated rows.

That is why the same query completed with memory_limit set to 200 MB — measured at 4.18 s, no slower than the unconstrained run. The engine was never holding the data in the first place.

Why GeoPandas wins on small joins

DuckDB pays fixed costs: parsing and planning the query, loading the extension, building the join structure. On 7,342 points against 4,596 polygons those costs dominate, and GeoPandas' in-memory R-tree does the whole job in 0.104 s against DuckDB's 0.168 s.

The crossover is somewhere in the low hundreds of thousands of rows, and it depends on the shapes involved. Below it, use whatever is convenient; above it, the difference stops being a matter of taste.

Why unmatched rows deserve a number

Every spatial join has an unmatched population, and it is never zero on real data. Measured: 1.9% of US-coded points fell outside the 51 state polygons, and 3.9% of global points fell outside every province.

Those points are offshore features, small islands the generalised boundaries omit, and slivers between polygons. None of that is an error, and all of it changes a total. Counting the unmatched rows converts an invisible discrepancy into a documented one.

Bar chart of matched and unmatched point shares in two spatial joins.
An inner join drops these silently; a left join with a filtered count reports them.

Edge cases or notes

  • ST_Intersects includes the boundary; ST_Within does not. A point on a shared edge matches two polygons with the first.
  • A CRS mismatch returns zero rows, not an error. Check the geometry column types.
  • ST_DWithin measures in the coordinate units — reproject before using it.
  • Invalid geometry can produce wrong results rather than an error; ST_IsValid before a join on untrusted data.
  • The output can be larger than either input. Aggregate if you need one row per feature.
  • An R-tree index does not speed up a table-to-table join — it speeds up a filter against a fixed geometry.
  • left join plus a count(... ) filter gives matched and unmatched in one pass.
  • Cross-check against GeoPandas once when porting a pipeline.

FAQ

How do I do a spatial join in DuckDB?

A normal JOIN with a spatial predicate: join areas a on st_intersects(st_point(p.lon, p.lat), a.geom). DuckDB plans a dedicated spatial join operator for it.

Is DuckDB faster than a GeoPandas sjoin?

At scale, yes. 13.5 million points against 4,596 polygons took 61.8 s and 305 MB in DuckDB against 117.5 s and 4,739 MB in GeoPandas. On a 7,342-point join GeoPandas was faster.

Why does my spatial join return zero rows?

Usually a CRS mismatch, which DuckDB does not report as an error. Check the geometry column types on both sides before believing an empty result.

Do I need a spatial index for the join?

No — the join uses a dedicated operator, not the R-tree. An index helps a bounding-box filter against a fixed geometry, not a table-to-table join.

Why does the join return more rows than my input?

Because a point on a shared boundary matches two polygons, and a line crossing several areas matches each. Aggregate if you need one row per input feature.

How do I find the points that matched nothing?

Use a left join and count the null side. In the measured US join, 42,691 of 2,241,395 points — 1.9% — fell outside every state polygon.