Spatial Joins in Constrained Environments

Within the Local Spatial Processing Patterns framework, this guide covers how to correlate streaming sensor coordinates with static reference geometry on gateways that have no server-class memory and no unbounded compute budget.

When an IoT gateway must relate live telemetry — GPS fixes, asset pings, NMEA tuples — to fixed reference layers such as zoning boundaries, asset footprints, or exclusion zones, the desktop spatial-join paradigm collapses. Loading both the sensor stream and the reference collection into RAM and handing them to a dynamic R-tree is exactly what triggers OOM kills on a Raspberry Pi 4 or an industrial Linux gateway. Reliable execution means trading monolithic, in-memory joins for a streaming, constraint-aware pipeline: a static index built once, a cheap candidate probe per event, and a precise predicate evaluated only on the survivors.

Streaming join: prune by bounding box, then confirm with a prepared GEOS predicate; matched pairs flush in chunks while misses advance to the next candidate.

Constraint Mapping #

Every decision in a constrained spatial join is dictated by a specific hardware limit, so it pays to name which limit drives which design choice before writing geometry code. These same ceilings are treated as first-class parameters in the device constraints and resource limits reference; the table below maps each one to the join behaviour it governs.

The recurring theme is that a join on the edge is I/O-bound and allocation-sensitive, not CPU-bound in the textbook sense. A standard point-in-polygon over a 500 MB GeoJSON file does not fail because the geometry math is slow; it fails because deserialization duplicates every coordinate array into Python objects and the garbage collector thrashes trying to keep up. The fix is structural: never hold the whole reference layer as live objects, and never let a geometry operation block the thread reading the sensor.

Implementation 1 — A Static Grid Index Built Once #

The primary technique replaces the dynamic tree with a uniform grid (a flat spatial hash) computed at startup from the bounding boxes of the reference features. Because the grid never grows during operation, there are no per-event allocations in the hot path and nothing for the generational collector to rescan. Each reference feature is registered into every grid cell its envelope overlaps; a sensor event then probes the single cell that contains its coordinate and tests only that short candidate list.

Grid cells are expressed in the same units as the incoming coordinates. If the stream arrives in WGS84 degrees, a CELL of 0.01° is roughly 1.1 km of latitude — keep the unit choice deliberate, because the relationship between degrees and metres varies with latitude, a subtlety covered in coordinate reference systems at the edge.

# grid_join.py — static uniform-grid spatial join for constrained gateways.
#
# Memory model: the grid is built ONCE at startup and is immutable thereafter.
# There are no allocations in the per-event hot path, so the generational GC has
# nothing to sweep. gc.freeze() moves the grid out of reach of later collections.
# Threading model: build runs at init; stream_join() runs on the asyncio loop
# thread and offloads only the precise predicate (see geos_ffi.py) to a worker.

import gc
from array import array

CELL = 0.01  # grid cell edge, in the coordinate units of the stream (degrees)

