Spatial Database Indexing & Performance

A zoning platform can be perfectly correct and still be unusable. The geometry is valid, the coordinate systems align, the taxonomy is clean — and then an analyst draws a study area over three census tracts, the overlay query touches a district table of two million multipolygons, and the request hangs for four minutes before the browser gives up. The data was never the problem; the access path was. When PostgreSQL cannot prove that a spatial predicate is selective, it falls back to a sequential scan and runs ST_Intersects against every row in the table, calling into GEOS millions of times for a result that should have touched a few hundred candidates. This topic sits inside Municipal Zoning Data Architecture & Compliance Frameworks and covers the discipline that keeps parcel and zoning queries fast as a dataset grows from one city to an entire state: choosing the right index type, shrinking pathological geometries, writing predicates the planner can actually accelerate, and reading the plan to prove it worked.

Prerequisites and operational context jump to heading

Indexing is the last mile of a spatial platform, not the first. The techniques below assume the surrounding architecture is already in place, and they misfire badly when it is not:

• A stable storage geometry type and SRID. Every geometry in a table must share one SRID and, ideally, one geometry type. A GiST index over a geometry(MultiPolygon, 4326) column behaves predictably; an index over a mixed geometry column that silently holds SRID 0 rows will produce candidate sets that skip half the table. This flows directly from the conventions set in municipal data structures.
• Geometries that survive validation. An invalid self-intersecting polygon can make ST_Intersects throw a TopologyException mid-query, aborting the whole statement. Run every feature through the schema validation and data quality checks layer — including ST_IsValid — before it lands, so the index is built over clean geometry.
• A working CRS in projected units where distance matters. Index acceleration is orthogonal to projection, but ST_DWithin radius queries and area computations are meaningless in degrees. Keep a projected copy (or use geography) so distance predicates are both correct and index-eligible.
• Realistic data volume to test against. A GiST index on ten thousand parcels looks fast no matter what you do. Load a full county — hundreds of thousands to millions of rows — before you trust any timing, because the planner’s cost model only starts preferring the index once the table is large enough that a sequential scan is genuinely expensive.

Get those right and indexing becomes a tuning exercise. Get them wrong and no index will save a query, because the planner will keep choosing a sequential scan for reasons that have nothing to do with the index itself.

Architecture: the two-phase filter and the index types that serve it jump to heading

Every fast spatial query in PostGIS is really two queries stacked together. The first phase is a cheap bounding-box test: does the minimum bounding rectangle of candidate geometry A overlap the bounding rectangle of query geometry B? That question is answered by the && operator, and a spatial index exists precisely to answer it without reading every row. The second phase is the exact predicate — ST_Intersects, ST_Contains, ST_Within — which runs the real, expensive GEOS computation but only on the small candidate set the first phase produced. The entire art of spatial performance is keeping phase one selective so phase two runs a few hundred times instead of a few million.

PostGIS ships three index access methods, and choosing between them is the first architectural decision:

• GiST (Generalized Search Tree) builds a balanced R-tree of bounding boxes. It is the default and the correct choice for almost all polygon and line work — parcels, zoning districts, overlay boundaries, rights-of-way. It handles arbitrary overlapping geometry, updates cleanly, and its bounding-box entries are exactly what the && operator consumes.
• SP-GiST (Space-Partitioned GiST) uses a quadtree/k-d-tree partitioning of space rather than overlapping rectangles. It can outperform GiST for point clouds and non-overlapping data where a partition boundary cleanly separates features — think address points or sensor locations — but it is a poor fit for large overlapping polygons.
• BRIN (Block Range Index) stores only the bounding box of each physical block range, making it tiny and nearly free to maintain. It is only useful when rows are physically sorted so that spatial locality matches on-disk order — a table clustered on a geohash or Hilbert curve. On unsorted parcel data it is worse than no index because its block summaries overlap everything.

The second architectural lever is geometry size. R-tree performance collapses when a single geometry has a huge bounding box, because that one box overlaps every query window and the index can never rule it out. A statewide “unincorporated area” overlay stored as one multipolygon spanning the whole state has a bounding box the size of the state; every query intersects it in phase one and then pays for an enormous exact intersection in phase two. The fix is ST_Subdivide, which splits any geometry into pieces with a bounded vertex count. Each fragment carries a tight bounding box, so the index regains its ability to discriminate, and the exact predicate operates on small pieces. This preprocessing is what makes the downstream spatial overlay analysis that consumes this table feasible at all.

The third lever is the planner’s knowledge of your data. PostgreSQL will not use an index it believes is unhelpful, and that belief comes from the statistics ANALYZE collects. On a freshly bulk-loaded parcel table the statistics are stale or absent, the planner estimates selectivity wrongly, and it chooses a sequential scan even though a perfect index sits right there. Running ANALYZE after every bulk load is not optional maintenance; it is part of making the index usable.

Production implementation jump to heading

