PostGIS Spatial Indexes Explained: GiST and the Two-Phase Filter

Problem statement

The same query, twice, on the same data — 32,683 HM Land Registry parcels along the River Aire in Leeds and the Environment Agency's 440 flood-zone polygons over them:

import time

def timed(sql):
    with engine.begin() as con:
        con.execute(text("SET LOCAL max_parallel_workers_per_gather = 0"))
        start = time.perf_counter()
        value = con.execute(text(sql)).scalar()
    return value, time.perf_counter() - start

JOIN = "SELECT COUNT(*) FROM parcels p JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)"
n, before = timed(JOIN)
print(f"no index:   {n:,} rows in {before:6.2f} s")

with engine.begin() as con:
    con.execute(text("CREATE INDEX parcels_geom_idx ON parcels USING GIST (geom)"))
    con.execute(text("CREATE INDEX flood_zones_geom_idx ON flood_zones USING GIST (geom)"))
    con.execute(text("ANALYZE parcels"))
    con.execute(text("ANALYZE flood_zones"))
n, after = timed(JOIN)
print(f"with index: {n:,} rows in {after:6.2f} s  ({before / after:.0f}× faster)")

no index:   3,679 rows in  33.86 s
with index: 3,679 rows in   0.23 s  (146× faster)

Nothing about the data or the query changed. Meanwhile this one gains nothing from the index at all:

STATION = "ST_SetSRID(ST_MakePoint(429900, 433200), 27700)"      # Leeds station
for label, predicate in [("ST_Distance < 500", f"ST_Distance(geom, {STATION}) < 500"),
                         ("ST_DWithin 500   ", f"ST_DWithin(geom, {STATION}, 500)")]:
    n, seconds = timed(f"SELECT COUNT(*) FROM parcels WHERE {predicate}")
    print(f"{label}  {n:,} parcels in {seconds * 1000:7.1f} ms")

ST_Distance < 500  610 parcels in    16.6 ms
ST_DWithin 500     610 parcels in     1.0 ms

Both behaviours come from the same mechanism. Understanding it tells you which queries an index will rescue, which it will not, and why.

Quick answer

Build a GiST index on every geometry column, and write predicates the planner can use:

CREATE INDEX parcels_geom_idx ON parcels USING GIST (geom);
ANALYZE parcels;                     -- the planner needs statistics too

The single most valuable substitution: replace ST_Distance(a, b) < d with ST_DWithin(a, b, d). They return the same rows; only one of them can use an index.

What the index stores, two phases, queries that use it

1. Understand what the index actually stores

A GiST index on a geometry column does not store the geometry. It stores each geometry's bounding box arranged in a tree, so that a box query can skip most of the table.

That substitution is what makes the index fast and what makes it approximate. A bounding box can overlap when the geometries do not — an L-shaped parcel and a river have overlapping boxes and no shared point. So a box test can produce false positives, but never false negatives: if the boxes do not overlap, the geometries certainly do not intersect.

That asymmetry is the whole design. A cheap test that can only over-report is exactly what you want as a first pass.

2. Follow the two phases

with engine.begin() as con:
    plan = con.execute(text(f"EXPLAIN (ANALYZE, COSTS OFF, TIMING OFF) {JOIN}")).scalars().all()

print("\n".join(line for line in plan if "Planning" not in line and "Execution" not in line))
print(q("""SELECT (SELECT COUNT(*) FROM parcels) * (SELECT COUNT(*) FROM flood_zones)          AS pairs,
                  (SELECT COUNT(*) FROM parcels p JOIN flood_zones f ON p.geom && f.geom)   AS boxes_overlap,
                  (SELECT COUNT(*) FROM parcels p JOIN flood_zones f
                     ON ST_Intersects(p.geom, f.geom))                                     AS intersect,
                  (SELECT MAX(ST_NPoints(geom)) FROM flood_zones)                          AS most_vertices
""").to_string(index=False))