def _key(lon, lat):
    # Integer cell coordinate. Integer keys hash faster than floats and never
    # drift, so the same point always lands in the same bucket.
    return (int(lon // CELL), int(lat // CELL))

def build_grid(features):
    """features: iterable of (feature_id, minx, miny, maxx, maxy) envelopes.
    Returns dict[cell] -> array('i', feature_ids). Call once, then it is frozen.
    Feature IDs are stored in a typed array('i'), not a Python list, to keep the
    index a few bytes per entry instead of a boxed-int object each."""
    grid = {}
    for fid, minx, miny, maxx, maxy in features:
        cx0, cy0 = _key(minx, miny)
        cx1, cy1 = _key(maxx, maxy)
        for cx in range(cx0, cx1 + 1):
            for cy in range(cy0, cy1 + 1):
                grid.setdefault((cx, cy), array('i')).append(fid)
    gc.collect()      # settle the build allocations
    gc.freeze()       # exclude the static grid from future GC generations
    return grid

def candidates(grid, lon, lat):
    # O(1) bucket probe. Returns a short candidate tuple, or () for empty cells.
    return grid.get(_key(lon, lat), ())

The join itself is a generator. It pulls one event at a time and yields matched (event_id, feature_id) pairs lazily, so at most one event and one candidate list are live at any moment. This is what keeps heap occupancy flat regardless of how long the stream runs.

def stream_join(grid, points, precise_contains):
    """points: generator of (event_id, lon, lat) from the sensor reader.
    precise_contains(fid, lon, lat) -> bool, the exact GEOS check.
    Yields (event_id, fid) lazily; one event is in flight at a time."""
    for event_id, lon, lat in points:
        for fid in candidates(grid, lon, lat):
            if precise_contains(fid, lon, lat):
                yield (event_id, fid)
                break          # first containing feature wins; stop scanning

Reference envelopes feed build_grid from a streaming parser — never from json.load. Building the (fid, minx, miny, maxx, maxy) tuples by walking the file with an incremental tokenizer is what keeps the build within the RAM ceiling; the mechanics of that parse are the subject of reducing RAM usage for GeoJSON parsing on Raspberry Pi.

Implementation 2 — Precise Predicates via GEOS FFI #

The grid probe is a coarse filter; it returns features whose envelope overlaps the event’s cell, not features that actually contain the point. The exact test is delegated to GEOS, but through direct foreign-function calls rather than a Python geometry wrapper. Bypassing the object-heavy layer and calling the reentrant C API through ctypes cuts per-event overhead by 60–80% on Cortex-A SoCs. The same FFI escalation underpins implementing polygon containment checks in C for the pure-geometry case.

Load the shared library once, hold a single reentrant context, and pre-build a GEOSPrepared version of each reference geometry so repeated point tests reuse the cached internal index. The exact signatures come from the official GEOS C API reference; declaring them explicitly to ctypes matters on aarch64, where a defaulted int restype truncates a 64-bit pointer.

# geos_ffi.py — minimal ctypes binding to GEOS prepared-geometry predicates.
# Load the library ONCE and reuse a single reentrant context per worker thread.
# The GEOS context handle is NOT thread-safe to share, so create one per thread
# in the executor (see Configuration & Tuning) rather than a single global.

import ctypes as C

_geos = C.CDLL("libgeos_c.so.1")
_geos.GEOS_init_r.restype = C.c_void_p
_geos.GEOSPrepare_r.argtypes = [C.c_void_p, C.c_void_p]
_geos.GEOSPrepare_r.restype = C.c_void_p
_geos.GEOSPreparedContains_r.argtypes = [C.c_void_p, C.c_void_p, C.c_void_p]
_geos.GEOSPreparedContains_r.restype = C.c_char  # 1 = true, 0 = false, 2 = error

class PreparedIndex:
    """Maps feature_id -> prepared GEOS geometry. Prepared geometries cache an
    internal monotone-chain index, so the second and later point tests against a
    polygon are far cheaper than the first."""
    def __init__(self):
        self.ctx = _geos.GEOS_init_r()
        self._prepared = {}     # fid -> GEOSPreparedGeometry*
        self._point = None      # reused mutable point; see contains()

    def add(self, fid, geom_ptr):
        # geom_ptr is a GEOSGeometry* built once from WKB at load time.
        self._prepared[fid] = _geos.GEOSPrepare_r(self.ctx, geom_ptr)

    def contains(self, fid, lon, lat):
        # Mutate one reused point geometry instead of allocating per event.
        # Building a fresh GEOSGeom_createPoint per call would allocate in C on
        # every fix; here we update the coordinate sequence of a single point.
        _set_point_xy(self.ctx, self._point, lon, lat)
        return _geos.GEOSPreparedContains_r(
            self.ctx, self._prepared[fid], self._point) == b"\x01"

The contains method deliberately mutates one reused point geometry rather than allocating a fresh GEOS point per fix. That single decision removes the largest remaining source of per-event allocation, which on a duty-cycled gateway is the difference between a flat RSS curve and a slow climb into the swap-or-die zone. When the precise test is too expensive to run inline — for example during a 100 Hz GNSS burst — fall back to bucket-only matching and let threshold-based event mapping decide which events deserve the exact check.

Configuration & Tuning #

The grid is only as good as its cell size, and cell size is the single most consequential knob. Too coarse and every bucket returns dozens of candidates, pushing the GEOS predicate count up until the CPU saturates; too fine and large polygons register into thousands of cells, inflating the index past the RAM ceiling. Tune CELL to the median feature size, not the largest.

Because the GEOS reentrant context is not safe to share across threads, give each worker in the pool its own PreparedIndex (or at least its own context handle), then offload the predicate from the asyncio loop with asyncio.to_thread() or a bounded concurrent.futures.ThreadPoolExecutor. This keeps the GIL released during the C call so the sensor reader keeps draining its socket. The structure of those offloads is detailed in async execution for spatial workloads; the broader event-loop rules are in the asyncio documentation. Pin the worker threads to physical cores with os.sched_setaffinity so a big.LITTLE scheduler does not migrate a hot predicate onto a throttled little core mid-burst.

If you build GEOS yourself for the target, compile the C wrapper that exposes the prepared predicates with -O2 -fno-exceptions -fvisibility=hidden to keep the symbol table small and the hot path branch-predictable; ship it as a single versioned .so (libgeos_c.so.1) loaded by the binding above.

Verification & Field Diagnostics #

A constrained join can be wrong in two silent ways: it can drop true matches because the grid registration missed a cell, or it can climb in memory because something in the hot path is still allocating. Both are caught with lightweight, always-available instrumentation rather than a profiler.

• Ground-truth replay. Before deployment, replay a recorded coordinate track with known-correct matches through stream_join and assert the yielded pairs match the reference set exactly. This catches CELL-registration bugs that only appear near feature boundaries.
• Resident memory watch. Sample /proc/self/status VmRSS, or psutil.Process().memory_info().rss, on a timer. A flat line confirms the no-allocation hot path is holding; a slow climb points to a per-event allocation that escaped, usually a GEOS point that is not being reused.
• Grid balance metric. Log the length of each candidates() result periodically. A healthy grid returns single-digit candidate counts; a p99 in the dozens means CELL is too coarse for the feature density.
• Staging-only leak hunt. Enable tracemalloc in staging to pin coordinate-buffer leaks to a line number, then disable it in production — its per-allocation bookkeeping is too costly for the gateway to carry continuously.
• Throughput correlation. When telemetry stalls, correlate gateway CPU throttling against join latency; a thermal trip and a latency spike landing together confirm the join is heat-bound, not logic-bound.

Failure Modes Specific to This Pattern #

The grid-plus-FFI join degrades in characteristic ways, and each has a detectable signature and a safe recovery path:

• Reference-layer OOM at startup. Materializing the GeoJSON to build the grid blows the RAM ceiling before the join even runs. Detection: the process is OOM-killed during init, never during streaming. Recovery: feed build_grid from the streaming parser only, and stage envelopes to disk if the feature count is very large.
• Candidate-list explosion. An undersized CELL (relative to feature size) makes buckets return long candidate lists, the GEOS predicate count spikes, and the CPU saturates. Detection: the grid-balance metric climbs while ingest rate is flat. Recovery: rebuild the grid with a larger CELL, or split oversized features.
• Cross-thread GEOS corruption. Sharing one reentrant context across pool workers produces non-deterministic crashes or wrong predicate results. Detection: failures that vanish at max_workers=1. Recovery: one GEOS context (and prepared index) per worker thread.
• Backpressure collapse. When the FFI queue cannot keep up, unbounded buffering bloats memory until the watchdog fires. Detection: queue depth and RSS rise together. Recovery: bound the executor queue, apply exponential backoff on sensor polling, and shed low-priority events via threshold mapping before they reach the join.
• Graceful-degradation handoff. When MemoryMax is approached (RSS past ~85% of cap), switch from exact GEOS containment to bucket-only matching so the service keeps producing approximate results instead of being killed. Flush matched events to the local broker in chunks; the queue and delivery side of that handoff is covered in message queue management at the edge.

Related #

• On-device geometry filtering — the upstream gatekeeper that discards irrelevant geometry before it reaches the join.
• Threshold-based event mapping — decide which matched events warrant the exact predicate and which to shed under load.
• Async execution for spatial workloads — how to offload synchronous FFI predicates without blocking the sensor loop.
• Reducing RAM usage for GeoJSON parsing on Raspberry Pi — the streaming parser that feeds envelopes into the grid build.
• Device constraints and resource limits — the hardware budgets every choice on this page is measured against.