My PostGIS Spatial Query Is Slow: How to Fix It

Problem statement

The query has run for three minutes and you have no idea whether it is nearly done or nearly hopeless:

SELECT COUNT(*)
FROM parcels p
JOIN flood_zones f ON ST_Intersects(p.geom, f.geom);

You already built the index. You already ran ANALYZE. It is still slow. Meanwhile a colleague runs an apparently harder query on the same tables and it comes back in a second.

Spatial slowness has a small number of causes, and they are distinguishable in about a minute using EXPLAIN. Guessing — adding indexes, raising work_mem, rewriting the join — usually changes nothing, because the actual cause is almost always one of six specific things, and five of them are visible in the query plan.

Quick answer

Read the plan before changing anything:

EXPLAIN (ANALYZE, BUFFERS)
SELECT p.id, w.ward_name
FROM parcels p JOIN wards w ON ST_Intersects(p.geom, w.geom);

The three single-line fixes that solve most cases:

ANALYZE parcels;                                        -- stale statistics
-- ST_Distance(a, b) < 500   →   ST_DWithin(a, b, 500)  -- index-aware predicate
CREATE INDEX ON parcels USING GIST (geom);              -- the index itself

Read the plan, check the index, keep the column bare

1. Read the plan, not the query

EXPLAIN (ANALYZE, BUFFERS, VERBOSE)
SELECT COUNT(*) FROM parcels p
JOIN flood_zones f ON ST_Intersects(p.geom, f.geom);

Finalize Aggregate  (actual time=197814 rows=1 loops=1)
  Buffers: shared hit=428690213 read=3776
  ->  Gather  (actual rows=3 loops=1)
        Workers Launched: 2
        ->  Partial Aggregate  (actual rows=1 loops=3)
              ->  Nested Loop  (actual rows=9960 loops=3)
                    Join Filter: st_intersects(p.geom, f.geom)
                    Rows Removed by Join Filter: 48375285
                    ->  Parallel Seq Scan on parcels p  (actual rows=40355 loops=3)
                    ->  Seq Scan on flood_zones f  (actual rows=1199 loops=121064)

That is the plan, 198 seconds long, for HM Land Registry's 121,064 Manchester parcels against 1,199 Environment Agency flood zones, before any spatial index (Example 1 reads the same plan). Four things matter here:

• Seq Scan on both tables — no index is being used.
• loops=121064 — the flood-zone scan runs once for every parcel.
• Rows Removed by Join Filter: 48375285 — per worker; across the three processes about 145 million geometry comparisons, to find 29,880 intersecting pairs.
• shared hit=428690213 — 429 million buffer pages visited, all but 3,776 of them from cache. The time went on geometry tests, not on disk.

actual time is what happened; cost is what the planner guessed. When those disagree wildly, the statistics are the problem.

2. Check that the index exists — and is used

SELECT indexname, indexdef FROM pg_indexes WHERE tablename = 'parcels';

ix_public_parcels_id  CREATE INDEX ix_public_parcels_id ON public.parcels USING btree (id)

No GiST index. Create it, and analyse:

CREATE INDEX parcels_geom_idx ON parcels USING GIST (geom);
ANALYZE parcels;

Then check it is being used rather than merely existing:

SELECT relname, indexrelname, idx_scan
FROM pg_stat_user_indexes WHERE indexrelname LIKE '%geom%';

parcels  parcels_geom_idx  1199

Here the flood-zone join used it 1,199 times — once per flood zone, the outer side of the loop — and the same count query took 1.1 s (Example 4) instead of 197.8 s (Example 1). idx_scan = 0 after running a query would mean the planner rejected the index. That is a different problem from not having one, and steps 3 and 4 cover its two causes.

3. Keep the indexed column bare

An index on geom can only be used when a query mentions geom directly. Wrap it in a function and the index is unreachable:

-- ❌ every row is transformed before comparison; no index
WHERE ST_Transform(geom, 4326) && ST_MakeEnvelope(-2.26, 53.46, -2.22, 53.48, 4326)

-- ✅ transform the constant instead; geom stays bare
WHERE geom && ST_Transform(ST_MakeEnvelope(-2.26, 53.46, -2.22, 53.48, 4326), 27700)

On the Manchester parcels the first took 0.52 s and the second 0.012 s (Example 4). They also returned different counts, 2,875 against 2,901: a rectangle in degrees is not a rectangle in British National Grid, so transforming the box changes which bounding boxes it touches. Check that a rewrite returns the same rows before trusting the speed-up.

-- ❌ a function of the column
WHERE ST_Distance(geom, %(pt)s) < 500

-- ✅ the same rows, index-aware
WHERE ST_DWithin(geom, %(pt)s, 500)

-- ❌ buffers every row
WHERE ST_Intersects(ST_Buffer(geom, 100), %(area)s)

-- ✅ buffer the constant
WHERE ST_Intersects(geom, ST_Buffer(%(area)s, 100))

ST_DWithin is the substitution worth memorising. It expands the reference geometry's bounding box by the distance, asks the index for candidates, and measures exactly on those — the same two-phase design described in PostGIS spatial indexes explained. ST_Distance(...) < d computes a distance for every row in the table before comparing.

When you genuinely need a function of the column in the predicate, index the function:

CREATE INDEX parcels_centroid_idx ON parcels USING GIST (ST_Centroid(geom));

4. Refresh the statistics

The planner chooses between an index scan and a sequential scan by estimating how many rows match. Bad estimates produce bad choices in both directions.

SELECT relname, n_live_tup, n_dead_tup, last_analyze, last_autoanalyze
FROM pg_stat_user_tables WHERE relname IN ('parcels', 'wards');

parcels  121064  24213

After an update that rewrote a fifth of the parcels, with autovacuum held off, 24,213 dead tuples — old row versions not yet reclaimed, which a sequential scan still has to read past (Example 4).

VACUUM ANALYZE parcels;

A plan whose estimated rows differ from actual rows by more than about an order of magnitude is almost always this.

5. Check whether the bounding boxes are any good

The index filters on bounding boxes and an exact test runs on the survivors. When the boxes are poor approximations, the exact test does nearly all the work — and the plan says so:

Index Scan using parcels_geom_idx on parcels p  (actual rows=25 loops=1199)
  Index Cond: (geom && f.geom)
  Filter: st_intersects(geom, f.geom)
  Rows Removed by Filter: 110        ← about four false positives per real hit

That is the indexed join on the Manchester data: for each flood zone the index returned about 135 candidate parcels and the exact test kept 25. Four false positives per hit is healthy; hundreds would mean the boxes are useless.

Find the offenders:

SELECT id,
       ST_NPoints(geom)                                     AS vertices,
       ROUND((ST_Area(ST_Envelope(geom)) / NULLIF(ST_Area(geom), 0))::numeric, 1)
                                                            AS box_ratio
FROM flood_zones
ORDER BY ST_NPoints(geom) DESC
LIMIT 5;

1086    15  10627.9
 519    44   7612.2
1042    10   6572.6
 937    15   6347.9
 558    73   2609.1

A box_ratio of 10,628 means the bounding box is ten thousand times the area of the geometry. In this layer the worst offenders are not large rivers but thin diagonal slivers with 10 to 73 vertices; a coastline or a river network does the same at far larger scale. Chop them up:

CREATE TABLE flood_zones_sub AS
SELECT id, zone_class, ST_Subdivide(geom, 256) AS geom
FROM flood_zones;

CREATE INDEX ON flood_zones_sub USING GIST (geom);
ANALYZE flood_zones_sub;

ST_Subdivide splits each geometry into pieces of at most 256 vertices, each with a small, tight box. Queries then join against flood_zones_sub and aggregate by id. On these flood zones, averaging 892 vertices with mostly compact boxes, it cut the join from 1.10 s to 0.82 s (Example 4), and found one more intersecting pair, 29,881 against 29,880, because the pieces' shared edges are computed afresh; on layers with huge, sprawling geometries the gain is far larger, which is why it is the most under-used PostGIS performance technique.

