Add spatial parquet points loading

.gitignore

@@ -13,11 +13,10 @@ packages/zarrextra/src/*.d.ts.map
 coverage/
 test-fixtures/
 validation-results/
+**/.pytest_cache/
 
-# Python virtual environments (version-specific, completely separate)
-python/v0.5.0/.venv/
-python/v0.6.1/.venv/
-python/v0.7.2/.venv/
+# Python virtual environments
+python/*/.venv/
 
 # Vendored from @cornerstonejs/codec-openjph (see scripts/vendor-openjph-for-python.mjs)
 python/spatialdata-codec-writer/src/spatialdata_codec_writer/vendor/openjph/

CONTEXT.md

@@ -40,6 +40,22 @@ _Avoid_: serializable prop, stack entry prop
 A small MDV/control-layer UI area that edits observable state directly while passing plain values through renderer and third-party boundaries.
 _Avoid_: MobX renderer contract
 
+**Points Render Resource**:
+The **Resource Resolver** output for a points **Spatial Entry**: a bundle `{ element, loader }` pairing the canonical `PointsElement` with a frozen **`PointsLoader`** facet.
+_Avoid_: treating `PointsLoader` alone as the full render resource, or storing the loader on/mutating the element
+
+**PointsElement**:
+The source identity and public source API for SpatialData points. It owns element metadata and source-level methods such as full or bounded point loading; deck render strategies receive the resolved **Points Render Resource** instead of mutating or calling the element directly.
+_Avoid_: render-time strategy object, mutable loader holder
+
+**PointsLoader**:
+The loader facet of a **Points Render Resource**: encoding capabilities plus a fetch API (`loadInBounds`, optional `loadAll`). Built by `@spatialdata/core` store-I/O factories; consumed by `@spatialdata/layers` render strategies — not by calling `PointsElement` methods directly from deck code.
+_Avoid_: conflating with Viv/image `loader` when discussing SpatialData element identity
+
+**Points Encoding**:
+The render-time points layout selected after resolver probing, e.g. `preloaded-columnar`, `morton-tiled`, or future `geoarrow-*` kinds. Distinct from persisted Parquet layout described in ADR 0002.
+_Avoid_: `experimentalOptimizations` as a synonym for encoding kind
+
 ## Relationships
 
 - A **Render Stack** contains zero or more ordered **Stack Entries**.
@@ -59,3 +75,12 @@ _Avoid_: MobX renderer contract
 - "Layer" was used for SpatialData elements, deck.gl layer instances, UI rows, and saved config entries. Resolved: use **Stack Entry** for saved/render order, **Spatial Entry** for SpatialData-backed entries, and deck.gl layer only for runtime renderer output.
 - "Snapshot" was used for both persisted config and temporary UI/render state. Resolved: persisted config is a **Render Stack**; live MDV direct-edit areas are **MobX Control Islands** and should not periodically snapshot the whole stack during interaction.
 - "Props" was used for both serialized renderer inputs and runtime callback objects. Resolved: `entry.props` must remain serializable renderer input; listeners, factories, portals, and raw deck integration points are **Runtime Attachments**.
+
+## In-progress work
+
+- **Points preload & feature filter** — runtime filter on preloaded scatter, row-group
+  geometry reads, catalog fallbacks for large dictionary-only transcripts.
+  Status and cleanup plan:
+  [`docs/plans/points-preload-feature-filter-status.md`](docs/plans/points-preload-feature-filter-status.md).
+  Related ADRs: [0002](docs/adr/0002-spatially-aware-vector-loading.md),
+  [0003](docs/adr/0003-points-render-resource.md).

docs/adr/0002-spatially-aware-vector-loading.md

@@ -0,0 +1,148 @@
+# Spatially-Aware Vector Loading
+
+SpatialData points and shapes can be large enough that whole-element Parquet
+loads are not a viable browser default. We will treat viewport-bounded vector
+loading as a first-class source API and keep persisted optimization artifacts in
+Parquet/GeoParquet rather than inventing a deck.gl-specific storage format.
+
+## Decision
+
+- Points v1 follows current Vitessce practice: a SpatialData Points Parquet
+  element may be sorted by 2D Morton order with a `morton_code_2d` column, a
+  feature-code column, controlled row-group sizes, and 2-4 leading sentinel rows
+  whose `morton_code_2d` is `0` and whose coordinates encode the full point
+  extent.
+- `@spatialdata/core` exposes bounded point loading through
+  `PointsElement.loadPointsInBounds()`. When the Parquet module supports
+  Vitessce's row-group APIs (`readMetadata` and `readParquetRowGroup`) and the
+  store supports range reads, the loader may fetch selected row groups. Otherwise
+  it degrades to the existing full-table read followed by bounds filtering.
+- Render-time code uses ADR 0003's **Points Render Resource** (`{ element,
+  loader }`) and calls the `PointsLoader` facet. `PointsElement` remains source
+  identity and public source API, not the deck strategy contract.
+- `@spatialdata/vis` may render compatible points through a deck.gl `TileLayer`.
+  The tile layer owns async viewport loads and abort signals; ordinary
+  `ScatterplotLayer` rendering remains the fallback for preloaded point data.
+- `points.experimental/<key>` and `shapes.experimental/<key>` are reserved as
+  top-level Experimental Optimization Collections. They link back to the source
+  element by key and metadata rather than modifying canonical SpatialData
+  element semantics.
+- GeoParquet is the durable shape optimization target. GeoArrow is a runtime
+  columnar layout / deck adapter option, not a duplicate persisted artifact.
+
+## Experimental Optimization Collections
+
+Use `points.experimental/<key>/` and `shapes.experimental/<key>/` only for
+persisted layouts that **standard SpatialData / Vitessce readers cannot correctly
+consume** — not for every browser optimization.
+
+| Layout | Where it lives | Why |
+|--------|----------------|-----|
+| **Morton v1** (`morton_code_2d`, sentinels, `{feature_key}_codes`, row groups) | **Canonical** `points/<key>/points.parquet` | Follows Vitessce practice. Extra columns are additive; Python `spatialdata` full-table reads still work. |
+| **Feature-primary sort** (Morton not primary key) | `points.experimental/<key>/` | Breaks Morton row-group bisect; needs a new tiling `kind` |
+| **Padua multiscale** (`__spatial_index__`, levels in schema metadata) | `points.experimental/<key>/` | Non-standard vs morton-points v1 |
+| **GeoParquet shapes tiling** | `shapes.experimental/<key>/` | Future |
+
+`experimentalOptimizations` in `@spatialdata/vis` means use TileLayer / row-group
+reads when **canonical** parquet schema supports morton tiling — not “look in
+`points.experimental/`”.
+
+The experimental writer defaults to **in-place** Morton sorting on
+`points/<key>/points.parquet`. Use `--experimental` only when writing a layout
+that must not replace the canonical element.
+
+## Multi-part Parquet (reader)
+
+Wild-type SpatialData points may store `points/<key>/points.parquet` as a
+**directory** with `part.0.parquet`, `part.1.parquet`, … The logical path remains
+`points/<key>/points.parquet`. `@spatialdata/core` supports both single-file and
+multipart layouts for metadata, schema, and row-group range reads. The
+experimental writer outputs a **single-file** Morton artifact by design; row-group
+range reads fetch only the byte ranges needed per viewport.
+
+## Feature / gene filtering
+
+Transcript and other feature-bearing points declare `feature_key` in element
+`spatialdata_attrs` (for example `"feature_name"` on xenium transcripts). This is
+distinct from `instance_key` (for example `"cell_id"`), which identifies the
+object a point belongs to.
+
+The Morton writer adds `{feature_key}_codes` (for example `feature_name_codes`)
+as `int32` categorical codes alongside the string feature column. Sorting is
+**spatial** (Morton on x/y) by default; row groups are spatial chunks.
+
+**Core API** — extend bounded loading with optional feature codes:
+
+```typescript
+interface PointsInBoundsOptions {
+  bounds: SpatialBounds;
+  /** Integer codes matching `{feature_key}_codes` in the parquet artifact */
+  featureCodes?: readonly number[];
+  signal?: AbortSignal;
+}
+```
+
+**Vis API** — extend `PointsLayerConfig` with `featureCodes?: number[]` and
+wire through TileLayer `updateTriggers.getTileData`.
+
+v1 applies feature filtering as a **read-time row predicate** after spatial
+bounds filtering (and after row-group fetch on the Morton path). It does not
+skip row groups by gene. String-based `features?: string[]` and a codebook
+artifact are deferred.
+
+**Implementation status** (preload vs runtime filter, catalog, workers):
+[`docs/plans/points-preload-feature-filter-status.md`](../plans/points-preload-feature-filter-status.md).
+
+Feature filtering is separate from **feature-primary sort** experiments
+(`[feature_codes, morton_code_2d]`), which may require a new tiling `kind` if
+promoted. Use `write-index-permutations` on a derivative Zarr store to benchmark
+sort strategies; see the writer README.
+
+A hypothetical **per-gene density map** (2D histogram / KDE for one feature) is
+out of scope for the Morton tile path and may be an offline aggregation or
+dedicated viz mode later.
+
+## Sort strategy experiments
+
+Default Morton v1 sort is spatial on `morton_code_2d` (optionally `z` when
+low-cardinality). Multi-key sorts under evaluation include
+`[morton_code_2d, feature_name_codes]` and `[feature_name_codes, morton_code_2d]`.
+The reader's row-group bisect assumes Morton is the **primary** sort key; do not
+silently swap sort order under the existing `morton-points` format id.
+
+Generate comparable permutations with:
+
+```bash
+spatialdata-experimental-writer write-index-permutations SOURCE_ZARR DEST_ZARR
+```
+
+The derivative store includes sibling `points/<condition>/` elements and
+`index-manifest.json` for benchmark tooling.
+
+## Prior Art
+
+- scverse Padua hackathon points work:
+  <https://github.com/scverse/2026_04_hackathon_padua/issues/17> and
+  <https://github.com/scverse/2026_04_hackathon_padua/issues/24>.
+- Padua branch prototype:
+  <https://github.com/scverse/2026_04_hackathon_padua/tree/viz/point_chunking/visualization>.
+- Vitessce tiled SpatialData Points:
+  <https://github.com/vitessce/vitessce/pull/2286>.
+- Vitessce sentinel bbox update:
+  <https://github.com/vitessce/vitessce/issues/2419> and
+  <https://github.com/vitessce/vitessce/pull/2489>.
+- Vitessce shapes format `0.3` compatibility:
+  <https://github.com/vitessce/vitessce/pull/2495> and
+  <https://github.com/vitessce/vitessce/releases/tag/v3.9.11>.
+
+## Consequences
+
+- Source loaders must expose typed/columnar batches and remain independent of
+  deck.gl. Rendering packages decide whether to use TileLayer, ScatterplotLayer,
+  or a future GeoArrow-aware layer.
+- Whole-table point loading is still supported and is the compatibility fallback,
+  but render paths can opt into experimental optimizations with a single
+  `experimentalOptimizations` switch.
+- Shapes format `0.3` remains on the current modern Parquet-backed path in
+  `VShapesSource`; large-shape spatial tiling still needs a separate GeoParquet
+  artifact/writer slice.

docs/adr/0003-points-render-resource.md

@@ -0,0 +1,97 @@
+# Points Render Resource
+
+ADR 0002 describes persisted Morton Parquet artifacts and bounded loading APIs on
+`PointsElement`. This ADR describes the **render-time** boundary between store
+I/O, the Resource Resolver, and the deck.gl `PointsLayer` composite.
+
+## Decision
+
+- A points **Spatial Entry** (`PointsElement`) remains the canonical spatial
+  identity handle. Deck layers stay associated with that element for picks,
+  tooltips, and Render Stack `elementKey`.
+- The **Resource Resolver** (today `resolvePointsRenderResource()` in
+  `@spatialdata/vis`) probes once and returns a **Points Render Resource**
+  bundle `{ element, loader }` with **frozen** encoding capabilities.
+- **`PointsLoader`** is the loader facet only: encoding kind, batch format,
+  bounds, and fetch methods. Render strategies call `loader.loadInBounds()` —
+  not `element.loadPointsInBounds()` directly from `@spatialdata/layers`.
+- **`PointsLayer`** (`@spatialdata/layers` `CompositeLayer`) takes
+  `resource: PointsRenderResource` plus cosmetic props. It delegates to
+  encoding-specific render strategies selected by `loader.capabilities.kind`.
+- **Store I/O loader factories** live in `@spatialdata/core` and close over
+  `PointsElement`. **Render strategies** and tile-debug overlay logic live in
+  `@spatialdata/layers`. The vis resolver associates element + loader.
+
+## Encoding selection (v1)
+
+| Condition | Encoding kind | Strategy |
+|-----------|---------------|----------|
+| Full table preloaded in resolver cache | `preloaded-columnar` | `ScatterplotLayer` |
+| Morton metadata with row-group range reads + bounds | `morton-tiled` | `TileLayer` + per-tile scatter |
+| Future GeoArrow batch from core | `geoarrow-binary` | stub → `GeoArrowScatterplotLayer` |
+| Future tiled Arrow/Parquet deck path | `geoarrow-tiled` | stub |
+
+Resolver probing is **eager**: capabilities do not change mid-session unless
+the element or resolver cache inputs change.
+
+## GeoArrow boundary
+
+- **Core** may later expose deck-free Apache Arrow `RecordBatch` batches from
+  Parquet row groups (x/y/z columns or geometry).
+- **Layers** owns [deck.gl-geoarrow](https://github.com/geoarrow/deck.gl-geoarrow)
+  integration: GeoArrow geometry shaping and `GeoArrowScatterplotLayer` /
+  future tiled deck paths.
+- Core must not import deck.gl or `@geoarrow/deck.gl-geoarrow`.
+
+## Batch contract
+
+`PointsBatch` is a tagged union:
+
+- `columnar-ndarray` — v1 Morton and preloaded paths
+- `arrow-record-batch` — reserved for GeoArrow strategies
+
+## Tile debug overlay
+
+When `showTileDebugOverlay` is enabled on a tiled encoding, the morton strategy
+emits a pickable `PolygonLayer` sublayer with per-tile status (pending, loading,
+loaded, empty, error, aborted). This is cosmetic for tile fetching and must not
+appear in `TileLayer.updateTriggers.getTileData`.
+
+## Relationship to ADR 0002
+
+- ADR 0002: persisted artifacts and source-level
+  `PointsElement.loadPointsInBounds()` API.
+- ADR 0003: render-time bundle, strategy registry, and deck composite ownership.
+
+## Consequences
+
+- Swapping encodings or deck.gl parquet layers requires new loader factories
+  and/or strategies — not changes to `PointsLayer` public props.
+- `PointsElement` does not grow a mutable `renderResource` attachment; the
+  resolver cache holds stable bundle references per element key.
+- Image precedent: `ImageElement` + Viv loader built in vis; points precedent:
+  `PointsElement` + `PointsLoader` built in vis, rendered by `PointsLayer`.
+
+## Future performance investigations
+
+These are documented follow-ups — not part of the v1 render bundle.
+
+### CPU / compute hot paths
+
+The scan+compact loops in `filterColumnarByFeatureCodes` /
+`filterPointsToBounds` (`packages/core/src/pointsTiling.ts`) are hot paths for
+large preloaded datasets. Candidates include WASM SIMD and WebGPU compute (e.g.
+[typegpu](https://github.com/software-mansion/typegpu)) for parallel index
+selection and column compaction. **Worker offload** (`@spatialdata/core/workers`)
+is the near-term mitigation; GPU/WASM is a follow-up benchmark task.
+
+### FBO-based render caching
+
+For viewport-stable layers (tiled points, filtered preloaded batches, static
+image tiles), cache rasterized sublayer output in **framebuffer objects (FBOs)**
+so pan/zoom and cosmetic prop changes do not re-draw the full payload every
+frame. This should integrate with the broader **Render Stack compositing**
+story (`Group Entry`, Viv/deck stacking) via shared FBO cache utilities —
+invalidation keyed on structural `updateTriggers`, composition order with host
+overlays — rather than as a points-only hack. Detail deferred until compositing
+utils exist.

docs/docs/core/error-handling.mdx

@@ -4,22 +4,41 @@ sidebar_position: 4
 
 # Error Handling
 
-The `@spatialdata/core` package uses a `Result` type pattern inspired by Rust for explicit error handling. This approach makes error cases visible in the type system and avoids unexpected exceptions.
+`@spatialdata/core` exports a `Result` type (inspired by Rust) for operations
+where failure is expected and should be handled explicitly. **Adoption is
+currently very limited** — a handful of APIs (notably coordinate-system
+lookups) return `Result`; most of the codebase still throws exceptions, returns
+`null`, or uses other ad hoc patterns.
 
-:::info Error types and safety
+Because error handling is not uniform, you **cannot treat `Result` as a
+system-wide contract**. When integrating with this library, assume exceptions
+unless a specific API's signature returns `Result`. The sections below document
+the type and patterns for the places that use it.
 
-The use of this pattern in some parts of the code should not be taken as an assertion that all possible exceptions are handled in this way.
+We find the `Result` approach appealing in principle, but wider use is an **open
+design decision**. If we lean in more thoroughly, we would likely adopt an
+**established library** (for example [`neverthrow`](https://github.com/supermacro/neverthrow)) rather than extend our small in-house definitions.
 
-This aspect of the design may be subject to review, and community feedback is encouraged.
+:::info Limited adoption — read this first
 
-:::
-
-:::note Result implementation
+The `Result` pattern is **not** adhered to consistently across packages. Do not
+infer from this page that most failures are typed, composable, or exception-free.
+Many code paths (including parquet I/O, element loading, and worker boundaries)
+still throw or use nullable returns.
 
-The `Result` type is implemented in `zarrextra` and re-exported from `@spatialdata/core` for convenience. This is currently a custom implementation for simplicity and to avoid dependencies. We may review using an existing Result library (such as `neverthrow`) in the future, but for now this provides a lightweight, dependency-free solution.
+Community feedback on whether and how to expand `Result` use is welcome.
 
 :::
 
+:::note Current implementation (may change)
+
+Today the `Result` type lives in `zarrextra` and is re-exported from
+`@spatialdata/core` — a minimal custom shape (`ok` / `value` / `error`) kept
+dependency-free while adoption remains narrow. **Do not build heavily on these
+definitions** as a long-term API promise; broader adoption would probably replace
+or wrap them with a maintained library and standard combinators (`map`, `andThen`,
+etc.).
+
 ## The Result Type
 
 ```ts
