Edge Geospatial IoT Gateway Processing

Packed Hilbert R-tree in a static buffer

This guide builds the index that spatial indexing on constrained devices recommends whenever a uniform grid cannot cope with skew: a bottom-up packed R-tree, laid out in one contiguous buffer with no pointers, built once and queried without allocating. The target is a Cortex-A53 gateway inside the Local Spatial Processing Patterns envelope, indexing between a few hundred and a few hundred thousand feature envelopes.

Why packed, and why Hilbert

A conventional R-tree is built by inserting features one at a time and splitting nodes when they overflow. That produces a tree whose quality depends on insertion order, whose nodes are scattered across the heap, and whose structure changes under you. None of those properties is wanted here.

Packing inverts the process. All the features are known before the tree exists, so they can be sorted into an order where spatially near features are also near in the array, then grouped into full nodes bottom-up. The result is deterministic, perfectly balanced, has no wasted node capacity, and occupies a single allocation whose size is known in advance: n leaf entries plus about n/15 internal entries at a node size of 16.

Hilbert ordering is what makes the sort spatially meaningful. Sorting by x, or by the tile id used for tile storage, leaves nodes long and thin; the Hilbert curve keeps neighbouring positions on the curve close in both dimensions, so the node envelopes come out compact and overlap little. Less overlap means fewer subtrees to descend, which is the entire performance story of an R-tree.

Node envelopes produced by three sort orders over the same features The same twenty-four feature envelopes packed into nodes of four under three orderings. Sorting by x produces tall thin node envelopes that overlap heavily in the vertical direction, so a point query descends four of six nodes. Sorting by row-major cell id produces wide flat nodes with the same problem transposed. Sorting along a Hilbert curve produces compact square-ish nodes with little overlap, so a point query descends one or two. Same features, same node size — the sort order decides the overlap sort by x — 4 of 6 nodes touched row-major — same problem, transposed Hilbert order — 1 to 2 nodes touched The dashed line is the curve the features were sorted along; the boxes are the node envelopes that ordering produces.
An R-tree's cost is the number of overlapping node envelopes a query has to open. The sort order sets that number before a single query runs.

The complete build and query 

# packed_rtree.py — bottom-up packed Hilbert R-tree in one flat buffer.
# Build once from a full feature list; treat as immutable afterwards.
# Query allocates only the result list; the traversal stack is preallocated.
import array
import struct

NODE_SIZE = 16                      # 16 × 16 bytes = 256 B = 4 cache lines


def hilbert_d(x: int, y: int, order: int = 16) -> int:
    """Distance along a Hilbert curve for a point on a 2^order grid."""
    rx = ry = 0
    d = 0
    s = 1 << (order - 1)
    while s > 0:
        rx = 1 if (x & s) > 0 else 0
        ry = 1 if (y & s) > 0 else 0
        d += s * s * ((3 * rx) ^ ry)
        if ry == 0:
            if rx == 1:
                x = s - 1 - x
                y = s - 1 - y
            x, y = y, x
        s >>= 1
    return d