Aggregate (actual rows=1 loops=1)
  ->  Nested Loop (actual rows=3679 loops=1)
        ->  Seq Scan on flood_zones f (actual rows=440 loops=1)
        ->  Index Scan using parcels_geom_idx on parcels p (actual rows=8 loops=440)
              Index Cond: (geom && f.geom)
              Filter: st_intersects(geom, f.geom)
              Rows Removed by Filter: 60
   pairs  boxes_overlap  intersect  most_vertices
14380080          30270       3679           6715

Three lines carry the whole story:

• Index Cond: (geom && ...) — the index phase. && means "bounding boxes overlap". This is what the tree can answer.
• Filter: st_intersects(...) — the exact phase. Real geometry comparison, run only on the survivors.
• Rows Removed by Filter — the false positives: candidates whose boxes overlapped without the geometries touching, shown per loop.

PostGIS adds the && condition automatically — you write ST_Intersects and get the two-phase plan for free. Which table is scanned in order and which through its index is the planner's choice; here it walks the smaller one.

3. Recognise the queries an index cannot help

The index answers exactly one kind of question: which rows have a bounding box overlapping this box? Anything that cannot be reduced to that form scans the table:

def uses_index(where):
    with engine.begin() as con:
        plan = "\n".join(con.execute(text(f"EXPLAIN SELECT COUNT(*) FROM parcels WHERE {where}")).scalars())
    return "Index" in plan

BOX_4326 = "ST_MakeEnvelope(-1.56, 53.79, -1.54, 53.80, 4326)"
for where in [f"ST_Distance(geom, {STATION}) < 500",
              f"ST_DWithin(geom, {STATION}, 500)",
              "ST_Area(geom) > 10000",
              f"ST_Transform(geom, 4326) && {BOX_4326}",
              f"geom && ST_Transform({BOX_4326}, 27700)"]:
    print(f"{'index' if uses_index(where) else 'scan ':5}  {where}")

scan   ST_Distance(geom, ST_SetSRID(ST_MakePoint(429900, 433200), 27700)) < 500
index  ST_DWithin(geom, ST_SetSRID(ST_MakePoint(429900, 433200), 27700), 500)
scan   ST_Area(geom) > 10000
scan   ST_Transform(geom, 4326) && ST_MakeEnvelope(-1.56, 53.79, -1.54, 53.80, 4326)
index  geom && ST_Transform(ST_MakeEnvelope(-1.56, 53.79, -1.54, 53.80, 4326), 27700)

ST_DWithin expands the reference geometry's box by 500 units, asks the index for candidates, then measures exactly. ST_Distance computes a distance for every row before comparing — the column is inside a function call, so no index applies. The transform pair is the same rule: transforming the column hides it from the index, transforming the constant leaves the column bare.

Keep the indexed column bare on one side of the operator. That single rule covers almost every case.

4. Build, and then tell the planner

CREATE INDEX parcels_geom_idx ON parcels USING GIST (geom);
ANALYZE parcels;

ANALYZE is not optional. The planner chooses between an index scan and a sequential scan by estimating how many rows will match, and without statistics its estimate for a fresh index is a guess. A common report of "I built the index and nothing changed" is a missing ANALYZE — especially after a bulk load, which is exactly when the table has just changed the most.

On a live table, build without locking writes. CONCURRENTLY cannot run inside a transaction block, so it needs an autocommit connection:

with engine.connect().execution_options(isolation_level="AUTOCOMMIT") as con:
    for concurrently in ("", "CONCURRENTLY"):
        con.execute(text("DROP INDEX IF EXISTS leeds_parcels_geom_idx"))
        start = time.perf_counter()
        con.execute(text(f"CREATE INDEX {concurrently} leeds_parcels_geom_idx ON leeds_parcels USING GIST (geom)"))
        print(f"CREATE INDEX {concurrently or 'plain':13} on {LEEDS_ROWS:,} rows: {time.perf_counter() - start:5.2f} s")

