How to Use DuckDB as the Engine in a GIS Pipeline

On this page

Problem statement

A GIS pipeline usually accretes: read a shapefile, filter it in pandas, join it, write a GeoPackage, read that back, aggregate it, write a CSV. Each stage materialises everything, each intermediate file has to be managed, and the whole thing runs at the speed of the slowest materialisation.

Using DuckDB as the pipeline's engine changes the structure rather than just the speed. The transformations become SQL views over the source files, nothing is materialised until the end, and the intermediate files disappear.

The scale difference is measurable: a 13,472,195-point join and aggregate ran in 50.9 s using 291 MB, against 121.8 s and 23,813 MB for the equivalent GeoPandas code (benchmark). But the more important change is that there is one artefact โ€” the query โ€” instead of six intermediate files.

Quick answer

Express the pipeline as chained SQL, materialise once:

import duckdb


def run_pipeline(config):
    con = duckdb.connect(config["database"])
    con.execute("install spatial; load spatial;")
    con.execute("set enable_progress_bar = false")
    con.execute(f"set memory_limit = '{config['memory_limit']}'")
    con.execute(f"set temp_directory = '{config['temp_dir']}'")

    con.execute(f"""
        create or replace view raw as
        select * from read_parquet('{config["input"]}')
    """)
    con.execute("""
        create or replace view cleaned as
        select * from raw
        where lat between -90 and 90 and lon between -180 and 180
          and lat is not null
    """)
    con.execute(f"""
        create or replace view joined as
        select c.*, a.name as area
        from cleaned c
        join st_read('{config["areas"]}') a
          on st_intersects(st_point(c.lon, c.lat), a.geom)
    """)
    con.execute(f"""
        copy (select area, count(*) as n from joined group by 1 order by n desc)
        to '{config["output"]}' (format parquet)
    """)


import time

t = time.perf_counter()
run_pipeline({"database": "pipeline.duckdb", "memory_limit": "1GB", "temp_dir": "duckdb_tmp",
              "input": "places.parquet", "areas": "provinces.shp", "output": "per_area.parquet"})
summary = duckdb.sql("select * from read_parquet('per_area.parquet')").df()
print(f"{len(summary):,} areas written in {time.perf_counter() - t:.1f} s; busiest:")
print(summary.head(3).to_string(index=False))
4,419 areas written in 52.7 s; busiest:
         area      n
Uttar Pradesh 172436
   California 126373
        Texas 113915

Views cost nothing โ€” they are query fragments. Only the final COPY touches the disk.

Chained views compose into one plan; materialised stages do not. six stages, two architecturesintermediate filesโœ•six writes and six readsโœ•six files to name and clean upโœ•six chances to go staleโœ•the optimiser sees one stage at a timechained viewsโœ“nothing written until the endโœ“one artefact: the queryโœ“filters push into the file readersโœ“the optimiser sees the whole planThe intermediate files were never the pipeline โ€” they were an artefact of tools that could not compose spatialworkflow.io
Materialise a stage only when it is expensive and reused several times.

Views as stages, SQL in files, validation between them

1. Make the stages views, not tables

A view is a named query. Chaining views composes into one plan that the optimiser sees whole, so filters push down into the file readers and columns nobody uses are never read.

Materialising each stage as a table does the opposite: it forces every intermediate result into storage and hides the optimisation opportunities.

Use a real table only when a stage is expensive and reused several times downstream.

2. Keep the SQL in files, not in strings

Pipeline SQL grows. Keeping it in .sql files makes it reviewable, diffable, syntax-highlighted and testable:

pipeline/
  01_raw.sql
  02_cleaned.sql
  03_joined.sql
  04_summary.sql

Then the Python is a thin runner that reads, parameterises and executes them in order โ€” which is also what makes the pipeline easy to run partially while developing.

3. Parameterise properly

String-formatting a path into SQL is fine for a trusted config and wrong for anything user-supplied. DuckDB supports parameters:

con = duckdb.connect()
path, country = "places.parquet", "IS"
print(con.execute("select count(*) from read_parquet(?) where country = ?", [path, country]).fetchone())
(16296,)

For paths that must be interpolated โ€” table names, file globs โ€” validate them against an allowlist rather than trusting them.

4. Validate between stages

The advantage of a database engine is that assertions are queries:

