diff --git a/ROADMAP.md b/ROADMAP.md index 1e37be80..003babfe 100644 --- a/ROADMAP.md +++ b/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) diff --git a/pipelines/aggregation_run.py b/pipelines/aggregation_run.py index 148bc55a..65c8f8b9 100644 --- a/pipelines/aggregation_run.py +++ b/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 process the frozen dirty[i::n] (matrix shard i of n) aggregation_run.py matrix print the shard matrix JSON (sized to the dirt) + aggregation_run.py step ... 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 | matrix ]") + sys.exit("usage: aggregation_run.py [freeze | shard | matrix | step ...]") if __name__ == "__main__": diff --git a/pipelines/smooth.py b/pipelines/smooth.py index 5849590b..05d87e3f 100644 --- a/pipelines/smooth.py +++ b/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 diff --git a/pipelines/source_datum.py b/pipelines/source_datum.py index d5a5d387..3a5e49d2 100644 --- a/pipelines/source_datum.py +++ b/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,7 +24,7 @@ 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) @@ -30,6 +32,8 @@ def transform_file(filepath, negate, offset, clamp_positive=False): 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") diff --git a/pipelines/source_polygonize.py b/pipelines/source_polygonize.py index 91592247..4ee310d8 100644 --- a/pipelines/source_polygonize.py +++ b/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) diff --git a/pipelines/source_smooth.py b/pipelines/source_smooth.py new file mode 100644 index 00000000..239c6adc --- /dev/null +++ b/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 ... + source_smooth [--sigma 4] + source_normalize ... + +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() diff --git a/pipelines/test_engine.py b/pipelines/test_engine.py index afeba98d..be2c337c 100644 --- a/pipelines/test_engine.py +++ b/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]) diff --git a/pipelines/tile_raster.py b/pipelines/tile_raster.py new file mode 100644 index 00000000..289b6b95 --- /dev/null +++ b/pipelines/tile_raster.py @@ -0,0 +1,105 @@ +"""Tile large source rasters in store/source// into smaller tiles, dropping all-nodata tiles. + +A recipe step for sources whose raw download is one raster too large for the per-file +`source_*` steps to hold in memory (`source_datum` reads a whole band via rasterio; +`source_polygonize` via gdal_calc). A ~23 GB Allen Coral Atlas regional GeoTIFF or the ~32 GB +GSC Pacific DEM would OOM. Run it after `source_download`, before `source_datum`: + + source_download + tile_raster # split each store/source//*.tif into tiles, drop empties + source_datum ... + +Memory-safe: each tile is cut with a windowed `gdal_translate -srcwin` (block streaming) and +emptiness is tested on a windowed mask read — the whole raster is never materialized. Source +CRS / nodata / dtype carry over (the values stay PRISTINE; `source_normalize` re-asserts CRS+nodata, +`source_datum` applies any unit/datum transform). Each input raster is replaced by its tiles +(`tile___.tif`); already-tiled outputs are skipped so a rerun is idempotent. + +ponytail: gdal_translate per tile (not gdal_retile) so emptiness is judged before a tile is +written — no all-nodata tiles created then pruned. +""" + +import argparse +import os +import subprocess +import sys +from glob import glob + +import rasterio +from rasterio.windows import Window + + +def tile_file(path, size): + """Split one raster into tile__RR_CC.tif beside it; return the kept tile names.""" + base = os.path.splitext(os.path.basename(path))[0] + out_dir = os.path.dirname(path) or "." + with rasterio.open(path) as src: + width, height = src.width, src.height + ncols = (width + size - 1) // size + nrows = (height + size - 1) // size + kept = [] + for r in range(nrows): + for c in range(ncols): + col_off, row_off = c * size, r * size + w = min(size, width - col_off) + h = min(size, height - row_off) + with rasterio.open(path) as src: # reopen per tile so only this window is read + if not src.read_masks(1, window=Window(col_off, row_off, w, h)).any(): + continue # all-nodata tile (open ocean) — skip + name = f"tile_{base}_{r:02d}_{c:02d}.tif" + subprocess.run( + ["gdal_translate", "-q", "-of", "GTiff", + "-srcwin", str(col_off), str(row_off), str(w), str(h), + "-co", "TILED=YES", "-co", "COMPRESS=ZSTD", "-co", "BIGTIFF=IF_NEEDED", + path, os.path.join(out_dir, name)], + check=True) + kept.append(name) + return kept + + +def tile_source(source, size): + paths = [p for p in sorted(glob(f"store/source/{source}/*.tif")) + if not os.path.basename(p).startswith("tile_")] # skip our own outputs (idempotent rerun) + print(f"{source}: tiling {len(paths)} file(s) at {size}px") + for p in paths: + kept = tile_file(p, size) + os.remove(p) # replace the monolith with its tiles + print(f" {os.path.basename(p)} -> {len(kept)} non-empty tile(s)") + + +def main(): + p = argparse.ArgumentParser(description="Tile large rasters in store/source//, dropping all-nodata tiles.") + p.add_argument("source") + p.add_argument("--size", type=int, default=16384, help="tile edge in pixels (default 16384)") + a = p.parse_args() + tile_source(a.source, a.size) + + +def _check(): + """Self-check: an all-nodata quadrant is dropped; data tiles keep their values.""" + import tempfile + import numpy as np + from rasterio.transform import from_origin + + d = tempfile.mkdtemp() + path = os.path.join(d, "big.tif") + arr = np.zeros((4, 4), dtype="int16") # 2x2 grid of 2x2 tiles, nodata=0 + arr[0:2, 0:2] = [[5, 10], [3, 0]] # only the top-left tile has valid data + with rasterio.open(path, "w", driver="GTiff", height=4, width=4, count=1, dtype="int16", + nodata=0, crs="EPSG:4326", transform=from_origin(0, 4, 1, 1)) as dst: + dst.write(arr, 1) + + kept = tile_file(path, size=2) + assert kept == ["tile_big_00_00.tif"], kept # other three tiles are all-nodata -> dropped + with rasterio.open(os.path.join(d, "tile_big_00_00.tif")) as t: + o = t.read(1) + assert t.nodata == 0, t.nodata # nodata carried over + assert o[0, 0] == 5 and o[0, 1] == 10 and o[1, 0] == 3, o + print("tile_raster.py self-check ok") + + +if __name__ == "__main__": + if sys.argv[1:2] == ["--check"]: + _check() + else: + main() diff --git a/sources/aca_ncaribbean/Justfile b/sources/aca_ncaribbean/Justfile new file mode 100644 index 00000000..c85c9994 --- /dev/null +++ b/sources/aca_ncaribbean/Justfile @@ -0,0 +1,17 @@ +# Allen Coral Atlas SDB (N. Caribbean, Florida & Bahamas) — 10 m, CC-BY 4.0. Run from pipelines/: +# just ../sources/aca_ncaribbean/ +# R2 MIRROR of the raw ACA download (account-gated behind an expiring signed GCS URL → not CI-fetchable; +# mirrored once like BATNAS). One ~23 GB Int16 GeoTIFF, too big for the per-file steps to hold in memory, +# so tile_raster splits it into per-tile rasters (dropping all-nodata ocean tiles) right after download. +# Tiles are pristine source values: Int16 CENTIMETRES, positive-down depth, NoData 0, EPSG:4326. +# SDB only resolves shallow water (~0–20 m where the bottom is visible); everything deeper is NoData → GEBCO. +[no-cd] +default: + uv run python source_download.py aca_ncaribbean + uv run python tile_raster.py aca_ncaribbean + uv run python source_datum.py aca_ncaribbean --negate --scale 0.01 # turns cm depth → metres + uv run python source_smooth.py aca_ncaribbean --sigma 8 --max-depth 18 # denoise + drop water >18 m (SDB's reliable limit) → GEBCO fills the deep via the merge + uv run python source_normalize.py aca_ncaribbean --crs EPSG:4326 + uv run python source_bounds.py aca_ncaribbean + uv run python source_polygonize.py aca_ncaribbean 8 + uv run python source_create_tarball.py aca_ncaribbean diff --git a/sources/aca_ncaribbean/file_list.txt b/sources/aca_ncaribbean/file_list.txt new file mode 100644 index 00000000..291e95d7 --- /dev/null +++ b/sources/aca_ncaribbean/file_list.txt @@ -0,0 +1,6 @@ +# Allen Coral Atlas SDB (N. Caribbean, Florida & Bahamas), 10 m, CC-BY 4.0 — our R2 MIRROR of the raw +# ACA download (account-gated behind an expiring signed GCS URL → not CI-fetchable; mirrored once like +# BATNAS). One ~23 GB Int16 GeoTIFF: CENTIMETRES, positive-down depth, NoData 0, EPSG:4326. The recipe +# tiles it (tile_raster — too big for the per-file steps to hold in memory), then source_datum --negate +# --scale 0.01 turns cm depth → metres elevation. Refresh = re-download from ACA, re-sync to R2 (rclone). +https://data.openwaters.io/bathymetry/mirror/aca_ncaribbean/ACA_NCaribbean_Florida_Bahamas_SDB_10m_cm.tif diff --git a/sources/aca_ncaribbean/metadata.json b/sources/aca_ncaribbean/metadata.json new file mode 100644 index 00000000..e00109b0 --- /dev/null +++ b/sources/aca_ncaribbean/metadata.json @@ -0,0 +1,8 @@ +{ + "name": "Allen Coral Atlas Satellite-Derived Bathymetry (N. Caribbean, Florida & Bahamas, 10 m)", + "producer": "Allen Coral Atlas / Arizona State University", + "website": "https://allencoralatlas.org/", + "license": "CC-BY 4.0", + "max_zoom": 11, + "datum": "~MSL (instantaneous sea surface; satellite-derived)" +} diff --git a/sources/cudem/metadata.json b/sources/cudem/metadata.json index bb1111e4..8f348436 100644 --- a/sources/cudem/metadata.json +++ b/sources/cudem/metadata.json @@ -4,5 +4,6 @@ "website": "https://www.ncei.noaa.gov/products/coastal-relief-model", "license": "public domain (U.S. Government work)", "max_zoom": 13, - "volatile": true + "volatile": true, + "priority": 1 } diff --git a/sources/cudem_third/metadata.json b/sources/cudem_third/metadata.json index fb036622..d7ac81f7 100644 --- a/sources/cudem_third/metadata.json +++ b/sources/cudem_third/metadata.json @@ -4,5 +4,6 @@ "website": "https://www.ncei.noaa.gov/products/coastal-relief-model", "license": "public domain (U.S. Government work)", "max_zoom": 12, - "volatile": true + "volatile": true, + "priority": 1 } diff --git a/sources/noaa_s102/metadata.json b/sources/noaa_s102/metadata.json index ebf9be8e..9be2f9a8 100644 --- a/sources/noaa_s102/metadata.json +++ b/sources/noaa_s102/metadata.json @@ -9,5 +9,5 @@ "volatile": true, "negate": true, "link_column": "S102V30", - "priority": 1 + "priority": 2 }