Point-in-Polygon Joins at Scale

Stamping every customer with the delivery zone that contains them is a spatial join, and it belongs to the family of patterns covered under spatial joins. The specific job here: a customers table typed geometry(Point, 4326) holds several million rows, a delivery_zones table typed geometry(Polygon, 4326) holds a few thousand service areas, and each customer needs the id of the zone they sit inside. Done wrong, this is a Cartesian product that tests every point against every polygon. Done right, it is a nested loop where a GiST index turns each point into a handful of candidate polygons, streamed through Python in bounded batches.

Problem Statement

The assignment must populate customers.zone_id for the whole table and keep it current as new signups arrive. The naive query joins the two tables on ST_Contains and pulls the result into the application. With five million customers and four thousand zones, the theoretical comparison space is twenty billion point-in-polygon tests. Even at a microsecond each that is hours of pure CPU, and materialising the join in Python exhausts memory long before it finishes. The scalable version leans on the && bounding-box operator backed by a GiST index to reduce each point to only the polygons whose envelopes overlap it, resolves overlapping zones deterministically, and processes the table in id-ordered slices.

Why the Naive Approach Fails

The first attempt usually pulls the entire join client-side:

# Anti-pattern: no index guarantee, unbounded result set in Python memory
cur.execute("""
    SELECT c.id AS customer_id, z.id AS zone_id
    FROM customers c, delivery_zones z
    WHERE ST_Contains(z.geom, c.geom)
""")
assignments = cur.fetchall()   # millions of tuples land in the Python heap at once
for customer_id, zone_id in assignments:
    ...

There are two independent problems. The comma join expresses a cross product, and while ST_Contains does inject a bounding-box prefilter, a stale or missing GiST index on delivery_zones.geom leaves the planner to evaluate the exact ST_Contains predicate against every polygon for every point — an O(n·m) sequential nested loop. On a cold cache this can run for hours with no visible progress.

The second problem is fetchall(). Even when the join is index-assisted and fast on the server, dragging every result tuple into Python at once inflates the process heap and, in a container with a fixed memory limit, triggers an OOM kill. And if two delivery zones overlap — a common situation where an express zone sits inside a standard zone — the join emits two rows for the customer in the overlap, so a downstream UPDATE ... SET zone_id either errors or silently keeps whichever row arrives last.

Production-Ready Implementation

The correct query is a LATERAL join: for each customer, look up the single best containing zone through the index, then batch the whole table by primary key. The LATERAL subquery collapses overlaps to one deterministic winner via a priority column.

-- Run once during schema setup; the polygon index is what makes the join scale.
CREATE INDEX IF NOT EXISTS idx_delivery_zones_geom_gist
    ON delivery_zones USING GIST (geom);
ANALYZE delivery_zones;

import psycopg
from psycopg.rows import dict_row

DSN = "postgresql://etl:secret@10.0.3.4:5432/retail"
BATCH = 5_000   # customers scanned per keyset slice


def assign_zones() -> int:
    """
    Assign every customer the id of the single delivery zone that contains
    them, resolving overlaps by delivery_zones.priority.

    Both geometry columns are SRID 4326. The GiST index on
    delivery_zones(geom) drives a nested-loop index scan; customers are
    walked in id order so the job is restartable and memory-bounded.
    """
    select_sql = """
        SELECT c.id AS customer_id, z.zone_id
        FROM customers c
        LEFT JOIN LATERAL (
            SELECT dz.id AS zone_id
            FROM delivery_zones dz
            WHERE dz.geom && c.geom                 -- GiST bbox prefilter
              AND ST_Contains(dz.geom, c.geom)      -- exact containment
            ORDER BY dz.priority                    -- deterministic overlap winner
            LIMIT 1
        ) z ON TRUE
        WHERE c.id > %(last_id)s
        ORDER BY c.id
        LIMIT %(batch)s;
    """
    update_sql = """
        UPDATE customers AS c
        SET zone_id = data.zone_id
        FROM (SELECT UNNEST(%(ids)s::bigint[])  AS id,
                     UNNEST(%(zones)s::bigint[]) AS zone_id) AS data
        WHERE c.id = data.id;
    """

    total = 0
    with psycopg.connect(DSN, row_factory=dict_row) as conn:
        last_id = 0
        while True:
            with conn.cursor() as cur:
                cur.execute(select_sql, {"last_id": last_id, "batch": BATCH})
                rows = cur.fetchall()

            if not rows:
                break

            ids = [r["customer_id"] for r in rows]
            zones = [r["zone_id"] for r in rows]   # None where no zone contains the point

            with conn.cursor() as cur:
                cur.execute(update_sql, {"ids": ids, "zones": zones})
            conn.commit()      # incremental commit → restartable, bounded WAL

            last_id = ids[-1]
            total += len(rows)
            print(f"assigned through customer id {last_id} ({total} total)")

    return total


if __name__ == "__main__":
    print(f"done: {assign_zones()} customers processed")

The dz.geom && c.geom clause is the cheap gate: the GiST index returns only zones whose bounding box overlaps the point, typically one or two candidates. ST_Contains then runs the exact test against that tiny set instead of all four thousand polygons. The LATERAL ... LIMIT 1 guarantees one row per customer even inside overlapping zones, and the LEFT JOIN keeps customers who fall outside every zone (their zone_id comes back NULL).

