Spatial Indexing on Constrained Devices

Within the Local Spatial Processing Patterns guide, this page covers the structure that decides whether a gateway can answer “which zones contain this point?” in microseconds or in milliseconds — the spatial index. Everything else in this section assumes one exists: the funnel in on-device geometry filtering needs a candidate set to filter, and the joins in spatial joins in constrained environments need somewhere to look features up. This page is about building and holding that index on hardware that cannot afford the desktop answer.

The desktop answer is a dynamic R-tree with insert, delete and rebalance, held in a language runtime’s heap. On a gateway that structure has three problems: its memory footprint is several times the data it indexes, its pointer chasing defeats a small cache, and rebalancing produces unpredictable pauses in a pipeline that has a deadline. The edge answer is almost always a static index — built once, laid out contiguously, and replaced wholesale rather than mutated.

Constraint mapping #

The row that surprises people is cache size. On a Cortex-A53 with a 512 KB L2, a pointer-based R-tree over 40 000 features misses cache on nearly every node traversal, and each miss costs 80–150 cycles. The same tree packed into a contiguous array traverses several levels within a handful of cache lines. The measured difference is routinely three to five times, which is larger than the difference between the algorithms themselves.

Implementation: a static grid built once at load #

For a zone set that fits a bounded area with reasonably even coverage, nothing beats a flat grid. The build is a single pass, the query is arithmetic, and the whole structure is two arrays.

# static_grid.py — build-once spatial grid over feature envelopes.
# No allocation in the query path; the two arrays are sized at build time.
# Threading: build on one thread, then treat as immutable and share freely.
import array
import math


class StaticGrid:
    __slots__ = ("min_x", "min_y", "cell", "nx", "ny", "starts", "items")

    def __init__(self, envelopes, cell_size):
        """envelopes: sequence of (min_x, min_y, max_x, max_y) in projected metres."""
        self.cell = float(cell_size)
        self.min_x = min(e[0] for e in envelopes)
        self.min_y = min(e[1] for e in envelopes)
        max_x = max(e[2] for e in envelopes)
        max_y = max(e[3] for e in envelopes)
        self.nx = max(1, int((max_x - self.min_x) / self.cell) + 1)
        self.ny = max(1, int((max_y - self.min_y) / self.cell) + 1)

        # Pass 1: count occupants per cell so the CSR arrays can be sized exactly.
        counts = array.array("i", [0]) * 0
        counts = array.array("i", bytes(4 * (self.nx * self.ny + 1)))
        spans = []
        for feature_id, (x0, y0, x1, y1) in enumerate(envelopes):
            cx0, cy0 = self._cell_of(x0, y0)
            cx1, cy1 = self._cell_of(x1, y1)
            spans.append((feature_id, cx0, cy0, cx1, cy1))
            for cy in range(cy0, cy1 + 1):
                base = cy * self.nx
                for cx in range(cx0, cx1 + 1):
                    counts[base + cx + 1] += 1

        # Prefix sum → start offsets (compressed sparse row layout).
        for i in range(1, len(counts)):
            counts[i] += counts[i - 1]
        self.starts = counts
        self.items = array.array("i", bytes(4 * counts[-1]))

        # Pass 2: place feature ids. cursor walks a copy of the start offsets.
        cursor = array.array("i", counts)
        for feature_id, cx0, cy0, cx1, cy1 in spans:
            for cy in range(cy0, cy1 + 1):
                base = cy * self.nx
                for cx in range(cx0, cx1 + 1):
                    slot = base + cx
                    self.items[cursor[slot]] = feature_id
                    cursor[slot] += 1

    def _cell_of(self, x, y):
        cx = min(self.nx - 1, max(0, int((x - self.min_x) / self.cell)))
        cy = min(self.ny - 1, max(0, int((y - self.min_y) / self.cell)))
        return cx, cy

    def query_point(self, x, y):
        """Candidate feature ids whose envelope covers this cell. No allocation
        beyond the returned memoryview slice."""
        cx, cy = self._cell_of(x, y)
        slot = cy * self.nx + cx
        return self.items[self.starts[slot]:self.starts[slot + 1]]

    def memory_bytes(self):
        return (len(self.starts) + len(self.items)) * 4

