Radius-1 IDW neighbor fallback and domain-gap fill for interpolate_in_space by cgrudz · Pull Request #16 · sofarocean/roguewave

cgrudz · 2026-08-12T00:12:35Z

Description

What kind of change is this?
Bug Fix / Addition

Scope

src/roguewave/interpolate/nd_interp.py — NdInterpolator
src/roguewave/wavewatch3/restart_file.py — RestartFile.interpolate_in_space, RestartFile._fancy_index, RestartFile.fill_missing_spectra (new)
src/roguewave/wavewatch3/io.py — write_restart_file
tests/interpolate/test_nd_interp.py, tests/restart_files/test_fill_missing_spectra.py (new), tests/restart_files/test_empty_index.py (new), tests/restart_files/test_io.py, tests/interpolate/__init__.py (new, test package scaffolding)

Current behavior:
Sofar's 0.25°/0.5°/1.0° WW3 grids share an origin and are exact integer-ratio refinements of one another, so a target point's bilinear stencil always collapses to 100% weight on a single coincident source grid point. When that coincident point is masked/land, interpolate_in_space returns a NaN spectrum for the target point, even when the point is a single source-grid step away from open ocean on every side. Measured against a real production 0.25° restart, this affects up to several percent of target sea points depending on target resolution, concentrated at real coastlines (Antarctic coast, Indonesian archipelago, etc.), invisible to skipna=True bulk diagnostics but present in the raw array written to restart files.

New behavior:

NdInterpolator gains an opt-in nan_fallback_radius constructor parameter (default 0, preserving current behavior for its other three callers: dataset.py, dataarray.py, cluster.py). When a point's primary bilinear lookup fails and nan_fallback_radius > 0, the interpolator checks that many source-grid-index steps around the coincident point — including the coincident point itself, which can still be the best available candidate for a general (non-grid-aligned) miss — inverse-distance-weights (haversine) whichever of those have valid data, and still returns NaN if none do. So a genuine domain gap (e.g. a source grid with no sea points nearby at all) is left alone rather than silently papered over. RestartFile.interpolate_in_space passes nan_fallback_radius=1, matching the fix's mechanism to the actual grid ratios in play. The sibling depth interpolator and RestartFileTimeStack are intentionally untouched — separate, already-tracked concerns.
For the residual points where even the radius-1 fallback correctly leaves NaN — genuine domain gaps like Hudson Bay, absent from the 0.25° source grid entirely, not a coastal artifact — RestartFile.fill_missing_spectra fills them with a WW3 cold-start-style spectrum instead, matching ww3_strt's ITYPE options: "calm" (zero energy, default, matches ITYPE 5), "user_defined" (broadcast a caller-supplied spectrum, matches ITYPE 4), "gaussian" (matches ITYPE 1), and "jonswap" (matches ITYPE 2). These are roguewave's own implementations of the same named spectral shapes/conventions, not a port of WW3's Fortran. Points that already have data are left untouched.
Fixed a separate, pre-existing latent bug surfaced by (1): RestartFile._fancy_index crashed on an empty index array — numpy.array([]) on an empty list of read spectra collapses to shape (0,) instead of (0, number_of_frequencies, number_of_directions), which xarray then rejected as a dimension-size conflict during Spectrum/Dataset construction. This is common whenever a radius-1 neighborhood is entirely land (e.g. exactly the Hudson Bay case above), so the fallback in (1) would not have worked on real data without it.
Fixed a separate latent bug in write_restart_file: it read .variance_density instead of .directional_variance_density on Spectrum/Dataset input, which broke composing interpolate_in_space's output directly into write_restart_file — exactly the composition this PR's regridding path relies on. Fixes write_restart_file reads .variance_density instead of .directional_variance_density -- breaks when given a Spectrum/Dataset with directional data #15. Added a regression test (test_clone_restart_file) covering the one in-repo path that hits this branch, previously untested.

Verified against the real production 0.25°→0.5° restart (2026-08-05 00Z cycle): before this PR, 885/156635 (0.565%) target sea points came back NaN; after the radius-1 fallback, 469/156635 (0.2994%) remain NaN, all genuine domain gaps (447 in Hudson Bay); after fill_missing_spectra(fill_type="calm"), 0 remain NaN. Separately, on the parallel cluster, a restart file produced by this full path (interpolate_in_space + fill_missing_spectra + write_restart_file) was confirmed to load correctly and let ww3_multi (the real production-pinned esm-box:0.5.7 build) step forward cleanly with no error, at production's exact 320-rank/10-node scale.

Links to context docs:

interpolate_in_space's depth interpolation treats masked/land corners as depth=0 instead of excluding them #14 (separate depth-interpolation bug, out of scope here)
write_restart_file reads .variance_density instead of .directional_variance_density -- breaks when given a Spectrum/Dataset with directional data #15 (write_restart_file bug, fixed by this PR)