The keyset predicate c.id > last_id walks the customer table in primary-key order rather than holding one giant transaction open. Each slice commits independently, so a crash resumes from the last committed id, work-in-flight memory stays flat regardless of table size, and WAL does not balloon into a single multi-gigabyte transaction. This mirrors the streaming discipline described in batch processing spatial joins in Python, applied to a write-back workload.

Configuration and Tuning Knobs

Because the polygon table is small and mostly static, the planner should choose it as the inner relation of the nested loop and index-probe it once per customer. If instead it hashes the customers, the LATERAL LIMIT 1 is doing its job of forcing a per-point lookup. Validate with EXPLAIN before trusting a full run.

Verification Steps

1. Confirm the plan is a nested loop over an index scan, not a sequential scan:

EXPLAIN (ANALYZE, BUFFERS)
SELECT c.id, z.zone_id
FROM customers c
LEFT JOIN LATERAL (
    SELECT dz.id AS zone_id
    FROM delivery_zones dz
    WHERE dz.geom && c.geom
      AND ST_Contains(dz.geom, c.geom)
    ORDER BY dz.priority
    LIMIT 1
) z ON TRUE
WHERE c.id > 0
ORDER BY c.id
LIMIT 5000;

Look for Index Scan using idx_delivery_zones_geom_gist inside a Nested Loop. A Seq Scan on delivery_zones means the SRIDs differ, the index is missing, or statistics are stale — run ANALYZE delivery_zones and re-check.

2. Reconcile assigned versus contained counts — the number of customers with a non-null zone must match a direct containment count:

def verify_assignment(conn):
    with conn.cursor() as cur:
        cur.execute("SELECT COUNT(*) FROM customers WHERE zone_id IS NOT NULL")
        assigned = cur.fetchone()[0]
        cur.execute("""
            SELECT COUNT(*) FROM customers c
            WHERE EXISTS (
                SELECT 1 FROM delivery_zones dz
                WHERE dz.geom && c.geom AND ST_Contains(dz.geom, c.geom)
            )
        """)
        containable = cur.fetchone()[0]
    assert assigned == containable, f"{assigned} assigned vs {containable} containable"
    print(f"verified: {assigned} customers assigned to a zone")

3. Spot-check a known point against ST_AsText of its assigned polygon to confirm the containment rule (strict interior vs boundary) behaves as intended for edge cases.

Gotchas Checklist

• SRID mismatch returns false, not an error. If customers.geom is SRID 4326 and delivery_zones.geom is SRID 3857, ST_Contains compares raw coordinates in different units and quietly matches nothing. Confirm with SELECT DISTINCT ST_SRID(geom) FROM delivery_zones before the run and reproject one side with ST_Transform if they differ.

• Overlapping zones duplicate points without the LATERAL LIMIT 1. A plain JOIN ... ON ST_Contains emits one row per matching polygon. Any customer inside two zones produces two rows, and the write-back becomes non-deterministic. The ORDER BY priority LIMIT 1 inside the lateral collapses that to a single, defined winner.

• Boundary points depend on the predicate. ST_Contains(polygon, point) is false for a point lying exactly on the polygon edge. If shared borders between adjacent zones matter — customers on a street that divides two areas — switch to ST_Intersects or ST_Covers and accept that a border point could match both zones (the priority tiebreak still resolves it).

• Invalid polygons silently drop matches. A self-intersecting or unclosed delivery-zone polygon can make ST_Contains return false for points that are visually inside it. Run SELECT id FROM delivery_zones WHERE NOT ST_IsValid(geom) and repair with ST_MakeValid before the join.

• fetchall() on the full join is the memory trap. Never materialise the whole assignment at once. The keyset loop above keeps in-flight rows to one BATCH and commits incrementally; see advanced GiST indexing and optimization for keeping the polygon index lean as zones are edited.

Frequently Asked Questions

Does ST_Contains use the GiST index on its own?

Yes. ST_Contains and ST_Intersects inject an internal bounding-box && operator that the GiST index resolves, so a separate && clause is not strictly required. Writing the && explicitly documents the intent and guarantees the prefilter even if a future refactor swaps the predicate for one that lacks the built-in index support.

Should I use ST_Contains or ST_Intersects for point-in-polygon assignment?

For a point that must fall strictly inside a zone, ST_Contains(polygon, point) is correct but returns false for a point exactly on the boundary. If a customer sitting on a shared border between two zones must still be assigned, use ST_Intersects or ST_Covers, which include the boundary. Pick one rule and apply it consistently or points on shared edges become duplicates.

How do I stop overlapping zones from producing duplicate rows?

A plain join emits one row per matching polygon, so a point inside two overlapping delivery zones appears twice. Wrap the polygon lookup in a LATERAL subquery with ORDER BY priority LIMIT 1 so each point resolves to exactly one zone, chosen deterministically by a priority column.

Related Topics

• Spatial Joins — parent topic: join strategies, predicate selection, and index interaction
• Batch Processing Spatial Joins in Python — server-side cursor streaming for read-heavy join workloads
• Bounding Box Filtering — how the && operator prunes candidates before the exact predicate
• Advanced GiST Indexing & Optimization — polygon index maintenance, bloat detection, and cost tuning