Spatial indexing & joins¶

A spatial index narrows a query to features whose bounding boxes could match, then applies the requested geometry predicate. gm.SpatialIndex builds an in-memory R-tree over geometry envelopes. It exposes crs and epoch for the shared frame of its rows.

query, candidates, nearest, and join accept a GeometryArray or an iterable of Geometry | None rows. Missing query rows remain None groups and keep their original row positions.

Candidates vs query¶

idx.candidates(q) returns bounding-box matches. idx.query(q, predicate=...) refines them with an exact geometry predicate. Both return source-row handles in ascending order; removing rows does not renumber the remaining handles.

import gometry as gm

query = gm.box(0.0, 0.0, 1.0, 1.0)
geoms = gm.GeometryArray([
    gm.box(0.2, 0.2, 0.8, 0.8),
    gm.Polygon([(0.5, 2.0), (2.0, 2.0), (2.0, 0.5)]),
    gm.box(0.9, 0.9, 1.6, 1.6),
    gm.box(3.0, 3.0, 4.0, 4.0),
])
idx = gm.SpatialIndex(geoms)

cands = idx.candidates(query)                     # bbox prefilter — candidate superset
hits = idx.query(query, predicate="intersects")   # exact refine — matching rows

assert cands.tolist() == [0, 1, 2]
assert hits.tolist() == [0, 2]

print("candidates (bbox):", cands.tolist())
print("exact intersects:", hits.tolist())

candidates (bbox): [0, 1, 2]
exact intersects: [0, 2]

In this example, geometry 1 is a false positive because its diagonal triangle's bounding box overlaps the query while the shape does not.

Array queries return one Groups row per query row. A missing query appears as None and is marked in Groups.is_missing; a present query with no matches appears as an empty row. Scalar queries keep their flat int64 result.

The prefilter admits rows 0, 1, 2; refinement keeps rows 0 and 2.

query
candidate envelopes
exact matches
bbox false positive
never a candidate

candidates supplies rows for caller-owned refinement or a coarse spatial key. query performs the exact predicate refinement. A positive-infinite candidate expansion admits every live non-empty indexed row for a query with bounds.

Periodic geographic bounds¶

For a geographic index, the candidate envelope uses the CRS-native longitude period and pole. Exact refinement then uses the same frame-aware topology as a direct predicate. EPSG:4807 is a public grad-unit example: one turn is 400 grads, and its north pole is 100 grads.

import gometry as gm

# EPSG:4807 uses grads, so the north pole is 100.
pole = 100.0
source = gm.Point(0.0, pole, crs=4807)
alias = gm.Point(200.0, pole, crs=4807)
index = gm.SpatialIndex([source])

print("periodic candidates:", index.candidates(alias).tolist())
print("exact intersects:", index.query(alias, predicate="intersects").tolist())
assert index.candidates(alias).tolist() == [0]
assert index.query(alias, predicate="intersects").tolist() == [0]

periodic candidates: [0]
exact intersects: [0]

The bounding-box search reaches the pole's alternate longitude spelling, and the exact predicate confirms that both spellings identify the same geographic point.

Predicates¶

query accepts any of the standard spatial predicates via predicate=. The default is "intersects".

The operand order is always predicate(query_geom, indexed_row). For example, SpatialIndex(polygons).query(point, predicate="within") finds polygons where gm.within(point, polygon) is true; predicate="contains" would ask whether the point contains each polygon.

import gometry as gm

pts = gm.points([21.0, 30.0, 21.5], [52.0, 52.0, 52.1], crs=4326)
idx = gm.SpatialIndex(pts)
box = gm.box(20.0, 51.0, 22.0, 53.0, crs=4326)

for pred in ("intersects", "contains", "contains_properly", "covers"):
    print(f"{pred:>17}:", idx.query(box, predicate=pred).tolist())

       intersects: [0, 2]
         contains: [0, 2]
contains_properly: [0, 2]
           covers: [0, 2]

The full vocabulary is the same as the top-level predicates — intersects (the default), contains, contains_properly, covers, within, covered_by, equals, touches, crosses, overlaps, plus the metric dwithin. Only disjoint is rejected: it matches everything far away, which no envelope index can accelerate.

For distance predicates, pass distance= and unit=:

import gometry as gm

pts = gm.points([21.0, 21.5, 30.0], [52.0, 52.1, 52.0], crs=4326)
idx = gm.SpatialIndex(pts)
probe = gm.Point(21.0, 52.0, crs=4326)

print("within 1 deg:", idx.query(probe, predicate="dwithin", distance=1.0, unit="planar").tolist())

within 1 deg: [0, 1]

Distance unit= follows the CRS metric rules in CRS, units & measurement. The unit='meters' override forces SI, and unit='planar' uses raw coordinate math. dwithin accepts a positive-infinite threshold and matches only non-empty geometries. On a geographic CRS, dwithin raises CRSError if the accumulated evaluation work limit cannot certify whether a present pair is within the radius; it does not report an unresolved pair as False. The CRS guide describes its numerical accuracy and failure behavior. Successful array results keep one boolean per input pair, and missing rows and present empty geometries remain in their positions as False.

Nearest neighbors¶

idx.nearest finds the k closest indexed geometries to a query geometry and returns their handles (nearest first). Distance uses the same unit= rule as query / free dwithin (CRS-natural by default).

import gometry as gm

stations = gm.points(
    [21.0, 21.05, 22.0, 19.0],
    [52.0, 52.02, 53.0, 51.0],
    crs=4326,
)
idx = gm.SpatialIndex(stations)
me = gm.Point(21.01, 52.0, crs=4326)

