How to Benchmark DuckDB Against GeoPandas Honestly

On this page

Problem statement

Benchmarks between GIS engines are usually wrong in the same few ways, and all of them favour whichever engine the author preferred:

The file read is inside one timer and outside the other. DuckDB reading Parquet in place looks instant beside GeoPandas' explicit read_file, when in reality the read is happening — just later.
The page cache is warm for the second run. Whichever engine goes second wins.
Only time is measured. The 82× memory difference measured on a real join is invisible in a stopwatch.
The results are never compared. A CRS mismatch in DuckDB returns zero rows without an error, so an "8× speed-up" can be an empty result.

A benchmark that fixes those four produces numbers you can act on — and they do not all say the same thing. Measured on the same machine and the same data (Example 2):

join                                  DuckDB              GeoPandas
7,342 pts × 4,596 polys       0.26 s / 215 MB      0.44 s / 318 MB
2,241,404 pts × 51 polys      3.66 s / 215 MB     13.83 s / 932 MB
13,472,195 pts × 4,596 polys 50.85 s / 291 MB   121.80 s / 23,813 MB

Quick answer

Time the whole task, measure memory, run each engine in its own process, and assert the results match:

import json
import subprocess
import sys


def benchmark(engine, script, *args):
    """Each engine in a clean process: no shared caches, real peak memory."""
    result = subprocess.run(
        [sys.executable, script, engine, *map(str, args)],
        capture_output=True, text=True, check=True)
    return json.loads(result.stdout.strip().splitlines()[-1])


duck = benchmark("duck", "bench.py", points, polygons)
gpd_ = benchmark("gpd", "bench.py", points, polygons)

assert duck["result"] == gpd_["result"], (
    f"results differ: {duck['result']} vs {gpd_['result']} — "
    f"check the CRS on both sides before comparing speed")

print(f"duckdb    {duck['secs']:6.2f}s  {duck['peak_mb']:7,.0f} MB")
print(f"geopandas {gpd_['secs']:6.2f}s  {gpd_['peak_mb']:7,.0f} MB")

Compare the outputs before comparing the times.

What you get — from The data used on this page

13,472,195 GeoNames points; places.parquet 306 MB, allCountries.txt 1.79 GB

Define the task, isolate engines, measure memory, compare results

1. Define the task, not the operation

Benchmark "count points per polygon from these files", not "sjoin". The task includes reading the data, because that is what a real workflow does and because the two engines do it at different points.

DuckDB reads inside the query; GeoPandas reads before it. Timing only the join hides half of GeoPandas' work and all of DuckDB's.

2. Run each engine in its own process

Peak resident memory is process-wide. Running both engines in one process reports the maximum of the two, attributed to whichever ran second.

A subprocess per engine gives a clean ru_maxrss, prevents warm caches inside a library from carrying over, and removes any chance of one import affecting the other.

3. Report memory as well as time

Time differences are often modest; memory differences are frequently decisive. The measured global join was 2.4× faster in DuckDB and used 82× less memory — 291 MB against 23,813 MB.

That second number is what decides whether the job runs on a laptop at all. A benchmark that omits it is measuring the less important axis.

4. Control the cache deliberately

Two policies, both defensible, and the important thing is to state which one you used:

Warm — run once to warm the page cache, then measure. Reflects a repeated interactive workload.
Cold — drop caches between runs (needs root on Linux) or use a file larger than RAM. Reflects a first run.

Also turn off DuckDB's own cache when measuring remote reads: set enable_external_file_cache = false, or the second query reads nothing and appears miraculous.

5. Assert that the results agree

This is the check that catches the embarrassing benchmark. A CRS mismatch in DuckDB returns zero rows with no error, and zero rows is very fast.

Compare the actual outputs: the matched-pair count, the group count, the top few groups. Measured on the global join, both engines returned 12,950,244 pairs across 4,419 groups with the same top three — which is what makes the timing comparison meaningful. On the first run they did not agree on the smallest join, and the cause was the aggregation, not the engines (Example 1).

6. Take the best of several runs, and say how many

Wall-clock timings on a shared machine are noisy. Take the minimum of three to five runs — the minimum is the least contaminated by other load — and report the run count alongside.