Deploy plan

None beyond a normal merge — pure library change, no infra side effects. Downstream consumers pick this up on their next roguewave version pin bump.

TODO

Update CHANGELOG.md — N/A, this repo has no CHANGELOG.md

…erpolate_in_space Sofar's 0.25/0.5/1.0 deg WW3 grids share an origin, so bilinear interpolation always collapses to a single coincident source point; when that point is masked/land, interpolate_in_space returned NaN instead of using nearby valid data. NdInterpolator now takes an opt-in nan_fallback_radius, defaulting off for its other callers, that inverse-distance-weights valid neighbors within that radius and still returns NaN when none exist (e.g. real domain gaps).

Copilot

Pull request overview

Adds radius-1 IDW fallback for masked WW3 grid points while preserving default interpolation behavior.

Changes:

Adds configurable Haversine-based neighbor fallback.
Enables fallback for restart-file spectrum interpolation.
Adds unit tests for fallback behavior.

Reviewed changes

Copilot reviewed 3 out of 4 changed files in this pull request and generated 2 comments.

File	Description
`src/roguewave/interpolate/nd_interp.py`	Implements neighbor fallback logic.
`src/roguewave/wavewatch3/restart_file.py`	Enables radius-1 fallback for spectra.
`tests/interpolate/test_nd_interp.py`	Tests fallback scenarios.
`tests/interpolate/__init__.py`	Marks the test package.

Suppressed comments (3)

src/roguewave/interpolate/nd_interp.py:376

