From cdbef5f0a7615ef395708c1e1fc6b474e6c1e1b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9rgio=20Souza=20Costa?= Date: Wed, 30 Sep 2026 16:46:14 -0300 Subject: [PATCH 1/3] Make DisSModel optional; cube as xarray Dataset with netCDF/GeoTIFF export --- CHANGELOG.md | 18 +++ CITATION.cff | 4 +- README.md | 21 +-- disscube/api/app.py | 2 +- disscube/cli.py | 6 +- disscube/client.py | 249 ++++++++++++++++++++++----------- disscube/export.py | 148 ++++++++++++++++++++ disscube/pipeline/runner.py | 85 ++--------- docs/architecture/catalog.md | 2 +- docs/examples.md | 4 +- docs/guides/mapbiomas.md | 4 +- examples/03_time_series.py | 37 +++-- examples/README.md | 2 +- pyproject.toml | 15 +- requirements.txt | 8 +- tests/test_error_fallbacks.py | 5 +- tests/test_export.py | 175 +++++++++++++++++++++++ tests/test_pipeline_options.py | 9 +- tests/test_temporal.py | 12 +- 19 files changed, 602 insertions(+), 204 deletions(-) create mode 100644 disscube/export.py create mode 100644 tests/test_export.py diff --git a/CHANGELOG.md b/CHANGELOG.md index e0410cc..aa6bb03 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,9 +9,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Changed +- **DisSModel is now optional.** DisSCube no longer depends on it: install + `disscube[dissmodel]` only to hand a cube to a model. `CubeClient.to_lucc_data()` + is replaced by `to_raster_backend()` (same result, needs the extra; raises an + `ImportError` that says how to install it). Nothing in the pipeline or the CLI + needs DisSModel any more, including exports. The package is described as + "Declarative spatial data cubes". +- `to_raster_backend()` / `to_dataset()` raise `ValueError` when a `period` + leaves no variable at all (the old backend came back empty). - **Slimmed `examples/` and migrated real-data cases**: Real-world datasets, bundled GIS assets (~12 MB), and TerraME parity benchmarks were moved to the dedicated [LambdaGeo/disscube-recipes](https://github.com/LambdaGeo/disscube-recipes) repository (`cases/terrame_fill`, `cases/ilha_maranhao`, `cases/prodes_br163`). The core `examples/` directory now contains strictly lightweight, offline, self-contained examples with tests for both the Python API (`01_quickstart.py`, `02_vector_drivers.py`, `03_time_series.py`) and declarative TOML pipelines (`examples/pipelines/quickstart.toml`). ### Added +- `CubeClient.to_dataset()`: the cube as an `xarray.Dataset` — `(y, x)` static and + `(time, y, x)` temporal variables, CRS and transform via rioxarray — the + primary, dependency-free output. +- `CubeClient.export_netcdf()` (extra `disscube[netcdf]`): CF-1.8 netCDF with a + real `time` axis, per-variable attributes (`spec_hash`, `operator`, ...) and + `mask` kept as a variable. +- `CubeClient.export_geotiff()`: one band per variable, and per year for temporal + ones (`_`), with `VARIABLE`/`YEAR`/`SPEC_HASH` band tags. + Pipelines and `disscube run/export` write netCDF when the output ends in `.nc` + or `[export] format = "netcdf"`. - **Declarative derivations.** A `Derivation` names a source, a target grid and an operator; operators register themselves (`OPERATOR_REGISTRY`) and cover zonal statistics (`mean`, `std`, `min`, `max`, `sum`, `majority`, diff --git a/CITATION.cff b/CITATION.cff index f4926f9..d4c895b 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -4,9 +4,9 @@ authors: - family-names: "Costa" given-names: "Sérgio Souza" orcid: "https://orcid.org/0000-0002-0232-4549" -title: "DisSCube: Declarative Spatial Layer for Dynamic Models" +title: "DisSCube: Declarative spatial data cubes" abstract: "An open-source, declarative engine for constructing cellular spatial data cubes for dynamic modeling, spatial simulation, and environmental analysis." -version: 0.3.0 +version: 0.4.0 date-released: "2026-09-29" url: "https://github.com/DisSModel/disscube" repository-code: "https://github.com/DisSModel/disscube" diff --git a/README.md b/README.md index 9cb8f91..b21cd37 100644 --- a/README.md +++ b/README.md @@ -100,16 +100,19 @@ da = cube.load("forest_pct", grid_id="AC/5km") print(da.shape) # (rows, cols) ``` -### 6. Hand off to DisSModel +### 6. Get the cube out ```python -backend = cube.to_lucc_data( - ["forest_pct", "dist_roads"], - grid_id="AC/5km", - period=("2015", "2020"), -) +ds = cube.to_dataset(["forest_pct", "dist_roads"], grid_id="AC/5km", period=("2015", "2020")) +# xarray.Dataset: (y, x) static and (time, y, x) temporal variables, CRS and transform via ds.rio + +cube.export_geotiff(["forest_pct"], "forest.tif", grid_id="AC/5km") # one band per variable and year +cube.export_netcdf(["forest_pct"], "cube.nc", grid_id="AC/5km") # CF-1.8; pip install "disscube[netcdf]" ``` +DisSCube does not need DisSModel. To hand a cube to a DisSModel model, install +`disscube[dissmodel]` and use `cube.to_raster_backend(...)`, which returns a `RasterBackend`. + ## Pipeline files (TOML) A whole data preparation — grid, sources, derived variables — can be declared @@ -280,7 +283,7 @@ DISSCUBE_CATALOG=./catalog.db DISSCUBE_STORE=./data/ uvicorn disscube.api.app:ap | `GET /catalog?grid=&role=` | List derived variables | | `GET /variables/{id}` | Metadata of one derived variable, including its Zarr `asset_url` | -The API does **not** serve raster data. Models load derived variables in-process with `CubeClient.load()` / `CubeClient.to_lucc_data()`, reading the same Zarr store (local, or S3 via fsspec) that the API writes to. Interactive docs are available at `/docs` once the server is running. +The API does **not** serve raster data. Models load derived variables in-process with `CubeClient.load()` / `CubeClient.to_dataset()`, reading the same Zarr store (local, or S3 via fsspec) that the API writes to. Interactive docs are available at `/docs` once the server is running. ## Known limitations @@ -311,9 +314,9 @@ If you use DisSCube in your research, dynamic modeling, or spatial data pipeline ```bibtex @software{costa_disscube_2026, author = {Costa, S{\'e}rgio Souza}, - title = {{DisSCube: Declarative Spatial Layer for Dynamic Models}}, + title = {{DisSCube: Declarative spatial data cubes}}, year = {2026}, - version = {0.3.0}, + version = {0.4.0}, url = {https://github.com/DisSModel/disscube} } ``` diff --git a/disscube/api/app.py b/disscube/api/app.py index c3b9b52..0f908e2 100644 --- a/disscube/api/app.py +++ b/disscube/api/app.py @@ -4,7 +4,7 @@ Scope: remote *orchestration* of the catalog — register grids and sources, trigger derivations, and query what has been derived. It does **not** serve raster data: models read derived variables in-process through -``CubeClient.load()`` / ``CubeClient.to_lucc_data()``, from the same Zarr +``CubeClient.load()`` / ``CubeClient.to_dataset()``, from the same Zarr store (local or S3 via fsspec) that the API writes to. Requires the optional ``api`` extra:: diff --git a/disscube/cli.py b/disscube/cli.py index c0fa654..bc21538 100644 --- a/disscube/cli.py +++ b/disscube/cli.py @@ -29,14 +29,14 @@ def main(argv: list[str] | None = None) -> int: p_run.add_argument("file") p_run.add_argument("--workspace", help="output folder (default: the file's 'workspace', " "else a folder named after the file)") - p_run.add_argument("--output", "-o", help="export derived variables to a multi-band GeoTIFF") + p_run.add_argument("--output", "-o", help="export derived variables (GeoTIFF; netCDF if the file ends in .nc)") p_run.add_argument("--dry-run", action="store_true", help="simulate plan execution without downloading or computing") p_run.add_argument("--json", action="store_true", help="output execution report as JSON") p_run.add_argument("-v", "--verbose", action="store_true", help="log each step") - p_exp = sub.add_parser("export", help="export derived variables from an existing data cube to GeoTIFF") + p_exp = sub.add_parser("export", help="export derived variables from an existing data cube to GeoTIFF or netCDF") p_exp.add_argument("file", help="pipeline TOML file") - p_exp.add_argument("--output", "-o", required=True, help="output GeoTIFF file path (e.g. data/cellspace.tif)") + p_exp.add_argument("--output", "-o", required=True, help="output file: GeoTIFF, or netCDF if it ends in .nc (e.g. data/cellspace.tif)") p_exp.add_argument("--workspace", help="workspace folder (default: data/cube or from pipeline)") p_exp.add_argument("--variables", nargs="*", help="specific variables to export (default: all derived)") p_exp.add_argument("--json", action="store_true", help="output export result as JSON") diff --git a/disscube/client.py b/disscube/client.py index cc1dfb6..537828b 100644 --- a/disscube/client.py +++ b/disscube/client.py @@ -19,7 +19,7 @@ if TYPE_CHECKING: # Type-only imports: Derivation is imported lazily at runtime to avoid a - # circular import; dissmodel is only needed by to_lucc_data(). + # circular import; dissmodel is optional and only needed by to_raster_backend(). from dissmodel.geo.raster.backend import RasterBackend from disscube.models import Derivation @@ -27,6 +27,22 @@ log = logging.getLogger("disscube.client.cube_client") +def _filter_period(da: xr.DataArray, period: tuple[str, str]) -> xr.DataArray: + """Keep the time slices of ``da`` whose value lies within ``period`` (inclusive).""" + start, end = period + time_vals = da.coords["time"].values + start_val: int | str + end_val: int | str + if len(time_vals) > 0 and isinstance(time_vals[0], (int, np.integer)): + try: + start_val, end_val = int(start), int(end) + except ValueError: + start_val, end_val = start, end + else: + start_val, end_val = start, end + return da.isel(time=(time_vals >= start_val) & (time_vals <= end_val)) + + class CubeClient: def __init__(self, catalog: str, store: str): self.catalog = SqliteCatalogStore(catalog) @@ -221,24 +237,73 @@ def _exists(d: DerivedVariable) -> bool: derived = static[0] return xr.open_zarr(derived.asset_url, consolidated=False)[derived.name] - def to_lucc_data( + # ------------------------------------------------------------------ + # Cube output: xarray Dataset (primary), exports, DisSModel adapter + # ------------------------------------------------------------------ + + def _load_variables( self, variables: list[str], grid_id: str | None = None, period: tuple[str, str] | None = None, - ) -> RasterBackend: + ) -> dict[str, xr.DataArray]: + """Load ``variables`` as DataArrays, applying the ``period`` filter. + + Static variables are ``(y, x)``; temporal ones are ``(time, y, x)`` with an + integer-year ``time`` coordinate. A temporal variable with no slice inside + ``period`` is skipped with a warning; static variables ignore ``period``. """ - Standard integration point for the DisSModel ecosystem. - Returns a RasterBackend containing all requested variables, with - the grid's CRS and affine transform (``backend.crs``, - ``backend.transform``), so it can be written back as a GeoTIFF. + loaded: dict[str, xr.DataArray] = {} + for var_name in variables: + da = self.load(var_name, grid_id=grid_id) - Static variables are stored as (y, x) arrays — identical to the - previous behaviour; existing executors require no changes. + if da.ndim == 3 and "time" in da.dims: + if period is not None: + da = _filter_period(da, period) + if da.sizes.get("time", 0) == 0: + log.warning("%s: no data in period %s, skipped", var_name, period) + continue + da = da.transpose("time", "y", "x") + else: + da = da.transpose("y", "x") + loaded[var_name] = da - Temporal variables are stored as (time, y, x) arrays with an - explicit time axis in ``backend.time_coords``. CA models retrieve - a 2D slice via ``backend.get(name, time=step)``. + if not loaded: + raise ValueError(f"No variables could be loaded: {variables}") + return loaded + + @staticmethod + def _detect_crs(arrays): + """The CRS of the first array that declares one (``crs`` attribute or ``spatial_ref``).""" + for var_name, da in arrays.items(): + crs = da.attrs.get("crs") + if not crs and "spatial_ref" in da.coords: + try: + crs = da.rio.crs + except Exception: # noqa: BLE001 — defensive fallback: the .rio accessor raises undocumented types (e.g. ValueError, CRSError) + log.debug("Could not read CRS from %s spatial_ref", var_name) + if crs: + return crs + return None + + def to_dataset( + self, + variables: list[str], + grid_id: str | None = None, + period: tuple[str, str] | None = None, + ) -> xr.Dataset: + """ + The cube as an ``xarray.Dataset`` — the primary, dependency-free output. + + Static variables have dimensions ``(y, x)``; temporal variables + ``(time, y, x)`` with an integer-year ``time`` coordinate shared by the + whole Dataset (a variable that lacks a year present in another is NaN + there). The grid's CRS and affine transform are written with rioxarray + (``ds.rio.crs``, ``ds.rio.transform()``). Variables stay lazy (Zarr-backed). + + Auxiliary quality layers stored as extra coordinates of a variable + (e.g. purity of a ``majority``) are not carried into the Dataset; + ``load()`` still returns them. Parameters ---------- @@ -248,92 +313,106 @@ def to_lucc_data( Restrict search to a specific grid. Required when the same variable name exists on multiple grids. period : tuple[str, str] | None - Optional ``(start, end)`` filter for temporal variables. - Only time slices whose value falls within [start, end] are loaded. - Static variables are unaffected. - Example: ``period=("2000", "2014")`` - - Notes - ----- - CONTRACT decisions (open — to be resolved before 1.0): - - 1. Canonical temporal type: ``valid_from`` / ``valid_until`` accept - year strings ("2020") or ISO dates ("2020-01-01"); - ``DerivedVariable.times`` stores ``list[int]`` (years only). - Open: validate year-only format at model construction time so - callers cannot silently store wrong temporal metadata. - - 2. Missing-time behavior: a temporal variable whose slices are all - outside ``period`` is skipped with ``log.warning`` and absent from - the returned backend. The caller cannot distinguish "variable was - static (period ignored)" from "existed but outside the range". - Open: raise ``ValueError``, return a NaN slice, or keep skip. - - 3. Empty-period backend: when every requested variable is filtered - out by ``period``, ``RasterBackend`` is initialized but holds no - data arrays. The caller receives a valid-looking but empty backend. - Open: raise before returning when no variable was stored. + Optional ``(start, end)`` filter for temporal variables; only time + slices whose value falls within [start, end] are kept (e.g. + ``period=("2000", "2014")``). Static variables are unaffected. + A temporal variable with no slice in the range is skipped with a + warning; if nothing is left, ``ValueError`` is raised. + """ + arrays = self._load_variables(variables, grid_id=grid_id, period=period) + + clean = {} + for name, da in arrays.items(): + extra = [c for c in da.coords if c not in da.dims and c != "spatial_ref"] + clean[name] = da.drop_vars(extra) if extra else da + ds = xr.Dataset(clean) + + crs = self._detect_crs(arrays) + if crs: + ds = ds.rio.write_crs(crs) + transform = self._grid_transform(grid_id, variables) + if transform is not None: + ds = ds.rio.write_transform(transform) + return ds + + def export_geotiff( + self, + variables: list[str], + output: str | os.PathLike[str], + grid_id: str | None = None, + period: tuple[str, str] | None = None, + ) -> None: """ - from dissmodel.geo.raster.backend import RasterBackend + Write variables to a multi-band GeoTIFF. - detected_crs = None - backend = None + Static variables give one band named after the variable; temporal ones + one band per year, named ``_`` and ordered by year. If a + ``mask`` variable is among ``variables`` (fraction of the cell inside the + territory), cells where it is 0 are set to NaN in every other band. + """ + from disscube.export import write_geotiff - for var_name in variables: - da = self.load(var_name, grid_id=grid_id) + arrays = self._load_variables(variables, grid_id=grid_id, period=period) + transform = self._grid_transform(grid_id, variables) + write_geotiff(arrays, output, crs=self._detect_crs(arrays), transform=transform) - # CRS detection — run once - if detected_crs is None: - detected_crs = da.attrs.get("crs") - if not detected_crs and "spatial_ref" in da.coords: - try: - detected_crs = da.rio.crs - except Exception: # noqa: BLE001 — defensive fallback: the .rio accessor raises undocumented types (e.g. ValueError, CRSError) - log.debug("Could not read CRS from %s spatial_ref", var_name) + def export_netcdf( + self, + variables: list[str], + output: str | os.PathLike[str], + grid_id: str | None = None, + period: tuple[str, str] | None = None, + ) -> None: + """ + Write the cube (see :meth:`to_dataset`) to a CF-1.8 netCDF file. - if backend is None: - rows, cols = da.sizes["y"], da.sizes["x"] - backend = RasterBackend(shape=(rows, cols)) + Requires the ``netcdf`` extra (``pip install disscube[netcdf]``). The + integer-year ``time`` axis is written as ``YYYY-01-01`` dates; variables + keep their attributes (``spec_hash``, ``operator``, ``source_id``...) and + ``mask`` is kept as a variable, not applied to the data. + """ + from disscube.export import write_netcdf - if da.ndim == 3 and "time" in da.dims: - # Temporal variable — optionally filter by period - if period is not None: - start, end = period - time_vals = da.coords["time"].values - start_val: int | str - end_val: int | str - if len(time_vals) > 0 and isinstance(time_vals[0], (int, np.integer)): - try: - start_val, end_val = int(start), int(end) - except ValueError: - start_val, end_val = start, end - else: - start_val, end_val = start, end - - mask = (time_vals >= start_val) & (time_vals <= end_val) - da = da.isel(time=mask) + write_netcdf(self.to_dataset(variables, grid_id=grid_id, period=period), output) - if da.sizes.get("time", 0) == 0: - log.warning("%s: no data in period %s, skipped", var_name, period) - continue + def to_raster_backend( + self, + variables: list[str], + grid_id: str | None = None, + period: tuple[str, str] | None = None, + ) -> RasterBackend: + """ + Adapter to the DisSModel ecosystem (requires ``pip install disscube[dissmodel]``). - time_coords = da.coords["time"].values - arr = da.transpose("time", "y", "x").values - backend.set(var_name, arr, time=time_coords) + Returns a ``RasterBackend`` with the requested variables, the grid's CRS and + affine transform (``backend.crs``, ``backend.transform``). Static variables + are ``(y, x)`` arrays; temporal ones ``(time, y, x)`` with the axis in + ``backend.time_coords``, and each keeps its own time axis (no NaN alignment). + CA models get a 2D slice via ``backend.get(name, time=step)``. + ``grid_id`` and ``period`` behave as in :meth:`to_dataset`. + """ + try: + from dissmodel.geo.raster.backend import RasterBackend + except ImportError as exc: + raise ImportError( + "to_raster_backend() needs DisSModel: pip install 'disscube[dissmodel]'" + ) from exc + + arrays = self._load_variables(variables, grid_id=grid_id, period=period) + first = next(iter(arrays.values())) + backend = RasterBackend(shape=(first.sizes["y"], first.sizes["x"])) + for var_name, da in arrays.items(): + if "time" in da.dims: + backend.set(var_name, da.values, time=da.coords["time"].values) else: - # Static variable — backward-compatible path - arr = da.transpose("y", "x").values - backend.set(var_name, arr) + backend.set(var_name, da.values) - if backend is None: - raise ValueError(f"No variables could be loaded: {variables}") - - if backend.crs is None and detected_crs: - backend.crs = detected_crs + crs = self._detect_crs(arrays) + if backend.crs is None and crs: + backend.crs = crs if backend.transform is None: backend.transform = self._grid_transform(grid_id, variables) - return backend def _grid_transform(self, grid_id: str | None, variables: list[str]): diff --git a/disscube/export.py b/disscube/export.py new file mode 100644 index 0000000..641f728 --- /dev/null +++ b/disscube/export.py @@ -0,0 +1,148 @@ +"""Writers for cube output: multi-band GeoTIFF and CF netCDF. + +Both work on plain xarray objects and need no DisSModel. They are normally +reached through :meth:`CubeClient.export_geotiff` and +:meth:`CubeClient.export_netcdf`. +""" + +from __future__ import annotations + +import logging +import os +from collections.abc import Mapping +from importlib.metadata import PackageNotFoundError, version +from pathlib import Path +from typing import Any, Literal + +import numpy as np +import xarray as xr + +log = logging.getLogger("disscube.export") + +#: Above this many bands a GeoTIFF is unwieldy; netCDF keeps the time axis explicit. +MANY_BANDS = 100 + + +def software_tag() -> str: + try: + return f"DisSCube {version('disscube')}" + except PackageNotFoundError: # not installed (running from a source tree) + return "DisSCube" + + +def geotiff_bands(arrays: Mapping[str, xr.DataArray]) -> list[tuple[str, xr.DataArray, int | None]]: + """Expand variables into bands: ``(band name, 2D array, year or None)``. + + A static variable is one band named after it; a temporal one gives one band per + year, ordered by year and named ``_``. + """ + bands: list[tuple[str, xr.DataArray, int | None]] = [] + for name, da in arrays.items(): + if "time" not in da.dims: + bands.append((name, da, None)) + continue + order = np.argsort(da.coords["time"].values) + for i in order: + year = int(da.coords["time"].values[i]) + bands.append((f"{name}_{year}", da.isel(time=int(i)), year)) + return bands + + +def write_geotiff(arrays: Mapping[str, xr.DataArray], output: str | os.PathLike[str], *, + crs=None, transform=None) -> Path: + """Write ``arrays`` (``(y, x)`` or ``(time, y, x)``) to a multi-band GeoTIFF. + + If ``arrays`` has a ``mask`` variable, cells where it is 0 become NaN in every + other band (and ``mask`` itself becomes 1 inside, NaN outside). Pixels are + written as float64 with ``nodata=NaN``. Each band carries ``VARIABLE``, + ``YEAR`` (when temporal) and ``SPEC_HASH`` tags. + """ + import rasterio + + if not arrays: + raise ValueError("nothing to export: no variables") + first = next(iter(arrays.values())) + height, width = first.sizes["y"], first.sizes["x"] + + if transform is None: + try: + transform = first.rio.transform() + except Exception as exc: # rio raises several undocumented types for degenerate coords + raise ValueError("cannot export a GeoTIFF: the grid transform is unknown") from exc + if crs is None: + crs = first.rio.crs + if crs is None: + raise ValueError("cannot export a GeoTIFF: the CRS is unknown") + + mask_arr = None + if "mask" in arrays: + mask_arr = np.asarray(arrays["mask"].values, dtype=np.float64) > 0.0 + + bands = geotiff_bands(arrays) + if len(bands) > MANY_BANDS: + log.warning("%d bands in one GeoTIFF; export_netcdf() keeps the time axis explicit", len(bands)) + + out_path = Path(output) + out_path.parent.mkdir(parents=True, exist_ok=True) + with rasterio.open(out_path, "w", driver="GTiff", height=height, width=width, count=len(bands), + dtype="float64", crs=crs, transform=transform, nodata=np.nan, + compress="deflate") as dst: + dst.update_tags(TIFFTAG_SOFTWARE=software_tag(), CONVENTIONS="CF-1.8", + BANDS=",".join(name for name, _, _ in bands)) + grid_ids = {da.attrs["grid_id"] for da in arrays.values() if "grid_id" in da.attrs} + if len(grid_ids) == 1: + dst.update_tags(GRID_ID=grid_ids.pop()) + for idx, (name, band, year) in enumerate(bands, start=1): + var = name if year is None else name[: -len(str(year)) - 1] + arr = np.asarray(band.values, dtype=np.float64).copy() + if mask_arr is not None: + arr = np.where(mask_arr, 1.0, np.nan) if var == "mask" else np.where(mask_arr, arr, np.nan) + dst.write(arr, idx) + dst.set_band_description(idx, name) + tags = {"VARIABLE": var} + if year is not None: + tags["YEAR"] = str(year) + if "spec_hash" in arrays[var].attrs: + tags["SPEC_HASH"] = str(arrays[var].attrs["spec_hash"]) + dst.update_tags(idx, **tags) + return out_path + + +def _netcdf_engine() -> Literal["h5netcdf", "netcdf4"]: + try: + import h5netcdf # noqa: F401 + import h5py # noqa: F401 + return "h5netcdf" + except ImportError: + pass + try: + import netCDF4 # noqa: F401 + return "netcdf4" + except ImportError: + raise ImportError("export_netcdf() needs a netCDF backend: pip install 'disscube[netcdf]'") from None + + +def write_netcdf(ds: xr.Dataset, output: str | os.PathLike[str]) -> Path: + """Write a cube Dataset to a compressed CF-1.8 netCDF file. + + An integer-year ``time`` axis becomes ``datetime64`` (``YYYY-01-01``). + """ + engine = _netcdf_engine() + ds = ds.copy() + if "time" in ds.coords and np.issubdtype(ds["time"].dtype, np.integer): + years = ds["time"].values.astype("int64") + ds = ds.assign_coords(time=np.array([f"{y:04d}-01-01" for y in years], dtype="datetime64[ns]")) + ds["time"].attrs.update(standard_name="time", long_name="time (year of the slice)") + ds.attrs.update(Conventions="CF-1.8", source=software_tag()) + + encoding: dict[str, dict[str, Any]] = { + str(name): {"zlib": True, "complevel": 4, "_FillValue": np.nan} + for name, da in ds.data_vars.items() if np.issubdtype(da.dtype, np.floating)} + if engine == "h5netcdf": # h5netcdf spells the compression options differently + encoding = {n: {"compression": "gzip", "compression_opts": 4, "_FillValue": np.nan} + for n in encoding} + + out_path = Path(output) + out_path.parent.mkdir(parents=True, exist_ok=True) + ds.to_netcdf(out_path, engine=engine, encoding=encoding) + return out_path diff --git a/disscube/pipeline/runner.py b/disscube/pipeline/runner.py index aece1df..5836c2c 100644 --- a/disscube/pipeline/runner.py +++ b/disscube/pipeline/runner.py @@ -241,71 +241,19 @@ class ExportReport: grid_id: str -def _save_geotiff_from_backend(backend, variables: list[str], grid: GridConfig, out_path: Path) -> None: - """Directly writes a multi-band GeoTIFF using rasterio from a RasterBackend instance, +def _export(cube, variables: list[str], out_path: Path, grid_id: str, fmt: str | None = None) -> None: + """Write ``variables`` to ``out_path``. - applying the territorial boundary mask and setting pixels outside Brazil to NaN. + ``fmt`` (``"geotiff"`` or ``"netcdf"``) wins when the caller gives one; otherwise + a ``.nc`` suffix means netCDF and anything else GeoTIFF. """ - import rasterio - from rasterio.transform import from_origin - - out_path.parent.mkdir(parents=True, exist_ok=True) - crs = getattr(backend, "crs", None) or grid.crs - transform = getattr(backend, "transform", None) - shape = getattr(backend, "shape", None) - - if transform is None or shape is None: - minx, miny, maxx, maxy = grid.bbox - res = grid.resolution - width = round((maxx - minx) / res) - height = round((maxy - miny) / res) - transform = from_origin(minx, maxy, res, res) + if fmt == "netcdf" or (fmt is None and out_path.suffix.lower() in (".nc", ".nc4", ".cdf")): + cube.export_netcdf(variables, out_path, grid_id=grid_id) + log.info("exported netCDF to %s (%d variables)", out_path, len(variables)) else: - height, width = shape + cube.export_geotiff(variables, out_path, grid_id=grid_id) + log.info("exported GeoTIFF to %s (%d variables)", out_path, len(variables)) - # 1. Recupera a máscara oficial do território brasileiro (se presente) - mask_arr = None - try: - raw_mask = backend.get("mask") - if raw_mask is not None: - # Considera célula ativa qualquer uma com fração de terra > 0 - mask_arr = np.asarray(raw_mask, dtype=np.float64) > 0.0 - except (KeyError, ValueError, LookupError) as exc: - log.debug("no 'mask' variable available for export: %s", exc) - - # 2. Converte os pixels fora do Brasil para NaN - arrays = [] - for var in variables: - arr = np.asarray(backend.get(var), dtype=np.float64).copy() - if mask_arr is not None: - if var == "mask": - arr = np.where(mask_arr, 1.0, np.nan) - else: - arr[~mask_arr] = np.nan # Mar, cantos e exterior viram NaN - arrays.append(arr) - - # 3. Grava o GeoTIFF declarando nodata=np.nan - with rasterio.open( - out_path, - "w", - driver="GTiff", - height=height, - width=width, - count=len(arrays), - dtype=arrays[0].dtype, - crs=crs, - transform=transform, - nodata=np.nan, - compress="deflate", - ) as dst: - dst.update_tags( - TIFFTAG_SOFTWARE="DisSCube 0.3.0", - GRID_ID=grid.name, - CONVENTIONS="CF-1.8", - ) - for idx, (var, arr) in enumerate(zip(variables, arrays), start=1): - dst.write(arr, idx) - dst.set_band_description(idx, var) def resolve_workspace(p: Plan, workspace: str | Path | None = None) -> Path: """ @@ -377,11 +325,10 @@ def run(pipeline: PipelineFile | Plan | str | Path, workspace: str | Path | None if isinstance(export, ExportConfig) and export.variables: vars_to_export = export.variables - backend = cube.to_lucc_data(vars_to_export, grid_id=grid_id) - out_tif = Path(target_export) - _save_geotiff_from_backend(backend, vars_to_export, p.grid, out_tif) - report.exported = out_tif - log.info("exported GeoTIFF to %s (%d bands)", out_tif, len(vars_to_export)) + out_path = Path(target_export) + fmt = export.format if isinstance(export, ExportConfig) and export_geotiff is None else None + _export(cube, vars_to_export, out_path, grid_id, fmt) + report.exported = out_path record = { "pipeline": pipeline_info, @@ -405,7 +352,7 @@ def export_cube(pipeline: PipelineFile | Plan | str | Path, output: str | Path, workspace: str | Path | None = None, variables: list[str] | None = None) -> ExportReport: - """Export derived variables from an existing data cube workspace to a multi-band GeoTIFF.""" + """Export derived variables from an existing data cube workspace (GeoTIFF, or netCDF for ``.nc`` outputs).""" from disscube import CubeClient p = pipeline if isinstance(pipeline, Plan) else plan(pipeline) @@ -428,9 +375,7 @@ def export_cube(pipeline: PipelineFile | Plan | str | Path, raise PipelineError("no variables found to export (pass --variables or declare [[derive]] in pipeline)") out_path = Path(output) - backend = cube.to_lucc_data(target_vars, grid_id=grid_id) - _save_geotiff_from_backend(backend, target_vars, p.grid, out_path) - log.info("exported GeoTIFF to %s (%d bands)", out_path, len(target_vars)) + _export(cube, target_vars, out_path, grid_id) return ExportReport(workspace=ws, output=out_path, variables=target_vars, grid_id=grid_id) diff --git a/docs/architecture/catalog.md b/docs/architecture/catalog.md index bb6d49c..9560161 100644 --- a/docs/architecture/catalog.md +++ b/docs/architecture/catalog.md @@ -79,7 +79,7 @@ The `times` field of `DerivedVariable` is a list of integers (years): `CubeClient.load()` automatically detects whether there are multiple slices and stacks them into `(time, y, x)`, ordered by the first value of `times`. -`CubeClient.to_lucc_data()` accepts `period=("2000", "2020")` to keep only the slices within the interval. +`CubeClient.to_dataset()` (and `to_raster_backend()`) accept `period=("2000", "2020")` to keep only the slices within the interval. ## Querying the catalog diff --git a/docs/examples.md b/docs/examples.md index 0e21d4b..732d609 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -21,7 +21,7 @@ Examples 01–03 generate their own inputs locally, allowing every calculated nu |---|---| | [`01_quickstart.py`](https://github.com/DisSModel/disscube/blob/main/examples/01_quickstart.py) | Grid, raster sources and declarative derivations: `percentage`, `majority`, `mean`; loading results; cache hits via `spec_hash` | | [`02_vector_drivers.py`](https://github.com/DisSModel/disscube/blob/main/examples/02_vector_drivers.py) | Drivers from points, lines and polygons: `min_distance`, `count`, `presence`, `attribute`; several variables per derivation | -| [`03_time_series.py`](https://github.com/DisSModel/disscube/blob/main/examples/03_time_series.py) | Time-stamped sources, `(time, y, x)` loading, and the hand-off to DisSModel with `to_lucc_data()` (including `period`) | +| [`03_time_series.py`](https://github.com/DisSModel/disscube/blob/main/examples/03_time_series.py) | Time-stamped sources, `(time, y, x)` loading, and the hand-off to DisSModel with `to_dataset()` and `to_raster_backend()` (including `period`) | ## Declarative Pipeline Files (TOML) @@ -56,6 +56,6 @@ Real-world datasets, historical reproductions and large-scale case studies are m ## Scope -DisSCube prepares data for models; it stops at `CubeClient.to_lucc_data()`. +DisSCube prepares data for models; it stops at `CubeClient.to_dataset()` (or `to_raster_backend()` for DisSModel). Examples that run simulations with the prepared data (BR-MANGUE, LUCC) belong to the model repositories, where those dependencies live. diff --git a/docs/guides/mapbiomas.md b/docs/guides/mapbiomas.md index b980ec8..385a234 100644 --- a/docs/guides/mapbiomas.md +++ b/docs/guides/mapbiomas.md @@ -28,7 +28,7 @@ for year in (2000, 2020): ``` Each year becomes a source with `time=year`, so derivations from several years -load as a `(time, y, x)` series and go to DisSModel with `to_lucc_data()`: +load as a `(time, y, x)` series and go to DisSModel with `to_raster_backend()` (`pip install "disscube[dissmodel]"`): ```python from disscube import Derivation @@ -39,7 +39,7 @@ for year in (2000, 2020): cube.derive_declarative(Derivation(target="urban_pct", source_id=f"lulc_{year}", operator="percentage", class_code=24), grid_id=grid.id) -backend = cube.to_lucc_data(["landuse", "urban_pct"], grid_id=grid.id) +backend = cube.to_raster_backend(["landuse", "urban_pct"], grid_id=grid.id) ``` ## Code 0 is nodata diff --git a/examples/03_time_series.py b/examples/03_time_series.py index 575a63f..56848f2 100644 --- a/examples/03_time_series.py +++ b/examples/03_time_series.py @@ -1,5 +1,5 @@ """ -03 — Time series and hand-off to DisSModel. +03 — Time series, exports and the optional DisSModel hand-off. LUCC models consume a mix of time-varying variables (land use observed in several years) and static drivers. This example derives: @@ -8,9 +8,11 @@ registered as a source with its own ``time`` - ``dist_road`` a static driver from a vector layer -``load()`` stacks the temporal slices into a ``(time, y, x)`` DataArray, and -``to_lucc_data()`` packs everything into a DisSModel ``RasterBackend``, -optionally restricted to a period. +``load()`` stacks the temporal slices into a ``(time, y, x)`` DataArray; +``to_dataset()`` returns the whole cube as an ``xarray.Dataset`` (optionally +restricted to a period) and ``export_geotiff()`` writes it out. If DisSModel is +installed (``pip install "disscube[dissmodel]"``), ``to_raster_backend()`` hands +the cube to a model. python examples/03_time_series.py # temporary workspace python examples/03_time_series.py ./scratch # keep the outputs @@ -94,22 +96,33 @@ def main(workspace: Path) -> None: for year in YEARS: print(f" forest cover {year}: {float(forest.sel(time=year).mean()):.0%}") - # Hand-off to DisSModel: a RasterBackend with all requested variables. - backend = cube.to_lucc_data(["forest_pct", "dist_road"], grid_id=grid_id) - print(f"\nRasterBackend : temporal={backend.temporal_band_names()} static={backend.static_band_names()}") - print(f" forest_pct time axis: {backend.time_axis('forest_pct').tolist()}") + # The cube as an xarray Dataset: (time, y, x) temporal and (y, x) static variables. + ds = cube.to_dataset(["forest_pct", "dist_road"], grid_id=grid_id) + print(f"\nDataset : {dict(ds.sizes)} variables={list(ds.data_vars)}") # A model reads each variable by name (and year, if temporal). - loss = backend.get("forest_pct", time=2010) - backend.get("forest_pct", time=2020) - dist_km = backend.get("dist_road") / 1_000 + loss = (ds["forest_pct"].sel(time=2010) - ds["forest_pct"].sel(time=2020)).values + dist_km = ds["dist_road"].values / 1_000 print(" mean forest loss 2010→2020 by distance to the road:") for lo, hi in ((0, 3), (3, 6), (6, 12)): band = (dist_km >= lo) & (dist_km < hi) print(f" {lo:>2}–{hi:<2} km : {loss[band].mean():.0%}") # `period` keeps only the slices inside the interval. - recent = cube.to_lucc_data(["forest_pct"], grid_id=grid_id, period=("2015", "2020")) - print(f"\nperiod 2015–2020: forest_pct time axis = {recent.time_axis('forest_pct').tolist()}") + recent = cube.to_dataset(["forest_pct"], grid_id=grid_id, period=("2015", "2020")) + print(f"\nperiod 2015–2020: forest_pct time axis = {recent['time'].values.tolist()}") + + # Export: one GeoTIFF band per variable and year, or a CF netCDF with the time axis. + cube.export_geotiff(["forest_pct", "dist_road"], workspace / "cube.tif", grid_id=grid_id) + print(f"\nGeoTIFF : {workspace / 'cube.tif'} (bands forest_pct_2010, ..., dist_road)") + + # Optional hand-off to DisSModel (pip install "disscube[dissmodel]"). + try: + backend = cube.to_raster_backend(["forest_pct", "dist_road"], grid_id=grid_id) + except ImportError: + print("\nDisSModel not installed: skipping the RasterBackend hand-off") + else: + print(f"\nRasterBackend : temporal={backend.temporal_band_names()} static={backend.static_band_names()}") if __name__ == "__main__": diff --git a/examples/README.md b/examples/README.md index a675eb0..6a4c19e 100644 --- a/examples/README.md +++ b/examples/README.md @@ -132,7 +132,7 @@ disscube export examples/pipelines/quickstart.toml --output outputs/cellspace.ti |---|---| | [`01_quickstart.py`](01_quickstart.py) | Grid, raster sources and declarative derivations: `percentage`, `majority`, `mean`; loading results; cache hits via `spec_hash` | | [`02_vector_drivers.py`](02_vector_drivers.py) | Drivers from vector layers: `min_distance`, `count`, `presence`, `attribute`; multiple variables per derivation | -| [`03_time_series.py`](03_time_series.py) | Time-stamped sources, `(time, y, x)` loading, and the hand-off to DisSModel with `to_lucc_data()` | +| [`03_time_series.py`](03_time_series.py) | Time-stamped sources, `(time, y, x)` loading, and the hand-off to DisSModel with `to_dataset()` and `to_raster_backend()` | All examples are tested automatically on CI (`tests/test_examples.py`). diff --git a/pyproject.toml b/pyproject.toml index be630c6..cf41a16 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,8 +4,8 @@ build-backend = "setuptools.build_meta" [project] name = "disscube" -version = "0.3.0" -description = "Declarative Spatial Layer for Dynamic Models" +version = "0.4.0" +description = "Declarative spatial data cubes: describe sources, grid and derived variables in TOML; get a cataloged, reproducible cube." readme = "README.md" requires-python = ">=3.11" authors = [ @@ -46,7 +46,6 @@ dependencies = [ "rioxarray", "affine", "pooch>=1.8.0", - "dissmodel>=0.6.0,<0.7.0", ] [project.scripts] @@ -68,6 +67,13 @@ bdc = [ s3 = [ "s3fs", ] +netcdf = [ # CubeClient.export_netcdf() + "h5netcdf", + "h5py", +] +dissmodel = [ # CubeClient.to_raster_backend(): hand-off to DisSModel models + "dissmodel>=0.6.0,<0.7.0", +] docs = [ # documentation site (mkdocs.yml, .github/workflows/docs.yml) "mkdocs>=1.6,<2", # MkDocs 2.0 drops the plugin/theme system Material relies on "mkdocs-material>=9.5", @@ -130,6 +136,9 @@ module = [ "fiona.*", # optional [bdc] extra "pystac_client.*", # optional [bdc] extra "s3fs.*", # optional [s3] extra + "h5netcdf.*", # optional [netcdf] extra + "h5py.*", # optional [netcdf] extra + "netCDF4.*", # optional netCDF backend (export_netcdf) "fastapi.*", # optional [api] extra "uvicorn.*", # optional [api] extra ] diff --git a/requirements.txt b/requirements.txt index e018412..a827164 100644 --- a/requirements.txt +++ b/requirements.txt @@ -14,7 +14,6 @@ toml rioxarray affine pooch>=1.8.0 -dissmodel>=0.6.0,<0.7.0 # Optional: BDC fiona @@ -23,6 +22,13 @@ pystac-client # Optional: S3 s3fs +# Optional: netCDF export +h5netcdf +h5py + +# Optional: DisSModel hand-off (to_raster_backend) +dissmodel>=0.6.0,<0.7.0 + # Optional: API fastapi uvicorn diff --git a/tests/test_error_fallbacks.py b/tests/test_error_fallbacks.py index 2b47f68..7a4d6ab 100644 --- a/tests/test_error_fallbacks.py +++ b/tests/test_error_fallbacks.py @@ -57,13 +57,14 @@ def boom(*args, **kwargs): assert aligned.shape == (grid.rows, grid.cols) -def test_to_lucc_data_ignores_malformed_spatial_ref(tmp_path): +def test_to_raster_backend_ignores_malformed_spatial_ref(tmp_path): + pytest.importorskip("dissmodel") da = _raster(n=4).drop_vars("spatial_ref").assign_coords(spatial_ref=0) da.spatial_ref.attrs["crs_wkt"] = "not a WKT string" cube = CubeClient(str(tmp_path / "catalog.db"), str(tmp_path / "store")) with patch.object(CubeClient, "load", return_value=da): - backend = cube.to_lucc_data(["v"], grid_id="g") + backend = cube.to_raster_backend(["v"], grid_id="g") assert backend.crs is None np.testing.assert_array_equal(backend.get("v"), da.values) diff --git a/tests/test_export.py b/tests/test_export.py new file mode 100644 index 0000000..b7b10a9 --- /dev/null +++ b/tests/test_export.py @@ -0,0 +1,175 @@ +"""The cube as an xarray Dataset, GeoTIFF/netCDF exports, and DisSModel being optional.""" + +import sys +from pathlib import Path + +import numpy as np +import pytest +import rasterio +import xarray as xr + +from disscube import CubeClient +from disscube.models import DerivedVariable, GridSpec +from disscube.pipeline import export_cube, run + +QUICKSTART = Path(__file__).resolve().parents[1] / "examples" / "pipelines" / "quickstart.toml" + +N = 6 +CRS = "EPSG:31982" + + +def _da(values, name, **attrs): + ys = 100.0 - 10.0 * np.arange(N) - 5.0 + xs = 10.0 * np.arange(N) + 5.0 + return xr.DataArray(values.astype("float32"), dims=("y", "x"), coords={"y": ys, "x": xs}, + attrs={"crs": CRS, "grid_id": "G1", **attrs}, name=name) + + +@pytest.fixture +def cube(tmp_path): + """Static ``elev`` and ``mask``; temporal ``forest`` for 2020 and 2010 (registered out of order).""" + cube = CubeClient(str(tmp_path / "cat.db"), str(tmp_path / "store")) + cube.register_grid(GridSpec(id="G1", type="local", crs=CRS, resolution=10, bbox=[0, 40, 60, 100])) + mask = np.ones((N, N)) + mask[:, 0] = 0 + layers = [ + ("elev", [], np.arange(N * N).reshape(N, N)), + ("mask", [], mask), + ("forest", [2020], np.full((N, N), 0.2)), + ("forest", [2010], np.full((N, N), 0.8)), + ] + for name, times, values in layers: + uid = f"{name}_{'_'.join(map(str, times))}" + path = tmp_path / "store" / f"{uid}.zarr" + _da(values, name, spec_hash=f"hash-{uid}").to_dataset(name=name).to_zarr(path, mode="w", consolidated=False) + cube.catalog.save_derived(DerivedVariable( + id=uid, name=name, grid_id="G1", role="driver", times=times, dtype="float32", + derivation_id=uid, spec_hash=f"hash-{uid}", tile_id=None, content_hash=None, asset_url=str(path))) + return cube + + +def test_to_dataset_mixes_static_and_temporal(cube): + ds = cube.to_dataset(["elev", "forest"], grid_id="G1") + assert isinstance(ds, xr.Dataset) + assert ds["elev"].dims == ("y", "x") + assert ds["forest"].dims == ("time", "y", "x") + assert ds["time"].values.tolist() == [2010, 2020] + assert float(ds["forest"].sel(time=2010).mean()) == pytest.approx(0.8) + assert ds.rio.crs.to_epsg() == 31982 + t = ds.rio.transform() + assert (t.a, t.e, t.c, t.f) == (10.0, -10.0, 0.0, 100.0) + + +def test_to_dataset_period_and_empty_period(cube): + ds = cube.to_dataset(["forest", "elev"], grid_id="G1", period=("2015", "2020")) + assert ds["time"].values.tolist() == [2020] + with pytest.raises(ValueError, match="No variables"): + cube.to_dataset(["forest"], grid_id="G1", period=("1990", "2000")) + + +def test_geotiff_one_band_per_variable_and_year(cube, tmp_path): + out = tmp_path / "out" / "cube.tif" + cube.export_geotiff(["elev", "forest"], out, grid_id="G1") + with rasterio.open(out) as src: + assert src.count == 3 + assert list(src.descriptions) == ["elev", "forest_2010", "forest_2020"] + assert src.crs.to_epsg() == 31982 + assert src.transform.f == 100.0 + assert src.tags(2)["YEAR"] == "2010" + assert src.tags(2)["SPEC_HASH"] == "hash-forest_2010" + assert src.tags()["BANDS"] == "elev,forest_2010,forest_2020" + assert src.read(2)[0, 1] == pytest.approx(0.8) + + +def test_geotiff_period_selects_years(cube, tmp_path): + out = tmp_path / "recent.tif" + cube.export_geotiff(["forest"], out, grid_id="G1", period=("2015", "2020")) + with rasterio.open(out) as src: + assert list(src.descriptions) == ["forest_2020"] + + +def test_geotiff_applies_mask_when_requested(cube, tmp_path): + out = tmp_path / "masked.tif" + cube.export_geotiff(["mask", "elev"], out, grid_id="G1") + with rasterio.open(out) as src: + mask, elev = src.read(1), src.read(2) + assert np.isnan(mask[:, 0]).all() and (mask[:, 1:] == 1).all() + assert np.isnan(elev[:, 0]).all() and not np.isnan(elev[:, 1:]).any() + + +def test_netcdf_roundtrip_keeps_time_crs_and_mask(cube, tmp_path): + pytest.importorskip("h5py") + out = tmp_path / "cube.nc" + cube.export_netcdf(["elev", "forest", "mask"], out, grid_id="G1") + with xr.open_dataset(out) as ds: + assert ds.attrs["Conventions"] == "CF-1.8" + assert [str(t)[:10] for t in ds["time"].values] == ["2010-01-01", "2020-01-01"] + assert ds["forest"].dims == ("time", "y", "x") + assert float(ds["forest"].isel(time=0).mean()) == pytest.approx(0.8) + assert ds["forest"].attrs["spec_hash"].startswith("hash-forest") + assert (ds["mask"].values[:, 0] == 0).all() # kept as a variable, not applied + assert ds.rio.crs.to_epsg() == 31982 + + +def test_netcdf_without_a_backend_says_how_to_install(cube, tmp_path, monkeypatch): + monkeypatch.setitem(sys.modules, "h5netcdf", None) + monkeypatch.setitem(sys.modules, "netCDF4", None) + monkeypatch.setitem(sys.modules, "h5py", None) + with pytest.raises(ImportError, match=r"disscube\[netcdf\]"): + cube.export_netcdf(["elev"], tmp_path / "x.nc", grid_id="G1") + + +# --- DisSModel is optional --------------------------------------------------- + +def _block_dissmodel(monkeypatch): + """Make every import of dissmodel fail, even if an earlier test already loaded it.""" + for name in ("dissmodel", "dissmodel.geo", "dissmodel.geo.raster", "dissmodel.geo.raster.backend"): + monkeypatch.setitem(sys.modules, name, None) + + +def test_exports_and_dataset_work_without_dissmodel(cube, tmp_path, monkeypatch): + _block_dissmodel(monkeypatch) + cube.to_dataset(["elev", "forest"], grid_id="G1") + cube.export_geotiff(["elev"], tmp_path / "a.tif", grid_id="G1") + ws = tmp_path / "ws" + run(QUICKSTART, workspace=ws, export_geotiff=tmp_path / "q.tif") + export_cube(QUICKSTART, output=tmp_path / "q2.tif", workspace=ws) + assert (tmp_path / "q.tif").exists() and (tmp_path / "q2.tif").exists() + + +def test_to_raster_backend_without_dissmodel_says_how_to_install(cube, monkeypatch): + _block_dissmodel(monkeypatch) + with pytest.raises(ImportError, match=r"disscube\[dissmodel\]"): + cube.to_raster_backend(["elev"], grid_id="G1") + + +def test_to_raster_backend_keeps_each_time_axis(cube): + pytest.importorskip("dissmodel") + backend = cube.to_raster_backend(["elev", "forest"], grid_id="G1") + assert backend.time_axis("forest").tolist() == [2010, 2020] + assert backend.transform.f == 100.0 + + +# --- pipeline export formats ------------------------------------------------- + +def test_run_exports_netcdf_by_suffix(tmp_path): + pytest.importorskip("h5py") + out = tmp_path / "q.nc" + run(QUICKSTART, workspace=tmp_path / "ws", export_geotiff=out) + with xr.open_dataset(out) as ds: + assert len(ds.data_vars) >= 1 + + +def test_export_table_format_overrides_suffix(tmp_path): + from disscube.pipeline import runner + + calls = [] + + class Fake: + def export_netcdf(self, *a, **k): calls.append("nc") + def export_geotiff(self, *a, **k): calls.append("tif") + + runner._export(Fake(), ["v"], tmp_path / "x.dat", "g", fmt="netcdf") + runner._export(Fake(), ["v"], tmp_path / "x.nc", "g") + runner._export(Fake(), ["v"], tmp_path / "x.tif", "g") + assert calls == ["nc", "nc", "tif"] diff --git a/tests/test_pipeline_options.py b/tests/test_pipeline_options.py index 8d0344f..631e5a3 100644 --- a/tests/test_pipeline_options.py +++ b/tests/test_pipeline_options.py @@ -3,7 +3,7 @@ a script: operator ``params`` (``distance`` in another CRS, ``subcells``), the ``area`` operator, ``fill = "nearest"``, file sources with ``read`` options, ``nodata`` and a NetCDF ``variable``, ``union`` sources, and the grid's -transform on ``to_lucc_data``. +transform on ``to_raster_backend``. """ from __future__ import annotations @@ -385,14 +385,15 @@ def test_unknown_derive_param_fails_the_plan(tmp_path): # --------------------------------------------------------------------------- -# to_lucc_data +# to_raster_backend # --------------------------------------------------------------------------- -def test_to_lucc_data_carries_the_grid_transform(cube, tmp_path): +def test_to_raster_backend_carries_the_grid_transform(cube, tmp_path): + pytest.importorskip("dissmodel") _vector(cube, tmp_path, "town", [Point(-49.55, -9.45)], "EPSG:4326") _derive(cube, "geo", "town", "distance") for grid_id in ("geo", None): - backend = cube.to_lucc_data(["v"], grid_id=grid_id) + backend = cube.to_raster_backend(["v"], grid_id=grid_id) assert backend.transform == GEO.transform diff --git a/tests/test_temporal.py b/tests/test_temporal.py index 471c8c6..05486bb 100644 --- a/tests/test_temporal.py +++ b/tests/test_temporal.py @@ -1,6 +1,6 @@ """ Tests for the temporal path: writer time extraction, load() shape consistency, -and to_lucc_data period-filter logging. +and to_raster_backend period-filter logging. """ import logging @@ -132,7 +132,7 @@ def test_load_single_temporal_slice_returns_3d(tmp_path): """load() must return (time, y, x) even when only one temporal slice exists. Regression guard: previously returned 2D (y, x) for single-slice variables, - causing to_lucc_data to treat them as static. + causing to_raster_backend to treat them as static. """ cube = CubeClient(str(tmp_path / "cat.db"), str(tmp_path / "store")) grid = _grid() @@ -145,7 +145,7 @@ def test_load_single_temporal_slice_returns_3d(tmp_path): result = cube.load("v", grid_id="G1") assert result.ndim == 3, ( f"Expected 3D (time, y, x) for a single temporal slice, got {result.ndim}D. " - "to_lucc_data would silently treat this as a static variable." + "to_raster_backend would silently treat this as a static variable." ) assert "time" in result.dims assert list(result.coords["time"].values) == [2020] @@ -176,10 +176,10 @@ def test_load_multiple_temporal_slices_are_sorted(tmp_path): # --------------------------------------------------------------------------- # -# to_lucc_data: period filter emits log.warning (not print) +# to_raster_backend: period filter emits log.warning (not print) # --------------------------------------------------------------------------- # -def test_to_lucc_data_period_skip_logs_warning(tmp_path, caplog): +def test_to_raster_backend_period_skip_logs_warning(tmp_path, caplog): """A temporal variable outside the requested period logs a warning. Verifies the fix from print() to log.warning() so the message is @@ -202,7 +202,7 @@ def test_to_lucc_data_period_skip_logs_warning(tmp_path, caplog): uid="hs_s") with caplog.at_level(logging.WARNING, logger="disscube.client.cube_client"): - cube.to_lucc_data(["v", "s"], grid_id="G1", period=("2010", "2020")) + cube.to_raster_backend(["v", "s"], grid_id="G1", period=("2010", "2020")) warning_messages = [r.message for r in caplog.records if r.levelno == logging.WARNING] assert any("v" in m for m in warning_messages), ( From db6e46bddda6f3d9571bf92df70af1858ce1f7b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9rgio=20Souza=20Costa?= Date: Wed, 30 Sep 2026 16:47:01 -0300 Subject: [PATCH 2/3] Remove the experimental HTTP API DisSCube is a library and a CLI. A server will live in its own package (depending on disscube), so disscube.api, its tests, the `api` extra (fastapi, uvicorn, python-multipart) and the dev-only httpx are dropped. CI and CONTRIBUTING install the `netcdf` and `dissmodel` extras instead. --- .github/workflows/ci.yml | 5 +- CHANGELOG.md | 3 +- CONTRIBUTING.md | 2 +- README.md | 20 ------- disscube/api/__init__.py | 5 -- disscube/api/app.py | 113 --------------------------------------- pyproject.toml | 8 --- requirements.txt | 6 --- tests/test_api.py | 105 ------------------------------------ 9 files changed, 6 insertions(+), 261 deletions(-) delete mode 100644 disscube/api/__init__.py delete mode 100644 disscube/api/app.py delete mode 100644 tests/test_api.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 79728ac..ab825b4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -47,7 +47,7 @@ jobs: - name: Install dependencies run: | python -m pip install --upgrade pip - pip install -e ".[dev,bdc,api]" + pip install -e ".[dev,bdc,netcdf,dissmodel]" # Coverage XML is produced but not uploaded anywhere yet. - name: Run tests with coverage @@ -61,7 +61,8 @@ jobs: test-minimal: # Core install without optional extras: checks that tests depending on - # optional packages (fiona via disscube[bdc], fastapi via disscube[api]) + # optional packages (fiona via disscube[bdc], h5py via disscube[netcdf], DisSModel via + # disscube[dissmodel]) # are skipped, not failed. name: test (core install, no extras) runs-on: ubuntu-latest diff --git a/CHANGELOG.md b/CHANGELOG.md index aa6bb03..1e8064c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 "Declarative spatial data cubes". - `to_raster_backend()` / `to_dataset()` raise `ValueError` when a `period` leaves no variable at all (the old backend came back empty). +- **Removed the experimental HTTP API** (`disscube.api` and the `api` extra). DisSCube is + a library and a CLI; a server belongs in its own package that depends on it. - **Slimmed `examples/` and migrated real-data cases**: Real-world datasets, bundled GIS assets (~12 MB), and TerraME parity benchmarks were moved to the dedicated [LambdaGeo/disscube-recipes](https://github.com/LambdaGeo/disscube-recipes) repository (`cases/terrame_fill`, `cases/ilha_maranhao`, `cases/prodes_br163`). The core `examples/` directory now contains strictly lightweight, offline, self-contained examples with tests for both the Python API (`01_quickstart.py`, `02_vector_drivers.py`, `03_time_series.py`) and declarative TOML pipelines (`examples/pipelines/quickstart.toml`). ### Added @@ -56,7 +58,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 and `export` (write the derived variables to a multi-band GeoTIFF straight from the cube). - `CubeClient.to_lucc_data()` hands a cube to DisSModel as a raster backend. -- An experimental HTTP API, available as the optional `disscube[api]` extra. - Self-contained examples, the TerraME `Fill` correspondence with a cell-by-cell parity test suite, and the MkDocs documentation site. - `CHANGELOG.md`, `CODE_OF_CONDUCT.md`, `CITATION.cff`, and a publish diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index a50076b..ed1e7ff 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -84,7 +84,7 @@ pip install -e ".[dev]" Tests that need optional packages are skipped, not failed, when the package is missing. To run all of them, install the extras too: ```bash -pip install -e ".[dev,bdc,api]" +pip install -e ".[dev,bdc,netcdf,dissmodel]" ``` ### 4. Running the Test Suite diff --git a/README.md b/README.md index b21cd37..9e2a2c4 100644 --- a/README.md +++ b/README.md @@ -234,7 +234,6 @@ disscube/ ├── pipeline/ Pipeline execution & planning (schema, runner) + internal stages ├── catalog/ CatalogStore (Protocol) + SQLite and JSON implementations ├── storage.py AssetStore (fsspec — local and S3) -├── api/ Experimental HTTP API (optional `api` extra) ├── cli.py `disscube validate` / `disscube run` / `disscube export` ├── sources/ Adapters that bring external data in as SpatialSources, │ │ each with a checksum and a provenance.json sidecar @@ -266,25 +265,6 @@ class WeightedMeanOperator(Operator): The operator is registered automatically and accepted by `Derivation` / `SpatialDerivation` with no other change. -## Experimental: HTTP API - -`disscube.api` exposes the catalog over HTTP for **remote orchestration**: registering grids and sources, triggering derivations and querying what has been derived. It is experimental and ships as an optional extra: - -```bash -pip install -e ".[api]" -DISSCUBE_CATALOG=./catalog.db DISSCUBE_STORE=./data/ uvicorn disscube.api.app:app -``` - -| Endpoint | Purpose | -|---|---| -| `GET` / `POST /grids` | List / register `GridSpec`s | -| `GET` / `POST /sources` | List / register `SpatialSource`s | -| `POST /derive` | Run a `SpatialDerivation` (errors → HTTP 400) | -| `GET /catalog?grid=&role=` | List derived variables | -| `GET /variables/{id}` | Metadata of one derived variable, including its Zarr `asset_url` | - -The API does **not** serve raster data. Models load derived variables in-process with `CubeClient.load()` / `CubeClient.to_dataset()`, reading the same Zarr store (local, or S3 via fsspec) that the API writes to. Interactive docs are available at `/docs` once the server is running. - ## Known limitations The limitations below are scope decisions for the current version, not bugs. They are documented so that users and reviewers understand what is implemented versus what is planned. diff --git a/disscube/api/__init__.py b/disscube/api/__init__.py deleted file mode 100644 index dc00786..0000000 --- a/disscube/api/__init__.py +++ /dev/null @@ -1,5 +0,0 @@ -"""Experimental HTTP API — requires ``pip install "disscube[api]"``.""" - -from .app import app, create_app - -__all__ = ["app", "create_app"] diff --git a/disscube/api/app.py b/disscube/api/app.py deleted file mode 100644 index 0f908e2..0000000 --- a/disscube/api/app.py +++ /dev/null @@ -1,113 +0,0 @@ -""" -Experimental HTTP API for DisSCube. - -Scope: remote *orchestration* of the catalog — register grids and sources, -trigger derivations, and query what has been derived. It does **not** serve -raster data: models read derived variables in-process through -``CubeClient.load()`` / ``CubeClient.to_dataset()``, from the same Zarr -store (local or S3 via fsspec) that the API writes to. - -Requires the optional ``api`` extra:: - - pip install "disscube[api]" - uvicorn disscube.api.app:app - -The catalog and store locations come from ``DISSCUBE_CATALOG`` and -``DISSCUBE_STORE`` (see :mod:`disscube.api.config`) and are opened when the -application starts, not when this module is imported. -""" - -from __future__ import annotations - -from collections.abc import AsyncIterator -from contextlib import asynccontextmanager -from typing import Annotated - -try: - from fastapi import Depends, FastAPI, HTTPException, Request -except ImportError as exc: # pragma: no cover - exercised only without the extra - raise ImportError( - "The DisSCube HTTP API requires the optional 'api' extra: pip install \"disscube[api]\"" - ) from exc - -import os - -from disscube.client import CubeClient -from disscube.models import DerivedVariable, GridSpec, SpatialDerivation, SpatialSource - -DEFAULT_CATALOG_PATH = os.getenv("DISSCUBE_CATALOG", "./catalog.db") -DEFAULT_STORE_PATH = os.getenv("DISSCUBE_STORE", "./data/") - - -def get_cube(request: Request) -> CubeClient: - """Dependency: the CubeClient opened in the application lifespan.""" - return request.app.state.cube - - -Cube = Annotated[CubeClient, Depends(get_cube)] - - -def create_app(catalog_path: str | None = None, store_path: str | None = None) -> FastAPI: - """ - Build the API application. - - Parameters default to ``DEFAULT_CATALOG_PATH`` / ``DEFAULT_STORE_PATH``. - The ``CubeClient`` is created in the application lifespan, so building - the app has no side effects on disk. - """ - - @asynccontextmanager - async def lifespan(app: FastAPI) -> AsyncIterator[None]: - app.state.cube = CubeClient( - catalog_path or DEFAULT_CATALOG_PATH, - store_path or DEFAULT_STORE_PATH, - ) - yield - - app = FastAPI( - title="DisSCube API (experimental)", - description="Catalog orchestration for DisSCube. Does not serve raster data.", - lifespan=lifespan, - ) - - @app.get("/grids", response_model=list[GridSpec]) - def list_grids(cube: Cube): - return cube.catalog.list_grids() - - @app.post("/grids") - def register_grid(grid: GridSpec, cube: Cube): - cube.register_grid(grid) - return {"status": "ok"} - - @app.get("/sources", response_model=list[SpatialSource]) - def list_spatial_sources(cube: Cube): - return cube.catalog.list_spatial_sources() - - @app.post("/sources") - def register_spatial_source(source: SpatialSource, cube: Cube): - cube.register_spatial_source(source) - return {"status": "ok"} - - @app.post("/derive", response_model=list[DerivedVariable]) - def derive(derivation: SpatialDerivation, cube: Cube): - try: - return cube.derive(derivation) - except Exception as e: # API boundary: report any derivation failure as a 400 - raise HTTPException(status_code=400, detail=str(e)) from e - - @app.get("/catalog", response_model=list[DerivedVariable]) - def get_catalog(cube: Cube, grid: str | None = None, role: str | None = None): - return cube.search(grid=grid, role=role) - - @app.get("/variables/{variable_id}", response_model=DerivedVariable) - def get_variable(variable_id: str, cube: Cube): - """Metadata of a derived variable (including its Zarr ``asset_url``), not its data.""" - for d in cube.search(): - if d.id == variable_id or d.name == variable_id: - return d - raise HTTPException(status_code=404, detail="Variable not found") - - return app - - -app = create_app() diff --git a/pyproject.toml b/pyproject.toml index cf41a16..567d243 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,7 +57,6 @@ dev = [ "pytest-cov", "mypy>=1.13", "ruff>=0.16,<0.17", - "httpx", "ipykernel", ] bdc = [ @@ -78,11 +77,6 @@ docs = [ # documentation site (mkdocs.yml, .github/workflows/ "mkdocs>=1.6,<2", # MkDocs 2.0 drops the plugin/theme system Material relies on "mkdocs-material>=9.5", ] -api = [ # experimental HTTP API (disscube.api) - "fastapi", - "uvicorn", - "python-multipart", -] [project.urls] Homepage = "https://github.com/DisSModel/disscube" @@ -139,8 +133,6 @@ module = [ "h5netcdf.*", # optional [netcdf] extra "h5py.*", # optional [netcdf] extra "netCDF4.*", # optional netCDF backend (export_netcdf) - "fastapi.*", # optional [api] extra - "uvicorn.*", # optional [api] extra ] ignore_missing_imports = true diff --git a/requirements.txt b/requirements.txt index a827164..e401613 100644 --- a/requirements.txt +++ b/requirements.txt @@ -29,17 +29,11 @@ h5py # Optional: DisSModel hand-off (to_raster_backend) dissmodel>=0.6.0,<0.7.0 -# Optional: API -fastapi -uvicorn -python-multipart - # Development & Testing pytest pytest-cov mypy>=1.13 ruff>=0.16,<0.17 -httpx ipykernel # Documentation diff --git a/tests/test_api.py b/tests/test_api.py deleted file mode 100644 index a51c37a..0000000 --- a/tests/test_api.py +++ /dev/null @@ -1,105 +0,0 @@ -""" -Tests for the experimental HTTP API (disscube.api). - -Skipped when the optional ``api`` extra (fastapi) is not installed. -""" - -import numpy as np -import pytest -import rasterio -from rasterio.transform import from_origin - -pytest.importorskip("fastapi") -from fastapi.testclient import TestClient - -from disscube.api import create_app - -GRID = { - "id": "T/100m", - "type": "local", - "crs": "EPSG:31983", - "resolution": 100.0, - "bbox": [0.0, 0.0, 400.0, 400.0], -} - - -@pytest.fixture -def client(tmp_path): - app = create_app(str(tmp_path / "catalog.db"), str(tmp_path / "store")) - with TestClient(app) as c: # the context manager runs the lifespan - yield c - - -@pytest.fixture -def raster_path(tmp_path): - """8x8 GeoTIFF at 50 m covering GRID's bbox; values 0..63.""" - path = tmp_path / "src.tif" - data = np.arange(64, dtype="float32").reshape(8, 8) - with rasterio.open( - path, "w", driver="GTiff", height=8, width=8, count=1, dtype="float32", - crs="EPSG:31983", transform=from_origin(0.0, 400.0, 50.0, 50.0), - ) as dst: - dst.write(data, 1) - return path - - -def test_create_app_has_no_side_effects(tmp_path): - catalog = tmp_path / "catalog.db" - create_app(str(catalog), str(tmp_path / "store")) - assert not catalog.exists() # the catalog is opened only at startup - - -def test_register_and_list_grids(client): - assert client.get("/grids").json() == [] - assert client.post("/grids", json=GRID).json() == {"status": "ok"} - grids = client.get("/grids").json() - assert [g["id"] for g in grids] == ["T/100m"] - - -def test_invalid_payload_is_rejected(client): - r = client.post("/grids", json={"id": "incomplete"}) - assert r.status_code == 422 - - -def test_unknown_variable_returns_404(client): - assert client.get("/variables/nope").status_code == 404 - - -def test_derive_failure_returns_400(client): - client.post("/grids", json=GRID) - r = client.post("/derive", json={ - "source_id": "missing_source", - "grid_id": "T/100m", - "role": "driver", - "variables": [{"name": "v", "operator": "mean"}], - }) - assert r.status_code == 400 - - -def test_derive_end_to_end(client, raster_path): - client.post("/grids", json=GRID) - client.post("/sources", json={ - "id": "src", - "name": "synthetic raster", - "format": "raster", - "asset_url": str(raster_path), - "crs": "EPSG:31983", - }) - assert [s["id"] for s in client.get("/sources").json()] == ["src"] - - r = client.post("/derive", json={ - "source_id": "src", - "grid_id": "T/100m", - "role": "driver", - "variables": [{"name": "elev", "operator": "mean"}], - }) - assert r.status_code == 200, r.text - derived = r.json() - assert [d["name"] for d in derived] == ["elev"] - - catalog = client.get("/catalog", params={"grid": "T/100m"}).json() - assert [d["name"] for d in catalog] == ["elev"] - - meta = client.get("/variables/elev").json() - assert meta["spec_hash"] == derived[0]["spec_hash"] - assert meta["asset_url"].endswith("elev.zarr") From b4e9d17a7ffa6d1032a0c30819b00e874546bd29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9rgio=20Souza=20Costa?= Date: Wed, 30 Sep 2026 16:47:58 -0300 Subject: [PATCH 3/3] Record provenance in exports - VariableWriter stores source_checksum with each new variable. - CubeClient.provenance(): spec_hash, content_hash, source_id and source_checksum per time slice; load() now shares its slice selection with it (_select). - netCDF: global `disscube_provenance` (JSON), plain attributes on single-slice variables, `history`. A temporal variable no longer exposes the first slice's spec_hash as if it were the whole variable's. - GeoTIFF: per-band SPEC_HASH / CONTENT_HASH / SOURCE_CHECKSUM / SOURCE_ID; the band of each year now gets that year's hash, not the first one's. - Pipeline exports record pipeline_file, pipeline_checksum, pipeline_name. --- CHANGELOG.md | 8 +++ README.md | 4 ++ disscube/client.py | 102 ++++++++++++++++++++++++++++++------ disscube/export.py | 36 ++++++++++--- disscube/pipeline/runner.py | 19 +++++-- disscube/pipeline/writer.py | 2 + tests/test_export.py | 49 +++++++++++++++-- 7 files changed, 188 insertions(+), 32 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1e8064c..58a60cb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - `CubeClient.export_netcdf()` (extra `disscube[netcdf]`): CF-1.8 netCDF with a real `time` axis, per-variable attributes (`spec_hash`, `operator`, ...) and `mask` kept as a variable. +- **Provenance in the exports.** The writer now stores `source_checksum` (the input the + slice came from) with each new variable, and `CubeClient.provenance()` lists, per time + slice, `spec_hash`, `content_hash`, `source_id` and `source_checksum`. netCDF files carry + it as the global JSON attribute `disscube_provenance` (and as plain attributes on + single-slice variables), plus `history`; GeoTIFF bands carry `SPEC_HASH`, `CONTENT_HASH`, + `SOURCE_CHECKSUM` and `SOURCE_ID` tags. Exports made by a pipeline also record + `pipeline_file`, `pipeline_checksum` and `pipeline_name`. Variables derived by earlier + versions lack `source_checksum` until they are derived again. - `CubeClient.export_geotiff()`: one band per variable, and per year for temporal ones (`_`), with `VARIABLE`/`YEAR`/`SPEC_HASH` band tags. Pipelines and `disscube run/export` write netCDF when the output ends in `.nc` diff --git a/README.md b/README.md index 9e2a2c4..a0255a5 100644 --- a/README.md +++ b/README.md @@ -110,6 +110,10 @@ cube.export_geotiff(["forest_pct"], "forest.tif", grid_id="AC/5km") # one band cube.export_netcdf(["forest_pct"], "cube.nc", grid_id="AC/5km") # CF-1.8; pip install "disscube[netcdf]" ``` +Exports carry their provenance: each GeoTIFF band and netCDF variable records the +`spec_hash`, the `content_hash` of the stored data and the `source_checksum` of the +input it came from (`cube.provenance("forest_pct")` lists them per year). + DisSCube does not need DisSModel. To hand a cube to a DisSModel model, install `disscube[dissmodel]` and use `cube.to_raster_backend(...)`, which returns a `RasterBackend`. diff --git a/disscube/client.py b/disscube/client.py index 537828b..e35a084 100644 --- a/disscube/client.py +++ b/disscube/client.py @@ -1,9 +1,11 @@ from __future__ import annotations +import json import logging import os import sys -from typing import TYPE_CHECKING +from collections.abc import Mapping +from typing import TYPE_CHECKING, Any import numpy as np import xarray as xr @@ -26,6 +28,9 @@ log = logging.getLogger("disscube.client.cube_client") +# Attributes that describe one stored slice: per-variable only when it has a single slice. +_SLICE_ATTRS = ("spec_hash", "content_hash", "source_id", "source_checksum") + def _filter_period(da: xr.DataArray, period: tuple[str, str]) -> xr.DataArray: """Keep the time slices of ``da`` whose value lies within ``period`` (inclusive).""" @@ -178,6 +183,23 @@ def load( For temporal variables (DerivedVariable.times is non-empty), returns a DataArray with dims (time, y, x). For static variables, returns (y, x). """ + selected = self._select(variable_id, tile_id, grid_id) + if selected[0].times: + # Stack temporal slices along the time axis + slices = [] + time_coords: list[int] = [] + for d in selected: + slices.append(xr.open_zarr(d.asset_url, consolidated=False)[d.name]) + time_coords.extend(d.times) + return xr.concat(slices, dim=xr.DataArray(time_coords, dims="time")) + + # Static + derived = selected[0] + return xr.open_zarr(derived.asset_url, consolidated=False)[derived.name] + + def _select(self, variable_id: str, tile_id: str | None, grid_id: str | None) -> list[DerivedVariable]: + """The catalog entries that make up a variable: its temporal slices sorted by + time, or its single static entry. Entries whose files are gone are skipped.""" matches = [] for d in self.catalog.search_derived_variables(tile_id=tile_id, grid_id=grid_id): if d.id == variable_id or d.name == variable_id: @@ -218,24 +240,49 @@ def _exists(d: DerivedVariable) -> bool: static = [d for d in matches if not d.times and _exists(d)] if temporal: - # Stack temporal slices along time axis sorted by first time value - temporal_sorted = sorted(temporal, key=lambda d: d.times[0]) - slices = [] - time_coords = [] - for d in temporal_sorted: - da = xr.open_zarr(d.asset_url, consolidated=False)[d.name] - slices.append(da) - time_coords.extend(d.times) - return xr.concat(slices, dim=xr.DataArray(time_coords, dims="time")) + return sorted(temporal, key=lambda d: d.times[0]) - # Static — original behaviour if not static: msg = f"Derived variable not found on disk: {variable_id}" if grid_id: msg += f" on grid {grid_id}" raise ValueError(msg) - derived = static[0] - return xr.open_zarr(derived.asset_url, consolidated=False)[derived.name] + return [static[0]] + + def provenance(self, variable_id: str, grid_id: str | None = None, + tile_id: str | None = None) -> list[dict[str, Any]]: + """Where each slice of a variable came from: one record per time slice (a single + record for a static variable), ordered by time. + + Each record has ``time`` (year or ``None``), ``spec_hash`` and ``content_hash`` + (SHA-256 of the stored Zarr), plus ``source_id`` and ``source_checksum`` (the + checksum of the input the slice was derived from) when they were recorded. + Variables derived before ``source_checksum`` was stored simply lack that key. + """ + records: list[dict[str, Any]] = [] + for d in self._select(variable_id, tile_id, grid_id): + attrs = xr.open_zarr(d.asset_url, consolidated=False)[d.name].attrs + record: dict[str, Any] = { + "time": d.times[0] if d.times else None, + "spec_hash": d.spec_hash, + "content_hash": d.content_hash, + "source_id": attrs.get("source_id"), + "source_checksum": attrs.get("source_checksum"), + } + records.append({k: v for k, v in record.items() if v is not None or k == "time"}) + return records + + def _provenance_of(self, arrays: Mapping[str, xr.DataArray], + grid_id: str | None) -> dict[str, list[dict[str, Any]]]: + """Provenance records for loaded arrays, keeping only the slices they hold.""" + out: dict[str, list[dict[str, Any]]] = {} + for name, da in arrays.items(): + records = self.provenance(name, grid_id=grid_id) + if "time" in da.dims: + years = {int(t) for t in da.coords["time"].values} + records = [r for r in records if r["time"] in years] + out[name] = records + return out # ------------------------------------------------------------------ # Cube output: xarray Dataset (primary), exports, DisSModel adapter @@ -321,11 +368,21 @@ def to_dataset( """ arrays = self._load_variables(variables, grid_id=grid_id, period=period) + provenance = self._provenance_of(arrays, grid_id) clean = {} for name, da in arrays.items(): extra = [c for c in da.coords if c not in da.dims and c != "spatial_ref"] - clean[name] = da.drop_vars(extra) if extra else da + da = da.drop_vars(extra) if extra else da + attrs = {k: v for k, v in da.attrs.items() if k not in _SLICE_ATTRS} + if "time" not in da.dims and provenance[name]: + # One slice: its provenance is the variable's. A temporal variable + # has one per year, so it lives only in ``disscube_provenance``. + attrs.update({k: v for k, v in provenance[name][0].items() if k in _SLICE_ATTRS}) + da = da.copy(deep=False) + da.attrs = attrs + clean[name] = da ds = xr.Dataset(clean) + ds.attrs["disscube_provenance"] = json.dumps({"variables": provenance}, sort_keys=True) crs = self._detect_crs(arrays) if crs: @@ -341,6 +398,7 @@ def export_geotiff( output: str | os.PathLike[str], grid_id: str | None = None, period: tuple[str, str] | None = None, + attrs: Mapping[str, str] | None = None, ) -> None: """ Write variables to a multi-band GeoTIFF. @@ -349,12 +407,17 @@ def export_geotiff( one band per year, named ``_`` and ordered by year. If a ``mask`` variable is among ``variables`` (fraction of the cell inside the territory), cells where it is 0 are set to NaN in every other band. + + Each band is tagged with the provenance of its slice (``SPEC_HASH``, + ``CONTENT_HASH``, ``SOURCE_CHECKSUM``, ``SOURCE_ID``). ``attrs`` adds + file-level tags (keys upper-cased), e.g. the pipeline that produced the cube. """ from disscube.export import write_geotiff arrays = self._load_variables(variables, grid_id=grid_id, period=period) transform = self._grid_transform(grid_id, variables) - write_geotiff(arrays, output, crs=self._detect_crs(arrays), transform=transform) + write_geotiff(arrays, output, crs=self._detect_crs(arrays), transform=transform, + provenance=self._provenance_of(arrays, grid_id), attrs=attrs) def export_netcdf( self, @@ -362,6 +425,7 @@ def export_netcdf( output: str | os.PathLike[str], grid_id: str | None = None, period: tuple[str, str] | None = None, + attrs: Mapping[str, str] | None = None, ) -> None: """ Write the cube (see :meth:`to_dataset`) to a CF-1.8 netCDF file. @@ -370,10 +434,16 @@ def export_netcdf( integer-year ``time`` axis is written as ``YYYY-01-01`` dates; variables keep their attributes (``spec_hash``, ``operator``, ``source_id``...) and ``mask`` is kept as a variable, not applied to the data. + + Provenance: a variable with a single slice carries ``spec_hash``, + ``content_hash``, ``source_id`` and ``source_checksum`` as attributes; for + every variable the global attribute ``disscube_provenance`` holds them per + slice as JSON. ``history`` records the export; ``attrs`` adds global + attributes (e.g. the pipeline file and its checksum). """ from disscube.export import write_netcdf - write_netcdf(self.to_dataset(variables, grid_id=grid_id, period=period), output) + write_netcdf(self.to_dataset(variables, grid_id=grid_id, period=period), output, attrs=attrs) def to_raster_backend( self, diff --git a/disscube/export.py b/disscube/export.py index 641f728..85d633f 100644 --- a/disscube/export.py +++ b/disscube/export.py @@ -10,6 +10,7 @@ import logging import os from collections.abc import Mapping +from datetime import UTC, datetime from importlib.metadata import PackageNotFoundError, version from pathlib import Path from typing import Any, Literal @@ -48,14 +49,22 @@ def geotiff_bands(arrays: Mapping[str, xr.DataArray]) -> list[tuple[str, xr.Data return bands +_BAND_TAGS = (("spec_hash", "SPEC_HASH"), ("content_hash", "CONTENT_HASH"), + ("source_checksum", "SOURCE_CHECKSUM"), ("source_id", "SOURCE_ID")) + + def write_geotiff(arrays: Mapping[str, xr.DataArray], output: str | os.PathLike[str], *, - crs=None, transform=None) -> Path: + crs=None, transform=None, + provenance: Mapping[str, list[dict[str, Any]]] | None = None, + attrs: Mapping[str, str] | None = None) -> Path: """Write ``arrays`` (``(y, x)`` or ``(time, y, x)``) to a multi-band GeoTIFF. If ``arrays`` has a ``mask`` variable, cells where it is 0 become NaN in every other band (and ``mask`` itself becomes 1 inside, NaN outside). Pixels are - written as float64 with ``nodata=NaN``. Each band carries ``VARIABLE``, - ``YEAR`` (when temporal) and ``SPEC_HASH`` tags. + written as float64 with ``nodata=NaN``. Each band carries ``VARIABLE`` and ``YEAR`` + (when temporal) tags, plus the ``SPEC_HASH``, ``CONTENT_HASH``, ``SOURCE_CHECKSUM`` + and ``SOURCE_ID`` of its slice found in ``provenance`` (``{variable: [records]}``, + as returned by ``CubeClient.provenance``). ``attrs`` become file-level tags. """ import rasterio @@ -88,7 +97,10 @@ def write_geotiff(arrays: Mapping[str, xr.DataArray], output: str | os.PathLike[ dtype="float64", crs=crs, transform=transform, nodata=np.nan, compress="deflate") as dst: dst.update_tags(TIFFTAG_SOFTWARE=software_tag(), CONVENTIONS="CF-1.8", + TIFFTAG_DATETIME=datetime.now(UTC).strftime("%Y:%m:%d %H:%M:%S"), BANDS=",".join(name for name, _, _ in bands)) + if attrs: + dst.update_tags(**{str(k).upper(): str(v) for k, v in attrs.items()}) grid_ids = {da.attrs["grid_id"] for da in arrays.values() if "grid_id" in da.attrs} if len(grid_ids) == 1: dst.update_tags(GRID_ID=grid_ids.pop()) @@ -102,8 +114,10 @@ def write_geotiff(arrays: Mapping[str, xr.DataArray], output: str | os.PathLike[ tags = {"VARIABLE": var} if year is not None: tags["YEAR"] = str(year) - if "spec_hash" in arrays[var].attrs: - tags["SPEC_HASH"] = str(arrays[var].attrs["spec_hash"]) + record = next((r for r in (provenance or {}).get(var, []) if r.get("time") == year), {}) + for key, tag in _BAND_TAGS: + if record.get(key) is not None: + tags[tag] = str(record[key]) dst.update_tags(idx, **tags) return out_path @@ -122,10 +136,13 @@ def _netcdf_engine() -> Literal["h5netcdf", "netcdf4"]: raise ImportError("export_netcdf() needs a netCDF backend: pip install 'disscube[netcdf]'") from None -def write_netcdf(ds: xr.Dataset, output: str | os.PathLike[str]) -> Path: +def write_netcdf(ds: xr.Dataset, output: str | os.PathLike[str], *, + attrs: Mapping[str, str] | None = None) -> Path: """Write a cube Dataset to a compressed CF-1.8 netCDF file. - An integer-year ``time`` axis becomes ``datetime64`` (``YYYY-01-01``). + An integer-year ``time`` axis becomes ``datetime64`` (``YYYY-01-01``). Global + attributes get ``Conventions``, ``source`` (the DisSCube version), ``history`` + (when it was written) and whatever is in ``attrs``. """ engine = _netcdf_engine() ds = ds.copy() @@ -133,7 +150,10 @@ def write_netcdf(ds: xr.Dataset, output: str | os.PathLike[str]) -> Path: years = ds["time"].values.astype("int64") ds = ds.assign_coords(time=np.array([f"{y:04d}-01-01" for y in years], dtype="datetime64[ns]")) ds["time"].attrs.update(standard_name="time", long_name="time (year of the slice)") - ds.attrs.update(Conventions="CF-1.8", source=software_tag()) + ds.attrs.update(Conventions="CF-1.8", source=software_tag(), + history=f"{datetime.now(UTC).isoformat(timespec='seconds')} {software_tag()}: exported") + if attrs: + ds.attrs.update({str(k): str(v) for k, v in attrs.items()}) encoding: dict[str, dict[str, Any]] = { str(name): {"zlib": True, "complevel": 4, "_FillValue": np.nan} diff --git a/disscube/pipeline/runner.py b/disscube/pipeline/runner.py index 5836c2c..339173f 100644 --- a/disscube/pipeline/runner.py +++ b/disscube/pipeline/runner.py @@ -241,17 +241,26 @@ class ExportReport: grid_id: str -def _export(cube, variables: list[str], out_path: Path, grid_id: str, fmt: str | None = None) -> None: +def _pipeline_attrs(pf: PipelineFile, name: str | None) -> dict[str, str]: + """File-level provenance for an export: which pipeline file (and content) made the cube.""" + attrs = {"pipeline_file": pf.path.name, "pipeline_checksum": pf.checksum} + if name: + attrs["pipeline_name"] = name + return attrs + + +def _export(cube, variables: list[str], out_path: Path, grid_id: str, fmt: str | None = None, + attrs: dict[str, str] | None = None) -> None: """Write ``variables`` to ``out_path``. ``fmt`` (``"geotiff"`` or ``"netcdf"``) wins when the caller gives one; otherwise a ``.nc`` suffix means netCDF and anything else GeoTIFF. """ if fmt == "netcdf" or (fmt is None and out_path.suffix.lower() in (".nc", ".nc4", ".cdf")): - cube.export_netcdf(variables, out_path, grid_id=grid_id) + cube.export_netcdf(variables, out_path, grid_id=grid_id, attrs=attrs) log.info("exported netCDF to %s (%d variables)", out_path, len(variables)) else: - cube.export_geotiff(variables, out_path, grid_id=grid_id) + cube.export_geotiff(variables, out_path, grid_id=grid_id, attrs=attrs) log.info("exported GeoTIFF to %s (%d variables)", out_path, len(variables)) @@ -327,7 +336,7 @@ def run(pipeline: PipelineFile | Plan | str | Path, workspace: str | Path | None out_path = Path(target_export) fmt = export.format if isinstance(export, ExportConfig) and export_geotiff is None else None - _export(cube, vars_to_export, out_path, grid_id, fmt) + _export(cube, vars_to_export, out_path, grid_id, fmt, _pipeline_attrs(pf, cfg.name)) report.exported = out_path record = { @@ -375,7 +384,7 @@ def export_cube(pipeline: PipelineFile | Plan | str | Path, raise PipelineError("no variables found to export (pass --variables or declare [[derive]] in pipeline)") out_path = Path(output) - _export(cube, target_vars, out_path, grid_id) + _export(cube, target_vars, out_path, grid_id, attrs=_pipeline_attrs(p.file, cfg.name)) return ExportReport(workspace=ws, output=out_path, variables=target_vars, grid_id=grid_id) diff --git a/disscube/pipeline/writer.py b/disscube/pipeline/writer.py index 9f9e7bf..7101b2f 100644 --- a/disscube/pipeline/writer.py +++ b/disscube/pipeline/writer.py @@ -49,6 +49,8 @@ def execute(self, ctx: PipelineContext) -> PipelineContext: da.attrs["operator"] = str(op_name) if getattr(derivation, "source_id", None): da.attrs["source_id"] = derivation.source_id + if getattr(derivation, "source_checksum", None): + da.attrs["source_checksum"] = derivation.source_checksum if tile_id: da.attrs["tile_id"] = tile_id if "spatial_ref" in da.coords: diff --git a/tests/test_export.py b/tests/test_export.py index b7b10a9..be6ece9 100644 --- a/tests/test_export.py +++ b/tests/test_export.py @@ -1,5 +1,6 @@ """The cube as an xarray Dataset, GeoTIFF/netCDF exports, and DisSModel being optional.""" +import json import sys from pathlib import Path @@ -41,10 +42,10 @@ def cube(tmp_path): for name, times, values in layers: uid = f"{name}_{'_'.join(map(str, times))}" path = tmp_path / "store" / f"{uid}.zarr" - _da(values, name, spec_hash=f"hash-{uid}").to_dataset(name=name).to_zarr(path, mode="w", consolidated=False) + _da(values, name, spec_hash=f"hash-{uid}", source_id=f"src-{uid}", source_checksum=f"sha256:src-{uid}").to_dataset(name=name).to_zarr(path, mode="w", consolidated=False) cube.catalog.save_derived(DerivedVariable( id=uid, name=name, grid_id="G1", role="driver", times=times, dtype="float32", - derivation_id=uid, spec_hash=f"hash-{uid}", tile_id=None, content_hash=None, asset_url=str(path))) + derivation_id=uid, spec_hash=f"hash-{uid}", tile_id=None, content_hash=f"sha256:zarr-{uid}", asset_url=str(path))) return cube @@ -77,6 +78,10 @@ def test_geotiff_one_band_per_variable_and_year(cube, tmp_path): assert src.transform.f == 100.0 assert src.tags(2)["YEAR"] == "2010" assert src.tags(2)["SPEC_HASH"] == "hash-forest_2010" + assert src.tags(3)["SPEC_HASH"] == "hash-forest_2020" # each year its own slice + assert src.tags(3)["SOURCE_CHECKSUM"] == "sha256:src-forest_2020" + assert src.tags(3)["CONTENT_HASH"] == "sha256:zarr-forest_2020" + assert src.tags(1)["SOURCE_ID"] == "src-elev_" assert src.tags()["BANDS"] == "elev,forest_2010,forest_2020" assert src.read(2)[0, 1] == pytest.approx(0.8) @@ -106,7 +111,14 @@ def test_netcdf_roundtrip_keeps_time_crs_and_mask(cube, tmp_path): assert [str(t)[:10] for t in ds["time"].values] == ["2010-01-01", "2020-01-01"] assert ds["forest"].dims == ("time", "y", "x") assert float(ds["forest"].isel(time=0).mean()) == pytest.approx(0.8) - assert ds["forest"].attrs["spec_hash"].startswith("hash-forest") + assert "spec_hash" not in ds["forest"].attrs # one per year: see provenance + assert ds["elev"].attrs["spec_hash"] == "hash-elev_" + assert ds["elev"].attrs["source_checksum"] == "sha256:src-elev_" + prov = json.loads(ds.attrs["disscube_provenance"])["variables"] + assert [(r["time"], r["source_checksum"]) for r in prov["forest"]] == [ + (2010, "sha256:src-forest_2010"), (2020, "sha256:src-forest_2020")] + assert prov["forest"][1]["content_hash"] == "sha256:zarr-forest_2020" + assert "exported" in ds.attrs["history"] and ds.attrs["source"].startswith("DisSCube") assert (ds["mask"].values[:, 0] == 0).all() # kept as a variable, not applied assert ds.rio.crs.to_epsg() == 31982 @@ -173,3 +185,34 @@ def export_geotiff(self, *a, **k): calls.append("tif") runner._export(Fake(), ["v"], tmp_path / "x.nc", "g") runner._export(Fake(), ["v"], tmp_path / "x.tif", "g") assert calls == ["nc", "nc", "tif"] + + +def test_provenance_lists_slices_in_time_order(cube): + records = cube.provenance("forest", grid_id="G1") + assert [r["time"] for r in records] == [2010, 2020] + assert records[0]["spec_hash"] == "hash-forest_2010" + (static,) = cube.provenance("elev", grid_id="G1") + assert static["time"] is None and static["source_id"] == "src-elev_" + + +def test_provenance_tolerates_products_without_source_checksum(cube, tmp_path): + # a zarr written before source_checksum was recorded + path = tmp_path / "store" / "old.zarr" + da = _da(np.ones((N, N)), "old", spec_hash="h-old") + da.to_dataset(name="old").to_zarr(path, mode="w", consolidated=False) + cube.catalog.save_derived(DerivedVariable( + id="old", name="old", grid_id="G1", role="driver", times=[], dtype="float32", + derivation_id="h-old", spec_hash="h-old", tile_id=None, content_hash=None, asset_url=str(path))) + (rec,) = cube.provenance("old", grid_id="G1") + assert rec == {"time": None, "spec_hash": "h-old"} + + +def test_pipeline_export_records_the_pipeline_and_source_checksums(tmp_path): + out = tmp_path / "q.tif" + run(QUICKSTART, workspace=tmp_path / "ws", export_geotiff=out) + with rasterio.open(out) as src: + tags = src.tags() + assert tags["PIPELINE_FILE"] == "quickstart.toml" + assert tags["PIPELINE_CHECKSUM"].startswith("sha256:") + assert src.tags(1)["SOURCE_CHECKSUM"].startswith("sha256:") + assert len(src.tags(1)["CONTENT_HASH"]) >= 32