Spatial Indexing in PostGIS with GiST

A spatial query over an unindexed PostGIS table scans every row; with a GiST index it jumps straight to candidates. This guide shows how to build, verify, and reason about GiST indexes so spatial predicates actually use them, driven from Python with SQLAlchemy 2.0. It is for anyone whose PostGIS queries are slower than expected. It sits under PostGIS Integration with Python in Spatial Analysis & Advanced Query Techniques.

Two-phase spatial predicate: GiST index scan versus sequential scan A query window enters both paths. With a GiST index, the bounding-box operator uses an index scan to narrow millions of rows down to a handful of candidates whose envelopes overlap, and only those candidates run the expensive exact ST_Intersects geometry test before reaching the result set. Without an index, a sequential scan reads every row and runs the exact test on all N geometries to reach the same result, far more slowly. Two-phase spatial predicate · index scan vs sequential scan With GiST index Query window && envelope GiST index scan bbox && filter Candidate rows a handful ST_Intersects exact geometry test Result set Phase 1 tree walk narrows millions of rows to the few whose boxes overlap → Phase 2 tests only those Without an index Query window && envelope Seq Scan read every row ST_Intersects on all N rows exact test, every geometry Result set No tree to descend: the expensive exact test runs on every geometry in the table Goal: EXPLAIN shows Index Scan or Bitmap Index Scan A Seq Scan on a large table means the index is missing, ANALYZE never ran, or the predicate wraps the indexed column
A spatial predicate always runs in two phases; the GiST index only accelerates the first, but that is what keeps the expensive exact test off millions of rows.

Why This Approach / What Goes Wrong

PostGIS spatial predicates run in two phases. First a fast bounding-box filter — the && operator — compares the axis-aligned envelope of each geometry against the query envelope; this is the phase a GiST index accelerates. Only rows that survive the box test proceed to the second phase: the exact geometry computation (ST_Intersects, ST_Contains, ST_DWithin, and friends), which is far more expensive because it walks vertices rather than comparing four numbers.

A GiST (Generalized Search Tree) index stores those bounding boxes in a balanced tree, so the planner can descend to the handful of rows whose envelopes overlap the query instead of box-testing all of them. Without the index, the planner has no choice but a sequential scan: it reads every row and box-tests it inline. On a few thousand rows that is fine; on a few million it is the difference between milliseconds and minutes.

A GiST descent pruning two of three subtrees The query envelope enters a root page that holds three child envelopes. Envelope A and envelope C do not overlap the query, so their whole subtrees are pruned and the millions of rows beneath them are never read. Envelope B overlaps, so the descent continues into its leaf entries, where two entries are bounding-box hits that go on to the exact geometry test and one is a box miss that is discarded. A footer notes that tree depth grows logarithmically, so doubling the table adds roughly one page read rather than doubling the work. Inside phase 1: the tree walk that never touches most of the table Query envelope the && operand Root page three child envelopes Envelope A no overlap · pruned Envelope B overlaps · descend Envelope C no overlap · pruned rows below are never read rows below are never read leaf · box hit exact test next leaf · box hit exact test next leaf · box miss discarded Phase 1 stops here — only surviving leaves are handed to ST_Intersects Depth grows with the logarithm of the row count, so doubling the table adds about one page read, not twice the work.
The tree walk is what makes the box phase cheap: whole subtrees whose envelopes miss the query are discarded before a single geometry is read.

Three mistakes stop the index from doing its job:

The goal is an EXPLAIN plan that shows an Index Scan or Bitmap Index Scan, not a Seq Scan.

Prerequisites

conda install -c conda-forge "sqlalchemy=2.0.*" "psycopg=3.1.*" "geopandas=0.14.*"

Step-by-Step Implementation

1. Create the GiST index on the geometry column. Use IF NOT EXISTS so the step is idempotent, and give the index a predictable name so you can spot it in EXPLAIN output.

from sqlalchemy import create_engine, text

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

with engine.begin() as conn:                         # begin() commits on exit
    conn.execute(text(
        "CREATE INDEX IF NOT EXISTS roads_gix ON roads USING GIST (geometry)"
    ))

For a table already carrying live traffic, build the index without taking a write lock:

# CREATE INDEX CONCURRENTLY cannot run inside a transaction block,
# so use an AUTOCOMMIT connection rather than engine.begin().
with engine.connect().execution_options(isolation_level="AUTOCOMMIT") as conn:
    conn.execute(text(
        "CREATE INDEX CONCURRENTLY IF NOT EXISTS roads_gix "
        "ON roads USING GIST (geometry)"
    ))

On a large table the build itself is worth tuning. GiST index creation is memory-bound, and the default maintenance_work_mem of 64 MB forces PostgreSQL to spill while assembling the tree. Raising it for the session routinely halves build time on a multi-million-row layer:

