Spatial Joins & Predicate Tuning

A spatial join is where most dbt geospatial projects first meet real pain. The model compiles, the SQL is legible, the results look plausible in a preview — and then the job that took eleven seconds against a sample takes forty minutes against production, or worse, it returns 4.1 million rows where the mart is supposed to hold 3.8 million. Both symptoms come from the same place: a join predicate the planner cannot serve from an index, joined against geometry whose bounding boxes overlap far more than the geometries themselves do. This topic covers how to write spatial joins that stay index-eligible, how to reshape the geometry so the index actually helps, and how to prove the join did not quietly change the grain of your model.

Spatial joins sit in the intermediate layer of the architecture described in core fundamentals and architecture for dbt geospatial — after staging has validated and normalized geometry, before marts serve it. That position matters for tuning: by the time a join runs, every input should already carry one canonical SRID and a valid geometry, so the join can be optimized purely as a query problem rather than as a data-cleanliness problem. If your join is slow and your inputs are dirty, fix the inputs first in geometry validation and data quality; a ST_MakeValid call inside a join condition defeats every index on the table.

Prerequisites checklist

• PostGIS 3.1+ or DuckDB spatial 0.10+. Both ship the bounding-box operators this topic depends on; PostGIS additionally exposes ST_Subdivide, which does the heavy lifting for oversized polygons.
• A GiST index on every geometry column that appears on the right-hand side of a join, created in a dbt post-hook. Index strategy is covered in index hints for spatial queries.
• A single canonical SRID across both join inputs, enforced per enforcing a canonical SRID across dbt models. A mixed-SRID join either errors or silently disables the index.
• ANALYZE after every full refresh. The planner’s choice between a nested loop with an index scan and a hash join comes from row-count and selectivity statistics that do not exist until the table is analyzed.
• Permission to create indexes in the target schema, since the fastest join is usually the one preceded by an index that dbt built as part of the model.

Architecture context: where the join belongs in the DAG

A spatial join should be a named intermediate model, not an inline subquery inside a mart. Materializing it gives you somewhere to hang an index, a row-count test, and an incremental strategy — and it means the expensive work runs once per build rather than once per downstream consumer.

The arrows are ref() edges, so this is also the lineage graph — the same structure discussed in spatial model dependency graphs. The important property is that geometry work stops at the intermediate layer. A mart that joins geometry again, to add one attribute, pays the whole cost a second time.

Which predicates the index can serve

PostGIS spatial indexes store bounding boxes, not geometries. Every index-accelerated predicate therefore works in two stages: an index scan that finds candidate rows whose boxes overlap, then an exact recheck that runs the real predicate on those candidates. A predicate is index-eligible when PostGIS can derive a box search from it. This is the single most useful table to internalize.

Two of these deserve emphasis because they are the mistakes that appear most often in review. ST_Distance(a, b) < 500 and ST_DWithin(a, b, 500) return identical rows and have wildly different plans: the first computes an exact distance for every pair in the cross product, the second asks the index for the boxes within 500 units and computes exact distances only for those. The second is the one to write, every time. Likewise, buffering inside a join condition looks harmless but converts an index scan into a sequential scan plus a per-row buffer construction — if you need a buffered geometry repeatedly, materialize it as a column in staging and index that column instead.

Configuration walkthrough

Give the join model its own configuration block. Two settings do most of the work: a materialization that persists the result, and a post-hook that indexes and analyzes it.

-- models/intermediate/int_pings_zoned.sql
{{ config(
    materialized = 'table',
    post_hook = [
      "CREATE INDEX IF NOT EXISTS {{ this.name }}_geom_idx ON {{ this }} USING GIST (ping_geom)",
      "ANALYZE {{ this }}"
    ]
) }}

select
    p.ping_id,
    p.trip_id,
    p.observed_at,
    p.geom as ping_geom,
    z.zone_id,
    z.zone_name
from {{ ref('stg_trip_pings') }} as p
join {{ ref('stg_service_zones') }} as z
  on st_intersects(p.geom, z.geom)

The ANALYZE in the post-hook is not decoration. Without fresh statistics the planner assumes a default selectivity for the spatial predicate and will frequently choose a hash join over the index-driven nested loop, which is exactly the wrong plan when one side is small and selective.

For the polygon side, the single highest-leverage preparation step is subdivision. A country boundary or a utility service area can be a polygon with tens of thousands of vertices whose bounding box covers half the map; every point in that box becomes a candidate, and every candidate pays an expensive exact test against a huge ring. ST_Subdivide cuts it into many small pieces whose boxes are tight.

-- 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,
    st_subdivide(geom, 256) as geom   -- max 256 vertices per piece
from {{ source('ops', 'service_zones') }}
where st_isvalid(geom)

Subdivision multiplies the row count of the polygon table, so the join now returns one row per point per polygon piece. Deduplicate downstream with distinct on the natural key, or aggregate with group by. That trade — more rows in staging, a distinct in the join model, and an order-of-magnitude faster exact test — is almost always worth taking for polygons above a few thousand vertices.

Core implementation: controlling fanout

