Merge pull request #29 from CU-ESIIL/codex/add-load_sentinel2_ndvi_cube-function

Add Sentinel-2 NDVI helper

Author: ttuff

README.md

@@ -4,7 +4,7 @@ CubeDynamics is a streaming-first climate cube math library with ggplot-style pi
 
 ## Features
 
-- **Streaming PRISM/gridMET/Sentinel-2 helpers** (`cd.load_prism_cube`, `cd.load_gridmet_cube`, `cd.load_s2_ndvi_cube`) for immediate analysis without bulk downloads.
+- **Streaming PRISM/gridMET/Sentinel-2 helpers** (`cd.load_prism_cube`, `cd.load_gridmet_cube`, `cd.load_s2_ndvi_cube`, `cd.load_sentinel2_ndvi_cube`) for immediate analysis without bulk downloads.
 - **Climate variance, correlation, trend, and synchrony cubes** that run on `xarray` objects and scale from laptops to clusters.
 - **Pipe + verb system** – build readable cube workflows with `pipe(cube) | v.month_filter(...) | v.variance(...)` syntax inspired by ggplot/dplyr.
 - **Verbs namespace (`cubedynamics.verbs`)** so transforms, stats, IO, and visualization live in focused modules.
@@ -95,10 +95,33 @@ pipe(cube) | v.month_filter([6, 7, 8]) | v.variance(dim="time")
 ### Sentinel-2 → NDVI Anomaly (z-score) Cube
 
 CubeDynamics works on remote-sensing image stacks in addition to climate
