Benchmarking Shapely vs PostGIS Predicate Checks

You need to know which engine evaluates a single spatial predicate faster — ST_Intersects, ST_DWithin, or ST_IsValid against their Shapely 2.x vectorized equivalents — and you want a measurement you can defend, not a number contaminated by cache warmth or a missing index. This page isolates each predicate, times it correctly with timeit on the Python side and EXPLAIN ANALYZE on the database side, and shows how warm-versus-cold cache and index presence change the result. It is the micro-benchmark counterpart to Comparing Spatial Validation Engines, zooming from whole-pipeline choice down to a single predicate.

Prerequisites

• Shapely 2.0+ with intersects, dwithin, is_valid, and STRtree. These vectorized predicates operate on NumPy arrays through one GEOS call; the pre-2.0 object-at-a-time API is not a fair opponent for PostGIS and must not be used.
• PostGIS 3.4+ on PostgreSQL 15+ with the ability to run EXPLAIN (ANALYZE, BUFFERS) and to toggle indexes. You need permission to DROP INDEX and CREATE INDEX to measure the index contribution.
• Matching geometries in one projected coordinate reference system (CRS). ST_DWithin thresholds are in CRS units; a benchmark that runs Shapely in degrees and PostGIS in metres compares two different questions. Fix both on the same projected CRS.
• A quiet host. Micro-benchmarks are sensitive to CPU contention and cache eviction; close other database sessions and background Python processes before measuring.
• Gotcha: the operating system page cache sits beneath the PostgreSQL buffer cache. To force a genuinely cold run you must restart PostgreSQL and drop OS caches, not merely open a new session — otherwise your “cold” number is still partly warm.

Step-by-Step Procedure

Step 1 — Isolate one predicate at a time

import numpy as np
import shapely
from shapely.geometry import Point

rng = np.random.default_rng(7)
n = 200_000
pts = shapely.points(rng.uniform(0, 10_000, n), rng.uniform(0, 10_000, n))
polys = shapely.buffer(
    shapely.points(rng.uniform(0, 10_000, 5_000), rng.uniform(0, 10_000, 5_000)),
    25.0,
)  # 5k circular zones to test predicates against

Verification: confirm pts and polys are NumPy arrays of Shapely geometries (pts.dtype == object, shapely.is_valid(pts).all()). Isolating a single predicate per benchmark is what makes the eventual index and cache effects attributable rather than tangled.

Step 2 — Time the Shapely predicates with timeit

import timeit

# ST_IsValid equivalent: pure per-feature, no index possible
t_isvalid = min(timeit.repeat(lambda: shapely.is_valid(pts), repeat=5, number=20)) / 20

# ST_Intersects equivalent WITH an STRtree (index built once, outside the timer)
tree = shapely.STRtree(polys)
def intersects_indexed():
    cand = tree.query(pts, predicate="intersects")   # bbox candidates + exact refine
    return cand
t_intersects = min(timeit.repeat(intersects_indexed, repeat=5, number=20)) / 20

# ST_DWithin equivalent: dwithin short-circuits at the threshold
def dwithin_check():
    return shapely.dwithin(pts, polys[0], 50.0)
t_dwithin = min(timeit.repeat(dwithin_check, repeat=5, number=20)) / 20

print(f"is_valid   {t_isvalid*1e3:8.2f} ms/call")
print(f"intersects {t_intersects*1e3:8.2f} ms/call")
print(f"dwithin    {t_dwithin*1e3:8.2f} ms/call")

Verification: build the STRtree outside the timed function — if tree = shapely.STRtree(polys) sits inside intersects_indexed, you are timing index construction on every call and the comparison against a persistent GiST index is meaningless. Take the min across repeats for the per-call figure.

Step 3 — Time the PostGIS predicates with EXPLAIN ANALYZE

-- ST_IsValid: per-row, no index involvement
EXPLAIN (ANALYZE, BUFFERS) SELECT count(*) FROM pts WHERE NOT ST_IsValid(geom);

-- ST_Intersects: index-assisted spatial join
EXPLAIN (ANALYZE, BUFFERS)
SELECT count(*)
FROM pts p JOIN zones z ON ST_Intersects(p.geom, z.geom);

-- ST_DWithin: index-assisted proximity, short-circuits at 50 units
EXPLAIN (ANALYZE, BUFFERS)
SELECT count(*)
FROM pts p JOIN zones z ON ST_DWithin(p.geom, z.geom, 50.0);