CREATE INDEX plain         on 293,729 rows:  0.30 s
CREATE INDEX CONCURRENTLY  on 293,729 rows:  0.35 s

CONCURRENTLY scans the table twice and waits for running transactions, so on a busy table it is slower — on this idle one the difference is small — but it does not block writes. PostGIS 3 builds a GiST index from sorted input, which is why a third of a million parcels index in well under a second. On a table being loaded, do the opposite — load first, index afterwards. Maintaining an index during a large insert costs more than building it once at the end.

5. Check that it is actually being used

print(q("""
    SELECT relname, indexrelname, idx_scan, idx_tup_read
    FROM pg_stat_user_indexes
    WHERE relname IN ('parcels', 'flood_zones', 'leeds_parcels')
    ORDER BY idx_scan
""").to_string(index=False))

      relname           indexrelname  idx_scan  idx_tup_read
  flood_zones   flood_zones_geom_idx         0             0
leeds_parcels leeds_parcels_geom_idx         0             0
      parcels       parcels_geom_idx      1761        121876

idx_scan = 0 on an index that has existed for weeks means nothing is using it. Either no query needs it, or every query that should is written in a form the planner cannot exploit. Both are worth knowing; an unused index still costs write time and disk.

Measuring the index and reading plans

The data used on this page

HM Land Registry's INSPIRE Index Polygons for Leeds and the Environment Agency's Flood Map for Planning (flood zones 2 and 3, rivers and sea) from its Web Feature Service, loaded into a local PostGIS 3.4 (the postgis/postgis:16-3.4 Docker image; connection details from PGPORT and PGPASSWORD). parcels and flood_zones are the 6 × 4 km stretch of the Aire valley from Kirkstall to the city centre; leeds_parcels is all of Leeds. Every table starts with no index — GeoPandas's to_postgis quietly builds a GiST index of its own, named idx_<table>_geom, so the setup drops it.

import io
import os
import time

import geopandas as gpd
import pandas as pd
import requests
from sqlalchemy import create_engine, text

# the server: PostGIS 3.4 in Docker on localhost:55432
if os.system("docker start guide-postgis-55432 >/dev/null 2>&1") != 0:
    os.system("docker run -d --name guide-postgis-55432 -e POSTGRES_USER=user -e POSTGRES_PASSWORD=spatialworkflow "
              "-e POSTGRES_DB=gis -p 127.0.0.1:55432:5432 postgis/postgis:16-3.4 >/dev/null")
for _ in range(60):                                  # the image restarts once while it installs PostGIS
    if os.system("docker exec guide-postgis-55432 pg_isready -h 127.0.0.1 -U user -d gis >/dev/null 2>&1") == 0:
        break
    time.sleep(2)
os.environ.setdefault("PGPORT", "55432")
os.environ.setdefault("PGPASSWORD", "spatialworkflow")
engine = create_engine("postgresql+psycopg://user@localhost/gis")


def q(sql, **params):
    return pd.read_sql(text(sql), engine, params=params)


reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Leeds.zip",
                               timeout=600)
reply.raise_for_status()
with open("leeds_inspire.zip", "wb") as f:
    f.write(reply.content)
leeds = gpd.read_file("zip://leeds_inspire.zip!Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"])
leeds = leeds.rename(columns={"INSPIREID": "id"}).rename_geometry("geom").set_crs(27700, allow_override=True)

zones = requests.get("https://environment.data.gov.uk/spatialdata/flood-map-for-planning-flood-zones/wfs",
                     timeout=300, params={
                         "service": "WFS", "version": "2.0.0", "request": "GetFeature",
                         "typeNames": "dataset-04532375-a198-476e-985e-0579a0a11b47:Flood_Zones_2_3_Rivers_and_Sea",
                         "bbox": "419000,425000,446000,450000,urn:ogc:def:crs:EPSG::27700",
                         "outputFormat": "application/json", "count": 100000})