6. Make sure it is the query that is slow

import time
import geopandas as gpd
from sqlalchemy import create_engine, text

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

SQL = "SELECT id, geom FROM parcels"

with engine.begin() as con:
    t0 = time.perf_counter()
    con.execute(text(f"EXPLAIN (ANALYZE) {SQL}")).fetchall()
    server = time.perf_counter() - t0

t0 = time.perf_counter()
gdf = gpd.read_postgis(SQL, engine, geom_col="geom")
client = time.perf_counter() - t0

print(f"server-side  {server:6.2f} s")
print(f"full read    {client:6.2f} s")
print(f"transfer     {client - server:6.2f} s  ({100*(client-server)/client:.0f}%)")

server-side    0.02 s
full read      1.50 s
transfer       1.48 s  (99%)

The query is fine. Nearly all of the time is serialising geometry, moving it over a socket, and building Shapely objects. No index will help — the fix is to select less, as covered in how to read a large PostGIS table.

Automatic plan diagnosis, before and after

The data used on this page

Real tables in a throwaway PostGIS container: HM Land Registry's INSPIRE Index Polygons for Manchester as parcels (each given its ONS ward code), the ONS ward boundaries for Manchester (December 2024) as wards, the Environment Agency's Flood Zones 2 and 3 over the city as flood_zones, and railway and tram stops from OpenStreetMap, as saved on 10 October 2026, as stations. The setup loads them unindexed but analysed — the state the page diagnoses — and sets PGPORT and PGPASSWORD so the connection strings below work unchanged.

import io
import os
import time
import zipfile

import geopandas as gpd
import requests
from sqlalchemy import create_engine, text

# OpenStreetMap as Overpass returned it on 10 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-postgis-query-is-slow-2026-10-10.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

os.environ.update(PGPORT="55433", PGPASSWORD="guide")          # the throwaway container below
INSPIRE = "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Manchester_City_Council.zip"
WARDS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
         "Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query"
         "?where=LAD24NM%3D%27Manchester%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson")
FLOOD = ("https://environment.data.gov.uk/geoservices/datasets/04532375-a198-476e-985e-0579a0a11b47/"
         "ogc/features/v1/collections/Flood_Zones_2_3_Rivers_and_Sea/items")

if os.system("docker inspect guide-postgis >/dev/null 2>&1") != 0:
    os.system("docker run -d --name guide-postgis -e POSTGRES_USER=user -e POSTGRES_PASSWORD=guide "
              "-e POSTGRES_DB=gis -p 127.0.0.1:55433:5432 postgis/postgis:16-3.4 >/dev/null")
    time.sleep(10)

if not os.path.exists("Manchester_City_Council.zip"):
    with requests.Session() as session:                        # the service sets a cookie, then redirects
        open("Manchester_City_Council.zip", "wb").write(session.get(INSPIRE, timeout=600).content)

engine = create_engine("postgresql+psycopg://user@localhost/gis")
with engine.begin() as con:
    loaded = con.execute(text("select to_regclass('parcels') is not null")).scalar()
if not loaded:
    parcels = gpd.read_file("/vsizip/Manchester_City_Council.zip/Land_Registry_Cadastral_Parcels.gml")
    parcels = parcels[["INSPIREID", "geometry"]].rename(columns={"INSPIREID": "inspire_id"}).set_crs(27700, allow_override=True)
    wards = gpd.read_file(WARDS).rename(columns={"WD24CD": "ward_code", "WD24NM": "ward_name"})
    centres = gpd.GeoDataFrame(parcels[["inspire_id"]], geometry=parcels.representative_point(), crs=27700)
    parcels["ward_code"] = gpd.sjoin(centres, wards[["ward_code", "geometry"]], how="left").ward_code.values
    box = wards.to_crs(4326).total_bounds
    params, features, url = {"bbox": ",".join(map(str, box)), "limit": 100, "f": "json"}, [], FLOOD
    while url:
        for attempt in range(5):                                  # the service sometimes answers with an error
            page = requests.get(url, params=params, timeout=120).json()
            if "features" in page:
                break
            time.sleep(10)
        features += page["features"]
        url, params = next((l["href"] for l in page["links"] if l["rel"] == "next"), None), None
    flood = gpd.GeoDataFrame.from_features(features, crs=4326).to_crs(27700)[["flood_zone", "geometry"]]
    import osmnx as ox
    ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"   # overpass-api.de mirror
    ox.settings.overpass_rate_limit = False
    stops = ox.features_from_polygon(wards.to_crs(4326).union_all(),
                                     {"railway": ["station", "halt", "tram_stop"]}).reset_index()
    stations = stops[stops.geom_type == "Point"][["name", "railway", "geometry"]].to_crs(27700)
    for name, frame in (("parcels", parcels), ("wards", wards), ("flood_zones", flood), ("stations", stations)):
        frame.rename_geometry("geom").to_postgis(name, engine, index=True, index_label="id")
    with engine.begin() as con:
        for name in ("parcels", "wards", "flood_zones", "stations"):   # start from the unindexed state the page diagnoses
            con.execute(text(f"DROP INDEX IF EXISTS idx_{name}_geom"))
            con.execute(text(f"ANALYZE {name}"))      # with statistics, as autovacuum would leave them
with engine.begin() as con:
    for name in ("parcels", "wards", "flood_zones", "stations"):
        print(name, con.execute(text(f"SELECT count(*) FROM {name}")).scalar())

parcels 121064
wards 32
flood_zones 1199
stations 110

Example 1: an automatic plan diagnosis

import json
from sqlalchemy import create_engine, text

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

def diagnose(sql, engine, params=None):
    with engine.begin() as con:
        raw = con.execute(
            text(f"EXPLAIN (ANALYZE, BUFFERS, FORMAT JSON) {sql}"), params or {}
        ).scalar()
    plan = (raw[0] if isinstance(raw, list) else json.loads(raw)[0])["Plan"]

    findings, nodes = [], []
    def walk(node, depth=0):
        nodes.append((depth, node))
        kind = node["Node Type"]
        rows = node.get("Actual Rows", 0)
        loops = node.get("Actual Loops", 1)
        removed = node.get("Rows Removed by Filter", 0)
        est = node.get("Plan Rows", 0)

        if kind == "Seq Scan" and rows * loops > 100_000:
            findings.append(
                f"Seq Scan on {node.get('Relation Name')} "
                f"({rows:,} rows × {loops} loops) — index it, or check the predicate")
        if removed and rows and removed / max(rows, 1) > 10:
            findings.append(
                f"{removed:,} rows removed by the exact filter against {rows:,} kept "
                f"— bounding boxes are poor; try ST_Subdivide")
        if est and rows and (est / max(rows, 1) > 10 or rows / max(est, 1) > 10):
            findings.append(
                f"{kind}: estimated {est:,} rows, got {rows:,} — run ANALYZE")
        for child in node.get("Plans", []):
            walk(child, depth + 1)
    walk(plan)

    hit = plan.get("Shared Hit Blocks", 0)
    read = plan.get("Shared Read Blocks", 0)
    if read and read > hit:
        findings.append(
            f"{read:,} blocks read from disk vs {hit:,} cached "
            f"(~{read * 8 / 1024:,.0f} MB) — cold cache, bloat, or too many columns")

    print(f"total time {plan['Actual Total Time']:,.0f} ms, buffers hit {hit:,} read {read:,}")
    for depth, node in nodes:
        removed = node.get("Rows Removed by Filter", 0) + node.get("Rows Removed by Join Filter", 0)
        print(f"  {'  ' * depth}{node['Node Type']:<22} "
              f"rows={node.get('Actual Rows', 0):>9,} loops={node.get('Actual Loops', 1)}"
              + (f" removed={removed:,}" if removed else ""))
    print()
    for f in findings or ["no obvious plan problem — check the transfer size"]:
        print(f"  → {f}")
    return findings

diagnose("""
    SELECT COUNT(*) FROM parcels p
    JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
""", engine)

total time 197,814 ms, buffers hit 428,690,213 read 3,776
  Aggregate              rows=        1 loops=1
    Gather                 rows=        3 loops=1
      Aggregate              rows=        1 loops=3
        Nested Loop            rows=    9,960 loops=3 removed=48,375,285
          Seq Scan               rows=   40,355 loops=3
          Seq Scan               rows=    1,199 loops=121064

  → Seq Scan on parcels (40,355 rows × 3 loops) — index it, or check the predicate
  → Seq Scan on flood_zones (1,199 rows × 121064 loops) — index it, or check the predicate

The thresholds are heuristics, not laws — a 10:1 estimate error is worth flagging, and a 3:1 one usually is not. What the function really buys is turning a wall of plan text into two or three sentences you can act on, which makes reading plans a habit rather than an ordeal.

Example 2: before and after, measured honestly

import time
from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://user@localhost/gis")
with engine.begin() as con:                       # the index first: step 2
    for table in ("parcels", "flood_zones", "wards", "stations"):
        con.execute(text(f"CREATE INDEX IF NOT EXISTS {table}_geom_idx ON {table} USING GIST (geom)"))
        con.execute(text(f"ANALYZE {table}"))

VARIANTS = {
    "ST_Distance < 500": """
        SELECT COUNT(*) FROM parcels p, stations s
        WHERE ST_Distance(p.geom, s.geom) < 500
    """,
    "ST_DWithin 500": """
        SELECT COUNT(*) FROM parcels p, stations s
        WHERE ST_DWithin(p.geom, s.geom, 500)
    """,
}

def bench(sql, engine, runs=3):
    times = []
    for _ in range(runs):
        with engine.begin() as con:
            t0 = time.perf_counter()
            n = con.execute(text(sql)).scalar()
            times.append(time.perf_counter() - t0)
    return n, min(times)

results = {}
for name, sql in VARIANTS.items():
    n, t = bench(sql, engine)
    results[name] = (n, t)
    print(f"{name:<22} {n:>10,} rows  {t:>8.2f} s")

counts = {n for n, _ in results.values()}
assert len(counts) == 1, f"variants disagree: {results}"

ST_Distance < 500          83,999 rows      1.26 s
ST_DWithin 500             83,999 rows      0.05 s

The assertion is the part that makes this a valid benchmark rather than a demonstration. A rewrite that is faster and returns different rows is not an optimisation, and ST_DWithin versus ST_Distance < is exactly the kind of substitution where an off-by-one on the boundary would be easy to miss. Taking min of several runs reduces the effect of an unlucky cache state.

Example 3: a query-performance health check

from sqlalchemy import create_engine, text
import pandas as pd

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

HEALTH = """
SELECT
    c.relname                                        AS table_name,
    pg_size_pretty(pg_total_relation_size(c.oid))    AS total_size,
    s.n_live_tup                                     AS live_rows,
    s.n_dead_tup                                     AS dead_rows,
    ROUND(100.0 * s.n_dead_tup / NULLIF(s.n_live_tup + s.n_dead_tup, 0), 1)
                                                     AS pct_dead,
    GREATEST(s.last_analyze, s.last_autoanalyze)     AS analysed,
    EXISTS (SELECT 1 FROM pg_index i
            JOIN pg_class ic ON ic.oid = i.indexrelid
            JOIN pg_am am ON am.oid = ic.relam
            WHERE i.indrelid = c.oid AND am.amname = 'gist') AS has_gist
FROM pg_class c
JOIN pg_namespace n ON n.oid = c.relnamespace
JOIN pg_stat_user_tables s ON s.relid = c.oid
WHERE n.nspname = 'public'
  AND EXISTS (SELECT 1 FROM geometry_columns g
              WHERE g.f_table_name = c.relname)
ORDER BY pg_total_relation_size(c.oid) DESC
"""

