Add optional xarray/rasterio support for indices and sampling

docs/concepts.md

@@ -36,6 +36,45 @@ regardless of which field an importer populated, `scene.cloud_mask`, `scene.clou
 (cloud masking, cloud filtering) read from this normalised schema rather than product-specific
 field names.
 
+## `Scene.data` backing stores
+
+`earthrs` has zero *required* runtime dependencies — `Scene.data` as a `dict[str, list]` mapping
+band names to nested lists (or any list-of-lists-like structure) is the always-available default
+and is what every test in this project exercises:
+
+```python
+from earthrs import Scene
+
+scene = Scene(data={"nir": [[0.8, 0.4]], "red": [[0.2, 0.4]]}, band_names=["nir", "red"])
+```
+
+`earthrs.indices` and `earthrs.sampling` additionally recognise two optional, geospatial-ecosystem
+backing stores, guarded behind `try`/`except ImportError` so importing `earthrs` never requires
+either package:
+
+- **`xarray.DataArray`** — must have a `band` dimension whose coordinate values are the band
+  names, in the same order as `scene.band_names`. Band access is `data.sel(band=band_name)`.
+- **`xarray.Dataset`** — one data variable per band; band access is `data[band_name]`, which
+  returns a 2D `DataArray`.
+- **A rasterio dataset handle** — anything exposing `.read(band_index)` with 1-indexed band
+  numbers (for example an open `rasterio.io.DatasetReader`). Band names map to rasterio band
+  indices positionally: `scene.band_names[i]` corresponds to rasterio band `i + 1`.
+
+Passing anything else raises a `TypeError` with a message listing the supported shapes.
+
+Install optional support for the xarray/rasterio-backed paths with:
+
+```bash
+pip install "earthrs[geo]"
+```
+
+This is a deliberate project-wide policy: xarray, rasterio, and similar packages are optional
+integrations that *improve* behaviour (native array types, richer metadata, lazy/on-disk reads)
+without ever being required — the zero-dependency `dict`-of-list path keeps working identically
+whether or not the `geo` extra is installed. See [Project status](status.md) for known
+limitations of the xarray/rasterio paths (for example, dimension-ordering conventions that are
+not yet recognised).
+
 ## Spectral indices
 
 `earthrs.indices` provides `ndvi`, `ndwi`, and `evi`, each taking a `Scene` and returning

docs/status.md

@@ -6,11 +6,18 @@ scaffolded but not yet built — see [Getting started](getting-started.md) and
 
 ## Known limitations
 
-- Sampling (`earthrs.sampling`) and spectral indices (`earthrs.indices`) currently only support
-  `Scene.data` as a mapping of band name to nested list/array data, not the xarray/rasterio-backed
-  data the `Scene` docstring describes as the typical case.
+- Sampling (`earthrs.sampling`) and spectral indices (`earthrs.indices`) support `Scene.data` as
+  a `dict[str, list]` of band name to nested list data (the always-available, zero-dependency
+  default), an `xarray.DataArray`/`Dataset`, or a rasterio dataset handle — see
+  [Core concepts](concepts.md#scenedata-backing-stores) for the exact conventions and how to
+  install optional support.
+- `xarray.DataArray` scene data must use a `band` dimension whose coordinate values match
+  `scene.band_names`; other dimension-ordering or naming conventions (for example a `variable`
+  dimension, or bands stored as separate leading axes without coordinates) are not recognised
+  and raise `TypeError`.
 - Spectral indices treat only Python `None` as a missing value; there is no nodata/fill-value
-  (for example `-9999`) handling yet.
+  (for example `-9999`) handling yet. This applies uniformly across dict, xarray, and rasterio
+  backed data — no backend currently reads a dataset's own nodata value.
 
 ## Scaffolded, not yet implemented
 

pyproject.toml

@@ -34,6 +34,10 @@ docs = [
   "mkdocs>=1.6",
   "mkdocs-material>=9.5",
 ]
+geo = [
+  "xarray",
+  "rasterio",
+]
 
 [project.urls]
 Homepage = "https://github.com/ysims/earthrs"

src/earthrs/indices.py

@@ -11,6 +11,16 @@
 
 from earthrs.scene import Scene
 
+try:
+    import rasterio
+except ImportError:  # pragma: no cover - exercised only when rasterio is absent
+    rasterio = None
+
+try:
+    import xarray
+except ImportError:  # pragma: no cover - exercised only when xarray is absent
+    xarray = None
+
 
 def ndvi(scene: Scene, *, nir_band: str = "nir", red_band: str = "red") -> Any:
     """Compute NDVI from a scene.
@@ -64,12 +74,48 @@ def evi(
 
 
 def _get_band(scene: Scene, band_name: str) -> Any:
-    if isinstance(scene.data, dict):
-        if band_name not in scene.data:
+    """Return band data as nested Python lists, regardless of the backing store.
+
+    Supports the following ``scene.data`` shapes:
+
+    - ``dict[str, list]`` (the always-available, zero-dependency default): returned as-is.
+    - ``xarray.DataArray`` with a ``band`` dimension whose coordinate values are band names
+      (matching the order of ``scene.band_names``), or ``xarray.Dataset`` with one data
+      variable per band: converted to nested lists via ``.values.tolist()``.
+    - A rasterio dataset handle (anything exposing ``.read(band_index)``, 1-indexed): band
+      names map to rasterio band indices by position, so ``scene.band_names[i]`` corresponds
+      to rasterio band ``i + 1``. Converted to nested lists via ``.tolist()``.
+
+    The result is always plain nested lists so downstream elementwise helpers (`_map_binary`,
+    `_map_ternary`) do not need to know about the original backing store.
+    """
+
+    data = scene.data
+    if isinstance(data, dict):
+        if band_name not in data:
+            raise ValueError(f"Band '{band_name}' not found in scene data.")
+        return data[band_name]
+    if xarray is not None and isinstance(data, xarray.Dataset):
+        if band_name not in data.data_vars:
+            raise ValueError(f"Band '{band_name}' not found in scene data.")
+        return data[band_name].values.tolist()
+    if xarray is not None and isinstance(data, xarray.DataArray):
+        if "band" not in data.dims or "band" not in data.coords:
+            raise TypeError(
+                "xarray.DataArray scene data must have a 'band' dimension with a matching "
+                "coordinate to support band access."
+            )
+        if band_name not in data.coords["band"].values.tolist():
+            raise ValueError(f"Band '{band_name}' not found in scene data.")
+        return data.sel(band=band_name).values.tolist()
+    if rasterio is not None and hasattr(data, "read") and hasattr(data, "count"):
+        if band_name not in scene.band_names:
             raise ValueError(f"Band '{band_name}' not found in scene data.")
-        return scene.data[band_name]
+        band_index = scene.band_names.index(band_name) + 1
+        return data.read(band_index).tolist()
     raise TypeError(
-        "Index calculations currently expect scene data as a mapping of band names to arrays."
+        "Index calculations currently expect scene data as a mapping of band names to arrays, "
+        "an xarray DataArray/Dataset, or a rasterio dataset."
     )
 
 

src/earthrs/sampling/core.py

@@ -7,6 +7,16 @@
 from earthrs.samples import Samples
 from earthrs.scene import Scene
 
+try:
+    import rasterio
+except ImportError:  # pragma: no cover - exercised only when rasterio is absent
+    rasterio = None
+
+try:
+    import xarray
+except ImportError:  # pragma: no cover - exercised only when xarray is absent
+    xarray = None
+
 
 def sample_points(scene: Scene, points: Any, *, method: str = "nearest") -> Samples:
     """Extract raster values at point locations.
@@ -187,15 +197,64 @@ def _coerce_transects(transects: Any) -> list[dict[str, list[tuple[int, int]]]]:
 
 
 def _get_band(scene: Scene, band_name: str) -> Any:
-    if isinstance(scene.data, dict):
-        if band_name not in scene.data:
+    """Return a single band's 2D data, regardless of the backing store.
+
+    Supports the following ``scene.data`` shapes:
+
+    - ``dict[str, list]`` (the always-available, zero-dependency default): returned as-is.
+    - ``xarray.DataArray`` with a ``band`` dimension whose coordinate values are band names
+      (matching the order of ``scene.band_names``): selected via ``.sel(band=band_name)``,
+      returning a 2D ``DataArray``.
+    - ``xarray.Dataset`` with one data variable per band: band access via ``data[band_name]``,
+      returning a 2D ``DataArray``.
+    - A rasterio dataset handle (anything exposing ``.read(band_index)``, 1-indexed): band
+      names map to rasterio band indices by position, so ``scene.band_names[i]`` corresponds
+      to rasterio band ``i + 1``. Returns a 2D NumPy array from ``.read(band_index)``.
+    """
+
+    data = scene.data
+    if isinstance(data, dict):
+        if band_name not in data:
+            raise ValueError(f"Band '{band_name}' not found in scene data.")
+        return data[band_name]
+    if xarray is not None and isinstance(data, xarray.Dataset):
+        if band_name not in data.data_vars:
+            raise ValueError(f"Band '{band_name}' not found in scene data.")
+        return data[band_name]
+    if xarray is not None and isinstance(data, xarray.DataArray):
+        if "band" not in data.dims or "band" not in data.coords:
+            raise TypeError(
+                "xarray.DataArray scene data must have a 'band' dimension with a matching "
+                "coordinate to support band access."
+            )
+        if band_name not in data.coords["band"].values.tolist():
             raise ValueError(f"Band '{band_name}' not found in scene data.")
-        return scene.data[band_name]
-    raise TypeError("Sampling currently expects scene data as a mapping of band names to arrays.")
+        return data.sel(band=band_name)
+    if rasterio is not None and hasattr(data, "read") and hasattr(data, "count"):
+        if band_name not in scene.band_names:
+            raise ValueError(f"Band '{band_name}' not found in scene data.")
+        band_index = scene.band_names.index(band_name) + 1
+        return data.read(band_index)
+    raise TypeError(
+        "Sampling currently expects scene data as a mapping of band names to arrays, an "
+        "xarray DataArray/Dataset, or a rasterio dataset."
+    )
 
 
 def _read_cell(data: Any, row: int, col: int) -> Any:
-    return data[row][col]
+    """Read a single value out of 2D band data, regardless of the backing store.
+
+    Nested Python lists/tuples (the zero-dependency default) are indexed exactly as before:
+    ``data[row][col]``. Anything else (an ``xarray.DataArray`` band slice, or a NumPy-like 2D
+    array as returned by rasterio's ``.read(band_index)``) is assumed to support NumPy-style
+    tuple indexing (``data[row, col]``); the result is unwrapped to a native Python scalar via
+    ``.item()`` when available, so callers get plain floats/ints regardless of backend.
+    """
+
+    if isinstance(data, (list, tuple)):
+        return data[row][col]
+    value = data[row, col]
+    return value.item() if hasattr(value, "item") else value
 
 
 def _reduce_values(values: list[Any], reducer: str) -> float:

src/earthrs/scene.py

@@ -17,7 +17,11 @@ class Scene:
     Parameters
     ----------
     data:
-        Raster data, typically an xarray ``DataArray`` or ``Dataset``.
+        Raster data. The always-available, zero-dependency default is a ``dict[str, list]``
+        mapping band names to nested-list arrays. ``earthrs.indices`` and ``earthrs.sampling``
+        also accept an xarray ``DataArray``/``Dataset`` or a rasterio dataset handle when the
+        optional ``geo`` extra is installed (``pip install "earthrs[geo]"``) — see
+        `docs/concepts.md`'s "Scene.data backing stores" section for the exact conventions.
     crs:
         Coordinate reference system for the scene.
     transform: