DuckDB Spatial Explained: A Spatial Database That Is Just a File
On this page
Problem statement
There is a gap in the Python GIS toolbox between a GeoDataFrame and PostGIS.
Below the gap, GeoPandas is ideal: a few hundred thousand features, everything in memory, and the whole pandas API available. Above it, PostGIS is ideal: a server, indexes, concurrent writers, and datasets that never fit in memory.
In between sits a very common case — a few million rows, one analyst, one question, no server. GeoPandas runs out of memory or takes minutes; standing up PostGIS to answer one question is disproportionate.
DuckDB fills that gap. It is an analytical database that runs inside your Python process, reads Parquet, shapefiles, GeoJSON and remote URLs directly, and executes spatial SQL. Measured on a real 13.5-million-point join against 4,596 polygons: 50.9 seconds and 291 MB of memory in DuckDB, against 121.8 seconds and 23,813 MB in GeoPandas — identical results (the benchmark).
Quick answer
Install it, load the spatial extension, and query files where they lie:
import duckdb
con = duckdb.connect() # in memory; or a path for a file
con.execute("install spatial; load spatial;")
top = con.execute("""
select v.name, count(*) as places
from read_parquet('geonames.parquet') g
join st_read('provinces.shp') v
on st_intersects(st_point(g.lon, g.lat), v.geom)
group by 1 order by 2 desc limit 10
""").df()
print(top.head(3).to_string(index=False)) name places
Uttar Pradesh 172436
California 126373
Texas 113915
No import step, no table creation, no server. The query reads the files, and the result comes back as a pandas DataFrame.
What DuckDB is, then install, read, query and hand over
1. Understand what DuckDB is and is not
It is an in-process OLAP engine: columnar storage, vectorised execution, and a query optimiser built for scans and aggregations over many rows.
It is not a transactional store. It has no concurrent-writer story comparable with Postgres, and it is poor at single-row lookups — measured, a point lookup against an indexed 5.2-million-row table took 0.46 ms, against 0.0093 ms in SQLite (measured here). That is a factor of 49, and it is a property of the design rather than a defect.
The distinction that matters: DuckDB is for questions about many rows at once.
2. Install it, and the spatial extension
import duckdb
con = duckdb.connect()
con.execute("install spatial") # downloads once, into ~/.duckdb
con.execute("load spatial") # every new connectioninstall is a one-off download; load runs per connection. In an environment that cannot reach the extension repository the install fails with IO Error: Failed to download extension "spatial" at URL … (Example 4), which is the usual cause of a spatial extension that will not load in CI or in a container.
3. Read the files you already have
The spatial extension wraps GDAL, so ST_Read opens anything GDAL can:
select count(*) from st_read('provinces.shp');
select count(*) from st_read('data.gpkg', layer='roads');
select count(*) from st_read('boundaries.geojson');and DuckDB's own readers handle the columnar formats:
select count(*) from read_parquet('geonames.parquet');
select count(*) from read_csv('points.csv');Measured on Natural Earth: ST_Read on a 7,342-feature shapefile took 0.073 s against GeoPandas' 0.178 s; on a 4,596-feature province file, 0.052 s against 0.143 s (Example 4).
4. Write spatial SQL
The function names follow the PostGIS convention, so most PostGIS SQL transfers:
select st_area(geom), st_length(geom), st_centroid(geom),
st_buffer(geom, 1000), st_intersection(a.geom, b.geom)Two important departures from PostGIS, both measured:
ST_Transformneeds explicit CRS strings and drops the CRS from the result type.ST_Distance_Sphereand the spheroid functions expectPOINT(latitude longitude), notPOINT(x y). With a normal lon/lat point, London to Paris comes back as 403,584 m instead of 343,556 m — 17.5% too far (Example 4).
5. Hand results to GeoPandas when you need it
DuckDB is the engine; GeoPandas is still the right place for plotting, .explore(), and the parts of the ecosystem built on GeoSeries:
import geopandas as gpd
con.execute("create or replace table prov as select name, geom from st_read('provinces.shp')")
df = con.execute("select name, st_aswkb(geom) as geom from prov").df()
gdf = gpd.GeoDataFrame(df.drop(columns="geom"), crs="EPSG:4326",
geometry=gpd.GeoSeries.from_wkb(df["geom"].map(bytes)))
print(len(gdf), gdf.geom_type.value_counts().to_dict())The boundary between the two is WKB, and it is cheap. One detail: DuckDB returns a BLOB column to pandas as bytearray, which shapely.wkb.loads rejects with TypeError: Expected bytes or string; map(bytes) converts it.
6. Know when it is the wrong tool
- Single-row lookups — SQLite or a dictionary, by two orders of magnitude.
- Concurrent writers — Postgres. DuckDB is single-writer.
- Small data — measured, DuckDB still won a 7,342 × 4,596 join, 0.26 s against 0.44 s, but by a fraction of a second. Below a few hundred thousand rows, choose by convenience, not speed.
- Anything needing the wider Python geospatial ecosystem in the middle of the query.
Replacing a GeoPandas pipeline over every GeoNames place
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
For the examples the setup also cuts the United States out of both: its 2.24 million GeoNames points and its 51 state polygons (with the District of Columbia), and Natural Earth's 7,342 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("populated_places.shp"):
with zipfile.ZipFile("populated_places.zip") as z:
for member in z.namelist():
if member.startswith("ne_10m_populated_places."):
open("populated_places" + os.path.splitext(member)[1], "wb").write(z.read(member))
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.gpkg"):
provinces = gpd.read_file("provinces.shp")
provinces[provinces.adm0_a3 == "USA"][["name", "geometry"]].to_file("us_states.gpkg")
print(duckdb.sql("select count(*) from read_parquet('us_points.parquet')").fetchone()[0], "US points,",
len(gpd.read_file("us_states.gpkg")), "US states")2241404 US points, 51 US states
Example 1 — the pattern that replaces a GeoPandas pipeline
import duckdb
def summarise_points_by_area(points_parquet, areas_path, area_name="name"):
"""Count points per polygon without loading either into memory."""
con = duckdb.connect()
con.execute("load spatial")
con.execute("set enable_progress_bar = false")
return con.execute(f"""
select a.{area_name} as area,
count(*) as n,
round(avg(g.population), 1) as mean_population
from read_parquet('{points_parquet}') g
join st_read('{areas_path}') a
on st_intersects(st_point(g.lon, g.lat), a.geom)
where g.lat is not null
group by 1
order by n desc
""").df()
import time
t = time.perf_counter()
states = summarise_points_by_area("us_points.parquet", "us_states.gpkg")
print(f"{len(states)} states, {states.n.sum():,} matched points in {time.perf_counter() - t:.2f} s")
print(states.head(3).to_string(index=False))51 states, 2,198,712 matched points in 3.63 s
area n mean_population
California 126298 2135.3
Texas 113843 1766.8
New York 77522 3061.9
Measured on 2,241,404 US points against 51 state polygons: 3.6 s (Example 1), against 13.6 s for the equivalent gpd.sjoin (Example 3), and 215 MB against 932 MB of peak memory when each runs in its own process (benchmark). Both produced 2,198,712 matched pairs.
Example 2 — a connection helper that sets the things you always want
import duckdb
def spatial_connection(path=":memory:", memory_limit=None, threads=None,
temp_dir=None):
"""One place for the settings every spatial session wants."""
con = duckdb.connect(path)
con.execute("install spatial; load spatial;")
con.execute("set enable_progress_bar = false")
if memory_limit:
con.execute(f"set memory_limit = '{memory_limit}'")
if threads:
con.execute(f"set threads = {threads}")
if temp_dir:
con.execute(f"set temp_directory = '{temp_dir}'")
version = con.execute(
"select extension_version from duckdb_extensions() "
"where extension_name = 'spatial'").fetchone()
print(f"duckdb {duckdb.__version__}, spatial {version[0] if version else '?'}")
return con
con = spatial_connection(memory_limit="2GB", threads=4)duckdb 1.5.5, spatial eb1e57c
Example 3 — a like-for-like comparison on your own data
import time
import duckdb
import geopandas as gpd
def compare_engines(points_parquet, polygons_path):
"""The only benchmark that matters is the one on your data."""
con = duckdb.connect()
con.execute("load spatial")
con.execute("set enable_progress_bar = false")
start = time.perf_counter()
duck_result = con.execute(f"""
select count(*) from read_parquet('{points_parquet}') p
join st_read('{polygons_path}') a
on st_intersects(st_point(p.lon, p.lat), a.geom)
""").fetchone()[0]
duck_secs = time.perf_counter() - start
start = time.perf_counter()
polygons = gpd.read_file(polygons_path)
frame = con.execute(f"select lon, lat from read_parquet('{points_parquet}')").df()
points = gpd.GeoDataFrame(frame,
geometry=gpd.points_from_xy(frame.lon, frame.lat),
crs=polygons.crs)
gpd_result = len(gpd.sjoin(points, polygons, predicate="intersects"))
gpd_secs = time.perf_counter() - start
print(f"duckdb {duck_result:>12,} pairs {duck_secs:6.2f} s")
print(f"geopandas {gpd_result:>12,} pairs {gpd_secs:6.2f} s")
assert duck_result == gpd_result, "the two engines disagree — investigate before trusting either"
compare_engines("us_points.parquet", "us_states.gpkg")duckdb 2,198,712 pairs 3.58 s
geopandas 2,198,712 pairs 13.57 s
Peak memory is deliberately left out: a process's peak resident size only ever grows, so measured in one process the second engine inherits the first one's peak. The benchmark guide runs each engine in its own process for that reason.
The assertion is the important line. A faster answer that differs from the slower one is not a speed-up.
Example 4 — the other measurements on this page
import shutil
for label, path in (("7,342 populated places", "populated_places.shp"), ("4,596 provinces", "provinces.shp")):
runs = []
for _ in range(3):
t0 = time.perf_counter(); con.execute(f"select count(*) from st_read('{path}')").fetchone()
runs.append(time.perf_counter() - t0)
gruns = []
for _ in range(3):
t0 = time.perf_counter(); gpd.read_file(path)
gruns.append(time.perf_counter() - t0)
print(f"{label:24} ST_Read {min(runs):.3f} s GeoPandas read_file {min(gruns):.3f} s")
london, paris = (51.5074, -0.1278), (48.8566, 2.3522)
right, swapped = con.execute(f"""select
st_distance_sphere(st_point({london[0]}, {london[1]}), st_point({paris[0]}, {paris[1]})),
st_distance_sphere(st_point({london[1]}, {london[0]}), st_point({paris[1]}, {paris[0]}))""").fetchone()
print(f"London–Paris, POINT(lat lon): {right:,.0f} m; POINT(lon lat): {swapped:,.0f} m ({swapped / right - 1:+.1%})")
offline = duckdb.connect(config={"extension_directory": "empty_extensions"})
try:
offline.execute("set custom_extension_repository = 'http://127.0.0.1:9'")
offline.execute("install spatial")
except duckdb.Error as error:
print(f"install with no reachable repository: {type(error).__name__}: {str(error).splitlines()[0][:90]}")
shutil.rmtree("empty_extensions", ignore_errors=True)7,342 populated places ST_Read 0.073 s GeoPandas read_file 0.178 s
4,596 provinces ST_Read 0.052 s GeoPandas read_file 0.143 s
London–Paris, POINT(lat lon): 343,556 m; POINT(lon lat): 403,584 m (+17.5%)
install with no reachable repository: IOException: IO Error: Failed to download extension "spatial" at URL "http://127.0.0.1:9/v1.5.5/linux_a
An in-process, columnar database
Why an in-process database changes the ergonomics
PostGIS is a service: it has to be installed, started, secured, loaded and kept in sync with the files on disk. That overhead is worth paying when many people share the data or when it never fits in memory.
DuckDB has no overhead at all. pip install duckdb, and the database is the file you already have. Nothing is loaded, so nothing is stale, and the query reads the current file every time.
That is what makes it fit the analyst's loop: ask a question, get an answer, change the question. The absence of an import step is not a convenience — it removes an entire category of "the database is out of date" bugs.
Why columnar storage wins on analytical queries
DuckDB stores columns together and processes them in vectors. A query touching three columns of a nineteen-column table reads three columns.
Measured on the same 13.5-million-row dataset in two formats (here): a grid aggregation took 0.42 s from a 1.79 GB TSV and 0.07 s from a 446 MB Parquet file — six times faster. Counting distinct values in one column: 0.56 s against 0.01 s. The engine was the same; the layout decided how much it had to read.
The row-oriented format has to parse every field of every row to reach one column. The columnar one reads only what the query names.
Why the memory numbers are the more important result
The timing difference on the global join was 2.4×. The memory difference was 82×: 291 MB against 23,813 MB for the same 13.5-million-point join.
That is the qualitative change. A machine with 8 GB of RAM can run the DuckDB version and cannot run the GeoPandas one, because DuckDB streams the join rather than materialising both sides. Below the memory ceiling the speed difference is a convenience; at the ceiling it is the difference between an answer and a MemoryError.
Why it is not a PostGIS replacement
DuckDB has one writer, no roles, no row-level security, no logical replication and no connection pooling. Its spatial join is a dedicated operator rather than an index lookup, which is excellent for a full join and useless for "give me this one feature by id" in a web request.
PostGIS is a system of record. DuckDB is an analysis engine you point at files. Projects frequently use both: PostGIS holds the data, and DuckDB reads a Parquet export of it to answer analytical questions without loading the server.
Per-connection loading and lat/lon argument order
load spatialis per connection;installis once per machine.- The extension is a separate download — an offline container needs it baked in.
ST_Readgoes through GDAL, so its format support is GDAL's, including its quirks.ST_Distance_SpheretakesPOINT(lat lon). With lon/lat input the London–Paris distance came back 17.5% too long.ST_Transformneeds explicit source and target CRS and returns a geometry with no CRS attached.- Set
enable_progress_bar = falsein scripts, or every query prints a progress bar. .df()needs pandas;.arrow()avoids the copy for large results.- It is single-writer. Two processes writing one database file is not supported.
Internal links
- How to install and load the DuckDB spatial extension — the setup, and what breaks
- DuckDB or PostGIS: which one a workflow needs — choosing between them
- Columnar or row storage: why DuckDB is fast on wide tables — the storage model
- How to read shapefiles, GeoJSON and GeoParquet in DuckDB — the readers
- How to run a spatial join in DuckDB — the operation this article measured
- How to move data between DuckDB and GeoPandas — the WKB boundary
- CRS in DuckDB: why ST_Transform moves your data to the wrong place — the axis-order trap
- Spatial SQL versus GeoPandas — the wider comparison
FAQ
What is DuckDB spatial?
An extension that adds geometry types, spatial SQL functions and GDAL-backed file readers to DuckDB, an analytical database that runs inside your Python process with no server.
Is DuckDB faster than GeoPandas?
For large operations, yes. Measured on a 13.5-million-point join against 4,596 polygons: 50.9 s and 291 MB against 121.8 s and 23,813 MB. On a small join — 7,342 points against 4,596 polygons — DuckDB was still faster, 0.26 s against 0.44 s, but by an amount nobody waits for.
Can DuckDB replace PostGIS?
No. It is single-writer, has no roles or security model, and is two orders of magnitude slower at single-row lookups. It complements PostGIS as an analysis engine over files.
Do I need to import my data first?
No. DuckDB queries Parquet, CSV, shapefiles, GeoPackages and GeoJSON in place, which removes the whole class of "the database copy is stale" problems.
How much memory does it need?
Far less than an in-memory library. The 13.5-million-point join peaked at 291 MB, and a 2.2-million-point join finished inside a 200 MB memory_limit at no cost in time, because DuckDB streams the aggregate (measured here).
What is the biggest gotcha?
CRS handling. ST_Transform requires explicit CRS strings and drops the CRS from its result, and the sphere and spheroid distance functions expect POINT(latitude longitude).