Benchmarking Spatial Containment in Async Python

Real-time telemetry pipelines for mobility, logistics, and IoT fleets routinely resolve coordinate streams at millions of points per second, and at the centre of dispatch routing, compliance enforcement, and dynamic geofencing sits one primitive: deciding whether a moving asset lies inside a defined polygon. When that primitive runs in Python, the collision between synchronous C-backed geometry engines and the cooperative asyncio event loop produces a specific, reproducible failure — p99 latency cliffs, loop starvation, and uncontrolled heap growth under moderate concurrency. This page is the measurement counterpart to async Python execution patterns for spatial math: the parent overview argues which offload boundary to choose, and this deep-dive shows how to benchmark each boundary so the choice is grounded in numbers rather than folklore. It sits within the wider latency contract defined in Core Architecture & Latency Constraints, where every stage of the pipeline is allocated a slice of a fixed per-event budget.

The reason a benchmark is non-optional: a naive polygon.contains(point) call benchmarks at sub-microsecond cost in isolation and then misses its SLA by two orders of magnitude in production. The gap is not the geometry — it is the event loop the geometry is blocking. Until you measure compute-bound stall separately from I/O wait, every tuning decision is a guess.

Concept & Specification #

Spatial containment is a point-in-polygon evaluation: given a test point and an ordered vertex ring, decide containment via ray casting or the winding-number rule. Both kernels are $O(V)$ in vertex count $V$, but their cost in an async runtime is dominated not by $V$ but by where they execute relative to the loop.

The quantity a benchmark must isolate is the blocking interval — the wall-clock time the kernel holds the Global Interpreter Lock and prevents the selector from polling sockets. The effective per-event latency observed by a client is

$$L_{\text{event}}=L_{\text{queue}}+L_{\text{compute}}+L_{\text{io}}$$

where $L_{\text{queue}}$ grows without bound the moment $L_{\text{compute}}$ runs inline, because every coroutine waiting behind the blocking call inherits the full kernel duration. The goal of the benchmark is to attribute the measured tail to one of these three terms. For an offloaded design, the amortized compute cost over a batch of size $B$ is

$$C_{\text{event}}=C_{\text{kernel}}+\frac{C_{\text{ipc}}}{B}$$

so the inter-process transfer term decays as $O(1/B)$ while the kernel term stays flat — the relationship the batch-size sweep below is designed to expose.

The parameters that materially move the numbers:

Step-by-Step Implementation #

Prerequisites: Python 3.11+, shapely>=2.0 (vectorized shapely.contains), numpy>=1.24, pytest-benchmark, pytest-asyncio. Coordinate batches are float64 arrays of shape (N, 2); polygons are pre-loaded into a read-only catalog so geometry never crosses a process boundary per call.

1. Pin a deterministic baseline kernel. Use a prepared geometry so repeated topology validation is not silently counted against the wrong term.

from __future__ import annotations

import numpy as np
import shapely
from numpy.typing import NDArray
from shapely import Polygon
from shapely.prepared import PreparedGeometry

def build_catalog(rings: dict[int, NDArray[np.float64]]) -> dict[int, PreparedGeometry]:
    """Prepare each polygon once; prepared geometry caches the GEOS index."""
    catalog: dict[int, PreparedGeometry] = {}
    for zone_id, ring in rings.items():
        catalog[zone_id] = shapely.prepare(Polygon(ring)) or Polygon(ring)
    return catalog

def evaluate_batch(poly: Polygon, coords: NDArray[np.float64]) -> NDArray[np.bool_]:
    """Vectorized containment; one GEOS call for the whole batch."""
    points = shapely.points(coords[:, 0], coords[:, 1])
    return shapely.contains(poly, points)

Gotcha: shapely.prepare mutates in place and returns None; never write geom = shapely.prepare(geom). The or fallback above keeps the reference valid.

2. Measure the blocking interval, not just throughput. Wrap the kernel so the loop-lag it induces is observable. Run the identical workload three ways — inline, via asyncio.to_thread, and via loop.run_in_executor against a ProcessPoolExecutor — to quantify thread-pool overhead against GIL contention.

import asyncio
from concurrent.futures import ProcessPoolExecutor
from time import perf_counter_ns

async def run_inline(poly: Polygon, coords: NDArray[np.float64]) -> float:
    t0 = perf_counter_ns()
    evaluate_batch(poly, coords)              # holds the GIL for the full kernel
    return (perf_counter_ns() - t0) / 1e6     # ms

async def run_threaded(poly: Polygon, coords: NDArray[np.float64]) -> float:
    t0 = perf_counter_ns()
    await asyncio.to_thread(evaluate_batch, poly, coords)
    return (perf_counter_ns() - t0) / 1e6

