DuckDB Spatial Explained: A Spatial Database That Is Just a File
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: 61.8 seconds and 305 MB of memory in DuckDB, against 117.5 seconds and 4,739 MB in GeoPandas โ identical results.
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;")
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()
No import step, no table creation, no server. The query reads the files, and the result comes back as a pandas DataFrame.
Step-by-step solution
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 13.4 ms, against 0.058 ms in SQLite. That is a factor of 230, 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 connection
install is a one-off download; load runs per connection. In an environment with no outbound network the install fails with IO Error: Extension "..." not found, 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.132 s against GeoPandas' 0.276 s; on a 4,596-feature province file, 0.162 s against 0.210 s.
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,552 m instead of the true 343,530 m โ 17.5% too far.
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
from shapely import wkb
df = con.execute("select name, st_aswkb(geom) as geom from prov").df()
gdf = gpd.GeoDataFrame(df, geometry=df["geom"].apply(wkb.loads), crs="EPSG:4326")
The boundary between the two is WKB, and it is cheap.
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, GeoPandas beat DuckDB on a 7,342 ร 4,596 join: 0.104 s against 0.168 s. Below a few hundred thousand rows the setup costs more than it saves.
- Anything needing the wider Python geospatial ecosystem in the middle of the query.
Code examples
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()
Measured on 2,241,395 US points against 51 state polygons: 4.27 s and 227 MB peak memory, against 15.5 s and 988 MB for the equivalent gpd.sjoin. Both produced 2,198,704 matched pairs and the same top three states.
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
Example 3 โ a like-for-like comparison on your own data
import time
import resource
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
duck_rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
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
gpd_rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
print(f"duckdb {duck_result:>12,} pairs {duck_secs:6.2f} s {duck_rss:7,.0f} MB")
print(f"geopandas {gpd_result:>12,} pairs {gpd_secs:6.2f} s {gpd_rss:7,.0f} MB")
assert duck_result == gpd_result, "the two engines disagree โ investigate before trusting either"
The assertion is the important line. A faster answer that differs from the slower one is not a speed-up.
Explanation
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: computing a grid aggregation from a 1.79 GB TSV took 1.01 s; from a 446 MB Parquet file, 0.06 s โ seventeen times faster. Counting distinct values in one column: 1.30 s against 0.03 s.
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 under 2ร. The memory difference was 15.5ร: 305 MB against 4,739 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.
Edge cases or notes
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: 61.8 s and 305 MB against 117.5 s and 4,739 MB. On a small join โ 7,342 points against 4,596 polygons โ GeoPandas was faster.
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 305 MB, and the same query completed with memory_limit set to 200 MB because DuckDB streams and spills.
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).