Reading EXPLAIN ANALYZE Output for Spatial Joins

This page addresses the exact scenario described under Query Plan Analysis with EXPLAIN: you have a PostGIS spatial join between two geometry tables, the query is slow, and you need to read the EXPLAIN ANALYZE output to pinpoint whether the GiST index is engaged, how much work the exact geometry predicate is doing, and whether the planner’s row estimates are accurate enough to produce a sensible join order. Three deterministic signals tell you the whole story: the presence of Index Cond, the presence of Recheck Cond, and the ratio of Plan Rows to Actual Rows.

Why Reading the Plain Output Fails

The naive approach is to run EXPLAIN ANALYZE in text format and scan visually for the word “Index”:

sql
-- Naive: text format loses structure and makes comparison error-prone
EXPLAIN ANALYZE
SELECT a.id, b.id
FROM parcels a
JOIN zones b ON ST_Intersects(a.geom, b.geom);

Text output looks convincing — you see Index Scan in the label — but it hides critical detail. The Index Cond and Filter lines are on separate rows with inconsistent indentation, making it easy to misread a Filter-only node (which runs the exact geometry test against every row) as an index-accelerated scan. You also cannot programmatically compare Plan Rows to Actual Rows across nodes without fragile string parsing. When the plan involves a Bitmap Heap Scan with Recheck Cond, the cost of bounding-box re-verification is invisible in aggregate timing.

The second mistake is omitting BUFFERS. Without it you cannot distinguish whether slow Index Scan nodes are spending time on cache misses (cold storage) versus genuine selectivity problems. A plan that looks like a GiST index win may be hammering shared buffers with random I/O on an unordered heap.

sql
-- Slow: missing BUFFERS, missing FORMAT JSON
EXPLAIN (ANALYZE, COSTS, TIMING)
SELECT a.id, b.id
FROM parcels a
JOIN zones b ON ST_Intersects(a.geom, b.geom);

The Two-Phase Execution Pipeline

Before reading any plan output, understand the two-phase architecture PostGIS uses for every spatial join. The diagram below shows how a query routes through the GiST index before the exact predicate runs.

Two-phase PostGIS spatial join execution pipeline Diagram showing Phase 1: GiST index bounding-box filter (Index Cond) reducing the candidate set, followed by Phase 2: exact geometry predicate (Filter) running only on surviving candidates before rows are returned to the application. Phase 1 GiST index scan Index Cond: a.geom && b.geom work_mem lossy? yes → Recheck Cond no Phase 2 Exact geometry test Filter: ST_Intersects(a.geom, b.geom) out All table rows → candidate set Candidates → matched rows Plan Rows ↔ Actual Rows divergence detected here

The key insight from this pipeline is that Index Cond marks the boundary where the GiST index is active. Everything inside Filter runs on the heap without index acceleration. A healthy plan minimises the rows reaching Filter by making Index Cond highly selective.

Production-Ready Implementation

Use FORMAT JSON with BUFFERS for all plan capture. The snippet below runs the spatial join with a complete plan extraction and surfaces only the nodes that matter for spatial join diagnosis.

python
import psycopg2
import json
from typing import Any

def analyze_spatial_join(
    conn_string: str,
    table_a: str,
    table_b: str,
    geom_col: str = "geom",
    srid: int = 4326,
) -> list[dict[str, Any]]:
    """
    Capture and parse EXPLAIN ANALYZE output for a PostGIS ST_Intersects join.

    Returns a flat list of plan nodes with the spatial-relevant fields extracted.
    Compatible with PostgreSQL 14+ / PostGIS 3.3+.

    Args:
        conn_string: libpq connection string, e.g. "dbname=geo host=localhost"
        table_a: fully-qualified name of the first geometry table
        table_b: fully-qualified name of the second geometry table
        geom_col: name of the geometry column (must exist on both tables)
        srid: expected SRID; used in the ST_Transform guard comment only
    """
    # Use parameterised identifiers via psycopg2.sql to prevent injection
    from psycopg2 import sql

    query = sql.SQL(
        """
        EXPLAIN (ANALYZE, COSTS, TIMING, BUFFERS, FORMAT JSON)
        SELECT a.id, b.id
        FROM {table_a} a
        JOIN {table_b} b
          ON ST_Intersects(a.{geom_col}, b.{geom_col})
        WHERE ST_SRID(a.{geom_col}) = {srid}   -- guard: fails fast on SRID mismatch
        """
    ).format(
        table_a=sql.Identifier(*table_a.split(".")),
        table_b=sql.Identifier(*table_b.split(".")),
        geom_col=sql.Identifier(geom_col),
        srid=sql.Literal(srid),
    )

    with psycopg2.connect(conn_string) as conn:
        # Autocommit avoids an implicit transaction wrapping EXPLAIN ANALYZE,
        # which would skew timing on the first run due to lock acquisition.
        conn.autocommit = True
        with conn.cursor() as cur:
            cur.execute(query)
            plan_json: list[dict] = cur.fetchone()[0]

    def _extract(node: dict, depth: int = 0) -> list[dict]:
        """Recursively flatten the JSON plan tree into a list of annotated nodes."""
        actual_rows = node.get("Actual Rows", 0)
        plan_rows = node.get("Plan Rows", 1)
        ratio = round(actual_rows / max(plan_rows, 1), 2)

        entry = {
            "depth": depth,
            "node_type": node.get("Node Type"),
            "actual_time_ms": round(node.get("Actual Total Time", 0), 3),
            "actual_rows": actual_rows,
            "plan_rows": plan_rows,
            "row_estimate_ratio": ratio,
            # Spatial-specific fields — None if absent
            "index_cond": node.get("Index Cond"),
            "filter": node.get("Filter"),
            "recheck_cond": node.get("Recheck Cond"),
            # Buffer stats reveal cache efficiency
            "shared_hit_blocks": node.get("Shared Hit Blocks", 0),
            "shared_read_blocks": node.get("Shared Read Blocks", 0),
            "loops": node.get("Actual Loops", 1),
        }

        rows = [entry]
        for child in node.get("Plans", []):
            rows.extend(_extract(child, depth + 1))
        return rows

    return _extract(plan_json[0]["Plan"])


