Add WSV DGM-W (German federal waterway bathymetry) source

pipelines/source_datum.py

@@ -1,11 +1,18 @@
-"""Apply the bathymetry value transform: ``negate`` then ``datum_offset_m``.
+"""Apply the bathymetry value transform: ``negate``, ``datum_offset_m``, then an
+optional ``offset_surface`` (a spatially-varying reference subtraction).
 
-Reads the knobs from ``metadata.json``:
+Reads the knobs from ``metadata.json`` / recipe args:
 
   - ``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
-    offset, not full VDatum — the Milestone 3 seam for a spatially-varying separation).
+    offset — the flat-datum case, e.g. a lake surface level).
+  - ``offset_surface``: a reference raster (height in the source's vertical frame)
+    subtracted per pixel, for sources whose datum separation varies in space rather
+    than being one constant. Turns ``bed`` into ``bed - reference`` — e.g. an NHN
+    riverbed minus the BSH SKN-Fläche (chart datum in NHN) → depth below chart datum.
+    Where the reference does not cover a pixel, that pixel becomes nodata (it cannot
+    be referenced), so partial-coverage references only touch what they cover.
 
 Operates per file in ``store/source/<id>/``, in the source's native CRS (the
 reprojection to Web Mercator happens later, in the aggregation stage), preserving
@@ -22,9 +29,27 @@
 
 import numpy as np
 import rasterio
-
-
-def transform_file(filepath, negate, offset, clamp_positive=False):
+from rasterio.warp import Resampling, reproject
+
+
+def _surface_on_grid(surface_path, profile):
+    """Resample a reference raster onto the tile grid (profile's crs/transform/shape),
+    bilinear (it's a smooth continuous surface). Returns (values float32, valid bool);
+    valid is False where the reference has no data over the tile."""
+    dst = np.full((profile["height"], profile["width"]), np.nan, dtype="float32")
+    with rasterio.open(surface_path) as ref:
+        reproject(
+            source=rasterio.band(ref, 1),
+            destination=dst,
+            src_transform=ref.transform, src_crs=ref.crs,
+            dst_transform=profile["transform"], dst_crs=profile["crs"],
+            src_nodata=ref.nodata, dst_nodata=np.nan,
+            resampling=Resampling.bilinear,
+        )
+    return dst, np.isfinite(dst)
+
+
+def transform_file(filepath, negate, offset, clamp_positive=False, offset_surface=None):
     with rasterio.open(filepath) as src:
         profile = src.profile
         data = src.read(1)
@@ -37,14 +62,27 @@ def transform_file(filepath, negate, offset, clamp_positive=False):
     if offset:
         valid = valid + np.float32(offset)
     data[mask] = valid
+
+    if offset_surface:
+        nodata = profile.get("nodata")
+        if nodata is None:
+            raise ValueError(f"{filepath}: --offset-surface needs a nodata value set")
+        ref, ref_valid = _surface_on_grid(offset_surface, profile)
+        # bed - reference where both cover; a bed pixel with no reference cannot be
+        # referenced to the datum, so drop it to nodata rather than leave it at NHN.
+        both = mask & ref_valid
+        data[mask & ~ref_valid] = np.float32(nodata)
+        data[both] = data[both] - ref[both]
+
     if clamp_positive:
-        # After the offset, 0 = water surface; anything > 0 is the surrounding terrain
-        # (a lake DEM's land fringe, or a topobathy playa) — drop it to nodata so it
-        # can't bleed into the water layer as false land.
+        # After the transform, 0 = the reference (water surface / chart datum); anything
+        # > 0 is above it (surrounding terrain, a lake DEM's land fringe, a topobathy
+        # playa) — drop it to nodata so it can't bleed into the water layer as false land.
         nodata = profile.get("nodata")
         if nodata is None:
             raise ValueError(f"{filepath}: --clamp-positive needs a nodata value set")
-        data[mask & (data > 0)] = np.float32(nodata)
+        cur = data != np.float32(nodata) if nodata is not None else mask
+        data[mask & cur & (data > 0)] = np.float32(nodata)
 
     profile.update(driver="GTiff", dtype="float32", tiled=True,
                    blockxsize=512, blockysize=512, compress="deflate")
@@ -54,44 +92,60 @@ def transform_file(filepath, negate, offset, clamp_positive=False):
     os.replace(tmp, filepath)
 
 
