Handling Polygon Overlaps in Quadtree Partitions

Real-time location pipelines powering IoT telemetry routing, ride-hailing dispatch, and logistics geofencing hit a hard correctness ceiling the moment a single geofence polygon spans more than one quadrant. A quadtree gives predictable logarithmic descent and excellent cache locality, but its rigid four-way decomposition fractures any geometry that crosses a partition edge — and the naive fix, copying the polygon into every node it touches, quietly converts a performance optimization into a duplicate-event generator. This page narrows the deterministic-descent model from the Quadtree vs R-Tree Performance Analysis down to one concrete sub-problem: keeping a polygon that straddles quadrant boundaries indexed exactly once. It sits inside the broader spatial index lookup contract, where every containment check must hold sub-10ms P99 no matter how many partitions a zone overlaps.

The trap is treating insertion as a clipping problem when it is really a reference routing problem. Clip-and-duplicate multiplies intersection work, fragments the heap with redundant GEOS handles, and breaks idempotency — one coordinate crossing a quadrant edge fires the same webhook twice because two nodes each evaluate their own copy of the geometry. The fix is to stop storing geometry in the tree at all.

Concept & Specification #

A quadtree node owns a square region of coordinate space and, on overflow, subdivides into four equal children. A polygon $P$ “overlaps” a partition when its bounding box intersects more than one leaf region. Let $b(P)$ be the axis-aligned bounding box of $P$ and let $L$ be the set of leaf regions. The set of leaves a polygon touches is

$$\text{cover}(P)=,\ell \in L:b(\ell )\cap b(P)\ne \varnothing ,$$

Under clip-and-duplicate, total stored geometry grows as ${\sum }_P|\text{cover}(P)|$, which is super-linear in zone size: a city-scale zone at leaf resolution $r$ can land in dozens of leaves, and each copy carries a full GEOS handle plus its WKB buffer. The architecture below caps stored geometry at exactly $|P|$ — one copy — by separating what the geometry is from where the tree routes to it.

The model has three invariants:

1. Reference-only nodes. A leaf stores (uuid, bbox) tuples — never raw geometry. Geometry lives once in a canonical registry keyed by a stable UUID.
2. Deterministic boundary snapping. Every vertex is snapped to a fixed tolerance grid before insertion, so floating-point drift cannot route the same polygon to different leaves across restarts.
3. Deferred exact evaluation. The tree returns candidate UUIDs via fast bounding-box overlap; the precise point-in-polygon test runs once, off the hot path, against the single canonical geometry.

Step-by-Step Implementation #

Prerequisites: Python 3.11+, shapely>=2.0 (vectorized, releases the GIL on bulk predicates), and coordinates as plain float tuples or a numpy.float64 array of shape (N, 2). Store the authoritative geometry as shapely Polygon/MultiPolygon once; everywhere else route by uuid.UUID and a 4-tuple bbox to avoid per-node object allocation.

1. Build the canonical registry. Each polygon is stored exactly once, keyed by a stable UUID, alongside a precomputed WKB buffer and prepared geometry for fast repeated tests.

from __future__ import annotations

import uuid
from dataclasses import dataclass

from shapely import Polygon, prepare, to_wkb
from shapely.geometry.base import BaseGeometry


@dataclass(frozen=True, slots=True)
class CanonicalZone:
    """A geofence geometry stored once and referenced by UUID."""

    zone_id: uuid.UUID
    geom: BaseGeometry
    wkb: bytes
    bbox: tuple[float, float, float, float]  # (minx, miny, maxx, maxy)


def register(geom: Polygon) -> CanonicalZone:
    prepare(geom)  # cache the GEOS intersection index in-place
    return CanonicalZone(
        zone_id=uuid.uuid4(),
        geom=geom,
        wkb=to_wkb(geom),
        bbox=geom.bounds,
    )

2. Snap vertices deterministically before insertion. Quantizing every coordinate to a fixed grid means a vertex sitting on a quadrant edge always rounds the same way, so the polygon’s cover() set is stable across restarts and replicas.

import numpy as np
from numpy.typing import NDArray

EPS = 1e-7  # ~1.1 cm at the equator


def snap(coords: NDArray[np.float64], eps: float = EPS) -> NDArray[np.float64]:
    # Round to the tolerance grid; identical input -> identical partition.
    return np.round(coords / eps) * eps

Gotcha: snap once, at the ingestion boundary, and store the snapped geometry as canonical. Snapping at query time but not at insert time reintroduces the drift you were trying to remove.

3. Insert references, never geometry. When a zone’s bbox overlaps several leaves, each leaf receives the same (uuid, bbox) card. The geometry is untouched and uncopied.

class QuadNode:
    __slots__ = ("region", "refs", "children")

    def __init__(self, region: tuple[float, float, float, float]) -> None:
        self.region = region
        self.refs: list[tuple[uuid.UUID, tuple[float, float, float, float]]] = []
        self.children: list[QuadNode] | None = None

    def insert(self, zone: CanonicalZone, cap: int = 128) -> None:
        if not _bbox_overlap(self.region, zone.bbox):
            return
        if self.children is None:
            self.refs.append((zone.zone_id, zone.bbox))
            if len(self.refs) > cap:
                self._subdivide(cap)
        else:
            for child in self.children:
                child.insert(zone, cap)

Gotcha: a zone legitimately appears in several leaves — that is reference fan-out, not duplication. The byte cost is one UUID (16 bytes) plus a 4-float bbox per leaf, not a whole geometry.

4. Defer the exact test. Query collects candidate UUIDs by bbox overlap, then runs one precise containment check per distinct candidate against the single canonical geometry — deduplicated by UUID, so a point near a boundary is evaluated once, not once per overlapping leaf.

