"""Bellhop bathymetry ``.bty`` reader (port of Matlab ``readbty.m`` core file I/O)."""
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
def _quoted_inner(line: str) -> str:
"""First single-quoted substring from a ``.bty`` header line."""
parts = line.split("'")
if len(parts) < 3:
raise ValueError(f"expected quoted bathymetry type on line: {line!r}")
return parts[1]
def _pad_bty_type(inner: str) -> str:
"""Normalize type token to two characters."""
return (inner + " ")[:2]
def _collect_floats_after(lines: list[str], start: int) -> list[float]:
"""Scan lines from ``start`` for comma/space-separated floats (skip ``!``, ``/``)."""
out: list[float] = []
for j in range(start, len(lines)):
s = lines[j].strip()
if not s:
continue
if "!" in s:
s = s.split("!", 1)[0].strip()
for tok in s.replace(",", " ").split():
if tok == "/":
continue
try:
out.append(float(tok))
except ValueError:
continue
return out
[docs]
@dataclass(frozen=True)
class BathymetryRead:
"""2D bathymetry from a ``.bty`` file (ranges in **m** after km→m)."""
interp_type: str
"""``L`` piecewise-linear or ``C`` curvilinear (first letter)."""
range_m: np.ndarray
depth_m: np.ndarray
tangent_u: np.ndarray
"""Unit tangent along each segment, shape ``(2, nseg)`` (range, depth)."""
normal_u: np.ndarray
"""Outward unit normal, shape ``(2, nseg)``."""
segment_length_m: np.ndarray
segment_curvature: np.ndarray
"""Curvature ``dφ/ds`` on each segment; zero for ``L``."""
tangent_node_u: np.ndarray | None
"""Node tangents for ``C`` (curvilinear); ``None`` for ``L``."""
normal_node_u: np.ndarray | None
"""Node normals for ``C``; ``None`` for ``L``."""
[docs]
def parse_bty(
text: str,
*,
max_depth_m: float | None = None,
) -> BathymetryRead:
"""Parse a 2D Bellhop ``.bty`` file (Matlab ``readbty``).
Ranges in the file are **km**; returned ``range_m`` uses **m** (×1000), matching
``readbty``. Endpoints at ±1e50 m are appended as in Matlab.
If ``max_depth_m`` is set (e.g. last SSP depth), bathymetry deeper than that raises
``ValueError``, matching ``readbty``.
"""
lines = text.splitlines()
i = 0
while i < len(lines) and not lines[i].strip():
i += 1
if i >= len(lines):
raise ValueError("empty .bty file")
bty_head = _pad_bty_type(_quoted_inner(lines[i]))
i += 1
interp = bty_head[0:1]
if interp not in ("L", "C"):
raise ValueError(f"unknown bathymetry interpolation type: {interp!r}")
while i < len(lines) and not lines[i].strip():
i += 1
if i >= len(lines):
raise ValueError(".bty: missing point count")
n_line = lines[i].strip()
if "!" in n_line:
n_line = n_line.split("!", 1)[0].strip()
n_pts = int(float(n_line.split()[0].rstrip(",")))
i += 1
floats = _collect_floats_after(lines, i)
need = 2 * n_pts
if len(floats) < need:
raise ValueError(f".bty: expected {need} floats for {n_pts} points, got {len(floats)}")
raw_r_km = np.asarray(floats[0::2][:n_pts], dtype=np.float64)
raw_z = np.asarray(floats[1::2][:n_pts], dtype=np.float64)
if max_depth_m is not None and np.any(raw_z > max_depth_m):
raise ValueError("bathymetry goes below last tabulated SSP value")
# Extend to ±infinity (km → m for ranges)
r_m = np.empty(n_pts + 2, dtype=np.float64)
z = np.empty(n_pts + 2, dtype=np.float64)
r_m[0] = -1e50
r_m[1:-1] = raw_r_km * 1000.0
r_m[-1] = 1e50
z[0] = raw_z[0]
z[1:-1] = raw_z
z[-1] = raw_z[-1]
n_ext = n_pts + 2
dr = np.diff(r_m)
dz = np.diff(z)
seg_len = np.sqrt(dr**2 + dz**2)
tx = dr / seg_len
tz = dz / seg_len
t_bot = np.vstack([tx, tz])
n_bot = np.vstack([-tz, tx])
n_seg = n_ext - 1
curv_seg = np.zeros(n_seg, dtype=np.float64)
t_node: np.ndarray | None = None
n_node: np.ndarray | None = None
if interp == "C":
t_node = np.zeros((2, n_ext), dtype=np.float64)
n_node = np.zeros((2, n_ext), dtype=np.float64)
for ii in range(1, n_ext - 1):
t_node[:, ii] = 0.5 * (t_bot[:, ii - 1] + t_bot[:, ii])
t_node[:, 0] = [1.0, 0.0]
t_node[:, -1] = [1.0, 0.0]
# readbty.m: normals at nodes from averaged tangents (not averaged normals)
n_node[0, :] = -t_node[1, :]
n_node[1, :] = t_node[0, :]
phi = np.arctan2(t_node[1, :], t_node[0, :])
curv_seg = np.diff(phi) / seg_len
return BathymetryRead(
interp_type=interp,
range_m=r_m,
depth_m=z,
tangent_u=t_bot,
normal_u=n_bot,
segment_length_m=seg_len,
segment_curvature=curv_seg,
tangent_node_u=t_node,
normal_node_u=n_node,
)
[docs]
def parse_bty_bytes(
data: bytes,
*,
encoding: str = "utf-8",
max_depth_m: float | None = None,
) -> BathymetryRead:
"""Like :func:`parse_bty` after decoding ``data`` with ``encoding``."""
return parse_bty(data.decode(encoding, errors="replace"), max_depth_m=max_depth_m)
def _fmt_bty_float(x: float) -> str:
"""Format a float for stable ``.bty`` round-trip."""
return f"{x:.17g}"