Point-in-Polygon Geofencing for Streaming Sensor Events with Shapely

Geofencing a sensor stream means answering one question for every event as fast as it arrives: which monitoring zone contains this coordinate? The naïve approach — testing each point against every zone polygon — is O(zones) per event and collapses once a network has a few hundred zones and thousands of sensors reporting per minute. The fix is a prebuilt spatial index: a shapely.STRtree constructed once over your zone polygons answers each query in expected O(log zones) time by discarding every polygon whose bounding box cannot contain the point. This page builds a self-contained, restartable geofencing class that handles the two things that break naïve implementations in production — boundary points and NaN coordinates — and slots directly into the spatial windowing and geofencing stage of the pipeline.

Why a Prebuilt Index Changes the Complexity

A point-in-polygon test is not cheap. Shapely’s exact predicate walks the polygon’s edges and runs a ray-crossing or winding-number computation whose cost scales with vertex count — and administrative or watershed boundaries routinely carry hundreds of vertices. Running that test against every zone for every event multiplies two large numbers together.

An STRtree (Sort-Tile-Recursive tree) attacks the problem in two phases. First it packs every zone’s bounding box into a balanced R-tree at construction time. Then, for a query point, it descends the tree and returns only the zones whose bounding box contains the point — typically one or two candidates instead of the whole registry. The exact, expensive geometry test then runs only on those few survivors. The bounding-box filter is what makes streaming geofencing tractable: it converts the per-event cost from “test every polygon” to “test the handful that could plausibly match.”

Two properties make the tree a natural fit for streams. It is immutable once built, so it is safe to share read-only across worker threads without locking. And because zone registries change on the scale of weeks, the tree’s construction cost is paid once at startup and amortized across millions of events — never per batch.

The bounding-box filter also explains why zone shape matters for throughput. A long, thin, diagonal boundary — a river-following watershed edge, say — has a bounding box that covers a large rectangle of empty space, so the tree returns it as a candidate for many points that ultimately fail the exact test. Compact, convex zones produce tight boxes and fewer false candidates. You cannot always change the geometry, but you can simplify dense boundaries and split pathologically elongated zones so the coarse filter stays selective. This is the geofencing analogue of choosing a good partition key: the index only helps if its coarse phase actually eliminates most of the work.

One subtlety that trips up first implementations: an STRtree gives you a fast candidate set, not a fast answer. The exact predicate still runs on every survivor, so a query point sitting inside three overlapping zones triggers three full geometry tests. That is fine — three is not three hundred — but it means the tie-break logic that picks a single winner is not an afterthought. It runs on the hot path for every ambiguous point, which is why the implementation below resolves overlaps with a cheap min over integer-indexed ids rather than a second round of geometry.

Production-Ready Implementation

The class below loads and repairs zone polygons, builds the tree once, guards each incoming coordinate, and resolves overlaps deterministically so the same coordinate always maps to the same zone across restarts.

# python 3.11 · shapely==2.0.4 · geopandas==0.14.3 · pyproj==3.6.1
from __future__ import annotations
import math
from shapely import STRtree, Point, make_valid
from shapely.geometry.base import BaseGeometry


class StreamGeofencer:
    """Assign streaming sensor events to monitoring zones via a prebuilt STRtree.

    Build once at process start; call ``locate`` per event. The tree is
    immutable and safe to share read-only across worker threads.
    """

    def __init__(
        self,
        polygons: list[BaseGeometry],
        zone_ids: list[str],
        predicate: str = "covers",
    ) -> None:
        if len(polygons) != len(zone_ids):
            raise ValueError("polygons and zone_ids must be the same length")

        repaired: list[BaseGeometry] = []
        clean_ids: list[str] = []
        for geom, zid in zip(polygons, zone_ids):
            fixed = make_valid(geom)          # repair self-intersections
            if fixed.is_empty or not fixed.is_valid:
                continue                       # drop unrecoverable geometry
            repaired.append(fixed)
            clean_ids.append(zid)

        self._zone_ids = clean_ids
        self._geoms = repaired
        self._predicate = predicate
        self._tree = STRtree(repaired)         # O(n log n), built exactly once

    @staticmethod
    def _valid_coord(lat: float | None, lon: float | None) -> bool:
        """Reject None, NaN, Inf, out-of-range, and null-island coordinates."""
        if lat is None or lon is None:
            return False
        if not (math.isfinite(lat) and math.isfinite(lon)):
            return False
        if not (-90.0 <= lat <= 90.0 and -180.0 <= lon <= 180.0):
            return False
        return not (lat == 0.0 and lon == 0.0)

    def locate(self, lat: float, lon: float) -> str | None:
        """Return the zone_id containing (lat, lon), or None.

        On overlapping zones, returns the lexicographically smallest zone_id
        so the mapping is deterministic across restarts. Invalid coordinates
        return None without raising.
        """
        if not self._valid_coord(lat, lon):
            return None

        pt = Point(lon, lat)                   # Shapely order is (x=lon, y=lat)
        idx = self._tree.query(pt, predicate=self._predicate)
        if len(idx) == 0:
            return None
        return min(self._zone_ids[i] for i in idx)