def flag_spatial_problems(nodes: list[dict]) -> list[str]:
    """
    Return a list of human-readable warnings based on the extracted plan nodes.
    Covers the four most common spatial join pathologies.
    """
    warnings = []
    for node in nodes:
        ntype = node["node_type"]

        # 1. Sequential scan instead of index scan
        if ntype == "Seq Scan" and node["actual_rows"] > 10_000:
            warnings.append(
                f"Seq Scan on large node ({node['actual_rows']} rows) — "
                "check GiST index existence and run ANALYZE."
            )

        # 2. Lossy bitmap causing Recheck overhead
        if node["recheck_cond"]:
            warnings.append(
                "Recheck Cond present — work_mem is forcing lossy bitmap storage. "
                "Run: SET work_mem = '64MB'; before this query."
            )

        # 3. Row-estimate divergence above 3x
        if node["row_estimate_ratio"] > 3 or node["row_estimate_ratio"] < 0.33:
            warnings.append(
                f"{ntype}: row estimate ratio {node['row_estimate_ratio']} "
                "(expected 0.33–3.0) — run ANALYZE on both tables."
            )

        # 4. High read blocks relative to hit blocks (cold cache / I/O bound)
        total_blocks = node["shared_hit_blocks"] + node["shared_read_blocks"]
        if total_blocks > 0:
            read_ratio = node["shared_read_blocks"] / total_blocks
            if read_ratio > 0.5:
                warnings.append(
                    f"{ntype}: {read_ratio:.0%} of block reads are from disk — "
                    "consider CLUSTER or increasing shared_buffers."
                )

    return warnings


# --- Usage example ---
# conn = "dbname=geo_prod host=db.internal user=app password=secret"
# nodes = analyze_spatial_join(conn, "public.parcels", "public.zones", srid=4326)
# for warning in flag_spatial_problems(nodes):
#     print(warning)
#
# # Print all index-level nodes for manual review:
# for n in nodes:
#     if n["node_type"] in ("Index Scan", "Bitmap Index Scan"):
#         print(
#             f"  [{n['node_type']}] {n['actual_time_ms']}ms | "
#             f"rows {n['actual_rows']}/{n['plan_rows']} "
#             f"(ratio {n['row_estimate_ratio']}) | "
#             f"index_cond={n['index_cond']}"
#         )

The ST_SRID guard in the query body is worth noting: it causes the query to fail immediately if the geometry column has an unexpected SRID rather than silently returning wrong results from an implicit coordinate system mismatch. For explicit reprojection patterns, see bounding box filtering and how it handles mixed-SRID datasets.

Configuration and Tuning Knobs

These four GUC settings directly affect what EXPLAIN ANALYZE shows for spatial joins:

GUC Recommended value Effect on plan
work_mem 32MB–128MB per session Eliminates Recheck Cond; allows lossless bitmap storage
enable_seqscan off (diagnostic only) Forces index path; reveals true index cost for comparison
join_collapse_limit 1 (diagnostic) Fixes join order; expose suboptimal planner choices
random_page_cost 1.1 on SSD, 4.0 on HDD Controls planner preference for index vs sequential scans

Set work_mem per session, not globally, unless all connections run spatial workloads:

sql
-- Session-scoped: safe for high-concurrency environments
SET work_mem = '64MB';

-- Then run your join — Recheck Cond should disappear from EXPLAIN output
EXPLAIN (ANALYZE, BUFFERS, FORMAT JSON)
SELECT a.id, b.id
FROM parcels a
JOIN zones b ON ST_Intersects(a.geom, b.geom);

