Tuning ST_Intersects joins with bounding-box prefilters

This page takes a spatial join that has degraded into a sequential scan and rewrites it, step by step, into a plan where the GiST index does the filtering and the exact geometry test runs on a handful of candidates.

When to use this approach

  • EXPLAIN shows a sequential scan on a geometry table you know is indexed. Something in the predicate is hiding the column from the planner — a transform, a buffer, or a function wrapper. Reach for this page rather than adding hardware.
  • The plan already uses the index, but the filter step dominates the run time. Here the index is doing its job and the exact recheck is the cost; the fix is subdivision, described in step 4, not a different predicate.
  • You need the join to stay correct while getting faster. Everything below preserves the result set exactly. If you are willing to change the answer — snapping points to a grid, for instance — the cheaper approach in discrete global grid macros may suit better.

The general rules behind these steps are collected in spatial joins and predicate tuning. Which step you need depends on what the plan is telling you, so read it first and enter the sequence at the right point.

Decision tree: which tuning step a query plan calls for Starting from the question of whether the plan shows an index condition, a no branch leads to removing per-row functions from the join condition, covered in step two. A yes branch asks whether rows removed by filter is large: yes leads to subdividing polygons in step four, no leads to shrinking the candidate set with a non-spatial predicate. Does the plan show Index Cond? Step 2 · unwrap the column move transforms into staging Is Rows Removed by Filter in the millions? Step 4 · subdivide the polygons the exact test is the bottleneck Shrink the candidate set instead date window, tenant filter, tighter extent no yes yes no

Prerequisites

  • PostGIS 3.1 or newer, with postgis installed in the target schema’s search path.
  • A GiST index on the geometry column of the joined table — see using spatial index hints in dbt materializations.
  • Both inputs stored in the same SRID; a mixed-SRID join is a different bug with the same symptom.
  • Permission to run EXPLAIN (ANALYZE, BUFFERS) against the target, and enough of a data volume that the planner does not choose a sequential scan simply because the table is tiny.

Step-by-step instructions

1. Capture the plan you are starting from

Never tune without a baseline. Compile the model and run its SQL under EXPLAIN.

bash
dbt compile --select int_pings_zoned
psql -X -f target/compiled/dbt_geospatial/models/intermediate/int_pings_zoned.sql \
     -c "EXPLAIN (ANALYZE, BUFFERS, TIMING)"
text
Nested Loop  (cost=0.00..8421934.11 rows=118422 width=104) (actual time=0.512..214883.9 rows=3822104 loops=1)
  Join Filter: st_intersects(p.geom, st_transform(z.geom, 4326))
  Rows Removed by Join Filter: 41155882
  ->  Seq Scan on stg_trip_pings p  ...
  ->  Seq Scan on stg_service_zones z  ...

Two lines matter: Rows Removed by Join Filter in the tens of millions means every pair was formed and tested, and the st_transform inside the join filter is why — the index on z.geom cannot answer a question about st_transform(z.geom, 4326).

2. Move per-row functions out of the join condition

Any function applied to a joined geometry column inside ON disables the index for that column. Push it upstream into staging, where it is computed once and can itself be indexed.

sql
-- models/staging/stg_service_zones.sql
{{ config(
    materialized = 'table',
    post_hook = "CREATE INDEX IF NOT EXISTS {{ this.name }}_geom_idx ON {{ this }} USING GIST (geom)"
) }}

select
    zone_id,
    zone_name,
    zone_priority,
    st_transform(geom, {{ var('canonical_srid') }}) as geom
from {{ source('ops', 'service_zones') }}

Verify the column is now stored in the canonical SRID, so no transform is needed at join time:

sql
select distinct st_srid(geom) from {{ ref('stg_service_zones') }};
-- Expect exactly one row: 4326
Moving ST_Transform out of the join condition restores index eligibility Before, on the left: the join condition wraps the zone geometry in ST_Transform, so the GiST index on that column cannot be used and the planner falls back to a sequential scan. After, on the right: the transform has moved into the staging model, the join condition references the stored column directly, and the index scan is available. before · transform in the join ON ST_Intersects(p.geom, ST_Transform(z.geom, 4326)) GiST on z.geom unused Seq Scan 41M pairs tested after · transform in staging ON ST_Intersects(p.geom, z.geom) GiST on z.geom used Index Scan candidates only The rewrite changes no rows — only which side of the pipeline pays for the projection.

3. Make the box prefilter explicit when the planner needs convincing

ST_Intersects already expands to && plus an exact test, so an explicit && is usually redundant. It earns its place in two situations: when the join condition also contains non-spatial clauses that confuse selectivity estimation, and when you are joining through a view or CTE that has blocked the rewrite. Writing it out costs nothing and makes the intent legible.

sql
select
    p.ping_id,
    z.zone_id
from {{ ref('stg_trip_pings') }} as p
join {{ ref('stg_service_zones') }} as z
  on p.geom && z.geom                       -- box overlap: index-served
 and st_intersects(p.geom, z.geom)          -- exact test: candidates only
where p.observed_at >= '{{ var("window_start") }}'

Verify that the plan now shows an index scan with a recheck condition:

text
Nested Loop  (actual time=0.084..2841.7 rows=3822104 loops=1)
  ->  Seq Scan on stg_trip_pings p
  ->  Index Scan using stg_service_zones_geom_idx on stg_service_zones z
        Index Cond: (p.geom && geom)
        Filter: st_intersects(p.geom, geom)

