Tuning GiST Indexes for Large-Scale Spatial Joins

This page solves one exact problem: an ST_Intersects or ST_DWithin join between two multi-million-row spatial tables that runs in minutes instead of milliseconds because the GiST index is missing, misbuilt, or silently ignored by the planner.

Index tuning is where the location intelligence architecture and data foundations layer earns its keep: a schema is only as fast as the indexes the planner can actually reach. A spatial join is a two-phase operation. PostGIS first runs a cheap bounding-box filter using the index — the && operator — and only then runs the expensive exact predicate (ST_Intersects, ST_DWithin) on the surviving candidates. If the index is healthy, phase one discards 99.9% of pairs before the exact geometry math ever runs. If the index is absent or defeated, the planner falls back to a sequential scan and evaluates the exact predicate against every row pair, which is where the minutes go. Everything below is about keeping phase one fast and making sure the planner actually uses it.

Prerequisites #

Before tuning, confirm the environment:

• PostGIS 3.x on PostgreSQL 13+. The USING GIST operator class for geometry ships with PostGIS; nothing extra to install. Confirm with SELECT postgis_full_version();.
• Two spatial tables with an enforced SRID. This task assumes the schema conventions from how to structure a geospatial database for multi-state retail chains — typed geometry(Point, 5070) and geometry(MultiPolygon, 5070) columns rather than a generic geometry. A mismatched or unconstrained SRID is the single most common reason an index is skipped (see Failure Modes).
• A metric, projected CRS. Store geometry in an equal-area or equal-distance projection such as EPSG:5070 (North America Albers) so that ST_DWithin distances are in metres. Running ST_DWithin on EPSG:4326 geometry measures in degrees, which is both wrong and un-indexable in the way you expect.
• psycopg (v3) and pandas if you want to run the benchmark harness below.
• ANALYZE privileges on both tables, so the planner has current statistics.

Index Choice: GiST, SP-GiST, or BRIN #

PostGIS offers three index access methods for geometry, and they are not interchangeable. The choice depends on the geometry type and, critically, on whether the table is physically sorted in space.

The default and correct first choice for a polygon-to-point join is GiST on both sides. Reach for SP-GiST only when profiling a pure-point table shows GiST index scans dominating; reach for BRIN only when the table is naturally sorted in space (for example, ingested region by region) and you can tolerate coarser pruning in exchange for an index that is a thousand times smaller.

Build Parameters That Matter #

The CREATE INDEX statement has a handful of knobs that change build time and steady-state query speed.

CLUSTER is the most under-used lever. After CLUSTER stores USING stores_geom_gist;, rows that are near in space sit near on disk, so a join that walks a region reads far fewer pages. The catch: CLUSTER is a one-time physical sort — new inserts are not kept in order, so re-cluster periodically (a nightly or weekly maintenance job) on tables that grow.

Annotated Build and Benchmark #

The SQL below builds tuned GiST indexes on both sides of the join, clusters the larger table, refreshes statistics, and inspects the plan. The comments call out why each step exists.

-- 1. Raise build memory for this session only (reset on disconnect).
SET maintenance_work_mem = '2GB';
SET max_parallel_maintenance_workers = 4;

-- 2. Build a GiST index on each geometry column participating in the join.
--    Both sides MUST be indexed; a join uses the index on the table it scans.
CREATE INDEX IF NOT EXISTS stores_geom_gist
    ON stores USING GIST (geom)
    WITH (fillfactor = 90);            -- static reference table: pack tightly

CREATE INDEX IF NOT EXISTS trade_areas_geom_gist
    ON trade_areas USING GIST (geom)
    WITH (fillfactor = 90);

-- 3. Physically order the larger table along the index for cache locality.
--    One-time sort; schedule a periodic re-CLUSTER as the table grows.
CLUSTER stores USING stores_geom_gist;

-- 4. Refresh planner statistics AFTER the index and CLUSTER, or the planner
--    estimates selectivity from stale samples and may skip the index.
ANALYZE stores;
ANALYZE trade_areas;

-- 5. Inspect the plan. We WANT to see an "Index Scan" / "Nested Loop" with the
--    && bbox filter pushed into the index, not a "Seq Scan".
EXPLAIN (ANALYZE, BUFFERS)
SELECT t.area_id, count(s.store_id) AS competitor_count
FROM   trade_areas t
JOIN   stores s
       ON ST_DWithin(t.geom, s.geom, 3000)   -- 3 km, metres because SRID 5070
GROUP  BY t.area_id;

A correct plan for step 5 shows a Nested Loop whose inner side is an Index Scan using stores_geom_gist, with an Index Cond referencing the && bounding-box operator that ST_DWithin expands to internally. If instead you see Seq Scan on stores, the index is being ignored — jump to Failure Modes.

To measure the real-world delta rather than trust the plan alone, drive the query from Python with a small psycopg harness that times the same statement with and without the index.