-def write_sidecar(source, negate, offset, clamp_positive):
+def write_sidecar(source, negate, offset, clamp_positive, offset_surface=None):
     """Record the applied transform in store/source/<id>/datum.json — the machine-readable
     provenance source_catalog folds into the catalog item (vertical-datum offset was invisible
     downstream when it lived only in this CLI arg). Written whenever the step runs, so a source
     whose recipe calls source_datum always leaves a sidecar."""
     os.makedirs(f"store/source/{source}", exist_ok=True)
     with open(f"store/source/{source}/datum.json", "w") as f:
         json.dump({"negate": bool(negate), "offset_m": float(offset),
-                   "clamp_positive": bool(clamp_positive)}, f, indent=2)
+                   "clamp_positive": bool(clamp_positive),
+                   "offset_surface": os.path.basename(offset_surface) if offset_surface else None}, f, indent=2)
+
+
+def _resolve_surface(source, arg):
+    """Accept a path or a filename cached under store/source/<id>/."""
+    if os.path.exists(arg):
+        return arg
+    candidate = f"store/source/{source}/{arg}"
+    if os.path.exists(candidate):
+        return candidate
+    sys.exit(f"--offset-surface: {arg} not found (also tried {candidate})")
 
 
 def main():
-    p = argparse.ArgumentParser(description="Apply negate + datum offset to a source's tifs.")
+    p = argparse.ArgumentParser(description="Apply negate + datum offset (scalar and/or surface) to a source's tifs.")
     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("--offset", type=float, default=0.0, help="metres added to reach ~MSL (flat datum)")
+    p.add_argument("--offset-surface", dest="offset_surface", default=None,
+                   help="reference raster (height in the source frame) subtracted per pixel — "
+                        "the spatially-varying datum case (e.g. BSH SKN-Fläche for chart datum). "
+                        "A path, or a filename under store/source/<id>/.")
     p.add_argument("--clamp-positive", action="store_true",
-                   help="after the offset, drop cells > 0 (above the water surface) to nodata — "
+                   help="after the transform, drop cells > 0 (above the reference) to nodata — "
                         "removes a lake DEM's land fringe / a topobathy playa")
     a = p.parse_args()
 
+    surface = _resolve_surface(a.source, a.offset_surface) if a.offset_surface else None
+
     # Record what this invocation applies even when it's a no-op, so the sidecar exists for
     # every source whose recipe runs source_datum (source_catalog's invariant).
-    write_sidecar(a.source, a.negate, a.offset, a.clamp_positive)
+    write_sidecar(a.source, a.negate, a.offset, a.clamp_positive, surface)
 
-    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 not a.clamp_positive and not surface:
+        print(f"{a.source}: no datum transform (negate=False, offset=0, no surface)")
         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}: negate={a.negate} offset={a.offset} surface={os.path.basename(surface) if surface else None} "
+          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, surface)
 
 
 def _check():
-    """Self-check the value transform on a synthetic raster (no GDAL CLI)."""
-    import os
+    """Self-check the value transform on synthetic rasters (no GDAL CLI)."""
     import tempfile
     from rasterio.transform import from_origin
 
@@ -107,7 +161,6 @@ def _check():
     transform_file(path, negate=True, offset=-1.0)  # depth->elev, then -1 m datum
     with rasterio.open(path) as src:
         out = src.read(1)
-    # valid pixels: -(v) - 1 ; nodata untouched
     assert out[0, 0] == -6.0 and out[0, 1] == -11.0, out
     assert out[1, 0] == -1.0 and out[1, 1] == -3.5 and out[1, 2] == -101.0, out
     assert out[0, 2] == nodata, out[0, 2]  # nodata not negated into +9999
@@ -124,6 +177,31 @@ 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
+
+    # offset_surface: bed (NHN) minus a chart-datum surface (NHN), partial coverage.
+    # bed grid: 2x2 at NHN; nodata one cell. Reference covers only the left column.
+    bedp = os.path.join(d, "bed.tif")
+    bed = np.array([[-10.0, -8.0], [-6.0, nodata]], dtype="float32")
+    with rasterio.open(bedp, "w", driver="GTiff", height=2, width=2, count=1,
+                       dtype="float32", nodata=nodata, crs="EPSG:4326",
+                       transform=from_origin(0, 2, 1, 1)) as dst:
+        dst.write(bed, 1)
+    # reference: SKN in NHN ~ -1.5 over the left half, nodata over the right half
+    refp = os.path.join(d, "ref.tif")
+    refnd = 99999.0
+    ref = np.array([[-1.5, refnd], [-1.5, refnd]], dtype="float32")
+    with rasterio.open(refp, "w", driver="GTiff", height=2, width=2, count=1,
+                       dtype="float32", nodata=refnd, crs="EPSG:4326",
+                       transform=from_origin(0, 2, 1, 1)) as dst:
+        dst.write(ref, 1)
+    transform_file(bedp, negate=False, offset=0.0, offset_surface=refp)
+    with rasterio.open(bedp) as src:
+        ob = src.read(1)
+    # left column referenced: bed - (-1.5) = bed + 1.5 ; right column no reference -> nodata
+    assert abs(ob[0, 0] - (-8.5)) < 1e-4, ob        # -10 - (-1.5)
+    assert abs(ob[1, 0] - (-4.5)) < 1e-4, ob        # -6  - (-1.5)
+    assert ob[0, 1] == nodata, ob                    # covered bed but no reference -> dropped
+    assert ob[1, 1] == nodata, ob                    # bed nodata stays nodata
     print("source_datum.py self-check ok")
 
 

