🧭 Concept DuckDB & analytic SQL Performance & scale DuckDB for Spatial Data · step 7 Notebook By · Updated

How DuckDB Uses (and Does Not Use) a Spatial Index

On this page

Problem statement

The spatial join is slow, so you do what you would do in PostGIS: create a spatial index.

create index areas_rtree on areas using rtree (geom);

The index builds in two hundredths of a second. The join runs at exactly the same speed. Checking the plan explains why:

┌───────────────────────────┐
│        SPATIAL_JOIN       │
│      Join Type: INNER     │
│ ST_Intersects(geom, geom) │
└───────────────────────────┘

There is no index scan anywhere in it. DuckDB's R-tree index is real, useful and does something quite specific — and joining two tables is not that thing.

Quick answer

The R-tree accelerates a filter against a fixed geometry. It does not participate in a table-to-table join, which uses a dedicated operator instead.

-- the index helps here: one constant geometry, many rows scanned
select count(*) from areas
where st_intersects(geom, st_makeenvelope(-74.05, 40.65, -73.85, 40.9));

-- the index does not participate here: two tables, a spatial join operator
select count(*) from points p join areas a on st_intersects(p.geom, a.geom);
query = "select count(*) from points p join areas a on st_intersects(p.geom, a.geom)"
plan = con.execute("explain " + query).fetchall()[0][1]   # look for RTREE_INDEX_SCAN vs SPATIAL_JOIN
print("SPATIAL_JOIN" in plan, "RTREE_INDEX_SCAN" in plan)
True False

Measured on 4,596 province polygons: the bounding-box query dropped from 9.2 ms to 0.33 ms with the index; a join of 2.2 million points against the same polygons took 4.2 s either way (Examples 1 and 4).

The query shapes a DuckDB R-tree index accelerates, and the one it does not. one-to-many yes, many-to-many nothe index is used✓filter against a fixed geometry✓a map viewport or study area✓repeated point-in-polygon lookups✓measured: 9.2 ms → 0.33 msthe index is not used✕table-to-table spatial joins✕ST_DWithin between two tables✕anything with no constant geometry✕the plan shows SPATIAL_JOINCreating an index to speed up a join is a wasted afternoon — read the plan first spatialworkflow.io
The R-tree answers "what is in this box?", not "join these two tables".

What an R-tree does, and what DuckDB does instead

1. Know what an R-tree does

An R-tree stores each geometry's bounding box in a balanced tree. A query with a fixed search box descends the tree, discards whole subtrees whose boxes cannot intersect, and returns a small candidate set for exact testing.

That is a one-to-many operation: one query geometry, many stored ones. It is the shape of "which features are in this view?" and "which polygon contains this point?", and it is what PostGIS uses for almost everything.

2. Know what DuckDB does for a join instead

A join is many-to-many. Probing an R-tree once per row of the other table works, and it is not the fastest strategy when both sides are large: each probe pays tree traversal, and there is no opportunity to process values in vectors.

DuckDB's SPATIAL_JOIN operator builds a spatial structure over one side and streams the other through it in batches, which suits its vectorised execution model. Measured: 13,472,195 points against 4,596 polygons in 50.9 s using 291 MB, without any index existing (benchmark).

3. Create the index when the query has a fixed geometry

create index areas_rtree on areas using rtree (geom);

Building it on 4,596 polygons took 0.024 s (Example 4). The queries it helps:

  • bounding-box filters — a map viewport, a study area
  • ST_Intersects(geom, <constant>), ST_Within(geom, <constant>)
  • repeated point-in-polygon lookups against the same table

The queries it does not help:

  • table-to-table joins
  • ST_DWithin between two tables
  • anything without a fixed geometry in the predicate

4. Read the plan rather than guessing

sql = "select count(*) from areas where st_intersects(geom, st_makeenvelope(-74.05, 40.65, -73.85, 40.9))"
plan = con.execute("explain " + sql).fetchall()[0][1]
print(plan)
print("uses the index:", "RTREE_INDEX_SCAN" in plan.upper())
print("spatial join:  ", "SPATIAL_JOIN" in plan.upper())
┌───────────────────────────┐
│    UNGROUPED_AGGREGATE    │
│    ────────────────────   │
│        Aggregates:        │
│        count_star()       │
└─────────────┬─────────────┘
┌─────────────┴─────────────┐
│          SEQ_SCAN         │
│    ────────────────────   │
│           Table:          │
│   index_demo.main.areas   │
│                           │
│   Type: Sequential Scan   │
│                           │
│          Filters:         │
│    ST_Intersects(geom,    │
│  POLYGON ((-74.05 40.65,  │
│ -74.05 40.9, -73.85 40.9, │
│ -73.85 40.65, -74.05 40.65│
│            )))            │
│                           │
│         ~919 rows         │
└───────────────────────────┘

uses the index: False
spatial join:   False

This takes seconds and settles the question for the specific query in front of you, which is more reliable than any general rule — including this article's.

5. Understand what an index costs

An R-tree is persisted with the table in a DuckDB database file, and it must be maintained on insert. For a static analytical table that is irrelevant; for a table you rewrite each run it is pure overhead.