zones.raise_for_status()
flood = gpd.read_file(io.BytesIO(zones.content)).set_crs(27700, allow_override=True).rename_geometry("geom")

AIRE = (425000, 432000, 431000, 436000)
with engine.begin() as con:
    con.execute(text("DROP TABLE IF EXISTS parcels, flood_zones, leeds_parcels CASCADE"))
leeds.to_postgis("leeds_parcels", engine, index=False, chunksize=50_000)
leeds.cx[AIRE[0]:AIRE[2], AIRE[1]:AIRE[3]].to_postgis("parcels", engine, index=False)
flood.cx[AIRE[0]:AIRE[2], AIRE[1]:AIRE[3]].to_postgis("flood_zones", engine, index=False)
with engine.begin() as con:
    for table in ("parcels", "flood_zones", "leeds_parcels"):        # to_postgis builds a GiST index
        con.execute(text(f"DROP INDEX IF EXISTS idx_{table}_geom"))
    con.execute(text("SELECT pg_stat_reset()"))
LEEDS_ROWS = len(leeds)
print(q("SELECT (SELECT COUNT(*) FROM parcels) AS parcels, (SELECT COUNT(*) FROM flood_zones) AS flood_zones, "
        "(SELECT COUNT(*) FROM leeds_parcels) AS leeds_parcels").to_string(index=False))

 parcels  flood_zones  leeds_parcels
   32682          440         293729

Example 1: measuring the index from Python

def time_query(sql, label, use_index=True):
    with engine.begin() as con:
        con.execute(text("SET LOCAL max_parallel_workers_per_gather = 0"))
        con.execute(text(f"SET LOCAL enable_indexscan = {use_index}"))
        con.execute(text(f"SET LOCAL enable_bitmapscan = {use_index}"))
        t0 = time.perf_counter()
        n = con.execute(text(sql)).scalar()
        dt = time.perf_counter() - t0
    print(f"{label:<20} {n:>10,} rows  {dt:>8.2f} s")
    return dt


slow = time_query(JOIN, "without index", use_index=False)
fast = time_query(JOIN, "with index",    use_index=True)
print(f"speedup: {slow / fast:.0f}×")

without index             3,679 rows     33.47 s
with index                3,679 rows      0.23 s
speedup: 145×

SET LOCAL enable_indexscan = off is a transaction-local planner hint, and it is the honest way to measure the difference — dropping and rebuilding the index changes the table's cached state as well.

The row count is identical in both runs. That is the point worth checking: the index changes how the answer is found, never what the answer is. If the counts differ, something else is wrong.

Example 2: reading the plan from Python and acting on it

import json


def explain(sql, params=None):
    with engine.begin() as con:
        raw = con.execute(
            text(f"EXPLAIN (ANALYZE, BUFFERS, FORMAT JSON) {sql}"), params or {}
        ).scalar()
    plan = raw[0]["Plan"] if isinstance(raw, list) else json.loads(raw)[0]["Plan"]

    seq_scans, index_scans, removed = [], [], 0

    def walk(node):
        nonlocal removed
        kind = node["Node Type"]
        loops = node.get("Actual Loops", 1)
        if kind == "Seq Scan":
            seq_scans.append((node["Relation Name"], node.get("Actual Rows", 0) * loops))
        elif "Index" in kind:
            index_scans.append((node.get("Index Name"), node.get("Actual Rows", 0) * loops))
        removed += node.get("Rows Removed by Filter", 0) * loops
        for child in node.get("Plans", []):
            walk(child)

    walk(plan)
    print(f"total time     {plan['Actual Total Time']:.1f} ms   (rows per loop are averaged, so totals are approximate)")
    for name, rows in index_scans:
        print(f"  index scan   {name} → {rows:,} rows")
    for name, rows in seq_scans:
        flag = "  ← candidate for an index" if rows > 10_000 else ""
        print(f"  seq scan     {name} → {rows:,} rows{flag}")
    print(f"  false positives removed by exact test: {removed:,}")
    return plan