ids = idx.nearest(me, k=2, unit="meters")
print("nearest 2 indices (closest first):", ids.tolist())

nearest 2 indices (closest first): [0, 1]

For geographic point pairs, nearest queries use ellipsoidal geodesic distance. For linework and filled geometries, they minimize that distance over straight chart-edge support, using a sound geodesic lower bound to prune candidates. Ordering follows the refined distances; the number of evaluations depends on the data and query.

max_distance= caps the search radius.
return_distance=True returns (indices, distances) for a scalar query or (matches, distances) for an array query; scalar indices is a read-only int64 NumPy array, array matches is CSR Groups, and distances is a read-only float64 NumPy array aligned with the ids.
For array queries, missing input rows are None groups; present queries with no indexed match are empty groups.
exclusive=True skips candidates structurally equal to the query and returns the nearest other features:

import gometry as gm

stations = gm.points([21.0, 21.05, 22.0], [52.0, 52.02, 53.0], crs=4326)
idx = gm.SpatialIndex(stations)

for i, station in enumerate(list(stations)):
    nearest_other = idx.nearest(station, exclusive=True)[0]
    print(f"station {i} → nearest other station {nearest_other}")

station 0 → nearest other station 1
station 1 → nearest other station 0
station 2 → nearest other station 1

The top-level function gm.nearest(geometry, values, k=..., unit=...) builds an index for one-shot queries.

Measure candidate selectivity¶

Compare the number of bounding-box candidates with the exact result for the same query and predicate:

import gometry as gm

q = gm.box(0.0, 0.0, 1.0, 1.0)
geoms = gm.GeometryArray([
    gm.box(0.2, 0.2, 0.8, 0.8),
    gm.Polygon([(0.5, 2.0), (2.0, 2.0), (2.0, 0.5)]),
    gm.box(0.9, 0.9, 1.6, 1.6),
    gm.box(3.0, 3.0, 4.0, 4.0),
])
idx = gm.SpatialIndex(geoms)

candidate_count = len(idx.candidates(q))
refined_count = len(idx.query(q, predicate="intersects"))
print(f"candidates: {candidate_count} → refined intersects: {refined_count}")

candidates: 3 → refined intersects: 2

A large candidate count relative to the refined count means the envelopes are not selective; large sparse multipolygons are a common cause.

Joins: many-to-many done right¶

gm.join returns every matching pair across two collections. It runs the bbox prefilter and exact refine in Rust and returns a (left, right) pair of read-only int64 NumPy columns. Zipping the columns produces row pairs.

import gometry as gm

points = gm.points([21.0, 30.0, 21.5], [52.0, 52.0, 52.1], crs=4326)
areas = gm.GeometryArray([
    gm.box(20.0, 51.0, 22.0, 53.0, crs=4326),
    gm.box(29.0, 51.0, 31.0, 53.0, crs=4326),
])

left, right = gm.join(points, areas, predicate="within")
print("matched (point, area) pairs:", list(zip(left.tolist(), right.tolist(), strict=True)))
for li, ri in zip(left, right, strict=True):
    print(f"point {li} within area {ri}")

matched (point, area) pairs: [(0, 0), (1, 1), (2, 0)]
point 0 within area 0
point 1 within area 1
point 2 within area 0

join refines to exact predicates. Use distance=/unit= for a predicate="dwithin" join. For a direct index query, compare len(idx.candidates(q)) with len(idx.query(q, predicate=...)) to measure the candidate superset and its exact refinement. Geographic dwithin refinement raises CRSError if the point-edge search exhausts its budget without proving the threshold result.

Many-to-many query cost

The brute-force for a in polys: for b in points: if gm.contains(a, b) pattern is O(N·M). gm.join builds and queries an index in Rust, then refines candidates exactly. Its cost is driven by bounding-box selectivity and candidate count; worst-case overlap can still approach O(N·M). Vectorized predicates reject mismatched non-scalar lengths; join returns the many-to-many result.

Insert and remove¶

An index is mutable, and its integer handles identify live rows. Lookup and removal accept non-negative integer-index handles and reject bool. idx.get(handle, default) returns default for an unknown or removed handle; negative and non-index inputs raise. Membership checks return False for non-index probes.

insert appends a geometry in the index's frame, and remove removes a row by handle. Queries include inserted rows and exclude removed ones without renumbering the remaining handles.

Batch insertion rejects empty geometries and GeometryArrays containing missing rows. It validates the complete input and frame before changing the index, so input-validation errors leave its rows and frame unchanged.

import gometry as gm

pts = gm.points([21.0, 30.0], [52.0, 52.0], crs=4326)
idx = gm.SpatialIndex(pts)
print("bulk rows:", len(idx))
handle = idx.insert(gm.Point(21.5, 52.1, crs=4326))
print("after insert:", len(idx), "| handle:", handle)
print("nearest:", idx.nearest(gm.Point(21.0, 52.0, crs=4326), k=2, unit="planar").tolist())
idx.remove(handle)
print("after remove handle", handle, "→ nearest still works:",
      idx.nearest(gm.Point(21.0, 52.0, crs=4326), k=2, unit="planar").tolist())

bulk rows: 2
after insert: 3 | handle: 2
nearest: [0, 2]
after remove handle 2 → nearest still works: [0, 1]

Equality and serialization¶

SpatialIndex compares like a mapping: equality uses its live handle-to-geometry contents, including after removals. It is mutable and unhashable.

Pickling preserves live handles and prior insertions and removals, rebuilding the search structure on load. For portable storage, keep the source geometries as WKB or EWKB and construct a new index.