RouteMatching

Building a Geofence Index for Fleet-Scale Point-in-Polygon Tests

Testing whether a fix is inside a geofence is trivial. Testing four million fixes against twelve thousand geofences is forty-eight billion tests, and a naive loop over polygon.contains(point) will still be running tomorrow.

The fix is the standard two-stage spatial query: a cheap bounding-box filter that eliminates almost everything, then an exact test on the handful that survive. Both stages have implementation details that matter at fleet scale, and this page covers them alongside the assertions that keep the index trustworthy as the fence set changes. It supports the matching work in location typing and POI matching for stops.

Compatibility and Configuration Requirements

Requirement Value Notes
shapely ≥ 2.0 STRtree with vectorised query, and prepared geometries
CRS One projected CRS for fences and fixes Mixing degrees and metres silently returns nothing
Fence validity All polygons valid An invalid polygon gives undefined containment results
Memory ~200 bytes per fence vertex A twelve-thousand-fence set is tens of megabytes
Rebuild cadence Whenever the fence set changes The tree is immutable once built

The CRS requirement is the one that bites. A fence set in EPSG:4326 queried with projected fixes produces zero matches and no error, because the coordinate ranges do not overlap. Asserting that both sides carry the same CRS before building the tree costs one line.

The Index

from __future__ import annotations

import numpy as np
import shapely
from shapely import STRtree
from shapely.geometry import Polygon


class GeofenceIndex:
    """Two-stage point-in-polygon over a fixed set of fences.

    Stage one is an R-tree over bounding boxes; stage two is an exact
    containment test against prepared geometries. Both fences and query
    points must be in the same projected CRS.
    """

    def __init__(self, fences: list[Polygon], fence_ids: list[str], crs_epsg: int):
        invalid = [i for i, g in enumerate(fences) if not g.is_valid]
        if invalid:
            raise ValueError(f"{len(invalid)} invalid fence geometries, e.g. index {invalid[0]}")
        self.crs_epsg = crs_epsg
        self.ids = np.asarray(fence_ids)
        self.geoms = np.asarray(fences, dtype=object)
        shapely.prepare(self.geoms)          # cache edge structure for fast contains
        self.tree = STRtree(self.geoms)

    def query(self, xs: np.ndarray, ys: np.ndarray) -> list[np.ndarray]:
        """Return the fence ids containing each point, as one array per point."""
        pts = shapely.points(xs, ys)
        # Stage one: bounding-box candidates for every point at once.
        pt_idx, fence_idx = self.tree.query(pts, predicate="intersects")
        # Stage two: exact test, only on the surviving pairs.
        hit = shapely.contains_xy(self.geoms[fence_idx], xs[pt_idx], ys[pt_idx])
        pt_idx, fence_idx = pt_idx[hit], fence_idx[hit]

        out = [np.empty(0, dtype=self.ids.dtype) for _ in range(len(xs))]
        for p in np.unique(pt_idx):
            out[p] = self.ids[fence_idx[pt_idx == p]]
        return out

Three choices do the work. shapely.prepare caches each polygon’s edge structure so repeated containment tests skip the setup cost. tree.query is called with the whole point array, not per point, which keeps the traversal in compiled code. And the exact test runs only on the surviving candidate pairs, which after the bounding-box stage is typically one or two per point rather than twelve thousand.

What each stage eliminates For one fix tested against twelve thousand fences, the R-tree bounding-box stage reduces the candidate set to three, and the exact containment test reduces it to one. The expensive test therefore runs three times per fix rather than twelve thousand. Candidates per fix, 12 000-fence set all fences 12 000 after R-tree 3 — bounding boxes overlap after exact test 1 — actually inside The R-tree removes 99.97 % of the work, and it removes it before any geometry is examined. The two survivors that fail the exact test are fences whose bounding box overlaps but whose shape does not. Long thin fences — a road corridor, a rail siding — have poor bounding boxes and inflate stage two.

Batching and Throughput

Calling query per point defeats the vectorisation. The tree traversal is compiled code, and the cost of crossing from Python into it once per fix dominates everything else at fleet scale.

BATCH = 500_000

def label_fixes(index: GeofenceIndex, fixes) -> list[np.ndarray]:
    out = []
    for start in range(0, len(fixes), BATCH):
        chunk = fixes[start : start + BATCH]
        out.extend(index.query(chunk["x"].to_numpy(), chunk["y"].to_numpy()))
    return out