It also does nothing at all when the table is a file being read in place. read_parquet('x.parquet') has no index and cannot have one — which is fine, because Parquet has its own pruning mechanism.

6. Use row-group pruning instead when reading files

For file-backed data, the equivalent of an index is sorting. Parquet stores per-row-group minimum and maximum values, so a filter on a sorted column skips whole groups.

Measured over HTTP with every byte counted, filtering the same 13.5-million-row dataset by one value (measured here):

                sorted file    shuffled file
bytes read         0.26 MB        10.26 MB
requests                 5             221
time                 0.01 s          1.50 s

Same rows, same query, 39× the bytes and 150× the time. For files, sorting is the index.

What DuckDB does instead of index probes for a table-to-table spatial join. the SPATIAL_JOIN operatorbuild sidea structure over onetableprobe side inbatchesvectorisedcandidate pairsbounding boxes firstexact predicateonly on survivorsPer-row index probes would defeat the vectorised execution the engine is built on spatialworkflow.io
13.5 million points against 4,596 polygons in 50.9 s, with no index existing.

Measuring when the index is used

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 load both into tables in a DuckDB database: the provinces as areas, and the 2.24 million GeoNames points in the United States as points, each with a real geometry column — an index needs a table, not a file read in place.

con = duckdb.connect("index_demo.duckdb")
con.execute("install spatial; load spatial; set enable_progress_bar = false")
con.execute("create or replace table areas as select name, geom from st_read('provinces.shp')")
con.execute("""create or replace table points as
               select name, st_point(lon, lat) as geom from read_parquet('places.parquet') where country = 'US'""")
print(con.execute("select (select count(*) from areas), (select count(*) from points)").fetchone())
(4596, 2241404)

Example 1 — measuring whether the index is used

import duckdb
import time


def index_benefit(con, table, geom_column="geom", box=None, repeats=5):
    """Time a bounding-box query with and without an index, and read the plan."""
    box = box or "st_makeenvelope(-74.05, 40.65, -73.85, 40.9)"
    query = f"select count(*) from {table} where st_intersects({geom_column}, {box})"

    def best_time():
        best = None
        for _ in range(repeats):
            start = time.perf_counter()
            con.execute(query).fetchone()
            elapsed = time.perf_counter() - start
            best = elapsed if best is None else min(best, elapsed)
        return best

    without = best_time()
    con.execute(f"create index if not exists {table}_rtree on {table} "
                f"using rtree ({geom_column})")
    with_index = best_time()

    plan = con.execute("explain " + query).fetchall()[0][1]
    print(f"without index {without * 1000:7.2f} ms")
    print(f"with index    {with_index * 1000:7.2f} ms   "
          f"({without / max(with_index, 1e-9):.1f}×)")
    print(f"plan uses RTREE_INDEX_SCAN: {'RTREE' in plan.upper()}")
    return without, with_index


con.execute("drop index if exists areas_rtree")
index_benefit(con, "areas")
without index    9.15 ms
with index       0.33 ms   (27.5×)
plan uses RTREE_INDEX_SCAN: True

Example 2 — the two query shapes, side by side

def compare_shapes(con, points_table, areas_table):
    """One query the index helps, one it cannot."""
    filter_sql = f"""
        select count(*) from {areas_table}
        where st_intersects(geom, st_makeenvelope(-74.05, 40.65, -73.85, 40.9))
    """
    join_sql = f"""
        select count(*) from {points_table} p
        join {areas_table} a on st_intersects(p.geom, a.geom)
    """
    for label, sql in (("fixed-geometry filter", filter_sql), ("table join", join_sql)):
        plan = con.execute("explain " + sql).fetchall()[0][1].upper()
        operator = ("RTREE_INDEX_SCAN" if "RTREE" in plan
                    else "SPATIAL_JOIN" if "SPATIAL_JOIN" in plan
                    else "SEQ_SCAN + FILTER")
        print(f"{label:24} → {operator}")


compare_shapes(con, "points", "areas")
fixed-geometry filter    → RTREE_INDEX_SCAN
table join               → SPATIAL_JOIN

Example 3 — sorting a Parquet file as the file-backed equivalent

def write_clustered(con, source, target, cluster_on, compression="zstd"):
    """Row-group pruning is the index you get for a file."""
    import os
    con.execute(f"""
        copy (select * from {source} order by {cluster_on})
        to '{target}' (format parquet, compression {compression})
    """)

    import pyarrow.parquet as pq
    meta = pq.read_metadata(target)
    print(f"{meta.num_rows:,} rows in {meta.num_row_groups} row groups, "
          f"{os.path.getsize(target) / 1e6:,.1f} MB")
    print(f"filters on {cluster_on} can now skip groups whose min/max exclude the value")


write_clustered(con, "read_parquet('places.parquet')", "places_by_country.parquet", "country")
13,472,195 rows in 110 row groups, 187.7 MB
filters on country can now skip groups whose min/max exclude the value

Sorting also shrinks the file, because clustered values compress better: measured, the same data was 446 MB sorted and 741 MB shuffled — 66% larger for the same rows.