plan = explain(JOIN)

total time     226.0 ms   (rows per loop are averaged, so totals are approximate)
  index scan   parcels_geom_idx → 3,520 rows
  seq scan     flood_zones → 440 rows
  false positives removed by exact test: 26,400

Two readings matter. A sequential scan over a few hundred flood zones is fine — with a table that small the planner is right to walk it in order and use the other table's index. A sequential scan over tens of thousands of rows in the same position is the finding.

And the false-positive count against the real hits says how well bounding boxes fit the geometry: here roughly seven candidates rejected for every parcel that really touches a flood zone. Flood zones are long, winding river corridors, so their boxes are poor approximations and the exact phase does much of the work. A much worse ratio would argue for subdividing the geometries, not another index.

Example 3: indexing a table you are loading from Python

def load_with_index(gdf, table, engine, srid=27700, chunksize=50_000):
    """Load, then index, then analyse — in that order."""
    gdf = gdf.to_crs(srid)
    gdf.to_postgis(table, engine, if_exists="replace",
                   index=False, chunksize=chunksize)
    with engine.begin() as con:
        con.execute(text(f'DROP INDEX IF EXISTS "idx_{table}_geom"'))       # GeoPandas made one; build ours
        con.execute(text(
            f'CREATE INDEX "{table}_geom_idx" ON "{table}" USING GIST (geom)'))
        con.execute(text(f'ANALYZE "{table}"'))
        size = con.execute(text(
            f"SELECT pg_size_pretty(pg_total_relation_size('{table}'))")).scalar()
        idx = con.execute(text(
            f"SELECT pg_size_pretty(pg_relation_size('{table}_geom_idx'))")).scalar()
    print(f"{table}: {len(gdf):,} rows, {size} total, {idx} index")


load_with_index(leeds, "leeds_parcels_loaded", engine)

leeds_parcels_loaded: 293,729 rows, 130 MB total, 12 MB index

The order is the point. to_postgis with if_exists="replace" drops the table and every index on it, so any index built beforehand is gone — and had it survived, maintaining it during the load would have cost more than building it afterwards. Note also that GeoPandas creates a GiST index of its own, idx_<table>_geom, after writing; the function drops it so the table ends with exactly one.

The index is 12 MB against 130 MB for the table and its index together, which is why "should I index this geometry column?" is almost always yes. Loading details are in how to write a GeoDataFrame to PostGIS.

Why spatial data needs its own index

A spatial index solves a problem that ordinary indexes cannot even express. A B-tree works because values have a total order: every number is before or after every other, so the tree can halve the search space at each step. Geometries have no such order. There is no meaningful sense in which one polygon is "less than" another, so there is nothing to sort by and nothing to bisect.

GiST — the Generalized Search Tree — sidesteps this by indexing something orderable that approximates the geometry: its bounding box. Boxes can be organised hierarchically, because a box can contain other boxes. The tree's internal nodes hold boxes that enclose all their children, so a query box that misses a node's box misses everything beneath it, and an entire subtree is skipped without examining a single geometry.

This is the same idea as the R-tree that GeoPandas builds in memory, and the two behave alike. The differences are operational rather than conceptual: the GiST index is stored on disk, maintained across transactions, shared by every connection, and available to the query planner as one option among several.

The approximation is what makes the two-phase design necessary. A bounding box overlap is a necessary condition for intersection, not a sufficient one. So the index phase produces a superset of the answer, and an exact phase filters it. PostGIS wires this up for you: writing ST_Intersects(a, b) produces a plan with a && b as the index condition and ST_Intersects as the filter. You get the optimisation without asking, which is why so few people ever need to write && by hand.

