How to Run a Spatial Join in DuckDB
On this page
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 a spatial join is one line that holds both datasets in memory; in DuckDB it is a SQL JOIN with a spatial predicate, streaming one table through a structure built over the other. Joining 13,472,195 GeoNames points to 4,596 provinces measures the difference (the benchmark):
engine time peak memory result
DuckDB 50.9 s 291 MB 12,950,244 pairs
GeoPandas 121.8 s 23,813 MB 12,950,244 pairs
Same answer, 2.4× the speed and 82× less memory. The memory figure is the one that decides whether the job runs at all on an ordinary machine.
DuckDB is not automatically worth the switch, though: on a small join — 7,342 points against the same 4,596 polygons — it was faster by less than a quarter of a second, 0.26 s against 0.44 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")
query = """select p.name, a.name as area
from read_parquet('places.parquet') p
join st_read('provinces.shp') a
on st_intersects(st_point(p.lon, p.lat), a.geom)"""
result = con.execute(query).df()
print(f"{len(result):,} point–area pairs")12,950,244 point–area pairs
Pick the predicate, check the CRS, aggregate in SQL
1. Choose the predicate deliberately
DuckDB's ST_ predicates agree exactly with GEOS, so a point on a shared state boundary intersects both states and is within neither, just as in GeoPandas. Which predicate asks your question is spatial predicates explained; here, ST_Intersects for points in states, and ST_DWithin in Example 2.
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 allA 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 shape ran in 3.66 s with 215 MB peak memory on 2,241,404 points against 51 polygons, against 13.83 s and 932 MB for the GeoPandas equivalent (benchmark).
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 (Example 1), 2,241,404 points produced 2,198,712 matches — 42,692 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. Aggregate before you trust a row count
The joined table can hold more rows than either input (spatial join cardinality); in SQL the cure is step 4's GROUP BY on the feature you want one row per.
A trustworthy join on every GeoNames point
The data used on this page
The points are GeoNames' full allCountries dump — every feature, about 13.5 million, rebuilt daily at download.geonames.org, so counts drift a little between runs — kept both as the 1.79 GB text file and as places.parquet (also linked as geonames.parquet) with name, latitude, longitude, feature class, country and population. The polygons are Natural Earth's 1:10m admin-1 states and provinces, 4,596 of them, as provinces.shp.
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)
fetch("https://download.geonames.org/export/dump/allCountries.zip", "allCountries.zip") # about 420 MB
if not os.path.exists("allCountries.txt"):
zipfile.ZipFile("allCountries.zip").extract("allCountries.txt")
if not os.path.exists("places.parquet"):
duckdb.sql("""
copy (select column01 as name, column04 as lat, column05 as lon, column06 as feature_class,
column08 as country, column14 as population
from read_csv('allCountries.txt', delim = '\\t', header = false, quote = '',
all_varchar = false, sample_size = -1))
to 'places.parquet' (format parquet)""")
if not os.path.exists("geonames.parquet"):
os.link("places.parquet", "geonames.parquet")
fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip", "provinces.zip")
if not os.path.exists("provinces.shp"):
with zipfile.ZipFile("provinces.zip") as z:
for member in z.namelist():
if member.startswith("ne_10m_admin_1_states_provinces."):
open("provinces" + os.path.splitext(member)[1], "wb").write(z.read(member))
n = duckdb.sql("select count(*) from read_parquet('places.parquet')").fetchone()[0]
print(f"{n:,} GeoNames points; places.parquet {os.path.getsize('places.parquet') / 1e6:.0f} MB, "
f"allCountries.txt {os.path.getsize('allCountries.txt') / 1e9:.2f} GB")13,472,195 GeoNames points; places.parquet 306 MB, allCountries.txt 1.79 GB
The examples use the United States cut out of both — its 2.24 million GeoNames points and its 51 state polygons, as us_states.shp — and, for the distance join, Great Britain's GeoNames populated places against the British towns in Natural Earth's populated places.
import geopandas as gpd
fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_populated_places.zip", "populated_places.zip")
if not os.path.exists("us_points.parquet"):
duckdb.sql("copy (select * from read_parquet('places.parquet') where country = 'US') "
"to 'us_points.parquet' (format parquet)")
if not os.path.exists("us_states.shp"):
provinces = gpd.read_file("provinces.shp")
provinces[provinces.adm0_a3 == "USA"][["name", "geometry"]].to_file("us_states.shp")
if not os.path.exists("uk_towns.parquet"):
towns = gpd.read_file("populated_places.zip")
towns = towns[towns.ADM0_A3 == "GBR"]
towns.assign(name=towns.NAME, lon=towns.geometry.x, lat=towns.geometry.y)[["name", "lon", "lat"]] \
.to_parquet("uk_towns.parquet")
print(duckdb.sql("select count(*) from read_parquet('us_points.parquet')").fetchone()[0], "US points;",
duckdb.sql("select count(*) from read_parquet('uk_towns.parquet')").fetchone()[0], "British towns")2241404 US points; 57 British towns
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()
states = spatial_join_summary("read_parquet('us_points.parquet')", "st_read('us_states.shp')")
print(states.head(3).to_string(index=False))polygon geometry type: GEOMETRY('EPSG:4326')
2,241,404 point-rows, 2,198,712 matched, 42,692 unmatched (1.9%)
area n
California 126298
Texas 113843
New York 77522
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()
con = duckdb.connect()
con.execute("load spatial")
con.execute("""create table gb as select name, st_point(lon, lat) as geom
from read_parquet('places.parquet') where country = 'GB' and feature_class = 'P'""")
con.execute("create table towns as select name, st_point(lon, lat) as geom from read_parquet('uk_towns.parquet')")
near = within_distance(con, "gb", "towns", 5000)
print(f"{len(near)} towns; the busiest:")
print(near.head(3).to_string(index=False))57 towns; the busiest:
site within_5000m nearest_m
London 185 208.9
Sheffield 102 469.4
Bradford 96 294.8
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")
cross_check(con, "read_parquet('places.parquet')", "provinces.shp")duckdb 48,539 geopandas 48,539 agree
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.
Example 4 — the silent CRS mismatch, reproduced
con.execute("create or replace table areas as select name, geom from st_read('us_states.shp')")
con.execute("""create or replace table areas_3857 as
select name, st_transform(geom, 'EPSG:4326', 'EPSG:3857', always_xy := true) as geom from areas""")
for table in ("areas", "areas_3857"):
kind = con.execute(f"select typeof(geom) from {table} limit 1").fetchone()[0]
n = con.execute(f"""select count(*) from read_parquet('us_points.parquet') p
join {table} a on st_intersects(st_point(p.lon, p.lat), a.geom)""").fetchone()[0]
print(f"{table:11} {kind:24} {n:,} matches")
print("a bare point:", con.execute("select typeof(st_point(-0.1, 51.5))").fetchone()[0])areas GEOMETRY('EPSG:4326') 2,198,712 matches
areas_3857 GEOMETRY 0 matches
a bare point: GEOMETRY
A join operator, not an index lookup
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 2.4× faster and 82× 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 US join completed with memory_limit set to 200 MB — measured at 3.53 s, no slower than the unconstrained run (here). The engine was never holding the data in the first place.
Why GeoPandas wins on small joins
DuckDB pays fixed costs — planning the query, loading the extension, building the join structure — and on the 7,342-point join they were most of the work.
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.
Invalid geometry
Invalid geometry can produce wrong results rather than an error; check ST_IsValid before joining untrusted data.
Internal links
- DuckDB spatial explained: a spatial database that is just a file — the engine
- How DuckDB uses (and does not use) a spatial index — why the index is not in the plan
- Fixing a DuckDB spatial join that never finishes — when it does not complete
- CRS in DuckDB: why ST_Transform moves your data to the wrong place — the zero-row failure
- How to join points to polygons with a spatial join — the GeoPandas equivalent
- Spatial join cardinality explained — why the row count changes
- How to aggregate millions of points into a grid with DuckDB — joining to a grid instead of polygons
- How to benchmark DuckDB against GeoPandas honestly — making the comparison fairly
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 50.9 s and 291 MB in DuckDB against 121.8 s and 23,813 MB in GeoPandas. On a 7,342-point join the difference was a fifth of a second.
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.
How do I find the points that matched nothing?
Use a left join and count the null side. In the measured US join, 42,692 of 2,241,404 points — 1.9% — fell outside every state polygon.