Example 4 — the join, with and without the index

JOIN = "select count(*) from points p join areas a on st_intersects(p.geom, a.geom)"
timings = {}
for state in ("without", "with"):
    con.execute("drop index if exists areas_rtree")
    if state == "with":
        start = time.perf_counter()
        con.execute("create index areas_rtree on areas using rtree (geom)")
        print(f"building the R-tree on {con.execute('select count(*) from areas').fetchone()[0]:,} polygons: "
              f"{time.perf_counter() - start:.3f} s")
    start = time.perf_counter()
    n = con.execute(JOIN).fetchone()[0]
    timings[state] = time.perf_counter() - start
    print(f"join {state:7} index: {n:,} pairs in {timings[state]:.2f} s")
plan = con.execute("explain " + JOIN).fetchall()[0][1]
print("operator in the plan:", "SPATIAL_JOIN" if "SPATIAL_JOIN" in plan else "other",
      "| index scan in the plan:", "RTREE" in plan.upper())
join without index: 2,200,310 pairs in 4.30 s
building the R-tree on 4,596 polygons: 0.024 s
join with    index: 2,200,310 pairs in 4.30 s
operator in the plan: SPATIAL_JOIN | index scan in the plan: False

Why a join probes no index

Why a join does not use an index

Consider probing an R-tree once per row of a 13-million-row table. Each probe descends a tree, chases pointers and returns a candidate list — a pointer-heavy, branch-heavy operation repeated 13 million times, with no opportunity to work on values in bulk.

DuckDB's whole design is vectorised: it processes columns in batches of a couple of thousand values, keeping the CPU's pipelines full. A dedicated spatial join operator that builds a structure once and streams batches through it fits that model; per-row index probes do not.

PostGIS makes the opposite choice because it is a row-oriented system where per-row index lookups are the natural motion, and because its planner has statistics and a cost model built around indexes.

Why the index still matters

The one-to-many shape is extremely common in an application: a map server asking "what is in this tile?", a lookup asking "which region contains this point?", a repeated query against a fixed study area.

For those, the R-tree turns a full scan into a tree descent. On the small measured table the absolute numbers are tiny — 9 ms to under 1 ms on 4,596 polygons — but the ratio grows with the table, and a web service doing this per request cares a great deal.

Why sorting is the index for file-backed data

There is no index on read_parquet('x.parquet'), and there does not need to be one. Parquet's row-group statistics provide coarse pruning for free, and their effectiveness depends entirely on whether the values are clustered.

The measured contrast is the whole argument: the same filter read 0.26 MB in 5 requests from a sorted file, and 10.26 MB in 221 requests from a shuffled one. Sorting on write is the analytical equivalent of choosing an index — it costs time once and pays back on every query.

Why reading the plan beats reasoning about it

Query planners change between versions. DuckDB's spatial support is developing quickly, and an operator that ignores the index today may use it in a later release.

EXPLAIN takes a second and tells you what this version does with this query on this data. It is the only claim in this article that will still be true in two years.

For file-backed data, sorting does what an index does for a table. the same filter, two file layoutssorted file · MB read0.26 MBshuffled file · MB read10.26 MBsorted · seconds0.01 sshuffled · seconds1.50 s011There is no index on read_parquet, and there does not need to be one spatialworkflow.io
Row-group statistics prune only when the values are clustered.

Persistent files, insert costs and the plan

  • create index ... using rtree (geom) requires the spatial extension loaded and a persistent database file.
  • Indexes need a table, not a file read in place. Create a table first if you want one.
  • Index maintenance costs on insert. For a rebuilt-every-run table, skip it.
  • EXPLAIN ANALYZE gives timings per operator — worth it when a query is unexpectedly slow.
  • Bounding-box filters benefit most; exact predicates still need a refinement step after the index.
  • Sorting on write helps every file-backed query, and shrinks the file.
  • Row-group size is a tuning knob for pruning granularity.
  • Do not create an index to fix a slow join. Look at the plan first.

FAQ

Does DuckDB have a spatial index?

Yes — an R-tree, created with create index ... using rtree (geom). It accelerates filters against a fixed geometry and is not used by table-to-table spatial joins.

Why did creating an index not speed up my join?

Because the join uses a dedicated SPATIAL_JOIN operator rather than index probes. EXPLAIN on the query shows the operator and no index scan.

When is an R-tree index worth creating?

When many queries filter the same table by a fixed geometry — a viewport, a study area, a point-in-polygon lookup. Measured, it took a bounding-box query from 9 ms to under 1 ms.

How do I speed up a slow spatial join then?

Filter first, simplify complex polygons, make sure both sides are in the same CRS, and set a memory limit so the engine spills rather than thrashes.

What is the equivalent of an index for a Parquet file?

Sorting. Row-group statistics let the reader skip groups: the same filter read 0.26 MB from a sorted file and 10.25 MB from a shuffled one.

How do I check whether my query uses the index?

Run EXPLAIN and look for RTREE_INDEX_SCAN in the plan. It takes a second and is more reliable than any rule of thumb.