def assert_rows(con, relation, minimum=1, maximum=None):
    n = con.execute(f"select count(*) from {relation}").fetchone()[0]
    if n < minimum or (maximum and n > maximum):
        raise ValueError(f"{relation}: {n:,} rows, expected {minimum:,}โ€“{maximum or 'โˆž'}")
    return n


print(assert_rows(con, "read_parquet('per_area.parquet')", minimum=100, maximum=5000))
try:
    assert_rows(con, "read_parquet('per_area.parquet')", minimum=5000)
except ValueError as error:
    print("ValueError:", error)
4419
ValueError: read_parquet('per_area.parquet'): 4,419 rows, expected 5,000โ€“โˆž

Cheap checks worth running after every stage: row count within an expected range, no nulls in key columns, no duplicate keys, and โ€” for spatial stages โ€” the count of features that matched nothing.

5. Record what ran

A pipeline that cannot say what it did is not reproducible. Record the inputs, their sizes and modification times, the DuckDB and extension versions, the row counts at each stage, and the wall-clock time:

con.execute("""
    create table if not exists _run_log (
        run_id varchar, stage varchar, rows bigint,
        seconds double, logged_at timestamp)
""")

Writing the log into the same database file keeps the provenance next to the output. Avoid at as a column name: it is a reserved word in DuckDB 1.5, and seconds double, at timestamp fails with Parser Error: syntax error at or near "at".

6. Make it idempotent

create or replace view and create or replace table are idempotent by construction. COPY ... TO overwrites. That means a rerun produces the same result and no duplicates, which is the property that makes a pipeline safe to retry.

The one thing to watch: a COPY ... TO into a partitioned directory needs overwrite_or_ignore, or a rerun fails on the existing directory.

A DuckDB pipeline as numbered SQL files run by a thin Python runner. SQL in files, Python as the runner101_raw.sqlread the files in place202_cleaned.sqlfilters as a view303_joined.sqlthe spatial join404_summary.sqlCOPY โ€” the only write5runner.pylog, validate, retryA 60-line query in a Python string is invisible to every tool that could help you spatialworkflow.io
The runner stops changing once it works; the SQL is where the pipeline lives.

A numbered-SQL pipeline with provenance

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

Example 1 โ€” a runner that executes numbered SQL files

import glob
import os
import time
import uuid
import duckdb


class SqlPipeline:
    """Runs numbered .sql files in order, with logging and validation."""

    def __init__(self, sql_dir, database=":memory:", settings=None, params=None):
        self.sql_dir = sql_dir
        self.params = params or {}
        self.run_id = uuid.uuid4().hex[:8]
        self.con = duckdb.connect(database)
        self.con.execute("install spatial; load spatial;")
        self.con.execute("set enable_progress_bar = false")
        for key, value in (settings or {}).items():
            self.con.execute(f"set {key} = '{value}'")
        self.con.execute("""
            create table if not exists _run_log (
                run_id varchar, stage varchar, rows bigint,
                seconds double, logged_at timestamp)
        """)

    def run(self, upto=None):
        files = sorted(glob.glob(os.path.join(self.sql_dir, "*.sql")))
        for path in files:
            stage = os.path.basename(path)
            if upto and stage > upto:
                break
            sql = open(path, encoding="utf-8").read().format(**self.params)

            started = time.perf_counter()
            self.con.execute(sql)
            elapsed = time.perf_counter() - started

            rows = self._rows_of(stage)
            self.con.execute(
                "insert into _run_log values (?, ?, ?, ?, now())",
                [self.run_id, stage, rows, elapsed])
            print(f"{stage:24} {rows if rows is not None else '-':>12}  {elapsed:6.2f}s")
        return self

    def _rows_of(self, stage):
        relation = os.path.splitext(stage)[0].split("_", 1)[-1]
        try:
            return self.con.execute(f"select count(*) from {relation}").fetchone()[0]
        except Exception:
            return None


os.makedirs("pipeline", exist_ok=True)
STAGES = {
    "01_raw.sql": "create or replace view raw as select * from read_parquet('{input}')",
    "02_cleaned.sql": ("create or replace view cleaned as select * from raw "
                       "where lat between -90 and 90 and lon between -180 and 180 and country = '{country}'"),
    "03_joined.sql": ("create or replace view joined as select c.*, a.name as area from cleaned c "
                      "left join st_read('{areas}') a on st_intersects(st_point(c.lon, c.lat), a.geom)"),
    "04_summary.sql": ("create or replace table summary as select area, count(*) as n from joined "
                       "group by 1 order by n desc"),
}
for name, sql in STAGES.items():
    open(os.path.join("pipeline", name), "w").write(sql)