Report the machine too: core count, RAM, and whether the data was on SSD or a network mount. Numbers without a machine are not reproducible.

State whether the cache was warm. Both choices are valid; silence is not.

A fair benchmark across data sizes

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

Three point–polygon pairs are prepared from those: Natural Earth's 7,342 populated places against all 4,596 provinces, the 2.24 million GeoNames points in the United States against the 51 US states (with the District of Columbia), and everything against everything.

import geopandas as gpd

fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_populated_places.zip", "populated_places.zip")
if not os.path.exists("towns.parquet"):
    towns = gpd.read_file("populated_places.zip")
    towns.assign(name=towns.NAME, lon=towns.geometry.x, lat=towns.geometry.y)[["name", "lon", "lat"]].to_parquet("towns.parquet")
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")
for name, sql in (("towns", "read_parquet('towns.parquet')"), ("US points", "read_parquet('us_points.parquet')")):
    print(f"{name}: {duckdb.sql(f'select count(*) from {sql}').fetchone()[0]:,}")
print("US states:", len(gpd.read_file("us_states.gpkg")))

towns: 7,342
US points: 2,241,404
US states: 51

Example 1 — the benchmark script, one engine per invocation

%%writefile bench.py
#!/usr/bin/env python
"""bench.py <engine> <points.parquet> <polygons.shp>"""
import json
import resource
import sys
import time


def run_duckdb(points, polygons):
    import duckdb
    con = duckdb.connect()
    con.execute("install spatial; load spatial;")
    con.execute("set enable_progress_bar = false")
    rows = con.execute(f"""
        select a.name, count(*) as n
        from read_parquet('{points}') p
        join st_read('{polygons}') a
          on st_intersects(st_point(p.lon, p.lat), a.geom)
        group by 1 order by n desc
    """).fetchall()
    return {"groups": len(rows), "total": sum(r[1] for r in rows),
            "top": [list(r) for r in rows[:3]]}


def run_geopandas(points, polygons):
    import duckdb
    import geopandas as gpd
    frame = duckdb.connect().execute(
        f"select lon, lat from read_parquet('{points}')").df()
    areas = gpd.read_file(polygons)
    pts = gpd.GeoDataFrame(frame,
                           geometry=gpd.points_from_xy(frame.lon, frame.lat),
                           crs=areas.crs)
    joined = gpd.sjoin(pts, areas, predicate="intersects")
    counts = joined.groupby("name", dropna=False).size().sort_values(ascending=False)   # keep unnamed areas, as SQL does
    return {"groups": len(counts), "total": int(counts.sum()),
            "top": [[k, int(v)] for k, v in list(counts.items())[:3]]}


if __name__ == "__main__":
    engine, points, polygons = sys.argv[1], sys.argv[2], sys.argv[3]
    started = time.perf_counter()
    result = run_duckdb(points, polygons) if engine == "duck" \
        else run_geopandas(points, polygons)
    elapsed = time.perf_counter() - started
    peak_mb = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
    print(json.dumps({"engine": engine, "secs": round(elapsed, 2),
                      "peak_mb": round(peak_mb, 1), "result": result}))

dropna=False is there because the first run of this benchmark failed its own assertion: DuckDB counted 7,324 town–province pairs and GeoPandas 7,321. Seven Natural Earth provinces have no name, and pandas' groupby silently drops rows whose key is missing, while SQL's group by keeps them as a NULL group. The join was identical; the aggregation was not.

Note that the GeoPandas path uses DuckDB to read the Parquet file. That is deliberate: it isolates the join difference from the Parquet-reading difference, and it is what a real GeoPandas user would do with a large Parquet input.

Example 2 — the runner, with repeats and assertions

import json
import statistics
import subprocess
import sys


def bench(engine, script, points, polygons, repeats=3):
    runs = []
    for _ in range(repeats):
        out = subprocess.run([sys.executable, script, engine, points, polygons],
                             capture_output=True, text=True, check=True)
        runs.append(json.loads(out.stdout.strip().splitlines()[-1]))
    return {
        "engine": engine,
        "best_secs": min(r["secs"] for r in runs),
        "median_secs": statistics.median(r["secs"] for r in runs),
        "peak_mb": max(r["peak_mb"] for r in runs),
        "result": runs[0]["result"],
        "runs": repeats,
    }


