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

Problem statement #

The same query, twice, on the same hardware:

SELECT COUNT(*)
FROM parcels p JOIN flood_zones f ON ST_Intersects(p.geom, f.geom);

Time: 412803.117 ms   (6 minutes 53 seconds)

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

Time: 2841.402 ms   (2.8 seconds)

145 times faster, and nothing about the data or the query changed. Meanwhile this one stays slow no matter how many indexes you build:

SELECT * FROM parcels WHERE ST_Distance(geom, %(point)s) < 500;   -- still minutes

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.

Step-by-step solution #

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 — four floats — 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 #

EXPLAIN ANALYZE
SELECT COUNT(*) FROM parcels p
JOIN flood_zones f ON ST_Intersects(p.geom, f.geom);

Nested Loop  (cost=0.42..184203.11 rows=41288 width=8)
  ->  Seq Scan on flood_zones f  (rows=1204)
  ->  Index Scan using parcels_geom_idx on parcels p  (rows=34)
        Index Cond: (geom && f.geom)
        Filter: st_intersects(geom, f.geom)
        Rows Removed by Filter: 11

Two lines carry the whole story:

• Index Cond: (geom && f.geom) — the index phase. && means "bounding boxes overlap". This is what the tree can answer.
• Filter: st_intersects(geom, f.geom) — the exact phase. Real geometry comparison, run only on the survivors.

Rows Removed by Filter: 11 is the false-positive count: 45 candidates per zone passed the box test, 34 were real intersections, 11 were boxes overlapping without geometries touching. PostGIS adds the && condition automatically — you write ST_Intersects and get the two-phase plan for free.

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.

-- ❌ a function of the column, not the column
WHERE ST_Distance(geom, %(pt)s) < 500

-- ✅ the same rows, as a box question
WHERE ST_DWithin(geom, %(pt)s, 500)

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 same rule explains the other cases:

-- ❌ no spatial index can help; consider a plain B-tree on a stored area column
WHERE ST_Area(geom) > 10000

-- ❌ transforms every row before comparing
WHERE ST_Transform(geom, 4326) && %(box4326)s

-- ✅ transform the constant instead, and the column stays bare
WHERE geom && ST_Transform(%(box4326)s, 27700)

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:

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

CONCURRENTLY takes roughly twice as long and cannot run inside a transaction block, but it does not block writes. On a table being loaded, do the opposite — load first, index afterwards:

DROP INDEX IF EXISTS parcels_geom_idx;
-- bulk load here
CREATE INDEX parcels_geom_idx ON parcels USING GIST (geom);
ANALYZE parcels;

Maintaining an index during a million-row insert costs far more than rebuilding it once at the end.

5. Check that it is actually being used #

SELECT relname, indexrelname, idx_scan, idx_tup_read
FROM pg_stat_user_indexes
WHERE indexrelname LIKE '%geom%'
ORDER BY idx_scan;

 relname     | indexrelname          | idx_scan | idx_tup_read
-------------+-----------------------+----------+--------------
 roads       | roads_geom_idx        |        0 |            0
 parcels     | parcels_geom_idx      |    18402 |     41288104

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.

Code examples #

Example 1: measuring the index from Python #

import time
import geopandas as gpd
from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://user@localhost/gis")

QUERY = """
SELECT COUNT(*) FROM parcels p
JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
"""

def time_query(sql, label, use_index=True):
    with engine.begin() as con:
        con.execute(text(f"SET enable_indexscan = {use_index}"))
        con.execute(text(f"SET 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(QUERY, "without index", use_index=False)
fast = time_query(QUERY, "with index",    use_index=True)
print(f"speedup: {slow / fast:.0f}×")

without index            41,288 rows    412.80 s
with index               41,288 rows      2.84 s
speedup: 145×

SET enable_indexscan = off is a session-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
from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://user@localhost/gis")

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"]
        if kind == "Seq Scan":
            seq_scans.append((node["Relation Name"], node.get("Actual Rows", 0)))
        elif "Index" in kind:
            index_scans.append((node.get("Index Name"), node.get("Actual Rows", 0)))
        removed += node.get("Rows Removed by Filter", 0)
        for child in node.get("Plans", []):
            walk(child)
    walk(plan)

    print(f"total time     {plan['Actual Total Time']:.1f} ms")
    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

explain("""
    SELECT COUNT(*) FROM parcels p
    JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
""")

total time     2841.4 ms
  index scan   parcels_geom_idx → 41,288 rows
  seq scan     flood_zones → 1,204 rows
  false positives removed by exact test: 13,244

Two readings matter. A sequential scan over 1,204 rows is fine — with a table that small the planner is right to skip the index. A sequential scan over 4 million rows in the same position is the finding.

And 13,244 false positives against 41,288 real hits is a healthy ratio of about 1:3. A ratio of 20:1 would say the bounding boxes are poor approximations of the geometries — long diagonal lines, sprawling multipolygons — and that the exact phase is doing most of the work. The fix there is subdividing the geometries, not another index.

Example 3: indexing a table you are loading from Python #

import geopandas as gpd
from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://user@localhost/gis")

def load_with_index(gdf, table, engine, srid=27700, chunksize=10_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'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(gpd.read_file("parcels.gpkg"), "parcels", engine)

parcels: 4,012,884 rows, 6284 MB total, 172 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.

The index is 172 MB against 6.3 GB of table: under 3%, for a 145× improvement on the queries that use it. That ratio is typical, and it is why "should I index this geometry column?" is almost always yes. Loading details are in how to write a GeoDataFrame to PostGIS.

Explanation #

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.

Edge cases or notes #

• 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 takes about twice as long and cannot run inside a transaction.
• 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. It typically costs under 3% of the table size and can make spatial queries a hundred 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.