Tuning PostGIS GiST indexes for parcel overlay queries

A nightly job joins Maricopa County’s 1.8 million parcels against a table of flood and historic-overlay districts to stamp each lot with the overlays it falls under. It used to finish in seconds on the staging extract; now, against the full county, it runs for four minutes and sometimes trips the statement timeout. Nothing about the query changed — only the row count did. This is the classic symptom of a spatial join that has fallen off its index: PostgreSQL is running ST_Intersects against every parcel-district pair because it no longer believes an index scan is cheaper. This guide is the focused rescue procedure for exactly that failure, a companion to the broader spatial database indexing and performance topic. We will read the plan, prove the seq scan, build and repair the GiST path, and confirm the timing collapses back to interactive speed.

Diagnosis: proving the sequential scan

Never tune from a hunch. Capture the plan first, because “it feels slow” and “it is running a sequential scan” are different problems with different fixes. Wrap the offending join in EXPLAIN (ANALYZE, BUFFERS):

EXPLAIN (ANALYZE, BUFFERS)
SELECT p.parcel_id, o.overlay_code
FROM parcels AS p
JOIN overlay_districts AS o
  ON ST_Intersects(p.geom, o.geom);

A healthy plan names an Index Scan or Bitmap Index Scan on a GiST index. The broken plan looks like this:

Nested Loop  (cost=0.00..48211900.55 rows=1823 width=41)
              (actual time=0.412..238190.771 rows=2144067 loops=1)
  Join Filter: st_intersects(p.geom, o.geom)
  Rows Removed by Join Filter: 3271884933
  ->  Seq Scan on overlay_districts o
        (actual time=0.020..3.114 rows=1794 loops=1)
  ->  Seq Scan on parcels p
        (actual time=0.004..96.517 rows=1823000 loops=1)
Planning Time: 0.203 ms
Execution Time: 238400.912 ms

Two lines are the smoking gun. Seq Scan on parcels means the index — if one even exists — is not being used, and Rows Removed by Join Filter: 3271884933 means GEOS ran ST_Intersects more than three billion times. That is the two-phase filter collapsed into one phase: with no usable bounding-box prefilter, every parcel is tested against every district. The diagram below maps the decision path from this plan back to a fast one.

Decision flow from a sequential-scan overlay join to an index scan EXPLAIN ANALYZE reveals a sequential scan with billions of rows removed by the join filter. The remedy runs left to right: build a GiST index on the parcel geometry, subdivide oversized overlay district polygons with ST_Subdivide, rewrite the join to use the bounding-box overlap operator plus ST_Intersects, then run ANALYZE to refresh statistics. A final re-check of the plan branches: if it now shows an index scan the job is fast and done; if it still shows a sequential scan, the flow loops back to the ANALYZE and rewrite steps. EXPLAIN ANALYZE Seq Scan · 3.2B rows removed by filter Build GiST index USING GIST (geom) on both tables ST_Subdivide split oversized overlay polygons Rewrite join && + ST_Intersects two-phase filter ANALYZE refresh stats planner sees rows Re-check plan EXPLAIN ANALYZE again index scan Done · seconds timing collapses still seq scan → re-ANALYZE, recheck predicate From a four-minute sequential scan to a sub-second index scan

Step-by-step implementation

Each step fixes one link in the chain. Run them in order and re-measure after the last.

Step 1 — Build a GiST index on both join sides

The parcel table almost certainly has no spatial index, or has one on the wrong column. Build a GiST index on the geometry of both tables. Use CONCURRENTLY so the nightly writers are not locked out, and confirm the geometry column is a typed, single-SRID column first.

-- Confirm the column is what the index expects before building.
SELECT ST_SRID(geom) AS srid, GeometryType(geom) AS gtype, count(*)
FROM parcels GROUP BY 1, 2;

-- Build the index without an ACCESS EXCLUSIVE lock on the live table.
CREATE INDEX CONCURRENTLY IF NOT EXISTS parcels_geom_gist
  ON parcels USING GIST (geom);

CREATE INDEX CONCURRENTLY IF NOT EXISTS overlay_districts_geom_gist
  ON overlay_districts USING GIST (geom);

If the SELECT returns more than one SRID, stop — an index over mixed SRIDs yields wrong candidate sets, and you must normalize the column before indexing.

Step 2 — Subdivide the oversized district polygons

Overlay districts are the usual culprit: a county-wide floodplain stored as a single multipolygon has a bounding box the size of the county, so the R-tree can never rule it out and every parcel still becomes a candidate. Split it so each fragment carries a tight bounding box.

-- Materialize subdivided overlays; cap each fragment at 256 vertices.
DROP TABLE IF EXISTS overlay_districts_sub;
CREATE TABLE overlay_districts_sub AS
SELECT district_id,
       overlay_code,
       ST_Subdivide(geom, 256)::geometry(Polygon, 4326) AS geom
FROM overlay_districts
WHERE ST_IsValid(geom);   -- skip invalid rings that ST_Subdivide would reject

CREATE INDEX CONCURRENTLY overlay_districts_sub_geom_gist
  ON overlay_districts_sub USING GIST (geom);

The join now targets overlay_districts_sub. One 40,000-vertex floodplain becomes a few hundred small polygons whose bounding boxes actually discriminate.

Step 3 — Rewrite the join for the two-phase filter

Make the bounding-box prefilter explicit. The && operator is what the GiST index accelerates; ST_Intersects then runs the exact test only on the survivors.