-archives. A typical Sentinel-2 workflow is:
+archives. The convenience helper `cd.load_sentinel2_ndvi_cube` streams
+Sentinel-2 Level-2A chips via [`cubo`](https://github.com/carbonplan/cubo),
+computes NDVI from B08 (NIR) and B04 (red), and standardizes the result across
+time. Install `cubo` alongside CubeDynamics to use the helper:
 
-1. Stream Level-2A chips with [`cubo`](https://github.com/carbonplan/cubo)
-   using bands B04 (red) and B08 (NIR).
+```python
+import cubedynamics as cd
+from cubedynamics import pipe, verbs as v
+
+ndvi_z = cd.load_sentinel2_ndvi_cube(
+    lat=40.0,
+    lon=-105.25,
+    start="2018-01-01",
+    end="2020-12-31",
+)
+
+pipe(ndvi_z) | v.show_cube_lexcube(title="Sentinel-2 NDVI z-score")
+```
+
+`load_sentinel2_ndvi_cube` returns a `(time, y, x)` NDVI z-score cube ready for
+the rest of the verbs API. Pass ``return_raw=True`` to also receive the raw
+Sentinel-2 reflectance stack and intermediate NDVI cube.
+
+If you prefer to customize every step, the helper simply wraps the manual
+pipeline below:
+
+1. Stream Level-2A chips with `cubo` using bands B04 (red) and B08 (NIR).
 2. Convert the reflectance cube to NDVI with `v.ndvi_from_s2(...)`.
 3. Standardize NDVI over time with `v.zscore(dim="time")` to highlight
    greenness anomalies.
@@ -119,7 +142,6 @@ LON = -102.18
 START = "2023-06-01"
 END = "2024-09-30"
 
-# 1. Stream Sentinel-2 reflectance without local downloads
 with warnings.catch_warnings():
     warnings.simplefilter("ignore")
     s2 = cubo.create(
@@ -134,19 +156,16 @@ with warnings.catch_warnings():
         query={"eo:cloud_cover": {"lt": 40}},
     )
 
-# 2. Pipe reflectance -> NDVI -> z-scores
 ndvi_z = (
     pipe(s2)
     | v.ndvi_from_s2(nir_band="B08", red_band="B04")
     | v.zscore(dim="time")
 ).unwrap()
 
-# 3. Optional: visualize and QA in notebooks
 (pipe(ndvi_z) | v.show_cube_lexcube(title="Sentinel-2 NDVI z-score", clim=(-3, 3)))
 median_series = ndvi_z.median(dim=("y", "x"))
 median_series.plot.line(x="time", ylabel="Median NDVI z-score")
 
-# 4. Optional: correlate with a PRISM anomaly cube at each pixel
 corr_cube = (
     pipe(ndvi_z)
     | v.correlation_cube(prism_anom_cube, dim="time")
@@ -189,7 +208,7 @@ Lexcube widgets require a live Python environment (Jupyter, Colab, Binder). They
 
 - `pipe` for wrapping any cube before piping.
 - `verbs` (``from cubedynamics import verbs as v``) exposes transforms, statistics, IO, and visualization helpers.
-- Streaming helpers: `cd.load_prism_cube`, `cd.load_gridmet_cube`, `cd.load_s2_cube`, `cd.load_s2_ndvi_cube`.
+- Streaming helpers: `cd.load_prism_cube`, `cd.load_gridmet_cube`, `cd.load_s2_cube`, `cd.load_s2_ndvi_cube`, `cd.load_sentinel2_ndvi_cube`.
 - Vegetation helper: `v.ndvi_from_s2` for direct NDVI calculation on Sentinel-2 cubes.
 - Stats verbs: `v.anomaly`, `v.month_filter`, `v.variance`, `v.zscore`, `v.correlation_cube`, and more under `cubedynamics.ops.*`.
 - IO verbs: `v.to_netcdf`, `v.to_zarr`, etc.

docs/streaming_sources.md

@@ -66,44 +66,31 @@ cube = cd.load_gridmet_cube(
 ## Sentinel-2 → NDVI anomaly (z-score) cube
 
 Remote-sensing chips stream into the same pipe + verbs grammar, so you can
-combine vegetation signals with climate anomalies. Sentinel-2 Level-2A imagery
-exposes B04 (red) and B08 (NIR) bands, which flow into
-`v.ndvi_from_s2(...)` and `v.zscore(...)` inside a pipe chain.
+combine vegetation signals with climate anomalies. Use
+`cd.load_sentinel2_ndvi_cube` (requires the `cubo` package) for a one-function
+workflow:
 
 ```python
-from __future__ import annotations
-
-import warnings
-
-import cubo
-
+import cubedynamics as cd
 from cubedynamics import pipe, verbs as v
 
-with warnings.catch_warnings():
-    warnings.simplefilter("ignore")
-    s2 = cubo.create(
-        lat=43.89,
-        lon=-102.18,
-        collection="sentinel-2-l2a",
-        bands=["B04", "B08"],
-        start_date="2023-06-01",
-        end_date="2024-09-30",
-        edge_size=512,
-        resolution=10,
-        query={"eo:cloud_cover": {"lt": 40}},
-    )
-
-ndvi_z = (
-    pipe(s2)
-    | v.ndvi_from_s2(nir_band="B08", red_band="B04")
-    | v.zscore(dim="time")
-).unwrap()
+ndvi_z = cd.load_sentinel2_ndvi_cube(
+    lat=40.0,
+    lon=-105.25,
+    start="2018-01-01",
+    end="2020-12-31",
+)
 
-(pipe(ndvi_z) | v.show_cube_lexcube(title="Sentinel-2 NDVI z-score", clim=(-3, 3)))
-median_series = ndvi_z.median(dim=("y", "x"))
-median_series.plot.line(x="time", ylabel="Median NDVI z-score")
+pipe(ndvi_z) | v.show_cube_lexcube(title="Sentinel-2 NDVI z-score")
 ```
 
+The helper streams Sentinel-2 Level-2A imagery via `cubo`, computes NDVI from
+bands B08 (NIR) and B04 (red), and runs `v.zscore(dim="time")` so the returned
+cube is standardized across time. Set ``return_raw=True`` to also grab the raw
+reflectance stack and intermediate NDVI cube. You can still manually reproduce
+the pipeline with `pipe(...) | v.ndvi_from_s2(...) | v.zscore(...)` if you need a
+different stacking order.
+
 The resulting cube highlights unusual greenness events (drought stress,
 disturbance, rapid recovery). Because every cube shares `(time, y, x)` axes, you
 can correlate NDVI anomalies with PRISM or gridMET cubes using

src/cubedynamics/__init__.py

@@ -28,6 +28,7 @@
 from .utils.chunking import coarsen_and_stride
 from .viz.lexcube_viz import show_cube_lexcube
 from .viz.qa_plots import plot_median_over_space
+from .sentinel import load_sentinel2_ndvi_cube
 
 # Streaming-first stubs for the new architecture ---------------------------------
 from .streaming.gridmet import stream_gridmet_to_cube
@@ -50,6 +51,7 @@
     "load_s2_ndvi_cube",
     "load_gridmet_cube",
     "load_prism_cube",
+    "load_sentinel2_ndvi_cube",
     "compute_ndvi_from_s2",
     "zscore_over_time",
     "temporal_anomaly",

src/cubedynamics/sentinel.py

@@ -0,0 +1,74 @@
+"""Sentinel-2 helper utilities for CubeDynamics."""
+
+from __future__ import annotations
+
+import warnings
+
+import cubo
+import xarray as xr
+
+from . import verbs as v
+from .piping import pipe
+
+
+def load_sentinel2_ndvi_cube(
+    lat: float,
+    lon: float,
+    start: str,
+    end: str,
+    *,
+    edge_size: int = 512,
+    resolution: int = 10,
+    max_cloud: int = 40,
+    return_raw: bool = False,
+) -> xr.DataArray | tuple[xr.DataArray, xr.DataArray, xr.DataArray]:
+    """Stream Sentinel-2 L2A data and compute NDVI and NDVI z-score cubes.
+
+    Parameters
+    ----------
+    lat, lon:
+        Center point for the Sentinel-2 spatial subset.
+    start, end:
+        Date range (``YYYY-MM-DD``) to request.
+    edge_size:
+        Spatial window size in pixels (default ``512``).
+    resolution:
+        Spatial resolution in meters (default ``10``).
+    max_cloud:
+        Maximum allowed cloud cover percentage (default ``40``).
+    return_raw:
+        If ``True`` return ``(s2, ndvi, ndvi_z)``; otherwise only the NDVI
+        z-score cube is returned.
+
+    Returns
+    -------
+    xarray.DataArray or tuple of DataArray
+        NDVI z-score cube or ``(s2, ndvi, ndvi_z)`` if ``return_raw`` is true.
+    """
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        s2 = cubo.create(
+            lat=lat,
+            lon=lon,
+            collection="sentinel-2-l2a",
+            bands=["B04", "B08"],
+            start_date=start,
+            end_date=end,
+            edge_size=edge_size,
+            resolution=resolution,
+            query={"eo:cloud_cover": {"lt": max_cloud}},
+        )
+
+    if "band" in s2.dims:
+        s2 = s2.transpose("time", "y", "x", "band")
+
+    ndvi = (pipe(s2) | v.ndvi_from_s2(nir_band="B08", red_band="B04")).unwrap()
+    ndvi_z = (pipe(ndvi) | v.zscore(dim="time")).unwrap()
+
+    if return_raw:
+        return s2, ndvi, ndvi_z
+    return ndvi_z
+
+
+__all__ = ["load_sentinel2_ndvi_cube"]

tests/conftest.py

@@ -2,11 +2,22 @@
 
 from __future__ import annotations
 
+import sys
+from pathlib import Path
+
 import numpy as np
 import pandas as pd
 import pytest
 import xarray as xr
 
+# Ensure ``src`` is importable even when cubedynamics is not installed in site-packages.
+PROJECT_ROOT = Path(__file__).resolve().parents[1]
+SRC_PATH = PROJECT_ROOT / "src"
+SRC_STR = str(SRC_PATH)
+if SRC_STR in sys.path:
+    sys.path.remove(SRC_STR)
+sys.path.insert(0, SRC_STR)
+
 try:
     import dask.array as da
 except ImportError:  # pragma: no cover

tests/test_sentinel_loader_helper.py

@@ -0,0 +1,80 @@
+"""Unit tests for the Sentinel-2 NDVI helper."""
+
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+import numpy as np
+import xarray as xr
+
+PROJECT_ROOT = Path(__file__).resolve().parents[1]
+SRC_PATH = PROJECT_ROOT / "src"
+SRC_STR = str(SRC_PATH)
+if SRC_STR in sys.path:
+    sys.path.remove(SRC_STR)
+sys.path.insert(0, SRC_STR)
+
+import cubedynamics
+from cubedynamics.sentinel import load_sentinel2_ndvi_cube
+
+
+class DummyCubo:
+    """Minimal cubo shim that returns a deterministic cube."""
+
+    @staticmethod
+    def create(**kwargs):
+        rng = np.random.default_rng(0)
+        time = np.arange(3)
+        y = np.arange(4)
+        x = np.arange(5)
+        band = ["B04", "B08"]
+        data = rng.random((len(time), len(y), len(x), len(band))).astype(np.float32)
+        return xr.DataArray(
+            data,
+            coords={"time": time, "y": y, "x": x, "band": band},
+            dims=("time", "y", "x", "band"),
+            name="reflectance",
+        )
+
+
+def test_public_api_exposes_loader():
+    assert hasattr(cubedynamics, "load_sentinel2_ndvi_cube")
+
+
+def test_load_sentinel2_ndvi_cube_returns_zscores(monkeypatch):
+    import cubedynamics.sentinel as sentinel_mod
+
+    monkeypatch.setattr(sentinel_mod, "cubo", DummyCubo)
+
+    ndvi_z = load_sentinel2_ndvi_cube(
+        lat=40.0,
+        lon=-105.25,
+        start="2018-01-01",
+        end="2018-12-31",
+    )
+
+    assert ndvi_z.dims == ("time", "y", "x")
+    assert ndvi_z.name == "ndvi_zscore"
+    mean_over_time = ndvi_z.mean(dim="time")
+    assert mean_over_time.max().item() < 1e-6
+    assert mean_over_time.min().item() > -1e-6
+
+
+def test_load_sentinel2_ndvi_cube_return_raw(monkeypatch):
+    import cubedynamics.sentinel as sentinel_mod
+
+    monkeypatch.setattr(sentinel_mod, "cubo", DummyCubo)
+
+    raw_s2, ndvi, ndvi_z = load_sentinel2_ndvi_cube(
+        lat=40.0,
+        lon=-105.25,
+        start="2018-01-01",
+        end="2018-12-31",
+        return_raw=True,
+    )
+
+    assert raw_s2.dims == ("time", "y", "x", "band")
+    assert set(raw_s2.coords["band"].values) == {"B04", "B08"}
+    assert ndvi.shape == raw_s2.sel(band="B08").shape
+    assert ndvi_z.shape == ndvi.shape