df = pd.read_sql(HEALTH, engine)
for r in df.itertuples():
    flags = []
    if not r.has_gist:                       flags.append("NO GIST INDEX")
    if r.pct_dead and r.pct_dead > 20:       flags.append(f"{r.pct_dead}% dead")
    if r.analysed is None:                   flags.append("never analysed")
    mark = "✗" if flags else "✓"
    print(f"  {mark} {r.table_name:<24} {r.total_size:>10}  {'; '.join(flags)}")

  ✓ parcels                       55 MB  
  ✓ flood_zones                   18 MB  
  ✓ wards                         64 kB  
  ✓ stations                      48 kB  

Run this before diagnosing an individual query. Two of the six causes — a missing index and table bloat — are properties of a table rather than of a query, so finding them here saves reading a plan at all. geometry_columns is PostGIS's own catalogue view, so the check only ever reports spatial tables.

Example 4 — every plan and table on this page, from the Manchester data

def rows(sql):
    with engine.begin() as con:
        return con.execute(text(sql)).fetchall()


def timed(sql):
    start = time.perf_counter()
    result = rows(sql)
    return result, time.perf_counter() - start


JOIN = "SELECT COUNT(*) FROM parcels p JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)"

with engine.begin() as con:
    con.execute(text("DROP INDEX IF EXISTS parcels_geom_idx"))
print("indexes before:", [r[0] for r in rows("SELECT indexname FROM pg_indexes WHERE tablename = 'parcels'")])
with engine.begin() as con:
    con.execute(text("CREATE INDEX parcels_geom_idx ON parcels USING GIST (geom)"))
    con.execute(text("ANALYZE parcels"))
result, secs = timed(JOIN)
n = result[0][0]
print(f"indexed join: {n:,} intersecting pairs in {secs:.2f} s")
print("index use:", rows("SELECT relname, indexrelname, idx_scan FROM pg_stat_user_indexes "
                         "WHERE indexrelname = 'parcels_geom_idx'"))
with engine.begin() as con:
    plan = [r[0] for r in con.execute(text(f"EXPLAIN (ANALYZE) {JOIN}"))]
print("\n".join(l for l in plan if "Index" in l or "Removed" in l))

ENV = "ST_MakeEnvelope(-2.26, 53.46, -2.22, 53.48, 4326)"
for label, sql in (("ST_Transform(geom) && box", f"SELECT count(*) FROM parcels WHERE ST_Transform(geom, 4326) && {ENV}"),
                   ("geom && ST_Transform(box)", f"SELECT count(*) FROM parcels WHERE geom && ST_Transform({ENV}, 27700)")):
    result, secs = timed(sql)
    print(f"{label:28} {result[0][0]:,} rows in {secs:.3f} s")

with engine.begin() as con:
    con.execute(text("ALTER TABLE parcels SET (autovacuum_enabled = false)"))
    con.execute(text("UPDATE parcels SET ward_code = ward_code WHERE id % 5 = 0"))
with engine.begin() as con:
    con.execute(text("SELECT pg_stat_force_next_flush()"))       # PostgreSQL 15+: publish the counters now
time.sleep(2)
with engine.begin() as con:
    con.execute(text("SELECT pg_stat_clear_snapshot()"))
    print("after updating a fifth:", con.execute(text(
        "SELECT relname, n_live_tup, n_dead_tup FROM pg_stat_user_tables WHERE relname = 'parcels'")).fetchall())
with engine.connect().execution_options(isolation_level="AUTOCOMMIT") as con:
    con.execute(text("VACUUM ANALYZE parcels"))
    con.execute(text("ALTER TABLE parcels RESET (autovacuum_enabled)"))

