Allen Coral Atlas SDB (Bahamas/N. Caribbean); multi-source pipeline improvements

ROADMAP.md

@@ -276,10 +276,13 @@ EMODnet already covers European seas **including the N. African Med shelf** (the
 | EMODnet DTM 2024 _(already ingested)_ | 115 m | European seas + **N. African Med shelf** | **LAT** ✓ | CC-BY 4.0 ✓ | z11 | covers the N-African Med coast as part of the existing product; EMODnet is European/NE-Atlantic so nothing reaches the Caribbean |
 | swIOBC | 250 m | SW Indian Ocean / E. Africa | ~MSL | CC-BY 3.0 ✓ | z9 | **BUILT** (B) — one ~711 MB GeoTIFF |
 | GMRT | ~100 m swaths | global multibeam | mixed | CC-BY 4.0 ✓ | z9–11 | OPPORTUNISTIC — patchy; "not for navigation" |
+| **Allen Coral Atlas (Caribbean/Bahamas)** | 10 m SDB | Bahamas + N. Caribbean shallow reef/bank (0–~25 m) | ~water-surface | **CC-BY 4.0** ✓ | z13 | **TABLED** (ingested as `sources/aca_ncaribbean`, recipe works, raw mirrored to R2) — but SDB is too noisy for a *chart*: ±1–2 m speckle, and false-shallow returns over km-deep water. A GEBCO-depth-gated feather fixes the terrain (deep→GEBCO, shallow→ACA), but 1 m shallow contours over the vast banks still explode to ~1 M features/macrotile (tippecanoe 45 min+), and coarsening them is unacceptable — fine shallow contours are essential for navigation. Revisit with ATL24 (ICESat-2, RMSE ~0.3 m, open) as depth control / a C-SHELPh re-derivation. Details: PR + `[[aca-bahamas-download-pending]]`. |
+| TCarta SDB (Caribbean GeoPortal) | 10 m / 100 m | Bahamas/Lucayan, Belize, Jamaica | n/a | **proprietary EULA — "NOT FOR NAVIGATION", no redistribution ✗** | — | SKIP — public download but all-rights-reserved. 10 m is RGB tiles only (no depth); 100 m `.lpkx` holds a real depth tif but is unlicensed for reuse. TNC Caribbean layers there are benthic/habitat, not depth. |
+| Bahamas Median DEM (IEEE DataPort) | 30 m | Andros I. only (ICESat-2 + Landsat-8) | m depth, datum unstated | **IEEE-subscription-gated; no stated data license ✗** | — | SKIP — "open source" = the method, not the data; gated, Andros-only ([DOI 10.21227/dwex-hm52](https://dx.doi.org/10.21227/dwex-hm52)) |
 | Red Sea / Strait of Tiran patches | 10–30 m | Red Sea rift + Tiran | unstated | CC-BY 4.0 ✓ | z12 | OPPORTUNISTIC (B) — Tiran coastal; rest deep curiosities |
 | Chilean fjord grids | 10–50 m | S. Chile fjords | SHOA-CD | per-record — **verify** | z12 | OPPORTUNISTIC (B) — license-gated scattered patches |
 | Brazil LEPLAC / de Wet SA shelf / Lesser Antilles / EOMAP / Israel / Mexico IBCCA | varies | — | — | **study-only / NC / no-license / commercial ✗** | — | SKIP |
-| Brazil DHN, Chile SHOA, Argentina SHN, Peru/Colombia/Ecuador, Caribbean states, SANHO, W/E Africa, Arabian Gulf | — | populated coasts | — | **closed / request-only ✗** | — | **GEBCO-only** — no open hi-res source exists |
+| Brazil DHN, Chile SHOA, Argentina SHN, Peru/Colombia/Ecuador, Caribbean states, SANHO, W/E Africa, Arabian Gulf | — | populated coasts | — | **closed / request-only ✗** | — | **GEBCO-only** for the national-HO grids — no open hi-res source exists (Caribbean/Bahamas shallow reefs are the exception: Allen Coral Atlas, above) |
 
 ### Inland waters — lakes & rivers (separate layer; Great Lakes covered elsewhere)
 

pipelines/aggregation_run.py

@@ -9,6 +9,9 @@
   aggregation_run.py freeze          write the dirty list into the covering (plan, once)
   aggregation_run.py shard <i> <n>   process the frozen dirty[i::n] (matrix shard i of n)
   aggregation_run.py matrix <max>    print the shard matrix JSON (sized to the dirt)
+  aggregation_run.py step <name>...  run only the named stage(s) over the covering (reproject |
+                                     merge | smooth | tile | contour) — stages run independently,
+                                     no skip flags; omit a stage by leaving it out
 
 `shard` takes a strided slice of the single dirty list `freeze` wrote, so every shard
 partitions the identical list — no overlap, no coordination, nothing recomputed per shard.
@@ -28,27 +31,46 @@
 import smooth
 import utils
 
-# Both forks share the merged DEM; set these to 1 for raster-only / no-smooth runs.
-SKIP_CONTOURS = os.environ.get("SKIP_CONTOURS", "")
-SKIP_SMOOTH = os.environ.get("SKIP_SMOOTH", "")
+# The aggregation stages, each operating on one tile's covering CSV. There are NO skip flags:
+# run the full sequence (production / `just planet` / a shard) or any subset independently via
+# the `step` CLI. To omit a stage (e.g. the slope blur in a value-exact test), just don't run it.
+def _smooth(filepath):
+    smooth.smooth_merged(filepath.replace("-aggregation.csv", "-tmp"))  # slope-selective blur
+
+
+STEPS = {
+    "reproject": aggregation_reproject.reproject,
+    "merge": aggregation_merge.merge,
+    "smooth": _smooth,
+    "tile": aggregation_tile.main,        # raster Terrain-RGB tiles
+    "contour": contour_run.generate,      # vector contours off the merged DEM
+}
+FULL = ["reproject", "merge", "smooth", "tile", "contour"]
 
 
 def run(filepath):
+    """Full pipeline for one tile: every stage, then drop the tmp DEM and mark it done."""
     item = filepath.split("/")[-1].replace("-aggregation.csv", "")
     print(f"{item} start")
-    aggregation_reproject.reproject(filepath)
-    aggregation_merge.merge(filepath)
-    tmp_folder = filepath.replace("-aggregation.csv", "-tmp")
-    if not SKIP_SMOOTH:
-        smooth.smooth_merged(tmp_folder)     # slope-selective blur, shared by both
-    aggregation_tile.main(filepath)          # raster Terrain-RGB tiles
-    if not SKIP_CONTOURS:
-        contour_run.generate(filepath)       # vector contours off the merged DEM
-    shutil.rmtree(tmp_folder)
+    for name in FULL:
+        STEPS[name](filepath)
+    shutil.rmtree(filepath.replace("-aggregation.csv", "-tmp"))
     utils.run_command(f'touch {filepath.replace("-aggregation.csv", "-aggregation.done")}')
     print(f"{item} end")
 
 
+def run_steps(names):
+    """Run only the named stage(s) over every covering tile, in order, leaving the tmp DEM in
+    place (no cleanup, no .done) so a later invocation can pick it up — for running stages
+    independently (dev / the engine test), with no skip flag."""
+    bad = [n for n in names if n not in STEPS]
+    if bad:
+        sys.exit(f"unknown step(s) {bad}; choose from {list(STEPS)}")
+    for filepath in covering_sorted():
+        for name in names:
+            STEPS[name](filepath)
+
+
 def covering_sorted():
     """Every aggregation CSV in the current covering, heaviest-first. Depends ONLY on the
     immutable covering, never on which tiles are already built — so every shard derives the
@@ -144,10 +166,12 @@ def main(argv):
     elif argv[:1] == ["shard"]:
         i, n = int(argv[1]), int(argv[2])
         run_all(work_list()[i::n])
+    elif argv[:1] == ["step"]:
+        run_steps(argv[1:])
     elif not argv:
         run_all(dirty_filepaths())
     else:
-        sys.exit("usage: aggregation_run.py [freeze | shard <i> <n> | matrix <max>]")
+        sys.exit("usage: aggregation_run.py [freeze | shard <i> <n> | matrix <max> | step <name>...]")
 
 
 if __name__ == "__main__":

pipelines/smooth.py

@@ -16,7 +16,7 @@
 Sigma is in merged-DEM pixels, so the physical blur scale tracks the tile's zoom
 (coarse base tiles blur more in metres, fine regional tiles less) — roughly what we
 want (coarse data is noisier). Revisit with a physical-scale sigma if it
-over/under-blurs. SKIP_SMOOTH=1 disables it.
+over/under-blurs. Run as the `smooth` stage of aggregation_run (omit that stage to skip it).
 """
 
 import glob

pipelines/source_datum.py

@@ -1,7 +1,9 @@
-"""Apply the bathymetry value transform: ``negate`` then ``datum_offset_m``.
+"""Apply the bathymetry value transform: ``scale`` then ``negate`` then ``datum_offset_m``.
 
 Reads the knobs from ``metadata.json``:
 
+  - ``scale``: multiply valid pixels for a unit conversion, applied first — e.g.
+    ``0.01`` for a source stored in centimetres (Allen Coral Atlas SDB) to reach metres.
   - ``negate``: flip positive-down depth sources (e.g. DDM, stored as +depth) to
     negative-down elevation.
   - ``datum_offset_m``: constant added to bring the source to ~MSL (a single
@@ -22,14 +24,16 @@
 import rasterio
 
 
-def transform_file(filepath, negate, offset, clamp_positive=False):
+def transform_file(filepath, negate, offset, clamp_positive=False, scale=1.0):
     with rasterio.open(filepath) as src:
         profile = src.profile
         data = src.read(1)
         mask = src.read_masks(1) != 0  # True where valid
 
     data = data.astype("float32")
     valid = data[mask]
+    if scale != 1.0:
+        valid = valid * np.float32(scale)
     if negate:
         valid = -valid
     if offset:
@@ -58,19 +62,22 @@ def main():
     p.add_argument("source")
     p.add_argument("--negate", action="store_true", help="flip positive-down depth to negative-down elevation")
     p.add_argument("--offset", type=float, default=0.0, help="metres added to reach ~MSL")
+    p.add_argument("--scale", type=float, default=1.0,
+                   help="multiply valid pixels (unit conversion, applied before negate) — "
+                        "e.g. 0.01 for centimetre depths to metres")
     p.add_argument("--clamp-positive", action="store_true",
                    help="after the offset, drop cells > 0 (above the water surface) to nodata — "
                         "removes a lake DEM's land fringe / a topobathy playa")
     a = p.parse_args()
 
-    if not a.negate and a.offset == 0 and not a.clamp_positive:
-        print(f"{a.source}: no datum transform (negate=False, offset=0)")
+    if not a.negate and a.offset == 0 and a.scale == 1.0 and not a.clamp_positive:
+        print(f"{a.source}: no datum transform (negate=False, offset=0, scale=1)")
         return
     filepaths = sorted(glob(f"store/source/{a.source}/*.tif"))
-    print(f"{a.source}: negate={a.negate} offset={a.offset} clamp_positive={a.clamp_positive} "
-          f"on {len(filepaths)} file(s)")
+    print(f"{a.source}: scale={a.scale} negate={a.negate} offset={a.offset} "
+          f"clamp_positive={a.clamp_positive} on {len(filepaths)} file(s)")
     for filepath in filepaths:
-        transform_file(filepath, a.negate, a.offset, a.clamp_positive)
+        transform_file(filepath, a.negate, a.offset, a.clamp_positive, a.scale)
 
 
 def _check():
@@ -108,6 +115,19 @@ def _check():
         o2 = src.read(1)
     assert o2[0, 0] == -50.0 and o2[0, 1] == -10.0, o2  # bed kept
     assert o2[1, 0] == nodata and o2[1, 1] == nodata, o2  # +5 land clamped; nodata untouched
+
+    # scale + negate: centimetre depths (ACA SDB) -> metres, flipped to elevation
+    path3 = os.path.join(d, "t3.tif")
+    arr3 = np.array([[100.0, 313.0, nodata]], dtype="float32")  # +depth in cm
+    with rasterio.open(path3, "w", driver="GTiff", height=1, width=3, count=1,
+                       dtype="float32", nodata=nodata, crs="EPSG:4326",
+                       transform=from_origin(0, 1, 1, 1)) as dst:
+        dst.write(arr3, 1)
+    transform_file(path3, negate=True, offset=0.0, scale=0.01)  # cm depth -> m elevation
+    with rasterio.open(path3) as src:
+        o3 = src.read(1)
+    assert abs(o3[0, 0] - (-1.0)) < 1e-6 and abs(o3[0, 1] - (-3.13)) < 1e-6, o3
+    assert o3[0, 2] == nodata, o3  # nodata untouched by scale+negate
     print("source_datum.py self-check ok")
 
 

pipelines/source_polygonize.py

@@ -14,16 +14,37 @@
 from multiprocessing import Pool
 import shutil
 
+import rasterio
+
 import utils
 
 SILENT = True
 
+# Cap the polygonized mask's longest side (pixels). Pixel-exact polygons of a speckly mask — e.g.
+# scattered satellite-derived reefs across a 16k-px Allen Coral Atlas tile — blow past gdal/sqlite
+# vertex limits ("sqlite3_bind_blob() failed: too big"). Larger masks are downsampled to this first,
+# only ever downscaling and taking the max per block so coverage never shrinks (it generalizes
+# outward, the conservative direction). The footprint feeds the covering (which aggregation tiles to
+# consider); ~native/N precision is ample there. Files already <= this are untouched.
+COVERAGE_MAX_PX = 1024
+
 
 def polygonize_tif(source, filename):
+    src_tif = f"store/source/{source}/{filename}"
     mask = f"store/polygon/{source}/{filename}"
     utils.run_command(
-        f'GDAL_CACHEMAX=1024 gdal_calc.py -A store/source/{source}/{filename} '
+        f'GDAL_CACHEMAX=1024 gdal_calc.py -A {src_tif} '
         f'--outfile={mask} --calc="A*0+1" --type=Byte --overwrite', silent=SILENT)
+    with rasterio.open(src_tif) as s:
+        width, height = s.width, s.height
+    factor = max(1, -(-max(width, height) // COVERAGE_MAX_PX))  # ceil-divide -> downscale factor
+    if factor > 1:
+        coarse = mask + ".coarse.tif"
+        utils.run_command(
+            f'GDAL_CACHEMAX=1024 gdalwarp -r max '
+            f'-ts {max(1, width // factor)} {max(1, height // factor)} '
+            f'-overwrite {mask} {coarse}', silent=SILENT)
+        os.replace(coarse, mask)
     utils.run_command(
         f'GDAL_CACHEMAX=1024 gdal_polygonize.py {mask} -b 1 -f "GPKG" '
         f'store/polygon/{source}/{filename}.gpkg -overwrite', silent=SILENT)

pipelines/source_smooth.py

@@ -0,0 +1,103 @@
+"""Denoise speckly sources with a gaussian blur over water (satellite-derived bathymetry).
+
+Allen Coral Atlas's raw 10 m SDB is per-pixel speckly. The merged-DEM smooth (smooth.py) can't
+remove it: that blur is slope-gated to PRESERVE steep gradients, and speckle *is* steep, so the
+noise survives into contour generation — a 2048 px all-water SDB window contours to ~1.55 M
+segments at 1 m intervals (gdal_contour then runs for tens of minutes and blows the sqlite blob
+limit). A small gaussian at the source removes the speckle — which is below SDB's ~±1-2 m vertical
+accuracy anyway, and real shoals span many pixels and survive it — so contours follow real
+bathymetry. Run after source_datum, before source_normalize:
+
+    source_datum <id> ...
+    source_smooth <id> [--sigma 4]
+    source_normalize <id> ...
+
+Measured on that window (25 fine levels): raw 1.55 M features -> sigma=4 ~3.3 k, sigma=8 ~0.6 k.
+A bias-shallow (one-sided) blur was tried first but barely helped (~740 k) — SDB speckle is
+symmetric, so the shallow half alone still explodes. Only water (elevation < 0) is blurred; land
+and nodata clamp to 0 so they don't drag the edge, and the nodata mask is preserved exactly, so
+bounds / coverage are unchanged.
+
+ponytail: whole-tile read (tile_raster pre-splits to <=16k px ~ a few GB); window it like smooth.py
+if a source ships larger single rasters.
+"""
+
+import argparse
+import os
+import sys
+from glob import glob
+
+import numpy as np
+import rasterio
+from scipy.ndimage import gaussian_filter
+
+
+def smooth_array(dem, sigma, nodata, max_depth=None):
+    water = (dem != nodata) & (dem < 0) if nodata is not None else (dem < 0)
+    work = np.where(water, dem, 0.0).astype("float32")  # clamp nodata/land to 0 (smooth.py trick)
+    blur = gaussian_filter(work, sigma=sigma, mode="nearest") if sigma else work  # sigma=0 → cutoff only
+    out = np.where(water, blur, dem).astype(dem.dtype)  # only water blurred; mask + land preserved
+    if max_depth and nodata is not None:
+        # SDB is only reliable to ~1.5x Secchi depth; drop water deeper than the cutoff to nodata so
+        # the merge feathers in the coarse-but-reliable base (GEBCO) below it — the "ACA shallow,
+        # GEBCO deep" blend. Gates on the smoothed value, so it misses deep noise that reads
+        # false-shallow (upgrade to a GEBCO-gated merge, or a depth-weighted feather, for that).
+        out = np.where(water & (out < -max_depth), out.dtype.type(nodata), out)
+    return out
+
+
+def smooth_file(path, sigma, max_depth=None):
+    with rasterio.open(path) as src:
+        profile = src.profile
+        dem = src.read(1)
+        nodata = src.nodata
+    out = smooth_array(dem.astype("float32"), sigma, nodata, max_depth).astype(profile["dtype"])
+    tmp = path + ".smooth.tif"
+    profile.update(driver="GTiff", tiled=True, blockxsize=512, blockysize=512, compress="deflate")
+    with rasterio.open(tmp, "w", **profile) as dst:
+        dst.write(out, 1)
+    os.replace(tmp, path)
+
+
+def main():
+    p = argparse.ArgumentParser(description="Gaussian-denoise (over water) a source's tifs; optionally mask deep water.")
+    p.add_argument("source")
+    p.add_argument("--sigma", type=float, default=4.0, help="gaussian sigma in pixels (default 4; 0 = cutoff only)")
+    p.add_argument("--max-depth", type=float, default=0.0,
+                   help="drop water deeper than this many metres to nodata (0 = off) — masks SDB "
+                        "past its reliable range so the merge fills the coarse base below it")
+    a = p.parse_args()
+    paths = sorted(glob(f"store/source/{a.source}/*.tif"))
+    print(f"{a.source}: gaussian denoise sigma={a.sigma} max_depth={a.max_depth or 'off'} on {len(paths)} file(s)")
+    for path in paths:
+        smooth_file(path, a.sigma, a.max_depth or None)
+
+
+def _check():
+    """Speckle removed over water; land/nodata untouched; mask preserved; deep water cut to nodata."""
+    rng = np.random.default_rng(0)
+    nodata = 0.0
+    dem = (-10 + rng.normal(0, 2, (64, 64))).astype("float32")  # noisy 10 m water
+    dem[:8, :] = nodata   # a nodata band
+    dem[-8:, :] = 5.0     # land (positive)
+    out = smooth_array(dem, sigma=3.0, nodata=nodata)
+    w = (dem != nodata) & (dem < 0)
+    assert out[w].std() < dem[w].std(), (out[w].std(), dem[w].std())  # water denoised
+    assert np.all(out[:8, :] == nodata), "nodata untouched"
+    assert np.all(out[-8:, :] == 5.0), "land untouched"
+    assert np.all((out != nodata) == (dem != nodata)), "nodata mask preserved"
+
+    # max_depth cutoff: a shallow patch is kept, a deep patch is dropped to nodata
+    d2 = np.full((16, 16), -8.0, dtype="float32")   # shallow water
+    d2[8:, :] = -40.0                               # deep water -> should be masked out
+    o2 = smooth_array(d2, sigma=0.0, nodata=nodata, max_depth=18.0)  # sigma=0 → cutoff only
+    assert np.all(o2[:8, :] == -8.0), "shallow water kept"
+    assert np.all(o2[8:, :] == nodata), "water deeper than 18 m dropped to nodata"
+    print("source_smooth.py self-check ok")
+
+
+if __name__ == "__main__":
+    if sys.argv[1:2] == ["--check"]:
+        _check()
+    else:
+        main()

pipelines/test_engine.py

@@ -29,10 +29,11 @@
 
 
 def run(tmp, *args):
-    # Small macrotile_z / num_overviews keep the synthetic rasters tiny.
-    # SKIP_CONTOURS/SKIP_SMOOTH: this is the raster priority test (both have/need none).
+    # Small macrotile_z / num_overviews keep the synthetic rasters tiny. No skip flags — the
+    # value-exact e2e composes only the stages it needs (`aggregation_run.py step reproject merge
+    # tile`, omitting the slope blur) rather than skipping a stage of the full run.
     env = {**os.environ, "SOURCES_DIR": "sources", "PYTHONPATH": PIPE,
-           "MACROTILE_Z": "10", "NUM_OVERVIEWS": "2", "SKIP_CONTOURS": "1", "SKIP_SMOOTH": "1"}
+           "MACROTILE_Z": "10", "NUM_OVERVIEWS": "2"}
     subprocess.run([sys.executable, os.path.join(PIPE, args[0]), *args[1:]],
                    cwd=tmp, env=env, check=True)
 
@@ -212,7 +213,7 @@ def main():
         run(tmp, "source_bounds.py", "base")
         run(tmp, "source_bounds.py", "fine")
         run(tmp, "aggregation_covering.py")
-        run(tmp, "aggregation_run.py")
+        run(tmp, "aggregation_run.py", "step", "reproject", "merge", "tile")
         run(tmp, "downsampling.py", "cover")
         run(tmp, "downsampling.py", "freeze")  # the CI path: shards + tail read this frozen list
         # Exercise the CI fan-out (deep shards + coarse tail), not just the single
@@ -245,7 +246,7 @@ def main():
         # pmtiles across the groups: a tile in no group is missing from the bundles (a
         # hole the Worker overzooms GEBCO into); a tile in two is double-bundled.
         cli_env = {**os.environ, "SOURCES_DIR": "sources", "PYTHONPATH": PIPE,
-                   "MACROTILE_Z": "10", "NUM_OVERVIEWS": "2", "SKIP_CONTOURS": "1", "SKIP_SMOOTH": "1"}
+                   "MACROTILE_Z": "10", "NUM_OVERVIEWS": "2"}
         names = json.loads(subprocess.run(
             [sys.executable, os.path.join(PIPE, "bundle.py"), "groups"],
             cwd=tmp, env=cli_env, check=True, capture_output=True, text=True).stdout.splitlines()[-1])