diff --git a/CHANGELOG.md b/CHANGELOG.md index f45f98f..e0410cc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,10 +8,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] -First public release (`0.2.0`). DisSCube is a declarative spatial layer for -dynamic models: it aligns heterogeneous geospatial sources onto a common grid, -derives variables with registered operators, and catalogs every result with a -`spec_hash` so a derived layer can be traced back to exactly how it was made. +### Changed +- **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 - **Declarative derivations.** A `Derivation` names a source, a target grid and diff --git a/README.md b/README.md index beb0738..9cb8f91 100644 --- a/README.md +++ b/README.md @@ -117,28 +117,31 @@ in one TOML file, the counterpart of a TerraME `fill` script, and run from the command line: ```bash -disscube validate examples/pipelines/itaituba_fill.toml # no downloads -disscube run examples/pipelines/itaituba_fill.toml --workspace outputs/itaituba +disscube validate examples/pipelines/quickstart.toml # no downloads +disscube run examples/pipelines/quickstart.toml --workspace outputs/quickstart ``` ```toml schema = 1 +name = "Quickstart" [grid] -name = "ilha_do_maranhao" -bbox = [-44.35, -2.62, -44.20, -2.47] # WGS84 → BDC Albers grid +name = "demo/300m" +crs = "EPSG:31983" +bbox = [570000.0, 9708000.0, 582000.0, 9720000.0] resolution = 300 [[source]] -id = "lulc_{year}" -type = "mapbiomas" -years = [2000, 2020] +id = "landuse" +type = "file" +path = "../data/quickstart/landuse.tif" +crs = "EPSG:31983" [[derive]] -target = "urban_pct" -source = "lulc_{year}" +target = "forest_pct" +source = "landuse" operator = "percentage" -class_code = 24 +class_code = 3 ``` See [`docs/guides/pipeline_files.md`](docs/guides/pipeline_files.md) and @@ -146,34 +149,19 @@ See [`docs/guides/pipeline_files.md`](docs/guides/pipeline_files.md) and ## Examples -[`examples/`](examples/) has runnable scripts; 01–06 need no downloads and run -in seconds: - -- **01–03, synthetic data** — a quickstart with raster operators, vector - drivers, and time series handed off to DisSModel. -- **04–06, real data** — the three *Fill* examples shipped with TerraME - (Itaituba, Emas National Park, Brazilian Amazon), derived with DisSCube and - compared cell by cell with TerraME's own output. Data are bundled in - [`examples/data/terrame/`](examples/data/terrame/). -- **07, Brazil Data Cube** — reads the Landsat 16-day data cube over Ilha do - Maranhão through the BDC STAC catalog (only the pixels of the area are - fetched) and derives NDVI, MNDWI and open-water drivers on a 300 m grid - snapped to the BDC Albers mesh. Needs network and `pip install -e ".[bdc]"`; - with `--offline` it runs on a synthetic stand-in. -- **08, MapBiomas** — land use of Ilha do Maranhão in 2000 and 2020 (Collection - 11, 30 m), read straight from the national files and aggregated on the same - grid as 07 (`majority`, urban/forest/mangrove `percentage`), then handed to - DisSModel with `to_lucc_data()`. Needs network; `--offline` runs a stand-in. -- **09, PRODES** — deforestation in the LuccME Lab15 area (south of Santarém, - Pará) at the end of 2008, 2016 and 2024: the PRODES edition is downloaded - once and cached, its legend read from the `.qml`, and the maps aggregated on - a 500 m BDC grid. Needs network; `--offline` runs a stand-in. +[`examples/`](examples/) has runnable, self-contained examples that run offline in seconds: + +- **01–03, synthetic data** — a quickstart with raster operators, vector drivers, and time series handed off to DisSModel. +- **quickstart.toml** — declarative pipeline equivalent to example 01. ```bash python examples/01_quickstart.py -python examples/04_terrame_fill_itaituba.py +disscube run examples/pipelines/quickstart.toml ``` +Real data workflows, TerraME parity benchmarks, and large-scale case studies are maintained in +[**LambdaGeo/disscube-recipes**](https://github.com/LambdaGeo/disscube-recipes) (`cases/terrame_fill`, `cases/ilha_maranhao`, `cases/prodes_br163`, `cases/luccme_br`). + See [`examples/README.md`](examples/README.md) for the full list, and [`docs/terrame_fill_correspondence.md`](docs/terrame_fill_correspondence.md) for how DisSCube's operators relate to TerraME's *Fill*. diff --git a/docs/examples.md b/docs/examples.md index eb4df86..0e21d4b 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -1,23 +1,21 @@ # Examples The [`examples/`](https://github.com/DisSModel/disscube/tree/main/examples) -folder has runnable scripts that need nothing beyond `pip install -e .`: each -one works in a temporary directory and finishes in a few seconds. Pass a -directory to keep the catalog, the inputs and the derived Zarr stores for -inspection (e.g. in QGIS): +folder contains runnable, self-contained examples demonstrating the two primary modes of using DisSCube: +the **Python API** (`CubeClient`, `Derivation`) and the **declarative TOML pipeline** (`disscube validate` / `disscube run`). + +All examples run offline, generate or consume tiny synthetic data (< 50 KB), and complete in seconds without downloading external assets: ```bash python examples/01_quickstart.py # temporary workspace -python examples/01_quickstart.py ./scratch # keep the outputs +python examples/01_quickstart.py ./scratch # keep the outputs for inspection ``` -All examples are executed by the test suite (`tests/test_examples.py`), so -they stay in sync with the API. +All examples are executed by the test suite (`tests/test_examples.py`), ensuring they always stay in sync with the codebase. -## Learning the API — synthetic data +## Python API Examples (Synthetic Data) -Examples 01–03 generate their own inputs, so every number they print can be -checked by hand. +Examples 01–03 generate their own inputs locally, allowing every calculated number to be verified directly: | Example | What it shows | |---|---| @@ -25,66 +23,36 @@ checked by hand. | [`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`) | -## Real data — TerraME's *Fill* examples - -Examples 04–06 use the three *Fill* examples shipped with TerraME's `gis` -package, bundled in -[`examples/data/terrame/`](https://github.com/DisSModel/disscube/tree/main/examples/data/terrame) -together with the cellular spaces TerraME produced. Each derives the same -attributes with DisSCube and prints a cell-by-cell comparison with TerraME's -own output. - -| Example | Cells | What it shows | Result | -|---|---|---|---| -| [`04_terrame_fill_itaituba.py`](https://github.com/DisSModel/disscube/blob/main/examples/04_terrame_fill_itaituba.py) | 620 × 5 km | TerraME's [Fill tutorial](https://github.com/TerraME/terrame/wiki/Fill): `mean`, `percentage` × `coverage_purity`, `min_distance` | averages and class coverage reproduce TerraME | -| [`05_terrame_fill_emas.py`](https://github.com/DisSModel/disscube/blob/main/examples/05_terrame_fill_emas.py) | 5 514 × 500 m | `presence` of lines, `max` / `min` of a raster; a study area defined by a limit polygon | 98.4–99.7 % of cells identical | -| [`06_terrame_fill_amazonia.py`](https://github.com/DisSModel/disscube/blob/main/examples/06_terrame_fill_amazonia.py) | 2 229 × 50 km | PRODES coverage with a declared nodata, distances to roads and ports | coverage identical wherever PRODES has data | - -What the differences mean — and which TerraME operations DisSCube does not -support yet — is discussed in -[TerraME Fill Cells Correspondence](terrame_fill_correspondence.md). - -## Real data — Brazil Data Cube - -Example 07 reads a satellite data cube straight from the Brazil Data Cube -STAC catalog; see [BDC Integration](guides/bdc.md#reading-data-cubes-via-stac). - -| Example | Cells | What it shows | -|---|---|---| -| [`07_bdc_cube.py`](https://github.com/DisSModel/disscube/blob/main/examples/07_bdc_cube.py) | 58 × 59 × 300 m | `LANDSAT-16D-1` over Ilha do Maranhão: windowed reads, dry-season median, `mean` of NDVI and MNDWI, `percentage` of open water on a grid snapped to the BDC Albers mesh | - -It needs network access and `pip install -e ".[bdc]"`. Without network, or -with `--offline`, it runs on a synthetic scene with the same grid and value -ranges and prints a warning that the numbers are not BDC data. +## Declarative Pipeline Files (TOML) -## Real data — MapBiomas +Pipelines can be declared in clean, version-controlled TOML files in +[`examples/pipelines/`](https://github.com/DisSModel/disscube/tree/main/examples/pipelines): -Example 08 reads the MapBiomas annual land-cover maps straight from the -national files; see [MapBiomas](guides/mapbiomas.md). - -| Example | Cells | What it shows | +| Pipeline | What it shows | Execution | |---|---|---| -| [`08_mapbiomas_land_use.py`](https://github.com/DisSModel/disscube/blob/main/examples/08_mapbiomas_land_use.py) | 58 × 59 × 300 m, 2 years | Collection 11 (30 m) over Ilha do Maranhão in 2000 and 2020: `majority` and `percentage` (urban, forest, mangrove) on the grid of example 07, `(time, y, x)` series handed to DisSModel with `to_lucc_data()` | - -It needs network access. Without it, or with `--offline`, it runs on a -synthetic scene with the same grid and class codes and says so. +| [`quickstart.toml`](https://github.com/DisSModel/disscube/blob/main/examples/pipelines/quickstart.toml) | Declarative counterpart of `01_quickstart.py`: derives `percentage`, `majority`, and `mean` on a 300 m grid | `disscube run examples/pipelines/quickstart.toml` | -## Real data — PRODES +Validate and execute: +```bash +disscube validate examples/pipelines/quickstart.toml +disscube run examples/pipelines/quickstart.toml --workspace outputs/quickstart +``` -Example 09 downloads the current PRODES edition once (≈130 MB, cached) and -reads its legend from the `.qml` in the ZIP; see -[Land-cover maps: PRODES and classified maps](guides/landcover.md). +See [Pipeline files (TOML)](guides/pipeline_files.md) for full syntax and options. -| Example | Cells | What it shows | -|---|---|---| -| [`09_prodes_deforestation.py`](https://github.com/DisSModel/disscube/blob/main/examples/09_prodes_deforestation.py) | 94 × 86 × 500 m, 3 years | The LuccME Lab15 area south of Santarém at the end of 2008, 2016 and 2024: forest / deforested / other, `percentage` of deforested and forest per cell | +## Real Data & Case Studies (External Repository) -## Pipeline files (TOML) +Real-world datasets, historical reproductions and large-scale case studies are maintained in the dedicated +[**DisSCube Recipes and Case Studies**](https://github.com/LambdaGeo/disscube-recipes) repository: -The same preparations can be declared as TOML files in -[`examples/pipelines/`](https://github.com/DisSModel/disscube/tree/main/examples/pipelines) -and run with `disscube run ` — including TerraME's Itaituba Fill -tutorial, offline. See [Pipeline files](guides/pipeline_files.md). +- **TerraME Fill Parity (`cases/terrame_fill`)**: Cell-by-cell numerical parity against TerraME's C++ `fillCellularSpace` across three reference areas: + - `itaituba`: 620 cells (5 km) with elevation, multi-class deforestation, and proximity drivers. + - `emas`: 5 514 cells (500 m) with firebreaks, rivers, and accumulation rasters. + - `amazonia`: 2 229 cells (50 km) with PRODES deforestation and indigenous territory area fractions. + - See [TerraME Fill Cells Correspondence](terrame_fill_correspondence.md). +- **Ilha do Maranhão (`cases/ilha_maranhao`)**: Brazil Data Cube (Landsat-16D STAC) and MapBiomas land-cover series (2000, 2020) aligned on a 300 m BDC Albers grid with Python hand-off to DisSModel. +- **PRODES BR-163 (`cases/prodes_br163`)**: Multi-year deforestation monitoring in the Mojuí dos Campos / BR-163 corridor (500 m grid). +- **LUCCME-BR (`cases/luccme_br`)**: National-scale land-use change reconstruction for Brazil (BigEarth). ## Scope diff --git a/docs/guides/bdc.md b/docs/guides/bdc.md index 0779c9e..af07ecb 100644 --- a/docs/guides/bdc.md +++ b/docs/guides/bdc.md @@ -89,7 +89,7 @@ or scale) produces a new file, a new checksum and a new product, while re-fetching the same data is a cache hit. For a layer computed from several assets — an index such as MNDWI — build the `Window2D` yourself and call `disscube.sources.register_raster(cube, source_id, window, out_dir, provenance)`; see -`examples/07_bdc_cube.py`. +the `ilha_maranhao` case study in the recipes repository (`cases/ilha_maranhao`). What the reader does for you: diff --git a/docs/guides/pipeline_files.md b/docs/guides/pipeline_files.md index e15e482..c475eea 100644 --- a/docs/guides/pipeline_files.md +++ b/docs/guides/pipeline_files.md @@ -8,14 +8,14 @@ before anything is downloaded, and recorded in the provenance of every product. ```bash -disscube validate examples/pipelines/itaituba_fill.toml # check the file, no downloads -disscube run examples/pipelines/itaituba_fill.toml --workspace outputs/itaituba +disscube validate examples/pipelines/quickstart.toml # check the file, no downloads +disscube run examples/pipelines/quickstart.toml --workspace outputs/quickstart ``` ## From `fill` to `[[derive]]` -TerraME's Fill tutorial (`examples/data/terrame/itaituba/itaituba.lua`) and -its pipeline file (`examples/pipelines/itaituba_fill.toml`), side by side: +TerraME's Fill tutorial (see `cases/terrame_fill` in the recipes repository) and +its DisSCube pipeline representation, side by side: ```lua -- TerraME @@ -177,10 +177,6 @@ file that declared it. | File | What it does | Needs | |---|---|---| -| `examples/pipelines/itaituba_fill.toml` | TerraME's Fill tutorial | nothing (bundled data) | -| `examples/pipelines/ilha_bdc.toml` | example 07: BDC NDVI and MNDWI | network, `.[bdc]` | -| `examples/pipelines/ilha_mapbiomas.toml` | example 08: MapBiomas 2000 and 2020 | network | -| `examples/pipelines/lab15_prodes.toml` | example 09: PRODES 2008, 2016, 2024 | network | +| `examples/pipelines/quickstart.toml` | Quickstart: synthetic land-use & elevation (`percentage`, `majority`, `mean`) | nothing (offline) | -All of them are validated by the test suite; the Itaituba one is also run and -compared with the same derivations made through the Python API. +For full real-world case studies and TerraME parity benchmarks, see the [DisSCube Recipes](https://github.com/LambdaGeo/disscube-recipes) repository (`cases/terrame_fill`, `cases/ilha_maranhao`, `cases/prodes_br163`, `cases/luccme_br`). diff --git a/docs/terrame_fill_correspondence.md b/docs/terrame_fill_correspondence.md index 80e082c..046f8db 100644 --- a/docs/terrame_fill_correspondence.md +++ b/docs/terrame_fill_correspondence.md @@ -35,8 +35,9 @@ how faithful they are — but the engineering around them: `PROJCS["unknown"]` round-trip problems. A TerraME fill script can also be written as a DisSCube pipeline file — -see [Pipeline files (TOML)](guides/pipeline_files.md) and -`examples/pipelines/itaituba_fill.toml`, the tutorial below as TOML. +see [Pipeline files (TOML)](guides/pipeline_files.md), +the [quickstart pipeline](https://github.com/DisSModel/disscube/blob/main/examples/pipelines/quickstart.toml), +and the parity cases in the recipes repository (`cases/terrame_fill`). ## Strategy → operator correspondence @@ -82,12 +83,13 @@ Amazônia — each with its input layers, the Lua script and the cellular space TerraME produced. Those outputs are the reference: the same inputs are derived with DisSCube on the same grid and compared cell by cell. -The data are bundled in `examples/data/terrame/` (provenance, license and -checksums in its README). The comparison runs in CI in -`tests/test_terrame_parity.py` — passing tests pin the parity below, strict -`xfail` tests record the known gaps — and examples -`04_terrame_fill_itaituba.py`, `05_terrame_fill_emas.py` and -`06_terrame_fill_amazonia.py` print the comparison for each dataset. +The full parity benchmark suite is maintained and executed on CI in the +[DisSCube Case Studies and Recipes](https://github.com/LambdaGeo/disscube-recipes) +repository (`cases/terrame_fill`), where input datasets are fetched from TerraME +with cryptographic checksum verification (`./run.sh all `). + +Below is the summary of the correspondences and cell-by-cell numerical parity +verified across the datasets. ### Itaituba — 620 cells, 5 km diff --git a/examples/04_terrame_fill_itaituba.py b/examples/04_terrame_fill_itaituba.py deleted file mode 100644 index 89a5dc8..0000000 --- a/examples/04_terrame_fill_itaituba.py +++ /dev/null @@ -1,120 +0,0 @@ -""" -04 — TerraME's *Fill* tutorial (Itaituba) reproduced with DisSCube. - -TerraME's ``gis`` package fills a 5 km cellular space over Itaituba (Pará) -from five layers; the script and its output ship with TerraME and are copied -in ``examples/data/terrame/itaituba/`` (see the README there). This example -derives the same attributes with DisSCube and compares them cell by cell with -TerraME's own output: - - TerraME fill DisSCube operator - ───────────────────────────────── ───────────────────────────────── - elevation average (raster) mean - defor_ coverage (raster) percentage (class_code=k) - distroad distance (lines) min_distance - distlocal distance (points) min_distance - population sum, area=true (polys) — not supported yet - -Coverage is where the two tools differ by design: TerraME divides by the -whole cell, DisSCube by the valid pixels, reporting the covered share as -``coverage_purity``. Their product reproduces TerraME's value. - - python examples/04_terrame_fill_itaituba.py # temporary workspace - python examples/04_terrame_fill_itaituba.py ./scratch # keep the outputs -""" - -import sys -import tempfile -from pathlib import Path - -import geopandas as gpd -import numpy as np - -from disscube import CubeClient, GridSpec, SpatialDerivation, SpatialSource, Variable - -DATA = Path(__file__).resolve().parent / "data" / "terrame" / "itaituba" -CRS = "EPSG:29191" # SAD69 / UTM zone 21S -RESOLUTION = 5_000.0 -# Class names as given in itaituba.lua (87 and 167); the others are unnamed there. -CLASSES = {7: "class 7", 87: "deforestation", 167: "river", 255: "class 255"} - - -def zipped(name: str) -> str: - """Shapefiles are stored zipped; GDAL reads them in place via zip://.""" - return f"zip://{DATA / name}.zip" - - -def main(workspace: Path) -> None: - # TerraME's output: 620 cells with the filled attributes. Its .prj is - # wrong (NAD83_Austin), so the real CRS is assigned explicitly. - terrame = gpd.read_file(zipped("itaituba")).set_crs(CRS, allow_override=True) - - workspace.mkdir(parents=True, exist_ok=True) - cube = CubeClient(catalog=str(workspace / "catalog.db"), store=str(workspace / "store")) - - # Same grid as TerraME: 5 km cells over the extent of its cellular space. - grid_id = "itaituba/5km" - cube.register_grid(GridSpec( - id=grid_id, type="local", crs=CRS, resolution=RESOLUTION, - bbox=terrame.total_bounds.tolist(), description="TerraME Fill tutorial grid", - )) - - sources = { - "elevation": ("raster", str(DATA / "itaituba-elevation.tif")), - "deforestation": ("raster", str(DATA / "itaituba-deforestation.tif")), - "roads": ("vector", zipped("itaituba-roads")), - "localities": ("vector", zipped("itaituba-localities")), - } - for sid, (fmt, url) in sources.items(): - cube.register_spatial_source(SpatialSource(id=sid, name=f"Itaituba {sid}", format=fmt, asset_url=url, crs=CRS)) - - def fill(source_id: str, *variables: Variable) -> None: - cube.derive(SpatialDerivation(source_id=source_id, grid_id=grid_id, role="driver", variables=list(variables))) - - fill("elevation", Variable(name="elevation", operator="mean")) - fill("deforestation", *(Variable(name=f"defor_{k}", operator="percentage", class_code=k) for k in CLASSES)) - fill("roads", Variable(name="distroad", operator="min_distance")) - fill("localities", Variable(name="distlocal", operator="min_distance")) - - # Match cells: TerraME numbers rows from the south, DisSCube grids are north-up. - xmin, _, _, ymax = terrame.total_bounds - centroids = terrame.geometry.centroid - cols = np.floor((centroids.x.to_numpy() - xmin) / RESOLUTION).astype(int) - rows = np.floor((ymax - centroids.y.to_numpy()) / RESOLUTION).astype(int) - - def ours(name: str, with_purity: bool = False) -> np.ndarray: - da = cube.load(name, grid_id=grid_id) - if with_purity: - da = da * da.coords["coverage_purity"] - return da.to_numpy()[rows, cols] - - print(f"{len(terrame)} cells, {RESOLUTION / 1000:.0f} km — DisSCube vs TerraME, cell by cell\n") - print(f"{'attribute':<12} {'DisSCube':<32} {'mean |diff|':>12} {'max |diff|':>11} within") - print("─" * 88) - - def row(attr, label, values, reference, tol, unit): - diff = np.abs(values - reference) - print(f"{attr:<12} {label:<32} {diff.mean():>10.2f}{unit:>2} {diff.max():>9.2f}{unit:>2}" - f" {np.mean(diff <= tol):>4.0%} of cells within {tol:g}{unit}") - - row("elevation", "mean", ours("elevation"), terrame["elevation"], 1, "m") - for k, label in CLASSES.items(): - ref = terrame[f"defor_{k}"].to_numpy() - row(f"defor_{k}", f"percentage ({label})", ours(f"defor_{k}") * 100, ref, 1, "pp") - row("", " × coverage_purity", ours(f"defor_{k}", with_purity=True) * 100, ref, 1, "pp") - row("distroad", "min_distance", ours("distroad"), terrame["distroad"], 100, "m") - row("distlocal", "min_distance", ours("distlocal"), terrame["distlocal"], 100, "m") - print(f"{'population':<12} {'— (area-weighted sum: not yet)':<32}") - - print("\nReading: averages and class coverage (with purity) reproduce TerraME;") - print("min_distance is a raster approximation between cell centres, while TerraME") - print("measures the exact distance from each cell polygon — see") - print("docs/terrame_fill_correspondence.md.") - - -if __name__ == "__main__": - if len(sys.argv) > 1: - main(Path(sys.argv[1])) - else: - with tempfile.TemporaryDirectory() as tmp: - main(Path(tmp)) diff --git a/examples/05_terrame_fill_emas.py b/examples/05_terrame_fill_emas.py deleted file mode 100644 index 0f9dea6..0000000 --- a/examples/05_terrame_fill_emas.py +++ /dev/null @@ -1,100 +0,0 @@ -""" -05 — TerraME's *Fill* example for Emas National Park. - -TerraME's ``gis`` package fills a 500 m cellular space over Emas National Park -(Goiás) from firebreak and river lines and a 30 m fire-accumulation raster; -the script and its output are copied in ``examples/data/terrame/emas/``. -This example derives the same attributes with DisSCube and compares them -cell by cell with TerraME's output: - - TerraME fill DisSCube operator - ───────────────────────────── ───────────────── - firebreak presence (lines) presence - river presence (lines) presence - maxcover maximum (raster) max - mincover minimum (raster) min - -TerraME kept only the cells that touch the park limit (``input = "limit"``), -so DisSCube derives the full rectangle and the comparison uses TerraME's cells. - - python examples/05_terrame_fill_emas.py # temporary workspace - python examples/05_terrame_fill_emas.py ./scratch # keep the outputs -""" - -import sys -import tempfile -from pathlib import Path - -import geopandas as gpd -import numpy as np - -from disscube import CubeClient, GridSpec, SpatialDerivation, SpatialSource, Variable - -DATA = Path(__file__).resolve().parent / "data" / "terrame" / "emas" -CRS = "EPSG:29192" # SAD69 / UTM zone 22S -RESOLUTION = 500.0 - - -def zipped(name: str) -> str: - """Shapefiles are stored zipped; GDAL reads them in place via zip://.""" - return f"zip://{DATA / name}.zip" - - -def main(workspace: Path) -> None: - workspace.mkdir(parents=True, exist_ok=True) - - # TerraME's output (5 514 cells). Its .prj is a generic UTM 22S, so the - # CRS of the input layers is assigned explicitly. - terrame = gpd.read_file(zipped("emas")).set_crs(CRS, allow_override=True) - - cube = CubeClient(catalog=str(workspace / "catalog.db"), store=str(workspace / "store")) - grid_id = "emas/500m" - cube.register_grid(GridSpec( - id=grid_id, type="local", crs=CRS, resolution=RESOLUTION, - bbox=terrame.total_bounds.tolist(), description="TerraME Emas example grid", - )) - - for sid, fmt, url in [ - ("firebreak", "vector", zipped("emas-firebreak")), - ("river", "vector", zipped("emas-river")), - ("cover", "raster", str(DATA / "emas-accumulation.tif")), - ]: - cube.register_spatial_source(SpatialSource(id=sid, name=f"Emas {sid}", format=fmt, asset_url=url, crs=CRS)) - - def fill(source_id: str, *variables: Variable) -> None: - cube.derive(SpatialDerivation(source_id=source_id, grid_id=grid_id, role="driver", variables=list(variables))) - - fill("firebreak", Variable(name="firebreak", operator="presence")) - fill("river", Variable(name="river", operator="presence")) - # One raster, two variables: a derivation may compute several at once. - fill("cover", Variable(name="maxcover", operator="max"), Variable(name="mincover", operator="min")) - - # Match cells by centroid (DisSCube grids are north-up). - xmin, _, _, ymax = terrame.total_bounds - centroids = terrame.geometry.centroid - cols = np.floor((centroids.x.to_numpy() - xmin) / RESOLUTION).astype(int) - rows = np.floor((ymax - centroids.y.to_numpy()) / RESOLUTION).astype(int) - - grid = cube.catalog.get_grid(grid_id) - print(f"{len(terrame)} TerraME cells (of a {grid.rows} × {grid.cols} rectangle), " - f"{RESOLUTION:.0f} m — DisSCube vs TerraME, cell by cell\n") - print(f"{'attribute':<11} {'DisSCube':<10} {'identical cells':>16} {'differing':>10}") - print("─" * 52) - for attr, op in [("firebreak", "presence"), ("river", "presence"), ("maxcover", "max"), ("mincover", "min")]: - ours = cube.load(attr, grid_id=grid_id).to_numpy()[rows, cols] - same = ours == terrame[attr].to_numpy() - print(f"{attr:<11} {op:<10} {same.mean():>15.1%} {int((~same).sum()):>10}") - - print("\nReading: the differing cells sit on cell borders. TerraME marks every cell") - print("a line touches and assigns each pixel to the cell containing its centre;") - print("DisSCube rasterizes lines through cell centres and, for min/max, also counts") - print("pixels that straddle a border (500 m / 30 m is not an integer). See") - print("docs/terrame_fill_correspondence.md.") - - -if __name__ == "__main__": - if len(sys.argv) > 1: - main(Path(sys.argv[1])) - else: - with tempfile.TemporaryDirectory() as tmp: - main(Path(tmp)) diff --git a/examples/06_terrame_fill_amazonia.py b/examples/06_terrame_fill_amazonia.py deleted file mode 100644 index de044fc..0000000 --- a/examples/06_terrame_fill_amazonia.py +++ /dev/null @@ -1,113 +0,0 @@ -""" -06 — TerraME's *Fill* example for the Brazilian Amazon. - -TerraME's ``gis`` package fills a 50 km cellular space over the Brazilian -Amazon from the PRODES deforestation map, roads, ports and indigenous lands; -the script and its output are copied in ``examples/data/terrame/amazonia/``. -This example derives the same attributes with DisSCube and compares them -cell by cell with TerraME's output: - - TerraME fill DisSCube operator - ────────────────────────────────── ───────────────────────────── - prodes_ coverage (raster) percentage × coverage_purity - distroads distance (lines) min_distance - distports distance (points) min_distance - protected area (polygons) — not supported yet - -PRODES is a 5 km raster with a declared nodata (255); cells with no PRODES -data at all are NaN in DisSCube (coverage_purity 0) and 0 in TerraME. - - python examples/06_terrame_fill_amazonia.py # temporary workspace - python examples/06_terrame_fill_amazonia.py ./scratch # keep the outputs -""" - -import sys -import tempfile -from pathlib import Path - -import geopandas as gpd -import numpy as np - -from disscube import CubeClient, GridSpec, SpatialDerivation, SpatialSource, Variable - -DATA = Path(__file__).resolve().parent / "data" / "terrame" / "amazonia" -CRS = "EPSG:29191" # SAD69 / UTM zone 21S -RESOLUTION = 50_000.0 -PRODES_CLASSES = (10, 208) - - -def zipped(name: str) -> str: - """Shapefiles are stored zipped; GDAL reads them in place via zip://.""" - return f"zip://{DATA / name}.zip" - - -def main(workspace: Path) -> None: - workspace.mkdir(parents=True, exist_ok=True) - - # TerraME's output (2 229 cells). Its .prj is wrong (NAD83_Austin), so the - # CRS of the input layers is assigned explicitly. - terrame = gpd.read_file(zipped("amazonia")).set_crs(CRS, allow_override=True) - - cube = CubeClient(catalog=str(workspace / "catalog.db"), store=str(workspace / "store")) - grid_id = "amazonia/50km" - cube.register_grid(GridSpec( - id=grid_id, type="local", crs=CRS, resolution=RESOLUTION, - bbox=terrame.total_bounds.tolist(), description="TerraME Amazonia example grid", - )) - - for sid, fmt, url in [ - ("prodes", "raster", str(DATA / "amazonia-prodes.tif")), - ("roads", "vector", zipped("amazonia-roads")), - ("ports", "vector", zipped("amazonia-ports")), - ]: - cube.register_spatial_source(SpatialSource(id=sid, name=f"Amazonia {sid}", format=fmt, asset_url=url, crs=CRS)) - - def fill(source_id: str, *variables: Variable) -> None: - cube.derive(SpatialDerivation(source_id=source_id, grid_id=grid_id, role="driver", variables=list(variables))) - - fill("prodes", *(Variable(name=f"prodes_{k}", operator="percentage", class_code=k) for k in PRODES_CLASSES)) - fill("roads", Variable(name="distroads", operator="min_distance")) - fill("ports", Variable(name="distports", operator="min_distance")) - - # Match cells by centroid (DisSCube grids are north-up). - xmin, _, _, ymax = terrame.total_bounds - centroids = terrame.geometry.centroid - cols = np.floor((centroids.x.to_numpy() - xmin) / RESOLUTION).astype(int) - rows = np.floor((ymax - centroids.y.to_numpy()) / RESOLUTION).astype(int) - - grid = cube.catalog.get_grid(grid_id) - print(f"{len(terrame)} TerraME cells (of a {grid.rows} × {grid.cols} rectangle), " - f"{RESOLUTION / 1000:.0f} km — DisSCube vs TerraME, cell by cell\n") - - print("coverage (percent of the cell):") - for k in PRODES_CLASSES: - da = cube.load(f"prodes_{k}", grid_id=grid_id) - ours = (da * da.coords["coverage_purity"]).to_numpy()[rows, cols] * 100 - theirs = terrame[f"prodes_{k}"].to_numpy() - no_data = np.isnan(ours) - diff = np.abs(ours[~no_data] - theirs[~no_data]) - print(f" prodes_{k:<4} percentage × coverage_purity: max |diff| {diff.max():.2f} pp in " - f"{(~no_data).sum()} cells with data; {no_data.sum()} cells without PRODES data " - f"(TerraME: 0, DisSCube: NaN)") - - print("\ndistance (m):") - for attr in ("distroads", "distports"): - ours = cube.load(attr, grid_id=grid_id).to_numpy()[rows, cols] - diff = ours - terrame[attr].to_numpy() - print(f" {attr:<10} min_distance: mean |diff| {np.abs(diff).mean():>8,.0f} " - f"mean bias {diff.mean():>+8,.0f} max |diff| {np.abs(diff).max():>8,.0f}") - - print("\nprotected area (fraction of the cell covered by indigenous lands): not supported yet") - - print("\nReading: class coverage reproduces TerraME wherever PRODES has data;") - print("min_distance is a raster approximation between 50 km cell centres, so it") - print("differs by tens of kilometres from TerraME's per-polygon distance. See") - print("docs/terrame_fill_correspondence.md.") - - -if __name__ == "__main__": - if len(sys.argv) > 1: - main(Path(sys.argv[1])) - else: - with tempfile.TemporaryDirectory() as tmp: - main(Path(tmp)) diff --git a/examples/07_bdc_cube.py b/examples/07_bdc_cube.py deleted file mode 100644 index 2000c5c..0000000 --- a/examples/07_bdc_cube.py +++ /dev/null @@ -1,195 +0,0 @@ -""" -07 — Brazil Data Cube: model drivers from a satellite data cube. - -Reads the Landsat 16-day data cube of the Brazil Data Cube (``LANDSAT-16D-1``) -over Ilha do Maranhão through the BDC STAC catalog, fetching only the pixels -of the area of interest, and turns one dry season into three drivers on a -300 m grid snapped to the BDC Albers mesh: - - - ``ndvi_mean`` mean NDVI (vegetation vigour) - - ``mndwi_mean`` mean MNDWI, (green − SWIR1) / (green + SWIR1) - - ``water_pct`` fraction of each cell with MNDWI > 0 (open water) - -The season is reduced to a per-pixel median of the 16-day composites, which -removes most of the remaining clouds. Each source is registered with the -SHA-256 of its file as checksum and a ``.provenance.json`` sidecar -(collection, items, period, reducer, scale), so re-running with another -period recomputes the drivers instead of returning cached ones. - - pip install -e ".[bdc]" # pystac-client for the search - python examples/07_bdc_cube.py # temporary workspace - python examples/07_bdc_cube.py ./scratch # keep the outputs - python examples/07_bdc_cube.py --offline # no network: synthetic stand-in - -Without network access (or with ``--offline`` / ``DISSCUBE_OFFLINE=1``, as in -the test suite) the script builds a synthetic scene with the same grid, CRS -and value ranges, and says so loudly: those numbers are not BDC data. -""" - -import json -import os -import sys -import tempfile -from pathlib import Path - -import numpy as np -import rasterio -from rasterio.crs import CRS -from rasterio.transform import from_origin - -from disscube import CubeClient, Derivation, SpatialSource -from disscube.sources import Window2D, register_raster -from disscube.utils import BDC_CRS, register_local_grid, sha256_file - -COLLECTION = "LANDSAT-16D-1" -PERIOD = "2020-07-01/2020-09-30" # dry season -BBOX = (-44.35, -2.62, -44.20, -2.47) # São Luís, Ilha do Maranhão (WGS84) -GRID_RES = 300.0 - - -# --------------------------------------------------------------------------- -# Inputs: real (BDC) or synthetic (offline) — both registered with provenance -# --------------------------------------------------------------------------- - -def sources_from_bdc(cube: CubeClient, raw: Path) -> Window2D: - """Register ``ndvi`` and ``mndwi`` from the BDC; return the MNDWI layer.""" - from disscube.sources import normalized_difference - from disscube.sources.bdc import ( - BDC_INDEX_SCALE, - BDC_STAC_URL, - items_provenance, - read_composite, - register_bdc_source, - search_items, - ) - - items = search_items(COLLECTION, BBOX, PERIOD) - if not items: - raise RuntimeError(f"no {COLLECTION} items over {BBOX} in {PERIOD}") - print(f"BDC items : {len(items)} ({items[0].id} … {items[-1].id})") - - # The catalog declares no scale; BDC stores indices as int16 × 10 000. - register_bdc_source(cube, "ndvi", COLLECTION, "NDVI", BBOX, PERIOD, raw, - items=items, scale=BDC_INDEX_SCALE) - - # MNDWI is a ratio, so a common scale factor on green and SWIR cancels out. - def season(asset): - return read_composite(COLLECTION, asset, BBOX, PERIOD, items=items) - - mndwi = normalized_difference(season("green"), season("swir16")) - register_raster( - cube, "mndwi", mndwi, raw, - provenance={ - "stac_url": BDC_STAC_URL, "collection": COLLECTION, "assets": ["green", "swir16"], - "expression": "(green - swir16) / (green + swir16)", "bbox_geo": list(BBOX), - "period": PERIOD, "reducer": "median", "items": items_provenance(items), - }, - name=f"{COLLECTION} MNDWI (median of {PERIOD})", time=2020, - tags=["bdc", f"collection:{COLLECTION}", f"period:{PERIOD}"], - ) - return mndwi - - -def sources_synthetic(cube: CubeClient, raw: Path) -> Window2D: - """Water to the north-west, an urban core, vegetation elsewhere — in BDC Albers.""" - from pyproj import Transformer - - to_albers = Transformer.from_crs("EPSG:4326", BDC_CRS, always_xy=True) - x0, y1 = to_albers.transform(BBOX[0], BBOX[3]) - x1, y0 = to_albers.transform(BBOX[2], BBOX[1]) - res = 30.0 - cols, rows = int((x1 - x0) / res), int((y1 - y0) / res) - transform = from_origin(x0, y1, res, res) - - rng = np.random.default_rng(7) - x = np.linspace(0, 1, cols)[None, :].repeat(rows, axis=0) - y = np.linspace(0, 1, rows)[:, None].repeat(cols, axis=1) - water = (x + (1 - y)) < 0.55 - urban = (x - 0.6) ** 2 + (y - 0.5) ** 2 < 0.04 - ndvi = np.clip(0.75 - 0.45 * urban + 0.05 * rng.standard_normal((rows, cols)), -1, 1) - ndvi[water] = -0.1 - mndwi = np.clip(-0.35 + 0.05 * rng.standard_normal((rows, cols)), -1, 1) - mndwi[water] = 0.45 - - crs = CRS.from_string(BDC_CRS) - layers = {name: Window2D(data=arr.astype("float32"), transform=transform, crs=crs) - for name, arr in (("ndvi", ndvi), ("mndwi", mndwi))} - for name, window in layers.items(): - register_raster(cube, name, window, raw, - provenance={"synthetic": True, "note": "offline stand-in, not BDC data"}, - name=f"{name} (synthetic stand-in)", time=2020, tags=["synthetic"]) - return layers["mndwi"] - - -def register_water_mask(cube: CubeClient, raw: Path, mndwi: Window2D) -> None: - """1 = open water (MNDWI > 0), 0 = land, 255 = no observation.""" - mask = np.where(mndwi.data > 0, 1, 0).astype("uint8") - mask[np.isnan(mndwi.data)] = 255 - rows, cols = mask.shape - path = raw / "water.tif" - with rasterio.open(path, "w", driver="GTiff", height=rows, width=cols, count=1, - dtype="uint8", crs=BDC_CRS, transform=mndwi.transform, nodata=255) as dst: - dst.write(mask, 1) - cube.register_spatial_source(SpatialSource( - id="water", name="open water (MNDWI > 0)", format="raster", asset_url=str(path), - crs=BDC_CRS, time=2020, checksum=sha256_file(path), tags=["derived:mndwi>0"], - )) - - -# --------------------------------------------------------------------------- -# Pipeline -# --------------------------------------------------------------------------- - -def main(workspace: Path, offline: bool) -> None: - raw = workspace / "raw" - raw.mkdir(parents=True, exist_ok=True) - cube = CubeClient(catalog=str(workspace / "catalog.db"), store=str(workspace / "store")) - - # 300 m grid snapped to the national BDC Albers mesh - grid = register_local_grid(cube, name="ilha_do_maranhao", bbox_geo=BBOX, resolution=GRID_RES) - - label = "synthetic stand-in" - mndwi = None - if not offline: - try: - mndwi = sources_from_bdc(cube, raw) - label = f"BDC {COLLECTION}, median of {PERIOD}" - except Exception as exc: # noqa: BLE001 — any network/catalog failure falls back, loudly - print(f"!! BDC not reachable ({type(exc).__name__}: {exc})") - if mndwi is None: - print("!! OFFLINE — using a synthetic scene; the numbers below are NOT BDC data") - mndwi = sources_synthetic(cube, raw) - register_water_mask(cube, raw, mndwi) - - for d in ( - Derivation(target="ndvi_mean", source_id="ndvi", operator="mean"), - Derivation(target="mndwi_mean", source_id="mndwi", operator="mean"), - Derivation(target="water_pct", source_id="water", operator="percentage", class_code=1), - ): - cube.derive_declarative(d, grid_id=grid.id) - - ndvi = cube.load("ndvi_mean", grid_id=grid.id) - mndwi_mean = cube.load("mndwi_mean", grid_id=grid.id) - water = cube.load("water_pct", grid_id=grid.id) - provenance = json.loads((raw / "ndvi.provenance.json").read_text(encoding="utf-8")) - - print(f"source : {label}") - print(f"grid : {grid.id} {ndvi.shape[-2]} × {ndvi.shape[-1]} cells") - print(f"NDVI (cell means) : {float(np.nanmin(ndvi)):.2f} … {float(np.nanmax(ndvi)):.2f}, " - f"mean {float(np.nanmean(ndvi)):.2f}") - print(f"MNDWI (cell means) : {float(np.nanmin(mndwi_mean)):.2f} … " - f"{float(np.nanmax(mndwi_mean)):.2f}") - print(f"open water : {float(np.nanmean(water)):.1%} of the area") - print(f"cells without data : {int(np.isnan(ndvi.values).sum())}") - print(f"ndvi provenance : {raw / 'ndvi.provenance.json'} " - f"({len(provenance.get('items', []))} items, {provenance['checksum'][:19]}…)") - - -if __name__ == "__main__": - args = [a for a in sys.argv[1:] if not a.startswith("--")] - offline = "--offline" in sys.argv or os.environ.get("DISSCUBE_OFFLINE") == "1" - if args: - main(Path(args[0]), offline) - else: - with tempfile.TemporaryDirectory() as tmp: - main(Path(tmp), offline) diff --git a/examples/08_mapbiomas_land_use.py b/examples/08_mapbiomas_land_use.py deleted file mode 100644 index 712f8ed..0000000 --- a/examples/08_mapbiomas_land_use.py +++ /dev/null @@ -1,153 +0,0 @@ -""" -08 — MapBiomas: land use of Ilha do Maranhão in two years, on the model grid. - -Reads the MapBiomas annual land-cover maps (Collection 11, 30 m) of 2000 and -2020 over Ilha do Maranhão — only the pixels of the area, straight from the -national files — and aggregates them onto the same 300 m grid as example 07 -(snapped to the BDC Albers mesh): - - - ``landuse`` dominant class of each cell (``majority``) - - ``urban_pct`` fraction of urban area (class 24) - - ``forest_pct`` fraction of forest formation (class 3) - - ``mangrove_pct`` fraction of mangrove (class 5) - -Each year is a separate source (``time`` = year), so every variable loads as a -``(time, y, x)`` series and goes to DisSModel with ``to_lucc_data()``. Each -source carries its checksum and a ``provenance.json`` (collection, year, URL, -legend, license). MapBiomas code 0 ("not observed", e.g. open sea) is nodata. - - python examples/08_mapbiomas_land_use.py # temporary workspace - python examples/08_mapbiomas_land_use.py ./scratch # keep the outputs - python examples/08_mapbiomas_land_use.py --offline # no network: synthetic stand-in - -Without network access (or with ``--offline`` / ``DISSCUBE_OFFLINE=1``, as in -the test suite) the script builds a synthetic scene with the same grid, CRS -and class codes, and says so loudly: those numbers are not MapBiomas data. -Data © MapBiomas (CC-BY-4.0). -""" - -import json -import os -import sys -import tempfile -from pathlib import Path - -import numpy as np -from rasterio.crs import CRS -from rasterio.transform import from_origin - -from disscube import CubeClient, Derivation -from disscube.sources import Window2D, register_raster -from disscube.sources.mapbiomas import CLASSES, MAPBIOMAS_NODATA, register_mapbiomas_source -from disscube.utils import register_local_grid - -YEARS = (2000, 2020) -BBOX = (-44.35, -2.62, -44.20, -2.47) # São Luís, Ilha do Maranhão (WGS84) -GRID_RES = 300.0 -FOREST, MANGROVE, PASTURE, URBAN, WATER = 3, 5, 15, 24, 33 - - -# --------------------------------------------------------------------------- -# Inputs: real (MapBiomas) or synthetic (offline) -# --------------------------------------------------------------------------- - -def sources_from_mapbiomas(cube: CubeClient, raw: Path) -> None: - for year in YEARS: - register_mapbiomas_source(cube, f"lulc_{year}", year, BBOX, raw) - print(f"MapBiomas {year} : read (collection 11, 30 m)") - - -def sources_synthetic(cube: CubeClient, raw: Path) -> None: - """Sea to the north-west, a town that grows between the years, mangrove on the coast.""" - res = 0.00026949458523585647 # MapBiomas 30 m pixel, in degrees - cols = int((BBOX[2] - BBOX[0]) / res) - rows = int((BBOX[3] - BBOX[1]) / res) - transform = from_origin(BBOX[0], BBOX[3], res, res) - x = np.linspace(0, 1, cols)[None, :].repeat(rows, axis=0) - y = np.linspace(0, 1, rows)[:, None].repeat(cols, axis=1) - rng = np.random.default_rng(3) - - for year, radius in zip(YEARS, (0.18, 0.28), strict=True): - lulc = np.full((rows, cols), FOREST, dtype="float32") - lulc[rng.random((rows, cols)) < 0.15] = PASTURE - coast = x + (1 - y) - lulc[(coast >= 0.45) & (coast < 0.55)] = MANGROVE - lulc[(x - 0.6) ** 2 + (y - 0.5) ** 2 < radius ** 2] = URBAN - lulc[coast < 0.45] = WATER - lulc[coast < 0.15] = np.nan # "not observed" (open sea) - register_raster( - cube, f"lulc_{year}", Window2D(lulc, transform, CRS.from_epsg(4326)), raw, - provenance={"synthetic": True, "note": "offline stand-in, not MapBiomas data"}, - name=f"land use {year} (synthetic stand-in)", time=year, tags=["synthetic"], - dtype="uint8", nodata=MAPBIOMAS_NODATA, - ) - - -# --------------------------------------------------------------------------- -# Pipeline -# --------------------------------------------------------------------------- - -def main(workspace: Path, offline: bool) -> None: - raw = workspace / "raw" - raw.mkdir(parents=True, exist_ok=True) - cube = CubeClient(catalog=str(workspace / "catalog.db"), store=str(workspace / "store")) - grid = register_local_grid(cube, name="ilha_do_maranhao", bbox_geo=BBOX, resolution=GRID_RES) - - label = "synthetic stand-in" - done = False - if not offline: - try: - sources_from_mapbiomas(cube, raw) - label = "MapBiomas collection 11 (30 m)" - done = True - except Exception as exc: # noqa: BLE001 — any network failure falls back, loudly - print(f"!! MapBiomas not reachable ({type(exc).__name__}: {exc})") - if not done: - print("!! OFFLINE — using a synthetic scene; the numbers below are NOT MapBiomas data") - sources_synthetic(cube, raw) - - for year in YEARS: - src = f"lulc_{year}" - for d in ( - Derivation(target="landuse", source_id=src, operator="majority"), - Derivation(target="urban_pct", source_id=src, operator="percentage", class_code=URBAN), - Derivation(target="forest_pct", source_id=src, operator="percentage", class_code=FOREST), - Derivation(target="mangrove_pct", source_id=src, operator="percentage", class_code=MANGROVE), - ): - cube.derive_declarative(d, grid_id=grid.id) - - landuse = cube.load("landuse", grid_id=grid.id) # (time, y, x) - print(f"source : {label}") - print(f"grid : {grid.id} {landuse.sizes['y']} × {landuse.sizes['x']} cells, " - f"years {[int(t) for t in landuse['time'].values]}") - print("share of the area : " + " ".join(f"{y:>6}" for y in YEARS)) - for var, code in (("urban_pct", URBAN), ("forest_pct", FOREST), ("mangrove_pct", MANGROVE)): - series = cube.load(var, grid_id=grid.id) - shares = [float(np.nanmean(series.sel(time=y))) for y in YEARS] - print(f" {CLASSES[code]:<17} : " + " ".join(f"{s:6.1%}" for s in shares)) - - first, last = (landuse.sel(time=y).values for y in YEARS) - valid = ~np.isnan(first) & ~np.isnan(last) - changed = int((first[valid] != last[valid]).sum()) - to_urban = int(((first != URBAN) & (last == URBAN) & valid).sum()) - print(f"dominant class : changed in {changed} of {int(valid.sum())} cells " - f"({to_urban} became urban)") - - # Hand-off to DisSModel: one RasterBackend with the (time, y, x) series. - backend = cube.to_lucc_data(["landuse", "urban_pct"], grid_id=grid.id) - times = [int(t) for t in backend.time_coords["landuse"]] - print(f"to_lucc_data : RasterBackend {backend.shape}, landuse and urban_pct at {times}") - - prov = json.loads((raw / f"lulc_{YEARS[-1]}.provenance.json").read_text(encoding="utf-8")) - print(f"provenance : {raw / f'lulc_{YEARS[-1]}.provenance.json'} " - f"({prov.get('url', prov.get('note'))})") - - -if __name__ == "__main__": - args = [a for a in sys.argv[1:] if not a.startswith("--")] - offline = "--offline" in sys.argv or os.environ.get("DISSCUBE_OFFLINE") == "1" - if args: - main(Path(args[0]), offline) - else: - with tempfile.TemporaryDirectory() as tmp: - main(Path(tmp), offline) diff --git a/examples/09_prodes_deforestation.py b/examples/09_prodes_deforestation.py deleted file mode 100644 index bb1fba9..0000000 --- a/examples/09_prodes_deforestation.py +++ /dev/null @@ -1,116 +0,0 @@ -""" -09 — PRODES: deforestation south of Santarém (LuccME Lab15 area), three years. - -Downloads the current PRODES edition of the Legal Amazon from TerraBrasilis -(≈130 MB, once — cached in ``$DISSCUBE_CACHE`` or ``~/.cache/disscube``), reads -its legend from the ``.qml`` shipped in the ZIP, and turns the area of the -LuccME Lab15 cellular space (along the BR-163, near Mojuí dos Campos, Pará) -into forest / deforested / other at the end of 2008, 2016 and 2024. The maps -are aggregated onto a 500 m grid snapped to the BDC Albers mesh: - - - ``deforested_pct`` fraction of each cell deforested (cumulative) - - ``forest_pct`` fraction still under native forest - - python examples/09_prodes_deforestation.py # temporary workspace - python examples/09_prodes_deforestation.py ./scratch # keep the outputs - python examples/09_prodes_deforestation.py --offline # no network: synthetic stand-in - -The PRODES edition distributed today separates deforestation by year from -2007 on, so earlier years are not available here (see -``disscube.sources.prodes``). Without network access (or with ``--offline`` / -``DISSCUBE_OFFLINE=1``, as in the test suite) the script uses a synthetic -PRODES-like ZIP and says so loudly: those numbers are not PRODES data. -""" - -import json -import os -import sys -import tempfile -import zipfile -from pathlib import Path - -import numpy as np -import rasterio -from rasterio.transform import from_origin - -from disscube import CubeClient, Derivation -from disscube.sources import prodes -from disscube.utils import register_local_grid - -YEARS = (2008, 2016, 2024) -BBOX = (-54.842, -3.587, -54.459, -3.168) # LuccME Lab15 cellular space (cs_moju) -GRID_RES = 500.0 - - -def synthetic_edition(folder: Path) -> prodes.ProdesFiles: - """A PRODES-like ZIP: clearings spreading east from a road, by year.""" - folder.mkdir(parents=True, exist_ok=True) - name = "prodes_synthetic_2025" - res = 0.00026949458523585647 - cols = int((BBOX[2] - BBOX[0]) / res) + 20 - rows = int((BBOX[3] - BBOX[1]) / res) + 20 - x = np.linspace(0, 1, cols)[None, :].repeat(rows, axis=0) - rng = np.random.default_rng(11) - noise = rng.random((rows, cols)) * 0.08 - data = np.full((rows, cols), 100, dtype="uint8") # native forest - for year in range(2007, 2026): # the frontier moves east - data[(x + noise < 0.05 + 0.02 * (year - 2007)) & (data == 100)] = year - 2000 - data[:, cols - 40:] = 91 # a river on the east edge - tif = folder / f"{name}.tif" - with rasterio.open(tif, "w", driver="GTiff", height=rows, width=cols, count=1, dtype="uint8", - crs="EPSG:4674", transform=from_origin(BBOX[0] - 10 * res, BBOX[3] + 10 * res, res, res), - nodata=255) as dst: - dst.write(data, 1) - entries = [(y - 2000, f"{y - 2000} d{y}") for y in range(2007, 2026)] - entries += [(91, "91 Hidrografia"), (100, "100 Vegetação nativa florestal")] - qml = "" + "".join(f'' for v, lab in entries) + "" - zpath = folder / f"{name}.zip" - with zipfile.ZipFile(zpath, "w") as zf: - zf.write(tif, tif.name) - zf.writestr(f"{name}.qml", qml) - return prodes.download(zpath.as_uri(), cache_dir=folder / "cache") - - -def main(workspace: Path, offline: bool) -> None: - raw = workspace / "raw" - raw.mkdir(parents=True, exist_ok=True) - cube = CubeClient(catalog=str(workspace / "catalog.db"), store=str(workspace / "store")) - grid = register_local_grid(cube, name="lab15", bbox_geo=BBOX, resolution=GRID_RES) - - files, label = None, "synthetic stand-in" - if not offline: - try: - files = prodes.download() - label = f"PRODES {files.zip_path.stem}" - except Exception as exc: # noqa: BLE001 — any network failure falls back, loudly - print(f"!! PRODES not reachable ({type(exc).__name__}: {exc})") - if files is None: - print("!! OFFLINE — using a synthetic PRODES-like map; the numbers below are NOT PRODES data") - files = synthetic_edition(workspace / "synthetic") - - for year in YEARS: - prodes.register_prodes_source(cube, f"prodes_{year}", year, BBOX, raw, files=files) - for target, code in (("deforested_pct", prodes.DEFORESTED), ("forest_pct", prodes.FOREST)): - cube.derive_declarative(Derivation(target=target, source_id=f"prodes_{year}", - operator="percentage", class_code=code), grid_id=grid.id) - - defor = cube.load("deforested_pct", grid_id=grid.id) - forest = cube.load("forest_pct", grid_id=grid.id) - print(f"source : {label}") - print(f"grid : {grid.id} {defor.sizes['y']} × {defor.sizes['x']} cells") - print("share of the area : " + " ".join(f"{y:>6}" for y in YEARS)) - print(" deforested : " + " ".join(f"{float(np.nanmean(defor.sel(time=y))):6.1%}" for y in YEARS)) - print(" native forest : " + " ".join(f"{float(np.nanmean(forest.sel(time=y))):6.1%}" for y in YEARS)) - prov = json.loads((raw / f"prodes_{YEARS[-1]}.provenance.json").read_text(encoding="utf-8")) - print(f"edition : {prov['edition']} (years separable from {prov['first_separable_year']})") - print(f"provenance : {raw / f'prodes_{YEARS[-1]}.provenance.json'}") - - -if __name__ == "__main__": - args = [a for a in sys.argv[1:] if not a.startswith("--")] - offline = "--offline" in sys.argv or os.environ.get("DISSCUBE_OFFLINE") == "1" - if args: - main(Path(args[0]), offline) - else: - with tempfile.TemporaryDirectory() as tmp: - main(Path(tmp), offline) diff --git a/examples/README.md b/examples/README.md index 7fe4037..a675eb0 100644 --- a/examples/README.md +++ b/examples/README.md @@ -1,62 +1,149 @@ # DisSCube — Examples -Self-contained, runnable examples. Examples 01–03 generate their own -synthetic input data; examples 04–06 use real data bundled in -[`data/terrame/`](data/terrame/); example 07 reads the Brazil Data Cube over -the network, 08 reads MapBiomas and 09 downloads PRODES (≈130 MB, cached). All of them work in a temporary directory; 01–06 finish in a few -seconds with no downloads: +Self-contained, offline and lightweight examples illustrating the two ways of using DisSCube: +the **Python API** and declarative **TOML pipelines**. + +All examples run in a few seconds without downloading external data: ```bash pip install -e . python examples/01_quickstart.py ``` -Pass a directory to keep the catalog, the raw inputs and the derived Zarr -stores for inspection (e.g. in QGIS): +Pass a directory as argument to keep the catalog and generated Zarr stores: ```bash python examples/01_quickstart.py ./scratch ``` -| Example | What it shows | -|---|---| -| [`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`; several 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()` (including `period`) | -| [`04_terrame_fill_itaituba.py`](04_terrame_fill_itaituba.py) | TerraME's *Fill* tutorial on real data (Itaituba, Pará, 5 km): `mean`, `percentage` × `coverage_purity`, `min_distance`, compared cell by cell with TerraME's own output | -| [`05_terrame_fill_emas.py`](05_terrame_fill_emas.py) | TerraME's Emas National Park example (500 m): `presence` of lines, `max` / `min` of a raster; a study area defined by a limit polygon | -| [`06_terrame_fill_amazonia.py`](06_terrame_fill_amazonia.py) | TerraME's Brazilian Amazon example (50 km): PRODES coverage with a declared nodata, distances to roads and ports | -| [`07_bdc_cube.py`](07_bdc_cube.py) | Brazil Data Cube: Landsat 16-day cube over Ilha do Maranhão via STAC (windowed reads), dry-season median, NDVI / MNDWI / open-water drivers on a 300 m BDC Albers grid. Needs network and `.[bdc]`; `--offline` runs a synthetic stand-in | -| [`08_mapbiomas_land_use.py`](08_mapbiomas_land_use.py) | MapBiomas Collection 11 (30 m) over Ilha do Maranhão in 2000 and 2020: windowed reads of the national files, code 0 as nodata, `majority` and `percentage` (urban, forest, mangrove) on the grid of 07, `(time, y, x)` series and `to_lucc_data()`. Needs network; `--offline` runs a synthetic stand-in | -| [`09_prodes_deforestation.py`](09_prodes_deforestation.py) | PRODES in the LuccME Lab15 area (south of Santarém) at the end of 2008, 2016 and 2024: edition downloaded once and cached, legend read from its `.qml`, forest / deforested / other, `percentage` on a 500 m BDC grid. Needs network; `--offline` runs a synthetic stand-in | +--- -All examples are executed by the test suite (`tests/test_examples.py`), so -they are kept in sync with the API. The suite sets `DISSCUBE_OFFLINE=1`, so -07, 08 and 09 run on their synthetic stand-ins there. +## Two Ways to Use DisSCube -## Pipeline files +DisSCube supports two equivalent, interoperable workflows for preparing spatial data cubes: -[`pipelines/`](pipelines/) has the same preparations declared as TOML files and -run with `disscube run `: TerraME's Itaituba Fill tutorial (offline, -bundled data) and the BDC, MapBiomas and PRODES examples. See -`docs/guides/pipeline_files.md`. +````carousel +```python +# (a) Python API (CubeClient & Derivation) +# examples/01_quickstart.py -## Scope +from disscube import CubeClient, Derivation, GridSpec, SpatialSource -DisSCube prepares data for models; it stops at `CubeClient.to_lucc_data()`. -Examples that run simulations with the prepared data (BR-MANGUE, LUCC) belong -to the model repositories, where those dependencies live. +# 1. Initialize client +cube = CubeClient("catalog.db", "store") -## Exporting to GeoTIFF +# 2. Register grid (300 m cells, SIRGAS 2000 / UTM 23S) +cube.register_grid(GridSpec( + id="demo/300m", type="local", crs="EPSG:31983", resolution=300.0, + bbox=[570000.0, 9708000.0, 582000.0, 9720000.0], +)) + +# 3. Register spatial sources +cube.register_spatial_source(SpatialSource( + id="landuse", format="raster", asset_url="landuse.tif", crs="EPSG:31983", +)) +cube.register_spatial_source(SpatialSource( + id="elevation", format="raster", asset_url="elevation.tif", crs="EPSG:31983", +)) + +# 4. Declarative derivations +cube.derive_declarative( + Derivation(target="forest_pct", source_id="landuse", operator="percentage", class_code=3), + grid_id="demo/300m", +) +cube.derive_declarative( + Derivation(target="landuse_major", source_id="landuse", operator="majority"), + grid_id="demo/300m", +) +cube.derive_declarative( + Derivation(target="elev_mean", source_id="elevation", operator="mean"), + grid_id="demo/300m", +) + +# 5. Load model-ready DataArrays +forest = cube.load("forest_pct", grid_id="demo/300m") +``` + +```toml +# (b) Declarative Pipeline (TOML) +# examples/pipelines/quickstart.toml -Derived data cubes can be exported to multi-band GeoTIFFs using the CLI: +schema = 1 +name = "Quickstart — Synthetic Land Use and Elevation" + +[grid] +name = "demo/300m" +crs = "EPSG:31983" +bbox = [570000.0, 9708000.0, 582000.0, 9720000.0] +resolution = 300 + +[[source]] +id = "landuse" +type = "file" +path = "../data/quickstart/landuse.tif" +crs = "EPSG:31983" + +[[source]] +id = "elevation" +type = "file" +path = "../data/quickstart/elevation.tif" +crs = "EPSG:31983" + +[[derive]] +target = "forest_pct" +source = "landuse" +operator = "percentage" +class_code = 3 + +[[derive]] +target = "landuse_major" +source = "landuse" +operator = "majority" + +[[derive]] +target = "elev_mean" +source = "elevation" +operator = "mean" +``` +```` +### Running the TOML Pipeline + +Validate without I/O or downloads: ```bash -disscube export examples/pipelines/itaituba_fill.toml --output data/cellspace.tif +disscube validate examples/pipelines/quickstart.toml ``` -Or directly during execution: +Execute and record provenance: +```bash +disscube run examples/pipelines/quickstart.toml --workspace outputs/quickstart +``` +Export directly to multi-band GeoTIFF: ```bash -disscube run examples/pipelines/itaituba_fill.toml --output data/cellspace.tif +disscube export examples/pipelines/quickstart.toml --output outputs/cellspace.tif ``` + +--- + +## Python API Examples + +| Example | What it shows | +|---|---| +| [`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()` | + +All examples are tested automatically on CI (`tests/test_examples.py`). + +--- + +## Real-World Cases and Recipes + +Real data workflows and parity benchmarks are maintained in the dedicated repository +[**LambdaGeo/disscube-recipes**](https://github.com/LambdaGeo/disscube-recipes): + +- `cases/terrame_fill`: Cell-by-cell numerical parity against TerraME's C++ `fillCellularSpace` (Itaituba, Emas, Amazônia). +- `cases/ilha_maranhao`: BDC Landsat-16D (STAC) and MapBiomas land-cover time series on a 300 m BDC Albers grid. +- `cases/prodes_br163`: Multi-year PRODES deforestation along the BR-163 corridor (Pará). +- `cases/luccme_br`: National LUCC reconstruction for Brazil (BigEarth). diff --git a/examples/data/quickstart/elevation.tif b/examples/data/quickstart/elevation.tif new file mode 100644 index 0000000..af24741 Binary files /dev/null and b/examples/data/quickstart/elevation.tif differ diff --git a/examples/data/quickstart/landuse.tif b/examples/data/quickstart/landuse.tif new file mode 100644 index 0000000..91e3ba9 Binary files /dev/null and b/examples/data/quickstart/landuse.tif differ diff --git a/examples/data/terrame/README.md b/examples/data/terrame/README.md deleted file mode 100644 index 4ab8816..0000000 --- a/examples/data/terrame/README.md +++ /dev/null @@ -1,72 +0,0 @@ -# TerraME *Fill* datasets - -Input layers and reference outputs of the three *Fill* examples shipped with -TerraME's `gis` package. They are used by examples -[`04_terrame_fill_itaituba.py`](../../04_terrame_fill_itaituba.py), -[`05_terrame_fill_emas.py`](../../05_terrame_fill_emas.py) and -[`06_terrame_fill_amazonia.py`](../../06_terrame_fill_amazonia.py), -and by the parity tests (`tests/test_terrame_parity.py`), which compare -DisSCube cell by cell with the cellular spaces TerraME itself produced. - -## Source - -- Repository: , directory `packages/gis/data/` -- Commit: `84faddc4e5149a89d7486444b61b8c50474e199c` (2023-08-08) -- License: GNU LGPL — version 3 according to the repository's `LICENSE`; - the headers of the `.lua` scripts say version 2.1. © INPE and TerraLAB/UFOP. - These files are redistributed unmodified with this attribution and are - **not** covered by DisSCube's MIT license. - -## Layout - -Each shapefile is stored as one zip (`.shp`, `.shx`, `.dbf`, `.prj`) and is -read directly, without unpacking, through a `zip://` path, e.g. -`zip:///…/itaituba/itaituba-roads.zip`. The members are byte-identical to -the files in the TerraME repository; the spatial-index `.qix` files and the -`.tme` window layouts were left out. Rasters are the original GeoTIFFs. Each -directory also keeps TerraME's script (`.lua`) that produced the -reference cells. - -| Dataset | TerraME script | Reference cells | Resolution | CRS | Fill operations | -|---|---|---|---|---|---| -| `itaituba/` | `itaituba.lua` | `itaituba.zip` (620) | 5 km | EPSG:29191 | `average`, `coverage`, `distance` (lines, points), `sum` with `area = true` | -| `emas/` | `emas.lua` | `emas.zip` (5 514) | 500 m | EPSG:29192 | `presence` (lines), `maximum`, `minimum` | -| `amazonia/` | `amazonia.lua` | `amazonia.zip` (2 229) | 50 km | EPSG:29191 | `coverage`, `distance` (lines, points), `area` | - -## Notes for reproducing TerraME's grids - -- **The `.prj` of the reference cells is wrong.** `itaituba` and `amazonia` - declare `NAD83_Austin` (a Lambert conformal projection for Texas) and - `emas` a generic UTM 22S; the cells are actually in the CRS of the input - layers (EPSG:29191 or EPSG:29192), which must be assigned explicitly. -- **Grid.** Cells are 5 km / 500 m / 50 km squares whose extent is the total - bounds of the reference layer. With `input = "limit"` (Emas, Amazônia) - TerraME keeps only the cells that touch the study-area limit, so the - reference is a subset of the full rectangle; match cells by centroid. -- **Rows.** TerraME numbers rows from the south (`row = 0` is the southern - row); DisSCube grids are north-up. - -## Checksums (SHA-256) - -| File | SHA-256 | -|---|---| -| `itaituba/itaituba-census.zip` | `24f1c5afcf4b17a1c9b4d30a721067354dc56697654c754f309471eb55a605f8` | -| `itaituba/itaituba-deforestation.tif` | `b31c9b9cf492e56609a859fa82566066aec34f51c977c922f341f1ef1be5dbf6` | -| `itaituba/itaituba-elevation.tif` | `f9fac6f61ed6ac29c1e543c8c6cfc239148eb6b308b2bcb0ff915f42abe2879c` | -| `itaituba/itaituba-localities.zip` | `0f50cd16bf2ff4c9aaa6347484eef9e848c5f2285741d199461db9d6831bf00e` | -| `itaituba/itaituba-roads.zip` | `f15b8ccc4ce484eb42393519e333607a8d93e4cbfa73055e3c1d775d4cf66749` | -| `itaituba/itaituba.lua` | `6ac758ce249779ee26f613f46098d557ffa7f382e0598263220c13bf079d8060` | -| `itaituba/itaituba.zip` | `76ac2b17b8aafbc90a0d7e3939e4602adbe76331c720f470009d255f09d09856` | -| `emas/emas-accumulation.tif` | `d3ef918861c36c46a7333d79346dbb54a3e48b80468787715ddd5ef317e2e58a` | -| `emas/emas-firebreak.zip` | `618859dbf0b7a3a22919d7180c7a7f3026712c2dc2eca2aed700d25acd5f77ba` | -| `emas/emas-limit.zip` | `66aaffc986d5ac2e6d14e543f0aaa27ab02fdf849005deaa6b73efbd3486e652` | -| `emas/emas-river.zip` | `f58b31befd67e716bf9a5c79502988b470b986673dd9e586924177f8ea4c6ab6` | -| `emas/emas.lua` | `b3fdfab361c9007691c893a7fb7de4ae33cee557fe26aa41b64ad1c64793805a` | -| `emas/emas.zip` | `cce329dc25d5d6243a0233b7d8adb993b94a637ac167c6b4bd553277dcb82483` | -| `amazonia/amazonia-indigenous.zip` | `46674013fd6b3bcfd9fee46595913a1d80ded53f4d0c82f5510ed332aa0c3574` | -| `amazonia/amazonia-limit.zip` | `d44ac5514302a09005b39215e859b53b353a08da2267fb841959f7dfc3e2674b` | -| `amazonia/amazonia-ports.zip` | `c7e537ca771ad1d7cfa0b17a050179ed5686825b685615f34d7304d975628b35` | -| `amazonia/amazonia-prodes.tif` | `28fef2ada04460ce05abf1ce40582cf08c928e6bc766610473c5ff7968ab365a` | -| `amazonia/amazonia-roads.zip` | `a5f503e80acb3328d9dd815125acbf32f2a9320a32438af748e28bf3cbfd3495` | -| `amazonia/amazonia.lua` | `0bbd3d0885fff24e58914c6205ac649347c0385a3f95c05360984c727dd34477` | -| `amazonia/amazonia.zip` | `97c01fa6efbee8e0775054176728e023e86b7ad47e31cf496ab944674f561893` | diff --git a/examples/data/terrame/amazonia/amazonia-indigenous.zip b/examples/data/terrame/amazonia/amazonia-indigenous.zip deleted file mode 100644 index 592abf7..0000000 Binary files a/examples/data/terrame/amazonia/amazonia-indigenous.zip and /dev/null differ diff --git a/examples/data/terrame/amazonia/amazonia-limit.zip b/examples/data/terrame/amazonia/amazonia-limit.zip deleted file mode 100644 index a4193ee..0000000 Binary files a/examples/data/terrame/amazonia/amazonia-limit.zip and /dev/null differ diff --git a/examples/data/terrame/amazonia/amazonia-ports.zip b/examples/data/terrame/amazonia/amazonia-ports.zip deleted file mode 100644 index 8b93fd8..0000000 Binary files a/examples/data/terrame/amazonia/amazonia-ports.zip and /dev/null differ diff --git a/examples/data/terrame/amazonia/amazonia-prodes.tif b/examples/data/terrame/amazonia/amazonia-prodes.tif deleted file mode 100644 index ccefa19..0000000 Binary files a/examples/data/terrame/amazonia/amazonia-prodes.tif and /dev/null differ diff --git a/examples/data/terrame/amazonia/amazonia-roads.zip b/examples/data/terrame/amazonia/amazonia-roads.zip deleted file mode 100644 index 9d6ed35..0000000 Binary files a/examples/data/terrame/amazonia/amazonia-roads.zip and /dev/null differ diff --git a/examples/data/terrame/amazonia/amazonia.lua b/examples/data/terrame/amazonia/amazonia.lua deleted file mode 100644 index 6b5616f..0000000 --- a/examples/data/terrame/amazonia/amazonia.lua +++ /dev/null @@ -1,130 +0,0 @@ -------------------------------------------------------------------------------------------- --- TerraME - a software platform for multiple scale spatially-explicit dynamic modeling. --- Copyright (C) 2001-2017 INPE and TerraLAB/UFOP -- www.terrame.org - --- This code is part of the TerraME framework. --- This framework is free software; you can redistribute it and/or --- modify it under the terms of the GNU Lesser General Public --- License as published by the Free Software Foundation; either --- version 2.1 of the License, or (at your option) any later version. - --- You should have received a copy of the GNU Lesser General Public --- License along with this library. - --- The authors reassure the license terms regarding the warranties. --- They specifically disclaim any warranties, including, but not limited to, --- the implied warranties of merchantability and fitness for a particular purpose. --- The framework provided hereunder is on an "as is" basis, and the authors have no --- obligation to provide maintenance, support, updates, enhancements, or modifications. --- In no event shall INPE and TerraLAB / UFOP be held liable to any party for direct, --- indirect, special, incidental, or consequential damages arising out of the use --- of this software and its documentation. --- -------------------------------------------------------------------------------------------- - --- @example Creates a database that can be used by the example deforestation of base package. - -import("gis") - -amazonia = Project{ - file = "amazonia.qgs", - clean = true, - title = "Amazonia database", - ports = filePath("amazonia-ports.shp", "gis"), - roads = filePath("amazonia-roads.shp", "gis"), - limit = filePath("amazonia-limit.shp", "gis") -} - -prodes = Layer{ - name = "prodes", - project = amazonia, - epsg = 29191, - file = filePath("amazonia-prodes.tif", "gis") -} - -protected = Layer{ - name = "protected", - project = amazonia, - epsg = 29191, - file = filePath("amazonia-indigenous.shp", "gis") -} - -amazoniaCells = Layer{ - project = amazonia, - file = "amazonia.shp", - clean = true, - input = "limit", - name = "cells", - resolution = 50000, - progress = false -} - -amazoniaCells:fill{ - operation = "coverage", - layer = "prodes", - attribute = "prodes", - progress = false -} - -amazoniaCells:fill{ - operation = "distance", - layer = "roads", - attribute = "distroads", - progress = false -} - -amazoniaCells:fill{ - operation = "distance", - layer = "ports", - attribute = "distports", - progress = false -} - -amazoniaCells:fill{ - operation = "area", - layer = "protected", - attribute = "protected", - progress = false -} - -cs = CellularSpace{ - project = amazonia, - layer = "cells", - as = { - forest = "prodes_208", - deforestation = "prodes_10" - } -} - -Map{ - target = cs, - select = "distroads", - slices = 10, - invert = true, - color = "YlOrBr" -} - -Map{ - target = cs, - select = "distports", - slices = 10, - invert = true, - color = "YlOrBr" -} - -Map{ - target = cs, - select = "protected", - slices = 10, - invert = true, - color = "PuBu" -} - -Map{ - target = cs, - select = "deforestation", - slices = 10, - invert = true, - color = "YlGn" -} - diff --git a/examples/data/terrame/amazonia/amazonia.zip b/examples/data/terrame/amazonia/amazonia.zip deleted file mode 100644 index 39234e4..0000000 Binary files a/examples/data/terrame/amazonia/amazonia.zip and /dev/null differ diff --git a/examples/data/terrame/emas/emas-accumulation.tif b/examples/data/terrame/emas/emas-accumulation.tif deleted file mode 100644 index 271cf17..0000000 Binary files a/examples/data/terrame/emas/emas-accumulation.tif and /dev/null differ diff --git a/examples/data/terrame/emas/emas-firebreak.zip b/examples/data/terrame/emas/emas-firebreak.zip deleted file mode 100644 index 0ea7ec4..0000000 Binary files a/examples/data/terrame/emas/emas-firebreak.zip and /dev/null differ diff --git a/examples/data/terrame/emas/emas-limit.zip b/examples/data/terrame/emas/emas-limit.zip deleted file mode 100644 index 1297aa1..0000000 Binary files a/examples/data/terrame/emas/emas-limit.zip and /dev/null differ diff --git a/examples/data/terrame/emas/emas-river.zip b/examples/data/terrame/emas/emas-river.zip deleted file mode 100644 index fe522b1..0000000 Binary files a/examples/data/terrame/emas/emas-river.zip and /dev/null differ diff --git a/examples/data/terrame/emas/emas.lua b/examples/data/terrame/emas/emas.lua deleted file mode 100644 index daecfaa..0000000 --- a/examples/data/terrame/emas/emas.lua +++ /dev/null @@ -1,118 +0,0 @@ -------------------------------------------------------------------------------------------- --- TerraME - a software platform for multiple scale spatially-explicit dynamic modeling. --- Copyright (C) 2001-2017 INPE and TerraLAB/UFOP -- www.terrame.org - --- This code is part of the TerraME framework. --- This framework is free software; you can redistribute it and/or --- modify it under the terms of the GNU Lesser General Public --- License as published by the Free Software Foundation; either --- version 2.1 of the License, or (at your option) any later version. - --- You should have received a copy of the GNU Lesser General Public --- License along with this library. - --- The authors reassure the license terms regarding the warranties. --- They specifically disclaim any warranties, including, but not limited to, --- the implied warranties of merchantability and fitness for a particular purpose. --- The framework provided hereunder is on an "as is" basis, and the authors have no --- obligation to provide maintenance, support, updates, enhancements, or modifications. --- In no event shall INPE and TerraLAB / UFOP be held liable to any party for direct, --- indirect, special, incidental, or consequential damages arising out of the use --- of this software and its documentation. --- -------------------------------------------------------------------------------------------- - --- @example Creates a database that can be used by the example fire-spread of base package. - -import("gis") - -project = Project{ - file = "emas.qgs", - clean = true, - author = "Almeida, R.", - title = "Emas database", - firebreak = filePath("emas-firebreak.shp", "gis"), - river = filePath("emas-river.shp", "gis"), - limit = filePath("emas-limit.shp", "gis") -} - -Layer{ - project = project, - name = "cover", - file = filePath("emas-accumulation.tif", "gis"), - epsg = 29192 -} - -cl = Layer{ - project = project, - file = "emas.shp", - clean = true, - input = "limit", - name = "cells", - resolution = 500, - progress = false -} - -cl:fill{ - operation = "presence", - attribute = "firebreak", - layer = "firebreak", - progress = false -} - -cl:fill{ - operation = "presence", - attribute = "river", - layer = "river", - progress = false -} - -cl:fill{ - operation = "maximum", - attribute = "maxcover", - layer = "cover", - progress = false -} - -cl:fill{ - operation = "minimum", - attribute = "mincover", - layer = "cover", - progress = false -} - -cs = CellularSpace{ - project = project, - layer = "cells" -} - -Map{ - target = cs, - select = "firebreak", - value = {0, 1}, - color = {"darkGreen", "black"}, - label = {"forest", "firebreak"} -} - -Map{ - target = cs, - select = "river", - value = {0, 1}, - color = {"darkGreen", "darkBlue"}, - label = {"forest", "river"} -} - -Map{ - target = cs, - select = "mincover", - slices = 5, - color = "YlGn" -} - -Map{ - target = cs, - select = "maxcover", - slices = 5, - color = "YlGn" -} - diff --git a/examples/data/terrame/emas/emas.zip b/examples/data/terrame/emas/emas.zip deleted file mode 100644 index 94468af..0000000 Binary files a/examples/data/terrame/emas/emas.zip and /dev/null differ diff --git a/examples/data/terrame/itaituba/itaituba-census.zip b/examples/data/terrame/itaituba/itaituba-census.zip deleted file mode 100644 index cdc21b1..0000000 Binary files a/examples/data/terrame/itaituba/itaituba-census.zip and /dev/null differ diff --git a/examples/data/terrame/itaituba/itaituba-deforestation.tif b/examples/data/terrame/itaituba/itaituba-deforestation.tif deleted file mode 100644 index b058c0a..0000000 Binary files a/examples/data/terrame/itaituba/itaituba-deforestation.tif and /dev/null differ diff --git a/examples/data/terrame/itaituba/itaituba-elevation.tif b/examples/data/terrame/itaituba/itaituba-elevation.tif deleted file mode 100644 index 5453971..0000000 Binary files a/examples/data/terrame/itaituba/itaituba-elevation.tif and /dev/null differ diff --git a/examples/data/terrame/itaituba/itaituba-localities.zip b/examples/data/terrame/itaituba/itaituba-localities.zip deleted file mode 100644 index 6676322..0000000 Binary files a/examples/data/terrame/itaituba/itaituba-localities.zip and /dev/null differ diff --git a/examples/data/terrame/itaituba/itaituba-roads.zip b/examples/data/terrame/itaituba/itaituba-roads.zip deleted file mode 100644 index 8111e1f..0000000 Binary files a/examples/data/terrame/itaituba/itaituba-roads.zip and /dev/null differ diff --git a/examples/data/terrame/itaituba/itaituba.lua b/examples/data/terrame/itaituba/itaituba.lua deleted file mode 100644 index 66ccf7b..0000000 --- a/examples/data/terrame/itaituba/itaituba.lua +++ /dev/null @@ -1,156 +0,0 @@ -------------------------------------------------------------------------------------------- --- TerraME - a software platform for multiple scale spatially-explicit dynamic modeling. --- Copyright (C) 2001-2017 INPE and TerraLAB/UFOP -- www.terrame.org - --- This code is part of the TerraME framework. --- This framework is free software; you can redistribute it and/or --- modify it under the terms of the GNU Lesser General Public --- License as published by the Free Software Foundation; either --- version 2.1 of the License, or (at your option) any later version. - --- You should have received a copy of the GNU Lesser General Public --- License along with this library. - --- The authors reassure the license terms regarding the warranties. --- They specifically disclaim any warranties, including, but not limited to, --- the implied warranties of merchantability and fitness for a particular purpose. --- The framework provided hereunder is on an "as is" basis, and the authors have no --- obligation to provide maintenance, support, updates, enhancements, or modifications. --- In no event shall INPE and TerraLAB / UFOP be held liable to any party for direct, --- indirect, special, incidental, or consequential damages arising out of the use --- of this software and its documentation. --- -------------------------------------------------------------------------------------------- - --- @example Creates a project for the Itaituba region, in Para state, Brazil. - -import("gis") - -itaituba = Project{ - file = "itaituba.tview", - clean = true, - localities = filePath("itaituba-localities.shp", "gis"), - roads = filePath("itaituba-roads.shp", "gis"), - census = filePath("itaituba-census.shp", "gis") -} - -Layer{ - project = itaituba, - name = "deforestation", - file = filePath("itaituba-deforestation.tif", "gis"), - epsg = 29191 -} - -Layer{ - project = itaituba, - name = "elevation", - file = filePath("itaituba-elevation.tif", "gis"), - epsg = 29191 -} - -itaitubaCells = Layer{ - project = itaituba, - name = "cells", - clean = true, - file = "itaituba.shp", - input = "census", - resolution = 5000, - progress = false -} - -itaitubaCells:fill{ - operation = "average", - layer = "elevation", - attribute = "elevation", - progress = false -} - -itaitubaCells:fill{ - operation = "coverage", - layer = "deforestation", - attribute = "defor", - progress = false -} - -itaitubaCells:fill{ - operation = "distance", - layer = "roads", - attribute = "distroad", - progress = false -} - -itaitubaCells:fill{ - operation = "distance", - layer = "localities", - attribute = "distlocal", - progress = false -} - -itaitubaCells:fill{ - operation = "sum", - layer = "census", - attribute = "population", - select = "population", - area = true, - progress = false -} - -cell = Cell{ - logpop = function(self) - return math.log(self.population) - end -} - -cs = CellularSpace{ - project = itaituba, - layer = "cells", - as = { - river = "defor_167", - deforestation = "defor_87" - }, - instance = cell -} - -Map{ - target = cs, - select = "elevation", - slices = 6, - color = "Blues" -} - -Map{ - target = cs, - select = "river", - slices = 6, - color = "Blues" -} - -Map{ - target = cs, - select = "deforestation", - slices = 6, - invert = true, - color = "Greens" -} - -Map{ - target = cs, - select = "distlocal", - slices = 10, - color = "Reds" -} - -Map{ - target = cs, - select = "distroad", - slices = 10, - color = "Reds" -} - -Map{ - target = cs, - select = "logpop", - slices = 10, - color = "Purples" -} - diff --git a/examples/data/terrame/itaituba/itaituba.zip b/examples/data/terrame/itaituba/itaituba.zip deleted file mode 100644 index cf290bf..0000000 Binary files a/examples/data/terrame/itaituba/itaituba.zip and /dev/null differ diff --git a/examples/pipelines/ilha_bdc.toml b/examples/pipelines/ilha_bdc.toml deleted file mode 100644 index c3c9a51..0000000 --- a/examples/pipelines/ilha_bdc.toml +++ /dev/null @@ -1,37 +0,0 @@ -# Example 07 as a pipeline file: Brazil Data Cube drivers over Ilha do Maranhão. -# Needs network and `pip install -e ".[bdc]"`. -# -# disscube run examples/pipelines/ilha_bdc.toml --workspace outputs/ilha_bdc - -schema = 1 -name = "Ilha do Maranhão — BDC Landsat drivers, dry season 2020" - -[grid] -name = "ilha_do_maranhao" -bbox = [-44.35, -2.62, -44.20, -2.47] # WGS84; grid in BDC Albers, snapped -resolution = 300 - -[[source]] -id = "ndvi" -type = "bdc" -collection = "LANDSAT-16D-1" -asset = "NDVI" -period = "2020-07-01/2020-09-30" -scale = 1e-4 # BDC stores indices as int16 × 10 000 - -[[source]] -id = "mndwi" -type = "bdc" -collection = "LANDSAT-16D-1" -normalized_difference = ["green", "swir16"] -period = "2020-07-01/2020-09-30" - -[[derive]] -target = "ndvi_mean" -source = "ndvi" -operator = "mean" - -[[derive]] -target = "mndwi_mean" -source = "mndwi" -operator = "mean" diff --git a/examples/pipelines/ilha_mapbiomas.toml b/examples/pipelines/ilha_mapbiomas.toml deleted file mode 100644 index 37820f4..0000000 --- a/examples/pipelines/ilha_mapbiomas.toml +++ /dev/null @@ -1,42 +0,0 @@ -# Example 08 as a pipeline file: MapBiomas land use of Ilha do Maranhão, 2000 and 2020. -# Needs network. -# -# disscube run examples/pipelines/ilha_mapbiomas.toml --workspace outputs/ilha_mapbiomas - -schema = 1 -name = "Ilha do Maranhão — MapBiomas land use" - -[grid] -name = "ilha_do_maranhao" -bbox = [-44.35, -2.62, -44.20, -2.47] -resolution = 300 - -[[source]] -id = "lulc_{year}" -type = "mapbiomas" -years = [2000, 2020] # one source per year: lulc_2000, lulc_2020 -collection = 11 -resolution = 30 - -[[derive]] -target = "landuse" -source = "lulc_{year}" # expanded over the source's years -operator = "majority" - -[[derive]] -target = "urban_pct" -source = "lulc_{year}" -operator = "percentage" -class_code = 24 - -[[derive]] -target = "forest_pct" -source = "lulc_{year}" -operator = "percentage" -class_code = 3 - -[[derive]] -target = "mangrove_pct" -source = "lulc_{year}" -operator = "percentage" -class_code = 5 diff --git a/examples/pipelines/itaituba_fill.toml b/examples/pipelines/itaituba_fill.toml deleted file mode 100644 index da3f64e..0000000 --- a/examples/pipelines/itaituba_fill.toml +++ /dev/null @@ -1,73 +0,0 @@ -# TerraME's Fill tutorial (Itaituba, Pará) as a DisSCube pipeline file. -# -# The same cellular space and layers as examples/data/terrame/itaituba/itaituba.lua, -# runnable offline with the bundled data: -# -# disscube validate examples/pipelines/itaituba_fill.toml -# disscube run examples/pipelines/itaituba_fill.toml --workspace outputs/itaituba -# -# TerraME DisSCube -# itaitubaCells:fill{operation = "average", [[derive]] operator = "mean" -# itaitubaCells:fill{operation = "coverage", [[derive]] operator = "percentage", class_code = … -# itaitubaCells:fill{operation = "distance", [[derive]] operator = "min_distance" -# (the census "sum, area = true" fill has no DisSCube operator yet) - -schema = 1 -name = "TerraME Fill tutorial — Itaituba" - -[grid] -name = "itaituba/5km" -crs = "EPSG:29191" # SAD69 / UTM 21S, as in the Lua script -bbox = [547177.3501989994, 9485214.18712, 702177.3501989994, 9585214.18712] -resolution = 5000 - -[[source]] -id = "elevation" -type = "file" -path = "../data/terrame/itaituba/itaituba-elevation.tif" -crs = "EPSG:29191" - -[[source]] -id = "deforestation" -type = "file" -path = "../data/terrame/itaituba/itaituba-deforestation.tif" -crs = "EPSG:29191" - -[[source]] -id = "roads" -type = "file" -path = "zip://../data/terrame/itaituba/itaituba-roads.zip" -crs = "EPSG:29191" - -[[source]] -id = "localities" -type = "file" -path = "zip://../data/terrame/itaituba/itaituba-localities.zip" -crs = "EPSG:29191" - -[[derive]] -target = "elevation" -source = "elevation" -operator = "mean" - -[[derive]] -target = "defor_87" # TerraME's "deforestation" class -source = "deforestation" -operator = "percentage" -class_code = 87 - -[[derive]] -target = "defor_167" # TerraME's "river" class -source = "deforestation" -operator = "percentage" -class_code = 167 - -[[derive]] -target = "distroad" -source = "roads" -operator = "min_distance" - -[[derive]] -target = "distlocal" -source = "localities" -operator = "min_distance" diff --git a/examples/pipelines/lab15_prodes.toml b/examples/pipelines/lab15_prodes.toml deleted file mode 100644 index 9e36094..0000000 --- a/examples/pipelines/lab15_prodes.toml +++ /dev/null @@ -1,29 +0,0 @@ -# Example 09 as a pipeline file: PRODES in the LuccME Lab15 area (south of Santarém). -# Downloads the PRODES edition once (~130 MB, cached in ~/.cache/disscube/prodes). -# -# disscube run examples/pipelines/lab15_prodes.toml --workspace outputs/lab15_prodes - -schema = 1 -name = "LuccME Lab15 area — PRODES deforestation" - -[grid] -name = "lab15" -bbox = [-54.842, -3.587, -54.459, -3.168] -resolution = 500 - -[[source]] -id = "prodes_{year}" -type = "prodes" -years = [2008, 2016, 2024] # the current edition separates years from 2007 - -[[derive]] -target = "deforested_pct" -source = "prodes_{year}" -operator = "percentage" -class_code = 2 # 1 forest, 2 deforested, 3 other - -[[derive]] -target = "forest_pct" -source = "prodes_{year}" -operator = "percentage" -class_code = 1 diff --git a/examples/pipelines/quickstart.toml b/examples/pipelines/quickstart.toml new file mode 100644 index 0000000..71fcaa0 --- /dev/null +++ b/examples/pipelines/quickstart.toml @@ -0,0 +1,48 @@ +# DisSCube — Quickstart pipeline (synthetic land use and elevation) +# +# Equivalent to examples/01_quickstart.py: +# - percentage: fraction of forest (class 3) in each 300 m cell +# - majority: dominant land-use class in each 300 m cell +# - mean: mean elevation in each 300 m cell +# +# Runs offline without external downloads: +# +# disscube validate examples/pipelines/quickstart.toml +# disscube run examples/pipelines/quickstart.toml --workspace outputs/quickstart + +schema = 1 +name = "Quickstart — Synthetic Land Use and Elevation" + +[grid] +name = "demo/300m" +crs = "EPSG:31983" +bbox = [570000.0, 9708000.0, 582000.0, 9720000.0] +resolution = 300 + +[[source]] +id = "landuse" +type = "file" +path = "../data/quickstart/landuse.tif" +crs = "EPSG:31983" + +[[source]] +id = "elevation" +type = "file" +path = "../data/quickstart/elevation.tif" +crs = "EPSG:31983" + +[[derive]] +target = "forest_pct" +source = "landuse" +operator = "percentage" +class_code = 3 + +[[derive]] +target = "landuse_major" +source = "landuse" +operator = "majority" + +[[derive]] +target = "elev_mean" +source = "elevation" +operator = "mean" diff --git a/examples/pipelines/test_aeroporto.toml b/tests/fixtures/pipelines/test_aeroporto.toml similarity index 100% rename from examples/pipelines/test_aeroporto.toml rename to tests/fixtures/pipelines/test_aeroporto.toml diff --git a/tests/test_examples.py b/tests/test_examples.py index 232a390..4d5800f 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -24,7 +24,7 @@ def test_example_runs(script, tmp_path): result = subprocess.run( [sys.executable, str(script), str(tmp_path / "workspace")], cwd=tmp_path, capture_output=True, text=True, timeout=300, check=False, - # examples that download data (07_bdc_cube.py) use their offline stand-in + # examples that support offline mode use DISSCUBE_OFFLINE=1 env={**os.environ, "DISSCUBE_OFFLINE": "1"}, ) assert result.returncode == 0, f"{script.name} failed:\n{result.stdout}\n{result.stderr}" diff --git a/tests/test_pipeline_file.py b/tests/test_pipeline_file.py index bcb2332..10caf18 100644 --- a/tests/test_pipeline_file.py +++ b/tests/test_pipeline_file.py @@ -16,14 +16,14 @@ from rasterio.crs import CRS from rasterio.transform import from_origin -from disscube import CubeClient, GridSpec, SpatialDerivation, SpatialSource, Variable +from disscube import CubeClient, Derivation, GridSpec, SpatialSource from disscube.cli import main as cli from disscube.pipeline import PipelineError, load, plan, run from disscube.sources import Window2D ROOT = Path(__file__).resolve().parents[1] PIPELINES = sorted((ROOT / "examples" / "pipelines").glob("*.toml")) -ITAITUBA = ROOT / "examples" / "pipelines" / "itaituba_fill.toml" +QUICKSTART = ROOT / "examples" / "pipelines" / "quickstart.toml" GRID = """ schema = 1 @@ -53,52 +53,52 @@ def _raster(path, value=1.0, dtype="float32"): # --------------------------------------------------------------------------- def test_example_pipelines_are_discovered(): - assert len(PIPELINES) >= 4 + assert len(PIPELINES) >= 1 @pytest.mark.parametrize("path", PIPELINES, ids=lambda p: p.name) def test_example_pipelines_validate(path): p = plan(path) - # Every example declares sources; a sources-only file (e.g. test_aeroporto.toml, + # Every example declares sources; a sources-only file (e.g. tests/fixtures/pipelines/test_aeroporto.toml, # which just fetches/validates one asset) is legitimate and has no [[derive]]. assert p.sources -def test_itaituba_pipeline_matches_the_python_api(tmp_path): - report = run(ITAITUBA, workspace=tmp_path / "toml") - assert len(report.derived) == 5 +def test_quickstart_pipeline_matches_the_python_api(tmp_path): + report = run(QUICKSTART, workspace=tmp_path / "toml") + assert len(report.derived) == 3 - data = ROOT / "examples" / "data" / "terrame" / "itaituba" + data = ROOT / "examples" / "data" / "quickstart" cube = CubeClient(catalog=str(tmp_path / "api.db"), store=str(tmp_path / "api")) - cfg = load(ITAITUBA).config.grid + cfg = load(QUICKSTART).config.grid cube.register_grid(GridSpec(id="api", type="local", crs=cfg.crs, resolution=cfg.resolution, bbox=cfg.bbox)) - cube.register_spatial_source(SpatialSource(id="elev", name="e", format="raster", crs=cfg.crs, - asset_url=str(data / "itaituba-elevation.tif"))) - cube.register_spatial_source(SpatialSource(id="roads", name="r", format="vector", crs=cfg.crs, - asset_url=f"zip://{data / 'itaituba-roads.zip'}")) - cube.derive(SpatialDerivation(source_id="elev", grid_id="api", role="driver", - variables=[Variable(name="elevation", operator="mean")])) - cube.derive(SpatialDerivation(source_id="roads", grid_id="api", role="driver", - variables=[Variable(name="distroad", operator="min_distance")])) + cube.register_spatial_source(SpatialSource(id="landuse", name="lu", format="raster", crs=cfg.crs, + asset_url=str(data / "landuse.tif"))) + cube.register_spatial_source(SpatialSource(id="elevation", name="el", format="raster", crs=cfg.crs, + asset_url=str(data / "elevation.tif"))) + cube.derive_declarative(Derivation(target="forest_pct", source_id="landuse", operator="percentage", class_code=3), + grid_id="api") + cube.derive_declarative(Derivation(target="elev_mean", source_id="elevation", operator="mean"), + grid_id="api") toml_cube = CubeClient(catalog=str(tmp_path / "toml" / "catalog.db"), store=str(tmp_path / "toml" / "store")) - for name in ("elevation", "distroad"): + for name in ("forest_pct", "elev_mean"): a = toml_cube.load(name, grid_id=report.grid_id).values b = cube.load(name, grid_id="api").values assert np.allclose(a, b, equal_nan=True), name def test_run_writes_record_and_annotates_provenance(tmp_path): - report = run(ITAITUBA, workspace=tmp_path / "ws") + report = run(QUICKSTART, workspace=tmp_path / "ws") record = json.loads(report.record.read_text()) - checksum = load(ITAITUBA).checksum + checksum = load(QUICKSTART).checksum assert record["pipeline"]["checksum"] == checksum - assert {d["target"] for d in record["derived"]} == {"elevation", "defor_87", "defor_167", - "distroad", "distlocal"} - prov = json.loads((tmp_path / "ws" / "raw" / "roads.provenance.json").read_text()) + assert {d["target"] for d in record["derived"]} == {"forest_pct", "landuse_major", "elev_mean"} + prov = json.loads((tmp_path / "ws" / "raw" / "landuse.provenance.json").read_text()) assert prov["pipeline"]["checksum"] == checksum and prov["checksum"].startswith("sha256:") + # --------------------------------------------------------------------------- # Years expansion # --------------------------------------------------------------------------- @@ -257,9 +257,9 @@ def fake_composite(collection, asset, bbox, period, **kw): # --------------------------------------------------------------------------- def test_cli_validate_and_run(tmp_path, capsys): - assert cli(["validate", str(ITAITUBA)]) == 0 + assert cli(["validate", str(QUICKSTART)]) == 0 assert "OK" in capsys.readouterr().out - assert cli(["run", str(ITAITUBA), "--workspace", str(tmp_path / "ws")]) == 0 + assert cli(["run", str(QUICKSTART), "--workspace", str(tmp_path / "ws")]) == 0 assert (tmp_path / "ws" / "run.json").exists() @@ -270,12 +270,12 @@ def test_cli_reports_errors(tmp_path, capsys): def test_cli_validate_json(capsys): - assert cli(["validate", str(ITAITUBA), "--json"]) == 0 + assert cli(["validate", str(QUICKSTART), "--json"]) == 0 out = capsys.readouterr().out data = json.loads(out) assert data["status"] == "ok" - assert data["grid"] == "itaituba/5km" - assert "elevation" in data["derives"] + assert data["grid"] == "demo/300m" + assert "elev_mean" in data["derives"] assert "elevation" in data["sources"] @@ -289,47 +289,47 @@ def test_cli_validate_json_error(tmp_path, capsys): def test_cli_run_dry_run(capsys): - assert cli(["run", str(ITAITUBA), "--dry-run"]) == 0 + assert cli(["run", str(QUICKSTART), "--dry-run"]) == 0 out = capsys.readouterr().out assert "[dry-run]" in out assert "Plan is valid" in out def test_cli_run_dry_run_json(capsys): - assert cli(["run", str(ITAITUBA), "--dry-run", "--json"]) == 0 + assert cli(["run", str(QUICKSTART), "--dry-run", "--json"]) == 0 out = capsys.readouterr().out data = json.loads(out) assert data["status"] == "ok" assert data["dry_run"] is True - assert data["grid"] == "itaituba/5km" + assert data["grid"] == "demo/300m" assert any(s["id"] == "elevation" for s in data["sources"]) - assert any(d["target"] == "elevation" for d in data["derives"]) + assert any(d["target"] == "elev_mean" for d in data["derives"]) def test_cli_run_json(tmp_path, capsys): ws = tmp_path / "ws_json" - assert cli(["run", str(ITAITUBA), "--workspace", str(ws), "--json"]) == 0 + assert cli(["run", str(QUICKSTART), "--workspace", str(ws), "--json"]) == 0 out = capsys.readouterr().out data = json.loads(out) assert data["status"] == "ok" assert data["workspace"] == str(ws) - assert data["grid_id"] == "itaituba/5km" - assert len(data["derived"]) == 5 + assert data["grid_id"] == "demo/300m" + assert len(data["derived"]) == 3 assert (ws / "run.json").exists() def test_cli_export_json(tmp_path, capsys): ws = tmp_path / "ws_export" - assert cli(["run", str(ITAITUBA), "--workspace", str(ws)]) == 0 + assert cli(["run", str(QUICKSTART), "--workspace", str(ws)]) == 0 capsys.readouterr() out_tif = tmp_path / "exported.tif" - assert cli(["export", str(ITAITUBA), "--workspace", str(ws), "--output", str(out_tif), "--json"]) == 0 + assert cli(["export", str(QUICKSTART), "--workspace", str(ws), "--output", str(out_tif), "--json"]) == 0 out = capsys.readouterr().out data = json.loads(out) assert data["status"] == "ok" assert data["workspace"] == str(ws) - assert data["grid_id"] == "itaituba/5km" - assert "elevation" in data["variables"] + assert data["grid_id"] == "demo/300m" + assert "elev_mean" in data["variables"] assert data["output"] == str(out_tif) assert out_tif.exists() diff --git a/tests/test_terrame_parity.py b/tests/test_terrame_parity.py deleted file mode 100644 index c87f21f..0000000 --- a/tests/test_terrame_parity.py +++ /dev/null @@ -1,178 +0,0 @@ -""" -Parity with TerraME's *Fill*: DisSCube vs TerraME's own output, cell by cell. - -The three Fill examples of TerraME's ``gis`` package (Itaituba, Emas, -Amazônia) ship with the cellular spaces TerraME produced; they are bundled in -``examples/data/terrame/``. Each test derives one attribute with DisSCube on -the same grid and compares it with TerraME's value in every cell. - -Passing tests pin the parity reached today; strict ``xfail`` tests record the -known gaps (see docs/terrame_fill_correspondence.md) and will start failing — -as a reminder to promote them — once a gap is closed. -""" - -from dataclasses import dataclass -from pathlib import Path - -import geopandas as gpd -import numpy as np -import pytest - -from disscube import CubeClient, GridSpec, SpatialDerivation, SpatialSource, Variable - -DATA = Path(__file__).resolve().parents[1] / "examples" / "data" / "terrame" - - -@dataclass -class Case: - cube: CubeClient - ref: gpd.GeoDataFrame - rows: np.ndarray - cols: np.ndarray - crs: str - folder: Path - - def source(self, sid: str, fmt: str, name: str) -> None: - url = f"zip://{self.folder / name}.zip" if fmt == "vector" else str(self.folder / name) - self.cube.register_spatial_source(SpatialSource(id=sid, name=sid, format=fmt, asset_url=url, crs=self.crs)) - - def derive(self, sid: str, name: str, operator: str, class_code: int | None = None, purity: bool = False): - self.cube.derive(SpatialDerivation( - source_id=sid, grid_id="g", role="driver", - variables=[Variable(name=name, operator=operator, class_code=class_code)], - )) - da = self.cube.load(name, grid_id="g") - if purity: - da = da * da.coords["coverage_purity"] - return da.to_numpy()[self.rows, self.cols] - - -def _case(tmp_path_factory, folder: str, epsg: int, resolution: float) -> Case: - path = DATA / folder - # The reference cells carry a wrong .prj; their real CRS is the inputs' one. - ref = gpd.read_file(f"zip://{path / folder}.zip").set_crs(epsg, allow_override=True) - xmin, _, _, ymax = ref.total_bounds - tmp = tmp_path_factory.mktemp(folder) - cube = CubeClient(str(tmp / "catalog.db"), str(tmp / "store")) - cube.register_grid(GridSpec(id="g", type="local", crs=f"EPSG:{epsg}", resolution=resolution, - bbox=ref.total_bounds.tolist())) - c = ref.geometry.centroid # TerraME rows grow northwards; DisSCube is north-up - cols = np.floor((c.x.to_numpy() - xmin) / resolution).astype(int) - rows = np.floor((ymax - c.y.to_numpy()) / resolution).astype(int) - return Case(cube, ref, rows, cols, f"EPSG:{epsg}", path) - - -def _share(ours, theirs, tol): - return float(np.mean(np.abs(ours - theirs) <= tol)) - - -# ═══ Itaituba — 620 cells, 5 km ═══════════════════════════════════════════════ - -@pytest.fixture(scope="module") -def itaituba(tmp_path_factory): - case = _case(tmp_path_factory, "itaituba", 29191, 5_000.0) - case.source("elev", "raster", "itaituba-elevation.tif") - case.source("defor", "raster", "itaituba-deforestation.tif") - case.source("roads", "vector", "itaituba-roads") - case.source("localities", "vector", "itaituba-localities") - case.source("census", "vector", "itaituba-census") - return case - - -def test_itaituba_average_elevation(itaituba): - ours = itaituba.derive("elev", "elevation", "mean") - err = np.abs(ours - itaituba.ref["elevation"].to_numpy()) - assert err.mean() < 1.0 # m (0.89 today) - assert err.max() < 11.0 # m (10.13 today) - - -@pytest.mark.parametrize("k", [7, 87, 167, 255]) -def test_itaituba_coverage(itaituba, k): - ours = itaituba.derive("defor", f"defor_{k}", "percentage", class_code=k, purity=True) * 100 - err = np.abs(ours - itaituba.ref[f"defor_{k}"].to_numpy()) - assert err.max() < 0.7 # pp, in every cell (0.64 today) - - -@pytest.mark.xfail(strict=True, reason="min_distance is a raster approximation; TerraME uses exact polygon distance") -@pytest.mark.parametrize(("sid", "attr"), [("roads", "distroad"), ("localities", "distlocal")]) -def test_itaituba_distance(itaituba, sid, attr): - ours = itaituba.derive(sid, attr, "min_distance") - assert _share(ours, itaituba.ref[attr].to_numpy(), 100) >= 0.75 - - -@pytest.mark.xfail(strict=True, reason="area-weighted sum over polygons (TerraME sum, area=true) not implemented") -def test_itaituba_area_weighted_sum(itaituba): - ours = itaituba.derive("census", "population", "sum") - assert _share(ours, itaituba.ref["population"].to_numpy(), 0.5) >= 0.95 - - -# ═══ Emas — 5 514 cells, 500 m ════════════════════════════════════════════════ - -@pytest.fixture(scope="module") -def emas(tmp_path_factory): - case = _case(tmp_path_factory, "emas", 29192, 500.0) - case.source("firebreak", "vector", "emas-firebreak") - case.source("river", "vector", "emas-river") - case.source("cover", "raster", "emas-accumulation.tif") - return case - - -@pytest.mark.parametrize("attr", ["firebreak", "river"]) -def test_emas_presence(emas, attr): - ours = emas.derive(attr, attr, "presence") - assert _share(ours, emas.ref[attr].to_numpy(), 0) >= 0.98 # 98.4 % / 99.7 % today - - -@pytest.mark.xfail(strict=True, reason="presence rasterizes line centres; TerraME marks every cell a line touches") -@pytest.mark.parametrize("attr", ["firebreak", "river"]) -def test_emas_presence_exact(emas, attr): - ours = emas.derive(attr, f"{attr}_exact", "presence") - assert _share(ours, emas.ref[attr].to_numpy(), 0) == 1.0 - - -@pytest.mark.parametrize(("attr", "op"), [("maxcover", "max"), ("mincover", "min")]) -def test_emas_min_max(emas, attr, op): - ours = emas.derive("cover", attr, op) - assert _share(ours, emas.ref[attr].to_numpy(), 0) >= 0.985 # 98.7 % / 99.0 % today - - -@pytest.mark.xfail(strict=True, reason="min/max count pixels partly inside a cell; TerraME uses pixel centres") -@pytest.mark.parametrize(("attr", "op"), [("maxcover", "max"), ("mincover", "min")]) -def test_emas_min_max_exact(emas, attr, op): - ours = emas.derive("cover", f"{attr}_exact", op) - assert _share(ours, emas.ref[attr].to_numpy(), 0) == 1.0 - - -# ═══ Amazônia — 2 229 cells, 50 km ════════════════════════════════════════════ - -@pytest.fixture(scope="module") -def amazonia(tmp_path_factory): - case = _case(tmp_path_factory, "amazonia", 29191, 50_000.0) - case.source("prodes", "raster", "amazonia-prodes.tif") - case.source("roads", "vector", "amazonia-roads") - case.source("ports", "vector", "amazonia-ports") - case.source("indigenous", "vector", "amazonia-indigenous") - return case - - -@pytest.mark.parametrize("k", [10, 208]) -def test_amazonia_coverage(amazonia, k): - ours = amazonia.derive("prodes", f"prodes_{k}", "percentage", class_code=k, purity=True) * 100 - theirs = amazonia.ref[f"prodes_{k}"].to_numpy() - no_data = np.isnan(ours) - # Cells without any PRODES pixel: DisSCube reports NaN, TerraME 0. - assert (theirs[no_data] == 0).all() - np.testing.assert_allclose(ours[~no_data], theirs[~no_data], atol=1e-6) - - -@pytest.mark.xfail(strict=True, reason="min_distance is a raster approximation; TerraME uses exact polygon distance") -@pytest.mark.parametrize(("sid", "attr"), [("roads", "distroads"), ("ports", "distports")]) -def test_amazonia_distance(amazonia, sid, attr): - ours = amazonia.derive(sid, attr, "min_distance") - assert _share(ours, amazonia.ref[attr].to_numpy(), 100) >= 0.75 - - -def test_amazonia_area_fraction(amazonia): - """TerraME's ``area`` (intersection area / cell area) is DisSCube's ``area``.""" - ours = amazonia.derive("indigenous", "protected", "area") - assert _share(ours, amazonia.ref["protected"].to_numpy(), 0.01) >= 0.99