print("worst boxes:", rows("""SELECT id, ST_NPoints(geom),
    ROUND((ST_Area(ST_Envelope(geom)) / NULLIF(ST_Area(geom), 0))::numeric, 1) AS ratio
    FROM flood_zones ORDER BY ratio DESC NULLS LAST LIMIT 5"""))
count, mean = rows("SELECT count(*), round(avg(ST_NPoints(geom))) FROM flood_zones")[0]
print(f"flood zones: {count:,}, mean vertices {mean}")
with engine.begin() as con:
    con.execute(text("DROP TABLE IF EXISTS flood_zones_sub"))
    con.execute(text("CREATE TABLE flood_zones_sub AS "
                     "SELECT id, flood_zone, ST_Subdivide(geom, 256) AS geom FROM flood_zones"))
    con.execute(text("CREATE INDEX ON flood_zones_sub USING GIST (geom)"))
    con.execute(text("ANALYZE flood_zones_sub"))
SUB = ("SELECT COUNT(DISTINCT (p.id, f.id)) FROM parcels p "
       "JOIN flood_zones_sub f ON ST_Intersects(p.geom, f.geom)")
for label, sql in (("whole flood zones", JOIN), ("subdivided", SUB)):
    best = min(timed(sql)[1] for _ in range(3))
    print(f"{label:18} {rows(sql)[0][0]:,} pairs, best of 3 {best:.2f} s")

indexes before: ['ix_public_parcels_id']
indexed join: 29,880 intersecting pairs in 1.12 s
index use: [('parcels', 'parcels_geom_idx', 1199)]
        ->  Index Scan using parcels_geom_idx on parcels p  (cost=0.28..172.53 rows=12 width=367) (actual time=0.500..0.914 rows=25 loops=1199)
              Index Cond: (geom && f.geom)
              Rows Removed by Filter: 110
ST_Transform(geom) && box    2,875 rows in 0.510 s
geom && ST_Transform(box)    2,901 rows in 0.012 s
after updating a fifth: [('parcels', 121064, 24213)]
worst boxes: [(1086, 15, Decimal('10627.9')), (519, 44, Decimal('7612.2')), (1042, 10, Decimal('6572.6')), (937, 15, Decimal('6347.9')), (558, 73, Decimal('2609.1'))]
flood zones: 1,199, mean vertices 892
whole flood zones  29,880 pairs, best of 3 1.10 s
subdivided         29,881 pairs, best of 3 0.82 s

Too many comparisons or expensive ones

Spatial queries are slow for one of two reasons: the database is doing more geometry comparisons than it needs to, or the geometry comparisons it does are expensive. Almost every cause reduces to one of those, and the plan tells you which.

The first family is about the index. A GiST index turns a spatial join from a product into something close to linear by testing bounding boxes first — a cheap test that can over-report but never under-report, so it safely eliminates most candidates. When there is no index, or the planner cannot use the one that exists, the query degrades to a nested loop over the whole table. The characteristic sign is enormous Rows Removed by Filter on a Seq Scan, and the ratio is brutal: in step 1, 1,199 flood zones tested against each of 121,064 parcels meant about 145 million comparisons where tens of thousands would do.

Within that family, the subtle case is an index that exists and is not used. Postgres can only use an index on geom if the query's predicate mentions geom directly. Any function wrapping it — ST_Transform(geom, …), ST_Buffer(geom, …), ST_Distance(geom, …) — produces a value the index has never seen, so the index is unreachable and the function runs per row. The rule that avoids the whole class of problem is one sentence: keep the indexed column bare and transform the constant instead.

The second family is about geometry complexity. A bounding box is a good stand-in for a compact shape and a terrible one for a sprawling one. When the index selects hundreds of candidates for every one the exact test keeps, the index has done little and the expensive phase is doing all the work. This is not an indexing problem — no index setting improves it — but a modelling one, and ST_Subdivide fixes it by making the boxes match the geometries again.

