🔧 Fix PostGIS & databases Performance & scale Spatial SQL & PostGIS · step 14 of 15 Updated

My PostGIS Spatial Query Is Slow: How to Fix It

Problem statement

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

SELECT p.id, w.ward_name
FROM parcels p
JOIN wards w ON ST_Intersects(p.geom, w.geom)
WHERE p.class = 'residential';

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 two seconds.

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);
Triage rows matching each query-plan symptom to its cause and fix.
Six causes. Five of them announce themselves in the plan.
What the plan shows Cause Fix
Seq Scan on a large table no index, or the predicate cannot use one index it; rewrite the predicate
no Index Cond: (geom && …) the geometry column is wrapped in a function keep the column bare
Rows Removed by Filter ≫ rows kept bounding boxes are poor approximations ST_Subdivide the big geometries
estimated rows ≫ or ≪ actual stale statistics ANALYZE
high read in BUFFERS data not cached, or the table is bloated VACUUM, more RAM, fewer columns
plan looks fine, wall clock is long the transfer, not the query reduce what you select

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

Step-by-step solution

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);
Aggregate  (actual time=412803.1..412803.1 rows=1 loops=1)
  Buffers: shared hit=1204 read=784112
  ->  Nested Loop  (actual time=0.9..411204.8 rows=41288 loops=1)
        ->  Seq Scan on flood_zones f  (actual rows=1204 loops=1)
        ->  Seq Scan on parcels p  (actual rows=34 loops=1204)
              Filter: st_intersects(p.geom, f.geom)
              Rows Removed by Filter: 4012850

Four numbers matter here:

  • Seq Scan on parcels — no index is being used.
  • loops=1204 — that scan runs once per flood zone.
  • Rows Removed by Filter: 4012850 — nearly the whole table, discarded, 1,204 times over. That is 4.8 billion geometry comparisons.
  • read=784112 — 784,112 buffer pages read from disk rather than found in cache: about 6 GB of I/O.

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';
      indexname       |                        indexdef
----------------------+---------------------------------------------------------
 parcels_pkey         | CREATE UNIQUE INDEX parcels_pkey ON 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%';
 relname | indexrelname     | idx_scan
---------+------------------+----------
 parcels | parcels_geom_idx |        0

idx_scan = 0 after running the query means the planner rejected it. 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.7, 53.3, -1.9, 53.9, 4326)

-- ✅ transform the constant instead; geom stays bare
WHERE geom && ST_Transform(ST_MakeEnvelope(-2.7, 53.3, -1.9, 53.9, 4326), 27700)
-- ❌ 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');
 relname | n_live_tup | n_dead_tup | last_analyze | last_autoanalyze
---------+------------+------------+--------------+------------------
 parcels |    4012884 |     882014 |              | 2026-06-02 03:14
 wards   |        215 |          0 |              | 2026-06-02 03:14

Statistics from June on a table that has since been reloaded. And 882,014 dead tuples — rows deleted or updated but not yet reclaimed, which the query still has to read past.

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

Scene comparing a compact polygon whose bounding box fits tightly with a sprawling multipolygon whose box covers everything.
The index tests boxes. When the box is nothing like the geometry, the index selects almost everything.

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  (actual rows=204 loops=1)
  Index Cond: (geom && f.geom)
  Filter: st_intersects(geom, f.geom)
  Rows Removed by Filter: 88412        ← 400 false positives per real hit

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;
 id  | vertices | box_ratio
-----+----------+-----------
 412 |   402113 |     184.2
 118 |   288044 |      92.7

A box_ratio of 184 means the bounding box is 184 times the area of the geometry — a river network or a coastline, where the box covers a whole region. 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. This routinely produces order-of-magnitude improvements and is the single 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, class, geom FROM parcels WHERE ward_code = 'E05011368'"

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.84 s
full read     18.42 s
transfer      17.58 s  (95%)

The query is fine. Ninety-five per cent 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.

Code examples

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("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")
    for depth, node in nodes:
        print(f"  {'  ' * depth}{node['Node Type']:<22} "
              f"rows={node.get('Actual Rows', 0):>9,} loops={node.get('Loops', 1)}")
    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 412,803 ms
  Aggregate              rows=        1 loops=1
    Nested Loop          rows=   41,288 loops=1
      Seq Scan           rows=    1,204 loops=1
      Seq Scan           rows=       34 loops=1204

  → Seq Scan on parcels (34 rows × 1204 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")

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          88,412 s   204.18 s
ST_DWithin 500             88,412 s     1.82 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                    6284 MB
  ✗ flood_zones                 412 MB  NO GIST INDEX
  ✗ incidents                    88 MB  34.2% dead
  ✓ wards                         4 MB

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.

Explanation

Panels contrasting predicates that wrap the indexed column with equivalent ones that transform the constant instead.
One rule removes most unused-index reports: transform the constant, never the column.

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: 4 million rows scanned per outer row means billions of comparisons where 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 88,000 candidates and the exact test keeps 204, the index has done nothing useful 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 4.9 GB of geometry to WKB, moving it over a socket, and constructing 4 million Shapely objects. When a plan looks healthy and the Python call takes eighteen seconds, 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.

Edge cases or notes

  • 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 by 4 million rows.
  • 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.

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 takes twenty seconds.

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.