Spatial Query And Join Assembly

Request Signals

  • sindex query

  • sjoin

  • dwithin

  • nearest

  • pandas semantics

  • join assembly

  • spatial index

  • cccl

Open First

  • docs/architecture/spatial-joins.md

  • src/vibespatial/spatial/query.py

  • src/vibespatial/api/sindex.py

  • src/vibespatial/api/tools/sjoin.py

  • docs/architecture/binary-predicates.md

Verify

  • uv run pytest tests/test_spatial_query.py

  • uv run pytest tests/upstream/geopandas/tests/test_sindex.py -q

  • uv run pytest tests/upstream/geopandas/tools/tests/test_sjoin.py -k "predicate or nearest"

  • uv run vsbench run spatial-query --rows 20000 --arg overlap_ratio=0.2

  • uv run python scripts/profile_fixture_spatial_query.py --fixture polygons-regular-grid-rows100000 --operation query --ensure

  • uv run python scripts/check_docs.py --check

Risks

  • Query/index speedups are meaningless if join assembly breaks index names, suffixes, or geometry-column retention.

  • A nearest implementation without a bounded-search fast path will not scale.

  • Falling back to Shapely/STRtree for unsupported inputs must stay explicit enough that future GPU work is not hidden.

Intent

Define how repo-owned spatial query primitives plug into GeoPandas-compatible sindex and sjoin surfaces while preserving pandas-compatible output rules without making pandas assembly the internal execution model.

Decision

Land a repo-owned spatial-query engine and make relation results the internal join boundary, with explicit GeoPandas export only at the public surface.

This means:

  • the repo-owned spatial-query engine defines how future GPU-friendly query, dwithin, and bounded-nearest work should be staged

  • regular-grid, axis-aligned rectangle polygon indexes may register an explicit GPU point-query specialization instead of paying the generic bbox-tiling cost

  • vendored SpatialIndex stays on STRtree by default for current host execution because the owned host path is not yet performance-competitive

  • sjoin and sjoin_nearest build low-level RelationJoinResult objects first, then wrap them in a deferred export result that owns join context until the explicit GeoDataFrame boundary so suffixes, index restoration, and geometry-column rules stay GeoPandas-compatible

  • public sindex.query indices, dense, and sparse outputs format NativeRelation pairs when the owned/native index path is admissible; dense and sparse remain public compatibility exports, not native carriers

  • public sindex.query_any returns one eager boolean per query row. Its general path reduces bounded candidate/refine tiles directly to a device selection. For intersects against isotropic point buffers, one nearest-point distance and the query bounding-box diagonal certify definite matches and misses; only the threshold-ambiguous rows enter the buffer index’s exact range-sliced intersects semijoin. Inactive capacity rows carry NaN bounds and therefore schedule no Morton work. Temporary storage is O(query rows plus one bounded candidate tile), never O(relation pairs), and public boolean filtering exports only the selection count before keeping row positions and geometry on device

  • public sindex.query_aggregate consumes native relation pairs into eager, input-sized pandas ExtensionArray columns; public Series arithmetic can combine those already-computed device columns before explicit host export

  • public sindex.query_pair_aggregate reduces two position-aligned indexes against the same query geometries into left, right, and shared match-count columns; bounded candidate tiles are consumed before public construction, so aligned relation comparison does not export either pair relation

  • public sindex.nearest(..., k=N, return_all=False) formats an exact bounded NativeRelation for owned device point inputs; k=1 preserves GeoPandas behavior and unsupported shapes keep the existing explicit fallback path

  • RelationJoinExportResult -> NativeTabularResult now builds a native attribute payload directly, so terminal Arrow-family and file sinks do not need to go back through the joined-frame materializer just to emit join rows

  • those deferred join exports can also lower into the shared NativeTabularResult boundary so terminal writers do not need to rebuild a public frame just to emit join results

  • dwithin uses distance-expanded bounds as the coarse pass, then exact shapely.dwithin on the compacted candidate set

  • bounded nearest uses exact-distance reduction over compacted candidates

  • unsupported geometry families fall back to the original Shapely/STRtree path

Query Strategy

The owned engine provides:

  • build and cache a flat spatial index over owned geometry buffers

  • allow the generic flat-index order build to use GPU morton-key sort when the runtime selection explicitly requests GPU execution

  • use bbox candidate generation for scalar and bulk query

  • refine candidate rows with the exact binary-predicate engine

  • support indices, dense, and sparse output formats

For benchmark-style joins where the tree is a regular grid of non-overlapping rectangles and the query side is points, the engine may skip generic candidate generation entirely:

  • detect the regular-grid rectangle layout while building the flat index

  • launch a point-to-cell GPU kernel that emits up to four polygon hits per point for edge and corner cases

  • return polygon row ids directly because point-vs-rectangle intersects semantics match the cell lookup exactly

  • compact duplicate polygon ids before any downstream dissolve/materialization work so repeated point hits do not dominate CPU join assembly

The planned integration seam makes this query stack prefer index-driven coarse filtering where possible and routes point-vs-polygon refinement through the dedicated point predicate pipeline instead of always using the generic binary-predicate path.

This maps cleanly onto a future GPU/CCCL path:

  • bounds compare: transform-style primitive

  • candidate compaction: DeviceSelect

  • exact predicate refine on compacted rows

  • scatter back into requested output representation

Profile public-path query changes against cached GeoParquet fixtures in .benchmark_fixtures/ before trusting synthetic-only kernel rails. The fixture workflow is the source of truth for host/device traversal because it exercises read_parquet -> sindex -> query through the public GeoPandas surface.

Nearest Strategy

Nearest uses a three-tier dispatch for point-point data:

  1. Zero-copy GPU grid nearest — extracts coords directly from Shapely arrays (bypassing _to_owned), builds a uniform-grid spatial hash on device, and runs ring-expansion search entirely on the GPU. Handles both bounded and unbounded nearest. Falls through for non-point inputs or degenerate grids.

  2. Indexed GPU nearest — uses _to_owned arrays with a sorted-X sweep and CCCL lower/upper bound for tie counting. Fallback when the grid path declines.

  3. STRtree host path — used for non-point or mixed-geometry inputs, or when the GPU is unavailable.

The bounded path (max_distance set) can also use distance-expanded bounds to prune candidates, then exact distance reduction per input geometry.

Pandas Semantics

GeoPandas-compatible result assembly is an explicit export step from native relation results:

  • relation pairs stay separate from join/export context until the boundary

  • how semantics for left, right, and inner

  • index restoration and unnamed-index handling

  • overlapping-column suffix policy

  • geometry-column retention rules

  • optional distance-column insertion for nearest joins

Consequences

  • o17.4.3 now has one owned spatial-query engine that future GPU work can target

  • targeted upstream sindex and sjoin families validate the result semantics

  • the default host adapter avoids a measured regression by staying on STRtree for now

  • future GPU work can replace query/nearest internals without re-implementing DataFrame join behavior

  • downstream native workflows can consume join relation results without paying immediate pandas materialization