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/disscube/client/cube_client.py b/disscube/client/cube_client.py index 04fee9e..4205aa0 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 @@ -206,6 +270,148 @@ 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, + time: Optional[int] = None, + ) -> 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. + 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 + ------- + 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 + - ``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, + 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}" + ) + + # 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: + layout.append({ + "tile_id": None, + "variable": d.name, + "url": d.asset_url, + "row_off": 0, + "col_off": 0, + "height": grid.rows, + "width": grid.cols, + "times": list(d.times or []), + }) + 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)), + "times": list(d.times or []), + }) + + return sorted(layout, key=lambda t: (t["row_off"], t["col_off"])) + def to_lucc_data( self, variables: List[str], 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/disscube/pipeline/aligner.py b/disscube/pipeline/aligner.py index 570d3ee..0094b5d 100644 --- a/disscube/pipeline/aligner.py +++ b/disscube/pipeline/aligner.py @@ -24,12 +24,15 @@ 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 rioxarray.exceptions import NoDataInBounds from shapely.geometry import box from disscube.operators.base import OPERATOR_REGISTRY @@ -97,6 +100,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 @@ -139,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)) @@ -186,12 +220,157 @@ 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) # ------------------------------------------------------------------ - 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 @@ -201,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 @@ -231,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/docs/architecture/tiling.md b/docs/architecture/tiling.md index 857fae3..4ebb70c 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`. @@ -73,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 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/examples/README.md b/examples/README.md index 30ea798..ab67ee9 100644 --- a/examples/README.md +++ b/examples/README.md @@ -20,7 +20,51 @@ 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. +- `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)`. +- `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); +- **`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. + +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 | 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 +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. - `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. @@ -32,10 +76,39 @@ Dois estudos sobre a mesma área geográfica e grade. | 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/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..6df3037 --- /dev/null +++ b/examples/case_studies/brmangue_dominio/01_dominio_estrutural.py @@ -0,0 +1,259 @@ +""" +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" + +# 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 +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 _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( + 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() + 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") + 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(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=CATALOGO, store=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), + )) + + # 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, 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( + 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} com dado válido " + f"({vazios} descartados por serem só nodata)") + + 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() 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..77774d8 --- /dev/null +++ b/examples/case_studies/brmangue_dominio/02_dominio_bdc.py @@ -0,0 +1,337 @@ +""" +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" +# 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 +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 _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( + 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() + 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") + 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(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") + 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 + + 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)) + if not candidatos: + raise RuntimeError("Nenhum tile BDC_SM cobre a extensão da fonte.") + 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) + 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=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], + 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() 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..230402e --- /dev/null +++ b/examples/case_studies/brmangue_dominio/03_carregar_no_haloexec.py @@ -0,0 +1,94 @@ +""" +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. + +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. 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 +""" + +from pathlib import Path + +from disscube.client import CubeClient + +try: + 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 + + +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 main() -> None: + cube = CubeClient(catalog="catalog.db", store="./data/") + 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/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']}") + + print(f"\n[2/2] montando o workspace em {WS_DIR}") + ws = MemmapRasterWorkspace.create( + WS_DIR, shape=(grid.rows, grid.cols), + arrays={nome: DTYPE for nome in VARIAVEIS}, + block_h=BLOCK_H, block_w=BLOCK_W, halo=HALO, + ) + for nome, tiles in layouts.items(): + load_zarr_tiles_into_workspace(ws, tiles, array=nome) + print(f" {nome}: carregado") + + 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__": + main() 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..e4b48b7 --- /dev/null +++ b/examples/case_studies/brmangue_dominio/04_serie_temporal_mangue.py @@ -0,0 +1,286 @@ +""" +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 + +# 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] + +# 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 _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( + 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/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}") + + # 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") + + # ── 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)") + 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() 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() 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) 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 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} diff --git a/tests/test_tile_layout.py b/tests/test_tile_layout.py new file mode 100644 index 0000000..a158401 --- /dev/null +++ b/tests/test_tile_layout.py @@ -0,0 +1,250 @@ +""" +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", + "times", + } + + +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) + + +# ── 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"] == [] 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 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())