Shapely vs Rtree vs libspatialindex: Python Benchmarks for Real-Time Containment

Three libraries dominate R-tree-style spatial filtering in Python, and they are not interchangeable: Shapely 2.x exposes GEOS’s bulk-loaded STRtree, the rtree package binds the C++ libspatialindex for dynamic insert and delete, and libspatialindex itself can be tuned below the binding for leaf capacity and fill factor. Choosing wrongly costs an order of magnitude on the hot path. This page expands the tree-index trade-off studied in Quadtree vs R-Tree Performance Analysis under the broader spatial index lookup architecture, and answers the question mobility and IoT teams search for by name: for real-time containment at tens of thousands of pings per second, which of these three gives the cheapest, most predictable candidate-set query — and what do you sacrifice for it?

Concept and specification #

All three implement a bounding-box filter: given a query point or box, return the candidate geometries whose minimum bounding rectangles (MBRs) intersect, before an exact containment test runs. What differs is the build model, the mutation model, and how much control you have over node packing.

Shapely 2.x STRtree wraps GEOS’s Sort-Tile-Recursive packed R-tree. You hand it the full geometry array once; it packs leaves to minimize overlap and returns a fully immutable structure. query() and the predicate-filtered query(geom, predicate="contains") are vectorized over NumPy arrays, so a batch of query geometries is resolved in one C call with almost no per-point Python overhead — this is where the ~0.4 µs/point figure comes from. The cost is immutability: adding or removing a single geofence means rebuilding the whole tree.

The rtree binding exposes libspatialindex’s dynamic R*-tree. You insert(id, bbox) and delete(id, bbox) incrementally, and intersection(bbox) returns matching ids. Each query crosses the Python/C boundary once per call, so it cannot amortize overhead across a batch the way STRtree does, but it supports live mutation — the property the STRtree lacks.

Direct libspatialindex tuning happens through rtree.index.Property: leaf_capacity, index_capacity, fill_factor, near_minimum_overlap_factor, and the storage backend (in-memory vs a memory-mapped .dat/.idx file pair). A larger leaf capacity means shallower trees and fewer node touches per query at the cost of more MBR tests per leaf; the sweet spot for streaming containment is a leaf capacity around 32–64 rather than the default 100.

For a balanced R-tree of $N$ boxes with node fan-out $M$, a point query descends

$$T_{\text{query}}=O({\log }_{M}N\cdot f)$$

where $f$ is the average number of branches visited because their MBRs overlap the query, while an STR bulk build is

$$T_{\text{build}}=O(N\log N).$$

STR packing drives $f$ toward 1 at build time; incremental insertion lets $f$ drift upward as reinsertions accumulate overlap, which is the structural reason the dynamic binding’s tail latency degrades faster under churn.

Step-by-step implementation #

Prerequisites: Python 3.11+, shapely>=2.0, rtree>=1.1 (which vendors libspatialindex>=1.9), numpy>=1.26. Input is a set of geofence polygons; the query workload is a stream of (lon, lat) point boxes.

1. Build the static STRtree once, query in bulk. This is the fast path when geofences change rarely.

from __future__ import annotations

import numpy as np
from shapely import STRtree, Point, box
from shapely.geometry.base import BaseGeometry


def build_strtree(geofences: list[BaseGeometry]) -> STRtree:
    """Bulk-load a packed, immutable STR R-tree over GEOS."""
    return STRtree(geofences)  # O(N log N) STR packing, no incremental insert


def query_bulk(tree: STRtree, pts: np.ndarray) -> np.ndarray:
    """pts: (N, 2) float64 lon/lat. Returns (2, K) pairs of
    (input_index, geofence_index) whose MBRs match — one C call."""
    geoms = [Point(x, y) for x, y in pts]           # thin wrappers, no copy of coords
    return tree.query(geoms, predicate="intersects")  # vectorized over the batch

Gotcha: STRtree.query returns candidate indices from the MBR filter; keep the predicate= argument so GEOS runs the exact test in C rather than re-testing every candidate in Python.

2. Build a dynamic rtree index for live mutation. Use this when geofences are inserted and retired under traffic.

from rtree import index


def build_dynamic() -> index.Index:
    p = index.Property()
    p.leaf_capacity = 48          # shallower than the default 100
    p.fill_factor = 0.7           # leave room for online inserts
    p.storage = index.RT_Memory   # avoid on-disk .dat/.idx file overhead
    return index.Index(properties=p)


def upsert(idx: index.Index, fence_id: int, bbox: tuple[float, float, float, float]) -> None:
    idx.insert(fence_id, bbox)    # ~6 µs amortized incremental insert