SELECT p.parcel_id, o.overlay_code
FROM parcels AS p
JOIN overlay_districts_sub AS o
  ON p.geom && o.geom              -- phase 1: index-accelerated bbox overlap
 AND ST_Intersects(p.geom, o.geom) -- phase 2: exact predicate on candidates
GROUP BY p.parcel_id, o.overlay_code;  -- collapse duplicate fragment hits

Because subdivision can make one district match a parcel through several fragments, GROUP BY (or SELECT DISTINCT) collapses those duplicate rows back to one overlay per parcel.

Step 4 — Run ANALYZE to refresh statistics

Even with a perfect index, the planner will keep choosing a sequential scan if its statistics predate the data. Refresh them.

ANALYZE parcels;
ANALYZE overlay_districts_sub;

ANALYZE re-samples the tables so the planner’s row estimates match reality and it starts costing the index scan correctly.

Step 5 — Re-check the plan and lock in the result

Re-run the EXPLAIN (ANALYZE, BUFFERS) from the diagnosis. Confirm the seq scan is gone and an index scan has taken its place.

import logging
import psycopg
from psycopg.rows import dict_row

logging.basicConfig(level=logging.INFO)
log = logging.getLogger("overlay_tuning")

CHECK = """
EXPLAIN (ANALYZE, BUFFERS)
SELECT p.parcel_id, o.overlay_code
FROM parcels p
JOIN overlay_districts_sub o
  ON p.geom && o.geom AND ST_Intersects(p.geom, o.geom);
"""

def verify(dsn: str) -> bool:
    """Return True only if the plan uses a spatial index scan."""
    try:
        with psycopg.connect(dsn, row_factory=dict_row, autocommit=True) as conn:
            plan = "\n".join(r["QUERY PLAN"] for r in conn.execute(CHECK).fetchall())
    except psycopg.OperationalError as exc:
        log.error("Cannot connect to verify plan: %s", exc)
        return False
    used_index = "Index Scan" in plan or "Bitmap Index Scan" in plan
    if not used_index:
        log.warning("Still no index scan:\n%s", plan)
    else:
        log.info("Index scan confirmed.")
    return used_index

if __name__ == "__main__":
    verify("postgresql://gis:gis@localhost:5432/zoning")

Verification & testing

Do not declare victory on timing alone; confirm the mechanism changed.

  • The plan names an index. The re-checked plan must show Index Scan or Bitmap Index Scan using overlay_districts_sub_geom_gist, and the Seq Scan on parcels line must be gone.
  • Rows removed by filter collapses. The Rows Removed by Join Filter count should drop from billions to a small multiple of the result size. That number is the clearest proof the prefilter is doing its job.
  • The counter moves. Query SELECT idx_scan FROM pg_stat_user_indexes WHERE indexrelname = 'parcels_geom_gist'; before and after; a real use increments it.
  • The result set is unchanged. Diff the parcel-to-overlay mapping against the slow query’s output on a small extract. Subdivision plus GROUP BY must produce identical assignments — same parcels, same overlay codes — just faster.
  • Timing lands in seconds. On a county-scale table the tuned join should complete in low single-digit seconds, not minutes.

Failure recovery

When the plan still shows a sequential scan after all five steps, work through the usual causes rather than adding more indexes.

  • Statistics never refreshed. If ANALYZE was skipped or ran before the load finished, the planner is still blind. Re-run ANALYZE, and for volatile tables set an aggressive autovacuum analyze_scale_factor so this does not recur.
  • The predicate hides the column. A join written as ST_Intersects(ST_Transform(p.geom, 3857), o.geom) cannot use an index on p.geom. Transform the other side into the stored SRID, or build a functional index on the exact transformed expression.
  • An INVALID index from a failed CONCURRENTLY build. Check SELECT indexrelid::regclass FROM pg_index WHERE NOT indisvalid;. Drop any invalid index and rebuild it; the planner ignores invalid indexes entirely.
  • The planner is mis-costing on a small dev extract. On ten thousand rows a sequential scan genuinely is cheaper, and forcing the index proves nothing. Verify against production-scale data; the index only wins once the table is large.
  • Downstream still slow despite a fast join. If the fast mapping feeds analytics, the bottleneck may have moved. The consuming spatial overlay analysis step has its own geometry-heavy operations, and a fast join into a slow overlay computation still looks slow end to end.

Frequently asked questions

Why did the query get slow only after loading the full county?

Because the planner’s cost model is data-dependent. On a small staging extract a sequential scan really is cheaper, so the plan you had was fine — it just did not use the index. Once the table grew to millions of rows the same plan became catastrophic, but nothing re-triggered a better one until you built the index and ran ANALYZE.

Do I have to keep the subdivided overlay_districts_sub table in sync?

Yes. It is a materialized derivative, so whenever the source overlay geometry changes you must re-run the ST_Subdivide step and rebuild its index. Wire it into the same job that loads the districts, or use a trigger, so the fast table never drifts from the authoritative one.

My plan uses the index but still returns duplicate overlay rows — why?

Subdivision splits one district into many fragments, and a parcel can intersect several fragments of the same district, producing one row per fragment hit. Collapse them with GROUP BY on the parcel and overlay code, or SELECT DISTINCT, to get one overlay assignment per parcel.

Should I index the original overlay table or only the subdivided one?

Index whichever table the join actually targets. Once you rewrite the join against overlay_districts_sub, that is the table that needs the GiST index. Keep an index on the original only if other queries still hit it directly; otherwise it is dead write cost.

Up: Spatial Database Indexing & Performance