@@ -212,13 +231,27 @@ if (errors.length > 0) {
 }
 ```
 
-## Why Result Instead of Exceptions?
+## Why use Result where we do?
+
+These are motivations for the pattern **where it is used today**, and for
+possible wider adoption if we decide to lean in:
 
 1. **Explicit error handling**: The type system shows you when operations can fail
 2. **No hidden control flow**: Errors don't jump up the call stack unexpectedly
 3. **Self-documenting**: Function signatures show error types
 4. **Composable**: Easy to chain, map, and aggregate results
 5. **Performance**: Avoids the overhead of exception creation when errors are common
 
-This is especially valuable for coordinate system lookups where it's common and expected that some systems may not be available for all elements.
+This is especially valuable for coordinate system lookups, where it is common
+and expected that some systems may not be available for all elements.
+
+## Parquet I/O (not yet migrated)
+
+Parquet loading in `VTableSource` / `VPointsSource` is a typical example of the
+inconsistency above: exceptions, `null`, and skip-on-miss policies coexist.
+Call sites disagree on whether a missing part should throw or be skipped. We are
+keeping this as-is for now. See
+[`docs/plans/parquet-io-error-handling.md`](../../plans/parquet-io-error-handling.md)
+for current behaviour, rationale, and a suggested migration order if we expand
+`Result` use.