diff --git a/CHANGELOG.md b/CHANGELOG.md index 298b5af..07a672b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,61 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 --- +## [0.6.0] - 2026-10-14 + +### Added +- **`disscube plot`**: a PNG of the bands of an exported GeoTIFF, to look at a result without a script per case. + `disscube plot cellspace.tif [-o look.png] [--bands NAME ...] [--year 2024] [--cols N]`. One panel per band; the scale + follows the values (shares in [0, 1] as 0 to 100 %, a few integer codes as classes, otherwise continuous between the 2nd + and 98th percentiles); cells without data are white. It is a quick look and not a publication figure. Needs matplotlib + (`pip install 'disscube[plot]'`). `--shared-colorbar` draws one bar for all panels when every panel is a share (0 to 100 %); + `--crop` trims the empty margin (bounding box of the cells with data in any drawn band). +- **PRODES: the local code of each category is declared in the pipeline file.** `[source.codes]` (or `codes=` in + `register_prodes_source`) gives the code of `forest`, `deforested`, `non_forest` and `water`; categories with the same code + are merged. The default is the map of before (1 forest, 2 deforested up to the year, 3 non-forest vegetation and water + together); `water = 4` keeps the hydrography apart, so a map of the classes shows rivers apart from non-forest vegetation. + The legend is still read by label (`dAAAA`, `rAAAA`, hidrografia, não floresta, floresta; clouds are no data), so the + categories do not depend on the numbers of an edition. See `docs/guides/landcover.md`. +- **`majority` and `minority` accept `params = {groups = {label = [codes]}}`**: the map is read at a coarser level and + the result is the label of the group that covers most (least) of the cell. One derivation gives the MapBiomas level-1 + class (`"1" = [3, 4, 5, 6, 49]`, …) from the leaf-class map. Codes in no group, and nodata, are ignored; a code in two + groups is an error. In a pipeline file the table is `[derive.params.groups]`. +- **`percentage` accepts a list of classes**: `class_code = [3, 4, 5, 6, 49]` is the share of the cell in any of + them, for classes that are a group of codes (a MapBiomas level-1 class is a group of leaf codes). Raster sources + only; a vector source still takes one class. +- **Sources and grids larger than memory.** The MapBiomas and PRODES adapters read the window and write the + registered GeoTIFF a block of rows at a time, in the file's own `uint8` (new `copy_window`, `register_file`), + instead of one float32 array of the whole window. `derive` aligns a raster source a block of grid cells (rows × columns) at + a time when it has more than `block_pixels` pixels (default 50 million, `DISSCUBE_BLOCK_PIXELS` or + `derive(..., block_pixels=)`), for the operators that work cell by cell (`Operator.block_safe`: mean, sum, std, + median, min, max, majority, minority, percentage). The aligner clips the source to the grid before + reprojecting. No new TOML syntax: the PRODES of the whole Legal Amazon (86,679 × 111,508 pixels) on a 5.28 km + grid, or MapBiomas over Brazil at 30 m, is one pipeline file. Registered MapBiomas files are byte for byte the + ones of 0.5.0. The block size is chosen from the measured footprint of each block in the source (a wide strip of a projected grid covers far more longitude/latitude than its height suggests, so rows alone do not bound the memory). Class shares, majority and median of a block equal those of the whole grid when GDAL's exact + transformer is used; with its default approximate one a cell on an edge can pick the neighbouring source + pixel (1 cell in 3,422 in the test). `mean` and `sum` agree to about 1e-12 with the exact transformer. The exact transformer can be selected only with rasterio 1.5 or later (Python 3.12 or later); with rasterio 1.4, on Python 3.11, a few cells of the categorical operators can differ between blocks and whole grid, and the equality test is skipped there. The BDC adapter still reads whole windows. +- **`run.json` records the DisSCube version** (`"disscube": "0.6.0"`, `null` when running from a source tree that is not installed), so a run can be tied to a release without noting the commit by hand. + +### Fixed +- **Invalid polygons in a vector source no longer stop the pipeline.** The reprojection can turn a valid polygon into a + self-intersecting one, and `clip` then raised `TopologyException` (found with the INCRA settlements layer of the + open-lucc-workflow, one bad feature among 10,000). The aligner now repairs the invalid features (`make_valid`) after the + reprojection and before the clip, and logs how many; valid features are not touched. +- **PRODES: the yearly classes d2010 to d2025 were dropped.** In the edition's legend every year from 2010 has a + yearly class (dYYYY, codes 10 to 25) and a residual class (rYYYY, codes 50 to 65); the lookup was built from a dictionary + keyed by year, so the residual class replaced the yearly one and the pixels of d2010 to d2025 came out as no data. A map of + forest or of accumulated deforestation lacked them (about a quarter of the deforestation of Acre, whose accumulated + deforestation rose by 1 thousand km² between 2008 and 2024 in place of 12 thousand). Both kinds of class are now in the + lookup. Maps made with an earlier version should be made again. +- **`Error in sys.excepthook:` at the end of `disscube run`.** The aligner opened the source raster and never + closed it; the open handle was finalised at interpreter exit. It is now closed once the aligned arrays are in memory. +- **SQLite connections of the catalog are closed.** `with sqlite3.connect(...)` only ends the transaction; + every catalog call left a connection open until interpreter exit (`ResourceWarning: unclosed database`). + +--- + +--- + ## [0.5.0] - 2026-10-05 ### Added diff --git a/CITATION.cff b/CITATION.cff index c53ee3b..c247e03 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -6,8 +6,8 @@ authors: orcid: "https://orcid.org/0000-0002-0232-4549" title: "DisSCube: Declarative spatial data cubes" abstract: "An open-source, declarative engine for constructing cellular spatial data cubes for dynamic modeling, spatial simulation, and environmental analysis." -version: 0.5.0 -date-released: "2026-10-05" +version: 0.6.0 +date-released: "2026-10-14" url: "https://github.com/DisSModel/disscube" repository-code: "https://github.com/DisSModel/disscube" license: MIT diff --git a/README.md b/README.md index e17d51e..71119c9 100644 --- a/README.md +++ b/README.md @@ -248,7 +248,7 @@ disscube/ ├── pipeline/ Pipeline execution & planning (schema, runner) + internal stages ├── catalog/ CatalogStore (Protocol) + SQLite and JSON implementations ├── storage.py AssetStore (fsspec — local and S3) -├── cli.py `disscube validate` / `disscube run` / `disscube export` +├── cli.py `disscube validate` / `disscube run` / `disscube export` / `disscube plot` ├── sources/ Adapters that bring external data in as SpatialSources, │ │ each with a checksum and a provenance.json sidecar │ ├── _raster.py Window2D, windowed reads, composites, mosaics, register_raster @@ -286,7 +286,7 @@ The operator is registered automatically and accepted by `Derivation` / `Spatial The limitations below are scope decisions for the current version, not bugs. They are documented so that users and reviewers understand what is implemented versus what is planned. **In-memory, single-tile processing** -Each call to `derive()` loads a tile's full data into memory. There is no lazy (Dask) or distributed processing. For continental-scale grids (e.g. `BR/1km`), use the tile loop — each tile is processed and saved independently. +Each call to `derive()` loads a tile's full data into memory, except raster sources of MapBiomas and PRODES, which are read, written and aligned by blocks of rows and columns (`DISSCUBE_BLOCK_PIXELS`, default 50 million pixels) and so can be larger than memory. Other sources are read in memory; for the Brazil Data Cube, the item windows of a BDC composite are all held before the reduction. Each `[[derive]]` of a block-read source re-reads and re-aligns it, one pass per variable and year. There is no lazy (Dask) or distributed processing. For continental-scale grids with other sources, use the tile loop; each tile is processed and saved independently. **Vector aggregation by rasterization (not area-weighted)** Operators over vector sources (`majority`, `percentage`, `attribute`, `presence`, `minority`) convert geometries to raster before aggregating pixels. For the share of each cell covered by polygons, use `area`, which intersects the polygons with the cells exactly. @@ -312,7 +312,7 @@ If you use DisSCube in your research, dynamic modeling, or spatial data pipeline author = {Costa, S{\'e}rgio Souza}, title = {{DisSCube: Declarative spatial data cubes}}, year = {2026}, - version = {0.5.0}, + version = {0.6.0}, url = {https://github.com/DisSModel/disscube} } ``` diff --git a/disscube/catalog/sqlite_store.py b/disscube/catalog/sqlite_store.py index bb32904..c36e268 100644 --- a/disscube/catalog/sqlite_store.py +++ b/disscube/catalog/sqlite_store.py @@ -1,3 +1,4 @@ +import contextlib import json import sqlite3 from pathlib import Path @@ -10,10 +11,17 @@ def __init__(self, path: str | Path): self.path = Path(path) self._init_db() + @contextlib.contextmanager def _get_connection(self): + # `with sqlite3.connect(...)` only ends the transaction (commit/rollback); it does not + # close the connection. Close it explicitly so none is left open until interpreter exit. conn = sqlite3.connect(self.path) conn.row_factory = sqlite3.Row - return conn + try: + with conn: + yield conn + finally: + conn.close() def _init_db(self): with self._get_connection() as conn: diff --git a/disscube/cli.py b/disscube/cli.py index 1a8fad7..deebd74 100644 --- a/disscube/cli.py +++ b/disscube/cli.py @@ -5,6 +5,7 @@ disscube fetch pipeline.toml disscube run pipeline.toml [--workspace DIR] [--output OUT.tif] [--dry-run] [--json] [-v] disscube export pipeline.toml --output OUT.tif [--workspace DIR] [--variables ...] [--json] [-v] + disscube plot OUT.tif [-o LOOK.png] [--bands NAME ...] [--year YYYY] [--cols N] [--shared-colorbar] [--crop] [--dpi N] """ from __future__ import annotations @@ -42,7 +43,26 @@ def main(argv: list[str] | None = None) -> int: p_exp.add_argument("--json", action="store_true", help="output export result as JSON") p_exp.add_argument("-v", "--verbose", action="store_true", help="log each step") + p_plot = sub.add_parser("plot", help="draw the bands of an exported GeoTIFF into a PNG, to look at the result") + p_plot.add_argument("file", help="GeoTIFF written by run --output or export") + p_plot.add_argument("--output", "-o", help="PNG (default: the GeoTIFF's name with .png)") + p_plot.add_argument("--bands", nargs="*", help="band or variable names (default: all)") + p_plot.add_argument("--year", type=int, help="only the bands of this year") + p_plot.add_argument("--cols", type=int, help="panels per row") + p_plot.add_argument("--shared-colorbar", action="store_true", help="one colour bar for all panels (when all are shares)") + p_plot.add_argument("--crop", action="store_true", help="trim the empty margin (bounding box of the cells with data)") + p_plot.add_argument("--dpi", type=int, default=150) + args = parser.parse_args(argv) + if args.command == "plot": + from disscube.plot import plot_geotiff + + try: + print(f"figure : {plot_geotiff(args.file, args.output, bands=args.bands, year=args.year, cols=args.cols, dpi=args.dpi, shared_colorbar=args.shared_colorbar, crop=args.crop)}") + except (ImportError, ValueError, OSError) as exc: + print(f"error: {exc}", file=sys.stderr) + return 2 + return 0 from disscube.pipeline import PipelineError, plan, run logging.basicConfig(level=logging.INFO if getattr(args, "verbose", False) else logging.WARNING, diff --git a/disscube/client.py b/disscube/client.py index e35a084..c50ef56 100644 --- a/disscube/client.py +++ b/disscube/client.py @@ -13,7 +13,7 @@ from disscube.catalog.sqlite_store import SqliteCatalogStore from disscube.models import DerivedVariable, GridSpec, SpatialDerivation, SpatialRelation, SpatialSource from disscube.pipeline import PipelineContext -from disscube.pipeline.aggregator import Aggregator +from disscube.pipeline.aggregator import Aggregator, write_georeference from disscube.pipeline.aligner import GridAligner from disscube.pipeline.normalizer import Normalizer from disscube.pipeline.writer import VariableWriter @@ -28,6 +28,82 @@ log = logging.getLogger("disscube.client.cube_client") +#: Source pixels aligned at once before a raster source is processed a block of grid cells at a time. +DEFAULT_BLOCK_PIXELS = 50_000_000 + + +def _source_pixels(bbox: list[float], grid_crs: str, ds: Any) -> float: + """Source pixels inside ``bbox`` (grid CRS): its footprint in the source CRS, clipped to the source.""" + from rasterio.warp import transform_bounds + + left, bottom, right, top = transform_bounds(grid_crs, ds.crs, *bbox, densify_pts=21) + b = ds.bounds + w = min(right, b.right) - max(left, b.left) + h = min(top, b.top) - max(bottom, b.bottom) + if w <= 0 or h <= 0: + return 0.0 + return (w / abs(ds.res[0])) * (h / abs(ds.res[1])) + + +def _blocks(source: SpatialSource, grid: GridSpec, derivation: SpatialDerivation, + block_pixels: int | None) -> list[list[GridSpec]] | None: + """ + Blocks of grid cells (a list of block rows, each a list of blocks), each covering at most about + ``block_pixels`` source pixels, or None to derive the grid whole. + + The footprint of a block in the source is measured, not assumed: a strip of grid rows that is + wide in a projected grid covers a far larger area of a longitude/latitude source than its + height suggests, so blocks are halved along their larger side until the budget is met. + None when the source is a vector, is small enough, or a variable uses an operator (or + ``fill = "nearest"``) whose cells depend on their neighbours. + """ + import rasterio + + from disscube.operators.base import OPERATOR_REGISTRY + + if source.format != "raster": + return None + for v in derivation.variables: + op = OPERATOR_REGISTRY.get(v.operator) + if op is None or not op.block_safe or getattr(v, "fill", None) == "nearest": + return None + budget = block_pixels or int(os.environ.get("DISSCUBE_BLOCK_PIXELS", DEFAULT_BLOCK_PIXELS)) + minx, _miny, _maxx, maxy = grid.bbox + res = grid.resolution + + def bbox_of(r0: int, r1: int, c0: int, c1: int) -> list[float]: + return [minx + c0 * res, maxy - r1 * res, minx + c1 * res, maxy - r0 * res] + + try: + with rasterio.open(source.asset_url) as ds: + def worst(nr: int, nc: int) -> float: + return max( + _source_pixels(bbox_of(r0, min(r0 + nr, grid.rows), c0, min(c0 + nc, grid.cols)), + grid.crs, ds) + for r0 in range(0, grid.rows, nr) for c0 in range(0, grid.cols, nc) + ) + + nr, nc = grid.rows, grid.cols + while worst(nr, nc) > budget and (nr > 1 or nc > 1): + if nr >= nc: + nr = max(1, (nr + 1) // 2) + else: + nc = max(1, (nc + 1) // 2) + except Exception: # noqa: BLE001 — not a plain file (NetCDF subdataset, URL…): derive whole + return None + if nr >= grid.rows and nc >= grid.cols: + return None + out = [] + for r0 in range(0, grid.rows, nr): + r1 = min(r0 + nr, grid.rows) + out.append([ + GridSpec(id=grid.id, type=grid.type, crs=grid.crs, resolution=res, + bbox=bbox_of(r0, r1, c0, min(c0 + nc, grid.cols)), description=grid.description) + for c0 in range(0, grid.cols, nc) + ]) + return out + + # Attributes that describe one stored slice: per-variable only when it has a single slice. _SLICE_ATTRS = ("spec_hash", "content_hash", "source_id", "source_checksum") @@ -68,7 +144,16 @@ def get_relations(self, grid_id: str) -> list[SpatialRelation]: def search(self, grid: str | None = None, role: str | None = None) -> list[DerivedVariable]: return self.catalog.search_derived_variables(grid_id=grid, role=role) - def derive(self, derivation: SpatialDerivation, tile_id: str | None = None) -> list[DerivedVariable]: + def derive(self, derivation: SpatialDerivation, tile_id: str | None = None, + block_pixels: int | None = None) -> list[DerivedVariable]: + """ + Derive the variables of ``derivation`` on its grid. + + ``block_pixels`` is the number of source pixels that may be aligned at once. A raster source + larger than that is processed a block of grid cells at a time (default + :data:`DEFAULT_BLOCK_PIXELS`, or the ``DISSCUBE_BLOCK_PIXELS`` environment variable), which + bounds the memory without changing the result; it needs operators that work cell by cell. + """ derivation = derivation.model_copy(deep=True) if not derivation.relations: derivation.relations = self.get_relations(derivation.grid_id) @@ -113,17 +198,35 @@ def derive(self, derivation: SpatialDerivation, tile_id: str | None = None) -> l description=f"Temporary tile grid for {tile_id}" ) - ctx = PipelineContext(source=source, grid=grid, derivation=derivation, tile_id=tile_id) - - pipeline = [ - Normalizer(), - GridAligner(), - Aggregator(), - VariableWriter(self.store, self.catalog) - ] - - for stage in pipeline: - ctx = stage.execute(ctx) + blocks = _blocks(source, grid, derivation, block_pixels) + if blocks is None: + ctx = PipelineContext(source=source, grid=grid, derivation=derivation, tile_id=tile_id) + for stage in (Normalizer(), GridAligner(), Aggregator()): + ctx = stage.execute(ctx) + else: + # A source too large to align in one piece: derive the grid a block of cells at a time + # and join the blocks. Every cell depends only on the source pixels inside it (see + # Operator.block_safe), so the result is the same as for the whole grid. + total = sum(len(row) for row in blocks) + n = 0 + band = [] + for row in blocks: + parts = [] + for sub in row: + n += 1 + log.info("block %d/%d of grid %s (%d rows x %d cols)", n, total, grid.id, sub.rows, sub.cols) + sctx = PipelineContext(source=source, grid=sub, derivation=derivation, tile_id=tile_id) + for stage in (Normalizer(), GridAligner(), Aggregator()): + sctx = stage.execute(sctx) + parts.append(sctx.data) + band.append(xr.concat(parts, dim="x", coords="minimal", compat="override") + if len(parts) > 1 else parts[0]) + joined = (xr.concat(band, dim="y", coords="minimal", compat="override") + if len(band) > 1 else band[0]) + joined = joined.assign_coords(y=grid.ys, x=grid.xs) + write_georeference(joined, grid) + ctx = PipelineContext(source=source, grid=grid, derivation=derivation, tile_id=tile_id, data=joined) + ctx = VariableWriter(self.store, self.catalog).execute(ctx) return [ d for d in self.catalog.search_derived_variables(tile_id=tile_id) @@ -135,6 +238,7 @@ def derive_declarative( derivation: Derivation, grid_id: str, tile_id: str | None = None, + block_pixels: int | None = None, ) -> list[DerivedVariable]: """ Thin convenience wrapper: build a ``SpatialDerivation`` from a @@ -151,7 +255,7 @@ def derive_declarative( tile_id : str | None Optional tile sub-identifier, forwarded to ``derive()``. """ - return self.derive(derivation.to_spatial_derivation(grid_id), tile_id=tile_id) + return self.derive(derivation.to_spatial_derivation(grid_id), tile_id=tile_id, block_pixels=block_pixels) def purge_stale(self) -> int: """ diff --git a/disscube/models/derivation.py b/disscube/models/derivation.py index b469925..52d4dc7 100644 --- a/disscube/models/derivation.py +++ b/disscube/models/derivation.py @@ -76,7 +76,7 @@ class Derivation(BaseModel): target: str source_id: str operator: str - class_code: int | None = None + class_code: int | list[int] | None = None role: str = "driver" valid_from: str | None = None valid_until: str | None = None diff --git a/disscube/models/variable.py b/disscube/models/variable.py index 55624ab..b99d5d7 100644 --- a/disscube/models/variable.py +++ b/disscube/models/variable.py @@ -10,7 +10,7 @@ class Variable(BaseModel): name: str operator: str - class_code: int | None = None + class_code: int | list[int] | None = None params: dict[str, Any] = {} """Operator options (e.g. ``{"crs": "EPSG:5880"}`` for ``distance``).""" fill: Literal["nearest"] | None = None diff --git a/disscube/operators/base.py b/disscube/operators/base.py index 1747d22..a32e4ac 100644 --- a/disscube/operators/base.py +++ b/disscube/operators/base.py @@ -64,6 +64,12 @@ class body is executed (via ``__init_subclass__``), so there is no # composition rather than a pre-collapsed mode. needs_fine_alignment: ClassVar[bool] = False + # True when each output cell depends only on the source pixels inside that cell (the zonal + # statistics and class shares over a raster). The pipeline may then derive the grid a block of + # rows at a time, so a source larger than memory can be processed; the result is the same. + # Operators that look across cells (distances, network costs) or at vector features stay False. + block_safe: ClassVar[bool] = False + # When False, GridAligner does not clip a vector source to the grid bbox # for this operator — needed by operators that measure to features lying # outside the grid (``distance``). A derivation is clipped only if every diff --git a/disscube/operators/zonal.py b/disscube/operators/zonal.py index ed06fef..ad8f9ba 100644 --- a/disscube/operators/zonal.py +++ b/disscube/operators/zonal.py @@ -73,6 +73,34 @@ def _fine_array(data: xr.DataArray) -> tuple[np.ndarray, float | None]: return arr, nodata +def _regroup(fine: np.ndarray, groups) -> np.ndarray: + """ + Replace the class codes of ``fine`` by the label of the group that lists them. + + ``groups`` is ``{label: [codes]}`` (labels may be strings, as TOML keys are). Pixels whose code is in no + group, and nodata, become NaN, which the reduction ignores. This is how a map of leaf classes is read + at a coarser level (MapBiomas level 1: ``{1 = [3, 4, 5, 6, 49], 3 = [15, 39, ...]}``). + """ + if not isinstance(groups, dict) or not groups: + raise ValueError("'groups' must be a table {label = [codes]}") + out = np.full(fine.shape, np.nan, dtype=np.float64) + seen: set[int] = set() + for label, codes in groups.items(): + try: + value = float(int(label)) + except (TypeError, ValueError): + raise ValueError(f"group label {label!r} must be an integer") from None + codes = [codes] if isinstance(codes, int) else list(codes) + if seen & set(codes): + raise ValueError(f"code(s) {sorted(seen & set(codes))} listed in more than one group") + seen |= set(codes) + out[np.isin(fine, np.asarray(codes, dtype=np.float64))] = value + return out + + +_GROUPS_DOC = "{label = [codes]}: read the map at a coarser level; the result is the label of the group (codes in no group are ignored)" + + def _windows(n_fine: int, n_target: int) -> list[tuple[int, int]]: """ Build per-target-cell (start, stop) index ranges along one axis. @@ -95,7 +123,7 @@ def _categorical_reduce( nodata: float | None, grid: GridSpec, mode: str, - class_code: int | None, + class_code: int | list[int] | None, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """ Reduce a fine categorical array into the target grid. @@ -110,8 +138,9 @@ def _categorical_reduce( Target grid (provides ``rows`` / ``cols``). mode : {"percentage", "majority", "minority"} Aggregation to compute for the primary output. - class_code : int | None - Target class for ``percentage``. + class_code : int | list[int] | None + Target class for ``percentage``; with a list, the share of the cell in any of them + (e.g. the leaf classes of one MapBiomas level-1 class). Returns ------- @@ -123,6 +152,7 @@ def _categorical_reduce( dominance_purity : np.ndarray count(dominant_class) / valid_pixels per cell (categorical metric). """ + wanted = np.atleast_1d(class_code if class_code is not None else []).astype(np.int64) fr, fc = fine.shape tr, tc = grid.rows, grid.cols @@ -161,7 +191,7 @@ def _categorical_reduce( dominance[ti, tj] = dom_count / n_valid if mode == "percentage": - hit = int(counts[classes == class_code].sum()) if class_code in classes else 0 + hit = int(counts[np.isin(classes, wanted)].sum()) value[ti, tj] = hit / n_valid elif mode == "majority": value[ti, tj] = dom_class @@ -269,6 +299,7 @@ def _continuous_reduce( # --------------------------------------------------------------------------- class MeanOperator(Operator): + block_safe = True name = "mean" _resampling = Resampling.average @@ -297,6 +328,7 @@ class SumOperator(Operator): clipped to the grid, so a polygon crossing the border keeps its full denominator and only the part inside the grid is distributed. """ + block_safe = True name = "sum" _resampling = Resampling.sum @@ -383,6 +415,7 @@ def _cell_boxes(grid: GridSpec) -> np.ndarray: class StdOperator(Operator): + block_safe = True name = "std" # Standard deviation has no rasterio resampling equivalent; it requires a # real per-cell window pass, so this operator uses fine alignment. @@ -412,6 +445,7 @@ class MedianOperator(Operator): aligned to a fine grid first: the median of per-cell means would not be the median of the pixels. """ + block_safe = True name = "median" _resampling = Resampling.nearest needs_fine_alignment = True @@ -427,6 +461,7 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: class MinOperator(Operator): + block_safe = True name = "min" _resampling = Resampling.min @@ -437,6 +472,7 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: class MaxOperator(Operator): + block_safe = True name = "max" _resampling = Resampling.max @@ -451,14 +487,20 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: # --------------------------------------------------------------------------- class MajorityOperator(Operator): + block_safe = True name = "majority" _resampling = Resampling.nearest needs_fine_alignment = True - params: ClassVar[dict[str, str]] = {"subcells": "at most this many fine pixels per cell along each axis (memory bound)"} + params: ClassVar[dict[str, str]] = { + "subcells": "at most this many fine pixels per cell along each axis (memory bound)", + "groups": _GROUPS_DOC, + } def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: if isinstance(data, xr.DataArray): fine, nodata = _fine_array(data) + if var.params.get("groups"): + fine, nodata = _regroup(fine, var.params["groups"]), None value, cov, dom = _categorical_reduce(fine, nodata, grid, "majority", None) return _attach_purity(value, cov, dom, grid) if isinstance(data, gpd.GeoDataFrame): @@ -469,14 +511,20 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: class MinorityOperator(Operator): + block_safe = True name = "minority" _resampling = Resampling.nearest needs_fine_alignment = True - params: ClassVar[dict[str, str]] = {"subcells": "at most this many fine pixels per cell along each axis (memory bound)"} + params: ClassVar[dict[str, str]] = { + "subcells": "at most this many fine pixels per cell along each axis (memory bound)", + "groups": _GROUPS_DOC, + } def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: if isinstance(data, xr.DataArray): fine, nodata = _fine_array(data) + if var.params.get("groups"): + fine, nodata = _regroup(fine, var.params["groups"]), None value, cov, dom = _categorical_reduce(fine, nodata, grid, "minority", None) return _attach_purity(value, cov, dom, grid) if isinstance(data, gpd.GeoDataFrame): @@ -487,6 +535,7 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: class PercentageOperator(Operator): + block_safe = True name = "percentage" requires_class_code = True _resampling = Resampling.nearest @@ -505,6 +554,8 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: if isinstance(data, gpd.GeoDataFrame): if var.class_code is None: raise ValueError("'percentage' requires class_code to be set.") + if isinstance(var.class_code, list): + raise TypeError("'percentage' over a vector source takes a single class_code, not a list.") shapes = [(g, var.class_code) for g in data.geometry if g is not None] return _wrap(_rasterize(shapes, grid), grid) raise TypeError(f"'percentage' got unexpected type {type(data).__name__}") diff --git a/disscube/pipeline/aggregator.py b/disscube/pipeline/aggregator.py index e5d893b..edde095 100644 --- a/disscube/pipeline/aggregator.py +++ b/disscube/pipeline/aggregator.py @@ -101,6 +101,12 @@ def _fill_nearest(result: xr.DataArray) -> xr.DataArray: return result.copy(data=values[rows, cols]) +def write_georeference(ds: xr.Dataset, grid) -> None: + """Write the CRS and the transform of ``grid`` on ``ds`` (in place).""" + ds.rio.write_crs(_crs_to_named_wkt(grid.crs), inplace=True) + ds.rio.write_transform(grid.transform, inplace=True) + + class Aggregator(PipelineStage): def execute(self, ctx: PipelineContext) -> PipelineContext: grid = ctx.grid @@ -144,9 +150,7 @@ def execute(self, ctx: PipelineContext) -> PipelineContext: # rioxarray.write_crs() sets grid_mapping="spatial_ref" on every existing # data variable — calling it before the loop would miss them all, causing # QGIS and other CF-compliant readers to not detect the projection. - named_wkt = _crs_to_named_wkt(grid.crs) - final_ds.rio.write_crs(named_wkt, inplace=True) - final_ds.rio.write_transform(grid.transform, inplace=True) + write_georeference(final_ds, grid) ctx.data = final_ds return ctx \ No newline at end of file diff --git a/disscube/pipeline/aligner.py b/disscube/pipeline/aligner.py index b1fb06f..caa3cab 100644 --- a/disscube/pipeline/aligner.py +++ b/disscube/pipeline/aligner.py @@ -28,6 +28,7 @@ import geopandas as gpd import numpy as np import rioxarray # registers the .rio accessor +import shapely import xarray as xr from pyproj import CRS as ProjCRS from pyproj.exceptions import CRSError @@ -102,6 +103,66 @@ def _with_non_colliding_nodata(band: xr.DataArray) -> xr.DataArray: return out.rio.write_nodata(np.iinfo(wider).max, encoded=False) +def _clip_to_grid(band: xr.DataArray, grid: GridSpec, pad: int = 2) -> xr.DataArray: + """ + Keep only the pixels of ``band`` that can reach the cells of ``grid``. + + The band is opened lazily, so slicing it reads just that window: a 9.7-billion-pixel raster + aligned to a small grid no longer has to be loaded whole. ``pad`` extra source pixels on each side + cover the pixels that only partly overlap the edge cells, so the aligned values are the same as + from the whole band. A grid with no overlap with the band gets a one-pixel window at the nearest + edge, which aligns to nodata. Bands that are not north-up are returned as they are. + """ + try: + from rasterio.warp import transform_bounds + from rasterio.windows import from_bounds + + crs, tr = band.rio.crs, band.rio.transform() + ny, nx = band.sizes["y"], band.sizes["x"] + if crs is None or tr.b != 0 or tr.d != 0 or tr.e >= 0 or tr.a <= 0: + return band + left, bottom, right, top = transform_bounds(grid.crs, crs, *grid.bbox, densify_pts=21) + w = from_bounds(left, bottom, right, top, transform=tr) + c0, c1 = int(np.floor(w.col_off)) - pad, int(np.ceil(w.col_off + w.width)) + pad + r0, r1 = int(np.floor(w.row_off)) - pad, int(np.ceil(w.row_off + w.height)) + pad + c0, c1, r0, r1 = max(c0, 0), min(c1, nx), max(r0, 0), min(r1, ny) + if c1 <= c0: + c0 = min(max(c0, 0), nx - 1); c1 = c0 + 1 + if r1 <= r0: + r0 = min(max(r0, 0), ny - 1); r1 = r0 + 1 + if (c0, c1, r0, r1) == (0, nx, 0, ny): + return band + return band.isel(y=slice(r0, r1), x=slice(c0, c1)) + except Exception: + log.debug("could not clip the source to the grid; using the whole band", exc_info=True) + return band + + + +def _repair_geometries(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame: + """Repair the invalid geometries of a vector source (self-intersections, often created by the reprojection). + + ``clip`` raises ``TopologyException`` on an invalid polygon, which stopped a pipeline for one bad feature in a layer + of thousands. Only the invalid features are changed; a polygon that becomes a collection keeps its polygonal parts. + """ + geoms = gdf.geometry + invalid = ~geoms.is_valid & geoms.notna() + if not invalid.any(): + return gdf + fixed = shapely.make_valid(geoms[invalid].values) + keep_polygons = geoms[invalid].geom_type.isin(["Polygon", "MultiPolygon"]).values + out = [] + for geom, polygonal in zip(fixed, keep_polygons): + if polygonal and geom.geom_type not in ("Polygon", "MultiPolygon"): + parts = [g for g in shapely.get_parts(geom) if g.geom_type in ("Polygon", "MultiPolygon")] + geom = shapely.union_all(parts) if parts else shapely.Polygon() + out.append(geom) + gdf = gdf.copy() + gdf.loc[invalid.values, gdf.geometry.name] = gpd.GeoSeries(out, index=gdf.index[invalid.values], crs=gdf.crs) + log.warning("repaired %d invalid geometries before clipping", int(invalid.sum())) + return gdf + + class GridAligner(PipelineStage): def execute(self, ctx: PipelineContext) -> PipelineContext: grid = ctx.grid @@ -125,6 +186,7 @@ def execute(self, ctx: PipelineContext) -> PipelineContext: needs_reproject = str(gdf.crs) != str(grid.crs) if needs_reproject: gdf = gdf.to_crs(grid.crs) + gdf = _repair_geometries(gdf) ops = [OPERATOR_REGISTRY.get(v.operator) for v in ctx.derivation.variables] if all(getattr(op, "clip_to_grid", True) for op in ops): gdf = gdf.clip(box(*grid.bbox)) @@ -170,105 +232,134 @@ def _align_raster( # file's time metadata, and units such as "years since 2000-1-1" do not # decode; the time of a source is set when it is registered. opened = rioxarray.open_rasterio(url, decode_times=False) - if not isinstance(opened, xr.DataArray): - raise TypeError( - f"{url!r} opened as {type(opened).__name__}, not a single raster; " - "a multi-subdataset file needs an explicit variable (e.g. NETCDF:file:var)" - ) - ds_src = opened - # GDAL names the band axis of a NetCDF variable after its own - # dimension (e.g. "time"); every band of a raster is a "band" here. - extra = [d for d in ds_src.dims if d not in ("y", "x")] - if len(extra) == 1 and extra[0] != "band": - ds_src = ds_src.rename({extra[0]: "band"}) - ds_src = ds_src.assign_coords(band=np.arange(1, ds_src.sizes["band"] + 1)) - if nodata is not None: - ds_src = ds_src.rio.write_nodata(nodata, encoded=False) - # Map of variable name -> aligned DataArray. A plain dict (not a - # Dataset) is used because fine-aligned categorical arrays have a - # different shape than the target grid; putting them in a Dataset - # keyed on grid coords would trigger coordinate realignment to NaN. - result: dict[str, xr.DataArray] = {} - - for i, var in enumerate(variables): - # ── Band selection ───────────────────────────────────────── - if "band" in ds_src.dims and ds_src.sizes["band"] > 1: - if band_map and var.name in band_map: - band_idx = band_map[var.name] - 1 # 1-based → 0-based - if not (0 <= band_idx < ds_src.sizes["band"]): + # Close the file when done: the aligned arrays are in memory, and a handle left open is + # finalised at interpreter exit, which printed `Error in sys.excepthook` after every run. + try: + if not isinstance(opened, xr.DataArray): + raise TypeError( + f"{url!r} opened as {type(opened).__name__}, not a single raster; " + "a multi-subdataset file needs an explicit variable (e.g. NETCDF:file:var)" + ) + ds_src = opened + # GDAL names the band axis of a NetCDF variable after its own + # dimension (e.g. "time"); every band of a raster is a "band" here. + extra = [d for d in ds_src.dims if d not in ("y", "x")] + if len(extra) == 1 and extra[0] != "band": + ds_src = ds_src.rename({extra[0]: "band"}) + ds_src = ds_src.assign_coords(band=np.arange(1, ds_src.sizes["band"] + 1)) + if nodata is not None: + ds_src = ds_src.rio.write_nodata(nodata, encoded=False) + # Map of variable name -> aligned DataArray. A plain dict (not a + # Dataset) is used because fine-aligned categorical arrays have a + # different shape than the target grid; putting them in a Dataset + # keyed on grid coords would trigger coordinate realignment to NaN. + result: dict[str, xr.DataArray] = {} + + for i, var in enumerate(variables): + # ── Band selection ───────────────────────────────────────── + if "band" in ds_src.dims and ds_src.sizes["band"] > 1: + if band_map and var.name in band_map: + band_idx = band_map[var.name] - 1 # 1-based → 0-based + if not (0 <= band_idx < ds_src.sizes["band"]): + raise ValueError( + f"Band index {band_idx + 1} for variable " + f"'{var.name}' is out of range; " + f"source has {ds_src.sizes['band']} bands." + ) + elif i < ds_src.sizes["band"]: + band_idx = i + else: raise ValueError( - f"Band index {band_idx + 1} for variable " - f"'{var.name}' is out of range; " - f"source has {ds_src.sizes['band']} bands." + f"No band available for variable '{var.name}' at " + f"index {i}; source has {ds_src.sizes['band']} bands " + "and no band_map was provided." ) - elif i < ds_src.sizes["band"]: - band_idx = i + band = ds_src.isel(band=band_idx) else: - raise ValueError( - f"No band available for variable '{var.name}' at " - f"index {i}; source has {ds_src.sizes['band']} bands " - "and no band_map was provided." + band = ds_src.isel(band=0) if "band" in ds_src.dims else ds_src + + whole_band = band + band = _clip_to_grid(band, grid) + + # ── Per-operator resampling method ───────────────────────── + op_cls = OPERATOR_REGISTRY.get(var.operator) + needs_fine = bool(getattr(op_cls, "needs_fine_alignment", False)) + + if needs_fine: + # Categorical operators must see sub-cell class composition. + # Reproject with NEAREST (never average a class code) onto a + # fine grid that shares the target grid origin, at a resolution + # that is an integer sub-multiple of the target cell size. + # The fine resolution is chosen from the whole source, not from the clipped window: + # the resolution GDAL proposes depends on the extent, and every strip of a grid must + # use the one the whole grid would. + aligned = self._align_fine(band, grid, var.params.get("subcells"), + src_res=self._source_resolution(whole_band, grid)) + result[var.name] = aligned + log.debug( + "fine-aligned '%s' via '%s' (nearest, fine shape=%s -> target=%s)", + var.name, var.operator, aligned.rio.shape, (grid.rows, grid.cols), ) - band = ds_src.isel(band=band_idx) - else: - band = ds_src.isel(band=0) if "band" in ds_src.dims else ds_src - - # ── Per-operator resampling method ───────────────────────── - op_cls = OPERATOR_REGISTRY.get(var.operator) - needs_fine = bool(getattr(op_cls, "needs_fine_alignment", False)) - - if needs_fine: - # Categorical operators must see sub-cell class composition. - # Reproject with NEAREST (never average a class code) onto a - # fine grid that shares the target grid origin, at a resolution - # that is an integer sub-multiple of the target cell size. - aligned = self._align_fine(band, grid, var.params.get("subcells")) - result[var.name] = aligned - log.debug( - "fine-aligned '%s' via '%s' (nearest, fine shape=%s -> target=%s)", - var.name, var.operator, aligned.rio.shape, (grid.rows, grid.cols), + continue + + resampling: Resampling = ( + op_cls.resampling() if op_cls else Resampling.nearest ) - continue - resampling: Resampling = ( - op_cls.resampling() if op_cls else Resampling.nearest - ) + # ── Reproject to target grid ─────────────────────────────── + # Continuous operators work in float64 with NaN as nodata (see + # _as_float_with_nan_nodata for why the source dtype is unsafe). + band = _as_float_with_nan_nodata(band) + aligned = band.rio.reproject( + grid.crs, + shape=(grid.rows, grid.cols), + transform=grid.transform, + resampling=resampling, + ) - # ── Reproject to target grid ─────────────────────────────── - # Continuous operators work in float64 with NaN as nodata (see - # _as_float_with_nan_nodata for why the source dtype is unsafe). - band = _as_float_with_nan_nodata(band) - aligned = band.rio.reproject( - grid.crs, - shape=(grid.rows, grid.cols), - transform=grid.transform, - resampling=resampling, - ) + # ── Alignment invariant ──────────────────────────────────── + actual = aligned.rio.shape + expected = (grid.rows, grid.cols) + if actual != expected: + raise ValueError( + f"GridAligner: variable '{var.name}' alignment produced " + f"shape {actual}, expected {expected} for grid '{grid.id}'. " + f"Source: {url}" + ) - # ── Alignment invariant ──────────────────────────────────── - actual = aligned.rio.shape - expected = (grid.rows, grid.cols) - if actual != expected: - raise ValueError( - f"GridAligner: variable '{var.name}' alignment produced " - f"shape {actual}, expected {expected} for grid '{grid.id}'. " - f"Source: {url}" + result[var.name] = aligned.transpose("y", "x") + log.debug( + "aligned '%s' via '%s' (resampling=%s, shape=%s)", + var.name, var.operator, resampling.name, actual, ) - result[var.name] = aligned.transpose("y", "x") - log.debug( - "aligned '%s' via '%s' (resampling=%s, shape=%s)", - var.name, var.operator, resampling.name, actual, - ) - - return result + return result + finally: + # open_rasterio returns a DataArray, a Dataset or a list of Datasets (multi-subdataset files) + for handle in opened if isinstance(opened, list) else [opened]: + handle.close() # ------------------------------------------------------------------ # Fine alignment for categorical operators # ------------------------------------------------------------------ + @staticmethod + def _source_resolution(band: xr.DataArray, grid: GridSpec) -> float: + """The source resolution in the grid CRS: the one GDAL would choose to reproject ``band``.""" + from rasterio.warp import calculate_default_transform + + try: + rows, cols = band.rio.shape + src_transform, _, _ = calculate_default_transform( + band.rio.crs, grid.crs, cols, rows, *band.rio.bounds() + ) + return abs(float(src_transform.a)) + except Exception: # noqa: BLE001 — defensive fallback: the .rio accessor raises undocumented types (e.g. ValueError, CRSError) + return grid.resolution + def _align_fine( - self, band: xr.DataArray, grid: GridSpec, max_subcells: int | None = None + self, band: xr.DataArray, grid: GridSpec, max_subcells: int | None = None, + src_res: float | None = None, ) -> xr.DataArray: """ Reproject ``band`` onto a fine grid snapped to the target grid origin. @@ -304,20 +395,9 @@ def _align_fine( """ from affine import Affine - # The source resolution in target CRS units: the resolution GDAL would - # choose to reproject it (the same one ``rio.reproject`` uses by - # default), computed without reprojecting the pixels. - from rasterio.warp import calculate_default_transform - band = _with_non_colliding_nodata(band) - try: - rows, cols = band.rio.shape - src_transform, _, _ = calculate_default_transform( - band.rio.crs, grid.crs, cols, rows, *band.rio.bounds() - ) - src_res = abs(float(src_transform.a)) - except Exception: # noqa: BLE001 — defensive fallback: the .rio accessor raises undocumented types (e.g. ValueError, CRSError) - src_res = grid.resolution + if src_res is None: + src_res = self._source_resolution(band, grid) target_res = grid.resolution if src_res <= 0 or src_res >= target_res: diff --git a/disscube/pipeline/runner.py b/disscube/pipeline/runner.py index e1b978d..4e75640 100644 --- a/disscube/pipeline/runner.py +++ b/disscube/pipeline/runner.py @@ -9,6 +9,7 @@ import tomllib from dataclasses import dataclass, field from datetime import UTC, datetime +from importlib.metadata import PackageNotFoundError, version from pathlib import Path from typing import Any, Literal @@ -67,7 +68,7 @@ class PlannedDerive: target: str source: str operator: str - class_code: int | None + class_code: int | list[int] | None role: str params: dict = field(default_factory=dict) fill: str | None = None @@ -279,6 +280,14 @@ def resolve_workspace(p: Plan, workspace: str | Path | None = None) -> Path: return pf.base_dir / cfg.workspace if cfg.workspace else pf.base_dir / pf.path.stem +def _installed_version() -> str | None: + """Version of the installed DisSCube, written to ``run.json``; ``None`` from a source tree that is not installed.""" + try: + return version("disscube") + except PackageNotFoundError: + return None + + def run(pipeline: PipelineFile | Plan | str | Path, workspace: str | Path | None = None, export_geotiff: str | Path | None = None) -> RunReport: """ @@ -342,6 +351,7 @@ def run(pipeline: PipelineFile | Plan | str | Path, workspace: str | Path | None report.exported = out_path record = { + "disscube": _installed_version(), "pipeline": pipeline_info, "started": started, "finished": datetime.now(UTC).isoformat(timespec="seconds"), @@ -449,6 +459,7 @@ def _register_source( raw, files=files, name=c.name, + codes=c.codes, ) if isinstance(c, OsmSource): from disscube.sources import osm diff --git a/disscube/pipeline/schema.py b/disscube/pipeline/schema.py index 3b61829..83f6ff8 100644 --- a/disscube/pipeline/schema.py +++ b/disscube/pipeline/schema.py @@ -118,6 +118,19 @@ class ProdesSource(_SourceBase): year: int | str | None = None url: str | None = None cache: str | None = None + codes: dict[str, int] | None = None + """Local code of each category (forest, deforested, non_forest, water); categories with one code are merged.""" + + @model_validator(mode="after") + def _check_codes(self): + if self.codes is not None: + from disscube.sources.prodes import resolve_codes + + try: + resolve_codes(self.codes) + except ValueError as e: + raise ValueError(f"source {self.id!r}: {e}") from e + return self class ClassifiedSource(_SourceBase): @@ -193,7 +206,7 @@ class DeriveConfig(_Strict): target: str source: str operator: str - class_code: int | None = None + class_code: int | list[int] | None = None role: str = "driver" years: list[int] | None = None params: dict[str, Any] = Field(default_factory=dict) diff --git a/disscube/plot.py b/disscube/plot.py new file mode 100644 index 0000000..4489be3 --- /dev/null +++ b/disscube/plot.py @@ -0,0 +1,136 @@ +""" +Quick look at the bands of an exported GeoTIFF: one panel per band, written as a PNG. + + disscube plot cellspace.tif [-o look.png] [--bands NAME ...] [--year 2024] [--cols 2] [--shared-colorbar] [--crop] [--dpi 150] + +A check of the result, not a publication figure. The scale follows the values: shares in [0, 1] are +drawn as 0 to 100 %, a band of a few integer codes as discrete classes, anything else as a continuous +scale clipped at the 2nd and 98th percentiles. Cells without data are left white. ``--shared-colorbar`` draws one bar for all panels when every +panel is a share (0 to 100 %), to leave more room for the maps. ``--crop`` trims the empty margin: the panels are cut to the +bounding box of the cells with data in any of the drawn bands. Needs matplotlib +(``pip install 'disscube[plot]'``). +""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np + +#: Longest side, in pixels, of a panel: a larger raster is read decimated. +MAX_SIDE = 1600 +#: A band of integer codes with at most this many distinct values is drawn as classes. +MAX_CLASSES = 12 + + +def _kind(values: np.ndarray) -> str: + """'share' (all in [0, 1], not integer codes), 'classes' or 'continuous'.""" + v = values[np.isfinite(values)] + if v.size == 0: + return "empty" + integer = bool(np.all(v == np.round(v))) + unique = np.unique(v) + if integer and unique.size <= MAX_CLASSES and (v.max() > 1 or unique.size > 2): + return "classes" + if v.min() >= 0.0 and v.max() <= 1.0: + return "share" + return "continuous" + + +def _crop(data: list[np.ndarray], extent: list[float]) -> tuple[list[np.ndarray], list[float]]: + """Cut the panels to the bounding box of the cells with data in any band (extent: left, right, bottom, top).""" + valid = np.any([np.isfinite(a) for a in data], axis=0) + rows, cols = np.where(valid.any(axis=1))[0], np.where(valid.any(axis=0))[0] + if rows.size == 0: + return data, extent + r0, r1, c0, c1 = rows[0], rows[-1] + 1, cols[0], cols[-1] + 1 + h, w = valid.shape + left, right, bottom, top = extent + dx, dy = (right - left) / w, (top - bottom) / h + cropped = [a[r0:r1, c0:c1] for a in data] + return cropped, [left + c0 * dx, left + c1 * dx, top - r1 * dy, top - r0 * dy] + + +def plot_geotiff(path: str | Path, output: str | Path | None = None, *, bands: list[str] | None = None, + year: int | None = None, cols: int | None = None, dpi: int = 150, shared_colorbar: bool = False, crop: bool = False) -> Path: + """Draw the bands of ``path`` (all, or those named in ``bands`` or of ``year``) into a PNG.""" + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + except ImportError: + raise ImportError("disscube plot needs matplotlib: pip install 'disscube[plot]'") from None + import rasterio + from rasterio.enums import Resampling + + path = Path(path) + out = Path(output) if output else path.with_suffix(".png") + with rasterio.open(path) as src: + names = [d or f"band {i}" for i, d in enumerate(src.descriptions, 1)] + tags = [src.tags(i) for i in range(1, src.count + 1)] + chosen = list(range(src.count)) + if bands: + chosen = [i for i in chosen if names[i] in bands or tags[i].get("VARIABLE") in bands] + if year is not None: + chosen = [i for i in chosen if tags[i].get("YEAR") == str(year)] + if not chosen: + raise ValueError(f"no band of {path.name} matches; bands: {', '.join(names)}") + scale = min(1.0, MAX_SIDE / max(src.width, src.height)) + shape = (max(1, round(src.height * scale)), max(1, round(src.width * scale))) + data = [src.read(i + 1, out_shape=shape, resampling=Resampling.nearest, masked=True) + .astype("float64").filled(np.nan) for i in chosen] + t = src.transform + extent = [t.c, t.c + t.a * src.width, t.f + t.e * src.height, t.f] + if crop: + data, extent = _crop(data, extent) + crs = src.crs.to_string() if src.crs else "unknown CRS" + res = abs(t.a) + unit = (src.crs.linear_units if src.crs and src.crs.is_projected else "degrees") if src.crs else "units" + unit = {"metre": "m", "meter": "m"}.get(unit, unit) + grid = src.tags().get("GRID_ID", "") + + ncols = cols or (1 if len(data) == 1 else 2 if len(data) <= 4 else 3) + nrows = -(-len(data) // ncols) + fig, axes = plt.subplots(nrows, ncols, figsize=(5.2 * ncols, 4.4 * nrows), squeeze=False, + constrained_layout=True) + cmaps = {"share": "viridis", "continuous": "magma"} + shared = shared_colorbar and all(_kind(a) == "share" for a in data) + last_im = None + for ax, i, arr in zip(axes.flat, chosen, data, strict=False): + kind = _kind(arr) + if kind == "empty": + ax.text(0.5, 0.5, "no data", ha="center", va="center", transform=ax.transAxes) + elif kind == "classes": + codes = np.unique(arr[np.isfinite(arr)]) + cmap = plt.get_cmap("tab20", max(len(codes), 2)) + index = np.full(arr.shape, np.nan) + for k, c in enumerate(codes): + index[arr == c] = k + im = ax.imshow(index, extent=extent, cmap=cmap, vmin=-0.5, vmax=len(codes) - 0.5, + interpolation="nearest") + cb = fig.colorbar(im, ax=ax, shrink=0.7, ticks=range(len(codes))) + cb.ax.set_yticklabels([f"{c:g}" for c in codes]) + cb.set_label("class code") + elif kind == "share": + im = ax.imshow(arr * 100, extent=extent, cmap=cmaps[kind], vmin=0, vmax=100) + last_im = im + if not shared: + fig.colorbar(im, ax=ax, shrink=0.7, label="share of the cell (%)") + else: + v = arr[np.isfinite(arr)] + lo, hi = np.percentile(v, [2, 98]) + im = ax.imshow(arr, extent=extent, cmap=cmaps[kind], vmin=lo, vmax=hi if hi > lo else lo + 1) + fig.colorbar(im, ax=ax, shrink=0.7, label="2nd to 98th percentile") + ax.set_title(names[i], fontsize=9) + ax.set_xticks([]) + ax.set_yticks([]) + for ax in list(axes.flat)[len(data):]: + ax.axis("off") + if shared and last_im is not None: + fig.colorbar(last_im, ax=list(axes.flat)[:len(data)], shrink=0.6, label="share of the cell (%)") + fig.suptitle(f"{path.name} {grid} {res:g} {unit} {crs}".replace(" ", " "), fontsize=9) + out.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(out, dpi=dpi) + plt.close(fig) + return out diff --git a/disscube/sources/__init__.py b/disscube/sources/__init__.py index 8e54381..7282505 100644 --- a/disscube/sources/__init__.py +++ b/disscube/sources/__init__.py @@ -16,11 +16,13 @@ from disscube.sources._raster import ( Window2D, composite, + copy_window, is_bdc_albers, mosaic, normalized_difference, portable_crs, read_window, + register_file, register_raster, software_versions, write_geotiff, @@ -29,6 +31,7 @@ __all__ = [ "Window2D", "composite", + "copy_window", "is_bdc_albers", "mosaic", "normalized_difference", @@ -37,6 +40,7 @@ "read_qml_legend", "read_window", "reclassify", + "register_file", "register_raster", "software_versions", "write_geotiff", diff --git a/disscube/sources/_raster.py b/disscube/sources/_raster.py index 63f1bce..a9b2905 100644 --- a/disscube/sources/_raster.py +++ b/disscube/sources/_raster.py @@ -69,6 +69,14 @@ def portable_crs(crs: CRS) -> CRS: return CRS.from_string(BDC_CRS) if is_bdc_albers(crs) else crs +def _window_of(ds, bbox_geo: Sequence[float]) -> Window: + """Pixel window of the open raster ``ds`` that covers ``bbox_geo`` (WGS84), clipped to the file.""" + bounds = transform_bounds("EPSG:4326", ds.crs, *bbox_geo, densify_pts=21) + win = from_bounds(*bounds, transform=ds.transform) + win = win.round_offsets().round_lengths() + return win.intersection(Window(0, 0, ds.width, ds.height)) + + def read_window( href: str, bbox_geo: Sequence[float], @@ -87,10 +95,7 @@ def read_window( value they use, such as MapBiomas' 0). """ with rasterio.open(href) as ds: - bounds = transform_bounds("EPSG:4326", ds.crs, *bbox_geo, densify_pts=21) - win = from_bounds(*bounds, transform=ds.transform) - win = win.round_offsets().round_lengths() - win = win.intersection(Window(0, 0, ds.width, ds.height)) + win = _window_of(ds, bbox_geo) raw = ds.read(1, window=win) transform = ds.window_transform(win) nodata = nodata if nodata is not None else ds.nodata @@ -216,6 +221,66 @@ def write_geotiff( +#: Rows of pixels read and written at a time by :func:`copy_window` (a 111 508-pixel-wide +#: uint8 strip of this height is about 114 MB). +COPY_BLOCK_ROWS = 1024 + + +def copy_window( + href: str, + bbox_geo: Sequence[float], + path: str | Path, + *, + nodata: float | None = None, + lookup: dict[int, int] | None = None, + out_nodata: int = 0, + block_rows: int = COPY_BLOCK_ROWS, +) -> tuple[Path, CRS]: + """ + Write the pixels of ``href`` that cover ``bbox_geo`` to a uint8 GeoTIFF, a block of rows at a time. + + Same result as ``read_window`` + ``write_geotiff(dtype="uint8")`` for an integer source, but memory + is bounded by ``block_rows`` instead of the size of the window (the window of PRODES over the Legal + Amazon is 9.7 billion pixels; as float32 it would take 39 GB). ``nodata`` overrides the file's; + those pixels are written as ``out_nodata``. ``lookup`` maps source codes to output codes + (``reclassify``'s job in the whole-window path); a code that is not in it becomes ``out_nodata``. + Returns the file and the source CRS. + """ + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + with rasterio.open(href) as ds: + win = _window_of(ds, bbox_geo) + nd = nodata if nodata is not None else ds.nodata + lut = None + if lookup is not None: + if ds.dtypes[0] != "uint8": + raise ValueError(f"a lookup table needs a uint8 source; {href} is {ds.dtypes[0]}") + lut = np.full(256, out_nodata, dtype="uint8") + for old, new in lookup.items(): + if 0 <= int(old) < 256: + lut[int(old)] = int(new) + if nd is not None and not np.isnan(nd) and 0 <= int(nd) < 256: + lut[int(nd)] = out_nodata + transform = ds.window_transform(win) + crs = ds.crs + with rasterio.open( + path, "w", driver="GTiff", height=int(win.height), width=int(win.width), count=1, dtype="uint8", + crs=portable_crs(crs), transform=transform, nodata=out_nodata, compress="deflate", + ) as dst: + for r in range(0, int(win.height), block_rows): + h = min(block_rows, int(win.height) - r) + sub = Window(win.col_off, win.row_off + r, win.width, h) + raw = ds.read(1, window=sub) + if lut is not None: + block = lut[raw] + else: + block = raw.astype("uint8") + if nd is not None and not np.isnan(nd): + block[raw == nd] = out_nodata + dst.write(block, 1, window=Window(0, r, int(win.width), h)) + return path, crs + + # --------------------------------------------------------------------------- # Registration with provenance # --------------------------------------------------------------------------- @@ -243,10 +308,27 @@ def register_raster( ``dtype``/``nodata`` are passed to :func:`write_geotiff` (integer types for categorical layers). Returns the ``SpatialSource``. """ - from disscube.models import SpatialSource - out_dir = Path(out_dir) tif = write_geotiff(window, out_dir / f"{source_id}.tif", dtype=dtype, nodata=nodata) + return register_file(cube, source_id, tif, window.crs, provenance, name=name, time=time, tags=tags) + + +def register_file( + cube, + source_id: str, + tif: str | Path, + crs: CRS, + provenance: dict, + *, + name: str | None = None, + time: int | None = None, + tags: Sequence[str] = (), +): + """Register a GeoTIFF already written by an adapter: checksum, provenance sidecar, catalog entry.""" + from disscube.models import SpatialSource + + tif = Path(tif) + out_dir = tif.parent checksum = sha256_file(tif) record = { **provenance, @@ -266,7 +348,7 @@ def register_raster( format="raster", asset_url=str(tif), checksum=checksum, - crs=BDC_CRS if is_bdc_albers(window.crs) else window.crs.to_string(), + crs=BDC_CRS if is_bdc_albers(crs) else crs.to_string(), time=time, tags=[*tags, f"provenance:{prov_path}"], ) diff --git a/disscube/sources/mapbiomas.py b/disscube/sources/mapbiomas.py index 2b7b33c..84a1518 100644 --- a/disscube/sources/mapbiomas.py +++ b/disscube/sources/mapbiomas.py @@ -28,7 +28,7 @@ class in ``percentage``/``majority``. Data © MapBiomas, CC-BY-4.0. from dataclasses import dataclass from pathlib import Path -from disscube.sources._raster import Window2D, read_window, register_raster +from disscube.sources._raster import Window2D, copy_window, read_window, register_file log = logging.getLogger(__name__) @@ -141,7 +141,6 @@ def register_mapbiomas_source( resolution, year, URL, bbox, legend and license). ``time`` is ``year``. """ href = url or coverage_url(year, collection, resolution) - window = read_coverage(year, bbox_geo, collection=collection, resolution=resolution, url=href) provenance = { "dataset": "MapBiomas Brazil — annual land cover and land use", "collection": collection, @@ -153,10 +152,13 @@ def register_mapbiomas_source( "legend": LEGEND_URL, "license": LICENSE, } - return register_raster( - cube, source_id, window, out_dir, provenance, + # Read and written a block of rows at a time, so the window can be larger than memory (MapBiomas + # over Brazil at 30 m is 9 billion pixels). Same file as read_coverage() + write_geotiff(). + tif, crs = copy_window(href, bbox_geo, Path(out_dir) / f"{source_id}.tif", + nodata=MAPBIOMAS_NODATA, out_nodata=MAPBIOMAS_NODATA) + return register_file( + cube, source_id, tif, crs, provenance, name=name or f"MapBiomas collection {collection} ({resolution} m), {year}", time=year, tags=["mapbiomas", f"collection:{collection}", f"resolution:{resolution}m", f"year:{year}"], - dtype="uint8", nodata=MAPBIOMAS_NODATA, ) diff --git a/disscube/sources/prodes.py b/disscube/sources/prodes.py index 3f1e859..c24db72 100644 --- a/disscube/sources/prodes.py +++ b/disscube/sources/prodes.py @@ -14,7 +14,9 @@ Codes of the local map: :data:`FOREST` (1), :data:`DEFORESTED` (2), :data:`OTHER` (3, non-forest vegetation and water), nodata 0 (clouds and -anything the legend does not describe). +anything the legend does not describe). The code of each category is +a choice (``codes=``, or ``[source.codes]`` in a pipeline file): categories given the same code are merged, so +``water = 4`` keeps the hydrography of the edition apart from the non-forest vegetation. **Years.** The edition distributed today labels deforestation by year from ``d2007`` on — the first of those codes holds everything deforested up to @@ -32,13 +34,13 @@ import tempfile import urllib.request import zipfile -from collections.abc import Sequence +from collections.abc import Mapping, Sequence from dataclasses import dataclass, field from pathlib import Path from urllib.parse import urljoin, urlparse from disscube.sources._categorical import read_qml_legend, reclassify, strip_code -from disscube.sources._raster import Window2D, read_window, register_raster +from disscube.sources._raster import Window2D, copy_window, read_window, register_file from disscube.utils import sha256_file log = logging.getLogger(__name__) @@ -51,6 +53,25 @@ "prodes_amazonia_legal_2025_v20260408.zip") FOREST, DEFORESTED, OTHER = 1, 2, 3 + +#: What the legend of an edition is read as, and the local code each category gets by default. +#: ``other`` (3) of the default map is ``non_forest`` and ``water`` together. +CATEGORIES = ("forest", "deforested", "non_forest", "water") +DEFAULT_CODES = {"forest": FOREST, "deforested": DEFORESTED, "non_forest": OTHER, "water": OTHER} + + +def resolve_codes(codes: Mapping[str, int] | None = None) -> dict[str, int]: + """The local code of every category: ``codes`` over :data:`DEFAULT_CODES`, checked.""" + codes = dict(codes or {}) + unknown = sorted(set(codes) - set(CATEGORIES)) + if unknown: + raise ValueError(f"unknown PRODES categories {unknown}; the categories are {list(CATEGORIES)}") + out = {**DEFAULT_CODES, **codes} + for name, code in out.items(): + if isinstance(code, bool) or not isinstance(code, int) or not 1 <= code <= 255: + raise ValueError(f"PRODES code of {name!r} must be an integer from 1 to 255 (0 is no data), got {code!r}") + return out + NODATA = 0 @@ -159,8 +180,14 @@ def first_year(self) -> int: raise ValueError("the legend has no yearly deforestation classes (dYYYY)") return min(self.deforestation) - def lookup(self, year: int) -> dict[int, int]: - """Code -> FOREST / DEFORESTED / OTHER for the situation at the end of ``year``.""" + def lookup(self, year: int, codes: Mapping[str, int] | None = None) -> dict[int, int]: + """ + Code -> local code for the situation at the end of ``year``. + + ``codes`` gives the local code of each category (:func:`resolve_codes`): by default forest 1, + deforested 2, non-forest vegetation and water 3. + """ + local = resolve_codes(codes) if year < self.first_year: raise ValueError( f"this PRODES edition separates deforestation from {self.first_year} on " @@ -169,11 +196,16 @@ def lookup(self, year: int) -> dict[int, int]: ) table: dict[int, int] = {} for code in self.forest: - table[code] = FOREST - for code in self.non_forest | self.water: - table[code] = OTHER - for yr, code in {**self.deforestation, **self.residual}.items(): - table[code] = DEFORESTED if yr <= year else FOREST + table[code] = local["forest"] + for code in self.non_forest: + table[code] = local["non_forest"] + for code in self.water: + table[code] = local["water"] + # Two separate loops: the yearly classes (d2015) and the residual ones (r2015) share their years, so + # merging the two dicts by year would drop every dYYYY that has an rYYYY. + for by_year in (self.deforestation, self.residual): + for yr, code in by_year.items(): + table[code] = local["deforested"] if yr <= year else local["forest"] return table @@ -204,13 +236,14 @@ def parse_legend(legend: dict[int, str]) -> Scheme: # Read and register # --------------------------------------------------------------------------- -def read_prodes(year: int, bbox_geo: Sequence[float], files: ProdesFiles) -> tuple[Window2D, Scheme]: - """Forest / deforested / other at the end of ``year`` over ``bbox_geo``.""" +def read_prodes(year: int, bbox_geo: Sequence[float], files: ProdesFiles, + codes: Mapping[str, int] | None = None) -> tuple[Window2D, Scheme]: + """The PRODES classes at the end of ``year`` over ``bbox_geo`` (``codes``: see :meth:`Scheme.lookup`).""" if files.qml is None: raise ValueError(f"{files.zip_path.name} has no .qml legend; cannot interpret the codes") scheme = parse_legend(read_qml_legend(files.qml)) raw = read_window(str(files.tif), bbox_geo) - return reclassify(raw, scheme.lookup(year)), scheme + return reclassify(raw, scheme.lookup(year, codes)), scheme def register_prodes_source( @@ -224,16 +257,23 @@ def register_prodes_source( cache_dir: str | Path | None = None, files: ProdesFiles | None = None, name: str | None = None, + codes: Mapping[str, int] | None = None, ): """ Register PRODES forest / deforested / other at the end of ``year`` as a source. + ``codes``: the local code of each category (``forest``, ``deforested``, ``non_forest``, ``water``); + categories with the same code are merged. Default: 1, 2, 3, 3. ``codes={"water": 4}`` keeps the + hydrography apart from the non-forest vegetation. + Downloads the edition on first use (see :func:`download`). The provenance records the edition URL, the ZIP's SHA-256, the legend as read from the ``.qml`` and the code mapping applied. """ files = files or download(url, cache_dir) - window, scheme = read_prodes(year, bbox_geo, files) + if files.qml is None: + raise ValueError(f"{files.zip_path.name} has no .qml legend; cannot interpret the codes") + scheme = parse_legend(read_qml_legend(files.qml)) provenance = { "dataset": "PRODES — Legal Amazon deforestation (INPE)", "edition": files.zip_path.stem, @@ -241,14 +281,17 @@ def register_prodes_source( "zip_checksum": files.zip_checksum, "year": year, "bbox_geo": list(bbox_geo), - "codes": {"forest": FOREST, "deforested": DEFORESTED, "other": OTHER, "nodata": NODATA}, + "codes": {**resolve_codes(codes), "nodata": NODATA}, "first_separable_year": scheme.first_year, "legend": {str(k): v for k, v in read_qml_legend(files.qml).items()} if files.qml else None, "source": "https://terrabrasilis.dpi.inpe.br/downloads/", } - return register_raster( - cube, source_id, window, out_dir, provenance, + # Read and written a block of rows at a time: the window can be the whole Legal Amazon (9.7 billion + # pixels), which does not fit in memory as one array. Same file as read_prodes() + write_geotiff(). + tif, crs = copy_window(str(files.tif), bbox_geo, Path(out_dir) / f"{source_id}.tif", + lookup=scheme.lookup(year, codes), out_nodata=NODATA) + return register_file( + cube, source_id, tif, crs, provenance, name=name or f"PRODES {files.zip_path.stem}, situation at the end of {year}", - time=year, tags=["prodes", f"edition:{files.zip_path.stem}", f"year:{year}"], - dtype="uint8", nodata=NODATA, + time=year, tags=["prodes", f"edition:{files.zip_path.stem}", f"year:{year}", *(["codes:custom"] if codes else [])], ) diff --git a/docs/guides/landcover.md b/docs/guides/landcover.md index 95a7b4f..1fdec18 100644 --- a/docs/guides/landcover.md +++ b/docs/guides/landcover.md @@ -28,15 +28,55 @@ for year in (2008, 2016, 2024): - **Legend from the file.** Codes are read from the `.qml` (`"7 d2007"`, `"91 Hidrografia"`, `"100 Vegetação nativa florestal"`, …), not hard-coded, so a new edition with another encoding is read correctly or fails loudly. -- **One map per year.** For a year Y the local map has `FOREST` (1), - `DEFORESTED` (2, any `dYYYY` or residual `rYYYY` up to Y) and `OTHER` (3, - non-forest vegetation and water); clouds and unknown codes are nodata. + Each label is read as one category, by its text (the first rule that fits wins): + + | Category | Label | + |---|---| + | `deforested` | `dAAAA` (deforestation of the year, e.g. `d2007`), and `rAAAA` (residual, e.g. `r2015`) | + | `water` | contains "hidrografia", "hydrograph" or "água" | + | clouds | contains "nuvem" or "cloud": no data | + | `non_forest` | contains "não florest", "nao florest", "non-forest" or "non forest" | + | `forest` | contains "florest" or "forest" | + | anything else | warning in the log; no data | + + For the edition `prodes_amazonia_legal_2025_v20260408` the legend is: 7–25 + `d2007`…`d2025`, 50–65 `r2010`…`r2025`, 91 Hidrografia, 100 Vegetação nativa + florestal, 101 Vegetação nativa não florestal. The legend of an edition is + kept in the provenance file. +- **One map per year, with the codes you choose.** For a year Y a pixel is + `deforested` when its label is `dAAAA` or `rAAAA` with AAAA ≤ Y (a later one is + still `forest`). The local code of each category is `codes`; the default is + `forest` 1, `deforested` 2, and `non_forest` and `water` together 3 + (non-forest vegetation and water). Categories given the same code are merged, so + `codes={"water": 4}` keeps the hydrography apart: + + ```python + prodes.register_prodes_source(cube, "prodes_2024", 2024, bbox, out_dir="raw", + files=files, codes={"water": 4}) + ``` + + In a pipeline file the same is a table of the source: + + ```toml + [[source]] + id = "prodes_{year}" + type = "prodes" + years = [2024] + + [source.codes] # local code of each category (1 to 255; 0 is no data) + forest = 1 # vegetação nativa florestal + deforested = 2 # desmatamento até o ano (dAAAA and rAAAA) + non_forest = 3 # vegetação nativa não florestal + water = 4 # hidrografia; 3 would merge it with the non-forest + ``` + + Only these four names are accepted. Clouds and unknown labels are always no data. - **Years.** The edition distributed today separates deforestation by year from 2007 on (`d2007` holds everything cleared up to then). Earlier years raise an error; use MapBiomas for them. The provenance records the edition, its URL, the ZIP's SHA-256, the legend and -the code mapping. +the code of each category. ## Any classified map, e.g. from SITS (`disscube.sources.classified`) diff --git a/docs/guides/pipeline_files.md b/docs/guides/pipeline_files.md index 27077b4..dbc7f2f 100644 --- a/docs/guides/pipeline_files.md +++ b/docs/guides/pipeline_files.md @@ -96,7 +96,7 @@ relative to the pipeline file, so a folder with its TOML and data is portable. | `file` | `path` (raster, vector, or `zip://…` shapefile), `crs`, `format`, `time`, `variable`, `nodata`, `read` | a local file, checksummed | | `bdc` | `collection`, `asset` **or** `normalized_difference = [a, b]`, `period`, `reducer`, `scale`, `offset`, `url`, `time` | `disscube.sources.bdc` | | `mapbiomas` | `year`, `collection` (11), `resolution` (30), `url` | `disscube.sources.mapbiomas` | -| `prodes` | `year`, `url`, `cache` | `disscube.sources.prodes` | +| `prodes` | `year`, `url`, `cache`, `codes` | `disscube.sources.prodes` | | `classified` | `path`, `legend` (table or `.qml`/`.json`/`.csv` file), `time`, `nodata`, `producer` | `disscube.sources.classified` | | `osm` | `query`, `margin`, `date`, `endpoints`, `cache`, `time` (see [OpenStreetMap](osm.md)) | `disscube.sources.osm` | | `dem` | `dem` or `tiles`, `product`, `margin`, `resolution`, `cache` (see [Elevation and slope](dem.md)) | `disscube.sources.dem` | @@ -172,7 +172,7 @@ other loops or conditions: what does not fit belongs in Python. ├── catalog.db the catalog ├── store/ derived variables (Zarr) ├── raw/ every source: GeoTIFF + .provenance.json -└── run.json this run: pipeline file and checksum, sources, products +└── run.json this run: DisSCube version, pipeline file and checksum, sources, products ``` Each provenance file gains a `pipeline` entry with the pipeline file's path, diff --git a/docs/terrame_fill_correspondence.md b/docs/terrame_fill_correspondence.md index cc3a0e7..b327cf7 100644 --- a/docs/terrame_fill_correspondence.md +++ b/docs/terrame_fill_correspondence.md @@ -44,7 +44,7 @@ and the parity cases in the recipes repository (`cases/terrame_fill`). | TerraME fill strategy | DisSCube operator (`name`) | Status | Notes | |---|---|---|---| | `presence` | `presence` | implemented; parity measured (Emas) | Binary mask: 1 where any feature is present. Matches TerraME in 98.4–99.7 % of cells: lines are rasterized through cell centres, while TerraME marks every cell a line touches. | -| `coverage` / `percentage` (raster) | `percentage` | implemented (window-based); **parity verified** (Itaituba, Amazônia) | Fraction (0..1) of the target class per cell, **over valid pixels** — the denominator TerraME 2.0.1 uses too, so the value reproduces its goldens with no correction (Itaituba: max 0.0064; Amazônia: identical). `coverage_purity` stays available as metadata. Requires `class_code`. | +| `coverage` / `percentage` (raster) | `percentage` | implemented (window-based); **parity verified** (Itaituba, Amazônia) | Fraction (0..1) of the target class per cell, **over valid pixels** — the denominator TerraME 2.0.1 uses too, so the value reproduces its goldens with no correction (Itaituba: max 0.0064; Amazônia: identical). `coverage_purity` stays available as metadata. Requires `class_code`, one code or a list of codes (the share in any of them). | | `area` (polygons) | `area` | implemented (exact); **parity verified** (Amazônia) | Fraction (0..1) of each cell covered by polygons (e.g. protected areas): intersection area / cell area, overlapping polygons counted once. Reproduces TerraME's `protected` within 0.01 in at least 99 % of the cells. | | `majority` / `mode` | `majority` | implemented (window-based) | Dominant class per cell; ties resolve to the smallest class value. | | `minority` | `minority` | implemented (window-based) | Least-frequent class per cell. | diff --git a/pyproject.toml b/pyproject.toml index 4755bb9..f6427f3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "disscube" -version = "0.5.0" +version = "0.6.0" description = "Declarative spatial data cubes: describe sources, grid and derived variables in TOML; get a cataloged, reproducible cube." readme = "README.md" requires-python = ">=3.11" @@ -66,6 +66,9 @@ bdc = [ s3 = [ "s3fs", ] +plot = [ # disscube plot: a PNG of the bands of an exported GeoTIFF + "matplotlib", +] netcdf = [ # CubeClient.export_netcdf() "h5netcdf", "h5py", diff --git a/tests/test_aligner.py b/tests/test_aligner.py index 073aa6e..c2ad789 100644 --- a/tests/test_aligner.py +++ b/tests/test_aligner.py @@ -109,3 +109,22 @@ def test_alignment_invariant_passes_on_correct_shape(tmp_path): GridAligner().execute(ctx) # must not raise assert ctx.data["v"].shape == (10, 10) + + +def test_invalid_polygon_is_repaired_before_clip(): + """A bow-tie polygon (invalid) must not stop the alignment: the aligner repairs it before ``clip``.""" + import geopandas as gpd + from shapely.geometry import Polygon + + from disscube.pipeline.aligner import _repair_geometries + + good = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]) + bowtie = Polygon([(0, 0), (4, 4), (4, 0), (0, 4)]) # self-intersecting + gdf = gpd.GeoDataFrame({"id": [1, 2]}, geometry=[good, bowtie], crs="EPSG:3857") + assert not gdf.geometry.is_valid.all() + fixed = _repair_geometries(gdf) + assert fixed.geometry.is_valid.all() + assert fixed.geometry.geom_type.isin(["Polygon", "MultiPolygon"]).all() + assert fixed.geometry.iloc[0].equals(good) # valid features are left alone + assert fixed.geometry.iloc[1].area > 0 + assert _repair_geometries(fixed) is fixed # nothing to repair: same object diff --git a/tests/test_blocks.py b/tests/test_blocks.py new file mode 100644 index 0000000..133ccb0 --- /dev/null +++ b/tests/test_blocks.py @@ -0,0 +1,166 @@ +""" +Sources and grids larger than memory: reading, writing and deriving a block at a time must give exactly +what the whole-window path gives. + +Offline: a small "national" GeoTIFF (WGS84, uint8) stands in for MapBiomas; PRODES' lookup path is +checked against ``reclassify``. +""" + +from __future__ import annotations + +import logging + +import numpy as np +import pytest +import rasterio +from rasterio.transform import from_origin + +from disscube import CubeClient, Derivation +from disscube.sources import mapbiomas, read_window, reclassify, write_geotiff +from disscube.sources._raster import copy_window +from disscube.utils.grids import register_local_grid + +BBOX = (-44.35, -2.62, -44.20, -2.47) +RES = 0.00026949458523585647 +CLASSES = np.array([0, 3, 15, 24, 33], dtype="uint8") + + +def _national(path): + rng = np.random.default_rng(7) + smooth = rng.integers(0, len(CLASSES), size=(100, 100)) + data = np.kron(smooth, np.ones((10, 10), dtype=int)) # patches of 10 x 10 pixels + data = CLASSES[data] + with rasterio.open(path, "w", driver="GTiff", height=1000, width=1000, count=1, dtype="uint8", + crs="EPSG:4326", transform=from_origin(BBOX[0] - 0.02, BBOX[3] + 0.02, RES, RES)) as dst: + dst.write(data, 1) + return str(path) + + +def test_copy_window_equals_read_window_and_write(tmp_path): + url = _national(tmp_path / "br.tif") + whole = write_geotiff(read_window(url, BBOX, nodata=0), tmp_path / "whole.tif", dtype="uint8", nodata=0) + blocks, _ = copy_window(url, BBOX, tmp_path / "blocks.tif", nodata=0, out_nodata=0, block_rows=37) + with rasterio.open(whole) as a, rasterio.open(blocks) as b: + assert (a.width, a.height, a.transform, a.nodata) == (b.width, b.height, b.transform, b.nodata) + assert np.array_equal(a.read(1), b.read(1)) + + +def test_copy_window_lookup_equals_reclassify(tmp_path): + url = _national(tmp_path / "br.tif") + lookup = {3: 1, 15: 2, 24: 3} # 33 is not in it: becomes nodata + whole = write_geotiff(reclassify(read_window(url, BBOX), lookup), tmp_path / "whole.tif", dtype="uint8", nodata=0) + blocks, _ = copy_window(url, BBOX, tmp_path / "blocks.tif", lookup=lookup, out_nodata=0, block_rows=50) + with rasterio.open(whole) as a, rasterio.open(blocks) as b: + assert np.array_equal(a.read(1), b.read(1)) + + +def _cube(tmp_path, name, url, resolution): + cube = CubeClient(catalog=str(tmp_path / f"{name}.db"), store=str(tmp_path / f"{name}_store")) + grid = register_local_grid(cube, name="g", bbox_geo=BBOX, resolution=resolution) + mapbiomas.register_mapbiomas_source(cube, "lulc", 2020, BBOX, tmp_path / f"{name}_raw", url=url) + return cube, grid + + +@pytest.fixture +def exact_transform(monkeypatch): + """Reproject with GDAL's exact coordinate transformer instead of its default approximate one. + + The approximate transformer (tolerance 0.125 source pixels) interpolates between points whose + positions depend on the extent being warped, so a strip and the whole grid can pick a neighbouring + pixel at a cell edge. With the exact transformer every pixel is mapped independently of the extent. + """ + import inspect + + import rasterio.warp + from rioxarray.raster_array import RasterArray + + # `tolerance` is a parameter of rasterio.warp.reproject only from rasterio 1.5 (Python >= 3.12); older + # versions ignore it, so the exact transformer cannot be selected and equality is not testable. + if "tolerance" not in inspect.signature(rasterio.warp.reproject).parameters: + pytest.skip("rasterio < 1.5 cannot reproject with GDAL's exact transformer") + + original = RasterArray.reproject + + def reproject(self, *args, **kwargs): + kwargs.setdefault("tolerance", 0.0) + return original(self, *args, **kwargs) + + monkeypatch.setattr(RasterArray, "reproject", reproject) + + +@pytest.mark.parametrize("resolution", [500.0, 300.0]) +def test_strips_give_the_same_cells_as_the_whole_grid(tmp_path, caplog, resolution, exact_transform): + url = _national(tmp_path / "br.tif") + whole, g1 = _cube(tmp_path, "whole", url, resolution) + strips, g2 = _cube(tmp_path, "strips", url, resolution) + derivations = [ + Derivation(target="major", source_id="lulc", operator="majority"), + Derivation(target="forest", source_id="lulc", operator="percentage", class_code=3), + Derivation(target="pasture", source_id="lulc", operator="percentage", class_code=15), + Derivation(target="mean", source_id="lulc", operator="mean"), + Derivation(target="total", source_id="lulc", operator="sum"), + Derivation(target="med", source_id="lulc", operator="median"), + ] + for d in derivations: + whole.derive_declarative(d, grid_id=g1.id) + with caplog.at_level(logging.INFO, logger="disscube.client.cube_client"): + caplog.clear() + strips.derive_declarative(d, grid_id=g2.id, block_pixels=40_000) # a window of 1M pixels: ~25 strips + assert any("block " in r.message for r in caplog.records), "the block path was not taken" + a, b = whole.load(d.target, grid_id=g1.id), strips.load(d.target, grid_id=g2.id) + assert a.shape == b.shape + if d.operator in ("mean", "sum"): # GDAL adds the pixels in a different order: ~1e-12 + assert np.allclose(a.values, b.values, rtol=1e-9, atol=1e-9, equal_nan=True), d.target + assert np.array_equal(np.isnan(a.values), np.isnan(b.values)) + else: # class shares, majority, median: exactly the same + assert np.array_equal(a.values, b.values, equal_nan=True), d.target + assert np.allclose(a["x"].values, b["x"].values) and np.allclose(a["y"].values, b["y"].values) + + +def test_small_sources_and_neighbourhood_operators_are_not_split(tmp_path, caplog): + url = _national(tmp_path / "br.tif") + cube, grid = _cube(tmp_path, "c", url, 500.0) + with caplog.at_level(logging.INFO, logger="disscube.client.cube_client"): + cube.derive_declarative(Derivation(target="m", source_id="lulc", operator="mean"), grid_id=grid.id) + assert not any("block " in r.message for r in caplog.records) # default budget: 50 M pixels + + +def test_strips_with_the_default_transformer_differ_in_at_most_a_few_edge_cells(tmp_path): + """With GDAL's default (approximate) transformer a strip may pick the neighbouring source pixel at + a cell edge: a small fraction of the cells, never a different layout.""" + url = _national(tmp_path / "br.tif") + whole, g1 = _cube(tmp_path, "whole", url, 300.0) + strips, g2 = _cube(tmp_path, "strips", url, 300.0) + d = Derivation(target="forest", source_id="lulc", operator="percentage", class_code=3) + whole.derive_declarative(d, grid_id=g1.id) + strips.derive_declarative(d, grid_id=g2.id, block_pixels=40_000) + a, b = whole.load("forest", grid_id=g1.id).values, strips.load("forest", grid_id=g2.id).values + assert np.array_equal(np.isnan(a), np.isnan(b)) or np.mean(np.isnan(a) != np.isnan(b)) < 0.005 + both = ~np.isnan(a) & ~np.isnan(b) + assert np.mean(np.abs(a[both] - b[both]) > 1e-9) < 0.15 + assert np.nanmax(np.abs(a - b)) < 0.2 + + +def test_blocks_are_split_by_columns_when_a_wide_projected_grid_covers_a_large_geographic_area(tmp_path): + """A strip of rows of a wide projected grid spans far more longitude/latitude than its height suggests + (curvature), so rows alone cannot bound the source window: the planner must also split columns.""" + from disscube.client import _blocks, _source_pixels + from disscube.models import GridSpec, SpatialDerivation, SpatialSource + + url = _national(tmp_path / "br.tif") + # Brazil Polyconic (EPSG:5880), 1 km cells, a grid of about 18 x 16 km around the source window + from pyproj import Transformer + x0, y0 = Transformer.from_crs("EPSG:4326", "EPSG:5880", always_xy=True).transform(BBOX[0], BBOX[1]) + grid = GridSpec(id="p", type="local", crs="EPSG:5880", resolution=1000.0, + bbox=[x0, y0, x0 + 16_000.0, y0 + 16_000.0]) + source = SpatialSource(id="s", name="s", format="raster", asset_url=url, crs="EPSG:4326") + derivation = SpatialDerivation.model_validate({ + "id": "d", "source_id": "s", "grid_id": "p", "role": "feature", + "variables": [{"name": "m", "operator": "mean", "role": "feature"}], + }) + blocks = _blocks(source, grid, derivation, block_pixels=20_000) + assert blocks is not None and len(blocks[0]) > 1, "columns were not split" + assert sum(b.rows * b.cols for row in blocks for b in row) == grid.rows * grid.cols + with rasterio.open(url) as ds: + assert all(_source_pixels(b.bbox, grid.crs, ds) <= 20_000 or (b.rows == 1 and b.cols == 1) + for row in blocks for b in row) diff --git a/tests/test_class_code_list.py b/tests/test_class_code_list.py new file mode 100644 index 0000000..969cd87 --- /dev/null +++ b/tests/test_class_code_list.py @@ -0,0 +1,62 @@ +"""`percentage` accepts a list of classes: the share of the cell in any of them (a MapBiomas level-1 class +is a group of leaf codes).""" +import numpy as np +import rasterio +from rasterio.transform import from_origin + +from disscube import CubeClient, Derivation +from disscube.sources import mapbiomas +from disscube.utils.grids import register_local_grid + +BBOX = (-44.35, -2.62, -44.20, -2.47) +RES = 0.00026949458523585647 + + +def test_percentage_of_a_list_is_the_sum_of_the_shares(tmp_path): + rng = np.random.default_rng(3) + data = np.array([0, 3, 4, 15, 24, 33], dtype="uint8")[rng.integers(0, 6, size=(100, 100))] + data = np.kron(data, np.ones((10, 10), dtype="uint8")) + url = tmp_path / "br.tif" + with rasterio.open(url, "w", driver="GTiff", height=1000, width=1000, count=1, dtype="uint8", crs="EPSG:4326", + transform=from_origin(BBOX[0] - 0.02, BBOX[3] + 0.02, RES, RES)) as dst: + dst.write(data, 1) + cube = CubeClient(catalog=str(tmp_path / "c.db"), store=str(tmp_path / "s")) + grid = register_local_grid(cube, name="g", bbox_geo=BBOX, resolution=500.0) + mapbiomas.register_mapbiomas_source(cube, "lulc", 2020, BBOX, tmp_path / "raw", url=str(url)) + for name, codes in (("a", 3), ("b", 4), ("both", [3, 4]), ("none", [])): + cube.derive_declarative(Derivation(target=name, source_id="lulc", operator="percentage", class_code=codes), + grid_id=grid.id) + a, b, both = (cube.load(n, grid_id=grid.id).values for n in ("a", "b", "both")) + assert np.allclose(both, a + b, equal_nan=True) + assert np.nanmax(both - a) > 0 + assert np.nanmax(cube.load("none", grid_id=grid.id).values) == 0.0 + + +def test_majority_with_groups_is_the_dominant_group(tmp_path): + import pytest + + rng = np.random.default_rng(5) + data = np.array([0, 3, 4, 15, 24, 33], dtype="uint8")[rng.integers(0, 6, size=(100, 100))] + data = np.kron(data, np.ones((10, 10), dtype="uint8")) + url = tmp_path / "br.tif" + with rasterio.open(url, "w", driver="GTiff", height=1000, width=1000, count=1, dtype="uint8", crs="EPSG:4326", + transform=from_origin(BBOX[0] - 0.02, BBOX[3] + 0.02, RES, RES)) as dst: + dst.write(data, 1) + cube = CubeClient(catalog=str(tmp_path / "c.db"), store=str(tmp_path / "s")) + grid = register_local_grid(cube, name="g", bbox_geo=BBOX, resolution=500.0) + mapbiomas.register_mapbiomas_source(cube, "lulc", 2020, BBOX, tmp_path / "raw", url=str(url)) + groups = {"1": [3, 4], "2": [15], "3": [24, 33]} # TOML keys are strings + cube.derive_declarative(Derivation(target="top", source_id="lulc", operator="majority", params={"groups": groups}), + grid_id=grid.id) + shares = [] + for i, codes in enumerate(groups.values()): + cube.derive_declarative(Derivation(target=f"s{i}", source_id="lulc", operator="percentage", class_code=codes), + grid_id=grid.id) + shares.append(cube.load(f"s{i}", grid_id=grid.id).values) + top = cube.load("top", grid_id=grid.id).values + expected = np.argmax(np.nan_to_num(np.stack(shares), nan=-1), axis=0) + 1.0 + valid = ~np.isnan(shares[0]) + assert np.array_equal(top[valid], expected[valid]) and set(np.unique(top[valid])) <= {1.0, 2.0, 3.0} + with pytest.raises(ValueError, match="more than one group"): + cube.derive_declarative(Derivation(target="bad", source_id="lulc", operator="majority", + params={"groups": {"1": [3], "2": [3, 4]}}), grid_id=grid.id) diff --git a/tests/test_clean_exit.py b/tests/test_clean_exit.py new file mode 100644 index 0000000..edf0d5d --- /dev/null +++ b/tests/test_clean_exit.py @@ -0,0 +1,40 @@ +"""The aligner must close the source raster once the aligned arrays are in memory. + +It used to leave the file open; the handle was finalised at interpreter exit and `disscube run` ended with +`Error in sys.excepthook:` on stderr (exit code 0, results fine). +""" +import numpy as np +import rasterio +import rioxarray +from rasterio.transform import from_origin + +from disscube import CubeClient, Derivation +from disscube.sources import mapbiomas +from disscube.utils.grids import register_local_grid + +BBOX = (-44.35, -2.62, -44.20, -2.47) +RES = 0.00026949458523585647 + + +def test_aligner_closes_the_source_file(tmp_path, monkeypatch): + url = tmp_path / "br.tif" + with rasterio.open(url, "w", driver="GTiff", height=1000, width=1000, count=1, dtype="uint8", crs="EPSG:4326", + transform=from_origin(BBOX[0] - 0.02, BBOX[3] + 0.02, RES, RES)) as dst: + dst.write(np.full((1000, 1000), 3, dtype="uint8"), 1) + cube = CubeClient(catalog=str(tmp_path / "c.db"), store=str(tmp_path / "s")) + grid = register_local_grid(cube, name="g", bbox_geo=BBOX, resolution=500.0) + mapbiomas.register_mapbiomas_source(cube, "lulc", 2020, BBOX, tmp_path / "raw", url=str(url)) + + closed, original = [], rioxarray.open_rasterio + + def spy(*args, **kwargs): + da = original(*args, **kwargs) + close = da._close + da.set_close(lambda: (closed.append(1), close())[1]) + return da + + monkeypatch.setattr("disscube.pipeline.aligner.rioxarray.open_rasterio", spy) + for op, kw in (("percentage", {"class_code": 3}), ("mean", {})): # fine and continuous paths + closed.clear() + cube.derive_declarative(Derivation(target=f"v_{op}", source_id="lulc", operator=op, **kw), grid_id=grid.id) + assert closed == [1], op diff --git a/tests/test_pipeline_file.py b/tests/test_pipeline_file.py index 10caf18..4cd6446 100644 --- a/tests/test_pipeline_file.py +++ b/tests/test_pipeline_file.py @@ -263,6 +263,13 @@ def test_cli_validate_and_run(tmp_path, capsys): assert (tmp_path / "ws" / "run.json").exists() +def test_run_json_records_the_disscube_version(tmp_path): + from importlib.metadata import version + assert cli(["run", str(QUICKSTART), "--workspace", str(tmp_path / "ws")]) == 0 + record = json.loads((tmp_path / "ws" / "run.json").read_text()) + assert record["disscube"] == version("disscube") + + def test_cli_reports_errors(tmp_path, capsys): path = _toml(tmp_path, '[[derive]]\ntarget="x"\nsource="nope"\noperator="mean"\n') assert cli(["validate", str(path)]) == 2 diff --git a/tests/test_plot.py b/tests/test_plot.py new file mode 100644 index 0000000..21782e2 --- /dev/null +++ b/tests/test_plot.py @@ -0,0 +1,92 @@ +"""disscube plot: a PNG of the bands of an exported GeoTIFF.""" + +from __future__ import annotations + +import numpy as np +import pytest +import rasterio +from rasterio.transform import from_origin + +pytest.importorskip("matplotlib") + +from disscube.cli import main +from disscube.plot import _kind, plot_geotiff + + +def _tif(path): + rng = np.random.default_rng(1) + share = rng.random((30, 40)) + codes = rng.choice([3.0, 15.0, 24.0], size=(30, 40)) + elev = rng.normal(500, 80, (30, 40)) + share[:3] = np.nan + names = ["forest_2015", "forest_2020", "lulc_2020", "elev"] + years = ["2015", "2020", "2020", None] + with rasterio.open(path, "w", driver="GTiff", height=30, width=40, count=4, dtype="float64", nodata=np.nan, + crs="EPSG:5880", transform=from_origin(0, 30_000, 1000, 1000)) as dst: + dst.update_tags(GRID_ID="g1") + for i, (name, year, arr) in enumerate(zip(names, years, (share, share, codes, elev), strict=True), 1): + dst.write(arr, i) + dst.set_band_description(i, name) + dst.update_tags(i, VARIABLE=name.split("_")[0], **({"YEAR": year} if year else {})) + return path + + +def test_kind_follows_the_values(): + assert _kind(np.array([0.0, 0.2, 0.9, np.nan])) == "share" + assert _kind(np.array([3.0, 15.0, 24.0, 3.0])) == "classes" + assert _kind(np.array([0.0, 1.0, 1.0])) == "share" # a binary mask is not a set of classes + assert _kind(np.array([210.0, 480.5, 760.2])) == "continuous" + assert _kind(np.array([np.nan, np.nan])) == "empty" + + +def test_plot_writes_one_panel_per_band(tmp_path): + out = plot_geotiff(_tif(tmp_path / "a.tif"), tmp_path / "look.png") + assert out.exists() and out.stat().st_size > 5_000 + + +def test_plot_selects_bands_by_year_name_or_variable(tmp_path): + tif = _tif(tmp_path / "a.tif") + assert plot_geotiff(tif, tmp_path / "y.png", year=2020).exists() + assert plot_geotiff(tif, tmp_path / "n.png", bands=["elev"]).exists() + assert plot_geotiff(tif, tmp_path / "v.png", bands=["forest"]).exists() # a variable name: both years + with pytest.raises(ValueError, match="no band"): + plot_geotiff(tif, tmp_path / "x.png", bands=["nothing"]) + + +def test_cli_plot_and_its_errors(tmp_path, capsys): + tif = _tif(tmp_path / "a.tif") + assert main(["plot", str(tif), "-o", str(tmp_path / "p.png")]) == 0 + assert (tmp_path / "p.png").exists() + assert main(["plot", str(tif), "--bands", "zzz"]) == 2 + assert "no band" in capsys.readouterr().err + + +def test_shared_colorbar_draws_one_bar_for_shares_only(tmp_path): + tif = _tif(tmp_path / "x.tif") + plain = plot_geotiff(tif, tmp_path / "a.png", bands=["forest_2015", "forest_2020"], cols=2) + one = plot_geotiff(tif, tmp_path / "b.png", bands=["forest_2015", "forest_2020"], cols=2, shared_colorbar=True) + assert one.stat().st_size > 0 and plain.stat().st_size > 0 + # with a class band in the mix the option is ignored (each panel keeps its own bar) + mixed = plot_geotiff(tif, tmp_path / "c.png", bands=["forest_2015", "lulc_2020"], shared_colorbar=True) + assert mixed.exists() + assert main(["plot", str(tif), "-o", str(tmp_path / "d.png"), "--shared-colorbar", "--bands", "forest_2015"]) == 0 + + +def test_crop_trims_the_empty_margin_and_keeps_the_extent(): + from disscube.plot import _crop + + a = np.full((10, 20), np.nan) + a[2:5, 4:9] = 0.5 # rows 2-4, cols 4-8 + b = np.full((10, 20), np.nan) + b[6:8, 10:12] = 0.2 # rows 6-7, cols 10-11 + (ca, cb), ext = _crop([a, b], [0.0, 200.0, 0.0, 100.0]) # 10 m cells + assert ca.shape == cb.shape == (6, 8) # rows 2-7, cols 4-11 + assert [float(v) for v in ext] == [40.0, 120.0, 20.0, 80.0] + empty = np.full((4, 4), np.nan) + assert _crop([empty], [0, 4, 0, 4])[0][0].shape == (4, 4) + + +def test_plot_with_crop_and_cli(tmp_path): + tif = _tif(tmp_path / "x.tif") # the first 3 rows have no data + assert plot_geotiff(tif, tmp_path / "c.png", bands=["forest_2015"], crop=True, shared_colorbar=True).exists() + assert main(["plot", str(tif), "-o", str(tmp_path / "d.png"), "--crop", "--bands", "forest_2015"]) == 0 diff --git a/tests/test_prodes.py b/tests/test_prodes.py index 50b06e9..9163f56 100644 --- a/tests/test_prodes.py +++ b/tests/test_prodes.py @@ -98,6 +98,18 @@ def test_legend_parsing(): assert scheme.first_year == 2007 +def test_lookup_keeps_the_yearly_class_when_a_residual_class_has_the_same_year(): + """In the real legend every year from 2010 has both dYYYY (code 10 to 25) and rYYYY (50 to 65).""" + legend = {9: "9 d2009", 10: "10 d2010", 15: "15 d2015", 24: "24 d2024", 25: "25 d2025", 50: "50 r2010", + 55: "55 r2015", 64: "64 r2024", 100: "100 Vegetação nativa florestal"} + scheme = prodes.parse_legend(legend) + t = scheme.lookup(2015) + assert t[10] == t[50] == t[15] == t[55] == prodes.DEFORESTED + assert t[24] == t[64] == t[25] == prodes.FOREST # not yet cleared in 2015 + assert set(legend) - {100} <= set(t) # no class is left out (it would be nodata) + assert scheme.lookup(2024)[24] == scheme.lookup(2024)[64] == prodes.DEFORESTED + + def test_lookup_by_year(): scheme = prodes.parse_legend(LEGEND) t2012 = scheme.lookup(2012) @@ -156,3 +168,39 @@ def test_register_before_first_year_fails(server, tmp_path): files = prodes.download(server["url"], cache_dir=tmp_path / "cache") with pytest.raises(ValueError, match="use MapBiomas"): prodes.register_prodes_source(cube, "p", 2004, BBOX, tmp_path / "raw", files=files) + + +def test_codes_choose_the_local_code_of_each_category(server, tmp_path): + scheme = prodes.parse_legend(LEGEND) + assert scheme.lookup(2012)[91] == scheme.lookup(2012)[101] == 3 # default: water and non-forest together + t = scheme.lookup(2012, {"water": 4}) + assert t[91] == 4 and t[101] == 3 and t[100] == prodes.FOREST + merged = scheme.lookup(2012, {"forest": 7, "non_forest": 7, "water": 7}) + assert merged[100] == merged[101] == merged[91] == 7 and merged[7] == prodes.DEFORESTED + with pytest.raises(ValueError, match="unknown PRODES categories"): + prodes.resolve_codes({"cloud": 9}) + for bad in (0, 256, 2.5, True): + with pytest.raises(ValueError, match="1 to 255"): + prodes.resolve_codes({"water": bad}) + + cube = CubeClient(catalog=str(tmp_path / "c.db"), store=str(tmp_path / "s")) + files = prodes.download(server["url"], cache_dir=tmp_path / "cache") + plain = prodes.register_prodes_source(cube, "p_plain", 2016, BBOX, tmp_path / "raw", files=files) + apart = prodes.register_prodes_source(cube, "p_water", 2016, BBOX, tmp_path / "raw", files=files, + codes={"water": 4}) + with rasterio.open(plain.asset_url) as a, rasterio.open(apart.asset_url) as b: + va, vb = a.read(1), b.read(1) + assert 4 not in np.unique(va) and 4 in np.unique(vb) + assert np.array_equal(np.where(vb == 4, 3, vb), va) # only the water pixels change + prov = json.loads((tmp_path / "raw" / "p_water.provenance.json").read_text()) + assert prov["codes"] == {"forest": 1, "deforested": 2, "non_forest": 3, "water": 4, "nodata": 0} + + +def test_codes_in_a_pipeline_file(tmp_path): + from pydantic import ValidationError + + from disscube.pipeline.schema import ProdesSource + + assert ProdesSource(id="p", type="prodes", codes={"water": 4}).codes == {"water": 4} + with pytest.raises(ValidationError, match="unknown PRODES categories"): + ProdesSource(id="p", type="prodes", codes={"agua": 4}) diff --git a/tests/test_sqlite_connections.py b/tests/test_sqlite_connections.py new file mode 100644 index 0000000..66885d7 --- /dev/null +++ b/tests/test_sqlite_connections.py @@ -0,0 +1,30 @@ +"""The catalog must not leave SQLite connections open (they showed up as ResourceWarning and, at interpreter +exit, as `Error in sys.excepthook` after `disscube run`).""" +import gc +import sqlite3 +import warnings + +from disscube.catalog.sqlite_store import SqliteCatalogStore + + +def test_connection_is_closed_after_use(tmp_path): + store = SqliteCatalogStore(tmp_path / "catalog.db") + with store._get_connection() as conn: + conn.execute("SELECT 1") + try: + conn.execute("SELECT 1") + except sqlite3.ProgrammingError: + pass + else: + raise AssertionError("the connection is still open after leaving the with block") + + +def test_no_unclosed_database_warning(tmp_path): + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + store = SqliteCatalogStore(tmp_path / "catalog.db") + for _ in range(5): + with store._get_connection() as conn: + conn.execute("SELECT 1") + gc.collect() + assert not [w for w in caught if "unclosed database" in str(w.message)]