pipe = SqlPipeline("pipeline", params={"input": "places.parquet", "areas": "provinces.shp", "country": "GB"}).run()
01_raw.sql                   13472195    0.00s
02_cleaned.sql                 109179    0.00s
03_joined.sql                  109179    0.00s
04_summary.sql                    247    4.13s

Example 2 โ€” validation as SQL

CHECKS = {
    "no null geometry": "select count(*) from {r} where geom is null",
    "no duplicate ids": ("select count(*) from (select id from {r} "
                         "group by 1 having count(*) > 1)"),
    "coordinates in range": ("select count(*) from {r} "
                             "where lat not between -90 and 90 "
                             "or lon not between -180 and 180"),
    "all rows matched an area": "select count(*) from {r} where area is null",
}


def validate(con, relation, checks=CHECKS, tolerate=None):
    """Every check returns a count of bad rows; zero is a pass."""
    tolerate = tolerate or {}
    failures = []
    for name, template in checks.items():
        bad = con.execute(template.format(r=relation)).fetchone()[0]
        allowed = tolerate.get(name, 0)
        status = "ok" if bad <= allowed else "FAIL"
        print(f"  {name:26} {bad:8,}  {status}")
        if bad > allowed:
            failures.append(f"{name}: {bad:,} rows (allowed {allowed:,})")
    if failures:
        raise ValueError(f"{relation} failed validation:\n  " + "\n  ".join(failures))


pipe.con.execute("""create or replace view checked as
                    select row_number() over () as id, lat, lon, st_point(lon, lat) as geom, area from joined""")
unmatched = pipe.con.execute("select count(*) from checked where area is null").fetchone()[0]
total = pipe.con.execute("select count(*) from checked").fetchone()[0]
print(f"{unmatched:,} of {total:,} British points outside every province ({unmatched / total:.1%})")
validate(pipe.con, "checked", tolerate={"all rows matched an area": unmatched})
8,684 of 109,179 British points outside every province (8.0%)
  no null geometry                  0  ok
  no duplicate ids                  0  ok
  coordinates in range              0  ok
  all rows matched an area      8,684  ok

The tolerate argument is what makes this usable on real data. A spatial join legitimately leaves some rows unmatched โ€” measured, 3.9% of the global point set (here) and 8.0% of the British points above fell outside every province polygon โ€” so the check should record the number and fail only when it moves.

Example 3 โ€” the pipeline's provenance record

def write_provenance(con, inputs, output, run_id):
    """What ran, against what, with which versions."""
    import os
    import datetime
    import duckdb

    rows = []
    for path in inputs:
        stat = os.stat(path)
        rows.append((run_id, path, stat.st_size,
                     datetime.datetime.fromtimestamp(stat.st_mtime).isoformat()))

    con.execute("""
        create table if not exists _inputs (
            run_id varchar, path varchar, bytes bigint, modified varchar)
    """)
    con.executemany("insert into _inputs values (?, ?, ?, ?)", rows)

    versions = con.execute("""
        select extension_name, extension_version from duckdb_extensions()
        where loaded
    """).fetchall()
    print(f"run {run_id}: duckdb {duckdb.__version__}, "
          f"{', '.join(f'{n} {v}' for n, v in versions)}")
    print(f"  {len(inputs)} inputs โ†’ {output}")


write_provenance(pipe.con, ["places.parquet", "provinces.shp"], "summary", pipe.run_id)
print(pipe.con.execute("select path, bytes from _inputs").fetchall())
run 92686a61: duckdb 1.5.5, core_functions v1.5.5, icu v1.5.5, json v1.5.5, parquet v1.5.5, spatial eb1e57c
  2 inputs โ†’ summary
[('places.parquet', 306039192), ('provinces.shp', 20998780)]

Example 4 โ€” rerunning a partitioned COPY

import shutil

shutil.rmtree("by_country", ignore_errors=True)
PART = ("copy (select country, name, lat, lon from read_parquet('places.parquet') where country in ('IS', 'MT')) "
        "to 'by_country' (format parquet, partition_by (country){extra})")