Index Cond is the line to look for. If it is absent, the index is still not in play and the cause is upstream of this step.

4. Subdivide the fat polygons

With the index working, the remaining cost is the exact test. A polygon with 90,000 vertices is expensive to test against and has a bounding box that catches candidates it will ultimately reject. ST_Subdivide replaces it with many small pieces.

sql
-- models/staging/stg_service_zones.sql (final form)
select
    zone_id,
    zone_name,
    zone_priority,
    st_subdivide(st_transform(geom, {{ var('canonical_srid') }}), 256) as geom
from {{ source('ops', 'service_zones') }}
where st_isvalid(geom)

Because subdivision emits several rows per original zone, the join must collapse them again:

sql
select distinct
    p.ping_id,
    z.zone_id
from {{ ref('stg_trip_pings') }} as p
join {{ ref('stg_service_zones') }} as z
  on st_intersects(p.geom, z.geom)

Verify the vertex count fell and the row count rose in the expected proportion:

sql
select
    count(*) as pieces,
    max(st_npoints(geom)) as max_vertices,
    count(distinct zone_id) as zones
from {{ ref('stg_service_zones') }};
-- Expect max_vertices <= 256 and pieces >> zones
Why subdivision tightens the bounding box of a complex polygon On the left, one irregular coastal polygon with a single large bounding box that covers a wide area of empty sea, so many points fall inside the box but outside the polygon. On the right, the same polygon after ST_Subdivide, split into nine pieces each with a small bounding box that hugs the shape, so far fewer points become candidates. one polygon · one loose box 4 candidates inside the box, 0 inside the shape nine pieces · nine tight boxes 1 candidate inside a box, 0 inside the shape

5. Re-analyze and re-measure

Statistics from before the rewrite will mislead the planner about the new tables.

bash
psql -c "ANALYZE analytics.stg_service_zones; ANALYZE analytics.int_pings_zoned;"
dbt build --select int_pings_zoned+

Verify by re-running the EXPLAIN from step 1 and comparing actual time on the top node. A join that started at 214 seconds and lands under 5 is a typical outcome for this sequence; if the number barely moves, the bottleneck is elsewhere — check for a missing index on the other side, or a distinct over a very wide row.

Configuration reference

Setting Where Accepted values Spatial note
max_vertices ST_Subdivide(geom, n) 8 – 10000, commonly 128 – 512 Lower means tighter boxes and more rows; below ~64 the row multiplication starts to dominate
fillfactor GiST index storage parameter 10 – 100, default 90 Lower it for tables that receive updates so the index has room to grow in place
max_parallel_workers_per_gather session or role 0 – 8 Above 0 lets the exact recheck spread across cores; irrelevant while the plan is index-bound
work_mem session or role e.g. 64MB Raises the ceiling before a hash join spills; a spatial sort spilling to disk is a common hidden cost
random_page_cost database 1.0 – 4.0 On SSD-backed storage, 1.1 makes the planner more willing to choose index scans
var('canonical_srid') dbt_project.yml any valid SRID Keeps the staging transform and the join agreeing on one value

Gotchas & edge cases

  • && is not a substitute for the exact test. It compares boxes only, so it returns rows whose geometries do not touch. Always pair it with ST_Intersects.
  • Subdivision changes row counts, and therefore tests. A unique test on zone_id in staging will start failing the moment you subdivide. Move that test to a model that has re-aggregated the pieces.
  • ST_Subdivide rejects invalid geometry. Run it after validation, never before; the pattern is in quarantining invalid geometries in staging.
  • A CTE can block the rewrite. If a materialized CTE hides the geometry column behind an expression, the planner sees an opaque column with no statistics. Reference the model directly in the join.
  • Indexes are not free on write. Each GiST index slows inserts into an incremental model. Create the index in a post-hook that runs after the load, not before it.

FAQ

Do I still need && if I already call ST_Intersects?

Usually not — PostGIS rewrites ST_Intersects(a, b) into a && b AND _ST_Intersects(a, b) itself. Write it explicitly when the plan shows the index is being skipped despite an index-eligible predicate, or when reviewers benefit from seeing the two stages spelled out. It never changes the result.

What value should I pass to ST_Subdivide?

Start at 256 vertices. The tuning curve is flat across a wide middle range: much lower and you pay for row multiplication and the downstream distinct; much higher and the boxes stay loose enough that the candidate set does not shrink. Measure the exact-test time with EXPLAIN (ANALYZE) at 128, 256 and 512 on your own data before settling.

Why did my row count change after subdividing?

Because the polygon table now holds one row per piece. The join emits a row per matching piece, so a point inside a zone that was split into four pieces near its location can match more than once. Collapse with select distinct on the natural key, or aggregate; the avoiding Cartesian blowups guide covers the grain tests that catch this in CI.

The plan uses the index but the join is still slow — what now?

Look at where the time is spent rather than which nodes appear. If Rows Removed by Filter is large, the exact test is the cost and subdivision is the answer. If the index scan itself dominates, the index may be bloated — rebuild it — or the candidate set genuinely is large, in which case reduce it with a non-spatial predicate such as a date window before the geometry work.

Up: Part of Spatial Joins & Predicate Tuning.