Point-in-polygon joins at scale in dbt

This page builds a point-in-polygon assignment model that stays predictable as the point table grows into the hundreds of millions: one row per point, one zone per point, an incremental strategy that only touches new points, and tests that fail the build if either invariant breaks.

When to use this approach

• Your point table grows continuously and your polygon table barely changes. That asymmetry is what makes incremental assignment worth the extra configuration; if both sides churn, a full table build is simpler and often faster.
• Downstream consumers assume one zone per point. Dashboards that count trips per zone break silently under fanout, so the grain policy below matters more than the raw speed.
• A full rebuild no longer fits the build window. Before reaching for incremental complexity, confirm the join itself is tuned — an untuned join is slow at every scale, and the fix is in tuning ST_Intersects joins with bounding-box prefilters.

The decision is really about which curve your build sits on. A full rebuild grows with the whole history; an incremental assignment grows with the daily slice, and pays a fixed overhead for the merge.

Prerequisites

• Points and polygons in one SRID, both validated in staging.
• A GiST index on the polygon geometry, created in a post-hook.
• A monotonic column on the point table — an ingest timestamp or an append-only id — for the incremental filter.
• dbt-utils installed, for the row-count and expression tests used at the end.
• A stable natural key on each polygon; the tie-break depends on it, and a key that changes between builds will reshuffle historical assignments.

Step-by-step instructions

1. Prepare the polygon side once

Assignment quality is decided here. Validate, transform to the canonical SRID, subdivide, and index — in that order.

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

with valid as (
    select
        zone_id,
        zone_name,
        zone_priority,
        st_transform(geom, {{ var('canonical_srid') }}) as geom
    from {{ source('ops', 'zones') }}
    where st_isvalid(geom)
)

select
    zone_id,
    zone_name,
    zone_priority,
    st_subdivide(geom, 256) as geom
from valid

Verify the pieces cover the same area as the originals — subdivision must not drop slivers:

select
    abs(sum(st_area(geom)) - (select sum(st_area(geom)) from {{ source('ops', 'zones') }}))
      / nullif((select sum(st_area(geom)) from {{ source('ops', 'zones') }}), 0) as area_drift
from {{ ref('stg_zones') }};
-- Expect a value near zero (floating point noise only)

2. Write the assignment model with an explicit grain

The join emits one row per point per matching piece. Collapse it to one row per point with a total ordering, so the answer is identical on every run.

-- models/intermediate/int_point_zone_assignment.sql
{{ config(
    materialized = 'incremental',
    unique_key = 'ping_id',
    incremental_strategy = 'delete+insert',
    post_hook = "ANALYZE {{ this }}"
) }}

with new_points as (
    select ping_id, trip_id, observed_at, geom
    from {{ ref('stg_trip_pings') }}
    {% if is_incremental() %}
    where observed_at > (select coalesce(max(observed_at), '1900-01-01') from {{ this }})
    {% endif %}
),

matched as (
    select distinct on (p.ping_id)
        p.ping_id,
        p.trip_id,
        p.observed_at,
        z.zone_id,
        z.zone_name
    from new_points as p
    join {{ ref('stg_zones') }} as z
      on st_intersects(p.geom, z.geom)
    order by p.ping_id, z.zone_priority asc, z.zone_id asc
)

select * from matched

distinct on with order by ping_id, zone_priority, zone_id is the whole grain policy: first row per point wins, ties broken by priority then by identifier, both deterministic. Dropping either tie-break column reintroduces run-to-run drift.

Verify the grain immediately after the first build:

select count(*) as duplicate_pings
from (select ping_id from {{ ref('int_point_zone_assignment') }} group by ping_id having count(*) > 1) d;
-- Expect 0

3. Handle points that match nothing

A left join is not the answer — it would reintroduce fanout. Instead, assign unmatched points a sentinel zone in a second pass so the model still holds one row per point, and so the count of unmatched points is measurable rather than invisible.

matched as ( ... as above ... ),

unmatched as (
    select
        p.ping_id, p.trip_id, p.observed_at,
        null::text as zone_id,
        'outside_service_area' as zone_name
    from new_points as p
    where not exists (
        select 1 from {{ ref('stg_zones') }} as z
        where st_intersects(p.geom, z.geom)
    )
)

select * from matched
union all
select * from unmatched

Verify that the two branches partition the input exactly:

select
    (select count(*) from {{ ref('stg_trip_pings') }}) as input_rows,
    (select count(*) from {{ ref('int_point_zone_assignment') }}) as output_rows;
-- Expect the two numbers to be equal on a full refresh

4. Add the tests that hold the grain

# models/intermediate/schema.yml
models:
  - name: int_point_zone_assignment
    description: One row per ping, with its assigned zone or an explicit outside-area marker.
    tests:
      - dbt_utils.equal_rowcount:
          compare_model: ref('stg_trip_pings')
    columns:
      - name: ping_id
        tests: [unique, not_null]
      - name: zone_name
        tests:
          - not_null
          - accepted_values:
              values: ['outside_service_area']
              quote: true
              config:
                where: "zone_id is null"

The equal_rowcount test is the important one: it is the assertion that the join did not change the grain, and it runs in seconds regardless of geometry size.

5. Refresh when boundaries move

Incremental assignment freezes history by design. When a zone boundary changes, every historical point in the affected area still carries its old assignment — correct if you want assignment-as-of-ingest, wrong if you want current boundaries applied retroactively. Make the choice explicit and automate it.

# Boundary release day: rebuild the assignment from scratch
dbt build --select stg_zones+ --full-refresh

Verify the rebuild changed what you expected and nothing else:

select zone_name, count(*) as pings
from {{ ref('int_point_zone_assignment') }}
group by zone_name
order by pings desc;
-- Compare against the previous release's counts; movement should be confined to changed zones

Configuration reference

Gotchas & edge cases

• A point on a shared boundary intersects both neighbours. That is correct behaviour, not a bug; the tie-break decides. Never rely on ST_Within to resolve it, because a boundary point is within neither polygon under some implementations.
• max(observed_at) in the incremental filter drops late arrivals. If pings can arrive out of order, use a lookback window (> max(observed_at) - interval '3 days') with delete+insert so re-processed rows replace rather than duplicate.
• The unmatched branch must use the same column list and order. A union all with mismatched types compiles and then fails at run time on the warehouse, often only in production where the sentinel branch first becomes non-empty.
• A distinct on over a wide row is expensive. Select only the columns the mart needs before deduplicating; carrying geometry through the sort is a common and avoidable cost.
• Zone attributes copied into the assignment freeze with it. If zone_name changes, historical rows keep the old label. Carry only zone_id and join names at the mart layer when labels are volatile.

FAQ

Related

• Spatial Joins & Predicate Tuning — predicate rules and grain policies in general.
• Avoiding Cartesian Blowups in Spatial Joins — the tests that catch fanout before it reaches a mart.
• Incremental Spatial Materializations — incremental strategies for geometry-heavy models.

Up: Part of Spatial Joins & Predicate Tuning.