async def run_process(
    pool: ProcessPoolExecutor, zone_id: int, coords: NDArray[np.float64]
) -> float:
    loop = asyncio.get_running_loop()
    t0 = perf_counter_ns()
    await loop.run_in_executor(pool, _worker_eval, zone_id, coords)
    return (perf_counter_ns() - t0) / 1e6

3. Probe loop responsiveness concurrently. A heartbeat coroutine measures the gap between scheduled and actual wake-ups; that gap is $L_{\text{queue}}$.

async def loop_lag_probe(samples: list[float], stop: asyncio.Event) -> None:
    """Each 5ms tick that arrives late reveals event-loop starvation."""
    interval = 0.005
    while not stop.is_set():
        scheduled = perf_counter_ns()
        await asyncio.sleep(interval)
        actual = (perf_counter_ns() - scheduled) / 1e6
        samples.append(actual - interval * 1e3)  # lag in ms beyond the sleep

4. Sweep batch size to amortize IPC. Drive the process-pool path at $B\in 1,16,64,128,256$ and record the per-event cost; the curve should flatten as $C_{\text{ipc}}/B$ collapses, confirming the model above.

5. Track allocation churn with tracemalloc. Snapshot every 10k evaluations and diff with compare_to to attribute heap growth to per-event Point objects versus the static catalog — the difference between a stable pipeline and one that triggers gen-2 GC pauses mid-evaluation.

Benchmark / Verification #

Figures below are a 4-core worker, CPython 3.11, a 200-vertex municipal polygon, NumPy-backed coordinate arrays, sustained 50k events/sec offered load, perf_counter_ns timing, prepared geometry, gen-2 GC frozen after warm-up.

The before/after story is the first and last compute rows: moving from inline evaluation to a batched process pool drops P99 from 180ms to 6.2ms and lifts sustained throughput roughly 12×, while the loop-lag probe confirms the event loop is no longer starved (174ms → 0.8ms). The starred lesson from the batch sweep: per-event process dispatch (B=1) is worse than threads on small payloads because pickling a coordinate plus geometry costs 15–40µs against a 120µs kernel; the win appears only once C_{ipc} is amortized across ≥64 events. The AABB row exists to justify the degraded-mode fallback — it is GIL-safe, cheap enough to run on the loop, and rejects ~92% of candidates before any GEOS call.

A regression gate worth wiring into CI: assert to_thread P95 stays under 3ms and batched-pool P99 under 8ms at the 50k/s fixture. Both thresholds catch the two real regressions — an accidental inline call, and a dropped batch buffer that reverts to per-event dispatch.

Failure Modes & Edge Cases #

Degenerate geometry. Self-intersecting rings make ray casting’s even-odd rule disagree with winding number on overlap regions; validate with shapely.is_valid and repair via shapely.make_valid before preparing, or the benchmark measures a kernel that ships wrong answers fast. Empty coordinate batches must short-circuit before dispatch — handing a zero-length array to a process pool still pays full IPC for no work.

NaN and out-of-range coordinates. A NaN latitude propagates through shapely.contains as False silently, so a corrupt GPS sample looks like a legitimate “outside” result. Filter with np.isfinite(coords).all(axis=1) at ingest and count the rejects as a metric; a rising reject rate is a sensor fault, not a containment miss.

GIL contention above core count. asyncio.to_thread recovers throughput only while the C kernel releases the GIL; the Python glue around it does not. Sizing the thread pool past the physical core count produces thread thrashing that raises P99 — the sweep will show throughput plateau then regress, the signature of oversubscription.

GC pressure masquerading as compute cost. The millisecond P99 outliers at high load are usually a gen-2 collection landing on top of an evaluation, not edge-scan cost. Confirm by correlating gc.get_stats() collection counts with the spike timestamps; if they align, call gc.freeze() after warm-up and keep the hot path allocation-free using the vectorized array API rather than per-point Point construction. This is the same heap discipline detailed in memory-constrained spatial processing.

Surge overload. When offered load exceeds the measured throughput ceiling, an unbounded executor queue converts a compute stall into a memory exhaustion crash. Bound the dispatch queue and trip a circuit breaker at 2× expected depth, routing overflow to the AABB fast-path fallback so the pipeline degrades to coarse accuracy rather than failing entirely.

Related #

• Async Python execution patterns for spatial math — the parent overview: which offload boundary to choose and why, with the dispatcher critical path these benchmarks measure.
• Point-in-polygon algorithm benchmarks — the kernel-level ray-casting vs winding-number figures that feed the compute term here.
• Memory-constrained spatial processing — the allocation-free polygon layout that keeps GC pauses out of these P99 numbers.
• Up one level: Core Architecture & Latency Constraints — the per-event latency budget every containment check must fit inside.