#!/usr/bin/env python3
"""
THE 60 m BAND AT THREE DEPTHS, side by side, plus the field with no erosion.
amplitude_law.py measured the law: incision is proportional to strength (log-log
slope +1.000) and to erosion_scale (+1.003), and at the map's operating point it
delivers 0.227 m rms / 0.392 m peak. So a target depth in metres is reachable by
scaling strength alone, at the band the map already asks for.
THE TARGET IS RMS, because that is the number every measurement in this thread
has quoted. Peak is printed beside it for every panel, so nothing here depends
on which one you mean by "depth" -- both are on the label.
THE OFF PANEL IS THE CONTROL and it is not decoration: at the map's shipped
0.227 m it should be nearly indistinguishable from the 3 m panel's neighbour,
and if OFF and 3 m look identical then the difference being judged is not the
erosion.
python depth_strip.py
"""
import argparse, os
import numpy as np
from sweep_assumed_slope import render, MAP_PARAMS, EROSION_H, GULLY_METRES
MAP = dict(MAP_PARAMS)
MEAS_RMS, MEAS_STRENGTH = 0.227, MAP["en_strength"] # amplitude_law.py
TARGETS = [3.0, 10.0, 30.0]
def args_with(**over):
d = dict(MAP_PARAMS); d.update(over)
out = []
for k, v in d.items():
out += ["--set", k, str(v)]
return out
def tint(h, lo, hi, step):
t = np.clip((h - lo) / max(hi - lo, 1e-6), 0, 1)
stops = np.array([[0.05,0.12,0.28],[0.20,0.35,0.45],[0.30,0.42,0.22],
[0.55,0.50,0.30],[0.60,0.55,0.52],[0.92,0.92,0.94]])
idx = t * (len(stops) - 1)
i0 = np.floor(idx).astype(int); i1 = np.minimum(i0 + 1, len(stops) - 1)
f = (idx - i0)[..., None]
rgb = stops[i0] * (1 - f) + stops[i1] * f
gy, gx = np.gradient(h, step)
n = 1.0 / np.sqrt(gx*gx + gy*gy + 1.0)
shade = np.clip(0.40 + 0.85 * ((-gx - gy) * n * 0.5 + 0.5), 0, 1.4)
return np.clip(rgb * shade[..., None], 0, 1)
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--bin", default="/tmp/erosion_render")
ap.add_argument("--ox", type=float, default=50000.0)
ap.add_argument("--oz", type=float, default=34000.0)
ap.add_argument("--n", type=int, default=768)
ap.add_argument("--mpp", type=float, default=8.0)
ap.add_argument("--outdir", default="/tmp/depthstrip")
a = ap.parse_args()
os.makedirs(a.outdir, exist_ok=True)
print(f"window {a.n*a.mpp/1000:.2f} km at ({a.ox:.0f}, {a.oz:.0f}) @ {a.mpp} m/px, "
f"gully band {GULLY_METRES:g} m ({GULLY_METRES/a.mpp:.1f} px per gully)")
off = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/ds_off.f32",
args_with(en_strength=0))
panels = [("erosion OFF", off, 0.0, 0.0, 0.0)]
shipped = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/ds_ship.f32", args_with())
e = shipped - off
panels.append((f"SHIPPED (asks 30 m)", shipped, MEAS_STRENGTH,
float(np.sqrt((e**2).mean())), float(np.abs(e).max())))
for tgt in TARGETS:
s = MEAS_STRENGTH * (tgt / MEAS_RMS)
h = render(a.bin, a.ox, a.oz, a.n, a.mpp, f"/tmp/ds_{tgt:g}.f32",
args_with(en_strength=s))
e = h - off
panels.append((f"{tgt:g} m", h, s,
float(np.sqrt((e**2).mean())), float(np.abs(e).max())))
print(f"\n{'panel':>20} {'strength':>12} {'rms m':>9} {'peak m':>9} {'d:w':>7}")
for name, _, s, rms, pk in panels:
dw = pk / GULLY_METRES if pk else 0.0
print(f"{name:>20} {s:>12.6g} {rms:>9.3f} {pk:>9.3f} {dw:>7.3f}")
lo = min(p[1].min() for p in panels); hi = max(p[1].max() for p in panels)
imgs = [tint(p[1], lo, hi, a.mpp) for p in panels]
strip = np.concatenate(imgs, axis=1)
from PIL import Image, ImageDraw
im = Image.fromarray((strip*255).astype(np.uint8))
d = ImageDraw.Draw(im)
for k, (name, _, s, rms, pk) in enumerate(panels):
x = k * a.n + 8
d.rectangle([x-6, 4, x + 250, 46], fill=(0, 0, 0))
d.text((x, 8), f"{name}", fill=(255, 255, 255))
d.text((x, 22), f"rms {rms:.2f} m peak {pk:.2f} m", fill=(210, 210, 210))
d.text((x, 34), f"strength {s:.5g}", fill=(160, 160, 160))
Image.fromarray((imgs[k]*255).astype(np.uint8)).save(
os.path.join(a.outdir, f"panel{k}_{name.split()[0]}.png"))
out = os.path.join(a.outdir, "strip.png")
im.save(out)
print(f"\nwrote {out} ({im.width}x{im.height}, shared tint {lo:.0f}..{hi:.0f} m)")
if __name__ == "__main__":
main()