#!/usr/bin/env python3
"""
BAND COMPARISON — terrain-diffusion's spectrum against our erosion filter's,
on the same physical axes.
hyper-games #46 §3, probe B. Both fields are sampled over the same physical
width at the same metres-per-pixel, so a given radial bin k is the same
wavelength in both and the two spectra can be read against each other directly.
The instrument is the same one, and its self-test runs first.
"""
import argparse, sys
import numpy as np
sys.path.insert(0, '/storage/projects/code/godot/terrain-lookdev')
from spectrum_probe import fetch_tile, radial_power_spectrum, selftest
def octave_table(name, field, mpp, kmin, kmax):
n = field.shape[0]
width = n * mpp
k, p = radial_power_spectrum(field)
rows = []
lo = kmin
while lo * 2 <= k.max():
hi = min(lo * 2, int(k.max()))
m = (k > lo) & (k <= hi) & (p > 0)
if m.sum() >= 3:
# local power-law slope inside this octave: the honest measure of
# "how rough is it HERE", with no extrapolation from a lower band.
s, _ = np.polyfit(np.log10(k[m]), np.log10(p[m]), 1)
rows.append((width / hi, width / lo, -s, m.sum()))
lo = hi
return k, p, rows
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--url", default="http://127.0.0.1:8817")
ap.add_argument("--erosion", default="/tmp/erosion-512.f32")
ap.add_argument("--n", type=int, default=512)
ap.add_argument("--mpp", type=float, default=30.0)
ap.add_argument("--i1", type=int, default=0)
ap.add_argument("--j1", type=int, default=2048)
args = ap.parse_args()
print("INSTRUMENT SELF-TEST")
if not selftest():
print(" FAILED — no number below is meaningful.")
return 2
print(" all pass\n")
n, mpp = args.n, args.mpp
print(f"Both fields: {n}x{n} px @ {mpp:g} m/px = {n*mpp/1000:.2f} km across, "
f"Nyquist lambda = {2*mpp:g} m\n")
ero = np.fromfile(args.erosion, dtype="<f4").reshape(n, n).astype(np.float64)
nan = int(np.isnan(ero).sum())
if nan:
print(f" NOTE: {nan} NaN cells in the erosion field; filled with the mean.")
ero = np.where(np.isnan(ero), np.nanmean(ero), ero)
dif, _ = fetch_tile(args.url, args.i1, args.j1, n, 1)
fields = [("terrain-diffusion (30m model, seed 1234)", dif),
("our erosion filter (map_dc_erosion, band-limited)", ero)]
for label, f in fields:
k, p, rows = octave_table(label, f, mpp, 4, n // 2)
print(f"{label}")
print(f" relief: min {f.min():.0f} m max {f.max():.0f} m "
f"sd {f.std():.0f} m")
print(f" {'octave (wavelength)':>26} {'local beta':>11} {'bins':>6}")
for w_lo, w_hi, beta, nb in rows:
unit = (f"{w_lo/1000:5.2f}-{w_hi/1000:5.2f} km" if w_hi >= 1000
else f"{w_lo:6.0f}-{w_hi:6.0f} m ")
print(f" {unit:>26} {beta:11.2f} {nb:6d}")
print()
print("READING IT: beta is the local 2-D radial power-law exponent inside")
print("that octave alone -- no extrapolation. Higher beta = smoother, i.e.")
print("less energy at the fine end of the octave than the coarse end.")
print("A self-similar generator holds beta constant down the column.")
print("Subtract 1 for the radially-INTEGRATED convention some papers use.")
return 0
if __name__ == "__main__":
sys.exit(main())