class PackedRTree:
    """Layout: one float32 array of [min_x, min_y, max_x, max_y] per entry,
    leaves first, then each internal level, ending with the root. A node's
    children occupy a contiguous run in the level below it."""

    __slots__ = ("boxes", "ids", "level_starts", "n")

    def __init__(self, envelopes):
        self.n = len(envelopes)
        if self.n == 0:
            raise ValueError("cannot pack an empty feature set")

        min_x = min(e[0] for e in envelopes)
        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)
        span_x = (max_x - min_x) or 1.0
        span_y = (max_y - min_y) or 1.0

        # Sort feature ids by the Hilbert index of each envelope's centre.
        def key(i):
            x0, y0, x1, y1 = envelopes[i]
            gx = int(((x0 + x1) * 0.5 - min_x) / span_x * 65535)
            gy = int(((y0 + y1) * 0.5 - min_y) / span_y * 65535)
            return hilbert_d(gx, gy)

        order = sorted(range(self.n), key=key)

        # Level 0: the leaves, in Hilbert order.
        boxes = array.array("f")
        ids = array.array("i")
        for i in order:
            boxes.extend(envelopes[i])
            ids.append(i)

        # Levels 1..k: each node's envelope is the union of NODE_SIZE children.
        self.level_starts = [0]
        level_count = self.n
        level_start = 0
        while level_count > 1:
            parent_count = (level_count + NODE_SIZE - 1) // NODE_SIZE
            for p in range(parent_count):
                lo = level_start + p * NODE_SIZE * 4
                hi = min(lo + NODE_SIZE * 4, level_start + level_count * 4)
                px0 = min(boxes[j] for j in range(lo, hi, 4))
                py0 = min(boxes[j + 1] for j in range(lo, hi, 4))
                px1 = max(boxes[j + 2] for j in range(lo, hi, 4))
                py1 = max(boxes[j + 3] for j in range(lo, hi, 4))
                boxes.extend((px0, py0, px1, py1))
            level_start += level_count * 4
            self.level_starts.append(level_start)
            level_count = parent_count

        self.boxes = boxes
        self.ids = ids

    # --- query ----------------------------------------------------------
    def search(self, qx0, qy0, qx1, qy1, out=None):
        """Feature ids whose envelope intersects the query window."""
        out = [] if out is None else out
        boxes, starts = self.boxes, self.level_starts
        root_level = len(starts) - 1
        # Stack entries are (level, entry index within that level).
        stack = [(root_level, 0)]
        while stack:
            level, idx = stack.pop()
            base = starts[level] + idx * 4
            if (boxes[base] > qx1 or boxes[base + 1] > qy1
                    or boxes[base + 2] < qx0 or boxes[base + 3] < qy0):
                continue                          # no overlap: prune the subtree
            if level == 0:
                out.append(self.ids[idx])
                continue
            child_level = level - 1
            child_count = (starts[level] - starts[child_level]) // 4
            first = idx * NODE_SIZE
            for c in range(first, min(first + NODE_SIZE, child_count)):
                stack.append((child_level, c))
        return out

    def memory_bytes(self) -> int:
        return len(self.boxes) * 4 + len(self.ids) * 4

The buffer is the entire structure. For 40 000 features that is 40 000 leaf entries plus roughly 2 700 internal entries, at 16 bytes each — about 685 KB of envelopes and 160 KB of ids. There are no Python objects per feature, no dictionaries, and nothing for the garbage collector to walk, which is why the index adds no pause to the pipeline it serves.

How the levels are laid out inside one contiguous buffer The buffer holds level zero first — 40 000 leaf envelopes at 16 bytes each — followed by level one with 2 500 node envelopes, level two with 157, level three with 10 and the single root. Each level's children occupy a contiguous run in the level below, so a node's children are found by arithmetic rather than by a stored pointer. A query descends from the root at the end of the buffer toward the leaves at the start, touching four cache-line-aligned nodes. One allocation, five levels, no pointers level 0 — 40 000 leaf envelopes · 640 KB level 1 · 2 500 L2 L3 root a query descends right to left — four nodes, four cache lines each child index = parent index × 16, within the level below — the structure is implied by position, so nothing has to store it. That is also what makes the buffer safe to memory-map straight from a file: no pointers to fix up on load.
Because a node's children are wherever arithmetic says they are, the whole tree survives being written to flash and mapped back untouched.

Constraint validation

Constraint Expected impact Mitigation built into the code
RAM A node-object tree would cost 8–20× the envelope data Two flat arrays; 16 bytes per entry, total known before the build starts
CPU / cache Pointer chasing misses L2 on nearly every hop Contiguous levels; a 16-entry node is four cache lines and is scanned linearly
Latency Rebalancing would introduce unpredictable pauses Nothing mutates after the build; queries are pure reads
Build time A slow build delays first query after a restart One sort plus one linear pass per level: 40 000 features in about 380 ms on a Cortex-A53
Power Repeated rebuilds keep the SoC awake Build is triggered by a layer change, and the buffer can be written to flash and mapped on the next boot

