waveray

Fast last-stage nearshore spectral wave transformation: precomputed backward ray-traced linear transfer operators over local bathymetry, with parametric depth-limited wave breaking at the target.

📖 Documentation: https://oceanum.github.io/waveray/

Built to downscale directional wave spectra from SWAN (or WW3) hindcasts and nowcasts through the final nearshore transformation to a target site — decades of hourly spectra in seconds, no wave model in the runtime loop.

Install

pip install waveray                # core
pip install "waveray[datamesh]"    # + Oceanum Datamesh bathymetry/spectra access

Documentation

Full docs: https://oceanum.github.io/waveray/

The pages are plain Markdown in docs/ and render on GitHub too.

Method

Setup, once per site:

1. Pull local bathymetry (e.g. from an Oceanum Datamesh datasource) into a LocalGrid on a local tangent plane.
2. For every (frequency, direction) bin, trace ray bundles backward from the target over the 2D bathymetry (O'Reilly & Guza style spectral refraction). Each sub-ray records where and in what direction it exits the domain, or whether it runs aground (island/headland sheltering).
3. Assemble a linear operator T[f, dir_t, bp, dir_b] mapping K boundary spectra (the SWAN output points around the domain perimeter — alongshore inhomogeneity is handled by exit-point interpolation) to the target spectrum. The per-ray coefficient is the exact linear invariant E(f, theta) * c * cg = const, which reproduces Snell refraction and energy-flux shoaling identically (validated against analytic plane-beach solutions in the test suite).

Along each ray path the coefficients also integrate the linear source terms: JONSWAP bottom friction (on by default) and, when a wind is supplied, SWAN's wind input — the Komen et al. (1984) exponential growth folded into the same path exponent, plus (optionally, agrow=True) the Cavaleri & Malanotte-Rizzoli (1981) linear growth accumulated into an additive spectrum that seeds locally generated wind sea. Wind can be a uniform (speed, direction) or a gridded snapshot, with u* from the Zijlema et al. (2012) drag law (SWAN's default).

Runtime, per hindcast:

4. E_t = einsum(T, E_b) over all timesteps at once.
5. Depth-limited breaking as a nonlinear post-step: transformed Hm0 is capped at a Miche-type (Battjes-Janssen) or gamma * d limit with optional tide-modulated depth; the cap scales the spectrum proportionally, matching how SWAN distributes surf-breaking dissipation over the spectrum.

The stationary, linear-propagation assumption (no temporal nonlinear evolution, no quadruplets/triads) is deliberate: over a last-mile nearshore domain the propagation time is minutes and the transformation is dominated by refraction, shoaling, sheltering and breaking.

Quickstart

import numpy as np
from waveray import SiteModel, fetch_datamesh_bathymetry

bathy = fetch_datamesh_bathymetry(
    "gebco_2025", bbox=(114.35, -28.95, 114.65, -28.60), positive="up"
)
model = SiteModel.build(
    bathy=bathy,
    target=(114.58, -28.77),                     # lon, lat
    boundary_points=[(114.35, -28.90), (114.35, -28.65)],  # SWAN output sites
    freqs=efth_boundary.freq.values,
    dirs=efth_boundary.dir.values,
)
efth_near = model.transform(efth_boundary, tide=tide_series)   # (time, freq, dir)
model.to_netcdf("geraldton_berth.nc")            # reuse without rebuilding

efth follows the wavespectra convention: dims (time, [site,] freq, dir), dir = coming-from nautical degrees, density per Hz per degree. The site dimension (K boundary points, same order as boundary_points) is required when K > 1.

Ray paths can be exported for inspection (QGIS, EIDOS, any web map) as a GeoJSON FeatureCollection — one MultiLineString per (freq, dir) bin, with period, direction, per-ray status and friction decay in the properties:

from waveray import ray_paths_geojson

x, y = bathy_grid.to_local(lon, lat)
ray_paths_geojson(bathy_grid, (x[0], y[0]), freqs=[0.06, 0.1],
                  dirs=np.arange(0, 360, 15), path="rays.geojson")

Notebooks

Executed, illustrated notebooks (plots included in the committed output):

uv sync --extra datamesh --extra notebooks
DATAMESH_TOKEN=... uv run jupyter lab notebooks/

Limitations (v0)

• No diffraction: accuracy degrades inside harbours / behind breakwaters, and island shadows are sharper than reality (rays block, they do not leak).
• Island / headland sheltering is included — a ray grounding on any land carries no energy — but blocking is binary and only as good as the bathymetry: a feature smaller than ~2 grid cells is smoothed away by the bilinear depth sampling and will not shelter (GEBCO at ~450 m cannot shelter behind a small reef or islet). See tests/test_island.py.
• Bottom friction IS included (JONSWAP, integrated along ray paths) and is ON by default with the SWAN swell coefficient cf_jonswap=0.038; pass cf_jonswap=None for pure refraction + shoaling. No triad interactions.
• Wind input (SWAN formulations) is available but stationary: the wind is baked into the operator at build time, like the bathymetry — build operators per wind condition if the wind matters and varies. Only the wind input term is represented; its nonlinear counterweights (whitecapping, quadruplets) cannot live in a linear operator, so wind input is unbounded and ships with a max_growth ceiling that warns when it binds. This is not a ray-method artefact — SWAN with its sinks disabled grows a 2 m swell to 77.8 m over the same 15 km fetch. An optional wind-sea saturation closure (off by default, empirical rather than derived) caps the energy the wind adds — never the total, so swell is exempt by construction — bringing the measured error against full-physics SWAN to within 13 % with swell present and ~28 % worst case generating a sea across 8–18 m/s, against an uncapped worst of 466 %.
• Breaking is an endpoint cap, not accumulated dissipation along the approach — appropriate at berths and outside the inner surf zone; tune gamma per site against observations.
• The operator is built at a fixed water level; tide only modulates the breaking cap. For strongly tide-dependent refraction, build operators per tide stage.
• Boundary spectra quality bounds achievable skill — nested-model literature consistently finds offshore directional accuracy, not nearshore resolution, is the limiting factor.

Development

uv sync                     # or: uv sync --extra datamesh
uv run pytest               # physics-validation suite
uv run ruff check . && uv run ruff format --check .

Comparison against a real SWAN run is opt-in — it pulls the official SWAN docker image and takes a few minutes:

docker pull delftwaves/swan:latest
uv run pytest -m swan -s    # stationary SWAN 41.51A vs waveray, same boundary spectra