The quality of the approximation determines how much work the exact phase does. For compact shapes — parcels, buildings, cells — a bounding box is a good stand-in and false positives are rare. For long diagonal lines, sprawling multipolygons or a single feature covering a whole country, the box is enormous relative to the geometry, and the index selects candidates that the exact phase then mostly rejects. This is why ST_Subdivide is a real performance technique: chopping a huge polygon into pieces with small, tight boxes makes the index selective again. It is not an index setting; it is a change to how well boxes model the data.

And the planner has the final say. An index is an option, not an instruction. Postgres compares the estimated cost of an index scan against a sequential scan and picks the cheaper. On a small table, or when the query matches most rows, a sequential scan genuinely is cheaper — reading a table in order beats random access to most of its pages. This is why ANALYZE matters as much as CREATE INDEX: without accurate row estimates, the planner makes that comparison with bad numbers and can choose wrongly in either direction.

ANALYZE, concurrent builds and bulk loads

• ANALYZE after every bulk load. Statistics drive the planner's choice, and a fresh index with stale statistics is often ignored.
• CREATE INDEX CONCURRENTLY avoids locking writes, but scans the table twice, waits for running transactions and cannot run inside a transaction block.
• Drop the index before a large load, rebuild after. Maintaining it per row is far more expensive.
• SP-GiST and BRIN exist for geometry. BRIN is tiny and useful when rows are physically ordered by location; GiST is the default and the right first choice.
• ST_Subdivide before indexing when a few geometries are enormous — it makes the bounding boxes selective again.
• A geography column indexes the same way, but distances are in metres and the calculations are more expensive.
• ST_DWithin is index-aware; ST_Distance < d is not. Also true for ST_Intersects versus ST_Distance = 0.
• Functional indexes work: CREATE INDEX ON parcels USING GIST (ST_Centroid(geom)) if you always query centroids.
• Indexes go stale in bloated tables. REINDEX after heavy update or delete traffic.
• pg_stat_user_indexes.idx_scan = 0 means the index is unused — it still costs write time and disk.

Internal links

• How PostGIS stores geometry: SRID, EWKB and the typed column — what is being indexed
• Spatial indexes explained: R-trees and why spatial joins are fast — the same idea in memory
• My PostGIS spatial query is slow — diagnosing when the index does not help
• How to do a spatial join in PostGIS with SQL — the operation that gains the most
• Spatial SQL or GeoPandas? — when the index is the deciding factor
• How to write a GeoDataFrame to PostGIS — load first, index second
• How to use the spatial index directly in GeoPandas (sindex) — the in-memory equivalent
• Spatial predicates explained — which predicates the two-phase filter supports

FAQ

Why did nothing get faster after I created the index?

Usually a missing ANALYZE, or a predicate the planner cannot use — ST_Distance(geom, x) < d instead of ST_DWithin(geom, x, d). Check with EXPLAIN for an Index Cond line.

What does the index actually store?

Bounding boxes, not geometries. That is why it is fast, and why an exact geometry test always runs afterwards on the candidates it returns.

What is && in a query plan?

The bounding-box overlap operator — the only spatial question the index can answer. PostGIS adds it automatically when you write ST_Intersects and friends.

Should every geometry column have a GiST index?

Effectively yes. On all of Leeds the index was 12 MB of a 130 MB table-plus-index, and it made the flood-zone join about 150 times faster.

Why is my spatial join still slow with indexes on both tables?

Often the geometries are large and their bounding boxes are poor approximations, so the exact phase does most of the work. ST_Subdivide the big geometries and reindex.

Does the index work for ST_Distance?

No. Use ST_DWithin(a, b, distance), which is index-aware and returns the same rows.

Should I index before or after loading data?

After. Building once at the end is much cheaper than maintaining the index through a bulk insert — and to_postgis(if_exists="replace") drops it anyway.