sources/README.md

@@ -22,6 +22,7 @@ Selection rule: **resolution sets the zoom cap, an openly-redistributable licens
 | [Vaklodingen](vaklodingen/)                                        | 20 m        | z12       | Netherlands                         | NAP (~MSL)              |
 | INFOMAR ([10 m](infomar_10m/), [25 m](infomar_25m/))               | 10 m / 25 m | z13 / z11 | Ireland inshore + shelf             | **LAT**                 |
 | [UK SurfZone](uk_surfzone/)                                        | 2 m         | z14       | England intertidal                  | ODN (~MSL)              |
+| [WSV DGM-W (Germany)](dgm_w/)                                       | 2 m         | z14       | German tidal estuaries (Elbe, Weser, Jade, Nordsee) | **SKN/LAT** (chart datum; inland reaches deferred) |
 | [GSC Atlantic](gsc_atlantic/)                                      | 100 m       | ~z10      | Scotian Shelf + NL                  | unverified              |
 | [GSC Pacific](gsc_pacific/)                                        | 10 m        | z13       | BC coast + Salish Sea               | unverified              |
 | [gbr30](gbr30/)                                                    | 30 m        | z12       | GBR + Coral Sea                     | MSL                     |

sources/dgm_w/Justfile

@@ -0,0 +1,34 @@
+# German Federal Waterway Bathymetry — WSV DGM-W (Digitales Geländemodell des Wasserlaufs), 2 m.
+# Run from pipelines/:  just ../sources/dgm_w/
+# Zipped GeoTIFF tiles from the WSV INSPIRE "Höhe" Atom download service (file_list.txt), one .tif
+# per zip, terrain + riverbed merged. Per-tile CRS (DHDN/Gauss-Krüger and UTM/ETRS89 zones) ->
+# mixed_crs: normalize WITHOUT --crs to keep each file's embedded CRS; aggregation reprojects each
+# tile from its own frame. Raw values are orthometric NHN elevation; nodata -32768.
+#
+# SCOPE: active reaches in file_list.txt — tidal estuaries (16 tiles: Nordsee, Außenelbe below
+# Elbe-km 586, Jade, Unter-/Außenweser), the free-flowing Rhein (3) and upper Elbe (25), and the
+# impounded canalised rivers Main (6), upper Rhein (16), Mosel (1), Saar (5), Lahn (3). The other
+# inland reaches stay commented out (bed above MSL, pending low-water referencing).
+#
+# DATUM: build_reference.py builds a per-reach low-water surface in NHN — SKN (chart datum ~= LAT)
+# for the tidal estuaries, a GlW/MNW ramp for the free-flowing Rhein + Elbe, per-pool Stauziel for
+# the impounded canalised rivers — into store/source/dgm_w/reference/. source_datum subtracts it (bed - datum =
+# depth) and --clamp-positive drops above-datum land/drying to nodata, so the synced COG is
+# low-water normalized. Licenses: DEM GeoNutzV (attribute WSV), BSH SKN-Fläche CC-BY 4.0 (attribute
+# BSH). Not for navigation.
+[no-cd]
+default: download normalize
+
+[no-cd]
+download:
+    uv run python source_download.py dgm_w
+
+[no-cd]
+normalize:
+    uv run python source_unzip.py dgm_w
+    uv run python ../sources/dgm_w/build_reference.py
+    uv run python source_datum.py dgm_w --offset-surface reference/reference.vrt --clamp-positive
+    uv run python source_normalize.py dgm_w --nodata -32768
+    uv run python source_bounds.py dgm_w
+    uv run python source_polygonize.py dgm_w 8
+    uv run python source_create_tarball.py dgm_w