The harness below is what a production ingestion job runs after loading a county’s parcels and the state’s overlay districts. It creates a GiST index, subdivides the oversized overlay geometries into an indexed working table, runs an overlay query wrapped in EXPLAIN (ANALYZE, BUFFERS), and then confirms the index was actually used by reading pg_stat_user_indexes. Every database call is wrapped so a single failure surfaces with context instead of a bare traceback.

"""Build spatial indexes for a county parcel/overlay workload and verify them."""
import logging
from contextlib import contextmanager
from typing import Optional

import psycopg
from psycopg import sql
from psycopg.rows import dict_row

logging.basicConfig(level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s")
logger = logging.getLogger("spatial_indexer")

DSN = "postgresql://gis:gis@localhost:5432/zoning"


@contextmanager
def connect(dsn: str = DSN):
    """Yield an autocommit connection; DDL like CREATE INDEX must not sit in an
    aborted transaction, so we commit each statement explicitly."""
    conn = None
    try:
        conn = psycopg.connect(dsn, row_factory=dict_row, autocommit=True)
        yield conn
    except psycopg.OperationalError as exc:
        logger.error("Cannot reach PostGIS at %s: %s", dsn, exc)
        raise
    finally:
        if conn is not None:
            conn.close()


def create_gist_index(conn, table: str, geom_col: str = "geom") -> None:
    """Create a GiST index if absent. CONCURRENTLY avoids locking writers on a
    live table; it cannot run inside a transaction block, hence autocommit."""
    idx = f"{table}_{geom_col}_gist"
    stmt = sql.SQL(
        "CREATE INDEX CONCURRENTLY IF NOT EXISTS {idx} ON {tbl} USING GIST ({col})"
    ).format(idx=sql.Identifier(idx), tbl=sql.Identifier(table), col=sql.Identifier(geom_col))
    try:
        conn.execute(stmt)
        logger.info("GiST index ready: %s", idx)
    except psycopg.errors.DuplicateTable:
        logger.info("Index %s already exists; skipping.", idx)
    except psycopg.Error as exc:
        # A failed CONCURRENTLY build leaves an INVALID index behind — flag it.
        logger.error("Index build failed on %s(%s): %s", table, geom_col, exc)
        raise


def subdivide_overlays(conn, source: str, target: str, max_vertices: int = 256) -> int:
    """Split oversized overlay multipolygons into an indexed working table so no
    single bounding box swamps the R-tree. Returns the fragment count."""
    stmt = sql.SQL("""
        DROP TABLE IF EXISTS {target};
        CREATE TABLE {target} AS
        SELECT district_id, zoning_code,
               ST_Subdivide(geom, %s)::geometry(Polygon, 4326) AS geom
        FROM {source}
        WHERE ST_IsValid(geom);
    """).format(target=sql.Identifier(target), source=sql.Identifier(source))
    try:
        conn.execute(stmt, (max_vertices,))
        row = conn.execute(
            sql.SQL("SELECT count(*) AS n FROM {t}").format(t=sql.Identifier(target))
        ).fetchone()
        n = row["n"]
        logger.info("Subdivided %s -> %s (%d fragments)", source, target, n)
        return n
    except psycopg.errors.InternalError_ as exc:
        # ST_Subdivide raises on some degenerate rings; the WHERE ST_IsValid
        # guard catches most, but log the survivors for the quarantine review.
        logger.error("Subdivision failed for %s: %s", source, exc)
        raise


def explain_overlay(conn, parcels: str, overlays: str) -> str:
    """Run the parcel/overlay join under EXPLAIN ANALYZE and return the plan text."""
    query = sql.SQL("""
        EXPLAIN (ANALYZE, BUFFERS, TIMING)
        SELECT p.parcel_id, o.zoning_code
        FROM {parcels} AS p
        JOIN {overlays} AS o
          ON p.geom && o.geom            -- phase 1: index-accelerated bbox test
         AND ST_Intersects(p.geom, o.geom)  -- phase 2: exact GEOS predicate
    """).format(parcels=sql.Identifier(parcels), overlays=sql.Identifier(overlays))
    try:
        rows = conn.execute(query).fetchall()
        plan = "\n".join(r["QUERY PLAN"] for r in rows)
        if "Seq Scan" in plan and "Index" not in plan:
            logger.warning("Plan still shows a sequential scan — run ANALYZE and recheck.")
        return plan
    except psycopg.errors.QueryCanceled as exc:
        logger.error("Overlay query exceeded statement_timeout: %s", exc)
        raise


def index_was_used(conn, table: str, geom_col: str = "geom") -> Optional[int]:
    """Read scan counters from pg_stat_user_indexes to prove the index earned its keep."""
    idx = f"{table}_{geom_col}_gist"
    row = conn.execute(
        "SELECT idx_scan FROM pg_stat_user_indexes WHERE indexrelname = %s", (idx,)
    ).fetchone()
    if row is None:
        logger.warning("No stats row for %s; index may not exist.", idx)
        return None
    logger.info("%s has served %d index scans.", idx, row["idx_scan"])
    return row["idx_scan"]


def main() -> None:
    with connect() as conn:
        # 1. Tame the oversized overlay geometries first.
        subdivide_overlays(conn, source="overlay_districts", target="overlay_districts_sub")
        # 2. Index both sides of the join.
        create_gist_index(conn, "parcels")
        create_gist_index(conn, "overlay_districts_sub")
        # 3. Refresh planner statistics so the index becomes eligible.
        conn.execute("ANALYZE parcels")
        conn.execute("ANALYZE overlay_districts_sub")
        # 4. Prove it: read the plan, then read the scan counters.
        plan = explain_overlay(conn, "parcels", "overlay_districts_sub")
        logger.info("Plan:\n%s", plan)
        index_was_used(conn, "overlay_districts_sub")


if __name__ == "__main__":
    try:
        main()
    except psycopg.Error:
        logger.exception("Indexing run aborted.")
        raise SystemExit(1)

Three decisions in that harness are deliberate. CREATE INDEX CONCURRENTLY runs against a live table without an ACCESS EXCLUSIVE lock, so ingestion can keep writing — at the cost of requiring autocommit and leaving an INVALID index if it fails, which the error branch flags for cleanup. ST_Subdivide runs before indexing so the fragments the index sees already have tight bounding boxes. And the run does not trust that the index helped; it reads pg_stat_user_indexes.idx_scan to confirm the counter actually moved, because a plan that “looks right” and an index that was truly used are two different claims.

Edge cases and gotchas jump to heading

Spatial indexes fail in ways that are invisible until the table gets large. Watch for these:

• The oversized-geometry trap. One statewide multipolygon in the overlay table silently defeats the index for every query, because its bounding box overlaps everything. Symptoms look like a planner problem but are a data-shape problem; ST_Subdivide is the fix, not more tuning.
• ST_Intersects without the && companion — usually fine, occasionally not. ST_Intersects internally emits an && and can use the index. But ST_Distance, ST_Contains on the wrong argument order, and buffer-based predicates often do not. When in doubt, write the && explicitly, or use ST_DWithin (which is index-aware) instead of ST_Distance(...) < r.
• Stale statistics after a bulk load. The single most common “the index isn’t working” report is a missing ANALYZE. A 2-million-row load with no fresh statistics leaves the planner estimating a handful of rows and choosing a nested loop over a sequential scan of the other table.
• Functional predicates that hide the column. ST_Intersects(ST_Transform(geom, 3857), window) cannot use an index on geom, because the indexed expression is geom, not ST_Transform(geom, 3857). Either transform the query window into the stored SRID, or build a functional index on the exact expression.
• Index bloat after heavy churn. Zoning tables that see frequent boundary corrections accumulate dead index entries. A periodic REINDEX CONCURRENTLY or VACUUM keeps the R-tree from ballooning and the plans from drifting.
• CLUSTER is a one-time sort, not a maintained order. CLUSTER parcels USING parcels_geom_gist physically reorders rows for locality and helps BRIN and buffer hit rates, but new inserts are not kept in order. Treat it as a batch operation after large loads, and re-run it if you adopt a BRIN strategy.

Index type reference jump to heading

Integration points jump to heading

Indexing does not stand alone — it is the performance contract between the tables above it and the analytics below it. Upstream, the geometry it indexes must already conform to the canonical types and SRIDs defined in municipal data structures; an index over a column with inconsistent SRIDs or mixed geometry types produces candidate sets you cannot trust. Alongside that, the schema validation and data quality checks layer is what guarantees every indexed geometry passes ST_IsValid, so the exact-predicate phase never aborts a query with a topology exception. Downstream, the whole reason these indexes exist is to make spatial overlay analysis — computing which parcels a zoning change actually touches — run in seconds rather than minutes, because that consumer issues exactly the parcel-versus-district joins this topic accelerates. The contract in both directions is a table whose geometry column is clean, single-SRID, and covered by a fresh, used index.

Compliance and audit artifacts jump to heading

For a spatial platform that underwrites decisions, a query result is only defensible if the access path that produced it is reproducible. Capture, per indexed workload:

• The index definitions themselves, dumped from pg_indexes, versioned in migration history so the exact index configuration behind any historical result is knowable.
• A representative EXPLAIN (ANALYZE, BUFFERS) plan for each critical query, stored with the schema version, so a reviewer can confirm the query used the index and see the row estimates the planner worked from.
• ANALYZE timestamps and statistics targets, because a result computed against stale statistics may have run a different plan than the same query does today.
• ST_Subdivide parameters and fragment counts, so the preprocessing that shaped every overlay bounding box is recorded rather than folded invisibly into a working table.

These pair with the lineage recorded upstream to give a continuous, provable chain from a raw municipal geometry to a fast, correct spatial answer.

FAQ jump to heading

Related jump to heading

• Tuning PostGIS GiST indexes for parcel overlay queries — a step-by-step rescue of a join stuck on a sequential scan
• Municipal data structures — the canonical geometry types and SRIDs the index depends on
• Schema validation & data quality checks — guarantees the geometry an index covers is valid
• Spatial overlay analysis — the downstream consumer whose joins these indexes make fast

Up: Municipal Zoning Data Architecture & Compliance Frameworks