Feat/out of core sources

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

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

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}
 }
 ```

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:

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,

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:
         """