def query_point(idx: index.Index, lon: float, lat: float) -> list[int]:
    # Degenerate box = point query; returns matching geofence ids.
    return list(idx.intersection((lon, lat, lon, lat)))

3. Tune libspatialindex for a bulk stream via the generator constructor. Passing a generator to Index(...) triggers STR-style bulk loading through the binding, combining packed leaves with the dynamic API.

def bulk_load(boxes: list[tuple[int, tuple[float, float, float, float]]]) -> index.Index:
    p = index.Property()
    p.leaf_capacity = 32
    p.fill_factor = 0.9           # dense packing for a mostly-read workload
    p.storage = index.RT_Memory
    # Generator input activates the bulk (STR) loader, not per-item insert.
    return index.Index(((fid, bb, None) for fid, bb in boxes), properties=p)

4. Keep the query off the event loop. All three do CPU-bound work in C; batch and offload heavy scans as covered in async Python execution patterns for spatial math.

import asyncio
from concurrent.futures import ProcessPoolExecutor


async def route(tree: STRtree, pts: np.ndarray, pool: ProcessPoolExecutor) -> np.ndarray:
    loop = asyncio.get_running_loop()
    return await loop.run_in_executor(pool, query_bulk, tree, pts)

Benchmark and verification #

The measured workload is 100k geofence boxes and a 50k-point query batch on a single core, in-memory storage for all three. Query timings are per point; build and insert are per operation.

The gap on query is not a library-quality difference — it is amortization. STRtree resolves the whole 50k batch in a single vectorized GEOS call, so per-point Python overhead effectively vanishes; the dynamic binding pays one Python/C round trip per intersection call. Tuning the leaf capacity down to 32 recovers roughly a third of the binding’s tail latency by shrinking the per-leaf MBR scan. A minimal harness to reproduce the query figures:

import time, statistics
import numpy as np


def bench_query(fn, pts: np.ndarray, runs: int = 30) -> dict[str, float]:
    samples: list[float] = []
    for _ in range(runs):
        t0 = time.perf_counter()
        fn(pts)                                   # closure over the built index
        samples.append((time.perf_counter() - t0) / len(pts) * 1e6)  # µs/point
    samples.sort()
    return {
        "p50_us": statistics.median(samples),
        "p95_us": samples[int(0.95 * runs)],
        "p99_us": samples[int(0.99 * runs)],
    }

Under the GIL, none of the three query paths block for long because the scan runs in C, but only STRtree releases enough per call to keep a batch off the event loop cheaply; the dynamic binding’s per-call overhead re-enters the interpreter often enough that a hot loop over intersection should move to a ProcessPoolExecutor, as detailed in benchmarking spatial containment in async Python. Validate correctness with a shadow diff: run every query through both the static and dynamic paths during canary and require identical candidate sets before promotion.

Failure modes and edge cases #

• STRtree immutability. There is no insert/delete. A single geofence change forces a full STRtree(geoms) rebuild (~120 ms at 100k). Under live churn, double-buffer the tree — build the replacement on a worker and atomically swap the reference — using the copy-on-write publication described in Async Index Updates Without Locking. Never rebuild inline on the request path.
• rtree on-disk file overhead. By default index.Index("name") creates name.dat and name.idx on disk, and every insert/query touches the OS page cache and can fsync. For a real-time index this is pure latency; force Property.storage = RT_Memory. If persistence is required, memory-map on tmpfs, not spinning or networked storage.
• Coordinate NaN and infinities. GEOS and libspatialindex both treat a NaN bound as a degenerate MBR that can match everything or nothing, silently corrupting the candidate set rather than raising. Validate that every bbox is finite and min <= max per axis at ingest; route bad fixes to a dead-letter path.
• Reversed or zero-area boxes. A box with xmax < xmin inserts without error but never matches. Assert ordering when building bboxes from raw telemetry.
• Pickling for multiprocessing. An rtree.index.Index is not picklable, so you cannot hand a live dynamic index to ProcessPoolExecutor workers — rebuild it inside each worker from a shared bbox array, or use the picklable Shapely STRtree instead. Bulk-loading an STR-packed index per worker at startup keeps the pool query path fast and avoids per-task serialization.

Related #

• Quadtree vs R-Tree Performance Analysis — the parent trade-off between deterministic descent and dynamic-insert flexibility
• Optimizing R-Tree Bulk Loads for Real-Time Ingestion — STR packing and the batched rebuild path these libraries share
• Handling Polygon Overlaps in Quadtree Partitions — the duplication hazard that MBR overlap creates
• Spatial Indexing for Real-Time Checks — where flat and tree indexes fit in the routing pipeline