Statistics sit underneath both. The planner's choice is a cost comparison, and a cost comparison with wrong row estimates picks wrongly. After a bulk load, n_live_tup and the column statistics can be badly stale, which is exactly when the table has changed most. ANALYZE is cheap and is the correct first response to any plan whose estimates and actuals disagree by an order of magnitude.

And then there is the case where the query is not the problem at all. EXPLAIN ANALYZE measures server-side execution. It does not measure serialising geometry to WKB, moving it over a socket, and constructing Shapely objects. Reading all 121,064 Manchester parcels, execution took hundredths of a second and the Python call over a second and a half; when a plan looks healthy and the call is still slow, no amount of tuning helps — the fix is to transfer less, which is a query-writing decision rather than a database one. That boundary is the subject of spatial SQL or GeoPandas, and knowing which side of it you are on is most of the diagnosis.

EXPLAIN against EXPLAIN ANALYZE

• EXPLAIN alone estimates; EXPLAIN ANALYZE executes. On a destructive statement, wrap it in a transaction and roll back.
• ANALYZE after every bulk load. Autovacuum's schedule is not tuned for a table that just changed wholesale.
• ST_DWithin uses the index; ST_Distance(...) < d does not. Same for ST_Intersects versus ST_Distance(...) = 0.
• ST_Subdivide needs an aggregation afterwards — one input geometry becomes many rows, so GROUP BY the original id.
• A geography column is slower than geometry because the calculations are spheroidal. Use geometry in a projected SRID when possible.
• work_mem affects sorts and hash joins, not spatial index scans. Raising it rarely helps a spatial query.
• Dead tuples above ~20% mean every scan reads past rows that no longer exist. VACUUM.
• Parallel query is disabled for many PostGIS functions unless they are marked PARALLEL SAFE; check pg_proc.proparallel if you expected parallelism.
• pg_stat_statements finds the slow queries in the first place, which is a different question from why one is slow.
• SET enable_seqscan = off is a diagnostic, not a fix — it proves the index would help, and should never be left on.

Internal links

• PostGIS spatial indexes explained — the mechanism behind five of these six causes
• How to do a spatial join in PostGIS with SQL — writing the join so the planner can use the index
• Spatial SQL or GeoPandas? — when the transfer, not the query, is the cost
• How to read a large PostGIS table into Python — fixing the transfer half
• How to write a GeoDataFrame to PostGIS — loading in a way that leaves the table healthy
• How to speed up GeoPandas: tips for large datasets — the client-side equivalent
• Spatial indexes explained: R-trees — the same two-phase idea in memory
• How to profile a slow Python GIS script — when it is neither the query nor the transfer

FAQ

I created the index and nothing got faster. Why?

Either you have not run ANALYZE, or the predicate wraps the geometry column in a function so the index cannot be used. Check for an Index Cond: (geom && …) line in the plan.

What does Rows Removed by Filter mean?

Rows that reached the exact geometry test and failed it. A small number is normal — that is the second phase doing its job. A number many times the rows kept means the bounding boxes are poor approximations.

Why is ST_Distance so much slower than ST_DWithin?

ST_Distance(geom, x) < d computes a distance for every row before comparing. ST_DWithin uses the index to find candidates first, then measures only those.

When should I use ST_Subdivide?

When a few geometries have very many vertices and bounding boxes much larger than the shapes — coastlines, river networks, national boundaries. Subdivide, index the pieces, and aggregate by the original id.

Does raising work_mem help?

Rarely for spatial queries. It affects sorts and hash joins, not GiST index scans. Read the plan before changing server settings.

My plan looks fine but the Python call is slow.

The time is in transfer and parsing, not the query. Select fewer columns, filter harder, or replace geometry with derived values such as ST_Area(geom).

How do I know whether the table is bloated?

n_dead_tup in pg_stat_user_tables. Above roughly 20% of live rows, every scan reads past rows that no longer exist. Run VACUUM ANALYZE.