Half a million points per call is a reasonable batch: large enough that the per-call overhead disappears, small enough that the intermediate candidate arrays stay comfortably in memory. The candidate arrays are the memory peak, and they are roughly the number of points times the average bounding-box overlap — a few million entries, not a few hundred million.

Where the fixes are already partitioned by vehicle-day, batching by partition rather than by a fixed count is simpler and gives similar throughput.

Batch size against throughput Points tested per second against batch size. At one point per call throughput is around eleven thousand per second, dominated by Python call overhead. It rises steeply to about two hundred thousand at ten thousand points per batch and is flat beyond a hundred thousand. Points tested per second, 12 000-fence index 240 k/s 120 k/s 0 flat — batch bigger buys nothing 1 1 k 100 k 5 M points per query call (log scale) A twenty-fold difference between the naive loop and a batched call, on identical geometry. Past a hundred thousand the only thing a larger batch changes is the peak memory of the candidate arrays.

Nesting and Priority

Fleet geofences nest by design: a loading bay inside a yard inside a site, or a customer zone inside a city zone. A fix inside the bay is inside all three, and the index correctly returns all three.

Collapsing to one at test time is the mistake. The hierarchy is information — a dwell in the bay and a dwell in the yard mean different things operationally — and it cannot be recovered once discarded. Return every hit and apply a priority rule downstream, where the consumer’s needs are known.

def most_specific(hits: np.ndarray, area_by_id: dict[str, float]) -> str | None:
    """Smallest containing fence — the usual 'where is it really' rule."""
    return min(hits, key=lambda fid: area_by_id[fid], default=None)

Smallest-area is the right default because a nested fence is almost always more specific than its parent. Where fences overlap without nesting — two adjacent customer zones with a shared access road — smallest-area is arbitrary, and an explicit priority column on the fence set is the only honest answer.

Store all the hits, not just the winner. A stop attributed to a bay is also in the site, and site-level reporting should not have to re-run the spatial test to find that out.

One fix, three correct answers A fix inside a loading bay, which is inside a yard, which is inside a customer site. All three fences contain it, and all three are useful: the bay for service-time analysis, the yard for access analysis, the site for customer reporting. customer site yard bay 3 all three returned bay 3 — service time yard — access and queueing site — customer reporting smallest area wins where one is needed Collapsing to one hit at test time discards a hierarchy that cannot be recovered later. Where fences overlap without nesting, smallest-area is arbitrary — set an explicit priority. Store the full hit list; site-level reporting should never need a second spatial pass.

Keeping the Index Honest

The index is rebuilt whenever the fence set changes, and three assertions catch the failures that matter.

Every fence is valid. An invalid polygon — self-intersecting, usually from a hand-drawn boundary — gives undefined containment results rather than an error. The constructor above rejects them, which is the right place because a build-time failure is far cheaper than silently wrong labels.

The CRS matches on both sides. Assert it rather than trusting it. The symptom of a mismatch is zero matches everywhere, which reads as “no vehicles visited any site” and is often blamed on the fixes.

Hit rate is stable across rebuilds. Compare the share of fixes matching at least one fence against the previous index. A step change means fences were added, removed or moved, and knowing which is much easier at rebuild time than a week later.

def assert_hit_rate_stable(new_rate: float, previous_rate: float, tol: float = 0.02) -> None:
    if abs(new_rate - previous_rate) > tol:
        raise SystemExit(
            f"geofence hit rate moved from {previous_rate:.3f} to {new_rate:.3f} — "
            "review the fence set change before promoting this index"
        )

Test that assertion by rebuilding with a fence deliberately removed; it fires immediately. A check nobody has seen fail is not evidence that the index is correct.

Common Pitfalls Specific to This Technique

Building the tree per batch. The R-tree construction is the expensive part and the fence set does not change between batches. Build once, query many times.

Testing in degrees. A bounding-box query in EPSG:4326 works, but any buffer or distance applied to a fence is then in degrees, and the resulting fence is a different size at each latitude.

Forgetting shapely.prepare on large polygons. Site boundaries traced from imagery routinely have hundreds of vertices, and the unprepared exact test on those is several times slower.

Assuming one fence per fix. Nested fences are normal, and code that takes hits[0] gets whichever one the tree happened to return first, which is not stable across rebuilds.

Related

Up to section Location Typing & POI Matching for Stops Related Matching GPS Stops to Commercial POI Databases in Python