import time
import statistics
import psycopg  # psycopg 3

DSN = "postgresql://analyst@localhost:5432/retail"

JOIN_SQL = """
    SELECT t.area_id, count(s.store_id) AS competitor_count
    FROM   trade_areas t
    JOIN   stores s
           ON ST_DWithin(t.geom, s.geom, 3000)
    GROUP  BY t.area_id;
"""


def assert_srid(conn, table: str, col: str, expected: int) -> None:
    """CRS discipline: refuse to benchmark if the column SRID is not what we
    expect. A mismatched SRID silently defeats the spatial index."""
    with conn.cursor() as cur:
        cur.execute(
            "SELECT DISTINCT ST_SRID(%s) FROM %s LIMIT 2;"
            % (col, table)  # table/col are trusted schema identifiers here
        )
        srids = [row[0] for row in cur.fetchall()]
    assert srids == [expected], f"{table}.{col} SRID {srids} != {expected}"


def time_query(conn, sql: str, runs: int = 5) -> float:
    """Return median wall-clock seconds over N runs (warm cache after run 1)."""
    timings = []
    with conn.cursor() as cur:
        for _ in range(runs):
            start = time.perf_counter()
            cur.execute(sql)
            cur.fetchall()
            timings.append(time.perf_counter() - start)
    return statistics.median(timings)


with psycopg.connect(DSN) as conn:
    # Guard: both join columns must share the SRID we projected to (EPSG:5070).
    assert_srid(conn, "stores", "geom", 5070)
    assert_srid(conn, "trade_areas", "geom", 5070)

    # Force a fair comparison: disable index scans, then re-enable.
    with conn.cursor() as cur:
        cur.execute("SET enable_indexscan = off; SET enable_bitmapscan = off;")
    seq = time_query(conn, JOIN_SQL)

    with conn.cursor() as cur:
        cur.execute("RESET enable_indexscan; RESET enable_bitmapscan;")
    idx = time_query(conn, JOIN_SQL)

    print(f"seq-scan median: {seq:8.3f}s")
    print(f"gist    median: {idx:8.3f}s")
    print(f"speedup        : {seq / idx:8.1f}x")

On a join between roughly two million store points and a hundred thousand trade-area polygons, a healthy GiST index typically turns a 40–90 second sequential scan into a sub-second index-driven nested loop — a two-to-three order of magnitude speedup. If the harness reports a speedup near 1x, the “index” run is not actually using the index, which points at one of the failures below.

Failure Modes and Debugging #

Two of these deserve emphasis because they are subtle. First, any function wrapping the indexed column in the predicate defeats the index: ST_DWithin(ST_Transform(s.geom, 4326), t.geom, 0.03) cannot use stores_geom_gist, because the index is built on geom, not on ST_Transform(geom, 4326). The fix is architectural — settle the join CRS once at ingest, store geometry already in that CRS, and never transform inside a hot predicate. Second, geography versus geometry is a genuine fork: geography gives correct great-circle distances globally but has different, often slower, index behaviour and forces metres regardless of projection. For continental retail work, projected geometry in an equal-distance CRS is the faster and more predictable choice.

The decision flow below shows how to diagnose a slow join from the EXPLAIN output downward.

Verification #

Confirm the tuning worked, not just that it ran:

1. The plan uses the index. EXPLAIN (ANALYZE, BUFFERS) shows Index Scan using …_gist with an && index condition and no Seq Scan on the large table.
2. Timing delta is real. The psycopg harness reports a speedup well above 1x (typically 100–1000x on multi-million-row joins). A speedup near 1x means the “indexed” run still seq-scanned.
3. Buffers dropped. In the BUFFERS line, shared read/shared hit counts fall by orders of magnitude between the seq-scan and index runs — the index reads far fewer pages.
4. Row counts are identical. The join result count must match between the two runs; an index that changes the answer signals an SRID or predicate bug, not a speedup.
5. Index is used, not bloated. SELECT pg_size_pretty(pg_relation_size('stores_geom_gist')); and pg_stat_user_indexes.idx_scan > 0 confirm the index exists, is a sane size, and is actually being hit.

-- Verification one-liner: has the index been used since the last stats reset?
SELECT relname, indexrelname, idx_scan
FROM   pg_stat_user_indexes
WHERE  indexrelname IN ('stores_geom_gist', 'trade_areas_geom_gist');

An idx_scan of zero after running the workload is the clearest possible signal that the planner is routing around your index — treat it as a failing gate and return to the diagnostic flow above.

Related #

• Setting Up PostGIS for Retail Analytics — the schema, SRID, and extension baseline this tuning assumes.
• How to Structure a Geospatial Database for Multi-State Retail Chains — table layout and SRID enforcement that keep indexes usable.
• Building a Competitor Proximity Index in PostGIS — a production ST_DWithin workload that lives or dies on the tuning here.

← Back to Setting Up PostGIS for Retail Analytics