This hard-codes the interpolated point dimension as axis 0. NdInterpolator supports passive dimensions before interp_index_coord_name, in which case neighbor_value and the output place the point axis elsewhere; constructing and boolean-indexing value_by_point this way produces incompatible shapes. Use the existing output-index helpers for both source and destination (which also avoids allocating another full output-sized array per neighbor).

            value_by_point = numpy.zeros(
                (number_points,) + neighbor_value.shape[1:], dtype=numpy.float64
            )
            value_by_point[usable_global_position] = neighbor_value[
                neighbor_value_is_valid

src/roguewave/interpolate/nd_interp.py:290

These are coordinates of the nearest source node, not of the requested interpolation point. For any non-coincident coastal target whose bilinear lookup fails, the fallback therefore computes IDW around the wrong location (for example, east/west neighbors become equally weighted even when the target is closer to one). Thread the requested latitude/longitude values into this helper and use them for the Haversine distances, or explicitly limit fallback eligibility to coincident targets.

        target_latitude = latitude_values[
            coincident_source_index_of_failed_points[latitude_axis_index]
        ]
        target_longitude = longitude_values[
            coincident_source_index_of_failed_points[longitude_axis_index]

src/roguewave/interpolate/nd_interp.py:347

neighbor_value_is_valid is sized only for the in-bounds subset, but it is applied here to arrays sized for all failed points. In a batch containing both an edge failure and an interior failure, an offset that excludes only the edge point causes a boolean-dimension mismatch and aborts interpolation. Index these coordinates with the already mapped usable_local_position.

            neighbor_latitude_value = latitude_values[
                neighbor_source_index_per_axis[latitude_axis_index][
                    neighbor_value_is_valid
                ]

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

_fancy_index crashed on an empty index array (numpy.array([]) collapses to shape (0,) instead of (0, num_freq, num_dir)), hit whenever a radius-1 neighborhood is entirely land -- common for domain gaps and the case that most needs to degrade to NaN gracefully. Also, per automated PR review: exclude out-of-domain points from the fallback (their clipped bracket indices aren't a real coincident node), measure fallback distances from the actually-requested point rather than the coincident source node, fix a sizing mismatch that could crash on a batch mixing in- and out-of-bounds points at the same offset, and fix a periodic-longitude test that never exercised the modulo wraparound it claimed to test.

cgrudz · 2026-08-12T01:30:51Z

Addressing the three additional findings from the review body's "suppressed comments" (these weren't rendered as separate inline threads, so replying here instead) — all fixed in 6bdd585:

nd_interp.py:376 (hard-coded point axis at 0): confirmed real — inherited from the same undocumented assumption already present in the primary bilinear loop above it (val[self.output_indexing_full(mask)] makes the identical assumption). Rather than silently extend that assumption further, documented it explicitly in a comment on _radius_neighbor_fallback per team discussion, since fixing it generically would mean auditing the whole class's positional-alignment convention across all 5 existing callers — out of scope for this opt-in feature, and not exercised by any current caller (output_index_coord_index == 0 holds for both interpolate_in_space call sites).
nd_interp.py:290 (target coordinates from the coincident source node instead of the requested point): confirmed real design gap, harmless only because of this coordinate pair's specific grid-nesting coincidence (verified: zero partial-weight cases exist for these grid ratios). Fixed by threading the actual requested points["latitude"]/points["longitude"] through into the distance calculation instead of re-deriving from the coincident index.
nd_interp.py:347 (neighbor_value_is_valid sized to the in-bounds subset but applied to full-size arrays): confirmed real bug, would raise on any batch mixing an in-bounds and out-of-bounds failed point at the same offset (e.g. near a pole). Fixed by subsetting neighbor_source_index_per_axis to the in-bounds candidates once, before any of the size-dependent lookups.

Re-ran the full unit test suite and the real production 0.25°→0.5° validation (469/156635 NaN, unchanged) after these fixes to confirm no regression.

Copilot

Pull request overview

Copilot reviewed 3 out of 4 changed files in this pull request and generated 1 comment.

Suppressed comments (1)

src/roguewave/interpolate/nd_interp.py:395

This hardcodes the interpolation-point axis as axis 0, even though NdInterpolator supports passive dimensions before that axis via output_index_coord_index. For data shaped like (frequency, latitude, longitude), neighbor_value is (frequency, failed_points), so this allocation has the wrong shape and the boolean index on the next line is applied to the frequency axis, causing an indexing error. Build the buffer in output_shape and use the existing axis-aware index helper.

            value_by_point = numpy.zeros(
                (number_points,) + neighbor_value.shape[1:], dtype=numpy.float64
            )
            value_by_point[usable_global_position] = neighbor_value[

_radius_neighbor_fallback previously allocated its accumulators (and a fresh scratch tensor per neighbor offset) at the full output shape -- number_points * frequency * direction -- even though fallback work only ever touches the failed-point subset, typically a small fraction of a production grid. For the shape in this PR (156635 points), that's gigabytes of unnecessary allocation. Size everything to the failed subset instead and scatter the small resolved result into the primary result at the end.

The radius-1 fallback correctly leaves a point NaN when no radius-1 neighbor has data at all -- a real domain gap (e.g. Hudson Bay, absent from the 0.25 source grid entirely) rather than a coastal artifact, per team discussion on the NASA OSTST Phase 2 dev plan. For those residual points, add a WW3 cold-start-style fill matching ww3_strt's ITYPE options: "calm" (zero energy, default), "user_defined" (broadcast a supplied spectrum), "gaussian", and "jonswap" (parametric shapes normalized/parameterized to match those ITYPEs' conventions, not a port of WW3's Fortran). Confirmed against the real production 0.25->0.5 restart: interpolate_in_space's 469 residual NaN points all resolve to 0 after fill_missing_spectra(fill_type="calm").

Copilot

Pull request overview

Copilot reviewed 4 out of 5 changed files in this pull request and generated 1 comment.

Suppressed comments (2)

src/roguewave/interpolate/nd_interp.py:323

The radius search excludes offset zero, which drops the nearest source node for non-grid-aligned queries. A primary lookup can fail with weights_sum <= 0.5 even when that nearest node is valid (for example, it is the only valid bilinear corner); the fallback then ignores the best available value and may return NaN or blend only farther nodes. Include the center candidate—the exact-coincident failure case remains unchanged because its center value is invalid.

        neighbor_offsets = [
            offset
            for offset in itertools.product(
                range(-radius, radius + 1), repeat=self.interp_ndims
            )
            if any(offset)
        ]

src/roguewave/wavewatch3/restart_file.py:353

The new empty-index path is not exercised by the added tests: their fake get_data bypasses RestartFile.__getitem__ and _fancy_index, although the real fallback reaches this branch whenever a neighbor batch contains no sea indices. Add a regression test using a real or minimally constructed RestartFile to verify that empty fancy indexing survives the conversion and Spectrum construction with shape (0, frequency, direction).

        if len(indices) == 0:
            # numpy.array([]) on an empty list of spectra collapses to shape
            # (0,) instead of (0, number_of_frequencies, number_of_directions),
            # which xarray then rejects as a dimension-size conflict.
            return numpy.empty(

Excluding the zero offset assumed the coincident node was always already-known-invalid whenever the fallback triggers. That holds for today's grid-nested use case (its weight is always exactly 0 or 1, never partial), but not in general: a non-grid-aligned bilinear miss can have a valid coincident node whose own weight is still <= 0.5 (e.g. the largest of four roughly-even corner weights), which the fallback would previously discard in favor of a worse, farther neighbor. Also adds a regression test for the empty fancy-index crash fixed earlier, using a fully synthetic, in-memory RestartFile (Grid and MetaData are plain dataclasses, no file I/O needed) -- the existing fallback tests use a synthetic get_data callback that never actually exercised RestartFile._fancy_index's real empty-array path.

cgrudz · 2026-08-12T20:56:33Z

Addressing the other finding from the latest review's "suppressed comments" (`restart_file.py:353`, no separate inline thread rendered for it) — fixed in 6ba6897:

Confirmed real gap: the fallback tests use a synthetic `get_data` callback that never calls the real `RestartFile.getitem`/`_fancy_index`, so the empty-index crash fix (from an earlier round) had zero regression coverage — only a manual rerun of the production validation script would have caught it breaking again.

Added `tests/restart_files/test_empty_index.py` with a fully synthetic, in-memory `RestartFile` — `Grid` and `MetaData` are plain dataclasses with no file I/O, and the empty-index path never touches `resource`, so this needed no S3 access or binary parsing. It exercises the real `getitem` → `_fancy_index` → `to_frequency_energy_density` → `Spectrum`/`Dataset` construction path end to end with an empty index array, asserting the correct (0, frequency, direction) shape survives all the way through.

Copilot

Pull request overview

Copilot reviewed 5 out of 6 changed files in this pull request and generated no new comments.

Suppressed comments (2)

src/roguewave/interpolate/nd_interp.py:398

This boolean index always targets axis 0, although neighbor_value_is_valid indexes the interpolated-point axis. If a passive dimension precedes the interpolation coordinates—for example data shaped (frequency, latitude, longitude, direction)—neighbor_value is (frequency, candidates, direction), so this raises a boolean-length mismatch (or selects frequencies). Move the configured point axis to the accumulator's axis 0 before filtering.

            usable_value = neighbor_value[neighbor_value_is_valid]

src/roguewave/interpolate/nd_interp.py:228

resolved_value is built with failed points on axis 0, but this destination keeps the interpolated-point axis at output_index_coord_index. For a supported layout such as (frequency, latitude, longitude, direction), that axis is 1, so assigning (failed, frequency, direction) into (frequency, failed, direction) fails or misorders values. Move the internal point axis back to the output position before scattering.

        result[self.output_indexing_full(failed_point_position)] = resolved_value

Its docstring described only the return shape, with no mention that masked/land coincident points are now filled from radius-1 neighbors when possible, or that fill_missing_spectra exists for the residual domain-gap points that still come back NaN.

…nal_variance_density Breaks composing interpolate_in_space's output directly with write_restart_file, which this branch's regridding path relies on. Fixes #15. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

cgrudz requested a balanced review from Copilot August 12, 2026 00:13

Copilot started reviewing on behalf of cgrudz August 12, 2026 00:14 View session

Copilot AI reviewed Aug 12, 2026

View reviewed changes

Comment thread src/roguewave/interpolate/nd_interp.py Outdated

Comment thread tests/interpolate/test_nd_interp.py Outdated

cgrudz requested a balanced review from Copilot August 12, 2026 01:32

Copilot started reviewing on behalf of cgrudz August 12, 2026 01:32 View session

Copilot AI reviewed Aug 12, 2026

View reviewed changes

Comment thread src/roguewave/interpolate/nd_interp.py Outdated

cgrudz added 2 commits August 12, 2026 09:58

cgrudz changed the title ~~Radius-1 IDW neighbor fallback for masked/land points in interpolate_in_space~~ Radius-1 IDW neighbor fallback and domain-gap fill for interpolate_in_space Aug 12, 2026

cgrudz requested a balanced review from Copilot August 12, 2026 19:19

Copilot started reviewing on behalf of cgrudz August 12, 2026 19:20 View session

Copilot AI reviewed Aug 12, 2026

View reviewed changes

Comment thread src/roguewave/interpolate/nd_interp.py

cgrudz requested a balanced review from Copilot August 12, 2026 21:01

Copilot started reviewing on behalf of cgrudz August 12, 2026 21:02 View session

Copilot AI reviewed Aug 12, 2026

View reviewed changes

cgrudz and others added 2 commits August 12, 2026 14:22

Navigation Menu

Uh oh!

Radius-1 IDW neighbor fallback and domain-gap fill for interpolate_in_space - #16

Radius-1 IDW neighbor fallback and domain-gap fill for interpolate_in_space#16
cgrudz wants to merge 7 commits into
mainfrom
cg/nan_filter_restart_interpolation

cgrudz commented Aug 12, 2026 •

edited

Loading

Uh oh!

Copilot AI left a comment

Uh oh!

Uh oh!

Uh oh!

cgrudz commented Aug 12, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Uh oh!

Copilot AI left a comment

Uh oh!

Uh oh!

cgrudz commented Aug 12, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

cgrudz commented Aug 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Deploy plan

TODO

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Uh oh!

cgrudz commented Aug 12, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Uh oh!

cgrudz commented Aug 12, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

cgrudz commented Aug 12, 2026 •

edited

Loading