The layout is compressed sparse row: one integer per cell holding the start offset, and one integer per (feature, cell) pair. For 40 zones over a 20 km area at 500 m cells that is 1 600 cells and roughly 260 occupancy entries — under 8 KB, entirely cache-resident, with a query that is two multiplications and a slice. The structure has no Python objects in the hot path at all, which is what keeps the collector out of it.

Implementation: when the grid stops fitting #

Grids fail on skew. A zone set that includes both a 50 m loading bay and a 200 km pipeline corridor cannot be served by a single cell size: fine cells make the corridor occupy thousands of them, coarse cells make the bay share with everything else. Two escapes exist.

The first is a two-level grid: a coarse grid whose dense cells hold a nested fine grid. This preserves the arithmetic query and handles moderate skew, at the cost of a branch per lookup and a build that has to decide which cells to subdivide.

The second is a packed Hilbert R-tree, which adapts to the data’s own distribution instead of imposing a uniform one. It is built by sorting features along a space-filling curve and packing them into fixed-size nodes bottom-up, producing a contiguous array with no pointers — a structure that suits flash and cache equally well and can be memory-mapped straight from a file. The construction and the query loop are covered in packed Hilbert R-tree in a static buffer.

Where the index lives: heap, file, or both #

An index has to exist somewhere, and on a gateway the choice between the heap and a memory-mapped file changes far more than it appears to. It decides how long a restart takes, whether the index survives a watchdog reset, and whether the kernel is allowed to evict it under pressure.

Holding the index on the heap is the default and the least flexible option. It costs a build at every start-up, it is destroyed by every restart, and it is unevictable — the kernel cannot reclaim it under memory pressure, so it is the pages that force everything else out. Its advantage is simplicity and the fastest possible query, since nothing can fault.

Building the index at provisioning and shipping it as a file the device memory-maps inverts every one of those properties. Start-up costs one open and one mmap regardless of index size. A watchdog reset does not rebuild anything. The pages are file-backed and clean, so under pressure the kernel evicts them and re-faults them later instead of killing the process, which converts a hard failure into a soft slowdown. The cost is that the first query after a restart faults pages in from flash, and that the index cannot reflect a reference layer the device received itself.

The hybrid that most mature deployments arrive at is a mapped base index for the layer shipped at provisioning, plus a small heap-resident overlay for features synced since. The query consults both and merges the candidate sets, which costs one extra traversal against a structure that is usually two orders of magnitude smaller than the base. Periodically — at a maintenance window, or when the overlay exceeds a size threshold — the two are merged into a new base file and the overlay is emptied.

That structure has a property worth naming: the expensive operation, rebuilding the base, happens on a schedule the device chooses, while the cheap operation, extending the overlay, happens whenever data arrives. Nothing in the ingestion path ever waits for a full rebuild, and the memory spike covered in the rebuild guide occurs at a moment when the device is idle rather than at the moment a sync lands.

Configuration and tuning #

• Cell size on a grid: start at the median feature diagonal, then measure candidates returned per probe. The band that works is usually one to four times the median feature size; the full trade is charted in the spatial joins guide.
• Node size on a packed R-tree: 16 entries is the usual sweet spot on ARM, because a 16-entry node of four floats each is 256 bytes — four cache lines — and the branching factor keeps the tree shallow.
• Build trigger: rebuild on reference-layer change, never on a timer. An index rebuilt every hour against unchanged data is pure cost and a recurring memory spike.
• Persistence: prefer building the index at provisioning and shipping it as a file the device memory-maps. It removes the build cost from every boot and makes the rebuild peak someone else’s problem.
• Fallback: keep the linear scan implementation and a switch to it. When an index build fails on a malformed layer, degrading to a slower correct answer beats failing to answer at all.

