From 7c2b51a01ebe70c94e006c6db82108284252df6e Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Thu, 27 Aug 2026 15:53:31 -0300 Subject: [PATCH 01/17] update --- disscube/pipeline/aligner.py | 139 +++++++++++++++++++++++++++++++ tests/test_aligner_resampling.py | 7 +- 2 files changed, 145 insertions(+), 1 deletion(-) diff --git a/disscube/pipeline/aligner.py b/disscube/pipeline/aligner.py index 570d3ee..db77024 100644 --- a/disscube/pipeline/aligner.py +++ b/disscube/pipeline/aligner.py @@ -24,12 +24,14 @@ import logging +import rasterio import rioxarray # noqa: F401 — registers the .rio accessor import numpy as np import xarray as xr import geopandas as gpd from pyproj import CRS as ProjCRS, Transformer from rasterio.warp import Resampling +from rasterio.windows import Window from shapely.geometry import box from disscube.operators.base import OPERATOR_REGISTRY @@ -97,6 +99,23 @@ def _align_raster( band_map : dict[str, int] Optional ``{variable_name: 1-based band index}`` from the source. """ + # ── Identity fast path ────────────────────────────────────────── + # If the source already shares the target grid's CRS, resolution and + # pixel origin (no resampling is mathematically needed — every source + # pixel maps onto exactly one target pixel), skip GDAL's warp + # machinery entirely: `.rio.reproject()` always allocates a full + # destination array via a real resampling pass, and `_align_fine` + # additionally burns a whole extra reproject just to *estimate* the + # source resolution. A plain windowed `rasterio.read()` produces the + # identical result — this is the same real-vs-needed-work gap + # confirmed against real MapBiomas/ANADEM data (2026-08-27): those + # sources already share one 30m EPSG:5880 grid, so every derive() + # against them was reprojecting via GDAL only to reproduce numbers a + # windowed copy already had. + offset = self._identity_offset(url, grid) + if offset is not None: + return self._align_raster_identity(url, grid, variables, band_map, offset) + ds_src = rioxarray.open_rasterio(url) # Map of variable name -> aligned DataArray. A plain dict (not a # Dataset) is used because fine-aligned categorical arrays have a @@ -186,6 +205,126 @@ def _align_raster( return result + # ------------------------------------------------------------------ + # Identity fast path — no resampling needed, plain windowed read + # ------------------------------------------------------------------ + + def _identity_offset( + self, url: str, grid: GridSpec, offset_tolerance_px: float = 1e-3 + ) -> tuple[int, int] | None: + """ + Return ``(row_off, col_off)`` — the integer pixel offset of + ``grid``'s origin inside the source raster's own pixel grid — if, + and only if, the source can supply ``grid`` with a plain windowed + read: same CRS, same resolution, no rotation, origin aligned to a + whole pixel (within ``offset_tolerance_px``), and the requested + window fully covered by the source extent. + + Returns ``None`` for any other case (different CRS/resolution, + sub-pixel misalignment, or the target grid falling partially or + fully outside the source) — the caller then falls back to the + general reproject path unchanged. + """ + try: + with rasterio.open(url) as ds: + src_crs = ds.crs + transform = ds.transform + src_h, src_w = ds.height, ds.width + except Exception: + return None + + if src_crs is None: + return None + + try: + if not ProjCRS.from_user_input(src_crs).equals(ProjCRS.from_user_input(grid.crs)): + return None + except Exception: + return None + + px_w, rot1, ox, rot2, px_h, oy = ( + transform.a, transform.b, transform.c, transform.d, transform.e, transform.f, + ) + if abs(rot1) > 1e-9 or abs(rot2) > 1e-9: + return None + if abs(abs(px_w) - grid.resolution) > 1e-6 or abs(abs(px_h) - grid.resolution) > 1e-6: + return None + + minx, _miny, _maxx, maxy = grid.bbox + col_f = (minx - ox) / px_w + row_f = (oy - maxy) / abs(px_h) + col_off, row_off = round(col_f), round(row_f) + if abs(col_f - col_off) > offset_tolerance_px or abs(row_f - row_off) > offset_tolerance_px: + return None + + if ( + row_off < 0 or col_off < 0 + or row_off + grid.rows > src_h or col_off + grid.cols > src_w + ): + # Grid falls partially/fully outside the source's own extent — + # a plain read can't pad with nodata, so fall back to reproject + # (rio.reproject handles the padding via the target transform). + return None + + return row_off, col_off + + def _align_raster_identity( + self, + url: str, + grid: GridSpec, + variables: list[Variable], + band_map: dict[str, int], + offset: tuple[int, int], + ) -> dict[str, xr.DataArray]: + """ + Build one ``(grid.rows, grid.cols)`` DataArray per variable via a + plain windowed ``rasterio`` read — no reprojection. Only called when + ``_identity_offset`` has already confirmed the source needs none. + """ + row_off, col_off = offset + window = Window(col_off, row_off, grid.cols, grid.rows) + result: dict[str, xr.DataArray] = {} + + with rasterio.open(url) as ds: + n_bands = ds.count + for i, var in enumerate(variables): + if n_bands > 1: + if band_map and var.name in band_map: + band_idx = band_map[var.name] # rasterio bands are 1-based + if not (1 <= band_idx <= n_bands): + raise ValueError( + f"Band index {band_idx} for variable " + f"'{var.name}' is out of range; " + f"source has {n_bands} bands." + ) + elif i < n_bands: + band_idx = i + 1 + else: + raise ValueError( + f"No band available for variable '{var.name}' at " + f"index {i}; source has {n_bands} bands " + "and no band_map was provided." + ) + else: + band_idx = 1 + + arr = ds.read(band_idx, window=window) + nodata = ds.nodatavals[band_idx - 1] if ds.nodatavals else ds.nodata + + da = xr.DataArray(arr, dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs}) + da.rio.write_crs(grid.crs, inplace=True) + if nodata is not None: + da.attrs["_disscube_nodata"] = nodata + da.rio.write_nodata(nodata, inplace=True) + + result[var.name] = da + log.debug( + "identity-aligned '%s' (windowed read, no reprojection; " + "window=%s)", var.name, (row_off, col_off, grid.cols, grid.rows), + ) + + return result + # ------------------------------------------------------------------ # Crop source to target grid extent (avoids reading/reprojecting the # whole source raster for a small target grid) diff --git a/tests/test_aligner_resampling.py b/tests/test_aligner_resampling.py index 5c23ca1..12b9ec7 100644 --- a/tests/test_aligner_resampling.py +++ b/tests/test_aligner_resampling.py @@ -192,8 +192,13 @@ def test_alignment_invariant_raises_on_shape_mismatch(tmp_path, monkeypatch): We force the failure by monkeypatching the .rio.reproject result to a wrong shape, isolating the invariant check itself. + + The source is deliberately at a different resolution (5m) than the grid + (10m) so GridAligner's identity fast path (same CRS/resolution/origin -> + plain windowed read, no reproject) does not intercept — this test only + means to exercise the general reproject path's shape invariant. """ - src = np.zeros((10, 10), dtype=np.float32) + src = np.zeros((20, 20), dtype=np.float32) # 100/20 = 5m pixels, vs grid's 10m url = _write_raster(tmp_path / "z.tif", src) grid = _grid(resolution=10) # expects 10x10 From 4f304594f38f036ee5d25c6b71794546c4a49759 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Thu, 27 Aug 2026 16:31:03 -0300 Subject: [PATCH 02/17] feat: add generic reclassify operator (value lookup table) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds a domain-agnostic per-pixel value-remapping operator, filling the gap between the zonal operators (which aggregate several source pixels into one target cell) and the need to simply remap class codes 1:1. The mechanism (apply a table) lives in disscube; the table itself is data supplied by the caller via Variable.mapping, so land-cover, soil-class and land-tenure regroupings all reuse this operator without disscube knowing anything about any of those domains. - Variable.mapping / Derivation.mapping carry {source_value: target_value}. Being a plain field, it is folded into spec_hash() automatically, so two derivations with different tables are always distinct products. - Operator.requires_mapping mirrors the existing requires_class_code fail-fast contract; Derivation validates it at construction time. - Lookup is vectorized via np.searchsorted over sorted keys rather than a dense 0..max array, so a sparse table with large codes costs the same as a dense one. - Unmapped values, source nodata and non-finite input all become NaN — the same "no valid result here" convention the zonal operators already use. - No purity coordinates are attached: with one source pixel per target cell there is no sub-cell composition to summarise, so coverage/dominance would be a constant carrying no information at the cost of two extra full-size arrays. Verified against real data: reproducing the BR-MANGUE papel_dominio_2024 classification through this operator (from band 3 of br_mangue_base_v2.vrt) matches the existing GDAL pipeline's output pixel-for-pixel over a dense 2048x2048 window containing all six classes, including the exact nodata count. Tests: 17 new, full suite 95 passing. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- disscube/derivation.py | 11 ++ disscube/models/variable.py | 6 ++ disscube/operators/__init__.py | 2 +- disscube/operators/base.py | 5 + disscube/operators/reclassify.py | 107 ++++++++++++++++++ tests/test_reclassify.py | 179 +++++++++++++++++++++++++++++++ 6 files changed, 309 insertions(+), 1 deletion(-) create mode 100644 disscube/operators/reclassify.py create mode 100644 tests/test_reclassify.py diff --git a/disscube/derivation.py b/disscube/derivation.py index 53866cc..b4cb32a 100644 --- a/disscube/derivation.py +++ b/disscube/derivation.py @@ -42,6 +42,11 @@ class Derivation(BaseModel): class_code : int | None Target class code used by class-aware operators (e.g. ``"percentage"``). Required when the operator demands it; reserved but optional otherwise. + mapping : dict[int, int] | None + Source-value -> target-value lookup table used by value-remapping + operators (e.g. ``"reclassify"``). Required when the operator demands + it. Folded into ``spec_hash()`` via ``Variable``, so changing the + table always yields a distinct product. role : str Semantic role of the variable (e.g. ``"driver"``, ``"state"``). Defaults to ``"driver"``. @@ -66,6 +71,7 @@ class Derivation(BaseModel): source_id: str operator: str class_code: int | None = None + mapping: dict[int, int] | None = None role: str = "driver" valid_from: str | None = None valid_until: str | None = None @@ -85,6 +91,10 @@ def _validate_operator(self) -> "Derivation": raise ValueError( f"Operator {self.operator!r} requires class_code to be set." ) + if meta.requires_mapping and not self.mapping: + raise ValueError( + f"Operator {self.operator!r} requires mapping to be set." + ) return self # ── Conversion helpers ──────────────────────────────────────────────────── @@ -103,6 +113,7 @@ def to_variable(self) -> Variable: name=self.target, operator=self.operator, class_code=self.class_code, + mapping=self.mapping, ) def to_spatial_derivation(self, grid_id: str) -> SpatialDerivation: diff --git a/disscube/models/variable.py b/disscube/models/variable.py index 67ec7d4..0130412 100644 --- a/disscube/models/variable.py +++ b/disscube/models/variable.py @@ -10,6 +10,12 @@ class Variable(BaseModel): operator: str class_code: int | None = None + # Explicit source-value -> target-value lookup table, used by value-remapping + # operators (e.g. "reclassify"). Kept as plain data (not a callable) so it + # travels through the catalog and is folded into spec_hash() — two + # derivations using different tables are always distinct products. + mapping: dict[int, int] | None = None + class SpatialSource(BaseModel): id: str diff --git a/disscube/operators/__init__.py b/disscube/operators/__init__.py index 357a01b..ea15afc 100644 --- a/disscube/operators/__init__.py +++ b/disscube/operators/__init__.py @@ -7,7 +7,7 @@ """ # Import submodules to trigger auto-registration of all operator classes. -from . import zonal, proximity # noqa: F401 +from . import zonal, proximity, reclassify # noqa: F401 from .zonal import ZonalAggregator from .proximity import ProximityAggregator diff --git a/disscube/operators/base.py b/disscube/operators/base.py index b7acbd1..9c8fc2f 100644 --- a/disscube/operators/base.py +++ b/disscube/operators/base.py @@ -52,6 +52,11 @@ class body is executed (via ``__init_subclass__``), so there is no name: ClassVar[str] requires_class_code: ClassVar[bool] = False + # When True, ``Derivation`` enforces that ``mapping`` (a source-value -> + # target-value lookup table) is set at construction time, the same + # fail-fast contract ``requires_class_code`` provides for class-aware + # operators. + requires_mapping: ClassVar[bool] = False _resampling: ClassVar[Resampling] = Resampling.nearest # When True, GridAligner must NOT pre-aggregate the band with this diff --git a/disscube/operators/reclassify.py b/disscube/operators/reclassify.py new file mode 100644 index 0000000..41ca328 --- /dev/null +++ b/disscube/operators/reclassify.py @@ -0,0 +1,107 @@ +""" +Value-remapping operators — per-pixel lookup, no spatial aggregation. + +Unlike the zonal operators, nothing here combines several source pixels into +one target cell: each pixel's new value depends only on its own old value, +via an explicit lookup table carried on ``Variable.mapping``. + +This is deliberately domain-agnostic. The *mechanism* (apply a table) lives +here; the *table* (which source code means what) is data supplied by the +caller, so a land-cover reclassification, a soil-class grouping and a +land-tenure regrouping all use this same operator with different tables and +no changes to disscube. +""" + +from __future__ import annotations + +import numpy as np +import xarray as xr +from rasterio.warp import Resampling + +from disscube.operators.base import Operator +from disscube.models.variable import Variable +from disscube.models.grid import GridSpec + + +def _lookup(arr: np.ndarray, mapping: dict[int, int], nodata: float | None) -> np.ndarray: + """ + Apply ``mapping`` to ``arr``, returning float64 with NaN where a value is + nodata or absent from the table. + + Vectorized via ``np.searchsorted`` over the sorted key array rather than a + dense 0..max lookup array, so a sparse table with large codes (e.g. + MapBiomas-style codes in the hundreds) costs the same as a dense one. + """ + keys = np.array(sorted(mapping), dtype=np.int64) + vals = np.array([mapping[int(k)] for k in keys], dtype=np.float64) + + finite = np.isfinite(arr) + valid = finite.copy() + if nodata is not None and np.isfinite(nodata): + valid &= arr != nodata + + # Round to integer codes only where the value is usable; searchsorted needs + # a clean integer array, and non-finite entries would poison the cast. + codes = np.where(valid, arr, 0).astype(np.int64) + + idx = np.searchsorted(keys, codes) + np.clip(idx, 0, len(keys) - 1, out=idx) + found = valid & (keys[idx] == codes) + + return np.where(found, vals[idx], np.nan) + + +class ReclassifyOperator(Operator): + """ + Remap source values to new values through an explicit lookup table. + + ``Variable.mapping`` supplies ``{source_value: target_value}``. Any pixel + whose value is source-nodata, non-finite, or simply absent from the table + becomes NaN — the same "no valid result here" convention the zonal + operators use, so downstream readers need no special case. Callers that + need a specific sentinel (255, -9999, …) convert NaN on write. + + Resampling is NEAREST: a class code must never be averaged. Note that + when the target grid is coarser than the source this samples one source + pixel per target cell rather than taking a majority — reclassify is a + per-pixel remap, not an aggregation. To reclassify *and* downsample, + derive the reclassified variable on a grid at source resolution, then + derive ``majority`` from it. + + No purity coordinates are attached (unlike the categorical zonal + operators): with one source pixel per target cell there is no sub-cell + composition to summarise, so coverage/dominance would be a constant 1.0 + carrying no information, at the cost of two extra full-size arrays. + """ + + name = "reclassify" + _resampling = Resampling.nearest + requires_mapping = True + + def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray: + if not isinstance(data, xr.DataArray): + raise TypeError( + f"'reclassify' requires a raster source, got {type(data).__name__}" + ) + if not var.mapping: + raise ValueError( + f"Operator 'reclassify' requires a non-empty mapping for " + f"variable {var.name!r}." + ) + + da = data.isel(band=0) if "band" in data.dims else data + da = da.transpose("y", "x") + + nodata = da.attrs.get("_disscube_nodata", None) + if nodata is None: + try: + nodata = da.rio.nodata + except Exception: + nodata = None + + arr = np.asarray(da.values, dtype=np.float64) + out = _lookup(arr, var.mapping, nodata) + + return xr.DataArray( + out, dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs} + ) diff --git a/tests/test_reclassify.py b/tests/test_reclassify.py new file mode 100644 index 0000000..69e3572 --- /dev/null +++ b/tests/test_reclassify.py @@ -0,0 +1,179 @@ +""" +Tests for the generic value-remapping operator (``reclassify``). + +Covers the operator itself (lookup semantics, nodata/unmapped handling, +sparse tables) and its integration with the declarative ``Derivation`` +front-end (fail-fast validation, spec_hash sensitivity to the table). +""" + +import numpy as np +import pytest +import xarray as xr + +from disscube.derivation import Derivation +from disscube.models import GridSpec, Variable +from disscube.operators.base import OPERATOR_REGISTRY + +CRS = "EPSG:31982" + + +def _grid(rows=2, cols=3, resolution=10): + return GridSpec( + id="G1", type="local", crs=CRS, resolution=resolution, + bbox=[0, 0, cols * resolution, rows * resolution], + ) + + +def _da(array, grid, nodata=None): + da = xr.DataArray( + np.asarray(array, dtype=np.float64), + dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs}, + ) + if nodata is not None: + da.attrs["_disscube_nodata"] = nodata + return da + + +def _compute(array, mapping, grid=None, nodata=None): + grid = grid or _grid(*np.shape(array)) + op = OPERATOR_REGISTRY["reclassify"]() + var = Variable(name="v", operator="reclassify", mapping=mapping) + return op.compute(_da(array, grid, nodata), var, grid).values + + +# ── Lookup semantics ───────────────────────────────────────────────────────── + +def test_maps_every_value_through_the_table(): + src = [[1, 2, 3], [4, 5, 6]] + mapping = {1: 10, 2: 20, 3: 30, 4: 40, 5: 50, 6: 60} + out = _compute(src, mapping) + assert np.array_equal(out, np.array([[10, 20, 30], [40, 50, 60]], dtype=np.float64)) + + +def test_many_to_one_grouping(): + """Several source codes collapsing onto one target code is the common case.""" + src = [[2, 3, 4], [5, 6, 1]] + mapping = {1: 1, 2: 2, 3: 2, 4: 2, 5: 2, 6: 2} + out = _compute(src, mapping) + assert np.array_equal(out, np.array([[2, 2, 2], [2, 2, 1]], dtype=np.float64)) + + +def test_unmapped_values_become_nan(): + src = [[1, 99, 2]] + out = _compute(src, {1: 10, 2: 20}, grid=_grid(rows=1, cols=3)) + assert out[0, 0] == 10 + assert np.isnan(out[0, 1]) + assert out[0, 2] == 20 + + +def test_nodata_becomes_nan_even_when_present_in_the_table(): + """An explicit nodata sentinel wins over a table entry for the same value.""" + src = [[1, -9999, 2]] + mapping = {1: 10, 2: 20, -9999: 77} + out = _compute(src, mapping, grid=_grid(rows=1, cols=3), nodata=-9999) + assert out[0, 0] == 10 + assert np.isnan(out[0, 1]) + assert out[0, 2] == 20 + + +def test_nan_input_stays_nan(): + src = [[1.0, np.nan, 2.0]] + out = _compute(src, {1: 10, 2: 20}, grid=_grid(rows=1, cols=3)) + assert np.isnan(out[0, 1]) + + +def test_sparse_table_with_large_codes(): + """Codes far apart must not require a dense 0..max lookup array.""" + src = [[3, 500, 33000]] + mapping = {3: 1, 500: 2, 33000: 3} + out = _compute(src, mapping, grid=_grid(rows=1, cols=3)) + assert np.array_equal(out, np.array([[1, 2, 3]], dtype=np.float64)) + + +def test_zero_is_a_mappable_code_not_treated_as_missing(): + src = [[0, 1]] + out = _compute(src, {0: 7, 1: 8}, grid=_grid(rows=1, cols=2)) + assert np.array_equal(out, np.array([[7, 8]], dtype=np.float64)) + + +def test_negative_codes_are_mappable(): + src = [[-3, 2]] + out = _compute(src, {-3: 1, 2: 5}, grid=_grid(rows=1, cols=2)) + assert np.array_equal(out, np.array([[1, 5]], dtype=np.float64)) + + +def test_output_is_grid_shaped_with_grid_coords(): + grid = _grid(rows=2, cols=3) + op = OPERATOR_REGISTRY["reclassify"]() + var = Variable(name="v", operator="reclassify", mapping={1: 1}) + result = op.compute(_da(np.ones((2, 3)), grid), var, grid) + assert result.dims == ("y", "x") + assert result.shape == (grid.rows, grid.cols) + assert np.array_equal(result.coords["x"].values, grid.xs) + assert np.array_equal(result.coords["y"].values, grid.ys) + + +def test_band_dimension_is_collapsed(): + grid = _grid(rows=1, cols=2) + da = xr.DataArray( + np.array([[[1.0, 2.0]]]), + dims=("band", "y", "x"), + coords={"band": [1], "y": grid.ys, "x": grid.xs}, + ) + op = OPERATOR_REGISTRY["reclassify"]() + var = Variable(name="v", operator="reclassify", mapping={1: 9, 2: 8}) + out = op.compute(da, var, grid).values + assert np.array_equal(out, np.array([[9, 8]], dtype=np.float64)) + + +# ── Guard rails ────────────────────────────────────────────────────────────── + +def test_empty_mapping_raises(): + grid = _grid(rows=1, cols=2) + op = OPERATOR_REGISTRY["reclassify"]() + var = Variable(name="v", operator="reclassify", mapping=None) + with pytest.raises(ValueError, match="requires a non-empty mapping"): + op.compute(_da([[1, 2]], grid), var, grid) + + +def test_vector_source_raises(): + grid = _grid(rows=1, cols=2) + op = OPERATOR_REGISTRY["reclassify"]() + var = Variable(name="v", operator="reclassify", mapping={1: 1}) + with pytest.raises(TypeError, match="requires a raster source"): + op.compute("not a raster", var, grid) + + +# ── Declarative front-end ──────────────────────────────────────────────────── + +def test_derivation_requires_mapping(): + with pytest.raises(ValueError, match="requires mapping"): + Derivation(target="papel", source_id="s1", operator="reclassify") + + +def test_derivation_with_mapping_ok(): + d = Derivation( + target="papel", source_id="s1", operator="reclassify", mapping={1: 1, 2: 2}, + ) + assert d.to_variable().mapping == {1: 1, 2: 2} + + +def test_mapping_changes_spec_hash(): + """Two derivations differing only in the table must be distinct products.""" + a = Derivation(target="p", source_id="s", operator="reclassify", mapping={1: 1}) + b = Derivation(target="p", source_id="s", operator="reclassify", mapping={1: 2}) + assert a.spec_hash() != b.spec_hash() + + +def test_same_mapping_same_spec_hash(): + a = Derivation(target="p", source_id="s", operator="reclassify", mapping={1: 1, 2: 5}) + b = Derivation(target="p", source_id="s", operator="reclassify", mapping={2: 5, 1: 1}) + assert a.spec_hash() == b.spec_hash() + + +def test_string_keys_are_coerced_to_int(): + """Tables loaded from JSON/TOML arrive with string keys.""" + d = Derivation( + target="p", source_id="s", operator="reclassify", mapping={"1": 10, "2": 20}, + ) + assert d.to_variable().mapping == {1: 10, 2: 20} From be4aa510c6dfe96f6992fc4b1a5578d644a8d244 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Thu, 27 Aug 2026 16:49:33 -0300 Subject: [PATCH 03/17] docs: add BR-MANGUE structural domain example (reclassify + tiling) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds a standalone case study showing the full path from raw tiles to a categorical domain raster, exercising the new generic `reclassify` operator together with disscube's existing tile-based derivation. The example is also the clearest demonstration of where the boundary between the cube and a project sits: nothing mangrove-specific lives in disscube — `reclassify` just applies a table that arrives as data — while the papeis_estados table and the two-band papel/elevation coupling stay in the example, since Variable derives one band per variable. Why tile by tile: a single derive() over the full 18352x21350 extent (391M cells) allocates tens of GB. derive(tile_id=...) keeps memory bounded; the whole extent completes in ~80s on a modest machine. Why the example writes the GeoTIFF itself: tools/zarr_to_tif.py converts one whole Zarr at a time, but reassembling N tiles into a single multi-band GeoTIFF needs windowed writes. Verified against the reference implementation it replaces (br_mangue's GDAL pipeline): identical output over all 391M cells in both bands, same CRS, transform and band descriptions. geomosaic is imported with a clear install hint — it is not a disscube dependency, and examples are not part of the installed package. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- examples/README.md | 17 +- .../brmangue_dominio/01_dominio_estrutural.py | 208 ++++++++++++++++++ 2 files changed, 224 insertions(+), 1 deletion(-) create mode 100644 examples/case_studies/brmangue_dominio/01_dominio_estrutural.py diff --git a/examples/README.md b/examples/README.md index 30ea798..dc7dbcb 100644 --- a/examples/README.md +++ b/examples/README.md @@ -20,7 +20,22 @@ Dois estudos sobre a mesma área geográfica e grade. - `python examples/case_studies/maranhao/02_brmangue_derive.py` — deriva uso, alt, solo para o modelo BR-MANGUE. - `python examples/case_studies/maranhao/03_brmangue_simulate.py` — executa BrmangueRasterExecutor. -### 4. Estudo de caso: Acre (AC/5km) +### 4. Estudo de caso: BR-MANGUE domínio estrutural (nacional, 30 m) +Standalone — não depende dos scripts de `setup/`; registra a própria grade e fontes. +- `python examples/case_studies/brmangue_dominio/01_dominio_estrutural.py` — reclassifica + a legenda de estados num raster de papéis, tile a tile, e escreve o GeoTIFF final. + +Demonstra três coisas de uso geral: +- o operador genérico **`reclassify`** (tabela `{valor_origem: valor_destino}` como dado); +- **`derive(tile_id=...)`** para processar uma extensão grande (391M células) com + memória limitada — ver `docs/architecture/tiling.md`; +- a fronteira entre o que é do cubo e o que é do projeto: a tabela de papéis e o + cruzamento entre duas bandas ficam no exemplo, não no pacote. + +Requer `geomosaic` (`pip install geomosaic`) e a variável `BRMANGUE_ENTRADA` +apontando para o diretório com os tiles ANADEM v2. + +### 5. Estudo de caso: Acre (AC/5km) - `python examples/drivers/02_acre_5km.py` — drivers regionais Acre 5 km. - `python examples/case_studies/lucc_acre/01_derive.py` — atributos de uso do solo de fonte vetorial. - `python examples/case_studies/lucc_acre/02_simulate.py` — executa LUCCRasterExecutor. diff --git a/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py b/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py new file mode 100644 index 0000000..b4a583c --- /dev/null +++ b/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py @@ -0,0 +1,208 @@ +""" +examples/case_studies/brmangue_dominio/01_dominio_estrutural.py + +BR-MANGUE — domínio estrutural nacional (30 m, EPSG:5880). + +Mostra o caminho completo para transformar tiles brutos num raster +categórico de domínio, usando o operador genérico `reclassify`: + + geomosaic tiles ANADEM -> VRT (uma banda por chamada) + disscube brm_state_2024 -> papel, tile a tile (derive(tile_id=...)) + rasterio escrita janelada do GeoTIFF de 2 bandas + +O que é genérico e o que é do projeto +------------------------------------- +Do disscube, nada aqui é específico de manguezal: `reclassify` só aplica +uma tabela `{valor_origem: valor_destino}` que chega como dado. O que é +BR-MANGUE são as ~10 linhas de PAPEIS_ESTADOS e o acoplamento +papel/elevação no fim do loop — ambos vivem neste exemplo, não no pacote. + +Por que tile a tile +------------------- +Um derive() único na extensão completa (18352x21350 = 391M células) +alocaria dezenas de GB e derruba máquinas modestas. `derive(tile_id=...)` +processa um recorte por vez com memória limitada — ver +docs/architecture/tiling.md. + +Por que a escrita é feita aqui +------------------------------ +`tools/zarr_to_tif.py` converte UM Zarr inteiro de uma vez; para remontar +N tiles num único GeoTIFF multibanda é preciso escrever janela a janela, +o que é feito no loop abaixo. + +Pré-requisitos +-------------- + - pip install geomosaic (não é dependência do disscube) + - tiles ANADEM v2 em $BRMANGUE_ENTRADA/anadem_v2/ + (3 bandas: elevation_m, lulc_mb_2024, brm_state_2024) + +Usage: + export BRMANGUE_ENTRADA=/caminho/para/pymangue/dados/entrada + python examples/case_studies/brmangue_dominio/01_dominio_estrutural.py +""" + +import os +import time +from pathlib import Path + +import numpy as np +import rasterio +from rasterio.windows import Window + +from disscube.client import CubeClient +from disscube.models import GridSpec, SpatialSource, SpatialDerivation, Variable + +try: + from geomosaic.core import build_mosaic_contract, write_vrt +except ImportError as exc: # pragma: no cover — exemplo, não faz parte do pacote + raise SystemExit( + "Este exemplo precisa do geomosaic (não é dependência do disscube):\n" + " pip install geomosaic\n" + "ou, a partir do repositório irmão:\n" + " pip install -e ../geomosaic" + ) from exc + + +# ── Configuração ────────────────────────────────────────────────────────────── +ENTRADA = Path(os.environ.get("BRMANGUE_ENTRADA", "./dados/entrada")) +TILES_DIR = ENTRADA / "anadem_v2" +TRABALHO = Path("./data/brmangue_dominio") +SAIDA = TRABALHO / "dominio_estrutural_2024.tif" + +GRID_ID = "brmangue/30m" +TILE = 4096 # lado do recorte processado por vez +BANDA_ESTADO = 3 # brm_state_2024 +BANDA_ELEV = 1 # elevation_m +NODATA_SAIDA = 255 +NODATA_ELEV = -9999.0 + +# ── Específico do BR-MANGUE — dado do projeto, não do disscube ──────────────── +# Agrupa a legenda de estados 0–12 nos seis papéis do modelo. +PAPEIS_ESTADOS = { + "excluidos": ([0, 12], 0), + "mangue_inicial": ([1], 1), + "acomodacao_natural": ([2, 3, 4, 5, 6], 2), + "uso_manejado_candidato": ([7], 3), + "restricao_ou_barreira": ([8, 9, 10], 4), + "agua_aberta_nao_validada": ([11], 5), +} +MAPPING = {c: papel for _, (codes, papel) in PAPEIS_ESTADOS.items() for c in codes} + + +def main() -> None: + if not TILES_DIR.is_dir(): + raise RuntimeError( + f"Diretório de tiles não encontrado: {TILES_DIR}\n" + "Defina BRMANGUE_ENTRADA apontando para .../pymangue/dados/entrada" + ) + tiles = sorted(TILES_DIR.glob("*.tif")) + if not tiles: + raise RuntimeError(f"Nenhum .tif em {TILES_DIR}") + + t0 = time.perf_counter() + TRABALHO.mkdir(parents=True, exist_ok=True) + + # ── 1. geomosaic: contrato + um VRT por banda ──────────────────────────── + print(f"\n[1/3] geomosaic: {len(tiles)} tiles -> contrato") + contract = build_mosaic_contract([str(t) for t in tiles]) + altura, largura = contract.mosaic_height, contract.mosaic_width + a, _b, ox, _d, e, oy = contract.mosaic_transform + resolution = a + print(f" grade mestra: {altura}x{largura} @ {contract.crs}") + + vrt_estado = write_vrt(contract, str(TRABALHO / "estado.vrt"), band=BANDA_ESTADO) + vrt_elev = write_vrt(contract, str(TRABALHO / "elevacao.vrt"), band=BANDA_ELEV) + + # ── 2. disscube: grade, fontes e tiles ─────────────────────────────────── + print("[2/3] disscube: registrando grade, fontes e tiles") + cube = CubeClient( + catalog=str(TRABALHO / "catalog.db"), store=str(TRABALHO / "store") + ) + cube.register_grid(GridSpec( + id=GRID_ID, type="reference", crs=str(contract.crs), resolution=resolution, + bbox=[ox, oy - altura * abs(e), ox + largura * resolution, oy], + )) + cube.register_spatial_source(SpatialSource( + id="estado_2024", name="brm_state_2024", format="raster", + asset_url=str(vrt_estado), crs=str(contract.crs), + )) + cube.register_spatial_source(SpatialSource( + id="elevacao", name="elevation_m", format="raster", + asset_url=str(vrt_elev), crs=str(contract.crs), + )) + + # Cada recorte vira um SpatialSource {grid_id}_{tile_id} carregando só o + # bbox — é assim que derive(tile_id=...) descobre a janela a processar. + janelas = [] + for r0 in range(0, altura, TILE): + for c0 in range(0, largura, TILE): + h, w = min(TILE, altura - r0), min(TILE, largura - c0) + tile_id = f"R{r0:05d}C{c0:05d}" + minx, maxy = ox + c0 * resolution, oy - r0 * abs(e) + cube.register_spatial_source(SpatialSource( + id=f"{GRID_ID}_{tile_id}", name=tile_id, format="raster", + asset_url=str(vrt_estado), crs=str(contract.crs), + bbox=[minx, maxy - h * abs(e), minx + w * resolution, maxy], + )) + janelas.append((tile_id, r0, c0, h, w)) + print(f" {len(janelas)} tiles de {TILE}x{TILE}") + + deriv_papel = SpatialDerivation( + source_id="estado_2024", grid_id=GRID_ID, role="state", + variables=[Variable(name="papel", operator="reclassify", mapping=MAPPING)], + ) + deriv_elev = SpatialDerivation( + source_id="elevacao", grid_id=GRID_ID, role="driver", + variables=[Variable(name="elevacao", operator="mean")], + ) + + # ── 3. derive por tile + escrita janelada ──────────────────────────────── + print(f"[3/3] derive() por tile -> {SAIDA}") + perfil = dict( + driver="GTiff", height=altura, width=largura, count=2, dtype="uint8", + nodata=NODATA_SAIDA, crs=contract.crs, + transform=rasterio.transform.Affine(a, 0.0, ox, 0.0, e, oy), + tiled=True, blockxsize=512, blockysize=512, + compress="deflate", zlevel=6, bigtiff="IF_SAFER", sparse_ok=True, + ) + + with rasterio.open(SAIDA, "w", **perfil) as dst: + dst.set_band_description(1, "papel_dominio_2024") + dst.set_band_description(2, "elevacao_valida") + + for i, (tile_id, r0, c0, h, w) in enumerate(janelas, start=1): + cube.derive(deriv_papel, tile_id=tile_id) + cube.derive(deriv_elev, tile_id=tile_id) + + papel = cube.load("papel", tile_id=tile_id).values + elev = cube.load("elevacao", tile_id=tile_id).values + + # ── Regra do projeto: acoplamento papel/elevação ────────────── + # `elevacao_valida` só é 0/1 onde o ESTADO é válido; onde o estado + # é nodata a banda 2 também é nodata. O modelo Variable do + # disscube deriva uma banda por variável, então este cruzamento + # entre duas bandas é responsabilidade de quem chama. + estado_valido = ~np.isnan(papel) + elev_valida = np.isfinite(elev) & (elev != NODATA_ELEV) + + janela = Window(c0, r0, w, h) + dst.write( + np.where(estado_valido, papel, NODATA_SAIDA).astype(np.uint8)[:h, :w], + 1, window=janela, + ) + dst.write( + np.where( + estado_valido, elev_valida.astype(np.uint8), NODATA_SAIDA + ).astype(np.uint8)[:h, :w], + 2, window=janela, + ) + + if i % 5 == 0 or i == len(janelas): + print(f" {i}/{len(janelas)} tiles ({100 * i / len(janelas):.0f}%)") + + print(f"\n=== domínio estrutural gerado em {time.perf_counter() - t0:.1f}s ===") + print(f" {SAIDA} ({altura}x{largura}, 2 bandas uint8)") + + +if __name__ == "__main__": + main() From 09f1c3391e2337ce752f8ca776681c22b90108b4 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Thu, 27 Aug 2026 17:25:07 -0300 Subject: [PATCH 04/17] fix: resolve BDC tile sources in derive(tile_id=...) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The documented BDC workflow did not work. bdc_importer registers tiles as `BDC_{LEVEL}_{tile}` (e.g. BDC_SM_009002), but derive() looked up only `{grid_id}_{tile_id}` (e.g. BR/5km_009002), so the exact example in docs/guides/bdc.md failed with: ValueError: Tile 009002 with valid bbox not found for grid BR/5km BDC tiles are grid-independent — the same envelope serves every grid sharing the BDC CRS — so they are registered once, not per grid. Tile lookup now tries, in order: the grid-scoped id, the tile_id as a fully qualified id, then each BDC level. Ambiguity is an error, not a precedence rule: tile ids are NOT unique across levels. In the real V2 grids 189 ids exist in both SM and MD, and they cover different areas (verified: BDC_SM_005004 is at 3152000,11425600 while BDC_MD_005004 is at 3680000,10897600). Resolving a bare ambiguous id by trying SM first would silently derive the wrong extent, so it raises and names the candidates, pointing at the fully-qualified form instead. Docs updated accordingly: the loop in guides/bdc.md now passes the full tile source id rather than splitting off the suffix, with a warning about cross-level collisions, and the lookup order is documented in both guides/bdc.md and architecture/tiling.md. Tests: 12 new covering both conventions, precedence, ambiguity and the not-found message. Full suite 107 passing. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- disscube/client/cube_client.py | 70 +++++++++++++++++++- docs/architecture/tiling.md | 5 +- docs/guides/bdc.md | 24 ++++++- tests/test_tile_resolution.py | 115 +++++++++++++++++++++++++++++++++ 4 files changed, 208 insertions(+), 6 deletions(-) create mode 100644 tests/test_tile_resolution.py diff --git a/disscube/client/cube_client.py b/disscube/client/cube_client.py index 04fee9e..64f5f65 100644 --- a/disscube/client/cube_client.py +++ b/disscube/client/cube_client.py @@ -64,9 +64,7 @@ def derive(self, derivation: SpatialDerivation, tile_id: Optional[str] = None) - raise ValueError(f"Grid not found: {derivation.grid_id}") if tile_id: - tile_source = self.catalog.get_spatial_source(f"{grid.id}_{tile_id}") - if not tile_source or not tile_source.bbox: - raise ValueError(f"Tile {tile_id} with valid bbox not found for grid {grid.id}") + tile_source = self._resolve_tile_source(grid.id, tile_id) grid = GridSpec( id=grid.id, @@ -94,6 +92,72 @@ def derive(self, derivation: SpatialDerivation, tile_id: Optional[str] = None) - if d.spec_hash == spec_hash ] + # BDC tile levels, coarsest-last. Order is only used for reporting; an + # ambiguous id is always an error, never resolved by precedence. + _BDC_TILE_LEVELS = ("SM", "MD", "LG") + + def _resolve_tile_source(self, grid_id: str, tile_id: str) -> SpatialSource: + """ + Find the ``SpatialSource`` whose ``bbox`` defines ``tile_id``. + + Two registration conventions are supported, tried in this order: + + 1. ``{grid_id}_{tile_id}`` — a tile mesh defined for one specific + grid (see ``docs/architecture/tiling.md``). + 2. ``BDC_{LEVEL}_{tile_id}`` — the national BDC tile grids, as + registered by ``disscube.utils.bdc_importer``. These are + grid-independent: the same tile envelope serves every grid that + shares the BDC CRS, so they are not duplicated per grid. + + A fully-qualified id (``"BDC_SM_027005"``) is also accepted directly, + which is the way to disambiguate a bare id that exists at more than + one BDC level. + + Raises + ------ + ValueError + If nothing matches, if the match has no ``bbox``, or if a bare + tile id exists at several BDC levels. The levels cover different + areas (SM ~1.5°, MD ~3°, LG ~6°) and share many ids, so guessing + a level would silently derive the wrong extent. + """ + scoped = self.catalog.get_spatial_source(f"{grid_id}_{tile_id}") + if scoped is not None: + if not scoped.bbox: + raise ValueError( + f"Tile source {scoped.id!r} has no bbox; a tile source must " + "carry the bbox of the partition to process." + ) + return scoped + + qualified = self.catalog.get_spatial_source(tile_id) + if qualified is not None and qualified.bbox: + return qualified + + matches = [] + for level in self._BDC_TILE_LEVELS: + candidate = self.catalog.get_spatial_source(f"BDC_{level}_{tile_id}") + if candidate is not None and candidate.bbox: + matches.append((level, candidate)) + + if len(matches) == 1: + return matches[0][1] + + if len(matches) > 1: + found = ", ".join(f"BDC_{lvl}_{tile_id}" for lvl, _ in matches) + raise ValueError( + f"Tile id {tile_id!r} is ambiguous: it exists at several BDC " + f"levels ({found}), which cover different areas. Pass the " + f"fully-qualified id instead, e.g. tile_id='BDC_" + f"{matches[0][0]}_{tile_id}'." + ) + + raise ValueError( + f"Tile {tile_id!r} with valid bbox not found for grid {grid_id!r}. " + f"Looked for {grid_id}_{tile_id}, {tile_id}, and BDC_" + f"{{{','.join(self._BDC_TILE_LEVELS)}}}_{tile_id}." + ) + def derive_declarative( self, derivation: "Derivation", # noqa: F821 — imported lazily to avoid circular refs diff --git a/docs/architecture/tiling.md b/docs/architecture/tiling.md index 857fae3..60b3997 100644 --- a/docs/architecture/tiling.md +++ b/docs/architecture/tiling.md @@ -23,7 +23,10 @@ cube.derive(derivation, tile_id="009002") Internamente: 1. Busca `GridSpec` da master grid. -2. Busca `SpatialSource` com id `{grid_id}_{tile_id}` para obter o `bbox` do tile. +2. Busca o `SpatialSource` que carrega o `bbox` do tile, tentando nesta ordem: + `{grid_id}_{tile_id}` (malha própria da grade), o `tile_id` como id completo, + e `BDC_{SM,MD,LG}_{tile_id}` (grades nacionais BDC). Um id simples que exista + em mais de um nível BDC levanta `ValueError` — ver `docs/guides/bdc.md`. 3. Cria um `GridSpec` temporário: mesmos CRS e resolução, bbox restrito ao tile. 4. Executa o pipeline nessa grade temporária. 5. Salva em `data/derived/{grid_id}/009002/{spec_hash}/{var}.zarr`. diff --git a/docs/guides/bdc.md b/docs/guides/bdc.md index 7e8b14c..89020cd 100644 --- a/docs/guides/bdc.md +++ b/docs/guides/bdc.md @@ -46,12 +46,32 @@ cube.derive(derivation, tile_id="009002") # Tiles são registrados com IDs no formato BDC_SM_ (ex: BDC_SM_009002) tiles = [s for s in cube.catalog.list_spatial_sources() if s.id.startswith("BDC_SM_")] for tile_source in tiles: - tile_id = tile_source.id.split("_")[-1] - cube.derive(derivation, tile_id=tile_id) + cube.derive(derivation, tile_id=tile_source.id) # ID completo: BDC_SM_009002 ``` Cada tile é processado de forma independente e pode ser paralelizado. +!!! warning "IDs de tile não são únicos entre níveis" + O mesmo número identifica tiles diferentes em níveis diferentes — nas grades + V2, 189 IDs existem em SM **e** MD, cobrindo áreas distintas. Um ID simples + (`"009002"`) só é aceito quando existe em um único nível; se existir em mais + de um, `derive()` levanta `ValueError` em vez de escolher por precedência. + **Em loops, prefira sempre o ID completo** (`tile_source.id`), como acima. + +## Como o tile é localizado + +`derive(tile_id=...)` procura o `SpatialSource` que carrega o `bbox`, nesta ordem: + +1. `{grid_id}_{tile_id}` — malha de tiles definida para uma grade específica + (ver `docs/architecture/tiling.md`); +2. o `tile_id` tratado como ID completo — é assim que se desambigua + (`"BDC_SM_009002"`); +3. `BDC_{SM,MD,LG}_{tile_id}` — as grades nacionais BDC, registradas uma única + vez e compartilhadas por todas as grades no CRS do BDC. + +Os tiles BDC são independentes de grade: a mesma envoltória serve qualquer grade +que compartilhe o CRS, por isso não são duplicados por grade. + ## Carregar resultado tileado ```python diff --git a/tests/test_tile_resolution.py b/tests/test_tile_resolution.py new file mode 100644 index 0000000..24ad9e2 --- /dev/null +++ b/tests/test_tile_resolution.py @@ -0,0 +1,115 @@ +""" +Tests for tile-source resolution in ``CubeClient._resolve_tile_source``. + +Two registration conventions must both work: + * ``{grid_id}_{tile_id}`` — a tile mesh defined for one specific grid; + * ``BDC_{LEVEL}_{tile_id}`` — the national BDC grids, registered once and + shared by every grid in the BDC CRS. + +BDC tile ids are NOT unique across levels (in the real V2 grids 189 ids exist +in both SM and MD, covering different areas), so a bare ambiguous id must be +an error rather than resolved by precedence. +""" + +import tempfile +from pathlib import Path + +import pytest + +from disscube.client import CubeClient +from disscube.models import SpatialSource + +GRID_ID = "BR/5km" +CRS = "EPSG:4326" + + +@pytest.fixture +def cube(): + with tempfile.TemporaryDirectory() as d: + yield CubeClient(catalog=str(Path(d) / "c.db"), store=str(Path(d) / "store")) + + +def _tile(cube, source_id, bbox=(0.0, 0.0, 1.0, 1.0)): + cube.register_spatial_source(SpatialSource( + id=source_id, name=source_id, format="raster", + asset_url="planned", crs=CRS, bbox=list(bbox), + )) + + +# ── Grid-scoped convention ─────────────────────────────────────────────────── + +def test_grid_scoped_tile_is_found(cube): + _tile(cube, f"{GRID_ID}_T01", bbox=(1, 2, 3, 4)) + assert cube._resolve_tile_source(GRID_ID, "T01").bbox == [1, 2, 3, 4] + + +def test_grid_scoped_wins_over_bdc(cube): + """A mesh registered for this grid takes precedence over a BDC tile.""" + _tile(cube, f"{GRID_ID}_027005", bbox=(9, 9, 10, 10)) + _tile(cube, "BDC_SM_027005", bbox=(1, 1, 2, 2)) + assert cube._resolve_tile_source(GRID_ID, "027005").bbox == [9, 9, 10, 10] + + +def test_tile_source_without_bbox_raises(cube): + cube.register_spatial_source(SpatialSource( + id=f"{GRID_ID}_T01", name="T01", format="raster", + asset_url="planned", crs=CRS, # sem bbox + )) + with pytest.raises(ValueError, match="no bbox"): + cube._resolve_tile_source(GRID_ID, "T01") + + +# ── BDC convention ─────────────────────────────────────────────────────────── + +def test_bare_bdc_tile_id_resolves(cube): + """The workflow documented in docs/guides/bdc.md: grid BR/5km + bare id.""" + _tile(cube, "BDC_SM_027005", bbox=(5, 6, 7, 8)) + assert cube._resolve_tile_source(GRID_ID, "027005").bbox == [5, 6, 7, 8] + + +@pytest.mark.parametrize("level", ["SM", "MD", "LG"]) +def test_each_bdc_level_resolves(cube, level): + _tile(cube, f"BDC_{level}_000123", bbox=(1, 1, 2, 2)) + assert cube._resolve_tile_source(GRID_ID, "000123").id == f"BDC_{level}_000123" + + +def test_fully_qualified_bdc_id_resolves(cube): + _tile(cube, "BDC_MD_000123", bbox=(3, 3, 4, 4)) + assert cube._resolve_tile_source(GRID_ID, "BDC_MD_000123").bbox == [3, 3, 4, 4] + + +# ── Ambiguity must fail loudly ─────────────────────────────────────────────── + +def test_ambiguous_bare_id_raises(cube): + """Same id at two levels = different areas; guessing would be silently wrong.""" + _tile(cube, "BDC_SM_005004", bbox=(1, 1, 2, 2)) + _tile(cube, "BDC_MD_005004", bbox=(50, 50, 60, 60)) + with pytest.raises(ValueError, match="ambiguous"): + cube._resolve_tile_source(GRID_ID, "005004") + + +def test_ambiguity_error_lists_candidates_and_suggests_qualified_id(cube): + _tile(cube, "BDC_SM_005004") + _tile(cube, "BDC_LG_005004") + with pytest.raises(ValueError) as exc: + cube._resolve_tile_source(GRID_ID, "005004") + msg = str(exc.value) + assert "BDC_SM_005004" in msg and "BDC_LG_005004" in msg + assert "tile_id='BDC_SM_005004'" in msg + + +def test_qualified_id_escapes_ambiguity(cube): + _tile(cube, "BDC_SM_005004", bbox=(1, 1, 2, 2)) + _tile(cube, "BDC_MD_005004", bbox=(50, 50, 60, 60)) + assert cube._resolve_tile_source(GRID_ID, "BDC_SM_005004").bbox == [1, 1, 2, 2] + assert cube._resolve_tile_source(GRID_ID, "BDC_MD_005004").bbox == [50, 50, 60, 60] + + +# ── Not found ──────────────────────────────────────────────────────────────── + +def test_missing_tile_raises_with_searched_ids(cube): + with pytest.raises(ValueError) as exc: + cube._resolve_tile_source(GRID_ID, "999999") + msg = str(exc.value) + assert f"{GRID_ID}_999999" in msg + assert "BDC_" in msg From 6c16c6b7ec2bf4004733e53617655bd2ac15674b Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Thu, 27 Aug 2026 17:40:00 -0300 Subject: [PATCH 05/17] fix: don't read the whole source when it doesn't overlap the grid MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `_crop_to_grid` caught every exception from `clip_box` and fell back to returning the unclipped band. One of those exceptions is `NoDataInBounds`, raised precisely when the source and the target grid do not overlap — so "there is nothing here" was handled as "read everything", then reprojected into an all-nodata window. Measured on real data (BDC_SM tile 027005 against the BR-MANGUE ANADEM mosaic, EPSG:5880 -> BDC Albers): the fallback materialised the full 18352x21350 source, 31.6x the 3520x3520 target window, peaking at 2.42GB RSS to produce 50MB of nodata. With the fix the same tile stays at 0.18GB. This is a normal case, not an error: tile meshes are selected by envelope, so some tiles fall outside the data. `_crop_to_grid` now returns None for it and `_align_raster` builds the empty result directly, carrying the source nodata in `_disscube_nodata` so categorical operators mark every cell invalid exactly as they would for a reprojected empty window. Other exceptions keep the previous fallback: an unexpected failure should degrade to a correct if expensive read, never to wrong output. Tests: 7 new covering both the crop contract and the aligned result, including a guard that the shortcut does not swallow the normal path. Full suite 114 passing. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- disscube/pipeline/aligner.py | 68 ++++++++++++++++++-- tests/test_aligner_no_overlap.py | 106 +++++++++++++++++++++++++++++++ 2 files changed, 168 insertions(+), 6 deletions(-) create mode 100644 tests/test_aligner_no_overlap.py diff --git a/disscube/pipeline/aligner.py b/disscube/pipeline/aligner.py index db77024..0094b5d 100644 --- a/disscube/pipeline/aligner.py +++ b/disscube/pipeline/aligner.py @@ -32,6 +32,7 @@ from pyproj import CRS as ProjCRS, Transformer from rasterio.warp import Resampling from rasterio.windows import Window +from rioxarray.exceptions import NoDataInBounds from shapely.geometry import box from disscube.operators.base import OPERATOR_REGISTRY @@ -158,6 +159,20 @@ def _align_raster( # this into a windowed read. band = self._crop_to_grid(band, grid) + if band is None: + # Source and grid do not overlap at all. Reading the source to + # reproject it would produce an all-nodata result at the cost + # of materialising the whole raster — build the empty result + # directly instead. This is a normal case in tiled workflows, + # where the tile mesh is selected by envelope and some tiles + # fall outside the data. + result[var.name] = self._empty_for_grid(grid, ds_src) + log.debug( + "'%s': source does not overlap grid %r; emitting empty result", + var.name, grid.id, + ) + continue + # ── Per-operator resampling method ───────────────────────── op_cls = OPERATOR_REGISTRY.get(var.operator) needs_fine = bool(getattr(op_cls, "needs_fine_alignment", False)) @@ -330,7 +345,32 @@ def _align_raster_identity( # whole source raster for a small target grid) # ------------------------------------------------------------------ - def _crop_to_grid(self, band: xr.DataArray, grid: GridSpec, buffer_px: int = 4) -> xr.DataArray: + def _empty_for_grid(self, grid: GridSpec, ds_src: xr.DataArray) -> xr.DataArray: + """ + Build an all-nodata ``(grid.rows, grid.cols)`` array for a source that + does not overlap ``grid``. + + Carries the source nodata in ``_disscube_nodata`` so the categorical + operators mark every cell invalid, exactly as they would for a + reprojected but empty window. + """ + try: + nodata = ds_src.rio.nodata + except Exception: + nodata = None + fill = float(nodata) if nodata is not None else np.nan + + da = xr.DataArray( + np.full((grid.rows, grid.cols), fill, dtype=np.float64), + dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs}, + ) + if nodata is not None: + da.attrs["_disscube_nodata"] = nodata + return da + + def _crop_to_grid( + self, band: xr.DataArray, grid: GridSpec, buffer_px: int = 4 + ) -> xr.DataArray | None: """ Crop ``band`` (still in its native CRS) to the region overlapping ``grid``'s bbox, with a small buffer for resampling kernels, before @@ -340,10 +380,20 @@ def _crop_to_grid(self, band: xr.DataArray, grid: GridSpec, buffer_px: int = 4) (when it differs from ``grid.crs``) so the crop is correct even when the source raster is not already in the target CRS. - Falls back to the unclipped ``band`` (same behaviour as before this - fix) if the CRS is unknown or the crop fails for any reason — e.g. - the source doesn't actually overlap the grid, which the existing - downstream shape/coverage checks already surface clearly. + Returns + ------- + xr.DataArray | None + The cropped band, or ``None`` when the source and the grid do not + overlap at all. ``None`` is a distinct outcome from a failed crop: + reading the whole source only to reproject it into an all-nodata + window costs the full raster in memory (measured at 31x the target + window for a real BDC tile) and produces nothing, so the caller + builds the empty result directly instead. + + Falls back to the unclipped ``band`` only when the crop cannot be + computed for some other reason (unknown CRS, transform failure) — the + pre-crop behaviour, kept so an unexpected failure degrades to a correct + if expensive read rather than to wrong output. """ try: src_crs = band.rio.crs @@ -370,9 +420,15 @@ def _crop_to_grid(self, band: xr.DataArray, grid: GridSpec, buffer_px: int = 4) try: return band.rio.clip_box(minx, miny, maxx, maxy, auto_expand=True) + except NoDataInBounds: + log.debug( + "crop-to-grid: source does not overlap grid %r; " + "skipping the read entirely", grid.id, + ) + return None except Exception: log.debug( - "crop-to-grid: clip_box failed (grid may not overlap source), " + "crop-to-grid: clip_box failed for an unexpected reason, " "falling back to unclipped read", exc_info=True, ) return band diff --git a/tests/test_aligner_no_overlap.py b/tests/test_aligner_no_overlap.py new file mode 100644 index 0000000..2340302 --- /dev/null +++ b/tests/test_aligner_no_overlap.py @@ -0,0 +1,106 @@ +""" +Tests for the no-overlap path in ``GridAligner``. + +A tile mesh is normally selected by envelope, so some tiles legitimately fall +outside the data. Reading the whole source only to reproject it into an +all-nodata window costs the entire raster in memory and produces nothing — +measured at 31x the target window for a real BDC tile against a real mosaic. +``_crop_to_grid`` must report the no-overlap case instead of falling back to +an unclipped read, and the alignment must still yield a correctly shaped, +fully invalid result. +""" + +import numpy as np +import pytest +import rasterio +import xarray as xr +from rasterio.transform import from_bounds + +from disscube.models import GridSpec, SpatialSource, SpatialDerivation, Variable +from disscube.pipeline import PipelineContext +from disscube.pipeline.aligner import GridAligner + +CRS = "EPSG:31982" + + +def _write_raster(path, array, bbox, nodata=-9999.0): + rows, cols = array.shape + with rasterio.open( + path, "w", driver="GTiff", height=rows, width=cols, count=1, + dtype="float32", crs=CRS, transform=from_bounds(*bbox, cols, rows), + nodata=nodata, + ) as dst: + dst.write(array.astype("float32"), 1) + return str(path) + + +def _grid(bbox, resolution=10, gid="G1"): + return GridSpec(id=gid, type="local", crs=CRS, resolution=resolution, bbox=list(bbox)) + + +def _align(url, grid, operator="mean"): + aligner = GridAligner() + source = SpatialSource(id="S1", name="S1", format="raster", asset_url=url, crs=CRS) + variables = [Variable(name="v", operator=operator)] + ctx = PipelineContext( + source=source, grid=grid, + derivation=SpatialDerivation( + source_id="S1", grid_id=grid.id, role="test", variables=variables + ), + ) + return aligner.execute(ctx).data["v"] + + +# ── _crop_to_grid contract ─────────────────────────────────────────────────── + +def test_crop_returns_none_when_grid_is_outside_source(tmp_path): + import rioxarray + url = _write_raster(tmp_path / "s.tif", np.ones((10, 10)), bbox=(0, 0, 100, 100)) + band = rioxarray.open_rasterio(url).isel(band=0) + # Grid far away from the source extent. + assert GridAligner()._crop_to_grid(band, _grid((500, 500, 600, 600))) is None + + +def test_crop_returns_data_when_grid_overlaps(tmp_path): + import rioxarray + url = _write_raster(tmp_path / "s.tif", np.ones((10, 10)), bbox=(0, 0, 100, 100)) + band = rioxarray.open_rasterio(url).isel(band=0) + cropped = GridAligner()._crop_to_grid(band, _grid((20, 20, 50, 50))) + assert cropped is not None + assert cropped.size < band.size # genuinely cropped, not the whole raster + + +# ── Alignment result for a non-overlapping grid ────────────────────────────── + +@pytest.mark.parametrize("operator", ["mean", "majority"]) +def test_non_overlapping_grid_yields_empty_result_of_grid_shape(tmp_path, operator): + url = _write_raster(tmp_path / "s.tif", np.ones((10, 10)), bbox=(0, 0, 100, 100)) + grid = _grid((500, 500, 600, 600)) + aligned = _align(url, grid, operator=operator) + assert aligned.shape == (grid.rows, grid.cols) + + +def test_non_overlapping_result_is_entirely_nodata(tmp_path): + url = _write_raster(tmp_path / "s.tif", np.ones((10, 10)), bbox=(0, 0, 100, 100), + nodata=-9999.0) + aligned = _align(url, _grid((500, 500, 600, 600))) + values = np.asarray(aligned.values) + assert np.all((values == -9999.0) | np.isnan(values)) + + +def test_non_overlapping_result_carries_source_nodata(tmp_path): + """Categorical operators rely on _disscube_nodata to mark cells invalid.""" + url = _write_raster(tmp_path / "s.tif", np.ones((10, 10)), bbox=(0, 0, 100, 100), + nodata=-9999.0) + aligned = _align(url, _grid((500, 500, 600, 600)), operator="majority") + assert aligned.attrs.get("_disscube_nodata") == -9999.0 + + +def test_overlapping_grid_still_reads_data(tmp_path): + """Guard: the no-overlap shortcut must not swallow the normal path.""" + src = np.arange(100, dtype="float32").reshape(10, 10) + url = _write_raster(tmp_path / "s.tif", src, bbox=(0, 0, 100, 100)) + aligned = _align(url, _grid((0, 0, 100, 100))) + values = np.asarray(aligned.values) + assert np.isfinite(values).any() + assert not np.all(values == -9999.0) From 44e8671dd709305be0efbe5fcb581531cba22863 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Thu, 27 Aug 2026 17:40:36 -0300 Subject: [PATCH 06/17] docs: add BR-MANGUE domain example on the BDC tile mesh MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Companion to 01_dominio_estrutural.py, producing the same product over the national BDC_SM mesh with the master grid in BDC Albers instead of an ad hoc pixel mesh in the source's own CRS. The pair makes the trade-off concrete rather than theoretical. Measured on the same data: the ad hoc mesh runs 30 tiles in 80s under 1GB and reproduces the source bit for bit, while the BDC mesh runs 32 tiles in 100s peaking at 3.7GB and loses ~1% of cells per class to edge resampling — the cost of being on a canonical, shareable grid. The README now states when each is the right choice. Tile selection densifies the source footprint before projecting it: taking min/max of the four transformed corners overestimates the area, because straight edges become curves in another projection. That over-selection picked up 3 tiles that do not touch the data (35 vs the correct 32). They would now derive cheaply as empty, but selecting them at all is wasted work. Passes the full tile id (BDC_SM_027005) rather than the bare number, since ids are not unique across BDC levels. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- examples/README.md | 23 +- .../brmangue_dominio/02_dominio_bdc.py | 259 ++++++++++++++++++ 2 files changed, 279 insertions(+), 3 deletions(-) create mode 100644 examples/case_studies/brmangue_dominio/02_dominio_bdc.py diff --git a/examples/README.md b/examples/README.md index dc7dbcb..107879f 100644 --- a/examples/README.md +++ b/examples/README.md @@ -24,16 +24,33 @@ Dois estudos sobre a mesma área geográfica e grade. Standalone — não depende dos scripts de `setup/`; registra a própria grade e fontes. - `python examples/case_studies/brmangue_dominio/01_dominio_estrutural.py` — reclassifica a legenda de estados num raster de papéis, tile a tile, e escreve o GeoTIFF final. +- `python examples/case_studies/brmangue_dominio/02_dominio_bdc.py` — o mesmo produto + sobre a malha nacional **BDC_SM**, com a grade mestra em BDC Albers. -Demonstra três coisas de uso geral: +Demonstram três coisas de uso geral: - o operador genérico **`reclassify`** (tabela `{valor_origem: valor_destino}` como dado); - **`derive(tile_id=...)`** para processar uma extensão grande (391M células) com memória limitada — ver `docs/architecture/tiling.md`; - a fronteira entre o que é do cubo e o que é do projeto: a tabela de papéis e o cruzamento entre duas bandas ficam no exemplo, não no pacote. -Requer `geomosaic` (`pip install geomosaic`) e a variável `BRMANGUE_ENTRADA` -apontando para o diretório com os tiles ANADEM v2. +A diferença entre os dois está só na partição e no CRS de trabalho: + +| | `01_` | `02_` | +|---|---|---| +| Malha | ad hoc, 4096² em pixel | **BDC_SM**, canônica | +| Grade | EPSG:5880 (a da fonte) | BDC Albers | +| Reprojeção | nenhuma (caminho de identidade) | 5880 → Albers | +| Medido | 30 tiles, 80 s, < 1 GB | 32 tiles, 100 s, 3,7 GB de pico | + +Use a malha BDC quando os derivados forem virar patrimônio do cubo — os tiles têm +significado fora do script e alinham com o resto do ecossistema BDC. Use a ad hoc +quando forem só um passo intermediário: é mais leve e reproduz a fonte bit a bit +(reprojetar reamostra, então `02_` perde ~1% de células por classe nas bordas). + +Requerem `geomosaic` (`pip install geomosaic`) e a variável `BRMANGUE_ENTRADA` +apontando para o diretório com os tiles ANADEM v2. O `02_` também precisa de +`pip install "disscube[bdc]"` e das grades em `data/bdc_grids/`. ### 5. Estudo de caso: Acre (AC/5km) - `python examples/drivers/02_acre_5km.py` — drivers regionais Acre 5 km. diff --git a/examples/case_studies/brmangue_dominio/02_dominio_bdc.py b/examples/case_studies/brmangue_dominio/02_dominio_bdc.py new file mode 100644 index 0000000..c3a5b59 --- /dev/null +++ b/examples/case_studies/brmangue_dominio/02_dominio_bdc.py @@ -0,0 +1,259 @@ +""" +examples/case_studies/brmangue_dominio/02_dominio_bdc.py + +BR-MANGUE — domínio estrutural na malha nacional BDC (tiles SM, 30 m). + +Variante do `01_dominio_estrutural.py`. A diferença está inteiramente na +definição da partição e do CRS de trabalho: + + 01_ malha ad hoc (4096x4096 em pixel), grade em EPSG:5880. + Sem reprojeção — a fonte já está na grade alvo. + 02_ malha BDC_SM (canônica, nacional), grade em BDC Albers. + Reprojeta 5880 -> Albers dentro do GridAligner. + +Quando usar cada uma +-------------------- +A malha ad hoc é mais rápida e reproduz a fonte bit a bit, mas os recortes +só têm significado dentro do script que os criou. A malha BDC dá tiles +canônicos, compartilháveis entre projetos e alinhados com o resto do +ecossistema BDC — ao custo de uma reprojeção real e de um produto que já +não é idêntico à fonte (Albers reamostrado, não Polyconic). + +Escolha a malha BDC quando os derivados forem virar patrimônio do cubo; +a ad hoc quando forem só um passo intermediário. + +Nota sobre `elevacao_valida` +---------------------------- +Com reprojeção, a máscara de validade da elevação é reamostrada junto +(operador `mean` -> `Resampling.average`), então as bordas entre válido e +nodata ficam aproximadas. Em `01_` isso não acontece porque não há +reamostragem. É consequência de mudar de grade, não um defeito do cubo. + +Pré-requisitos +-------------- + - pip install geomosaic + - pip install "disscube[bdc]" (fiona, para ler os shapefiles BDC) + - grades BDC em data/bdc_grids/BDC_{SM,MD,LG}_V2.zip + - tiles ANADEM v2 em $BRMANGUE_ENTRADA/anadem_v2/ + +Usage: + export BRMANGUE_ENTRADA=/caminho/para/pymangue/dados/entrada + python examples/case_studies/brmangue_dominio/02_dominio_bdc.py +""" + +import os +import time +from pathlib import Path + +import numpy as np +import rasterio +from pyproj import CRS as ProjCRS, Transformer +from rasterio.windows import Window +from shapely.geometry import MultiPoint, shape + +from disscube.client import CubeClient +from disscube.models import GridSpec, SpatialSource, SpatialDerivation, Variable +from disscube.utils.grids import BDC_CRS + +try: + from geomosaic.core import build_mosaic_contract, write_vrt +except ImportError as exc: # pragma: no cover — exemplo, fora do pacote + raise SystemExit( + "Este exemplo precisa do geomosaic: pip install geomosaic" + ) from exc + +try: + import fiona +except ImportError as exc: # pragma: no cover + raise SystemExit( + 'Este exemplo precisa do fiona: pip install "disscube[bdc]"' + ) from exc + + +# ── Configuração ────────────────────────────────────────────────────────────── +ENTRADA = Path(os.environ.get("BRMANGUE_ENTRADA", "./dados/entrada")) +TILES_DIR = ENTRADA / "anadem_v2" +BDC_SM = "zip://data/bdc_grids/BDC_SM_V2.zip" +TRABALHO = Path("./data/brmangue_dominio_bdc") +SAIDA = TRABALHO / "dominio_estrutural_2024_bdc.tif" + +GRID_ID = "brmangue/30m_bdc" +RESOLUCAO = 30.0 +BANDA_ESTADO = 3 +BANDA_ELEV = 1 +NODATA_SAIDA = 255 +NODATA_ELEV = -9999.0 + +# ── Específico do BR-MANGUE — dado do projeto, não do disscube ──────────────── +PAPEIS_ESTADOS = { + "excluidos": ([0, 12], 0), + "mangue_inicial": ([1], 1), + "acomodacao_natural": ([2, 3, 4, 5, 6], 2), + "uso_manejado_candidato": ([7], 3), + "restricao_ou_barreira": ([8, 9, 10], 4), + "agua_aberta_nao_validada": ([11], 5), +} +MAPPING = {c: papel for _, (codes, papel) in PAPEIS_ESTADOS.items() for c in codes} + + +def main() -> None: + if not TILES_DIR.is_dir(): + raise RuntimeError( + f"Diretório de tiles não encontrado: {TILES_DIR}\n" + "Defina BRMANGUE_ENTRADA apontando para .../pymangue/dados/entrada" + ) + tiles_src = sorted(TILES_DIR.glob("*.tif")) + if not tiles_src: + raise RuntimeError(f"Nenhum .tif em {TILES_DIR}") + + t0 = time.perf_counter() + TRABALHO.mkdir(parents=True, exist_ok=True) + + # ── 1. geomosaic: VRTs na projeção nativa da fonte (EPSG:5880) ─────────── + print(f"\n[1/4] geomosaic: {len(tiles_src)} tiles -> contrato") + contract = build_mosaic_contract([str(t) for t in tiles_src]) + a, _b, ox, _d, e, oy = contract.mosaic_transform + fonte_bbox = ( + ox, oy - contract.mosaic_height * abs(e), + ox + contract.mosaic_width * a, oy, + ) + print(f" fonte: {contract.mosaic_height}x{contract.mosaic_width} @ {contract.crs}") + + vrt_estado = write_vrt(contract, str(TRABALHO / "estado.vrt"), band=BANDA_ESTADO) + vrt_elev = write_vrt(contract, str(TRABALHO / "elevacao.vrt"), band=BANDA_ELEV) + + # ── 2. Selecionar os tiles BDC_SM que cobrem a fonte ───────────────────── + print("[2/4] BDC: selecionando tiles SM sobre a extensão da fonte") + para_bdc = Transformer.from_crs( + ProjCRS.from_user_input(str(contract.crs)), + ProjCRS.from_user_input(BDC_CRS), + always_xy=True, + ) + # O footprint é densificado antes de projetar: transformar só os quatro + # cantos e tomar min/max devolve uma envoltória maior que a área real + # (as bordas viram curvas em outra projeção), o que seleciona tiles que + # não encostam no dado. Eles seriam derivados como vazios — correto, mas + # trabalho à toa. + passo = 200 + contorno = [] + x0, y0, x1, y1 = fonte_bbox + for i in range(passo + 1): + f = i / passo + contorno += [ + (x0 + (x1 - x0) * f, y0), (x0 + (x1 - x0) * f, y1), + (x0, y0 + (y1 - y0) * f), (x1, y0 + (y1 - y0) * f), + ] + px, py = para_bdc.transform([p[0] for p in contorno], [p[1] for p in contorno]) + alvo = MultiPoint(list(zip(px, py))).convex_hull + + selecionados = [] + with fiona.open(BDC_SM) as src: + for rec in src: + geom = shape(rec["geometry"]) + if geom.intersects(alvo): + selecionados.append((rec["properties"]["tile"], geom.bounds)) + if not selecionados: + raise RuntimeError("Nenhum tile BDC_SM cobre a extensão da fonte.") + selecionados.sort() + print(f" {len(selecionados)} tiles BDC_SM selecionados") + + # Extensão de saída = envoltória dos tiles escolhidos (já alinhada à malha). + out_minx = min(b[0] for _t, b in selecionados) + out_miny = min(b[1] for _t, b in selecionados) + out_maxx = max(b[2] for _t, b in selecionados) + out_maxy = max(b[3] for _t, b in selecionados) + altura = int(round((out_maxy - out_miny) / RESOLUCAO)) + largura = int(round((out_maxx - out_minx) / RESOLUCAO)) + print(f" saída: {altura}x{largura} @ BDC Albers, {RESOLUCAO:.0f}m") + + # ── 3. disscube: grade mestra BDC + fontes + tiles ─────────────────────── + print("[3/4] disscube: registrando grade, fontes e tiles") + cube = CubeClient( + catalog=str(TRABALHO / "catalog.db"), store=str(TRABALHO / "store") + ) + cube.register_grid(GridSpec( + id=GRID_ID, type="reference", crs=BDC_CRS, resolution=RESOLUCAO, + bbox=[out_minx, out_miny, out_maxx, out_maxy], + description="BR-MANGUE 30 m sobre a malha BDC (Albers)", + )) + cube.register_spatial_source(SpatialSource( + id="estado_2024", name="brm_state_2024", format="raster", + asset_url=str(vrt_estado), crs=str(contract.crs), + )) + cube.register_spatial_source(SpatialSource( + id="elevacao", name="elevation_m", format="raster", + asset_url=str(vrt_elev), crs=str(contract.crs), + )) + + # Tiles BDC registrados com o id canônico. São independentes de grade — + # o mesmo envelope serve qualquer grade no CRS do BDC. + for tile, bounds in selecionados: + cube.register_spatial_source(SpatialSource( + id=f"BDC_SM_{tile}", name=f"BDC SM Tile {tile}", format="raster", + asset_url="planned", crs=BDC_CRS, bbox=list(bounds), + )) + + deriv_papel = SpatialDerivation( + source_id="estado_2024", grid_id=GRID_ID, role="state", + variables=[Variable(name="papel", operator="reclassify", mapping=MAPPING)], + ) + deriv_elev = SpatialDerivation( + source_id="elevacao", grid_id=GRID_ID, role="driver", + variables=[Variable(name="elevacao", operator="mean")], + ) + + # ── 4. derive por tile BDC + escrita janelada ──────────────────────────── + print(f"[4/4] derive() por tile BDC -> {SAIDA}") + perfil = dict( + driver="GTiff", height=altura, width=largura, count=2, dtype="uint8", + nodata=NODATA_SAIDA, crs=BDC_CRS, + transform=rasterio.transform.Affine( + RESOLUCAO, 0.0, out_minx, 0.0, -RESOLUCAO, out_maxy + ), + tiled=True, blockxsize=512, blockysize=512, + compress="deflate", zlevel=6, bigtiff="IF_SAFER", sparse_ok=True, + ) + + with rasterio.open(SAIDA, "w", **perfil) as dst: + dst.set_band_description(1, "papel_dominio_2024") + dst.set_band_description(2, "elevacao_valida") + + for i, (tile, bounds) in enumerate(selecionados, start=1): + # id COMPLETO: ids simples não são únicos entre níveis BDC. + tile_id = f"BDC_SM_{tile}" + cube.derive(deriv_papel, tile_id=tile_id) + cube.derive(deriv_elev, tile_id=tile_id) + + papel = cube.load("papel", tile_id=tile_id).values + elev = cube.load("elevacao", tile_id=tile_id).values + + # Regra do projeto: acoplamento papel/elevação (duas bandas). + estado_valido = ~np.isnan(papel) + elev_valida = np.isfinite(elev) & (elev != NODATA_ELEV) + + c0 = int(round((bounds[0] - out_minx) / RESOLUCAO)) + r0 = int(round((out_maxy - bounds[3]) / RESOLUCAO)) + h, w = papel.shape + janela = Window(c0, r0, w, h) + + dst.write( + np.where(estado_valido, papel, NODATA_SAIDA).astype(np.uint8), + 1, window=janela, + ) + dst.write( + np.where( + estado_valido, elev_valida.astype(np.uint8), NODATA_SAIDA + ).astype(np.uint8), + 2, window=janela, + ) + + if i % 5 == 0 or i == len(selecionados): + print(f" {i}/{len(selecionados)} tiles " + f"({100 * i / len(selecionados):.0f}%)") + + print(f"\n=== domínio na malha BDC gerado em {time.perf_counter() - t0:.1f}s ===") + print(f" {SAIDA} ({altura}x{largura}, 2 bandas uint8, BDC Albers)") + + +if __name__ == "__main__": + main() From 820a5fc66ca6e0b85a9d0324844c440873db8463 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Thu, 27 Aug 2026 17:54:15 -0300 Subject: [PATCH 07/17] perf: select BDC tiles by valid-data mask, not by extent MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The mosaic's extent is a rectangle but the valid data is a coastal strip, so most tiles intersecting the extent cover nothing but nodata: 17 of 32 derived to fully empty rasters. Tile selection now takes one decimated pass over the source to build a low-resolution validity mask, and keeps only the tiles that touch it. Memory for that pass is bounded by the target width, not by the source size. Measured on the same data: 32 -> 15 tiles, 99.9s -> 56.9s, peak RSS 3.70GB -> 2.81GB, with class counts identical to the previous run (zero valid cells lost). Disk was never the problem — the empty Zarrs compressed to almost nothing, so the store only went 257MB -> 252MB. The win is compute, not storage. The mask is dilated by one cell: nearest-neighbour decimation can drop thin features, and the mangrove strip is a few pixels wide in places. The costs are asymmetric — a false positive derives one tile for nothing, a false negative silently loses data from the product. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- .../brmangue_dominio/02_dominio_bdc.py | 85 +++++++++++++++++-- 1 file changed, 80 insertions(+), 5 deletions(-) diff --git a/examples/case_studies/brmangue_dominio/02_dominio_bdc.py b/examples/case_studies/brmangue_dominio/02_dominio_bdc.py index c3a5b59..8049c85 100644 --- a/examples/case_studies/brmangue_dominio/02_dominio_bdc.py +++ b/examples/case_studies/brmangue_dominio/02_dominio_bdc.py @@ -96,6 +96,59 @@ MAPPING = {c: papel for _, (codes, papel) in PAPEIS_ESTADOS.items() for c in codes} +def _mascara_valida(vrt: Path, largura_alvo: int = 2048): + """ + Máscara booleana, em baixa resolução, de onde a fonte tem dado válido. + + Uma única leitura decimada — a memória fica limitada por ``largura_alvo``, + não pelo tamanho da fonte. A máscara é dilatada em uma célula porque a + decimação por vizinho mais próximo pode perder feições finas (a faixa de + mangue tem poucos pixels de largura em muitos trechos), e aqui um falso + positivo custa um tile derivado à toa enquanto um falso negativo custa + dado faltando no produto final. + + Returns + ------- + (mascara, transform, resolucao_da_mascara) + """ + from scipy.ndimage import binary_dilation + + with rasterio.open(vrt) as ds: + fator = max(1, int(np.ceil(ds.width / largura_alvo))) + altura = max(1, ds.height // fator) + largura = max(1, ds.width // fator) + baixa = ds.read( + 1, out_shape=(altura, largura), resampling=rasterio.enums.Resampling.nearest + ) + nodata = ds.nodata + transform = ds.transform + resolucao = abs(transform.a) * (ds.width / largura) + + valida = np.isfinite(baixa) + if nodata is not None and np.isfinite(nodata): + valida &= baixa != nodata + return binary_dilation(valida), transform, resolucao + + +def _tem_dado(bounds, mascara, transform, resolucao, para_fonte) -> bool: + """True se o tile (bounds em BDC Albers) toca algum pixel válido da fonte.""" + minx, miny, maxx, maxy = bounds + xs, ys = para_fonte.transform( + [minx, minx, maxx, maxx], [miny, maxy, miny, maxy] + ) + c0 = int(np.floor((min(xs) - transform.c) / resolucao)) + c1 = int(np.ceil((max(xs) - transform.c) / resolucao)) + r0 = int(np.floor((transform.f - max(ys)) / resolucao)) + r1 = int(np.ceil((transform.f - min(ys)) / resolucao)) + + h, w = mascara.shape + r0, r1 = max(0, r0), min(h, r1) + c0, c1 = max(0, c0), min(w, c1) + if r0 >= r1 or c0 >= c1: + return False + return bool(mascara[r0:r1, c0:c1].any()) + + def main() -> None: if not TILES_DIR.is_dir(): raise RuntimeError( @@ -146,16 +199,38 @@ def main() -> None: px, py = para_bdc.transform([p[0] for p in contorno], [p[1] for p in contorno]) alvo = MultiPoint(list(zip(px, py))).convex_hull - selecionados = [] + candidatos = [] with fiona.open(BDC_SM) as src: for rec in src: geom = shape(rec["geometry"]) if geom.intersects(alvo): - selecionados.append((rec["properties"]["tile"], geom.bounds)) - if not selecionados: + candidatos.append((rec["properties"]["tile"], geom.bounds)) + if not candidatos: raise RuntimeError("Nenhum tile BDC_SM cobre a extensão da fonte.") - selecionados.sort() - print(f" {len(selecionados)} tiles BDC_SM selecionados") + candidatos.sort() + + # A extensão do mosaico é um retângulo, mas o dado válido é só a faixa + # costeira — mais da metade dos tiles que a intersectam cobrem apenas + # nodata. Derivá-los é trabalho e disco jogados fora (medido: 17 de 32 + # tiles saíam 100% vazios). Uma passada decimada sobre a fonte dá a + # máscara de onde há dado de verdade, e a seleção passa a ser por ela. + print(f" {len(candidatos)} tiles intersectam a extensão; " + "lendo máscara de dado válido") + mascara, m_transform, m_res = _mascara_valida(vrt_estado) + para_fonte = Transformer.from_crs( + ProjCRS.from_user_input(BDC_CRS), + ProjCRS.from_user_input(str(contract.crs)), + always_xy=True, + ) + + selecionados = [ + (tile, bounds) for tile, bounds in candidatos + if _tem_dado(bounds, mascara, m_transform, m_res, para_fonte) + ] + if not selecionados: + raise RuntimeError("Nenhum tile BDC_SM contém dado válido.") + print(f" {len(selecionados)} tiles com dado válido " + f"({len(candidatos) - len(selecionados)} descartados por serem só nodata)") # Extensão de saída = envoltória dos tiles escolhidos (já alinhada à malha). out_minx = min(b[0] for _t, b in selecionados) From 14c6d59852e5d98a2b1e29b091b8e57db37db840 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Thu, 27 Aug 2026 19:57:26 -0300 Subject: [PATCH 08/17] perf: select tiles by valid-data mask in the ad hoc mesh example too MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 820a5fc fixed this for the BDC example but left 01_dominio_estrutural.py deriving the full rectangular mesh, where half the tiles cover nothing but nodata — the same waste, reported by the user after running it. Same approach: one decimated pass over the source builds a low-resolution validity mask, and only tiles touching it are derived. Simpler here than in 02_, since the tiles are pixel windows in the source itself, so the mask is indexed directly with no CRS transform. Measured: 30 -> 15 tiles, 80s -> 36.3s, and the output is still bit for bit identical to the reference GDAL pipeline over all 391M cells in both bands. Skipping is safe precisely because the mask comes from the state band: the reference implementation writes 255 to BOTH bands wherever state is nodata, so a tile with no valid state contributes nothing but the nodata the output already holds. README's comparison table updated with the new figures. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- examples/README.md | 7 ++- .../brmangue_dominio/01_dominio_estrutural.py | 49 ++++++++++++++++++- 2 files changed, 53 insertions(+), 3 deletions(-) diff --git a/examples/README.md b/examples/README.md index 107879f..fef47c1 100644 --- a/examples/README.md +++ b/examples/README.md @@ -41,7 +41,12 @@ A diferença entre os dois está só na partição e no CRS de trabalho: | Malha | ad hoc, 4096² em pixel | **BDC_SM**, canônica | | Grade | EPSG:5880 (a da fonte) | BDC Albers | | Reprojeção | nenhuma (caminho de identidade) | 5880 → Albers | -| Medido | 30 tiles, 80 s, < 1 GB | 32 tiles, 100 s, 3,7 GB de pico | +| Medido | 15 tiles, 36 s, < 1 GB | 15 tiles, 57 s, 2,8 GB de pico | + +Ambos selecionam os recortes por uma máscara de dado válido, não pela extensão: +o mosaico é um retângulo mas o dado é uma faixa costeira, então metade dos +recortes cobriria só nodata (30→15 e 32→15). Uma leitura decimada da fonte +resolve isso, e o resultado é bit a bit o mesmo — verificado. Use a malha BDC quando os derivados forem virar patrimônio do cubo — os tiles têm significado fora do script e alinham com o resto do ecossistema BDC. Use a ad hoc diff --git a/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py b/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py index b4a583c..68c735b 100644 --- a/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py +++ b/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py @@ -89,6 +89,35 @@ MAPPING = {c: papel for _, (codes, papel) in PAPEIS_ESTADOS.items() for c in codes} +def _mascara_valida(vrt: Path, largura_alvo: int = 2048): + """ + Máscara booleana, em baixa resolução, de onde a fonte tem dado válido, + junto com o fator de decimação usado. + + Uma única leitura decimada — a memória fica limitada por ``largura_alvo``, + não pelo tamanho da fonte. A máscara é dilatada em uma célula porque a + decimação por vizinho mais próximo pode perder feições finas (a faixa de + mangue tem poucos pixels de largura em muitos trechos), e aqui um falso + positivo custa um recorte derivado à toa enquanto um falso negativo custa + dado faltando no produto final. + """ + from scipy.ndimage import binary_dilation + + with rasterio.open(vrt) as ds: + fator = max(1, int(np.ceil(ds.width / largura_alvo))) + baixa = ds.read( + 1, + out_shape=(max(1, ds.height // fator), max(1, ds.width // fator)), + resampling=rasterio.enums.Resampling.nearest, + ) + nodata = ds.nodata + + valida = np.isfinite(baixa) + if nodata is not None and np.isfinite(nodata): + valida &= baixa != nodata + return binary_dilation(valida), fator + + def main() -> None: if not TILES_DIR.is_dir(): raise RuntimeError( @@ -131,12 +160,27 @@ def main() -> None: asset_url=str(vrt_elev), crs=str(contract.crs), )) + # A extensão do mosaico é um retângulo, mas o dado válido é só a faixa + # costeira: boa parte dos recortes cobre apenas nodata, e derivá-los é + # tempo jogado fora. Uma passada decimada sobre a fonte dá a máscara de + # onde há dado de verdade, e só esses recortes entram no loop. + mascara, fator = _mascara_valida(vrt_estado) + # Cada recorte vira um SpatialSource {grid_id}_{tile_id} carregando só o # bbox — é assim que derive(tile_id=...) descobre a janela a processar. - janelas = [] + janelas, vazios = [], 0 for r0 in range(0, altura, TILE): for c0 in range(0, largura, TILE): h, w = min(TILE, altura - r0), min(TILE, largura - c0) + + # Os recortes são janelas em pixel na própria fonte, então basta + # indexar a máscara na mesma proporção — sem transformar CRS. + sub = mascara[r0 // fator:-(-(r0 + h) // fator), + c0 // fator:-(-(c0 + w) // fator)] + if not sub.any(): + vazios += 1 + continue + tile_id = f"R{r0:05d}C{c0:05d}" minx, maxy = ox + c0 * resolution, oy - r0 * abs(e) cube.register_spatial_source(SpatialSource( @@ -145,7 +189,8 @@ def main() -> None: bbox=[minx, maxy - h * abs(e), minx + w * resolution, maxy], )) janelas.append((tile_id, r0, c0, h, w)) - print(f" {len(janelas)} tiles de {TILE}x{TILE}") + print(f" {len(janelas)} tiles de {TILE}x{TILE} com dado válido " + f"({vazios} descartados por serem só nodata)") deriv_papel = SpatialDerivation( source_id="estado_2024", grid_id=GRID_ID, role="state", From fd1221592a8aab82a91a542c2c827ee465d0b242 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Thu, 27 Aug 2026 20:53:24 -0300 Subject: [PATCH 09/17] fix: use the shared catalog and store in the brmangue examples MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Both examples created their own catalog.db and store under a per-example directory, unlike the other twelve examples, which all use CubeClient(catalog="catalog.db", store="./data/"). That was not merely cosmetic: with a private catalog the derived variables never showed up in cube.search() alongside the rest, and could not be combined through to_lucc_data() — losing the reason to catalog them at all. Each example was an island. Now both share the standard catalog and store, and follow the repo's directory conventions: VRTs go to data/raw/brmangue/ next to the other raw inputs, the GeoTIFFs to data/. Verified by running both: one catalog holds the two grids (brmangue/30m in EPSG:5880 and brmangue/30m_bdc in BDC Albers) with 60 DerivedVariable entries — papel and elevacao, 15 tiles each, on both grids. The two meshes are now comparable within the same cube. Note that both examples write the same VRTs and register the same SpatialSource ids, deriving them from the same ANADEM tiles. That is idempotent today, but the second run would silently overwrite the first if they ever diverge. Reported by the user. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- .../brmangue_dominio/01_dominio_estrutural.py | 18 ++++++++++++------ .../brmangue_dominio/02_dominio_bdc.py | 19 +++++++++++-------- 2 files changed, 23 insertions(+), 14 deletions(-) diff --git a/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py b/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py index 68c735b..6df3037 100644 --- a/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py +++ b/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py @@ -66,8 +66,14 @@ # ── Configuração ────────────────────────────────────────────────────────────── ENTRADA = Path(os.environ.get("BRMANGUE_ENTRADA", "./dados/entrada")) TILES_DIR = ENTRADA / "anadem_v2" -TRABALHO = Path("./data/brmangue_dominio") -SAIDA = TRABALHO / "dominio_estrutural_2024.tif" + +# Catálogo e store compartilhados, como em todos os demais exemplos: os +# derivados entram no MESMO cubo e aparecem em cube.search() junto com as +# outras variáveis — que é o ponto de catalogá-los. +CATALOGO = "catalog.db" +STORE = "./data/" +VRT_DIR = Path("./data/raw/brmangue") +SAIDA = Path("./data/brmangue_dominio_2024.tif") GRID_ID = "brmangue/30m" TILE = 4096 # lado do recorte processado por vez @@ -129,7 +135,7 @@ def main() -> None: raise RuntimeError(f"Nenhum .tif em {TILES_DIR}") t0 = time.perf_counter() - TRABALHO.mkdir(parents=True, exist_ok=True) + VRT_DIR.mkdir(parents=True, exist_ok=True) # ── 1. geomosaic: contrato + um VRT por banda ──────────────────────────── print(f"\n[1/3] geomosaic: {len(tiles)} tiles -> contrato") @@ -139,13 +145,13 @@ def main() -> None: resolution = a print(f" grade mestra: {altura}x{largura} @ {contract.crs}") - vrt_estado = write_vrt(contract, str(TRABALHO / "estado.vrt"), band=BANDA_ESTADO) - vrt_elev = write_vrt(contract, str(TRABALHO / "elevacao.vrt"), band=BANDA_ELEV) + vrt_estado = write_vrt(contract, str(VRT_DIR / "estado.vrt"), band=BANDA_ESTADO) + vrt_elev = write_vrt(contract, str(VRT_DIR / "elevacao.vrt"), band=BANDA_ELEV) # ── 2. disscube: grade, fontes e tiles ─────────────────────────────────── print("[2/3] disscube: registrando grade, fontes e tiles") cube = CubeClient( - catalog=str(TRABALHO / "catalog.db"), store=str(TRABALHO / "store") + catalog=CATALOGO, store=STORE ) cube.register_grid(GridSpec( id=GRID_ID, type="reference", crs=str(contract.crs), resolution=resolution, diff --git a/examples/case_studies/brmangue_dominio/02_dominio_bdc.py b/examples/case_studies/brmangue_dominio/02_dominio_bdc.py index 8049c85..77774d8 100644 --- a/examples/case_studies/brmangue_dominio/02_dominio_bdc.py +++ b/examples/case_studies/brmangue_dominio/02_dominio_bdc.py @@ -74,8 +74,13 @@ ENTRADA = Path(os.environ.get("BRMANGUE_ENTRADA", "./dados/entrada")) TILES_DIR = ENTRADA / "anadem_v2" BDC_SM = "zip://data/bdc_grids/BDC_SM_V2.zip" -TRABALHO = Path("./data/brmangue_dominio_bdc") -SAIDA = TRABALHO / "dominio_estrutural_2024_bdc.tif" +# Catálogo e store compartilhados, como em todos os demais exemplos: os +# derivados entram no MESMO cubo e aparecem em cube.search() junto com as +# outras variáveis — que é o ponto de catalogá-los. +CATALOGO = "catalog.db" +STORE = "./data/" +VRT_DIR = Path("./data/raw/brmangue") +SAIDA = Path("./data/brmangue_dominio_2024_bdc.tif") GRID_ID = "brmangue/30m_bdc" RESOLUCAO = 30.0 @@ -160,7 +165,7 @@ def main() -> None: raise RuntimeError(f"Nenhum .tif em {TILES_DIR}") t0 = time.perf_counter() - TRABALHO.mkdir(parents=True, exist_ok=True) + VRT_DIR.mkdir(parents=True, exist_ok=True) # ── 1. geomosaic: VRTs na projeção nativa da fonte (EPSG:5880) ─────────── print(f"\n[1/4] geomosaic: {len(tiles_src)} tiles -> contrato") @@ -172,8 +177,8 @@ def main() -> None: ) print(f" fonte: {contract.mosaic_height}x{contract.mosaic_width} @ {contract.crs}") - vrt_estado = write_vrt(contract, str(TRABALHO / "estado.vrt"), band=BANDA_ESTADO) - vrt_elev = write_vrt(contract, str(TRABALHO / "elevacao.vrt"), band=BANDA_ELEV) + vrt_estado = write_vrt(contract, str(VRT_DIR / "estado.vrt"), band=BANDA_ESTADO) + vrt_elev = write_vrt(contract, str(VRT_DIR / "elevacao.vrt"), band=BANDA_ELEV) # ── 2. Selecionar os tiles BDC_SM que cobrem a fonte ───────────────────── print("[2/4] BDC: selecionando tiles SM sobre a extensão da fonte") @@ -243,9 +248,7 @@ def main() -> None: # ── 3. disscube: grade mestra BDC + fontes + tiles ─────────────────────── print("[3/4] disscube: registrando grade, fontes e tiles") - cube = CubeClient( - catalog=str(TRABALHO / "catalog.db"), store=str(TRABALHO / "store") - ) + cube = CubeClient(catalog=CATALOGO, store=STORE) cube.register_grid(GridSpec( id=GRID_ID, type="reference", crs=BDC_CRS, resolution=RESOLUCAO, bbox=[out_minx, out_miny, out_maxx, out_maxy], From f57056e8c46d604162bd544eeddc5f3727ae70ad Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Fri, 28 Aug 2026 08:34:06 -0300 Subject: [PATCH 10/17] novo exemplo usando halo --- .../03_carregar_no_haloexec.py | 262 ++++++++++++++++++ 1 file changed, 262 insertions(+) create mode 100644 examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py diff --git a/examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py b/examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py new file mode 100644 index 0000000..1fb136f --- /dev/null +++ b/examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py @@ -0,0 +1,262 @@ +""" +examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py + +BR-MANGUE — carrega as variáveis derivadas num MemmapRasterWorkspace do +haloexec, prontas para um modelo rodar em disco. + +É o par em disco do `maranhao/03_brmangue_simulate.py`: aquele entrega o +cubo em RAM via `to_lucc_data()` -> `RasterBackend`; este entrega em disco +via workspace memmap, para domínios que não cabem na memória. + +NÃO roda simulação — só deixa o workspace pronto e prova que os dados +chegaram corretos, inclusive nas costuras entre tiles. + +Por padrão usa a grade **BDC Albers** (`brmangue/30m_bdc`, do exemplo 02), +por ser a malha canônica: tiles BDC_SM compartilháveis entre projetos. +Troque `GRID_ID` para `brmangue/30m` para usar a malha ad hoc do 01. + +Pré-requisitos: + - python examples/case_studies/brmangue_dominio/02_dominio_bdc.py + - pip install "haloexec[zarr]" + +Usage: + python examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py + +Por que o carregamento é feito à mão aqui +------------------------------------------ +As variáveis foram derivadas TILE A TILE (`derive(tile_id=...)`), então +existe um Zarr por tile, não um só. Nenhum dos dois lados resolve isso hoje: + + * `CubeClient.load(nome, grid_id=...)` levanta ValueError quando há mais + de um tile — é preciso pedir um tile específico (comportamento + documentado em docs/architecture/tiling.md: "mosaico automático não + está implementado"); + * `haloexec.load_zarr_into_workspace()` recebe UM store e exige que o + shape dele seja igual ao do workspace inteiro — não há equivalente + multi-tile como `load_geotiffs_into_workspace` é para GeoTIFF. + +Então este exemplo faz a costura: descobre no catálogo quais tiles existem, +usa o bbox de cada um (registrado como SpatialSource) para achar sua posição +na grade mestra — a mesma ideia de `build_mosaic_contract` do geomosaic, mas +partindo da geometria catalogada em vez do nome do arquivo — e monta cada +bloco do workspace lendo dos tiles que o cobrem. + +Buracos na malha são esperados (o MapBiomas não exporta onde não há costa) +e viram nodata, não erro. +""" + +from pathlib import Path + +import numpy as np + +from disscube.client import CubeClient +from disscube.models import GridSpec + +try: + import zarr + from haloexec import Block, MemmapRasterWorkspace +except ImportError as exc: # pragma: no cover — exemplo, fora do pacote + raise SystemExit( + 'Este exemplo precisa do haloexec: pip install "haloexec[zarr]"' + ) from exc + + +# ── Configuração ────────────────────────────────────────────────────────────── +CATALOGO = "catalog.db" +STORE = "./data/" +GRID_ID = "brmangue/30m_bdc" # malha BDC (exemplo 02); ou "brmangue/30m" +VARIAVEIS = ["papel", "elevacao"] +WS_DIR = Path("./data/workspace_brmangue_bdc") + +# float32 em vez de float64: o workspace mantém DOIS slots (double buffer), +# então cada variável ocupa 2x o tamanho da grade em disco. +DTYPE = "float32" +BLOCK_H = BLOCK_W = 512 +HALO = 2 + + +def _tiles_da_variavel(cube: CubeClient, nome: str, grid: GridSpec) -> list[dict]: + """Descobre os tiles de ``nome`` e onde cada um cai na grade mestra. + + A posição vem do bbox do SpatialSource do tile — o mesmo bbox que + `derive(tile_id=...)` usou para recortar — convertido para deslocamento + em pixel contra a origem da grade. O id do SpatialSource pode ser + ``{grid_id}_{tile_id}`` (malha própria) ou o próprio ``tile_id`` quando + já é um id canônico BDC. + """ + achados = [] + for d in cube.search(grid=grid.id): + if d.name != nome or not d.tile_id: + continue + src = (cube.catalog.get_spatial_source(f"{grid.id}_{d.tile_id}") + or cube.catalog.get_spatial_source(d.tile_id)) + if src is None or not src.bbox: + raise RuntimeError( + f"tile {d.tile_id} sem SpatialSource com bbox — não dá para " + "saber onde ele cai na grade mestra" + ) + arr = zarr.open(d.asset_url, mode="r")[nome] + achados.append({ + "tile_id": d.tile_id, + "url": d.asset_url, + "r0": int(round((grid.bbox[3] - src.bbox[3]) / grid.resolution)), + "c0": int(round((src.bbox[0] - grid.bbox[0]) / grid.resolution)), + "h": arr.shape[0], + "w": arr.shape[1], + }) + return sorted(achados, key=lambda t: (t["r0"], t["c0"])) + + +def _carregar(ws, nome: str, tiles: list[dict]) -> tuple[int, int]: + """Preenche ``nome`` no workspace, bloco a bloco, a partir de N tiles. + + Percorre os blocos do workspace (e não os tiles) porque um bloco pode + cair sobre dois tiles vizinhos, ou sobre um buraco da malha: montá-lo a + partir de tudo que o cobre é o que faz a costura ficar correta. + """ + abertos = {t["tile_id"]: zarr.open(t["url"], mode="r")[nome] for t in tiles} + com = sem = 0 + + for blk in ws.blocks(): + buf = np.full((blk.r1 - blk.r0, blk.c1 - blk.c0), np.nan, dtype=DTYPE) + tocou = False + for t in tiles: + r0 = max(blk.r0, t["r0"]); r1 = min(blk.r1, t["r0"] + t["h"]) + c0 = max(blk.c0, t["c0"]); c1 = min(blk.c1, t["c0"] + t["w"]) + if r0 >= r1 or c0 >= c1: + continue + buf[r0 - blk.r0:r1 - blk.r0, c0 - blk.c0:c1 - blk.c0] = np.asarray( + abertos[t["tile_id"]][r0 - t["r0"]:r1 - t["r0"], + c0 - t["c0"]:c1 - t["c0"]] + ) + tocou = True + ws.write_block_to_read_slot(blk, nome, buf) + com += tocou + sem += not tocou + + ws.flush() + return com, sem + + +def _provar_costura(ws, nome: str, esq: dict, dir_: dict) -> tuple[bool, str]: + """Prova que a janela com halo sobre a costura entre DOIS Zarr está certa. + + A comparação é feita contra cada arquivo Zarr de origem separadamente — + o lado esquerdo contra o Zarr de ``esq``, o direito contra o de ``dir_``, + célula a célula. Note o reindexamento: a primeira coluna do vizinho é a + coluna 0 DELE, não a coluna global — é exatamente aí que um erro de + offset apareceria. + + A linha é escolhida onde os dois lados DIFEREM. Sem isso a prova é fraca: + numa região toda 0.0 ou toda NaN a comparação passa sem provar nada. + """ + za = zarr.open(esq["url"], mode="r")[nome] + zb = zarr.open(dir_["url"], mode="r")[nome] + col = dir_["c0"] # coluna global onde o vizinho começa + + # linhas em que ambos os tiles existem + r_ini = max(esq["r0"], dir_["r0"]) + HALO + 1 + r_fim = min(esq["r0"] + esq["h"], dir_["r0"] + dir_["h"]) - HALO - 9 + if r_fim <= r_ini: + return True, "sem sobreposição vertical suficiente" + + ult = np.asarray(za[r_ini - esq["r0"]:r_fim - esq["r0"], esq["w"] - 1]) + pri = np.asarray(zb[r_ini - dir_["r0"]:r_fim - dir_["r0"], 0]) + dif = np.where(np.isfinite(ult) & np.isfinite(pri) & (ult != pri))[0] + if dif.size == 0: + return True, "nenhuma linha com lados diferentes (prova fraca, ignorada)" + + linha = r_ini + int(dif[len(dif) // 2]) + blk = Block(r0=linha - 4, r1=linha + 4, c0=col - 4, c1=col + 4) + jan = ws.read_block_with_halo(blk, boundary_value=np.nan)[nome] + + jr0, jc0 = blk.r0 - HALO, blk.c0 - HALO + for jr in range(jan.shape[0]): + for jc in range(jan.shape[1]): + gr, gc = jr0 + jr, jc0 + jc + if gc < col: + fonte, t = za, esq + else: + fonte, t = zb, dir_ + lr, lc = gr - t["r0"], gc - t["c0"] + if not (0 <= lr < t["h"] and 0 <= lc < t["w"]): + continue + esperado = np.float32(np.asarray(fonte[lr, lc])) + obtido = jan[jr, jc] + if not ((np.isnan(esperado) and np.isnan(obtido)) or esperado == obtido): + return False, (f"janela[{jr},{jc}] (global r={gr} c={gc}, tile " + f"{t['tile_id']} local r={lr} c={lc}): " + f"esperado {esperado}, obtido {obtido}") + + n = len(np.unique(jan[np.isfinite(jan)])) + return True, f"linha r={linha}, {n} valores distintos na janela" + + +def main() -> None: + cube = CubeClient(catalog=CATALOGO, store=STORE) + grid = cube.catalog.get_grid(GRID_ID) + if grid is None: + raise SystemExit( + f"Grade {GRID_ID!r} não encontrada. Rode antes:\n" + " python examples/case_studies/brmangue_dominio/02_dominio_bdc.py" + ) + + print(f"\n[1/3] catálogo: {GRID_ID} = {grid.rows}x{grid.cols} @ {grid.crs[:40]}") + por_var = {} + for nome in VARIAVEIS: + tiles = _tiles_da_variavel(cube, nome, grid) + if not tiles: + raise SystemExit(f"nenhum tile derivado para {nome!r} em {GRID_ID}") + por_var[nome] = tiles + print(f" {nome}: {len(tiles)} tiles de {tiles[0]['h']}x{tiles[0]['w']}") + + try: + cube.load(VARIAVEIS[0], grid_id=GRID_ID) + print(" (load() sem tile_id funcionou — havia um tile só)") + except ValueError: + print(" load() sem tile_id recusa multi-tile, como esperado") + + gib = grid.rows * grid.cols * np.dtype(DTYPE).itemsize * len(VARIAVEIS) * 2 / 1024**3 + print(f"\n[2/3] workspace {grid.rows}x{grid.cols} {DTYPE} " + f"({len(VARIAVEIS)} arrays x 2 slots = {gib:.1f} GB em disco)") + ws = MemmapRasterWorkspace.create( + WS_DIR, shape=(grid.rows, grid.cols), + arrays={n: DTYPE for n in VARIAVEIS}, + block_h=BLOCK_H, block_w=BLOCK_W, halo=HALO, + ) + for nome, tiles in por_var.items(): + com, sem = _carregar(ws, nome, tiles) + print(f" {nome}: {com} blocos com dado, {sem} inteiramente vazios") + + # ── prova de costura entre arquivos Zarr distintos ────────────────────── + print("\n[3/3] costuras entre tiles (cada lado conferido contra o SEU Zarr)") + ws2 = MemmapRasterWorkspace(WS_DIR) + nome = VARIAVEIS[0] + tiles = por_var[nome] + ok = fracas = falhas = 0 + for a in tiles: + b = next((o for o in tiles + if o["r0"] == a["r0"] and o["c0"] == a["c0"] + a["w"]), None) + if b is None: + continue + passou, msg = _provar_costura(ws2, nome, a, b) + rotulo = f"{a['tile_id']} | {b['tile_id']}" + if not passou: + falhas += 1 + print(f" FALHA {rotulo}: {msg}") + elif "fraca" in msg or "sem sobreposição" in msg: + fracas += 1 + print(f" pulada {rotulo}: {msg}") + else: + ok += 1 + print(f" OK {rotulo}: {msg}") + print(f"\n {ok} costuras provadas, {fracas} sem contraste, {falhas} FALHAS") + + print(f"\n=== workspace pronto em {WS_DIR} ===") + print(f" shape={ws2.shape} arrays={list(ws2.metadata['arrays'])} " + f"blocos={len(ws2.blocks())} halo={ws2.halo}") + print(" (nenhum modelo foi executado — o workspace está pronto para receber um)") + + +if __name__ == "__main__": + main() From 0cdf099d9d271c075fdefa6a98e593bcd439e284 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Fri, 28 Aug 2026 09:21:58 -0300 Subject: [PATCH 11/17] feat: add CubeClient.tile_layout() and use it in the haloexec example MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A variable derived tile by tile exists as N Zarr stores, and load() refuses that case because there is no automatic mosaic (docs/architecture/tiling.md). Until now the package stated the problem and offered no way out: the caller hit a ValueError listing the tiles and was on their own. tile_layout() answers the question load() cannot — where each piece sits on the master grid — without assembling anything and without imposing a destination. It returns plain data (path plus position), not an object from another package, so whoever consumes it decides: load into memory, write to a disk workspace, or just inspect coverage. Same arrangement geomosaic already uses when it hands back tile_offsets. Dimensions come from the registered bbox, not from the file: it is the same bbox derive(tile_id=...) cropped with, so it states the intended position. A store whose contents disagree is an inconsistency to catch, not to paper over. A variable with no tiles returns a single item covering the whole grid, so callers can treat both cases the same way. This makes the example 03 collapse from ~100 lines of offset arithmetic to six, and it now demonstrates use rather than working around a limitation. Its seam checks moved to haloexec as tests (f22057d) — they were proving the assembly code, not the data, and re-running them on every load protects nothing. Verified after the rewrite: the same 10 seams on the real BDC grid still check out against the source Zarrs, 0 divergences. Tests: 13, aimed at the position arithmetic, which is where an error hides — a wrong sign on row_off flips the mosaic vertically and nothing raises. Also covers both tile registration conventions, the global case, and the errors. Full suite 127 passing. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- disscube/client/cube_client.py | 109 ++++++++ .../03_carregar_no_haloexec.py | 234 +++--------------- tests/test_tile_layout.py | 175 +++++++++++++ 3 files changed, 317 insertions(+), 201 deletions(-) create mode 100644 tests/test_tile_layout.py diff --git a/disscube/client/cube_client.py b/disscube/client/cube_client.py index 64f5f65..91a97ac 100644 --- a/disscube/client/cube_client.py +++ b/disscube/client/cube_client.py @@ -270,6 +270,115 @@ def _exists(d: DerivedVariable) -> bool: derived = static[0] return xr.open_zarr(derived.asset_url, consolidated=False)[derived.name] + def tile_layout( + self, + variable_id: str, + grid_id: str, + ) -> List[dict]: + """ + Onde cada pedaço de uma variável derivada cai na grade mestra. + + Responde a pergunta que ``load()`` não responde: uma variável + derivada tile a tile existe como N stores Zarr, e ``load()`` recusa + o caso multi-tile porque não há mosaico automático (ver + ``docs/architecture/tiling.md``). Este método devolve a informação + necessária para montá-la — sem montar nada, e sem impor um destino. + + O retorno é **dado puro**, não um objeto de outro pacote: uma lista + de dicionários com caminho e posição. Quem consome decide o que + fazer — carregar em memória, escrever num workspace em disco, + inspecionar a cobertura. É o mesmo arranjo que o ``geomosaic`` usa + ao devolver ``tile_offsets`` sem conhecer quem vai lê-los. + + Uma variável derivada sem tiles (``tile_id`` nulo, partição + ``global``) devolve uma lista de UM elemento cobrindo a grade + inteira, para que quem consome trate os dois casos igual. + + Parameters + ---------- + variable_id : str + Nome da variável derivada (ex.: ``"papel"``). + grid_id : str + Grade mestra sobre a qual os pedaços se posicionam. + + Returns + ------- + list[dict] + Ordenada por ``(row_off, col_off)``. Cada item tem: + + - ``tile_id`` — identificador do tile, ou ``None`` se global + - ``variable`` — nome da variável dentro do store Zarr + - ``url`` — caminho do store + - ``row_off`` — linha, em pixel, onde o pedaço começa na grade + - ``col_off`` — coluna, em pixel + - ``height`` — altura do pedaço, em pixel + - ``width`` — largura do pedaço, em pixel + + As dimensões vêm do ``bbox`` registrado, não do arquivo: é o + mesmo ``bbox`` que ``derive(tile_id=...)`` usou para recortar, + então descrevem a posição pretendida. Um store cujo conteúdo + divirja disso é inconsistência a detectar, não a mascarar. + + Raises + ------ + ValueError + Se a grade não existir, se a variável não tiver nenhum derivado + nela, ou se um tile não tiver ``SpatialSource`` com ``bbox`` — + sem o bbox não há como saber onde o pedaço cai. + """ + grid = self.catalog.get_grid(grid_id) + if grid is None: + raise ValueError(f"Grid not found: {grid_id}") + + derived = [ + d for d in self.catalog.search_derived_variables(grid_id=grid_id) + if d.name == variable_id + ] + if not derived: + raise ValueError( + f"Derived variable not found: {variable_id} on grid {grid_id}" + ) + + layout: List[dict] = [] + for d in derived: + if not d.tile_id: + layout.append({ + "tile_id": None, + "variable": d.name, + "url": d.asset_url, + "row_off": 0, + "col_off": 0, + "height": grid.rows, + "width": grid.cols, + }) + continue + + # Duas convenções de registro de tile, na mesma ordem que + # derive() usa para resolvê-las. + source = ( + self.catalog.get_spatial_source(f"{grid_id}_{d.tile_id}") + or self.catalog.get_spatial_source(d.tile_id) + ) + if source is None or not source.bbox: + raise ValueError( + f"Tile {d.tile_id!r} of {variable_id!r} has no SpatialSource " + f"with a bbox, so its position on grid {grid_id!r} is unknown. " + f"Looked for {grid_id}_{d.tile_id} and {d.tile_id}." + ) + + minx, miny, maxx, maxy = source.bbox + layout.append({ + "tile_id": d.tile_id, + "variable": d.name, + "url": d.asset_url, + "row_off": int(round((grid.bbox[3] - maxy) / grid.resolution)), + "col_off": int(round((minx - grid.bbox[0]) / grid.resolution)), + "height": int(round((maxy - miny) / grid.resolution)), + "width": int(round((maxx - minx) / grid.resolution)), + }) + + return sorted(layout, key=lambda t: (t["row_off"], t["col_off"])) + def to_lucc_data( self, variables: List[str], diff --git a/examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py b/examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py index 1fb136f..230402e 100644 --- a/examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py +++ b/examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py @@ -8,12 +8,22 @@ cubo em RAM via `to_lucc_data()` -> `RasterBackend`; este entrega em disco via workspace memmap, para domínios que não cabem na memória. -NÃO roda simulação — só deixa o workspace pronto e prova que os dados -chegaram corretos, inclusive nas costuras entre tiles. +NÃO roda simulação — só deixa o workspace pronto. + +Como as duas metades se encaixam +-------------------------------- +As variáveis foram derivadas tile a tile (`derive(tile_id=...)`), então +existe um Zarr por tile. `cube.tile_layout()` responde ONDE cada pedaço +cai na grade mestra — caminho e posição, dado puro — e o loader do +haloexec monta o workspace a partir disso, bloco a bloco. + +Nenhum dos dois pacotes importa o outro: o contrato entre eles é a lista +de dicionários, do mesmo jeito que o geomosaic entrega `tile_offsets` sem +conhecer quem vai lê-los. Por padrão usa a grade **BDC Albers** (`brmangue/30m_bdc`, do exemplo 02), -por ser a malha canônica: tiles BDC_SM compartilháveis entre projetos. -Troque `GRID_ID` para `brmangue/30m` para usar a malha ad hoc do 01. +por ser a malha canônica. Troque `GRID_ID` para `brmangue/30m` para usar a +malha ad hoc do 01. Pré-requisitos: - python examples/case_studies/brmangue_dominio/02_dominio_bdc.py @@ -21,179 +31,33 @@ Usage: python examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py - -Por que o carregamento é feito à mão aqui ------------------------------------------- -As variáveis foram derivadas TILE A TILE (`derive(tile_id=...)`), então -existe um Zarr por tile, não um só. Nenhum dos dois lados resolve isso hoje: - - * `CubeClient.load(nome, grid_id=...)` levanta ValueError quando há mais - de um tile — é preciso pedir um tile específico (comportamento - documentado em docs/architecture/tiling.md: "mosaico automático não - está implementado"); - * `haloexec.load_zarr_into_workspace()` recebe UM store e exige que o - shape dele seja igual ao do workspace inteiro — não há equivalente - multi-tile como `load_geotiffs_into_workspace` é para GeoTIFF. - -Então este exemplo faz a costura: descobre no catálogo quais tiles existem, -usa o bbox de cada um (registrado como SpatialSource) para achar sua posição -na grade mestra — a mesma ideia de `build_mosaic_contract` do geomosaic, mas -partindo da geometria catalogada em vez do nome do arquivo — e monta cada -bloco do workspace lendo dos tiles que o cobrem. - -Buracos na malha são esperados (o MapBiomas não exporta onde não há costa) -e viram nodata, não erro. """ from pathlib import Path -import numpy as np - from disscube.client import CubeClient -from disscube.models import GridSpec try: - import zarr - from haloexec import Block, MemmapRasterWorkspace + from haloexec import MemmapRasterWorkspace, load_zarr_tiles_into_workspace except ImportError as exc: # pragma: no cover — exemplo, fora do pacote raise SystemExit( 'Este exemplo precisa do haloexec: pip install "haloexec[zarr]"' ) from exc -# ── Configuração ────────────────────────────────────────────────────────────── -CATALOGO = "catalog.db" -STORE = "./data/" GRID_ID = "brmangue/30m_bdc" # malha BDC (exemplo 02); ou "brmangue/30m" VARIAVEIS = ["papel", "elevacao"] WS_DIR = Path("./data/workspace_brmangue_bdc") -# float32 em vez de float64: o workspace mantém DOIS slots (double buffer), +# float32 em vez de float64: o workspace mantém dois slots (double buffer), # então cada variável ocupa 2x o tamanho da grade em disco. DTYPE = "float32" BLOCK_H = BLOCK_W = 512 HALO = 2 -def _tiles_da_variavel(cube: CubeClient, nome: str, grid: GridSpec) -> list[dict]: - """Descobre os tiles de ``nome`` e onde cada um cai na grade mestra. - - A posição vem do bbox do SpatialSource do tile — o mesmo bbox que - `derive(tile_id=...)` usou para recortar — convertido para deslocamento - em pixel contra a origem da grade. O id do SpatialSource pode ser - ``{grid_id}_{tile_id}`` (malha própria) ou o próprio ``tile_id`` quando - já é um id canônico BDC. - """ - achados = [] - for d in cube.search(grid=grid.id): - if d.name != nome or not d.tile_id: - continue - src = (cube.catalog.get_spatial_source(f"{grid.id}_{d.tile_id}") - or cube.catalog.get_spatial_source(d.tile_id)) - if src is None or not src.bbox: - raise RuntimeError( - f"tile {d.tile_id} sem SpatialSource com bbox — não dá para " - "saber onde ele cai na grade mestra" - ) - arr = zarr.open(d.asset_url, mode="r")[nome] - achados.append({ - "tile_id": d.tile_id, - "url": d.asset_url, - "r0": int(round((grid.bbox[3] - src.bbox[3]) / grid.resolution)), - "c0": int(round((src.bbox[0] - grid.bbox[0]) / grid.resolution)), - "h": arr.shape[0], - "w": arr.shape[1], - }) - return sorted(achados, key=lambda t: (t["r0"], t["c0"])) - - -def _carregar(ws, nome: str, tiles: list[dict]) -> tuple[int, int]: - """Preenche ``nome`` no workspace, bloco a bloco, a partir de N tiles. - - Percorre os blocos do workspace (e não os tiles) porque um bloco pode - cair sobre dois tiles vizinhos, ou sobre um buraco da malha: montá-lo a - partir de tudo que o cobre é o que faz a costura ficar correta. - """ - abertos = {t["tile_id"]: zarr.open(t["url"], mode="r")[nome] for t in tiles} - com = sem = 0 - - for blk in ws.blocks(): - buf = np.full((blk.r1 - blk.r0, blk.c1 - blk.c0), np.nan, dtype=DTYPE) - tocou = False - for t in tiles: - r0 = max(blk.r0, t["r0"]); r1 = min(blk.r1, t["r0"] + t["h"]) - c0 = max(blk.c0, t["c0"]); c1 = min(blk.c1, t["c0"] + t["w"]) - if r0 >= r1 or c0 >= c1: - continue - buf[r0 - blk.r0:r1 - blk.r0, c0 - blk.c0:c1 - blk.c0] = np.asarray( - abertos[t["tile_id"]][r0 - t["r0"]:r1 - t["r0"], - c0 - t["c0"]:c1 - t["c0"]] - ) - tocou = True - ws.write_block_to_read_slot(blk, nome, buf) - com += tocou - sem += not tocou - - ws.flush() - return com, sem - - -def _provar_costura(ws, nome: str, esq: dict, dir_: dict) -> tuple[bool, str]: - """Prova que a janela com halo sobre a costura entre DOIS Zarr está certa. - - A comparação é feita contra cada arquivo Zarr de origem separadamente — - o lado esquerdo contra o Zarr de ``esq``, o direito contra o de ``dir_``, - célula a célula. Note o reindexamento: a primeira coluna do vizinho é a - coluna 0 DELE, não a coluna global — é exatamente aí que um erro de - offset apareceria. - - A linha é escolhida onde os dois lados DIFEREM. Sem isso a prova é fraca: - numa região toda 0.0 ou toda NaN a comparação passa sem provar nada. - """ - za = zarr.open(esq["url"], mode="r")[nome] - zb = zarr.open(dir_["url"], mode="r")[nome] - col = dir_["c0"] # coluna global onde o vizinho começa - - # linhas em que ambos os tiles existem - r_ini = max(esq["r0"], dir_["r0"]) + HALO + 1 - r_fim = min(esq["r0"] + esq["h"], dir_["r0"] + dir_["h"]) - HALO - 9 - if r_fim <= r_ini: - return True, "sem sobreposição vertical suficiente" - - ult = np.asarray(za[r_ini - esq["r0"]:r_fim - esq["r0"], esq["w"] - 1]) - pri = np.asarray(zb[r_ini - dir_["r0"]:r_fim - dir_["r0"], 0]) - dif = np.where(np.isfinite(ult) & np.isfinite(pri) & (ult != pri))[0] - if dif.size == 0: - return True, "nenhuma linha com lados diferentes (prova fraca, ignorada)" - - linha = r_ini + int(dif[len(dif) // 2]) - blk = Block(r0=linha - 4, r1=linha + 4, c0=col - 4, c1=col + 4) - jan = ws.read_block_with_halo(blk, boundary_value=np.nan)[nome] - - jr0, jc0 = blk.r0 - HALO, blk.c0 - HALO - for jr in range(jan.shape[0]): - for jc in range(jan.shape[1]): - gr, gc = jr0 + jr, jc0 + jc - if gc < col: - fonte, t = za, esq - else: - fonte, t = zb, dir_ - lr, lc = gr - t["r0"], gc - t["c0"] - if not (0 <= lr < t["h"] and 0 <= lc < t["w"]): - continue - esperado = np.float32(np.asarray(fonte[lr, lc])) - obtido = jan[jr, jc] - if not ((np.isnan(esperado) and np.isnan(obtido)) or esperado == obtido): - return False, (f"janela[{jr},{jc}] (global r={gr} c={gc}, tile " - f"{t['tile_id']} local r={lr} c={lc}): " - f"esperado {esperado}, obtido {obtido}") - - n = len(np.unique(jan[np.isfinite(jan)])) - return True, f"linha r={linha}, {n} valores distintos na janela" - - def main() -> None: - cube = CubeClient(catalog=CATALOGO, store=STORE) + cube = CubeClient(catalog="catalog.db", store="./data/") grid = cube.catalog.get_grid(GRID_ID) if grid is None: raise SystemExit( @@ -201,61 +65,29 @@ def main() -> None: " python examples/case_studies/brmangue_dominio/02_dominio_bdc.py" ) - print(f"\n[1/3] catálogo: {GRID_ID} = {grid.rows}x{grid.cols} @ {grid.crs[:40]}") - por_var = {} - for nome in VARIAVEIS: - tiles = _tiles_da_variavel(cube, nome, grid) - if not tiles: - raise SystemExit(f"nenhum tile derivado para {nome!r} em {GRID_ID}") - por_var[nome] = tiles - print(f" {nome}: {len(tiles)} tiles de {tiles[0]['h']}x{tiles[0]['w']}") - - try: - cube.load(VARIAVEIS[0], grid_id=GRID_ID) - print(" (load() sem tile_id funcionou — havia um tile só)") - except ValueError: - print(" load() sem tile_id recusa multi-tile, como esperado") + print(f"\n[1/2] {GRID_ID}: {grid.rows}x{grid.cols}") + layouts = {nome: cube.tile_layout(nome, GRID_ID) for nome in VARIAVEIS} + for nome, tiles in layouts.items(): + print(f" {nome}: {len(tiles)} tiles de " + f"{tiles[0]['height']}x{tiles[0]['width']}") - gib = grid.rows * grid.cols * np.dtype(DTYPE).itemsize * len(VARIAVEIS) * 2 / 1024**3 - print(f"\n[2/3] workspace {grid.rows}x{grid.cols} {DTYPE} " - f"({len(VARIAVEIS)} arrays x 2 slots = {gib:.1f} GB em disco)") + print(f"\n[2/2] montando o workspace em {WS_DIR}") ws = MemmapRasterWorkspace.create( WS_DIR, shape=(grid.rows, grid.cols), - arrays={n: DTYPE for n in VARIAVEIS}, + arrays={nome: DTYPE for nome in VARIAVEIS}, block_h=BLOCK_H, block_w=BLOCK_W, halo=HALO, ) - for nome, tiles in por_var.items(): - com, sem = _carregar(ws, nome, tiles) - print(f" {nome}: {com} blocos com dado, {sem} inteiramente vazios") - - # ── prova de costura entre arquivos Zarr distintos ────────────────────── - print("\n[3/3] costuras entre tiles (cada lado conferido contra o SEU Zarr)") - ws2 = MemmapRasterWorkspace(WS_DIR) - nome = VARIAVEIS[0] - tiles = por_var[nome] - ok = fracas = falhas = 0 - for a in tiles: - b = next((o for o in tiles - if o["r0"] == a["r0"] and o["c0"] == a["c0"] + a["w"]), None) - if b is None: - continue - passou, msg = _provar_costura(ws2, nome, a, b) - rotulo = f"{a['tile_id']} | {b['tile_id']}" - if not passou: - falhas += 1 - print(f" FALHA {rotulo}: {msg}") - elif "fraca" in msg or "sem sobreposição" in msg: - fracas += 1 - print(f" pulada {rotulo}: {msg}") - else: - ok += 1 - print(f" OK {rotulo}: {msg}") - print(f"\n {ok} costuras provadas, {fracas} sem contraste, {falhas} FALHAS") + for nome, tiles in layouts.items(): + load_zarr_tiles_into_workspace(ws, tiles, array=nome) + print(f" {nome}: carregado") - print(f"\n=== workspace pronto em {WS_DIR} ===") - print(f" shape={ws2.shape} arrays={list(ws2.metadata['arrays'])} " - f"blocos={len(ws2.blocks())} halo={ws2.halo}") + print(f"\n=== workspace pronto ===") + print(f" shape={ws.shape} arrays={list(ws.metadata['arrays'])} " + f"blocos={len(ws.blocks())} halo={ws.halo}") print(" (nenhum modelo foi executado — o workspace está pronto para receber um)") + print("\n Para rodar um modelo sobre ele, componha o mixin de disco:") + print(" class FloodModelDiskHalo(DiskChunkedSyncRasterModel, FloodModel): pass") + print(" FloodModelDiskHalo(workspace=ws, ...)") if __name__ == "__main__": diff --git a/tests/test_tile_layout.py b/tests/test_tile_layout.py new file mode 100644 index 0000000..3bbd3a8 --- /dev/null +++ b/tests/test_tile_layout.py @@ -0,0 +1,175 @@ +""" +Tests for ``CubeClient.tile_layout()``. + +O layout é o contrato entre o cubo e quem monta o dado — dado puro, sem +objeto de outro pacote. Estes testes fixam o formato e a aritmética de +posição, que é onde um erro passa despercebido: um offset trocado produz +um mosaico com tiles no lugar errado, sem erro nenhum. +""" + +import tempfile +from pathlib import Path + +import pytest + +from disscube.client import CubeClient +from disscube.models import DerivedVariable, GridSpec, SpatialSource + +CRS = "EPSG:31982" +RES = 10.0 +# grade 100x100 px: bbox de 1000x1000 unidades +GRID_BBOX = [0.0, 0.0, 1000.0, 1000.0] + + +@pytest.fixture +def cube(): + with tempfile.TemporaryDirectory() as d: + c = CubeClient(catalog=str(Path(d) / "c.db"), store=str(Path(d) / "s")) + c.register_grid(GridSpec( + id="G", type="local", crs=CRS, resolution=RES, bbox=list(GRID_BBOX), + )) + yield c + + +def _tile_source(cube, source_id, bbox): + cube.register_spatial_source(SpatialSource( + id=source_id, name=source_id, format="raster", + asset_url="planned", crs=CRS, bbox=list(bbox), + )) + + +def _derived(cube, name, tile_id, url="x.zarr"): + cube.catalog.save_derived(DerivedVariable( + id=f"{name}_{tile_id or 'global'}", name=name, grid_id="G", role="test", + times=[], dtype="float64", derivation_id="d", spec_hash="h", + tile_id=tile_id, asset_url=url, + )) + + +# ── formato do contrato ────────────────────────────────────────────────────── + +def test_layout_item_has_the_documented_keys(cube): + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + _derived(cube, "v", "T1") + item = cube.tile_layout("v", "G")[0] + assert set(item) == { + "tile_id", "variable", "url", "row_off", "col_off", "height", "width" + } + + +def test_variable_name_travels_with_the_tile(cube): + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + _derived(cube, "papel", "T1") + assert cube.tile_layout("papel", "G")[0]["variable"] == "papel" + + +# ── aritmética de posição ──────────────────────────────────────────────────── + +def test_top_left_tile_sits_at_origin(cube): + """bbox no canto superior-esquerdo da grade -> offset (0, 0).""" + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + _derived(cube, "v", "T1") + t = cube.tile_layout("v", "G")[0] + assert (t["row_off"], t["col_off"]) == (0, 0) + assert (t["height"], t["width"]) == (50, 50) + + +def test_row_offset_grows_downward_not_upward(cube): + """y cresce para cima no CRS, mas row cresce para baixo — o sinal aqui + é o erro clássico, e inverteria o mosaico verticalmente.""" + # tile na METADE DE BAIXO da grade: y de 0 a 500 + _tile_source(cube, "G_low", (0.0, 0.0, 500.0, 500.0)) + _derived(cube, "v", "low") + t = cube.tile_layout("v", "G")[0] + assert t["row_off"] == 50, "tile inferior deve começar na linha 50, não 0" + + +def test_column_offset_grows_rightward(cube): + _tile_source(cube, "G_right", (500.0, 500.0, 1000.0, 1000.0)) + _derived(cube, "v", "right") + t = cube.tile_layout("v", "G")[0] + assert (t["row_off"], t["col_off"]) == (0, 50) + + +def test_four_quadrants_tile_the_grid_without_gap_or_overlap(cube): + quadrantes = { + "NO": (0.0, 500.0, 500.0, 1000.0), + "NE": (500.0, 500.0, 1000.0, 1000.0), + "SO": (0.0, 0.0, 500.0, 500.0), + "SE": (500.0, 0.0, 1000.0, 500.0), + } + for tid, bbox in quadrantes.items(): + _tile_source(cube, f"G_{tid}", bbox) + _derived(cube, "v", tid) + + layout = cube.tile_layout("v", "G") + assert len(layout) == 4 + coberto = sum(t["height"] * t["width"] for t in layout) + assert coberto == 100 * 100, "os quatro quadrantes devem cobrir a grade" + cantos = {(t["row_off"], t["col_off"]) for t in layout} + assert cantos == {(0, 0), (0, 50), (50, 0), (50, 50)} + + +def test_layout_is_ordered_by_position(cube): + for tid, bbox in [ + ("SE", (500.0, 0.0, 1000.0, 500.0)), + ("NO", (0.0, 500.0, 500.0, 1000.0)), + ("NE", (500.0, 500.0, 1000.0, 1000.0)), + ]: + _tile_source(cube, f"G_{tid}", bbox) + _derived(cube, "v", tid) + ordem = [(t["row_off"], t["col_off"]) for t in cube.tile_layout("v", "G")] + assert ordem == sorted(ordem) + + +# ── convenções de registro de tile ─────────────────────────────────────────── + +def test_bdc_style_tile_id_resolves(cube): + """Tiles BDC são registrados pelo id canônico, sem prefixo de grade.""" + _tile_source(cube, "BDC_SM_009002", (0.0, 500.0, 500.0, 1000.0)) + _derived(cube, "v", "BDC_SM_009002") + assert cube.tile_layout("v", "G")[0]["tile_id"] == "BDC_SM_009002" + + +def test_grid_scoped_wins_over_bare_id(cube): + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + _tile_source(cube, "T1", (500.0, 0.0, 1000.0, 500.0)) + _derived(cube, "v", "T1") + t = cube.tile_layout("v", "G")[0] + assert (t["row_off"], t["col_off"]) == (0, 0) + + +# ── variável global (sem tiles) ────────────────────────────────────────────── + +def test_global_variable_yields_one_item_covering_the_grid(cube): + """Quem consome deve poder tratar global e tileado do mesmo jeito.""" + _derived(cube, "v", None) + layout = cube.tile_layout("v", "G") + assert len(layout) == 1 + t = layout[0] + assert t["tile_id"] is None + assert (t["row_off"], t["col_off"]) == (0, 0) + assert (t["height"], t["width"]) == (100, 100) + + +# ── erros ──────────────────────────────────────────────────────────────────── + +def test_unknown_grid_raises(cube): + with pytest.raises(ValueError, match="Grid not found"): + cube.tile_layout("v", "inexistente") + + +def test_variable_without_derivations_raises(cube): + with pytest.raises(ValueError, match="Derived variable not found"): + cube.tile_layout("nao_existe", "G") + + +def test_tile_without_bbox_raises_naming_what_was_searched(cube): + """Sem bbox não há posição — falhar alto é melhor que empilhar em (0,0).""" + cube.register_spatial_source(SpatialSource( + id="G_T1", name="T1", format="raster", asset_url="planned", crs=CRS, + )) + _derived(cube, "v", "T1") + with pytest.raises(ValueError) as exc: + cube.tile_layout("v", "G") + assert "G_T1" in str(exc.value) From 88c6357e4fba828bcf819c1c3f41af524a4963f8 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Fri, 28 Aug 2026 10:09:23 -0300 Subject: [PATCH 12/17] fix: tile_layout() conflated the time slices of a temporal variable MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A variable derived with valid_from/valid_until has one set of pieces PER SLICE, all at the same tile positions. tile_layout() returned them together, so every position appeared once per year — and assembling from that would let one year overwrite another with no error at all. Found by running the real thing: the MapBiomas mangrove series over five marker years came back as 140 pieces in 28 positions, each repeated five times with a different spec_hash in the URL. The existing tests missed it because they only covered static variables. tile_layout() now takes an optional `time`, and refuses a temporal variable when it is omitted, naming the available slices — the same discipline load() already applies to multi-tile. Each item also carries `times`, empty for static variables, so a caller can tell which slice it is holding. Verified on the series: without time it raises listing [1985, 1995, 2005, 2015, 2024]; with time= each slice yields 28 pieces in 28 distinct positions, no overlap. Static variables are unaffected and still need no time. Unrelated but worth recording, since it was the run that found this: the derived tiles were checked one by one against the source VRT window — 140/140 identical, and the 2024 count (7,527,925 mangrove cells) matches exactly what the structural domain reports for class 1, by an independent route. Tests: 7 new for the temporal cases. Full suite 134 passing. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- disscube/client/cube_client.py | 33 +++++++++++++++ tests/test_tile_layout.py | 77 +++++++++++++++++++++++++++++++++- 2 files changed, 109 insertions(+), 1 deletion(-) diff --git a/disscube/client/cube_client.py b/disscube/client/cube_client.py index 91a97ac..4205aa0 100644 --- a/disscube/client/cube_client.py +++ b/disscube/client/cube_client.py @@ -274,6 +274,7 @@ def tile_layout( self, variable_id: str, grid_id: str, + time: Optional[int] = None, ) -> List[dict]: """ Onde cada pedaço de uma variável derivada cai na grade mestra. @@ -300,6 +301,16 @@ def tile_layout( Nome da variável derivada (ex.: ``"papel"``). grid_id : str Grade mestra sobre a qual os pedaços se posicionam. + time : int, optional + Fatia temporal a devolver, para variáveis derivadas com janela + de validade (``valid_from``/``valid_until``). Uma variável + temporal tem um conjunto de pedaços POR ANO, todos nas mesmas + posições — devolvê-los juntos daria um layout em que cada + posição aparece várias vezes, e quem montasse a partir dele + sobrescreveria um ano com outro sem perceber. Por isso, se a + variável tiver mais de uma fatia e ``time`` não for informado, + o método falha em vez de escolher (mesma disciplina de + ``load()`` diante de multi-tile). Returns ------- @@ -313,6 +324,8 @@ def tile_layout( - ``col_off`` — coluna, em pixel - ``height`` — altura do pedaço, em pixel - ``width`` — largura do pedaço, em pixel + - ``times`` — anos cobertos por este pedaço (lista, vazia + quando a variável é estática) As dimensões vêm do ``bbox`` registrado, não do arquivo: é o mesmo ``bbox`` que ``derive(tile_id=...)`` usou para recortar, @@ -339,6 +352,24 @@ def tile_layout( f"Derived variable not found: {variable_id} on grid {grid_id}" ) + # Uma variável temporal repete cada posição uma vez por fatia; sem + # escolher a fatia, o layout descreveria a mesma célula N vezes. + fatias = sorted({t for d in derived for t in (d.times or [])}) + if time is not None: + derived = [d for d in derived if time in (d.times or [])] + if not derived: + raise ValueError( + f"Derived variable {variable_id!r} on grid {grid_id!r} has no " + f"slice for time {time}. Available: {fatias}" + ) + elif len(fatias) > 1: + raise ValueError( + f"Derived variable {variable_id!r} on grid {grid_id!r} is temporal " + f"and spans several slices: {fatias}. Each slice repeats the same " + f"tile positions, so a combined layout would describe every cell " + f"more than once — pass time= to pick one." + ) + layout: List[dict] = [] for d in derived: if not d.tile_id: @@ -350,6 +381,7 @@ def tile_layout( "col_off": 0, "height": grid.rows, "width": grid.cols, + "times": list(d.times or []), }) continue @@ -375,6 +407,7 @@ def tile_layout( "col_off": int(round((minx - grid.bbox[0]) / grid.resolution)), "height": int(round((maxy - miny) / grid.resolution)), "width": int(round((maxx - minx) / grid.resolution)), + "times": list(d.times or []), }) return sorted(layout, key=lambda t: (t["row_off"], t["col_off"])) diff --git a/tests/test_tile_layout.py b/tests/test_tile_layout.py index 3bbd3a8..a158401 100644 --- a/tests/test_tile_layout.py +++ b/tests/test_tile_layout.py @@ -53,7 +53,8 @@ def test_layout_item_has_the_documented_keys(cube): _derived(cube, "v", "T1") item = cube.tile_layout("v", "G")[0] assert set(item) == { - "tile_id", "variable", "url", "row_off", "col_off", "height", "width" + "tile_id", "variable", "url", "row_off", "col_off", "height", "width", + "times", } @@ -173,3 +174,77 @@ def test_tile_without_bbox_raises_naming_what_was_searched(cube): with pytest.raises(ValueError) as exc: cube.tile_layout("v", "G") assert "G_T1" in str(exc.value) + + +# ── variáveis temporais ────────────────────────────────────────────────────── +# Uma variável derivada com valid_from/valid_until tem um conjunto de pedaços +# POR FATIA, todos nas mesmas posições. Juntá-los daria um layout em que cada +# célula aparece N vezes — e quem montasse a partir dele sobrescreveria uma +# fatia com outra sem erro nenhum. Foi o que dado real revelou (série mangue +# 1985-2024): 140 pedaços em 28 posições. + +def _derived_t(cube, name, tile_id, times, url): + cube.catalog.save_derived(DerivedVariable( + id=f"{name}_{tile_id}_{times[0]}", name=name, grid_id="G", role="test", + times=times, dtype="float64", derivation_id="d", + spec_hash=f"h{times[0]}", tile_id=tile_id, asset_url=url, + )) + + +def test_temporal_variable_without_time_raises(cube): + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + for ano in (1985, 1995): + _derived_t(cube, "v", "T1", [ano], f"v_{ano}.zarr") + with pytest.raises(ValueError, match="temporal"): + cube.tile_layout("v", "G") + + +def test_error_lists_the_available_slices(cube): + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + for ano in (1985, 1995, 2005): + _derived_t(cube, "v", "T1", [ano], f"v_{ano}.zarr") + with pytest.raises(ValueError) as exc: + cube.tile_layout("v", "G") + assert "1985" in str(exc.value) and "2005" in str(exc.value) + + +def test_time_selects_one_slice(cube): + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + for ano in (1985, 1995): + _derived_t(cube, "v", "T1", [ano], f"v_{ano}.zarr") + layout = cube.tile_layout("v", "G", time=1995) + assert len(layout) == 1 + assert layout[0]["url"] == "v_1995.zarr" + assert layout[0]["times"] == [1995] + + +def test_each_slice_covers_every_position_exactly_once(cube): + """O ponto todo: com a fatia escolhida, nenhuma posição se repete.""" + for tid, bbox in [("A", (0.0, 500.0, 500.0, 1000.0)), + ("B", (500.0, 500.0, 1000.0, 1000.0))]: + _tile_source(cube, f"G_{tid}", bbox) + for ano in (1985, 1995): + _derived_t(cube, "v", tid, [ano], f"v_{tid}_{ano}.zarr") + layout = cube.tile_layout("v", "G", time=1985) + posicoes = [(t["row_off"], t["col_off"]) for t in layout] + assert len(posicoes) == 2 and len(set(posicoes)) == 2 + + +def test_unknown_time_raises_naming_what_exists(cube): + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + _derived_t(cube, "v", "T1", [1985], "v.zarr") + with pytest.raises(ValueError, match="1985"): + cube.tile_layout("v", "G", time=2020) + + +def test_single_slice_needs_no_time(cube): + """Uma variável com uma fatia só não é ambígua — não deve exigir time.""" + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + _derived_t(cube, "v", "T1", [1985], "v.zarr") + assert len(cube.tile_layout("v", "G")) == 1 + + +def test_static_variable_reports_empty_times(cube): + _tile_source(cube, "G_T1", (0.0, 500.0, 500.0, 1000.0)) + _derived(cube, "v", "T1") + assert cube.tile_layout("v", "G")[0]["times"] == [] From c39c32997e4bd3cf2de5069a4395ea92ffda5f83 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Fri, 28 Aug 2026 10:18:42 -0300 Subject: [PATCH 13/17] docs: add temporal mangrove series example, and call gc.collect() in it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Brings the MapBiomas binary mangrove presence into the cube as a TEMPORAL variable: one slice per year, so load() returns (time, y, x). The earlier examples only derive static variables, so nothing yet exercised the slice machinery — the SpatialSource time=, the valid_from/valid_until on the derivation, or tile_layout()'s per-slice selection. Source shape drives the structure: the mangue_* products hold one DECADE per file with one year per band, and write_vrt mosaics a single band per call, so each year becomes its own VRT. The loop calls gc.collect() after every derive, which is not incidental. Each derive() leaves ~50 objects in reference cycles holding ~380 MB. Reference counting cannot break a cycle; only the generational collector can, and it almost never fires here because its trigger is the NUMBER of allocations while numpy concentrates hundreds of MB in very few objects — the collector never sees the memory pressure. Measured on a cold run of all 140 derivations: without the collect, RSS grows ~380 MB per tile, passes 7 GB and the run was killed at tile 21 of year 2015; with it, peak RSS is 1.07 GB and the run finishes in 100.8s. Seven times less memory for about 7% more time (0.05s per collect against ~0.7s per derive). The call is in the example rather than inside derive() because imposing that cost on every caller is the package's decision, not this example's. The reasoning is written at the bottom of the file so it does not read as a stray defensive line. README updated with 03 and 04, which were both missing from the listing. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- examples/README.md | 5 + .../04_serie_temporal_mangue.py | 241 ++++++++++++++++++ 2 files changed, 246 insertions(+) create mode 100644 examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py diff --git a/examples/README.md b/examples/README.md index fef47c1..93f0350 100644 --- a/examples/README.md +++ b/examples/README.md @@ -26,6 +26,11 @@ Standalone — não depende dos scripts de `setup/`; registra a própria grade e a legenda de estados num raster de papéis, tile a tile, e escreve o GeoTIFF final. - `python examples/case_studies/brmangue_dominio/02_dominio_bdc.py` — o mesmo produto sobre a malha nacional **BDC_SM**, com a grade mestra em BDC Albers. +- `python examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py` — carrega + as variáveis derivadas num `MemmapRasterWorkspace` do haloexec (par em disco do + `maranhao/03_brmangue_simulate.py`, que entrega em RAM). Requer `haloexec[zarr]`. +- `python examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py` — série + **temporal** mangue/não-mangue: cada ano é uma fatia, `load()` devolve `(time, y, x)`. Demonstram três coisas de uso geral: - o operador genérico **`reclassify`** (tabela `{valor_origem: valor_destino}` como dado); diff --git a/examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py b/examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py new file mode 100644 index 0000000..2786d63 --- /dev/null +++ b/examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py @@ -0,0 +1,241 @@ +""" +examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py + +BR-MANGUE — série temporal binária mangue/não-mangue (MapBiomas). + +Traz para o cubo a presença de mangue ano a ano, como variável TEMPORAL: +cada ano é uma fatia, e `load()` devolve `(time, y, x)`. + +O que este exemplo mostra, além do 01/02 +---------------------------------------- +Os exemplos anteriores derivam variáveis estáticas. Aqui a mesma variável +existe em várias fatias de tempo, o que muda três coisas: + + * cada ano é um `SpatialSource` com `time=` e uma derivação com + `valid_from`/`valid_until` — é isso que faz o catálogo tratá-los como + fatias de uma série, e não como variáveis diferentes; + * `tile_layout()` exige `time=`, porque cada fatia repete as mesmas + posições de tile — sem escolher uma, o layout descreveria cada célula + várias vezes; + * o loop chama `gc.collect()` (ver nota no fim do arquivo). + +Estrutura da fonte +------------------ +Os produtos `mangue_*` do MapBiomas trazem uma DÉCADA por arquivo, com um +ano por banda (`mangue_1985` ... `mangue_1994`). Como `write_vrt` mosaica +uma banda por chamada, cada ano vira um VRT próprio. + +Pré-requisitos: + - pip install geomosaic + - tiles em $BRMANGUE_ENTRADA/mapbiomas_historico/ + +Usage: + export BRMANGUE_ENTRADA=/caminho/para/pymangue/dados/entrada + python examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py +""" + +import gc +import os +import time as _time +from pathlib import Path + +import numpy as np +import rasterio + +from disscube.client import CubeClient +from disscube.models import GridSpec, SpatialSource, SpatialDerivation, Variable + +try: + from geomosaic.core import build_mosaic_contract, write_vrt +except ImportError as exc: # pragma: no cover — exemplo, fora do pacote + raise SystemExit( + "Este exemplo precisa do geomosaic: pip install geomosaic" + ) from exc + + +# ── Configuração ────────────────────────────────────────────────────────────── +ENTRADA = Path(os.environ.get("BRMANGUE_ENTRADA", "./dados/entrada")) +MB_DIR = ENTRADA / "mapbiomas_historico" + +CATALOGO = "catalog.db" +STORE = "./data/" +VRT_DIR = Path("./data/raw/brmangue_serie") + +GRID_ID = "brmangue/30m" +VARIAVEL = "mangue" +TILE = 4096 + +# Cinco anos-marco, um por década. Trocar por range(1985, 2025) traz a série +# completa — 40 anos x ~28 tiles, o que leva bem mais tempo. +ANOS = [1985, 1995, 2005, 2015, 2024] + +# Em qual produto (década) e banda cada ano está. +DECADAS = [ + ("mangue_1985_1994", 1985), + ("mangue_1995_2004", 1995), + ("mangue_2005_2014", 2005), + ("mangue_2015_2024", 2015), +] + + +def _produto_e_banda(ano: int) -> tuple[str, int]: + """Década que contém o ano, e a banda dele dentro dela (1-based).""" + for produto, inicio in DECADAS: + if inicio <= ano <= inicio + 9: + return produto, ano - inicio + 1 + raise ValueError(f"Ano {ano} fora da cobertura {DECADAS[0][1]}–{DECADAS[-1][1]+9}") + + +def _mascara_valida(vrt: Path, largura_alvo: int = 2048): + """Onde a fonte tem dado, em baixa resolução — ver nota no exemplo 01.""" + from scipy.ndimage import binary_dilation + + with rasterio.open(vrt) as ds: + fator = max(1, int(np.ceil(ds.width / largura_alvo))) + baixa = ds.read( + 1, + out_shape=(max(1, ds.height // fator), max(1, ds.width // fator)), + resampling=rasterio.enums.Resampling.nearest, + ) + nodata = ds.nodata + + valida = np.isfinite(baixa) + if nodata is not None and np.isfinite(nodata): + valida &= baixa != nodata + return binary_dilation(valida), fator + + +def main() -> None: + if not MB_DIR.is_dir(): + raise SystemExit( + f"Diretório não encontrado: {MB_DIR}\n" + "Defina BRMANGUE_ENTRADA apontando para .../pymangue/dados/entrada" + ) + + t0 = _time.perf_counter() + VRT_DIR.mkdir(parents=True, exist_ok=True) + + # ── 1. geomosaic: um contrato por década, um VRT por ano ──────────────── + print(f"\n[1/4] geomosaic: {len(ANOS)} anos") + contratos: dict[str, object] = {} + vrts: dict[int, Path] = {} + for ano in ANOS: + produto, banda = _produto_e_banda(ano) + if produto not in contratos: + tiles = sorted(str(p) for p in MB_DIR.glob(f"*{produto}_serie*.tif")) + if not tiles: + raise SystemExit(f"Nenhum tile de {produto} em {MB_DIR}") + contratos[produto] = build_mosaic_contract(tiles) + c = contratos[produto] + print(f" {produto}: {len(tiles)} tiles -> " + f"{c.mosaic_height}x{c.mosaic_width}") + vrts[ano] = write_vrt( + contratos[produto], str(VRT_DIR / f"{VARIAVEL}_{ano}.vrt"), band=banda + ) + + formas = {(c.mosaic_height, c.mosaic_width) for c in contratos.values()} + if len(formas) > 1: + raise SystemExit(f"As décadas têm extensões diferentes: {formas}") + + # ── 2. grade e fontes ─────────────────────────────────────────────────── + print("\n[2/4] disscube: grade e fontes") + cube = CubeClient(catalog=CATALOGO, store=STORE) + contrato = next(iter(contratos.values())) + a, _b, ox, _d, e, oy = contrato.mosaic_transform + bbox = [ox, oy - contrato.mosaic_height * abs(e), + ox + contrato.mosaic_width * a, oy] + + grid = cube.catalog.get_grid(GRID_ID) + if grid is None: + grid = GridSpec(id=GRID_ID, type="reference", crs=str(contrato.crs), + resolution=a, bbox=bbox) + cube.register_grid(grid) + print(f" grade {GRID_ID} registrada: {grid.rows}x{grid.cols}") + else: + if not (np.allclose(grid.bbox, bbox) and grid.resolution == a): + raise SystemExit( + f"A grade {GRID_ID} já existe com outra extensão. " + "Use outro GRID_ID para esta série." + ) + print(f" grade {GRID_ID} reutilizada: {grid.rows}x{grid.cols}") + + # time= no SpatialSource é o que marca a fatia temporal. + for ano, vrt in vrts.items(): + cube.register_spatial_source(SpatialSource( + id=f"{VARIAVEL}_{ano}", name=f"MapBiomas mangue {ano}", + format="raster", asset_url=str(vrt), crs=str(contrato.crs), time=ano, + )) + print(f" {len(vrts)} fontes registradas, uma por ano") + + # ── 3. tiles com dado + derivação por ano ─────────────────────────────── + print("\n[3/4] derive por tile") + mascara, fator = _mascara_valida(vrts[ANOS[0]]) + tiles_com_dado = [] + for r0 in range(0, grid.rows, TILE): + for c0 in range(0, grid.cols, TILE): + h, w = min(TILE, grid.rows - r0), min(TILE, grid.cols - c0) + sub = mascara[r0 // fator:-(-(r0 + h) // fator), + c0 // fator:-(-(c0 + w) // fator)] + if not sub.any(): + continue + tile_id = f"R{r0:05d}C{c0:05d}" + minx, maxy = ox + c0 * a, oy - r0 * abs(e) + cube.register_spatial_source(SpatialSource( + id=f"{GRID_ID}_{tile_id}", name=tile_id, format="raster", + asset_url=str(vrts[ANOS[0]]), crs=str(contrato.crs), + bbox=[minx, maxy - h * abs(e), minx + w * a, maxy], + )) + tiles_com_dado.append(tile_id) + print(f" {len(tiles_com_dado)} tiles com dado válido") + + for ano in ANOS: + t = _time.perf_counter() + derivacao = SpatialDerivation( + source_id=f"{VARIAVEL}_{ano}", grid_id=GRID_ID, role="land_use", + variables=[Variable(name=VARIAVEL, operator="majority")], + valid_from=str(ano), valid_until=str(ano), + ) + for tile_id in tiles_com_dado: + cube.derive(derivacao, tile_id=tile_id) + # Ver a nota no fim do arquivo: sem isto, a memória cresce a cada + # derive e um loop longo acumula vários GB. + gc.collect() + print(f" {ano}: {len(tiles_com_dado)} tiles em " + f"{_time.perf_counter() - t:.1f}s") + + # ── 4. o que ficou no catálogo ────────────────────────────────────────── + print("\n[4/4] catálogo") + derivados = [d for d in cube.search(grid=GRID_ID) if d.name == VARIAVEL] + fatias = sorted({t for d in derivados for t in (d.times or [])}) + print(f" {len(derivados)} pedaços em {len(fatias)} fatias: {fatias}") + + # tile_layout() precisa da fatia: cada ano repete as mesmas posições. + for ano in ANOS: + layout = cube.tile_layout(VARIAVEL, GRID_ID, time=ano) + posicoes = {(t["row_off"], t["col_off"]) for t in layout} + print(f" {ano}: {len(layout)} pedaços em {len(posicoes)} posições") + + print(f"\n=== série no cubo em {_time.perf_counter() - t0:.1f}s ===") + print(f" Para carregar um ano num workspace do haloexec:") + print(f" tiles = cube.tile_layout({VARIAVEL!r}, {GRID_ID!r}, time={ANOS[0]})") + print(f" load_zarr_tiles_into_workspace(ws, tiles)") + print(f" Para o backend em RAM do DisSModel:") + print(f" cube.to_lucc_data([{VARIAVEL!r}], grid_id={GRID_ID!r})") + + +# ── Nota: por que gc.collect() no loop ──────────────────────────────────────── +# Cada derive() deixa ~50 objetos presos em ciclos de referência, retendo +# ~380 MB. Contagem de referências não desfaz ciclo — só o coletor geracional +# desfaz, e ele quase nunca dispara aqui: o gatilho é o NÚMERO de alocações, e +# o numpy concentra centenas de MB em pouquíssimos objetos, então o coletor não +# enxerga a pressão de memória. +# +# Medido em 28 tiles seguidos: sem collect, o RSS cresce ~380 MB por tile e +# passa de 7 GB; com collect, fica estável em ~0,35 GB. O custo é ~0,05 s por +# chamada contra ~0,7 s do derive — cerca de 7%. +# +# A chamada está aqui, e não dentro do disscube, porque impor esse custo a todo +# derive() é decisão do pacote, não deste exemplo. + +if __name__ == "__main__": + main() From 42d42d252a21949a8924f0e7a703227c859e9b65 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Fri, 28 Aug 2026 10:21:35 -0300 Subject: [PATCH 14/17] docs: record the derive() memory finding, and correct a stale limitation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two things, both found by running the pipeline on real data. The memory one is new and practical: each derive() leaves ~50 objects in reference cycles holding ~380 MB, and the generational collector almost never fires to reclaim them, because its trigger is the NUMBER of allocations while NumPy concentrates hundreds of MB in very few objects. A long loop therefore grows until it dies — measured at 140 tiles, RSS passed 7 GB and the run was killed; with gc.collect() per call it stayed near 1 GB, for about 7% more time. Documented under known limitations, with the workaround and a pointer to example 04, until the collection is done internally. The stale one: both the README and docs/architecture/tiling.md still said load() "silently returns the first result" for a multi-tile variable. It raises ValueError naming the tiles — the explicit error they described as planned is already there. Corrected, and tiling.md now documents tile_layout() as the way to actually consume a multi-tile variable, including the temporal case where a slice must be chosen. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- README.md | 11 +++++++++-- docs/architecture/tiling.md | 36 ++++++++++++++++++++++++++++++++---- 2 files changed, 41 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 53ea4b8..40fcd4a 100644 --- a/README.md +++ b/README.md @@ -196,8 +196,15 @@ Cada chamada a `derive()` carrega o dado completo de um tile em memória. Não h **Agregação vetorial por rasterização (não área-ponderada)** Operadores sobre fontes vetoriais (`majority`, `percentage`, `attribute`, `presence`, `minority`) convertem geometrias em raster antes de agregar pixels. A fração de cobertura de cada célula é estimada por contagem de pixels, não por cálculo de área de interseção. Para cobertura proporcional mais precisa, use uma fonte raster em resolução substancialmente maior que a célula-alvo. -**Desambiguação de tiles em `load()`** -`CubeClient.load(name)` sem `tile_id` retorna silenciosamente o primeiro resultado quando múltiplos tiles da mesma variável existem na mesma grade. Erro explícito ou mosaico automático estão planejados. **Especifique sempre `tile_id` em workloads multi-tile.** +**Sem mosaico automático em `load()`** +`CubeClient.load(name)` sem `tile_id` levanta `ValueError` quando múltiplos tiles da mesma variável existem na mesma grade — remontá-los num array único não está implementado. Para consumir uma variável multi-tile, use `tile_layout(name, grid_id)`: ele devolve onde cada pedaço cai na grade (caminho e posição, como dado puro) e deixa a montagem a cargo de quem consome. Para variáveis temporais é preciso escolher a fatia, com `time=`, já que cada fatia repete as mesmas posições de tile. + +**Memória acumula ao longo de um loop de `derive()`** +Cada `derive()` deixa cerca de 50 objetos presos em ciclos de referência, retendo aproximadamente 380 MB. Contagem de referências não desfaz ciclo, e o coletor geracional quase nunca dispara aqui: o gatilho dele é o *número* de alocações, enquanto o NumPy concentra centenas de MB em pouquíssimos objetos — o coletor não enxerga a pressão de memória. + +Na prática, um loop longo cresce até estourar. Medido ao derivar 140 tiles: o RSS sobe ~380 MB por tile e passa de 7 GB, sendo morto pelo sistema; com `gc.collect()` após cada chamada, fica estável em ~1 GB. O custo é cerca de 7% (0,05 s por coleta contra ~0,7 s por `derive()`). + +Até que a coleta seja feita internamente, **quem deriva muitos tiles ou muitas fatias em sequência deve chamar `gc.collect()` no loop** — ver `examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py`. Rodar cada tile em subprocesso separado também resolve, ao custo de reabrir o catálogo a cada vez. **`SpatialRelation` não atua no pipeline** O modelo `SpatialRelation` é persistido no catálogo, mas nenhum estágio do pipeline usa as relações durante a derivação — e por isso elas são **excluídas do `spec_hash`**. Incluí-las tornaria a chave de cache sensível a metadados que não afetam o resultado, quebrando a garantia de reprodutibilidade. A integração com estratégias hierárquicas de grades está reservada para versão futura. diff --git a/docs/architecture/tiling.md b/docs/architecture/tiling.md index 60b3997..4ebb70c 100644 --- a/docs/architecture/tiling.md +++ b/docs/architecture/tiling.md @@ -76,10 +76,38 @@ da = cube.load("dist_road", tile_id="009002") da = cube.load("dist_road", grid_id="BR/5km") ``` -> **Limitação atual:** `load()` sem `tile_id` retorna silenciosamente o primeiro resultado -> quando múltiplos tiles da mesma variável existem na mesma grade. A desambiguação -> automática (mosaico ou erro explícito) está planejada mas não implementada. -> Especifique sempre `tile_id` em workloads multi-tile. +> **Limitação atual:** `load()` sem `tile_id` levanta `ValueError` quando múltiplos +> tiles da mesma variável existem na mesma grade — remontá-los num array único não +> está implementado. + +## Consumir uma variável multi-tile + +`tile_layout()` responde o que `load()` não responde: onde cada pedaço cai na grade +mestra. Devolve dado puro — caminho e posição — sem montar nada e sem impor um +destino. + +```python +tiles = cube.tile_layout("dist_road", "BR/5km") +# [{"tile_id": "009002", "variable": "dist_road", "url": ".../dist_road.zarr", +# "row_off": 0, "col_off": 3520, "height": 3520, "width": 3520, "times": []}, ...] +``` + +Uma variável **temporal** tem um conjunto de pedaços por fatia, todos nas mesmas +posições; juntá-los descreveria cada célula várias vezes. Por isso `tile_layout()` +exige a fatia nesse caso: + +```python +tiles = cube.tile_layout("uso", "BR/5km", time=2015) +``` + +Sem `time`, uma variável com mais de uma fatia levanta `ValueError` nomeando as +disponíveis — a mesma disciplina que `load()` aplica a multi-tile. + +Uma variável sem tiles (partição `global`) devolve um único item cobrindo a grade +inteira, para que quem consome trate os dois casos igual. + +O destino fica a critério de quem chama. Para montar em disco, por exemplo, +`haloexec.load_zarr_tiles_into_workspace()` aceita essa lista diretamente. ## Vantagens From ba7e730bb7ce6640464892b86c768ce5f60f4615 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Fri, 28 Aug 2026 10:37:26 -0300 Subject: [PATCH 15/17] docs: export a multiband GeoTIFF from the temporal series MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit One band per year, written window by window from the Zarr tiles of each slice. It is not needed by anything in the cube — a multi-tile Zarr is not directly viewable, so this exists to open in QGIS and confirm the slices landed correctly. tools/zarr_to_tif.py does not cover this: it converts one whole store at a time, so it cannot reassemble N tiles, and it currently picks data_vars[0], which in a DisSCube Zarr is spatial_ref rather than the data. Verified against the same source the series was checked against: the five bands report 2.91% to 2.97% mangrove, matching what the Zarrs hold. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- .../04_serie_temporal_mangue.py | 47 ++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py b/examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py index 2786d63..e4b48b7 100644 --- a/examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py +++ b/examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py @@ -65,6 +65,12 @@ VARIAVEL = "mangue" TILE = 4096 +# Exportação para conferência: um GeoTIFF com um ano por banda. Não é +# necessário para nada no cubo — serve para abrir no QGIS e ver que as fatias +# estão certas, já que Zarr multi-tile não é diretamente visualizável. +SAIDA_TIF = Path("./data/brmangue_serie_mangue.tif") +NODATA_SAIDA = 255 + # Cinco anos-marco, um por década. Trocar por range(1985, 2025) traz a série # completa — 40 anos x ~28 tiles, o que leva bem mais tempo. ANOS = [1985, 1995, 2005, 2015, 2024] @@ -105,6 +111,40 @@ def _mascara_valida(vrt: Path, largura_alvo: int = 2048): return binary_dilation(valida), fator +def _exportar_tif(cube, grid, anos: list[int]) -> None: + """Escreve um GeoTIFF com uma banda por ano, a partir dos Zarr do cubo. + + Vai janela a janela, montando cada banda dos tiles daquela fatia — o Zarr + de uma variável multi-tile não é um arquivo só, então não há conversão + direta. (`tools/zarr_to_tif.py` converte UM store inteiro, o que não + cobre este caso.) + """ + import zarr + + perfil = dict( + driver="GTiff", height=grid.rows, width=grid.cols, count=len(anos), + dtype="uint8", nodata=NODATA_SAIDA, crs=grid.crs, + transform=grid.transform, + tiled=True, blockxsize=512, blockysize=512, + compress="deflate", zlevel=6, bigtiff="IF_SAFER", sparse_ok=True, + ) + SAIDA_TIF.parent.mkdir(parents=True, exist_ok=True) + + with rasterio.open(SAIDA_TIF, "w", **perfil) as dst: + for banda, ano in enumerate(anos, start=1): + dst.set_band_description(banda, f"{VARIAVEL}_{ano}") + for t in cube.tile_layout(VARIAVEL, GRID_ID, time=ano): + arr = np.asarray(zarr.open(t["url"], mode="r")[t["variable"]]) + dst.write( + np.where(np.isnan(arr), NODATA_SAIDA, arr).astype(np.uint8), + banda, + window=rasterio.windows.Window( + t["col_off"], t["row_off"], t["width"], t["height"] + ), + ) + print(f" banda {banda}: {VARIAVEL}_{ano}") + + def main() -> None: if not MB_DIR.is_dir(): raise SystemExit( @@ -204,7 +244,7 @@ def main() -> None: f"{_time.perf_counter() - t:.1f}s") # ── 4. o que ficou no catálogo ────────────────────────────────────────── - print("\n[4/4] catálogo") + print("\n[4/5] catálogo") derivados = [d for d in cube.search(grid=GRID_ID) if d.name == VARIAVEL] fatias = sorted({t for d in derivados for t in (d.times or [])}) print(f" {len(derivados)} pedaços em {len(fatias)} fatias: {fatias}") @@ -215,7 +255,12 @@ def main() -> None: posicoes = {(t["row_off"], t["col_off"]) for t in layout} print(f" {ano}: {len(layout)} pedaços em {len(posicoes)} posições") + # ── 5. exportação para conferência visual ─────────────────────────────── + print(f"\n[5/5] GeoTIFF multibanda -> {SAIDA_TIF}") + _exportar_tif(cube, grid, ANOS) + print(f"\n=== série no cubo em {_time.perf_counter() - t0:.1f}s ===") + print(f" {SAIDA_TIF} — {len(ANOS)} bandas, uma por ano, para abrir no QGIS") print(f" Para carregar um ano num workspace do haloexec:") print(f" tiles = cube.tile_layout({VARIAVEL!r}, {GRID_ID!r}, time={ANOS[0]})") print(f" load_zarr_tiles_into_workspace(ws, tiles)") From 2a5607164d39b03131bf8cd75327ca1c7ffce520 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Fri, 28 Aug 2026 10:58:21 -0300 Subject: [PATCH 16/17] docs: add the temporal series on the BDC mesh, with a windowed TIF export MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Combines what 04 and 02 do separately: the temporal mangrove series, but on the national BDC_SM mesh with the master grid in BDC Albers, so every tile goes through a real 5880 -> Albers reprojection rather than the identity path. 120 derivations, 170s cold, peaking at 1.47 GB — the gc.collect() from 04 is what keeps that bounded. The export deserves a note, because a first version of it was wrong in a way that produced invented data. Writing tile by tile leaves the uncovered part of each border block filled with GDAL's default zero rather than the declared nodata, and here zero is a VALID value (non-mangrove). The result had 19.6 million cells outside the study area reading as observed non-mangrove, in two of the five bands. The cause is that BDC tiles are 3520x3520 while the GeoTIFF blocks are 512, and 3520 is not a multiple of 512, so every tile edge lands mid-block. Examples 01 and 02 never hit this because their tiles are 4096, an exact multiple. The fix is the same lesson the Zarr loader already encodes: walk the blocks of the DESTINATION, not the tiles, and build each block from whatever covers it. Verified after the fix: all five bands report 206,007,691 valid cells, matching the Zarrs exactly, with zero cells outside the tiled area. Also records what the reprojection costs — class counts land ~1% below the EPSG:5880 series of example 04, which does not resample. That is the price of a canonical mesh, not an error. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- examples/README.md | 2 + .../brmangue_dominio/05_serie_temporal_bdc.py | 327 ++++++++++++++++++ 2 files changed, 329 insertions(+) create mode 100644 examples/case_studies/brmangue_dominio/05_serie_temporal_bdc.py diff --git a/examples/README.md b/examples/README.md index 93f0350..15b13a3 100644 --- a/examples/README.md +++ b/examples/README.md @@ -31,6 +31,8 @@ Standalone — não depende dos scripts de `setup/`; registra a própria grade e `maranhao/03_brmangue_simulate.py`, que entrega em RAM). Requer `haloexec[zarr]`. - `python examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py` — série **temporal** mangue/não-mangue: cada ano é uma fatia, `load()` devolve `(time, y, x)`. +- `python examples/case_studies/brmangue_dominio/05_serie_temporal_bdc.py` — a mesma + série sobre a malha **BDC_SM**, com reprojeção 5880 → Albers. Requer `disscube[bdc]`. Demonstram três coisas de uso geral: - o operador genérico **`reclassify`** (tabela `{valor_origem: valor_destino}` como dado); diff --git a/examples/case_studies/brmangue_dominio/05_serie_temporal_bdc.py b/examples/case_studies/brmangue_dominio/05_serie_temporal_bdc.py new file mode 100644 index 0000000..ea29ec2 --- /dev/null +++ b/examples/case_studies/brmangue_dominio/05_serie_temporal_bdc.py @@ -0,0 +1,327 @@ +""" +examples/case_studies/brmangue_dominio/05_serie_temporal_bdc.py + +BR-MANGUE — a série temporal mangue/não-mangue sobre a malha BDC (Albers). + +Combina o que os dois exemplos anteriores fazem separadamente: + + 04_ série temporal, malha ad hoc, grade EPSG:5880 (a da fonte) + 05_ série temporal, malha BDC_SM, grade em BDC Albers ← este + +A diferença não é cosmética. Em `04` a fonte já está na grade alvo, então o +GridAligner usa o caminho de identidade — leitura janelada, sem reamostrar. +Aqui cada tile passa por uma reprojeção real (EPSG:5880 → BDC Albers), o que +custa tempo e memória, e reamostra: o resultado não é bit a bit igual à fonte. + +Quando usar cada um +------------------- +Use a malha BDC quando os derivados forem virar patrimônio do cubo — os tiles +têm significado fora do script e alinham com o resto do ecossistema BDC. Use a +ad hoc quando forem passo intermediário: é mais leve e preserva a fonte. + +Pré-requisitos: + - pip install geomosaic + - pip install "disscube[bdc]" (fiona, para ler os shapefiles BDC) + - grades BDC em data/bdc_grids/BDC_{SM,MD,LG}_V2.zip + - tiles em $BRMANGUE_ENTRADA/mapbiomas_historico/ + +Usage: + export BRMANGUE_ENTRADA=/caminho/para/pymangue/dados/entrada + python examples/case_studies/brmangue_dominio/05_serie_temporal_bdc.py +""" + +import gc +import os +import time as _time +from pathlib import Path + +import numpy as np +import rasterio +from pyproj import CRS as ProjCRS, Transformer +from shapely.geometry import MultiPoint, shape + +from disscube.client import CubeClient +from disscube.models import GridSpec, SpatialSource, SpatialDerivation, Variable +from disscube.utils.grids import BDC_CRS + +try: + from geomosaic.core import build_mosaic_contract, write_vrt +except ImportError as exc: # pragma: no cover — exemplo, fora do pacote + raise SystemExit( + "Este exemplo precisa do geomosaic: pip install geomosaic" + ) from exc + +try: + import fiona +except ImportError as exc: # pragma: no cover + raise SystemExit( + 'Este exemplo precisa do fiona: pip install "disscube[bdc]"' + ) from exc + + +# ── Configuração ────────────────────────────────────────────────────────────── +ENTRADA = Path(os.environ.get("BRMANGUE_ENTRADA", "./dados/entrada")) +MB_DIR = ENTRADA / "mapbiomas_historico" +BDC_SM = "zip://data/bdc_grids/BDC_SM_V2.zip" + +CATALOGO = "catalog.db" +STORE = "./data/" +VRT_DIR = Path("./data/raw/brmangue_serie") +SAIDA_TIF = Path("./data/brmangue_serie_mangue_bdc.tif") + +GRID_ID = "brmangue/30m_bdc" +VARIAVEL = "mangue" +RESOLUCAO = 30.0 +NODATA_SAIDA = 255 + +ANOS = [1985, 1995, 2005, 2015, 2024] + +DECADAS = [ + ("mangue_1985_1994", 1985), + ("mangue_1995_2004", 1995), + ("mangue_2005_2014", 2005), + ("mangue_2015_2024", 2015), +] + + +def _produto_e_banda(ano: int) -> tuple[str, int]: + for produto, inicio in DECADAS: + if inicio <= ano <= inicio + 9: + return produto, ano - inicio + 1 + raise ValueError(f"Ano {ano} fora da cobertura") + + +def _mascara_valida(vrt: Path, largura_alvo: int = 2048): + from scipy.ndimage import binary_dilation + + with rasterio.open(vrt) as ds: + fator = max(1, int(np.ceil(ds.width / largura_alvo))) + baixa = ds.read( + 1, + out_shape=(max(1, ds.height // fator), max(1, ds.width // fator)), + resampling=rasterio.enums.Resampling.nearest, + ) + nodata, transform = ds.nodata, ds.transform + resolucao = abs(transform.a) * (ds.width / baixa.shape[1]) + + valida = np.isfinite(baixa) + if nodata is not None and np.isfinite(nodata): + valida &= baixa != nodata + return binary_dilation(valida), transform, resolucao + + +def _tem_dado(bounds, mascara, transform, resolucao, para_fonte) -> bool: + minx, miny, maxx, maxy = bounds + xs, ys = para_fonte.transform([minx, minx, maxx, maxx], [miny, maxy, miny, maxy]) + c0 = int(np.floor((min(xs) - transform.c) / resolucao)) + c1 = int(np.ceil((max(xs) - transform.c) / resolucao)) + r0 = int(np.floor((transform.f - max(ys)) / resolucao)) + r1 = int(np.ceil((transform.f - min(ys)) / resolucao)) + h, w = mascara.shape + r0, r1 = max(0, r0), min(h, r1) + c0, c1 = max(0, c0), min(w, c1) + if r0 >= r1 or c0 >= c1: + return False + return bool(mascara[r0:r1, c0:c1].any()) + + +BLOCO_TIF = 512 + + +def _exportar_tif(cube, grid, anos: list[int]) -> None: + """GeoTIFF com uma banda por ano, montado dos Zarr — conferência visual. + + Percorre os BLOCOS do GeoTIFF, não os tiles, e é preciso que seja assim. + Os tiles BDC são 3520x3520, que não é múltiplo do bloco de 512, então cada + borda de tile cai no meio de um bloco. Escrevendo tile a tile, a parte não + coberta desses blocos de borda fica com o preenchimento padrão do GDAL — + zero, não o nodata declarado — e zero é um valor VÁLIDO aqui (não-mangue). + O resultado seriam milhões de células fora da área de estudo aparecendo + como "não-mangue" observado, o que é dado inventado. + + Montando cada bloco a partir dos tiles que o cobrem, com nodata no resto, + o problema não existe. (Nos exemplos 01/02 os tiles são 4096, múltiplo de + 512, e por isso a questão nunca apareceu lá.) + """ + import zarr + + perfil = dict( + driver="GTiff", height=grid.rows, width=grid.cols, count=len(anos), + dtype="uint8", nodata=NODATA_SAIDA, crs=grid.crs, transform=grid.transform, + tiled=True, blockxsize=BLOCO_TIF, blockysize=BLOCO_TIF, + compress="deflate", zlevel=6, bigtiff="IF_SAFER", sparse_ok=True, + ) + SAIDA_TIF.parent.mkdir(parents=True, exist_ok=True) + + with rasterio.open(SAIDA_TIF, "w", **perfil) as dst: + for banda, ano in enumerate(anos, start=1): + dst.set_band_description(banda, f"{VARIAVEL}_{ano}") + tiles = cube.tile_layout(VARIAVEL, GRID_ID, time=ano) + abertos = {t["tile_id"]: zarr.open(t["url"], mode="r")[t["variable"]] + for t in tiles} + + for r0 in range(0, grid.rows, BLOCO_TIF): + r1 = min(r0 + BLOCO_TIF, grid.rows) + for c0 in range(0, grid.cols, BLOCO_TIF): + c1 = min(c0 + BLOCO_TIF, grid.cols) + buf = None + for t in tiles: + tr0 = max(r0, t["row_off"]); tr1 = min(r1, t["row_off"] + t["height"]) + tc0 = max(c0, t["col_off"]); tc1 = min(c1, t["col_off"] + t["width"]) + if tr0 >= tr1 or tc0 >= tc1: + continue + if buf is None: + buf = np.full((r1 - r0, c1 - c0), NODATA_SAIDA, np.uint8) + trecho = np.asarray(abertos[t["tile_id"]][ + tr0 - t["row_off"]:tr1 - t["row_off"], + tc0 - t["col_off"]:tc1 - t["col_off"], + ]) + buf[tr0 - r0:tr1 - r0, tc0 - c0:tc1 - c0] = np.where( + np.isnan(trecho), NODATA_SAIDA, trecho + ).astype(np.uint8) + if buf is not None: + dst.write(buf, banda, + window=rasterio.windows.Window(c0, r0, c1 - c0, r1 - r0)) + print(f" banda {banda}: {VARIAVEL}_{ano}") + + +def main() -> None: + if not MB_DIR.is_dir(): + raise SystemExit( + f"Diretório não encontrado: {MB_DIR}\n" + "Defina BRMANGUE_ENTRADA apontando para .../pymangue/dados/entrada" + ) + + t0 = _time.perf_counter() + VRT_DIR.mkdir(parents=True, exist_ok=True) + + # ── 1. geomosaic: um VRT por ano, na projeção nativa da fonte ─────────── + print(f"\n[1/5] geomosaic: {len(ANOS)} anos") + contratos: dict[str, object] = {} + vrts: dict[int, Path] = {} + for ano in ANOS: + produto, banda = _produto_e_banda(ano) + if produto not in contratos: + tiles = sorted(str(p) for p in MB_DIR.glob(f"*{produto}_serie*.tif")) + if not tiles: + raise SystemExit(f"Nenhum tile de {produto} em {MB_DIR}") + contratos[produto] = build_mosaic_contract(tiles) + vrts[ano] = write_vrt( + contratos[produto], str(VRT_DIR / f"{VARIAVEL}_{ano}.vrt"), band=banda + ) + contrato = next(iter(contratos.values())) + print(f" fonte: {contrato.mosaic_height}x{contrato.mosaic_width} " + f"@ {contrato.crs}") + + # ── 2. selecionar os tiles BDC com dado ───────────────────────────────── + print("[2/5] BDC: selecionando tiles SM com dado válido") + a, _b, ox, _d, e, oy = contrato.mosaic_transform + fonte_bbox = (ox, oy - contrato.mosaic_height * abs(e), + ox + contrato.mosaic_width * a, oy) + + para_bdc = Transformer.from_crs( + ProjCRS.from_user_input(str(contrato.crs)), + ProjCRS.from_user_input(BDC_CRS), always_xy=True, + ) + # Footprint densificado antes de projetar: só os quatro cantos dariam uma + # envoltória maior que a área real, porque as bordas viram curvas. + passo, contorno = 200, [] + x0, y0, x1, y1 = fonte_bbox + for i in range(passo + 1): + f = i / passo + contorno += [(x0 + (x1 - x0) * f, y0), (x0 + (x1 - x0) * f, y1), + (x0, y0 + (y1 - y0) * f), (x1, y0 + (y1 - y0) * f)] + px, py = para_bdc.transform([p[0] for p in contorno], [p[1] for p in contorno]) + alvo = MultiPoint(list(zip(px, py))).convex_hull + + candidatos = [] + with fiona.open(BDC_SM) as src: + for rec in src: + geom = shape(rec["geometry"]) + if geom.intersects(alvo): + candidatos.append((rec["properties"]["tile"], geom.bounds)) + candidatos.sort() + + mascara, m_transform, m_res = _mascara_valida(vrts[ANOS[0]]) + para_fonte = Transformer.from_crs( + ProjCRS.from_user_input(BDC_CRS), + ProjCRS.from_user_input(str(contrato.crs)), always_xy=True, + ) + selecionados = [ + (tile, bounds) for tile, bounds in candidatos + if _tem_dado(bounds, mascara, m_transform, m_res, para_fonte) + ] + if not selecionados: + raise SystemExit("Nenhum tile BDC_SM contém dado válido.") + print(f" {len(candidatos)} intersectam a extensão, " + f"{len(selecionados)} têm dado válido") + + out_minx = min(b[0] for _t, b in selecionados) + out_miny = min(b[1] for _t, b in selecionados) + out_maxx = max(b[2] for _t, b in selecionados) + out_maxy = max(b[3] for _t, b in selecionados) + + # ── 3. grade BDC + fontes ─────────────────────────────────────────────── + print("[3/5] disscube: grade Albers, fontes e tiles") + cube = CubeClient(catalog=CATALOGO, store=STORE) + bbox = [out_minx, out_miny, out_maxx, out_maxy] + + grid = cube.catalog.get_grid(GRID_ID) + if grid is None: + grid = GridSpec(id=GRID_ID, type="reference", crs=BDC_CRS, + resolution=RESOLUCAO, bbox=bbox, + description="BR-MANGUE 30 m sobre a malha BDC (Albers)") + cube.register_grid(grid) + print(f" grade {GRID_ID} registrada: {grid.rows}x{grid.cols}") + else: + if not np.allclose(grid.bbox, bbox): + raise SystemExit( + f"A grade {GRID_ID} já existe com outra extensão. " + "Use outro GRID_ID para esta série." + ) + print(f" grade {GRID_ID} reutilizada: {grid.rows}x{grid.cols}") + + for ano, vrt in vrts.items(): + cube.register_spatial_source(SpatialSource( + id=f"{VARIAVEL}_{ano}", name=f"MapBiomas mangue {ano}", + format="raster", asset_url=str(vrt), crs=str(contrato.crs), time=ano, + )) + for tile, bounds in selecionados: + cube.register_spatial_source(SpatialSource( + id=f"BDC_SM_{tile}", name=f"BDC SM Tile {tile}", format="raster", + asset_url="planned", crs=BDC_CRS, bbox=list(bounds), + )) + print(f" {len(vrts)} fontes + {len(selecionados)} tiles BDC") + + # ── 4. derive por tile BDC, por ano ───────────────────────────────────── + print("[4/5] derive por tile (COM reprojeção 5880 -> Albers)") + for ano in ANOS: + t = _time.perf_counter() + derivacao = SpatialDerivation( + source_id=f"{VARIAVEL}_{ano}", grid_id=GRID_ID, role="land_use", + variables=[Variable(name=VARIAVEL, operator="majority")], + valid_from=str(ano), valid_until=str(ano), + ) + for tile, _bounds in selecionados: + # id completo: ids simples não são únicos entre níveis BDC. + cube.derive(derivacao, tile_id=f"BDC_SM_{tile}") + # Sem isto a memória cresce a cada derive — ver a nota no exemplo 04. + gc.collect() + print(f" {ano}: {len(selecionados)} tiles em " + f"{_time.perf_counter() - t:.1f}s") + + # ── 5. exportação para conferência ────────────────────────────────────── + print(f"[5/5] GeoTIFF multibanda -> {SAIDA_TIF}") + _exportar_tif(cube, grid, ANOS) + + derivados = [d for d in cube.search(grid=GRID_ID) if d.name == VARIAVEL] + print(f"\n=== série na malha BDC em {_time.perf_counter() - t0:.1f}s ===") + print(f" {len(derivados)} pedaços, {len(ANOS)} fatias, grade {grid.rows}x{grid.cols}") + print(f" {SAIDA_TIF} — {len(ANOS)} bandas em Albers") + print(f"\n A reprojeção reamostra: as contagens por classe ficam ~1% abaixo") + print(f" das do exemplo 04 (EPSG:5880), que não reamostra. É o custo de") + print(f" estar numa malha canônica, não um erro.") + + +if __name__ == "__main__": + main() From 04f0fa030c5b049c24b9a4d02e8bf21cebeaf67a Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Fri, 28 Aug 2026 11:23:03 -0300 Subject: [PATCH 17/17] feat: add tools/tiles_to_tif.py to inspect a tile mesh MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Checking whether tiles landed in the right place meant opening the whole mosaic, and where the values on both sides of a boundary agree, that shows nothing at all. This crops only the requested tiles to their own envelope and adds a band carrying the tile index, which makes the boundaries visible even there. Styled as categorical in QGIS, the seams show up directly. It complements tools/zarr_to_tif.py rather than replacing it: that one converts a single store, this one assembles several tiles — which is the case a multi-tile variable actually presents, and the one zarr_to_tif cannot cover. Writes by walking the blocks of the GeoTIFF, not the tiles. When the tile side is not a multiple of the block (BDC tiles are 3520, blocks are 512), every tile edge lands mid-block and the uncovered part of those blocks keeps GDAL's default zero rather than the declared nodata. Where zero is a valid value that is invented data — the same trap example 05 fell into before being fixed. Resolves the catalog's relative asset_url against the catalog's own directory, so the tool works when run from anywhere, not only from the repository root. Errors speak the tool's language rather than the API's: an unknown tile is refused listing what exists, and a temporal variable without --tempos names the available slices and shows the flag to use. README documents it with both modes and explains the block-walking, so the reason survives longer than this session. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB --- examples/README.md | 31 +++++- tools/tiles_to_tif.py | 218 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 248 insertions(+), 1 deletion(-) create mode 100644 tools/tiles_to_tif.py diff --git a/examples/README.md b/examples/README.md index 15b13a3..ab67ee9 100644 --- a/examples/README.md +++ b/examples/README.md @@ -76,10 +76,39 @@ apontando para o diretório com os tiles ANADEM v2. O `02_` também precisa de | Script | Uso | |---|---| -| `tools/zarr_to_tif.py` | Converte Zarr derivado para GeoTIFF | +| `tools/zarr_to_tif.py` | Converte **um** Zarr derivado para GeoTIFF | +| `tools/tiles_to_tif.py` | Recorta **alguns tiles** de uma variável num GeoTIFF, para inspecionar a malha | | `tools/import_bdc_tiles.py` | Importa tiles BDC SM/MD/LG no catálogo (one-time, lento) | ```bash python tools/zarr_to_tif.py data/derived/.../var.zarr output.tif python tools/import_bdc_tiles.py ``` + +### Inspecionar a malha de tiles + +Abrir o mosaico inteiro para conferir se os tiles estão no lugar certo é +desconfortável e, quando os valores dos dois lados de uma fronteira coincidem, +não mostra nada. `tools/tiles_to_tif.py` recorta só os tiles pedidos e +acrescenta uma banda com o **índice do tile**, que torna as fronteiras visíveis +mesmo aí: + +```bash +# três tiles vizinhos, duas fatias temporais +python tools/tiles_to_tif.py --grid brmangue/30m_bdc --variavel mangue \ + --tiles 029006 030006 030007 --tempos 1985 2024 + +# variável estática: basta omitir --tempos +python tools/tiles_to_tif.py --grid BR/5km --variavel slope --tiles 009002 +``` + +Saída: uma banda por fatia pedida, mais `indice_do_tile` (1, 2, 3… na ordem em +que foram pedidos). No QGIS, estilize essa última como categórica e os limites +aparecem. Um tile inexistente é recusado listando os disponíveis, e uma variável +temporal sem `--tempos` avisa quais fatias existem. + +> **Por que ele monta por blocos, e não tile a tile:** quando o lado do tile não +> é múltiplo do bloco do GeoTIFF (os tiles BDC são 3520, o bloco é 512), cada +> borda de tile cai no meio de um bloco, e a parte não coberta fica com o +> preenchimento padrão do GDAL — zero, não o nodata declarado. Onde zero é um +> valor válido, isso vira dado inventado. diff --git a/tools/tiles_to_tif.py b/tools/tiles_to_tif.py new file mode 100644 index 0000000..0ae8abf --- /dev/null +++ b/tools/tiles_to_tif.py @@ -0,0 +1,218 @@ +#!/usr/bin/env python3 +""" +tools/tiles_to_tif.py + +Gera um GeoTIFF só de alguns tiles, recortado à envoltória deles. + +Serve para inspecionar a malha: em vez de abrir o mosaico inteiro, sai um +arquivo pequeno cobrindo exatamente os tiles pedidos, com uma banda extra +marcando de qual tile veio cada pixel — assim as fronteiras ficam visíveis +no QGIS mesmo onde os valores dos dois lados coincidem. + +Bandas de saída: + 1..N uma por fatia temporal pedida (ou uma só, se a variável é estática) + N+1 índice do tile (1, 2, 3... na ordem em que foram pedidos; 255 = fora) + +A banda de índice é o que torna as fronteiras visíveis: onde os valores dos +dois lados coincidem, só ela distingue de qual tile veio cada pixel. + +Diferente de `tools/zarr_to_tif.py`, que converte UM store inteiro, aqui a +saída é montada de vários tiles — e percorrendo os blocos do GeoTIFF, não os +tiles: quando o lado do tile não é múltiplo do bloco (tiles BDC são 3520, +blocos são 512), escrever tile a tile deixa a parte não coberta dos blocos de +borda com zero em vez do nodata declarado. + +Uso: + python tools/tiles_to_tif.py --tiles 029006 030006 030007 + python tools/tiles_to_tif.py --tiles 029006 --anos 1985 2024 + python tools/tiles_to_tif.py --grid brmangue/30m --tiles R00000C00000 \ + --variavel papel --anos 0 +""" +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import numpy as np +import rasterio +import zarr +from rasterio.transform import Affine +from rasterio.windows import Window + +from disscube.client import CubeClient + +# Raiz do repositório — tools/ fica um nível abaixo. +REPO = Path(__file__).resolve().parents[1] +NODATA = 255 +BLOCO = 512 + + +def _resolver(url: str, raiz: Path) -> str: + """O catálogo guarda asset_url RELATIVO ao diretório de onde o disscube foi + usado. Rodando de outro lugar, é preciso resolver contra a raiz do repo — + senão o caminho não existe e a leitura falha.""" + p = Path(url) + return str(p if p.is_absolute() else (raiz / url).resolve()) + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--catalogo", default=str(REPO / "catalog.db")) + ap.add_argument("--store", default=str(REPO / "data")) + ap.add_argument("--grid", required=True) + ap.add_argument("--variavel", required=True) + ap.add_argument("--tiles", nargs="+", required=True, + help="ids dos tiles; para tiles BDC o prefixo BDC_SM_ é opcional") + ap.add_argument("--tempos", nargs="+", type=int, default=None, + help="fatias temporais a incluir, uma banda cada; " + "omita para variável estática") + ap.add_argument("--saida", default=None) + args = ap.parse_args() + + cube = CubeClient(catalog=args.catalogo, store=args.store) + grid = cube.catalog.get_grid(args.grid) + if grid is None: + raise SystemExit(f"Grade {args.grid!r} não encontrada em {args.catalogo}") + + # ── localizar os tiles pedidos ─────────────────────────────────────────── + pedidos, vistos = [], set() + for t in args.tiles: + if t in vistos: + print(f" aviso: tile {t!r} repetido no pedido — usando uma vez só") + continue + vistos.add(t) + pedidos.append(t) + + tempos = args.tempos + try: + layout = cube.tile_layout(args.variavel, args.grid, + time=tempos[0] if tempos else None) + except ValueError as erro: + # tile_layout fala na linguagem da API (time=); aqui o usuário + # tem um flag, então vale traduzir em vez de repassar o traceback. + if "temporal" in str(erro): + fatias = sorted({t for d in cube.search(grid=args.grid) + if d.name == args.variavel for t in (d.times or [])}) + raise SystemExit( + f"{args.variavel!r} é uma variável temporal com as fatias {fatias}.\n" + f"Escolha uma ou mais com --tempos, por exemplo:\n" + f" --tempos {fatias[0]}" + + (f" {fatias[-1]}" if len(fatias) > 1 else "") + ) from erro + raise SystemExit(str(erro)) from erro + if tempos is None: + tempos = [None] # uma banda só, sem fatia + por_id = {} + for item in layout: + por_id[item["tile_id"]] = item + por_id[item["tile_id"].replace("BDC_SM_", "")] = item + + escolhidos = [] + for t in pedidos: + if t not in por_id: + disp = sorted({i["tile_id"].replace("BDC_SM_", "") for i in layout}) + raise SystemExit( + f"Tile {t!r} não existe em {args.grid!r}.\nDisponíveis: {disp}" + ) + escolhidos.append(por_id[t]) + + # ── envoltória dos tiles escolhidos ───────────────────────────────────── + r0 = min(t["row_off"] for t in escolhidos) + c0 = min(t["col_off"] for t in escolhidos) + r1 = max(t["row_off"] + t["height"] for t in escolhidos) + c1 = max(t["col_off"] + t["width"] for t in escolhidos) + altura, largura = r1 - r0, c1 - c0 + + print(f"\ngrade {args.grid}: {grid.rows}x{grid.cols}") + print(f"recorte: linhas {r0}–{r1}, colunas {c0}–{c1} ({altura}x{largura})") + for i, t in enumerate(escolhidos, start=1): + print(f" {i}. {t['tile_id']:<16} em ({t['row_off']},{t['col_off']}) " + f"-> local ({t['row_off']-r0},{t['col_off']-c0})") + + # transform deslocado para a origem do recorte + base = grid.transform + transform = Affine(base.a, base.b, base.c + c0 * base.a, + base.d, base.e, base.f + r0 * base.e) + + saida = Path(args.saida) if args.saida else Path( + f"./{args.variavel}_{'_'.join(pedidos)}.tif" + ) + + perfil = dict( + driver="GTiff", height=altura, width=largura, + count=len(tempos) + 1, dtype="uint8", nodata=NODATA, + crs=grid.crs, transform=transform, + tiled=True, blockxsize=BLOCO, blockysize=BLOCO, + compress="deflate", zlevel=6, + ) + + # ── escrever, percorrendo os BLOCOS do destino ────────────────────────── + # Os tiles não são múltiplos do bloco do GeoTIFF, então escrever tile a + # tile deixaria a parte não coberta dos blocos de borda com zero em vez + # do nodata — e zero é valor válido aqui. + with rasterio.open(saida, "w", **perfil) as dst: + for banda, tempo in enumerate(tempos, start=1): + rotulo = args.variavel if tempo is None else f"{args.variavel}_{tempo}" + dst.set_band_description(banda, rotulo) + lay_ano = {i["tile_id"]: i for i in + cube.tile_layout(args.variavel, args.grid, time=tempo)} + raiz = Path(args.catalogo).resolve().parent + abertos = { + t["tile_id"]: zarr.open( + _resolver(lay_ano[t["tile_id"]]["url"], raiz), mode="r" + )[args.variavel] + for t in escolhidos + } + _escrever_banda(dst, banda, escolhidos, abertos, r0, c0, altura, largura, + valor=None) + print(f" banda {banda}: {rotulo}") + + # banda extra: índice do tile, para ver as fronteiras + idx = len(tempos) + 1 + dst.set_band_description(idx, "indice_do_tile") + _escrever_banda(dst, idx, escolhidos, None, r0, c0, altura, largura, + valor="indice") + print(f" banda {idx}: indice_do_tile (1..{len(escolhidos)})") + + with rasterio.open(saida) as ds: + print(f"\n{saida}") + print(f" {ds.height}x{ds.width}, {ds.count} bandas, crs={ds.crs.to_string()[:24]}") + for b in range(1, ds.count + 1): + a = ds.read(b) + v = a != NODATA + u = np.unique(a[v]) + print(f" {ds.descriptions[b-1]:<18} válidos={int(v.sum()):>10,} " + f"valores={u[:6]}") + return 0 + + +def _escrever_banda(dst, banda, escolhidos, abertos, r0, c0, altura, largura, valor): + for br in range(0, altura, BLOCO): + br1 = min(br + BLOCO, altura) + for bc in range(0, largura, BLOCO): + bc1 = min(bc + BLOCO, largura) + buf = None + for i, t in enumerate(escolhidos, start=1): + # posição do tile dentro do recorte + tr, tc = t["row_off"] - r0, t["col_off"] - c0 + sr0, sr1 = max(br, tr), min(br1, tr + t["height"]) + sc0, sc1 = max(bc, tc), min(bc1, tc + t["width"]) + if sr0 >= sr1 or sc0 >= sc1: + continue + if buf is None: + buf = np.full((br1 - br, bc1 - bc), NODATA, np.uint8) + if valor == "indice": + buf[sr0 - br:sr1 - br, sc0 - bc:sc1 - bc] = i + else: + trecho = np.asarray(abertos[t["tile_id"]][ + sr0 - tr:sr1 - tr, sc0 - tc:sc1 - tc]) + buf[sr0 - br:sr1 - br, sc0 - bc:sc1 - bc] = np.where( + np.isnan(trecho), NODATA, trecho).astype(np.uint8) + if buf is not None: + dst.write(buf, banda, window=Window(bc, br, bc1 - bc, br1 - br)) + + +if __name__ == "__main__": + raise SystemExit(main())