def compare(script, points, polygons, repeats=3):
    duck = bench("duck", script, points, polygons, repeats)
    gpd_ = bench("gpd", script, points, polygons, repeats)

    assert duck["result"]["total"] == gpd_["result"]["total"], (
        f"engines disagree: {duck['result']['total']:,} vs "
        f"{gpd_['result']['total']:,} matched pairs. Check the CRS on both sides.")

    print(f"{'engine':10} {'best':>8} {'median':>8} {'peak MB':>9}  ({repeats} runs)")
    for row in (duck, gpd_):
        print(f"{row['engine']:10} {row['best_secs']:8.2f} "
              f"{row['median_secs']:8.2f} {row['peak_mb']:9,.0f}")
    print(f"\nspeed  {gpd_['best_secs'] / duck['best_secs']:.1f}× "
          f"memory {gpd_['peak_mb'] / duck['peak_mb']:.1f}×  in DuckDB's favour")
    print(f"matched pairs {duck['result']['total']:,} in {duck['result']['groups']:,} groups; "
          f"same top three: {duck['result']['top'] == gpd_['result']['top']}\n")


for label, points, polygons, repeats in (("7,342 towns × 4,596 provinces", "towns.parquet", "provinces.shp", 3),
                                         ("US points × 51 states", "us_points.parquet", "us_states.gpkg", 3),
                                         ("all points × 4,596 provinces", "places.parquet", "provinces.shp", 1)):
    print(label)
    compare("bench.py", points, polygons, repeats)

7,342 towns × 4,596 provinces
engine         best   median   peak MB  (3 runs)
duck           0.26     0.26       215
gpd            0.44     0.44       318

speed  1.7× memory 1.5×  in DuckDB's favour
matched pairs 7,324 in 2,709 groups; same top three: True

US points × 51 states
engine         best   median   peak MB  (3 runs)
duck           3.66     3.67       215
gpd           13.83    13.97       932

speed  3.8× memory 4.3×  in DuckDB's favour
matched pairs 2,198,712 in 51 groups; same top three: True

all points × 4,596 provinces
engine         best   median   peak MB  (1 runs)
duck          50.85    50.85       291
gpd          121.80   121.80    23,813

speed  2.4× memory 81.7×  in DuckDB's favour
matched pairs 12,950,244 in 4,419 groups; same top three: True

Example 3 — a size sweep, because the answer changes with scale

def scale_sweep(script, points_parquet, polygons, sizes=(10_000, 100_000,
                                                         1_000_000, 10_000_000)):
    """The crossover is the interesting result, not any single number."""
    import duckdb
    import tempfile
    import os

    con = duckdb.connect()
    print(f"{'rows':>12} {'duckdb':>9} {'geopandas':>10} {'winner':>10}")
    for n in sizes:
        with tempfile.NamedTemporaryFile(suffix=".parquet", delete=False) as handle:
            subset = handle.name
        con.execute(f"""copy (select * from read_parquet('{points_parquet}')
                              limit {n}) to '{subset}' (format parquet)""")
        duck = bench("duck", script, subset, polygons, repeats=3)
        gpd_ = bench("gpd", script, subset, polygons, repeats=3)
        winner = "duckdb" if duck["best_secs"] < gpd_["best_secs"] else "geopandas"
        print(f"{n:12,} {duck['best_secs']:9.2f} {gpd_['best_secs']:10.2f} {winner:>10}")
        os.unlink(subset)


scale_sweep("bench.py", "places.parquet", "provinces.shp")

        rows    duckdb  geopandas     winner
      10,000      0.20       0.41     duckdb
     100,000      0.62       0.98     duckdb
   1,000,000      2.64       9.49     duckdb
  10,000,000     48.15      90.96     duckdb

On this machine there was no crossover: DuckDB was faster at every size, from 0.20 s against 0.41 s at 10,000 rows to 48 s against 91 s at ten million. Below a few hundred thousand rows the difference is a fraction of a second, which is convenience territory; the size at which the memory difference starts to decide whether the job runs is the useful output.