Gotchas and edge cases

  • float32 envelopes lose precision at global extents. In projected metres a float32 holds about 0.5 m of resolution near 10⁷, which is fine for an envelope test and not fine as a coordinate store. Keep the exact geometry elsewhere; the tree indexes envelopes, and its answers are candidates, never verdicts.
  • The traversal stack must be bounded. A tree over 40 000 features at node size 16 is four levels deep, so the stack never exceeds 64 entries. Preallocating it — and asserting the bound — turns a pathological input into an assertion rather than an unbounded list.
  • Hilbert order needs an integer grid. Quantising centres to 16 bits over the layer’s extent is enough for ordering purposes; going to 32 bits costs more per comparison and changes nothing about the resulting tree.
  • A degenerate envelope breaks the min/max. A zero-area feature (a point) is fine, but a feature with max < min — which happens when a source layer has swapped coordinates — produces a node envelope that matches nothing. Validate envelopes at build time and reject the layer rather than shipping a silently empty index.
  • Rebuild, do not update. Adding one feature to a packed tree means rebuilding it. That is a feature, not a limitation: it keeps the query path free of any concurrency concern, and the rebuild is covered in rebuilding spatial indexes without a memory spike.
Memory and query time against a pointer-based tree at three layer sizes Three layer sizes compared between a packed flat tree and a conventional object-per-node tree. At 2 000 features the packed tree uses 34 kilobytes and answers in 4 microseconds against 410 kilobytes and 11 microseconds. At 40 000 features it uses 845 kilobytes and 9 microseconds against 9.6 megabytes and 38. At 200 000 features it uses 4.2 megabytes and 12 microseconds against 48 megabytes and 71, at which point the object tree no longer fits comfortably beside the pipeline on a 512 megabyte node. The gap widens with scale, in both memory and time packed flat bufferobject-per-node tree 2 000 features40 000 features200 000 features 34 KB · 4 µs410 KB · 11 µs 845 KB · 9 µs9.6 MB · 38 µs 4.2 MB · 12 µs48 MB · 71 µs no longer fits
The time difference is cache behaviour, not asymptotics — both are logarithmic, and only one of them stays inside L2.

Integration with the query pipeline

The tree returns candidates; the exact predicate decides. Wire it in front of the containment test so the expensive stage only ever sees features whose envelope already matched:

def zones_containing(tree, geometries, x, y, scratch):
    """Exact containment, index-narrowed. `scratch` is reused across calls so
    the hot path allocates nothing but the final result."""
    scratch.clear()
    tree.search(x, y, x, y, out=scratch)          # degenerate window = point
    return [fid for fid in scratch if geometries[fid].contains_xy(x, y)]

That is the same funnel described in on-device geometry filtering, with the tree standing in for the bounding-box stage. Persisting the two arrays to a file and memory-mapping them on the next boot removes the build entirely from the start-up path, which on a device that reboots after every watchdog event is worth more than any query-side optimisation.

Validating a build before it is published

A packed tree is easy to build slightly wrongly and hard to notice, because a subtly broken tree returns some candidates and the exact predicate filters the rest — producing answers that are correct most of the time and quietly incomplete near the mistakes. Three assertions, run once after every build, catch essentially all of it.

Containment: every node’s envelope must contain all of its children’s envelopes. One pass over the levels, comparing each parent against its 16-entry run below. A violation means the union was computed over the wrong range, which is the most common indexing slip.

Completeness: querying the root’s own envelope must return exactly n leaf ids, each once. That single query exercises the whole traversal and catches an off-by-one in the child-count arithmetic that a spot check would miss.

Round trip: for a random sample of a few hundred features, query the feature’s own envelope and assert its id is in the result. That catches coordinate ordering mistakes — a swapped min and max — that the first two assertions pass over happily.

The three together cost a few hundred milliseconds on a 40 000-feature tree and run only at build time, which is exactly where the cost belongs.

Related in Spatial Indexing on Constrained Devices