The case for keeping the linear scan #

Every index in this section replaces a loop over the features with something cleverer, and every deployment should keep that loop. Not as a comment about what the index replaced — as running, tested, switchable code.

Three situations make it the correct choice rather than a fallback. Below roughly two hundred features a scan over a contiguous array of envelopes beats every index on this page, because the whole array fits in L1 and the branch predictor learns it; the index’s traversal costs more than the comparisons it avoids. During an index rebuild, a scan lets queries continue against the new layer immediately rather than waiting for a build to finish. And when an index build fails — a malformed layer, a memory budget refused — a scan is the difference between a device that answers slowly and one that does not answer.

The implementation is short enough that keeping it costs nothing:

def scan_candidates(envelopes, x, y, out):
    """Linear envelope test. `envelopes` is a flat array of 4 floats each;
    `out` is reused across calls so the hot path allocates nothing."""
    out.clear()
    for i in range(0, len(envelopes), 4):
        if (envelopes[i] <= x <= envelopes[i + 2]
                and envelopes[i + 1] <= y <= envelopes[i + 3]):
            out.append(i >> 2)
    return out

What makes it valuable is that it is also the reference implementation. Any index in this section can be validated against it directly: run both over the same query set, and assert the candidate sets agree. That test finds coordinate-ordering mistakes, off-by-one errors in cell arithmetic and traversal bugs, all of which otherwise present as intermittently missing results rather than as failures.

Wire the switch as configuration rather than as a code path chosen at build time, and export which mode is active in the health snapshot. A device that has quietly fallen back to scanning because its index build failed is a device with a latency problem and a working answer, and the operator needs to know which of those they are looking at.

Verification and field diagnostics #

Instrument three numbers and the index stops being opaque. Candidates per probe is the index’s actual selectivity — if it drifts upward after a layer update, the new features are shaped differently than the cell size assumes. Build time and build peak tell you whether a restart will meet its deadline and whether the rebuild will fit in memory. Query time percentiles catch the case where a well-behaved median hides a tail caused by one enormous cell.

A useful field command dumps the occupancy histogram: number of cells holding zero, one, two to four, five to sixteen, and more than sixteen features. A healthy grid is dominated by zero and one; a long tail past sixteen means the cell size no longer matches the data, and it means it before any latency alert fires.

Failure modes specific to this pattern #

A short checklist before shipping an index #

Five questions, asked once, prevent most of the failures catalogued above.

Have you measured candidates per probe against real traffic? Not against a synthetic uniform sample — against a recorded hour from a device in the field. Real traffic concentrates, and an index tuned against uniform queries over-returns for the corridor the fleet actually drives.

Do you know the build peak, and is it inside the budget on the smallest device in the fleet? The rebuild is the largest allocation the process makes, and the smallest device is the one that will meet it first.

Can the index be replaced while queries run? If the answer involves a lock held across the rebuild, queries stall for the duration; if it involves mutating in place, some query will see a half-built structure.

What happens if the build fails? A device with no index and no fallback answers nothing. A device that silently keeps a stale index answers wrongly. Only the fallback-to-scan path degrades honestly.

Is the active mode visible from the console? Grid or tree, base-plus-overlay or single, indexed or scanning — a field diagnosis that has to infer this from latency is a diagnosis that takes a day instead of a minute.

None of the five requires new code beyond what this section already describes. They require deciding the answers before the deployment rather than during the incident.

Related #

• Packed Hilbert R-tree in a static buffer — the contiguous, pointer-free tree and its query loop.
• H3 vs geohash cell indexing for zone lookups — when the index key has to leave the process.
• Rebuilding spatial indexes without a memory spike — the replacement path and its peak.
• Spatial Joins in Constrained Environments — what the candidate sets are used for.
• On-Device Geometry Filtering — the exact predicates that run after the index narrows the field.