For partial GiST indexes on active subsets, random_page_cost tuning is particularly impactful because smaller indexes fit in shared_buffers and have a fundamentally different I/O profile than full-table indexes.

Verification Steps

After applying any tuning, confirm the improvements with this checklist query. It returns a single diagnostic row you can assert against in CI or a monitoring script:

python
import psycopg2
import json

def verify_spatial_join_uses_index(
    conn_string: str,
    table_a: str,
    table_b: str,
    geom_col: str = "geom",
) -> dict[str, bool]:
    """
    Returns a dict of pass/fail assertions for the four critical plan properties.
    Raises AssertionError if any assertion fails (use in test suites or CI).
    """
    from psycopg2 import sql

    query = sql.SQL(
        "EXPLAIN (ANALYZE, FORMAT JSON) "
        "SELECT a.id, b.id FROM {a} a JOIN {b} b "
        "ON ST_Intersects(a.{g}, b.{g}) LIMIT 1000"
    ).format(
        a=sql.Identifier(*table_a.split(".")),
        b=sql.Identifier(*table_b.split(".")),
        g=sql.Identifier(geom_col),
    )

    with psycopg2.connect(conn_string) as conn:
        conn.autocommit = True
        with conn.cursor() as cur:
            cur.execute(query)
            plan = cur.fetchone()[0]

    plan_text = json.dumps(plan)

    results = {
        "uses_index_scan": (
            "Index Scan" in plan_text or "Bitmap Index Scan" in plan_text
        ),
        "has_index_cond_bbox": "&&" in plan_text and "Index Cond" in plan_text,
        "no_recheck_cond": "Recheck Cond" not in plan_text,
        "no_seq_scan": "Seq Scan" not in plan_text,
    }

    failed = [k for k, v in results.items() if not v]
    if failed:
        raise AssertionError(
            f"Spatial join plan assertions failed: {failed}\n"
            f"Full plan: {json.dumps(plan, indent=2)}"
        )

    return results

Run this in staging after every schema change or bulk data load. A passing result confirms the GiST index is active, bounding-box pruning is happening in the index layer, and there is no lossy bitmap overhead.

Gotchas Checklist

  • Index Cond can appear without the GiST index being useful. If the index has low selectivity — for example, a geometry column where all bounding boxes overlap heavily — the planner may use the index but still scan nearly every row through Filter. Check actual_rows at the Index Scan node versus actual_rows at its parent to see how much Filter reduces the set.

  • loops > 1 multiplies all costs. Actual Total Time in JSON format is per-loop, not cumulative. Multiply actual_time_ms by loops to get the true contribution of nested-loop inner nodes. This catches join order problems where a small table drives a large indexed table with many outer iterations.

  • SRID mismatches produce silent wrong results, not errors, in some predicates. ST_Intersects on mismatched SRIDs may return rows but compare geometries in incompatible coordinate spaces. Always verify SRID consistency with SELECT DISTINCT ST_SRID(geom) FROM table_name before trusting join results.

  • Statistics staleness after bulk loads skews Plan Rows. PostGIS stores bounding-box histograms in pg_statistic. After a bulk insert or ST_Transform migration, run ANALYZE table_name immediately. The planner’s row estimates for spatial joins degrade faster than for scalar columns because geometry distributions are inherently less uniform. See composite spatial indexes for how compound statistics affect multi-column spatial predicates.

  • enable_indexscan = off for diagnostics only, never in production. When you toggle SET enable_indexscan = off to force a sequential scan baseline, do it inside an explicit transaction you roll back, or in a separate psycopg2 connection with autocommit = False. Accidentally leaving it off in an application connection pool silently degrades every subsequent spatial query on that connection.

FAQ

Why is Index Cond missing from my spatial join EXPLAIN output?

The planner has fallen back to a sequential scan. The most common causes are: no GiST index on the geometry column, stale statistics making the planner underestimate selectivity, or a table small enough that the planner prefers a full scan. Run ANALYZE on both tables, confirm the index exists with \d table_name, and check enable_indexscan is on.

What does Recheck Cond mean in a spatial join plan?

Recheck Cond appears when work_mem forces PostgreSQL to store the bitmap in lossy mode, recording which pages rather than which individual tuples matched. The heap scan must then re-test each tuple’s bounding box. Increase work_mem per session (SET work_mem = '64MB') to eliminate this overhead on large ST_DWithin radius searches and ST_Intersects joins.

How large a row-estimate divergence is acceptable?

A Plan Rows to Actual Rows ratio up to 3x is generally acceptable. Beyond that, the planner is likely choosing suboptimal join orders or overestimating index selectivity. Run ANALYZE, and if geometry is clustered or skewed, consider CLUSTER table_name USING idx_geom to physically co-locate nearby geometries on disk, which also improves buffer hit rates visible in the BUFFERS output.