A spatial join is a many-to-many relation. Where zones overlap, or where a point falls on a shared boundary, one ping matches several zones and the model’s grain silently changes from “one row per ping” to “one row per ping per zone”. Downstream aggregates then double-count without any error being raised.

Decide the grain explicitly. The three usual policies:

-- Policy A: keep every match, and say so in the model name and its unique key.
select p.ping_id, z.zone_id, ...
from pings p join zones z on st_intersects(p.geom, z.geom)

-- Policy B: one row per ping, choose the zone with the largest overlap.
select distinct on (p.ping_id)
       p.ping_id, z.zone_id,
       st_area(st_intersection(p.geom, z.geom)) as overlap_area
from pings p join zones z on st_intersects(p.geom, z.geom)
order by p.ping_id, overlap_area desc, z.zone_id

-- Policy C: one row per ping, deterministic tie-break by priority then id.
select distinct on (p.ping_id)
       p.ping_id, z.zone_id
from pings p join zones z on st_intersects(p.geom, z.geom)
order by p.ping_id, z.zone_priority asc, z.zone_id asc

Policy B is the honest choice for overlapping administrative areas but only makes sense for areal geometry — ST_Area of a point intersection is always zero, so for point inputs use Policy C. Whichever you pick, the tie-break must be deterministic: distinct on without a total ordering produces a different answer on every run, which turns into an unexplainable diff in the mart and a failing snapshot test.

For point-in-polygon specifically, prefer ST_Intersects over ST_Within even when containment is what you mean. They use the index identically, but boundary semantics differ: a point exactly on a shared edge intersects both neighbouring polygons and is within neither in some implementations. Pair ST_Intersects with an explicit tie-break and you get a defined answer instead of a silent drop.

Validation and testing

Two dbt tests catch the overwhelming majority of spatial-join regressions, and neither requires a spatial extension to run.

# models/intermediate/schema.yml
models:
  - name: int_pings_zoned
    description: One row per ping per matched service zone.
    tests:
      - dbt_utils.equal_rowcount:
          compare_model: ref('stg_trip_pings')
          config:
            severity: warn   # only equal when the join is guaranteed 1:1
    columns:
      - name: ping_id
        tests:
          - not_null
      - name: zone_id
        tests:
          - not_null
          - relationships:
              to: ref('stg_service_zones')
              field: zone_id

For a join you intend to be one-to-one, assert the grain directly with a uniqueness test on the key; for a join you intend to fan out, assert an upper bound instead, so an accidental cross join fails the build:

-- tests/assert_zone_fanout_bounded.sql
select ping_id, count(*) as zone_matches
from {{ ref('int_pings_zoned') }}
group by ping_id
having count(*) > 8

A third check worth adding is a coverage assertion — the count of input rows that matched nothing. Silent non-matching is how a CRS mistake hides: reproject one side by a few metres and the join still runs, returns fewer rows, and nobody notices until a dashboard drops.

-- tests/assert_ping_zone_coverage.sql
select count(*) as unmatched
from {{ ref('stg_trip_pings') }} p
left join {{ ref('int_pings_zoned') }} z using (ping_id)
where z.ping_id is null
having count(*) > 1000

Advanced patterns

Prefilter with the box operator when you must use a non-indexable predicate. If a business rule genuinely needs ST_Distance, keep the index in play by adding an ST_DWithin guard first; the planner uses the index for the guard and evaluates the exact distance only on survivors.

where st_dwithin(a.geom, b.geom, 1000)      -- index-eligible, cuts the candidate set
  and st_distance(a.geom, b.geom) < 1000    -- exact, runs on candidates only

Join on a grid key before joining on geometry. Attaching an H3 or geohash cell to both sides turns the first pass into an equality join the planner can hash, with geometry used only to resolve cells that straddle a boundary. That pattern is developed in discrete global grid macros and pairs well with the partitioning approach in partitioning geospatial tables with H3 and dbt.

Make the join incremental on the point side. Points accumulate; zones change rarely. An incremental model that processes only new pings against the full zone table turns a nightly full join into a minutes-long append, following the pattern in incremental spatial materializations. Remember to trigger a full refresh when zone boundaries change, or historical rows keep their old assignment.

Watch the parallel plan. PostGIS marks most spatial functions as parallel-safe, but a nested loop over a GiST index will not be parallelized unless max_parallel_workers_per_gather is above zero and the outer relation is large enough to justify it. When a join is CPU-bound on the exact recheck rather than I/O-bound, raising the worker count is often a bigger win than any rewrite.

Troubleshooting

FAQ

Related

• Tuning ST_Intersects Joins with Bounding-Box Prefilters — the prefilter pattern in full, with plans before and after.
• Point-in-Polygon Joins at Scale in dbt — subdivision, grain policy and incremental assignment for the most common join of all.
• Avoiding Cartesian Blowups in Spatial Joins — how fanout happens and the tests that catch it.
• Index Hints for Spatial Queries — building and verifying the indexes these joins depend on.
• Optimizing Proximity Joins — nearest-neighbour work, where the tuning rules differ.

Up: Part of Core Fundamentals & Architecture for dbt Geospatial.