Fixing - #22
Merged
Merged
Conversation
Adds a domain-agnostic per-pixel value-remapping operator, filling the gap
between the zonal operators (which aggregate several source pixels into one
target cell) and the need to simply remap class codes 1:1.
The mechanism (apply a table) lives in disscube; the table itself is data
supplied by the caller via Variable.mapping, so land-cover, soil-class and
land-tenure regroupings all reuse this operator without disscube knowing
anything about any of those domains.
- Variable.mapping / Derivation.mapping carry {source_value: target_value}.
Being a plain field, it is folded into spec_hash() automatically, so two
derivations with different tables are always distinct products.
- Operator.requires_mapping mirrors the existing requires_class_code
fail-fast contract; Derivation validates it at construction time.
- Lookup is vectorized via np.searchsorted over sorted keys rather than a
dense 0..max array, so a sparse table with large codes costs the same as
a dense one.
- Unmapped values, source nodata and non-finite input all become NaN — the
same "no valid result here" convention the zonal operators already use.
- No purity coordinates are attached: with one source pixel per target cell
there is no sub-cell composition to summarise, so coverage/dominance
would be a constant carrying no information at the cost of two extra
full-size arrays.
Verified against real data: reproducing the BR-MANGUE papel_dominio_2024
classification through this operator (from band 3 of br_mangue_base_v2.vrt)
matches the existing GDAL pipeline's output pixel-for-pixel over a dense
2048x2048 window containing all six classes, including the exact nodata
count.
Tests: 17 new, full suite 95 passing.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
Adds a standalone case study showing the full path from raw tiles to a categorical domain raster, exercising the new generic `reclassify` operator together with disscube's existing tile-based derivation. The example is also the clearest demonstration of where the boundary between the cube and a project sits: nothing mangrove-specific lives in disscube — `reclassify` just applies a table that arrives as data — while the papeis_estados table and the two-band papel/elevation coupling stay in the example, since Variable derives one band per variable. Why tile by tile: a single derive() over the full 18352x21350 extent (391M cells) allocates tens of GB. derive(tile_id=...) keeps memory bounded; the whole extent completes in ~80s on a modest machine. Why the example writes the GeoTIFF itself: tools/zarr_to_tif.py converts one whole Zarr at a time, but reassembling N tiles into a single multi-band GeoTIFF needs windowed writes. Verified against the reference implementation it replaces (br_mangue's GDAL pipeline): identical output over all 391M cells in both bands, same CRS, transform and band descriptions. geomosaic is imported with a clear install hint — it is not a disscube dependency, and examples are not part of the installed package. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
The documented BDC workflow did not work. bdc_importer registers tiles as
`BDC_{LEVEL}_{tile}` (e.g. BDC_SM_009002), but derive() looked up only
`{grid_id}_{tile_id}` (e.g. BR/5km_009002), so the exact example in
docs/guides/bdc.md failed with:
ValueError: Tile 009002 with valid bbox not found for grid BR/5km
BDC tiles are grid-independent — the same envelope serves every grid
sharing the BDC CRS — so they are registered once, not per grid. Tile
lookup now tries, in order: the grid-scoped id, the tile_id as a fully
qualified id, then each BDC level.
Ambiguity is an error, not a precedence rule: tile ids are NOT unique
across levels. In the real V2 grids 189 ids exist in both SM and MD, and
they cover different areas (verified: BDC_SM_005004 is at 3152000,11425600
while BDC_MD_005004 is at 3680000,10897600). Resolving a bare ambiguous id
by trying SM first would silently derive the wrong extent, so it raises
and names the candidates, pointing at the fully-qualified form instead.
Docs updated accordingly: the loop in guides/bdc.md now passes the full
tile source id rather than splitting off the suffix, with a warning about
cross-level collisions, and the lookup order is documented in both
guides/bdc.md and architecture/tiling.md.
Tests: 12 new covering both conventions, precedence, ambiguity and the
not-found message. Full suite 107 passing.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
`_crop_to_grid` caught every exception from `clip_box` and fell back to returning the unclipped band. One of those exceptions is `NoDataInBounds`, raised precisely when the source and the target grid do not overlap — so "there is nothing here" was handled as "read everything", then reprojected into an all-nodata window. Measured on real data (BDC_SM tile 027005 against the BR-MANGUE ANADEM mosaic, EPSG:5880 -> BDC Albers): the fallback materialised the full 18352x21350 source, 31.6x the 3520x3520 target window, peaking at 2.42GB RSS to produce 50MB of nodata. With the fix the same tile stays at 0.18GB. This is a normal case, not an error: tile meshes are selected by envelope, so some tiles fall outside the data. `_crop_to_grid` now returns None for it and `_align_raster` builds the empty result directly, carrying the source nodata in `_disscube_nodata` so categorical operators mark every cell invalid exactly as they would for a reprojected empty window. Other exceptions keep the previous fallback: an unexpected failure should degrade to a correct if expensive read, never to wrong output. Tests: 7 new covering both the crop contract and the aligned result, including a guard that the shortcut does not swallow the normal path. Full suite 114 passing. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
Companion to 01_dominio_estrutural.py, producing the same product over the national BDC_SM mesh with the master grid in BDC Albers instead of an ad hoc pixel mesh in the source's own CRS. The pair makes the trade-off concrete rather than theoretical. Measured on the same data: the ad hoc mesh runs 30 tiles in 80s under 1GB and reproduces the source bit for bit, while the BDC mesh runs 32 tiles in 100s peaking at 3.7GB and loses ~1% of cells per class to edge resampling — the cost of being on a canonical, shareable grid. The README now states when each is the right choice. Tile selection densifies the source footprint before projecting it: taking min/max of the four transformed corners overestimates the area, because straight edges become curves in another projection. That over-selection picked up 3 tiles that do not touch the data (35 vs the correct 32). They would now derive cheaply as empty, but selecting them at all is wasted work. Passes the full tile id (BDC_SM_027005) rather than the bare number, since ids are not unique across BDC levels. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
The mosaic's extent is a rectangle but the valid data is a coastal strip, so most tiles intersecting the extent cover nothing but nodata: 17 of 32 derived to fully empty rasters. Tile selection now takes one decimated pass over the source to build a low-resolution validity mask, and keeps only the tiles that touch it. Memory for that pass is bounded by the target width, not by the source size. Measured on the same data: 32 -> 15 tiles, 99.9s -> 56.9s, peak RSS 3.70GB -> 2.81GB, with class counts identical to the previous run (zero valid cells lost). Disk was never the problem — the empty Zarrs compressed to almost nothing, so the store only went 257MB -> 252MB. The win is compute, not storage. The mask is dilated by one cell: nearest-neighbour decimation can drop thin features, and the mangrove strip is a few pixels wide in places. The costs are asymmetric — a false positive derives one tile for nothing, a false negative silently loses data from the product. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
820a5fc fixed this for the BDC example but left 01_dominio_estrutural.py deriving the full rectangular mesh, where half the tiles cover nothing but nodata — the same waste, reported by the user after running it. Same approach: one decimated pass over the source builds a low-resolution validity mask, and only tiles touching it are derived. Simpler here than in 02_, since the tiles are pixel windows in the source itself, so the mask is indexed directly with no CRS transform. Measured: 30 -> 15 tiles, 80s -> 36.3s, and the output is still bit for bit identical to the reference GDAL pipeline over all 391M cells in both bands. Skipping is safe precisely because the mask comes from the state band: the reference implementation writes 255 to BOTH bands wherever state is nodata, so a tile with no valid state contributes nothing but the nodata the output already holds. README's comparison table updated with the new figures. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
Both examples created their own catalog.db and store under a per-example directory, unlike the other twelve examples, which all use CubeClient(catalog="catalog.db", store="./data/"). That was not merely cosmetic: with a private catalog the derived variables never showed up in cube.search() alongside the rest, and could not be combined through to_lucc_data() — losing the reason to catalog them at all. Each example was an island. Now both share the standard catalog and store, and follow the repo's directory conventions: VRTs go to data/raw/brmangue/ next to the other raw inputs, the GeoTIFFs to data/. Verified by running both: one catalog holds the two grids (brmangue/30m in EPSG:5880 and brmangue/30m_bdc in BDC Albers) with 60 DerivedVariable entries — papel and elevacao, 15 tiles each, on both grids. The two meshes are now comparable within the same cube. Note that both examples write the same VRTs and register the same SpatialSource ids, deriving them from the same ANADEM tiles. That is idempotent today, but the second run would silently overwrite the first if they ever diverge. Reported by the user. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
A variable derived tile by tile exists as N Zarr stores, and load() refuses that case because there is no automatic mosaic (docs/architecture/tiling.md). Until now the package stated the problem and offered no way out: the caller hit a ValueError listing the tiles and was on their own. tile_layout() answers the question load() cannot — where each piece sits on the master grid — without assembling anything and without imposing a destination. It returns plain data (path plus position), not an object from another package, so whoever consumes it decides: load into memory, write to a disk workspace, or just inspect coverage. Same arrangement geomosaic already uses when it hands back tile_offsets. Dimensions come from the registered bbox, not from the file: it is the same bbox derive(tile_id=...) cropped with, so it states the intended position. A store whose contents disagree is an inconsistency to catch, not to paper over. A variable with no tiles returns a single item covering the whole grid, so callers can treat both cases the same way. This makes the example 03 collapse from ~100 lines of offset arithmetic to six, and it now demonstrates use rather than working around a limitation. Its seam checks moved to haloexec as tests (f22057d) — they were proving the assembly code, not the data, and re-running them on every load protects nothing. Verified after the rewrite: the same 10 seams on the real BDC grid still check out against the source Zarrs, 0 divergences. Tests: 13, aimed at the position arithmetic, which is where an error hides — a wrong sign on row_off flips the mosaic vertically and nothing raises. Also covers both tile registration conventions, the global case, and the errors. Full suite 127 passing. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
A variable derived with valid_from/valid_until has one set of pieces PER SLICE, all at the same tile positions. tile_layout() returned them together, so every position appeared once per year — and assembling from that would let one year overwrite another with no error at all. Found by running the real thing: the MapBiomas mangrove series over five marker years came back as 140 pieces in 28 positions, each repeated five times with a different spec_hash in the URL. The existing tests missed it because they only covered static variables. tile_layout() now takes an optional `time`, and refuses a temporal variable when it is omitted, naming the available slices — the same discipline load() already applies to multi-tile. Each item also carries `times`, empty for static variables, so a caller can tell which slice it is holding. Verified on the series: without time it raises listing [1985, 1995, 2005, 2015, 2024]; with time=<year> each slice yields 28 pieces in 28 distinct positions, no overlap. Static variables are unaffected and still need no time. Unrelated but worth recording, since it was the run that found this: the derived tiles were checked one by one against the source VRT window — 140/140 identical, and the 2024 count (7,527,925 mangrove cells) matches exactly what the structural domain reports for class 1, by an independent route. Tests: 7 new for the temporal cases. Full suite 134 passing. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
Brings the MapBiomas binary mangrove presence into the cube as a TEMPORAL variable: one slice per year, so load() returns (time, y, x). The earlier examples only derive static variables, so nothing yet exercised the slice machinery — the SpatialSource time=<year>, the valid_from/valid_until on the derivation, or tile_layout()'s per-slice selection. Source shape drives the structure: the mangue_* products hold one DECADE per file with one year per band, and write_vrt mosaics a single band per call, so each year becomes its own VRT. The loop calls gc.collect() after every derive, which is not incidental. Each derive() leaves ~50 objects in reference cycles holding ~380 MB. Reference counting cannot break a cycle; only the generational collector can, and it almost never fires here because its trigger is the NUMBER of allocations while numpy concentrates hundreds of MB in very few objects — the collector never sees the memory pressure. Measured on a cold run of all 140 derivations: without the collect, RSS grows ~380 MB per tile, passes 7 GB and the run was killed at tile 21 of year 2015; with it, peak RSS is 1.07 GB and the run finishes in 100.8s. Seven times less memory for about 7% more time (0.05s per collect against ~0.7s per derive). The call is in the example rather than inside derive() because imposing that cost on every caller is the package's decision, not this example's. The reasoning is written at the bottom of the file so it does not read as a stray defensive line. README updated with 03 and 04, which were both missing from the listing. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
Two things, both found by running the pipeline on real data. The memory one is new and practical: each derive() leaves ~50 objects in reference cycles holding ~380 MB, and the generational collector almost never fires to reclaim them, because its trigger is the NUMBER of allocations while NumPy concentrates hundreds of MB in very few objects. A long loop therefore grows until it dies — measured at 140 tiles, RSS passed 7 GB and the run was killed; with gc.collect() per call it stayed near 1 GB, for about 7% more time. Documented under known limitations, with the workaround and a pointer to example 04, until the collection is done internally. The stale one: both the README and docs/architecture/tiling.md still said load() "silently returns the first result" for a multi-tile variable. It raises ValueError naming the tiles — the explicit error they described as planned is already there. Corrected, and tiling.md now documents tile_layout() as the way to actually consume a multi-tile variable, including the temporal case where a slice must be chosen. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
One band per year, written window by window from the Zarr tiles of each slice. It is not needed by anything in the cube — a multi-tile Zarr is not directly viewable, so this exists to open in QGIS and confirm the slices landed correctly. tools/zarr_to_tif.py does not cover this: it converts one whole store at a time, so it cannot reassemble N tiles, and it currently picks data_vars[0], which in a DisSCube Zarr is spatial_ref rather than the data. Verified against the same source the series was checked against: the five bands report 2.91% to 2.97% mangrove, matching what the Zarrs hold. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
…port Combines what 04 and 02 do separately: the temporal mangrove series, but on the national BDC_SM mesh with the master grid in BDC Albers, so every tile goes through a real 5880 -> Albers reprojection rather than the identity path. 120 derivations, 170s cold, peaking at 1.47 GB — the gc.collect() from 04 is what keeps that bounded. The export deserves a note, because a first version of it was wrong in a way that produced invented data. Writing tile by tile leaves the uncovered part of each border block filled with GDAL's default zero rather than the declared nodata, and here zero is a VALID value (non-mangrove). The result had 19.6 million cells outside the study area reading as observed non-mangrove, in two of the five bands. The cause is that BDC tiles are 3520x3520 while the GeoTIFF blocks are 512, and 3520 is not a multiple of 512, so every tile edge lands mid-block. Examples 01 and 02 never hit this because their tiles are 4096, an exact multiple. The fix is the same lesson the Zarr loader already encodes: walk the blocks of the DESTINATION, not the tiles, and build each block from whatever covers it. Verified after the fix: all five bands report 206,007,691 valid cells, matching the Zarrs exactly, with zero cells outside the tiled area. Also records what the reprojection costs — class counts land ~1% below the EPSG:5880 series of example 04, which does not resample. That is the price of a canonical mesh, not an error. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
Checking whether tiles landed in the right place meant opening the whole mosaic, and where the values on both sides of a boundary agree, that shows nothing at all. This crops only the requested tiles to their own envelope and adds a band carrying the tile index, which makes the boundaries visible even there. Styled as categorical in QGIS, the seams show up directly. It complements tools/zarr_to_tif.py rather than replacing it: that one converts a single store, this one assembles several tiles — which is the case a multi-tile variable actually presents, and the one zarr_to_tif cannot cover. Writes by walking the blocks of the GeoTIFF, not the tiles. When the tile side is not a multiple of the block (BDC tiles are 3520, blocks are 512), every tile edge lands mid-block and the uncovered part of those blocks keeps GDAL's default zero rather than the declared nodata. Where zero is a valid value that is invented data — the same trap example 05 fell into before being fixed. Resolves the catalog's relative asset_url against the catalog's own directory, so the tool works when run from anywhere, not only from the repository root. Errors speak the tool's language rather than the API's: an unknown tile is refused listing what exists, and a temporal variable without --tempos names the available slices and shows the flag to use. README documents it with both modes and explains the block-walking, so the reason survives longer than this session. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MHivBmXr3hkX1H7SayXXkB
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.