Correct US coastal depths to chart datum (MLLW) - #117
Open
bkeepers wants to merge 30 commits into
Open
Conversation
There was a problem hiding this comment.
🟡 Not ready to approve
There are a few concrete correctness/robustness issues in the new datum-surface build and provenance recording (stale bundle extraction, shared work dir, and non-canonical offset_surface sidecar values) that should be addressed before merging.
Once you've addressed the issues Copilot identified, you can request another Copilot review.
This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.
Pull request overview
This PR introduces a spatially-varying vertical datum correction pipeline to reference US CUDEM coastal depths from NAVD88 to chart datum (MLLW) using NOAA VDatum separation grids, and updates source metadata to reflect verified vertical datums.
Changes:
- Add
offset_surfacesupport to the source prep/datum pipeline so per-pixel datum separations can be applied during staging. - Introduce
pipelines/datum_grid.pyand a Snakemakedatum_surfacerule to compose/cache a NAVD88→MLLW reference surface from the pinned NOAA VDatum bundle with benchmark gating. - Update source metadata/docs to reflect verified datums (EMODnet=LAT, GSC Pacific=LLWLT) and change CUDEM sources from raw-mirrored to processed/prepared.
File summaries
| File | Description |
|---|---|
| sources/README.md | Updates source materialization/datum documentation for CUDEM processing and verified datums. |
| sources/gsc_pacific/metadata.json | Records verified vertical datum (LLWLT). |
| sources/gsc_atlantic/metadata.json | Records that vertical datum is undocumented/mixed. |
| sources/emodnet/metadata.json | Records verified vertical datum (LAT). |
| sources/cudem/metadata.json | Switches CUDEM to prepared source with offset_surface and datum notes. |
| sources/cudem/file_list.txt | Updates CUDEM ingest notes to reflect per-pixel datum correction. |
| sources/cudem_third/metadata.json | Switches CUDEM 1/3 to prepared source with offset_surface and datum notes. |
| sources/cudem_third/file_list.txt | Updates CUDEM 1/3 ingest notes to reflect per-pixel datum correction and AK limits. |
| Snakefile | Adds DATUM_SURFACES discovery and datum_surface build rule; wires surface into prep inputs. |
| pipelines/common.smk | Adds offset_surface() input resolver for prep rules. |
| pipelines/source_prep.py | Plumbs offset_surface into datum sidecar + prep-time transform, with checks and a self-test. |
| pipelines/source_datum.py | Implements striped per-pixel reference subtraction (--offset-surface) and expands self-checks. |
| pipelines/source_catalog.py | Carries offset_surface provenance into catalog properties (seascape:datum_surface). |
| pipelines/publish.smk | Prevents sweeping old raw mirrors when a source transitions raw→processed by excluding objects/**. |
| pipelines/datum_grid.py | New: composes NAVD88→MLLW reference COG from VDatum bundle; benchmark-gated. |
| pipelines/aggregation_reproject.py | Updates comments to reflect datum correction is now baked at prep. |
| Justfile | Adds source_datum.py --check and datum_grid.py --check to the sources test target. |
| docs/plans/2026-08-03-vdatum-inventory.md | New supporting inventory/ground-truth doc for VDatum bundle and sign/formula verification. |
| docs/plans/2026-08-03-datum-vdatum.md | New planning doc describing the approach and design decisions for VDatum-based correction. |
| CONTRIBUTING.md | Updates workflow docs to reflect CUDEM is no longer a raw source. |
Review details
- Files reviewed: 20/20 changed files
- Comments generated: 3
- Review effort level: Lite
We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.
There was a problem hiding this comment.
🟡 Human review recommended
It changes core geospatial datum semantics and overview/decimation behavior across multiple pipeline stages, so it warrants final human review despite only minor nits found in this pass.
Review details
Suppressed comments (2)
pipelines/source_prep.py:12
- The docstring example surface name "navd88_chart" doesn’t match the surface this PR introduces ("navd88_mllw"), which can mislead future source authors when setting offset_surface.
offset_surface reference raster subtracted per pixel for a spatially-varying datum
separation (source_datum --offset-surface); names a raster in the datum
store, e.g. "navd88_chart" (built by datum_grid.py)
pipelines/source_datum.py:23
- This docstring says source_normalize produces a "LERC" COG, but source_normalize writes ZSTD-compressed COGs; the wording here is now inaccurate.
tiled GeoTIFF; source_normalize makes the final LERC COG.
- Files reviewed: 30/30 changed files
- Comments generated: 0 new
- Review effort level: Lite
We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.
There was a problem hiding this comment.
🟡 Changes recommended
pipelines/terrain.py sets READ_RESAMPLE = "nearest" but gdalwarp -r expects near, which will fail terrain window reads at runtime.
Once you've addressed the issues Copilot identified, you can request another Copilot review.
This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.
Review details
- Files reviewed: 36/36 changed files
- Comments generated: 1
- Review effort level: Lite
We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.
| # _read_window), so this kernel only ever replicates a coarser neighbour tile — it must not reduce. | ||
| # A reducing kernel here would decimate a second time, class-blind, over whatever the pyramid | ||
| # already decided. Carried into _config so changing it reruns every terrain stem. | ||
| READ_RESAMPLE = "nearest" |
CUDEM ships on NAVD88, which sits between -0.53 m (Columbia estuary) and +1.7 m and more (Boston, Alaska) off MLLW, so its depths read too deep against chart datum and drying flats classify as navigable water. No scalar offset fixes a separation that varies by metres along the coast. - datum_grid.py composes a NAVD88-to-MLLW reference surface from the pinned VDatum bundle (vdatum_all_20250917). The composition branches on each region's .met horz frame (NAD83 / IGS14+geoid crossing / PRVI-as-LMSL); --check reproduces published CO-OPS separations at 14 benchmark stations spanning all branches and both signs, worst residual 3.4 cm. Full inventory in docs/plans/2026-08-03-vdatum-inventory.md. - source_datum gains --offset-surface: subtract a reference raster (chart-datum height in the source frame) resampled onto each file, striped so the largest CUDEM tile peaks under 400 MB; reference nodata leaves pixels uncorrected rather than dropping them. - cudem/cudem_third convert from raw (mirrored verbatim) to prepared sources: the correction is preparation, so prep applies it while staging, presence-on-disk is the incremental state, and the reference builds and caches on the box (rule datum_surface). mixed_crs declared (HI/Guam tiles carry different frames). Registration records seascape:datum_surface as provenance. - EMODnet/GSC datums verified against product docs (refs #26): EMODnet is LAT (GeoTIFF product), GSC Pacific is LLWLT, GSC Atlantic declares no vertical datum and stays fallback-tier. Plan and rationale: docs/plans/2026-08-03-datum-vdatum.md. Refs #16.
build() now runs the 14-station CO-OPS check against the freshly composed COG before it replaces the published surface, so a formula, region-selection, or input-layout regression fails the producing rule instead of silently mis-correcting every CUDEM depth. The valid-pixel count moves into fill_holes' existing block pass, dropping the post-publish rescan of the full multi-billion-pixel COG.
- harvest() clears the extracted vdatum/ tree before extraction, so a bundle-edition change can't mix regions from two editions into the composed surface. - build() keys its scratch dir to the output path, so concurrent surface builds don't race on one shared work tree. - The datum sidecar records the surface's canonical name (stem), never the invocation path, so seascape:datum_surface compares equal however source_datum was invoked.
Hawaii, Guam, CNMI, and American Samoa tiles carry genuinely different horizontal frames (WGS 84, NAD83(MA11)) from the NAD83 mainline, and 159 of the 163 are compound-CRS. Splitting them out (cudem_pacific, cudem_pacific_third; mixed_crs, no VDatum surface) leaves cudem and cudem_third CRS-homogeneous with Puerto Rico/USVI - verified by reading all 1,545 tile headers - so the mainline pair drops mixed_crs and takes the single-VRT warp path, measured 2x cheaper per aggregation tile. - Compound CRSs are flattened unconditionally at prep (copy-edit- replace, never in place: bare rasters stage as hard-links to raw/, and an in-place header edit corrupts the verbatim download - pinned by a byte-compare in source_prep's self-check). A materialized warp applies the geoid separation to any compound source it meets, so no staged raster may carry one. - create_virtual_raster now fails loudly when gdalbuildvrt drops off-CRS inputs, the silent mechanism that cost production the 2024-vintage CUDEM tiles.
A coarse pixel must never read deeper than the finest data under it, and AVERAGE pyramids violated that wherever depth decimates: a -38 m pinnacle charted at -187 m by z6, a -4.8 m Chesapeake shoal at -24 m by z10. Overview content IS the charted depth below native zoom (a planted wrong overview reaches the warp output), so this was live. - utils.shoal_overviews(): per-level gdalwarp -r max cascade (max is associative; GDAL ships no max kernel for gdaladdo/COG overviews) attached via gdal raster overview add. Internal overviews only: a .ovr sidecar is silently dropped by the COG translate, which then rebuilds an AVERAGE pyramid. - mosaic tile COGs + planet-z8, terrain's read resample, and depare's uniform coarsen move to the max kernel; prepared-source COGs gain a shoal pyramid (they had none), halving the aggregation warp read on a real CUDEM tile and lifting the shoals it stopped drowning +5.9 m. - mosaic_tile version=2: every mosaic tile rebuilds so no AVERAGE pyramid survives in served tiles. Land dilation from max is 0.1-2.4% of water pixels across the tested tiles; the 32768-px tile translate goes 32 s -> 145 s (half of that is the warp cascade, reducible to a striped numpy block-max if it hurts).
The pyramid previously cost three full-res passes (base GTiff, overview attach, COG re-translate) - a 4.5x slowdown on a 32768px tile. Now the levels come from a striped numpy 2x2 nodata-aware block max cascading off the level above, attached zero-copy as VRT <Overview> elements, and one COG translate with FORCE_USE_EXISTING reads the band once: +20% on the 32768px tile, and level values are exact block maxima rather than warp-kernel output. The decode pool is bounded at 4 threads - the box deliberately oversubscribes cores across jobs, so ALL_CPUS would thrash.
NOAA charts the Columbia against CRD (Columbia River Datum), not MLLW, and NAVD88-to-MLLW there runs to -0.53 m - correcting the river to MLLW charted it DEEPER than leaving it alone. The reference now repaints the Columbia reach from the bundle's crd.gtx, whose values are CRD height above NAVD88 verbatim (horz=NAD83, no geoid chain; verified against three CO-OPS CRD_OFFSET anchors to 2-3 cm, with Astoria confirming the mouth transition where CRD merges into MLLW). The surface renames to navd88_chart since it is no longer MLLW everywhere. Also fixes a composer defect found en route: src_nodata passed with a rasterio.band source is silently ignored (rasterio 1.4.4/GDAL 3.10.3), so the GTX sentinel -88.8888 bilinearly bled into the composed grid near region edges; sampling now goes through a masked array. All 16 benchmark stations (13 MLLW + 3 CRD) pass within 5 cm.
The benchmark stations validate formula and sign but cannot see a missing coast: out of coverage the reference is nodata and pixels pass through silently. The build now writes a per-region pixel-count sidecar beside the surface and refuses to publish when any of the bundle's 52 regions contributed nothing (naming it) or the total drifts more than 2% from the published grid's measured count. At prep, each file reports corrected/valid; warnings fire only in the 0.995-1.0 band (a sub-half-percent remainder is a hole, not a boundary tile - VDatum's edge costs boundary tiles percent-scale, and Alaska / the Pacific territories are documented zero-coverage no-ops). The source aggregate publishes as seascape:datum_corrected so partial coverage is catalog provenance, not a silent mix.
The shoal max was class-blind: one positive land sample in a 2x2 block turned a coastal water pixel positive at coarse zooms, and the terrain classifier then minted false drying or land in navigable water - compounding to 8.2% of water-bearing blocks charted as land at x64 on a real Chesapeake tile. The reduction now splits on DRYING_CAP: a block with any water-domain child takes the shoal max of those children alone; only all-land blocks stay land. The rule is associative, so levels still cascade. Zero water-bearing blocks read land at any level after; shoal preservation and nodata semantics are byte-identical. Two latent full-res-read defects fixed en route, both measured by planting wrong overview values: gdalbuildvrt -resolution highest over mixed resolutions ignores overviews entirely (terrain's coastal windows were full-res reads decimated by a class-blind max kernel), so window VRTs now pin -tr to the target resolution and read each tile's exact pyramid level 1:1; and the planet-z8 build warped the whole native-res mosaic through -r max (128x128 mixed-class blocks), now assembled from the tiles' own 512px z8 levels. depare's coarsen uses the same block reduce; -r max survives nowhere that mixes classes. Cost: +25-34% on the pyramid pass (~+6% of a tile write), with an all-water/all-land stripe fast path keeping open ocean at the old cost.
The served raster, the contours, and the depth-area polygons reduce the same mosaic through different paths, and nothing asserted they agree. test_consistency.py pushes one fixture (channel through land, lobed shoal, drying flat with sub-legible pond, noisy deep basin) through the real render/smooth/contour/depare stages at native and coarse zoom, then asserts seven bounded properties: raster==fork window up to the shoal-ward quantization step, raster class vs depare coverage within one pixel-centre disagreement, band bracketing per the ladder, the channel open in both layers at both zooms, and isobaths separating the raster to measured margins. Eight negative controls prove each assertion bites. One real defect is bounded rather than hidden, named by ZOOM_SMOOTHING_SLACK_M = 3.0: per-zoom sigma smoothing plus water widening lets narrow water chart up to 2.88 m deeper one zoom out than the shoalest native pixel under it (and 2 m deeper than its own depare band). The constant retires when smoothing moves to a physical-scale sigma or a smooth-once-decimate shape; any worsening fails the test.
Prep ran one file at a time: the live cudem job sat at load average 1.00 for 95 minutes on an 8-vCPU box while ~966 tiles transformed serially. The per-file unit (datum transform, CRS flatten, normalize + pyramid) is independent and already stripe-bounded in memory, so it now runs through a thread pool. Threads, not processes, measured on 16 real CUDEM tiles: 70 s vs 84 s at 4 workers and 0.7 GB less peak RSS, because rasterio, numpy, and the COG translate subprocess all drop the GIL and one shared block cache beats one per worker. End to end 187.9 s -> 72.8 s (2.58x) with a byte-identical output digest; a source with no datum knobs gets 1.52x, since the win concentrates in the transform. Worker count and GDAL's own thread count are now one arithmetic: DEFAULT_WORKERS = cpu_count//2 against GDAL_WORKER_THREADS = 4 keeps the product at the box's declared 2x oversubscription, and the COG write stops asking for ALL_CPUS inside a fan-out. Peak RSS at 4 workers is 2.5 GB against the rule's 8 GB reservation, so concurrency on the box is unchanged. Aggregates assemble in input order from completion-order results, so the coverage report and corrected fraction are identical however the pool finishes; a failing file cancels the rest and raises naming itself.
Prepared COGs gained a class-aware shoal pyramid, but only the sources that happened to re-prep this cycle got one; the rest keep reading full resolution and decimating through a blending kernel. A version param is the declared way to invalidate a rule (the profile's rerun-triggers carry params), so the next sources run re-preps all 25 with no force override.
Hawaii, Guam and CNMI shipped uncorrected on island datums that are MSL realizations (GUVD04-MSL and NMVD03-MSL both measure +3 mm at their origin stations), so their depths read deep by the local MSL-MLLW separation. Each territory becomes its own source carrying a scalar datum_offset_m: Hawaii +0.351, Guam +0.418, CNMI +0.393, verified live against the CO-OPS datums API across a full station sweep of the region. The conservative pick across an island's disagreeing stations is the MAXIMUM, not the minimum as the plan and inventory both said: the scalar is added to raise the bed, so an undersized one leaves part of the original deep bias while an oversized one charts shallower than truth. Both docs are corrected. American Samoa stays uncorrected and documented: Pago Pago publishes no ASVD02 row and its datums sit on a modified 2011-2019 epoch after the 2009 earthquake re-determination, so the sign of the offset can't be bounded. Per-island offset surfaces were rejected on measurement - within-island spread already matches between-island spread, so island polygons would cut Hawaii's worst residual only 0.101 m to 0.079 m, inside the transform's own uncertainty. Re-sweeping all 163 Pacific tile headers found every NCEI dataset internally single-CRS, so mixed_crs is dropped and all ten CUDEM sources take the single-VRT aggregation path.
Bare rasters staged as <source>_<enumeration-index>.tif, so one upstream insertion renamed every later tile and re-keyed the whole registration - a full re-aggregation on every NCEI re-survey, not once. Names now derive from the URL's filename plus the same 16 hex that names raw/<hash>, so the staged-to-raw pair is readable and insertion leaves neighbours untouched. The hash suffix is not belt-and-braces: nz_coastal carries 23 duplicate basenames because LINZ tiles a national grid, so stem-only naming would hard-fail that source. Collisions stay impossible by construction, and a duplicated item URL still errors naming the URL. Also corrects the shoal-pyramid follow-up. Measured: a re-prep with unchanged filenames, bounds, dims, datum and recipe writes a byte-identical catalog.json, and no downstream key carries pixel content - so corrected COGs sit inert until a registration change or a mosaic_tile bump invalidates the merge. And within one invocation there is no content curing, so the prep_source version bump must be absorbed by a sources run before any build, or the build re-preps inside itself and re-merges the planet unconditionally.
warp_mixed wrote a full materialized GTiff and re-read it to make the COG - two passes over a transient that reaches ~1.8 GB at z14. Each input now warps to its own lazy -of vrt -tap warped VRT, and gdalbuildvrt -te -tr recombines them (legal: they are all 3857 by then), so both branches of reproject() converge on one translate. The target grid is already an integer multiple of the target resolution from the mercator origin, so -tap lands every VRT on it with no geometry code, and assert_vrt_complete still gates silent drops. Measured against the single-CRS path: 1.47x CPU -> 0.60x, transient 480 MB -> 110 KB, peak RSS 861 -> 641 MB. The win is not overview selection (a planted-sentinel probe refutes that - CUDEM's mercator y resolution makes the x2 overview coarser than the target, so AUTO rightly declines it); it is deleting the materialize-then-recompress round trip. A warped VRT is also the only mosaic shape that cannot shift values, because a VRT cannot express a pixel-value transform. Both a materialized warp and a GTI apply the -21.9 m NAVD88-to-ellipsoid shift on a compound-CRS input; this path cannot.
CUDEM was bisected by datum AND zoom into ten sources, but the zoom split was a pipeline artifact: aggregation_covering derives each FILE's build depth from its own bounds and pixel dims, and the merge groups by (priority, maxzoom, source), so a source's finer files already form the earlier group and win. Verified on real headers - across 133 Hawaii tiles the 1/9" band resolves to z15 and the 1/3" band to z13, two zooms apart inside one source - and pinned by new self-checks in aggregation_covering and test_engine.check_priority. Five sources remain, split only by datum treatment: cudem on the VDatum surface, and one per Pacific island group on its own scalar. The cap stays at 13 because it is a build budget, not an overzoom guard - uncapping the 1/9 band to its native z15 is the native-resolution plan's own step. Consequence: the 22 macrotiles covered by 1/3 but not 1/9 move from child_z 12 to 13, +13% merge pixels over CUDEM's 141 macrotiles, at the 1/3 band's own native depth so nothing is upsampled. Guam keeps mixed_crs (its bands are 6325 and 4269 and now meet in one VRT); that path costs 0.60x the single-VRT path since 20be674.
The covering registers every file within 2*macrotile_buffer_3857 of a macrotile, which is wider than the reproject's buffer_pixels halo, so an edge file can land wholly outside the warp window. It contributes no pixel, and gdalbuildvrt -te drops it silently at exit 0 - which assert_vrt_complete then reads as the heterogeneous-CRS drop it exists to catch, failing a tile that is perfectly fine. Filter those inputs by their own warped extent before the combine, and fail loudly only if nothing reaches the window at all.
bkeepers
force-pushed
the
chart-datum-vdatum
branch
from
August 6, 2026 02:20
ea0d3a0 to
7b5382d
Compare
The lazy warped-VRT path lost its small-input win at S-102 scale, and the mechanism is /vsicurl dataset-open latency (~1.15 s wall against 0.09 s CPU per open), paid once per input per pass across three serial passes. With the network removed the lazy path beats the materialized one in both wall and CPU, so nothing about its pixel path is slow; the per-block-lookup, spatial-index and cache-pressure hypotheses were all probed and refuted, and GDAL's curl/pool knobs measured as noise. Three levers, each measured on the real 131-input z14 S-102 tile: gdalwarp -srcband warps the object directly, deleting the band-select VRT pass (two opens per input becomes one); the per-input warps run through a thread pool (latency-bound at ~8% CPU duty, so it does not multiply against GDAL_WORKER_THREADS); and the combine builds 16-wide chunk VRTs through the same pool, with assert_vrt_complete at chunk and root level and contiguous slices preserving last-input-wins. The combine also moves onto _run so a transient /vsicurl failure retries - 131 serial remote opens in one gdalbuildvrt is fragile, and a warp that died without retry is indistinguishable from a heterogeneous-CRS drop. Against the materialized path at 131 inputs: CPU 0.87x, peak RSS 0.49x, transient 610 MB to 0.7 MB, wall 1.09x - the residual is the single lazy translate reading destination-block-ordered, and the only way back to source-ordered reads is materializing, which reintroduces the transient and the NAVD88-to-ellipsoid shift a VRT structurally cannot apply. Output is bit-identical to the prior lazy shape.
CUDEM's gridder plants smooth ~70 m cones peaking just above datum in dredged channels - 7 in Thimble Shoal Channel, uncharted, rendering as drying dots in a federal fairway. They are counted, never filtered: nothing gets deeper, and the count feeds the confidence grade. Detector: an isolated at-or-above-datum patch under 25 m across whose entire 50 m ring is deeper than 8 m. The ring separates artifact from signal - a charted rock belongs to a reef or platform, so something within 50 m of it is shallow, while a dome stands alone; ring-median variants flag real Monterey rocks and were rejected on measurement, and -8 m sits in the measured gap between real drying rocks (3.4-7.0 m of surrounding water) and the artifacts (10.8-14.5 m). Runs inside transform_file's existing stripe loop at 0.4 s per 81122 tile (inside run variance): each stripe is scored one iteration late with its neighbours as context and owns patches whose first row it holds, so a dome straddling a stripe boundary is neither split nor double-counted. Ground truth: 7/7 on Thimble Shoal, 0 false positives across 40 real tiles (2.6 Gpx) including 8 New York Harbour tiles full of dredged channels. Publishes as seascape:dome_candidates, null (not 0) where no transform ran - "not measured" must not read as "clean".
The feather's factor came from zoom-level arithmetic (1 << levels), but resolved_maxzoom rounds toward overzoom, so a group levels down holds cells anywhere in (1<<levels, 2<<levels] - Jacksonville's 16 m S-102 file measured 7.8 render px and got factor 4, leaving 73% of its boundary in staircase runs, and at cz14 the same cells (3.9 px) fell under the integer gate and got no feather at all. group_cell_px now derives the group's coarsest cell straight from the registration rows (3857 bounds / pixel dims - already in memory, no I/O; 7.832 px computed vs 7.788 px measured by lattice phase, 0.6% apart). The gate becomes FEATHER_MIN_CELL_PX = 3.0: an integer levels gate can only sit on a bracket edge with a factor-2 blind spot, the original rationale (bounded displacement; no coastline redraw by a very coarse source) is enforced by sigma and the ceiling rather than the gate, and the same physical shoreline must be drawn the same way at neighbouring zooms - a gate at 4 splits cz15 (feathered) from cz14 (not). On the real 16 m file at the z15 grid: staircase runs 98.5% raw -> 73.1% at the old factor 4 -> 39.2% at the ratio, with valid area growing shoal-ward (+0.36 pp) and below the over-smoothing knee at 16. Topology preservation is untouched - only the scale input changed. build.smk's MERGE_CFG gains feather_max_factor, which was missing from the rerun key, so tuning the ceiling now invalidates merged tiles.
Zooming out could chart narrow water deeper than the shoalest native pixel under it - 2.88 m on the consistency fixture, 48 m in real data at a source-coverage boundary - because every zoom smooths its own window independently and smooth_array's shoal clamp pins raw values over a region whose physical extent doubles per zoom out. A physical-scale sigma was the documented fix and is refuted by measurement (up to 130 m worse): the class-aware reduction hands any mixed block to its water, so narrow water doubles in physical width per zoom-out and the reduction is non-extensive in the smoothing's input space - no kernel scale can reconcile the zooms. Instead each rendered tile is clamped shoal-ward against the shoal reduction of the already-served zoom below (water only, before classification). Each zoom keeps its own generalization; only monotonicity is enforced, and it holds by construction: serve(cz) = max(own, R(serve(cz+1))) with a shoal-ward quantizer, recursively up the cascade. Measured: 0 px deepening on the fixture and over the engine's full 69-tile cascade; depare's deep-side bracket closes to 0.00 m; native zoom byte-identical. ZOOM_SMOOTHING_SLACK_M retires; the consistency test asserts no deepening with no budget. The remaining COARSE_SHOAL_ALLOWANCE_M covers only the licensed shallow side. terrain._config gains finer_clamped so every stem re-renders once; terrain renders now chain by zoom, so the critical path grows by cascade depth while job count and per-job cost are unchanged.
- DRYING_CAP joins MERGE_CFG and prep_source's params: the class-aware
reduction branches on the cap, so it shapes source and mosaic pyramid
VALUES the aggregation warp reads back - a cap change previously
re-rendered terrain and re-cut depare against pyramids still built at
the old cap, the exact cross-layer split the consistency test exists
to catch and cannot see across processes.
- The planet z8 VRT gets assert_vrt_complete: gdalbuildvrt skips an
unopenable tile at exit 0, and a dropped tile there becomes a
shoreline-coloured nodata hole in every z0-z7 render with nothing
naming it.
- The land nudge widens from ==0 to <=0: the reduction hands a mixed
block to its water, so a land feature narrower than a block (sea
stack, cliff toe) reads negative at coarse zooms and charted as
depth; OSM land is authoritative here exactly as it already is for
the vector layers. The reduction's erosion of isolated above-water
features over open water is now asserted in the self-check as a
documented trade-off rather than a surprise.
- datum_grid stops inheriting pyramid builds from COG_OPTS (nothing
reads a reference pyramid), refuses an --out name it does not
compose (the datum_surface rule templates {name}, so a second
registered surface would silently receive this one's content), and
transform_file refuses a reference whose nodata is not float32-exact
- reference_on's src_nodata handling matches by value, and a
non-exact sentinel would interpolate into the correction as height.
* origin/main: Report no benchmark rather than a fabricated one Report build progress in one line per job # Conflicts: # Justfile
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
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.
US coastal depths (CUDEM) ship on NAVD88, a land survey datum that sits anywhere from half a metre below MLLW chart datum in the Columbia River estuary to 1.7 m and more above it in New England and Alaska. Uncorrected, Atlantic-coast depths read about a metre too deep and drying flats classify as navigable water, which is the unsafe direction for a chart. This PR references CUDEM to MLLW using NOAA's VDatum separation grids, applied per pixel.
How
pipelines/datum_grid.pycomposes a NAVD88-to-MLLW reference surface from the pinned VDatum bundle (vdatum_all_20250917.zip). The composition branches on each region's frame (NAD83; IGS14 with a geoid crossing; Puerto Rico/USVI, whose datums are local MSL realizations). The build refuses to publish a surface that fails to reproduce the published CO-OPS separations at 14 benchmark stations; the current worst residual is 3.4 cm. Ground truth for the formulas is in docs/plans/2026-08-03-vdatum-inventory.md.source_datumgains--offset-surface: it subtracts a reference raster (chart-datum height in the source's vertical frame) resampled onto each file, processed in stripes so the largest tile peaks under 400 MB.cudemandcudem_thirdchange from raw (mirrored verbatim) to prepared sources. The correction is preparation, so it runs while staging, and the reference surface builds and caches on the sources box. Registration recordsseascape:datum_surfacefor provenance.crd.gtxrepainted over the MLLW composition. Correcting the river to MLLW would have charted it up to 0.53 m deeper than leaving it alone. The reference surface is namednavd88_chartaccordingly, and the gate now covers 16 stations (13 MLLW plus 3 CRD anchors, with Astoria confirming the mouth transition), worst residual 3.4 cm.cudem_pacificandcudem_pacific_third: their horizontal frames genuinely differ from the NAD83 mainline, and 159 of 163 are compound-CRS. That leavescudem/cudem_thirdCRS-homogeneous with Puerto Rico/USVI (verified by reading all 1,545 tile headers), so the mainline pair takes the single-VRT warp path, measured 2x cheaper per aggregation tile.gdalbuildvrtdropping an off-CRS input is now a loud failure instead of a silent one, and compound CRSs are flattened unconditionally at prep.seascape:datum_corrected, so partial coverage is provenance rather than a silent mix.seascape:dome_candidates(7/7 on the known Thimble Shoal artifacts, 0 false positives over 40 tiles).Validation
Compound CRS
Benchmarking this PR surfaced a latent CRS problem in the CUDEM collection. The 2024-vintage tiles carry a compound CRS (NAD83 + NAVD88 height, EPSG:5498) while 2017-vintage tiles are plain NAD83, and both vintages occur inside single NOAA regions. Two consequences, measured on a Puget Sound aggregation tile with 22 CUDEM inputs:
gdalbuildvrtskips inputs whose CRS differs from the first file's, so 6 of the 22 tiles never entered the mosaic and their coverage fell through to coarser sources. Production has been missing that data since the 2024 tiles landed.This PR closes both paths off at prep: after the datum correction the staged COGs are rewritten with their 2D horizontal CRS (the compound vertical would be wrong metadata anyway, since the values are MLLW rather than NAVD88), and CUDEM aggregates through the per-tile-VRT
mixed_crspath, which keeps every tile regardless of its declared frame. The corrected benchmark tile carries all 22 inputs with no vertical shift. Measurements: the mixed-CRS bench notes in the session scratchpad; the drop and the shift both reproduce with stock GDAL against the raw NOAA files.Rollout
The first sources run builds the reference (about an hour, cached on the volume) and preps the four CUDEM sources (~220 GB of raws already cached). The next build is a planet-wide mosaic rebuild: the tile rule's version bump retires every AVERAGE pyramid in served tiles, and the US coast picks up the corrected CUDEM in the same pass. The retired mirror data is gone from the volume; the R2 copy (
source/cudem*/objects/, ~216 GB) is excluded from sync and awaits a deliberaterclone purge. The other 21 prepared sources carry no pyramids until deliberately re-prepped, which is a lost speedup rather than a regression. Open follow-ups are recorded in the plan doc's Follow-ups section.Plan: docs/plans/2026-08-03-datum-vdatum.md. Refs #16.
Updated through c52831e: per-island corrections, source consolidation, prep fan-out, dome signal, lazy mixed-CRS warp, ratio-sized feather, the zoom monotonicity clamp, and the independent review's fixes.