Verification: in each plan confirm the intended access path — a Seq Scan where you expected an Index Scan on the GiST index means the predicate is running unindexed and the number is not what you think. Record Execution Time and the shared hit/shared read split to label the cache state.

Step 4 — Toggle the index to quantify its contribution

-- Cold run: drop the index, force a full scan
DROP INDEX IF EXISTS zones_gix;
EXPLAIN (ANALYZE, BUFFERS)
SELECT count(*) FROM pts p JOIN zones z ON ST_Intersects(p.geom, z.geom);

-- Indexed run: rebuild and re-measure
CREATE INDEX zones_gix ON zones USING GIST (geom);
ANALYZE zones;
EXPLAIN (ANALYZE, BUFFERS)
SELECT count(*) FROM pts p JOIN zones z ON ST_Intersects(p.geom, z.geom);

Verification: the indexed run should show a nested-loop or index scan with dramatically fewer shared buffers touched than the sequential-scan run. The ratio between the two execution times is the index contribution — report it explicitly so nobody mistakes an index win for an engine win.

Interpreting Results

Read the three predicates as three different stories. ST_IsValid and shapely.is_valid are both purely per-feature — no index applies to either — so this is the cleanest engine-versus-engine comparison, and Shapely’s single vectorized GEOS pass over an in-memory array is typically very competitive because there is no query planning or result serialisation. ST_Intersects is where the index dominates: an indexed PostGIS join or an STRtree-backed Shapely query both leave an unindexed elementwise comparison far behind, so the honest comparison is indexed-against-indexed. ST_DWithin adds short-circuiting: because it stops as soon as a feature is within the threshold, it usually beats computing an exact distance then filtering, in both engines.

The most important reading is which lever moved the number. If enabling the GiST index cut execution time by fifty-fold, then the headline is “use an index,” not “PostGIS beats Shapely.” If a warm cache halved the time, the headline is a cache effect, not throughput. Micro-benchmarks earn their keep precisely by attributing the change to the right lever, which whole-pipeline timings in PostGIS vs GeoPandas for Validation at Scale cannot isolate because there many effects move at once.

Finally, remember that the STRtree query returns bounding-box candidates and must be refined by the exact predicate to equal what ST_Intersects computes. If your Shapely intersects benchmark looks impossibly fast, check that you passed predicate="intersects" to refine candidates rather than counting raw bounding-box overlaps. This candidate-versus-exact distinction is the single most common way a predicate benchmark lies.

Gotchas & Edge Cases

Counting bounding-box candidates instead of exact hits. tree.query(pts) without a predicate argument returns every geometry whose bounding box overlaps — a superset of true intersections. Benchmarking that against ST_Intersects, which returns exact hits, compares two different computations. Always pass the predicate to refine.

“Cold” that is only session-cold. Opening a new psql session does not clear the buffer cache. A true cold run needs a PostgreSQL restart and an OS cache drop; otherwise your cold number is warm and the cache effect vanishes into noise.

timeit including array construction. If geometry arrays or the STRtree are built inside the timed callable, you measure setup, not the predicate. Construct all inputs before the timed block and reference them as closures.

ST_DWithin argument order and units. The threshold is in the geometry’s CRS units, and swapping to a geographic CRS silently changes metres to degrees — a threshold of 50 becomes 50 degrees, matching almost everything. Confirm the CRS before trusting a proximity benchmark.

JIT and query plan caching in PostgreSQL. For large joins PostgreSQL may spend time on just-in-time compilation on the first execution, inflating the cold number beyond pure I/O. Check the JIT section of the plan; if present, decide whether to disable it with SET jit = off for a cleaner comparison.

When to Escalate

Once a single predicate is characterised, the decision usually returns to a larger frame. If the predicate you benchmarked is a per-feature validity check, the production pattern for running it at volume in Python is Implementing Shapely Geometry Checks in Python, which wraps these predicates in a repair-aware pipeline. If the predicate is a set-level join whose scaling you now need to project across feature counts, move up to PostGIS vs GeoPandas for Validation at Scale for the throughput and memory curves. When the micro-benchmark shows the two engines genuinely close on your predicate, let developer familiarity and pipeline fit decide, using the weighted matrix in Comparing Spatial Validation Engines.

Frequently Asked Questions

Related

• Implementing Shapely Geometry Checks in Python — the production pattern for running these Shapely predicates with repair
• PostGIS vs GeoPandas for Validation at Scale — full-pipeline throughput and memory once a predicate is characterised
• Comparing Spatial Validation Engines — the weighted decision matrix these micro-benchmarks feed

Back to Comparing Spatial Validation Engines