Source: https://gometry.monicz.dev/guide/indexing/

# Spatial indexing & joins

A spatial index narrows a query to features whose bounding boxes could match,
then applies the requested geometry predicate.
[`gm.SpatialIndex`](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex "SpatialIndex") builds an in-memory
[R-tree](https://en.wikipedia.org/wiki/R-tree) over geometry envelopes. It exposes
[`crs`](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.crs "crs property") and [`epoch`](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.epoch "epoch property") 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)`](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.candidates "candidates") returns bounding-box
matches. [`idx.query(q, predicate=...)`](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.query "query") refines
them with an exact geometry predicate. Both return source-row handles in
ascending order; removing rows does not renumber the remaining handles.

```python
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())
```

```text
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.

The envelopes of rows 0, 1, and 2 overlap the query envelope, so all three become candidates. Exact intersects keeps rows 0 and 2 and rejects the diagonal triangle in row 1; row 3 never becomes a candidate.

- 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.

```python
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]
```

```text
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.

```python
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())
```

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

The full vocabulary is the same as the [top-level predicates](https://gometry.monicz.dev/guide/predicates/)
— [`intersects`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.intersects "intersects") (the default), [`contains`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.contains "contains"), [`contains_properly`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.contains_properly "contains_properly"), [`covers`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.covers "covers"),
[`within`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.within "within"), [`covered_by`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.covered_by "covered_by"), [`equals`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.equals "equals"), [`touches`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.touches "touches"), [`crosses`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.crosses "crosses"), [`overlaps`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.overlaps "overlaps"), plus the
metric [`dwithin`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.dwithin "dwithin"). Only [`disjoint`](https://gometry.monicz.dev/api/toplevel/predicates/#gometry.disjoint "disjoint") is rejected: it matches everything far
away, which no envelope index can accelerate.

For distance predicates, pass `distance=` and `unit=`:

```python
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())
```

```text
within 1 deg: [0, 1]
```

Distance `unit=` follows the CRS metric rules in [CRS, units & measurement](https://gometry.monicz.dev/guide/crs/).
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](https://gometry.monicz.dev/guide/crs/#geographic-point-to-segment-accuracy) 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`](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.nearest "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).

```python
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())
```

```text
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`](https://gometry.monicz.dev/api/geometry/coordinates/#gometry.Groups "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:

```python
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}")
```

```text
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=...)`](https://gometry.monicz.dev/api/toplevel/index-join/#gometry.nearest "nearest") 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:

```python
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}")
```

```text
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`](https://gometry.monicz.dev/api/toplevel/index-join/#gometry.join "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.

```python
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}")
```

```text
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`](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.insert "insert") appends a geometry in the index's frame,
and [`remove`](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.remove "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.

```python
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())
```

```text
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](https://gometry.monicz.dev/ecosystem/text-formats/) and construct a new index.

## See also

- [Predicates](https://gometry.monicz.dev/guide/predicates/) — the exact relationships [`query`](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.query "query") refines with.
- [Discrete grids](https://gometry.monicz.dev/guide/grids/) — cell coverages as another candidate layer.
- [Arrays & performance](https://gometry.monicz.dev/guide/arrays/) — array query → [`Groups`](https://gometry.monicz.dev/api/geometry/coordinates/#gometry.Groups "Groups"), vectorized refine.
- [API: SpatialIndex](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex "SpatialIndex") ·
  [SpatialIndex.candidates](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.candidates "candidates") ·
  [SpatialIndex.query](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.query "query") ·
  [SpatialIndex.nearest](https://gometry.monicz.dev/api/spatial-index/#gometry.SpatialIndex.nearest "nearest") ·
  [join](https://gometry.monicz.dev/api/toplevel/index-join/#gometry.join "join")