con.execute(PART.format(extra=""))
print("first run:", sorted(os.listdir("by_country")))
try:
    con.execute(PART.format(extra=""))
except duckdb.Error as error:
    print(f"second run: {type(error).__name__}: {str(error).splitlines()[0][:100]}")
con.execute(PART.format(extra=", overwrite_or_ignore"))
print("with overwrite_or_ignore: ok,", sorted(os.listdir("by_country")))
first run: ['country=IS', 'country=MT']
second run: IOException: IO Error: Directory "by_country" is not empty! Enable OVERWRITE option to overwrite files
with overwrite_or_ignore: ok, ['country=IS', 'country=MT']

Views over intermediate files

Why views beat intermediate files

Each materialised intermediate is a write, a read, a file to name, a file to clean up and a file that can go stale. Six stages means six of each.

Views collapse the chain into one plan. The optimiser can then push the final filter all the way down into the Parquet reader, so a query whose output is a few thousand rows may read a fraction of the input columns โ€” measured, as little as 0.06% of a file for a selective query.

The intermediate files were never the pipeline; they were an artefact of using tools that could not compose.

Why SQL in files is better than SQL in strings

A 60-line query inside a Python triple-quoted string is invisible to every tool: no syntax highlighting, no linting, awkward diffs, and no way to run it by hand while debugging.

The same query in a .sql file can be opened in any client, run against the database interactively, reviewed properly in a pull request, and formatted. The Python that runs it becomes a dozen lines that do not change.

Why validation belongs between stages

A pipeline that validates only its output tells you that something is wrong. A pipeline that validates between stages tells you where.

The checks are cheap because they are queries against data the engine already has open. Counting nulls in a key column of a 13-million-row view is a single scan of one column โ€” measured at hundredths of a second on a columnar file.

Why idempotence is the property that makes it operable

Reruns happen: a stage fails, a machine restarts, somebody re-runs to check. If a rerun duplicates rows or fails on an existing file, every rerun becomes a manual cleanup.

create or replace and COPY ... TO with overwrite make the pipeline safe to run twice. That single property is what allows retries, and retries are what allow the pipeline to run unattended.

Cheap validation queries to run after each pipeline stage. assertions are just queriesrow count inside an expected rangecatches a silent empty readno nulls in key columnsone scan of one columnno duplicate keysa group-by with a havingunmatched rows within tolerance3.9% was normal in one joinvalidating only the outputtells you that, not whereA tolerance is required: real spatial joins always leave some rows unmatched spatialworkflow.io
The checks are cheap because the engine already has the data open.

Re-evaluated views and column type changes

  • Views are re-evaluated on every reference. If a stage is expensive and used three times, materialise it as a table.
  • create or replace view cannot change a view's column types while something depends on it; drop dependents first.
  • COPY ... TO a partitioned directory needs overwrite_or_ignore.
  • Set memory_limit and temp_directory in the runner, not in each SQL file.
  • A DuckDB database file is a single-writer artefact โ€” two pipeline runs must not share one.
  • Keep the run log in the database so provenance travels with the output.
  • explain a slow stage before optimising it; the plan usually names the problem.
  • Pin the DuckDB version. Extension binaries are version-specific and plans can change.

FAQ

Should each pipeline stage be a table or a view?

A view, unless the stage is expensive and referenced several times. Views compose into one plan the optimiser can see whole; tables force materialisation at every step.

Where should the SQL live?

In numbered .sql files. They are reviewable, runnable by hand and diffable, and the Python runner stops changing once it works.

How do I validate between stages?

With queries returning a count of bad rows โ€” nulls, duplicates, out-of-range coordinates, unmatched joins. Allow a tolerance where the data legitimately has one.

Is a DuckDB pipeline faster than a GeoPandas one?

At scale, and with much less memory. A 13.5-million-point join and aggregate ran in 50.9 s and 291 MB against 121.8 s and 23,813 MB.

How do I make it safe to rerun?

create or replace for views and tables, overwriting COPY ... TO for outputs, and overwrite_or_ignore for partitioned writes. Then a rerun is a no-op rather than a cleanup.

Can two pipeline runs share one database file?

No. DuckDB is single-writer. Give each run its own file, or use an in-memory database and write only the outputs.