Measured in Example 3. The ratio moves with size in both directions.

Reads inside the timer, results that must agree

Why the file read belongs inside the timer

DuckDB's whole design point is that it reads lazily inside the query. Excluding the read from its timer measures nothing, and excluding GeoPandas' explicit read measures a different task.

Timing the complete task — files in, answer out — is the only comparison that corresponds to what a user experiences. It also means the benchmark automatically credits DuckDB for column pruning, which is a real advantage that a join-only timer hides.

Why memory is the more decisive axis

The measured global join was 2.4× faster in DuckDB and used an eighty-second of the memory. On a machine with 16 GB of RAM the GeoPandas version simply does not run, which is a qualitative difference that no speed ratio expresses.

Memory also determines what else can run at the same time. A pipeline stage that peaks at 23.8 GB cannot run at all on a 16 GB machine; one that peaks at 291 MB can run dozens of times over.

Why the equality assertion matters more than any timing

A DuckDB spatial join across mismatched coordinate systems returns zero rows and no error. Zero rows is extremely fast.

Every benchmark should therefore compare outputs before comparing times. In the measured runs both engines produced 12,950,244 matched pairs across 4,419 groups with identical leaders — and that agreement is what licenses the timing comparison. Without it, the fast number could be measuring nothing at all.

Why the scale, not one ratio, is the useful output

"DuckDB is 2× faster" is not actionable, because the ratio moves with scale. The measured speed advantage was 1.7× at 7,342 points, 3.8× at 2.2 million and 2.4× at 13.5 million; the memory advantage went 1.5×, 4.3×, 82×.

The number worth extracting from a benchmark is where your workload sits on that curve. While both finish in seconds within the memory you have, use whichever library is more convenient; beyond that, the choice starts to determine whether the job completes.

Time differed by 1.7–3.8×; memory by up to 82×.

Machines, repeats and units of `ru_maxrss`

Report the machine. Core count, RAM and storage type change every number here.
Take the minimum of several runs, not the mean; it is least polluted by other load.
ru_maxrss is kilobytes on Linux and bytes on macOS. Normalise it.
Disable DuckDB's external file cache when measuring remote reads.
A subprocess per engine is the only reliable way to attribute peak memory.
Warm the cache or state that you did not. Both are valid; silence is not.
Check thread counts. DuckDB uses all cores by default; GeoPandas is largely single-threaded.
Compare the outputs first. A fast wrong answer is the most common benchmark result.

Internal links

DuckDB spatial explained: a spatial database that is just a file — where the headline numbers come from
How to run a spatial join in DuckDB — the operation benchmarked here
DuckDB's memory model: streaming, spilling and limits — why the memory differs so much
CRS in DuckDB: why ST_Transform moves your data to the wrong place — the silent zero-row failure
How to profile a slow Python GIS script — profiling rather than benchmarking
Why GeoPandas is slow, explained — the other side of the comparison
How to move data between DuckDB and GeoPandas — using both
DuckDB or PostGIS: which one a workflow needs — the third engine

FAQ

How should I time a GIS engine comparison?

Time the whole task — files in, answer out — in a separate process per engine, take the minimum of several runs, and record the machine.

Should I measure memory too?

Yes, and it is often the more important number. The measured global join was 2.4× faster in DuckDB and used 82× less memory: 291 MB against 23,813 MB.

Why do I need a separate process per engine?

Peak resident memory is process-wide, so running both in one process attributes the maximum to whichever ran second.

What is the most common benchmark mistake?

Not comparing the results. A DuckDB join across mismatched CRSs returns zero rows with no error, and zero rows is very fast.

Is DuckDB always faster than GeoPandas?

Not by much at small sizes. On 7,342 points against 4,596 polygons DuckDB took 0.26 s and GeoPandas 0.44 s — a difference nobody waits for. The gap that matters opens with size, and in memory first.

Should the cache be warm or cold?

Either, as long as you say which. Warm reflects repeated interactive use; cold reflects a first run. Silence about it makes the numbers uninterpretable.