from shapely import Point, intersects


def locate(
    root: QuadNode, x: float, y: float, registry: dict[uuid.UUID, CanonicalZone]
) -> list[uuid.UUID]:
    seen: set[uuid.UUID] = set()
    pt = Point(x, y)
    hits: list[uuid.UUID] = []
    for zid in _candidates(root, x, y):  # bbox-filtered UUIDs across leaves
        if zid in seen:
            continue  # idempotency: evaluate each zone exactly once
        seen.add(zid)
        if intersects(registry[zid].geom, pt):  # prepared geom -> fast path
            hits.append(zid)
    return hits

5. Run exact tests in a batched worker pool. Keep ingestion threads doing only bbox routing; hand deferred intersects batches to a pool so GEOS work overlaps with I/O and never blocks the descent path.

import asyncio
from concurrent.futures import ProcessPoolExecutor

from shapely import from_wkb, intersects, points


def _batch_contains(wkbs: list[bytes], xs: list[float], ys: list[float]) -> list[bool]:
    geoms = [from_wkb(b) for b in wkbs]
    pts = points(xs, ys)  # vectorized; releases the GIL inside GEOS
    return [bool(intersects(g, p)) for g, p in zip(geoms, pts)]


async def resolve(pool: ProcessPoolExecutor, work: list[bytes], xs, ys) -> list[bool]:
    loop = asyncio.get_running_loop()
    return await loop.run_in_executor(pool, _batch_contains, work, xs, ys)

Gotcha: pass WKB across the process boundary, not shapely objects. GEOS handles do not pickle cleanly, and re-hydrating from WKB in the child is cheaper than the alternative.

Benchmark / Verification #

The figures below come from a 4-core CPython 3.11 worker, 12k municipal geofence zones (median 180 vertices, ~9% of zones spanning 2+ quadrants), and a sustained 50k coordinate events/sec. The “clip-and-duplicate” row copies geometry into every overlapping leaf; the “reference + deferred” row is the design above.

The duplicate rate collapses to zero because UUID deduplication guarantees one evaluation per zone regardless of cover() size, and the P99 drops ~7× because the descent path no longer instantiates or intersects redundant geometry. To verify the index is genuinely storing references only, assert that total stored geometry equals the registry size, and confirm boundary determinism by re-snapping and re-inserting a known straddling zone twice:

total_refs = sum(len(n.refs) for n in _walk(root))
assert sum(1 for _ in _walk(root)) >= 1
# Each zone is one registry entry no matter how many leaves reference it.
assert len(registry) <= total_refs
# Determinism: identical input snaps to an identical cover set.
assert cover(zone_a) == cover(register(Polygon(snap(coords_a))))

A py-spy dump during peak ingest should show descent threads in _bbox_overlap/_candidates, never in GEOS intersects on the hot path — exact tests belong in the pool. If RSS climbs past ~3× the registry, geometry is leaking into nodes; the bulk-load discipline that keeps that residency bounded is covered in Optimizing R-Tree Bulk Loads for Real-Time Ingestion.

Failure Modes & Edge Cases #

• Floating-point edge drift. A vertex landing within $\epsilon$ of a quadrant edge routes to different leaves on different runs without snapping, breaking cache locality and invalidating precomputed candidate sets. Snap deterministically at ingest (step 2) and treat the snapped geometry as canonical — never snap on one path but not the other.
• GEOS handle leaks under churn. In CPython, shapely 1.x geometries that outlive their C-API context fail to release the underlying GEOS handle, so RSS grows even as Python refcounts drop. Migrate to shapely>=2.0, operate on WKB buffers across thread/process boundaries, and correlate gc.get_stats() collection counts with RSS spikes; the broader allocation discipline lives in memory-constrained spatial processing.
• GIL contention on synchronous tests. Calling intersects inline from several descent threads serializes them on the GIL and the GEOS lock, producing a runaway tail. Keep exact tests in the batched pool (step 5) so the C work releases the GIL; the boundary trade-offs mirror those in async Python execution patterns for spatial math.
• Degenerate geometry. Self-intersecting rings, NaN/infinite coordinates, and empty MultiPolygon parts poison bounding boxes and silently break overlap tests. Validate at the registry boundary — geom.is_valid and np.isfinite(coords).all() — and make_valid or drop bad zones before insertion, never after publishing.
• Pathological cover fan-out. A single very large zone overlapping hundreds of leaves inflates candidate sets even though it is one geometry. Cap leaf references at $c$ to bound per-query fan-out, and for super-zones fall back to a coarse cell grid (H3 or S2) so descent does not enumerate every leaf the zone touches.
• Stale snapshots pinning memory. A long-running reader holding an index snapshot across many rebuilds prevents reference-count reclamation and drifts peak memory past the expected 2–3×. Cap retained versions and bound reader lifetime — the same swap-and-reclaim model detailed in Thread-Safe Spatial Index Updates in Python.

Related #

• Quadtree vs R-Tree Performance Analysis — parent reference: how the two tree families fail in opposite ways and when quadrant descent beats MBR overlap.
• Optimizing R-Tree Bulk Loads for Real-Time Ingestion — the sibling write-path problem: building a balanced index without per-insert churn.
• Thread-Safe Spatial Index Updates in Python — publishing rebuilt partitions atomically while readers stay lock-free.
• Memory Footprint of Streaming Polygon Indexes — keeping the canonical registry and its WKB buffers inside the container memory budget.
• Up one level: Spatial Indexing for Real-Time Checks — the index-structure and profiling contract every partitioning strategy must satisfy.