with engine.begin() as conn:
    conn.execute(text("SET maintenance_work_mem = '1GB'"))   # session-scoped, not global
    conn.execute(text("CREATE INDEX IF NOT EXISTS roads_gix ON roads USING GIST (geometry)"))

There is also a version cliff here. PostGIS 3.2 on PostgreSQL 14 or newer builds GiST indexes on geometry by sorting along a space-filling curve instead of inserting rows one at a time, which is several times faster and produces a better-packed tree. If an index build on a big table is taking an hour on an older stack, the upgrade is a bigger win than any parameter.

2. Run ANALYZE so the planner gets row and selectivity statistics. PostGIS ships a custom estimator for spatial columns, but it only has data to work with after ANALYZE samples the table.

with engine.begin() as conn:
    conn.execute(text("ANALYZE roads"))

The estimator builds a two-dimensional histogram of the column's extent and uses it to guess how many rows a given query window will return. That guess is what decides between an index scan and a sequential scan, so it degrades in two predictable situations: when the data is extremely clustered (a national table where 90% of rows sit in three cities) and when the extent has outliers (one stray geometry at 0°, 0° stretches the histogram over the whole planet and makes every real window look tiny). Deleting null-island rows before indexing is not tidiness — it measurably improves plan quality. For a heavily skewed column you can also raise the sample size with ALTER TABLE roads ALTER COLUMN geometry SET STATISTICS 1000; and re-run ANALYZE.

3. Write predicates the index can use — keep the indexed column bare and the SRIDs matched. The && operator and index-aware functions such as ST_Intersects and ST_DWithin recheck the exact geometry internally but let the planner reach for the GiST index for the box phase.

import geopandas as gpd

# ST_Intersects hits the GiST index for the bbox phase, then does the exact test.
# ST_MakeEnvelope's last argument is the SRID — it MUST match the column's SRID.
window_sql = """
SELECT r.road_id, r.geometry
FROM roads r
WHERE ST_Intersects(
    r.geometry,
    ST_MakeEnvelope(7.6, 45.0, 7.8, 45.1, 4326)
)
"""
window = gpd.read_postgis(window_sql, engine, geom_col="geometry")

4. Avoid index-defeating patterns. The single most common way to lose the index is to transform the indexed column. Buffer or transform the constant side instead, or switch to a function that is index-aware by construction.

# BAD: ST_Buffer wraps the indexed column, so the GiST index can't be used.
#   WHERE ST_Intersects(ST_Buffer(r.geometry, 250), :depot)
#
# GOOD: ST_DWithin is index-aware — it uses the index and applies the
# distance internally, so nothing wraps r.geometry.
#   WHERE ST_DWithin(r.geometry, :depot, 250)

Note the distance unit: ST_DWithin measures in the column's own units. On a projected column (metres) 250 means 250 m; on a geographic geometry column in EPSG:4326 it means 250 degrees, which is almost never what you want. For metric radius searches, either store the column in a projected CRS such as the appropriate UTM zone for your data or cast to geography (ST_DWithin(geometry::geography, :depot, 250)), which measures in metres on the ellipsoid and still uses a GiST index. This is exactly the pattern behind nearest-neighbour and radius search inside the database.

Which predicate shapes a GiST index can serve A six-row matrix comparing predicate forms. The double-ampersand envelope test, ST_Intersects against a constant envelope, ST_DWithin on a bare column, and ST_DWithin on a geography cast all keep the indexed column unwrapped and are served by the index; the distance forms measure in the column's own units, which means degrees on a geographic column and ellipsoidal metres after a geography cast. Wrapping the column in ST_Buffer or in ST_Transform makes the expression uncomputable from the stored envelopes, so the planner falls back to a sequential scan with no error raised. Keep the indexed column bare — everything else is negotiable Predicate written in the WHERE clause Planner Distance unit Why geometry && ST_MakeEnvelope(...) index scan n/a Pure box test — the cheapest filter there is ST_Intersects(geometry, envelope) index scan n/a Box phase indexed, then the exact test ST_DWithin(geometry, depot, 250) index scan column units Index-aware, but 250 means degrees on 4326 ST_DWithin(geometry::geography, depot, 250) index scan metres Ellipsoidal metres and still indexed ST_Intersects(ST_Buffer(geom, 250), depot) seq scan column units A function wraps the column, so no index ST_Intersects(ST_Transform(geom, 25832), win) seq scan metres Reprojects every row — transform the window Neither of the last two rows raises an error; the only symptom is a plan that says Seq Scan.
The stored envelopes describe the column as it sits on disk, so any function applied to that column makes the index unusable — transform or buffer the constant side instead.

5. Reach for a variant index only when the plain GiST index is provably the bottleneck. PostGIS supports three operator-class families on geometry, and they are not interchangeable:

Two index shapes matter more often than the operator class. A partial index skips rows you never query, which shrinks the tree and speeds the build:

with engine.begin() as conn:
    conn.execute(text(
        "CREATE INDEX IF NOT EXISTS roads_active_gix ON roads USING GIST (geometry) "
        "WHERE status = 'active'"
    ))

The planner only uses that index for queries whose WHERE clause implies status = 'active', so the predicate has to appear in the query too. And when a spatial filter is nearly always paired with an attribute filter, a multi-column index over both columns lets one index scan do both jobs — but this needs the btree_gist extension, because a plain GiST index cannot hold a text or integer column:

CREATE EXTENSION IF NOT EXISTS btree_gist;
CREATE INDEX roads_class_geom_gix ON roads USING GIST (road_class, geometry);

Order matters: put the equality-tested column first. Measure before and after, because two separate indexes combined by a bitmap AND are frequently just as fast and far more flexible.

Verification

Use EXPLAIN ANALYZE to confirm the planner actually chose an index scan. Reading the plan from Python lets you assert on it in a test.

from sqlalchemy import create_engine, text

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

with engine.connect() as conn:
    plan = conn.execute(text(
        "EXPLAIN ANALYZE SELECT road_id FROM roads "
        "WHERE geometry && ST_MakeEnvelope(7.6, 45.0, 7.8, 45.1, 4326)"
    )).fetchall()

plan_text = "\n".join(row[0] for row in plan)
print(plan_text)
# Index Scan using roads_gix on roads  (cost=0.28..8.30 rows=1 ...)
#   Index Cond: (geometry && '0103...'::geometry)
# ...
assert "Index Scan" in plan_text or "Bitmap Index Scan" in plan_text

Both Index Scan and Bitmap Index Scan mean the index is working — Postgres picks the bitmap variant when it expects to return many rows. Seeing Seq Scan instead means the index is missing, ANALYZE was not run, or the predicate is not index-eligible.

You can also confirm the index physically exists and see how large it grew:

with engine.connect() as conn:
    info = conn.execute(text("""
        SELECT indexname, pg_size_pretty(pg_relation_size(indexrelid)) AS size
        FROM pg_indexes
        JOIN pg_class ON pg_class.relname = indexname
        JOIN pg_index ON pg_index.indexrelid = pg_class.oid
        WHERE tablename = 'roads' AND indexdef ILIKE '%USING gist%'
    """)).fetchall()
print(info)   # [('roads_gix', '12 MB')]

A one-off EXPLAIN proves the index can be used; production tells you whether it is. pg_stat_user_indexes counts every scan since the last statistics reset, so an index with idx_scan = 0 after a week of traffic is dead weight — it slows every write and occupies disk for nothing:

with engine.connect() as conn:
    usage = conn.execute(text("""
        SELECT indexrelname, idx_scan, idx_tup_read
        FROM pg_stat_user_indexes
        WHERE relname = 'roads'
        ORDER BY idx_scan DESC
    """)).fetchall()
print(usage)
# [('roads_gix', 148213, 2914558), ('roads_active_gix', 0, 0)]

Zero scans on a spatial index almost always means the application is filtering in Python after fetching, not that the index is broken — pushing the predicate into SQL is the fix, as covered in connecting GeoPandas to PostGIS.

Edge Cases & Debugging

Frequently Asked Questions

Should I create the index before or after a bulk load? After, without exception. Loading into an indexed table maintains the tree on every insert, which can triple the load time and leaves a more fragmented index than a single build over the finished table. The one caveat is that the table is unindexed and therefore slow for readers during the load — if that matters, load into a staging table, index it there, and rename it into place.

Does an index on geometry also accelerate geography queries? Only if the column and the predicate use the same type. A GiST index on a geometry column cannot serve ST_DWithin(geom::geography, ...), because the cast produces a value the stored envelopes do not describe — the plan quietly becomes a sequential scan. Either store the column as geography and index that, or add a functional index on the cast: CREATE INDEX ON sensors USING GIST ((geom::geography));.

How much disk does a GiST index cost? Roughly 10–20% of the table's size for typical point and polygon layers, since each entry stores only a bounding box and a row pointer regardless of vertex count. That means the index over a table of 40 000-vertex coastlines is proportionally tiny, while an index over a table of points can approach the size of the data. Check the real number with pg_relation_size rather than estimating.

Why does ORDER BY geom <-> point LIMIT 5 beat ORDER BY ST_Distance(...) LIMIT 5? The <-> operator is index-aware: GiST walks the tree in order of increasing distance and stops after five rows. ST_Distance is an ordinary function, so the planner has to compute it for every row and sort the lot. On a million-row table that is the difference between a millisecond and several seconds — the same distinction between a bounded and an unbounded search that drives the choices in Nearest Neighbor & KD-Tree Search.

Related