Usage is a single construction followed by a hot per-event call:

from shapely.geometry import Polygon

zones = [
    Polygon([(-0.15, 51.50), (-0.05, 51.50), (-0.05, 51.55), (-0.15, 51.55)]),
    Polygon([(-0.05, 51.50), ( 0.05, 51.50), ( 0.05, 51.55), (-0.05, 51.55)]),
]
fencer = StreamGeofencer(zones, zone_ids=["central", "east"])

# Boundary point shared by both squares -> deterministic 'central'
print(fencer.locate(lat=51.52, lon=-0.05))   # "central"
print(fencer.locate(lat=51.52, lon=-0.10))   # "central"
print(fencer.locate(lat=float("nan"), lon=-0.10))  # None
print(fencer.locate(lat=51.52, lon=12.34))   # None (outside all zones)

The query call returns integer positions into the geometry list, which is Shapely 2.0 behaviour; the min over the matched zone_ids gives a stable tie-break for the boundary case where a point sits on the shared edge of two squares.

Parameter Tuning Guide

The predicate and how you pre-process zones matter more than any numeric knob. Pick per sensor deployment:

Sensor / zone type Predicate Zone pre-processing Tie-break
Air quality by district covers make_valid, dissolve slivers lowest zone_id
Water quality by watershed covers simplify to ~5 m tolerance lowest zone_id
Industrial fenceline exclusion covers buffer 0 to close gaps any match = alert
Overlapping regulatory + alert zones intersects keep all layers emit all matches
Point-source proximity (metres) buffered covers reproject to metric CRS nearest source

Two tuning notes. Simplifying dense boundaries with shapely.simplify to a few metres of tolerance cuts vertex count — and therefore per-test cost — with negligible effect on which zone a sensor falls in. And for proximity checks measured in metres, do not buffer in degrees; reproject the point and source to a metric CRS with pyproj first, because one degree of longitude is a different distance at every latitude.

Verification and Testing

Geofencing bugs are silent — a swapped axis or an excluded boundary produces valid-looking output that is simply wrong. Test the two failure modes explicitly.

import math
import pytest
from shapely.geometry import Polygon


def make_fencer():
    a = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
    b = Polygon([(10, 0), (20, 0), (20, 10), (10, 10)])
    return StreamGeofencer([a, b], zone_ids=["west", "east"])


def test_interior_point_matches_expected_zone():
    f = make_fencer()
    assert f.locate(lat=5, lon=5) == "west"
    assert f.locate(lat=5, lon=15) == "east"


def test_shared_boundary_is_deterministic():
    f = make_fencer()
    # x=10 is the shared edge; covers matches both, tie-break -> 'east'? no: min('east','west')='east'
    assert f.locate(lat=5, lon=10) == "east"
    # repeated calls never flip
    assert all(f.locate(5, 10) == "east" for _ in range(100))


def test_invalid_coordinates_return_none():
    f = make_fencer()
    assert f.locate(lat=math.nan, lon=5) is None
    assert f.locate(lat=math.inf, lon=5) is None
    assert f.locate(lat=0.0, lon=0.0) is None       # null island
    assert f.locate(lat=95.0, lon=5) is None        # out of range


def test_point_outside_all_zones_returns_none():
    f = make_fencer()
    assert f.locate(lat=5, lon=99) is None

The boundary test is the important one: covers matches both squares along x=10, and min("east", "west") resolves to "east" every time. Assert the exact winner so a future change to the tie-break rule is caught rather than shrugged off.

Gotchas

Latitude/longitude axis swap. Point(lat, lon) type-checks, runs, and returns confidently wrong zones. Shapely is (x, y) = (lon, lat). Build Point(lon, lat) and keep the assertion at the call site — this is the most common geofencing defect and it never raises on its own.

contains drops boundary points. With contains, a sensor sitting on the exact edge between two zones matches neither and reports as unassigned. Use covers so boundary points are always claimed, then break ties deterministically.

Rebuilding the tree per batch. The tree’s whole value is amortized construction. Rebuild only when the zone registry version changes; a per-batch rebuild reintroduces the O(z log z) cost on every window and erases the benefit over a plain loop.

Invalid polygons fail silently. A self-intersecting zone makes covers return arbitrary results for that zone with no exception. Run make_valid at load — as the constructor above does — and drop any geometry still invalid afterward instead of trusting it.

GPS jitter flips a stationary sensor between adjacent zones. A fixed sensor whose reported position wobbles a few metres around a shared boundary will match central one minute and east the next, producing phantom zone transitions that pollute per-zone counts. For sensors at known fixed locations, snap the event to the registered coordinate before the geofence query rather than trusting each transmitted fix. Reserve live coordinates for genuinely mobile platforms — drifting buoys, vehicle-mounted air monitors, wildlife trackers.

Metric buffers built in degrees are wrong everywhere but the equator. A “within 500 m of this outfall” check implemented by buffering the source point by 0.0045 degrees assumes a fixed degrees-to-metres ratio that only holds near the equator. At 55° latitude a degree of longitude is barely half its equatorial length, so the buffer becomes an east-west-squashed ellipse. Reproject both the point and the source to a local metric CRS with pyproj before any distance-based geofence.