interpolate_in_space's depth interpolation treats masked/land corners as depth=0 instead of excluding them

[cgrudz]

Summary

RestartFile.interpolate_in_space's inner _get_depth function (src/roguewave/wavewatch3/restart_file.py) initializes its output buffer with numpy.zeros(...) and never marks masked/land corners as invalid:

def _get_depth(indices, _dummy):
    index = self.grid.index(
        latitude_index=indices[0], longitude_index=indices[1]
    )
    output = numpy.zeros(len(indices[0]))
    mask = index >= 0
    output[mask] = self.depth[index[mask]]
    return output

This is inconsistent with the sibling _get_data in the same method, which explicitly sets output[~mask, :, :] = numpy.nan for masked corners — and inconsistent with RestartFileTimeStack's structurally identical _get_depth_data (restart_file_time_stack.py), which does correctly set output[~mask] = numpy.nan. Same codebase, same pattern, correct in one place and not the other.

Impact

NdInterpolator._data_interpolator only excludes a bilinear corner from the weighted average when its value is NaN. Because _get_depth never produces NaN, every masked land corner is silently treated as a valid depth reading of 0 m and included in the average — biasing interpolated depth toward zero at every coastal point.

Confirmed with a reproduction against the real NdInterpolator class: a synthetic stencil with one masked/land corner and three valid sea corners all at 100 m depth returns 75 m today; the corrected pattern (matching RestartFileTimeStack) returns 100 m.

Investigated downstream impact — conclusion: latent, not currently corrupting any written restart file, given how the pipeline is wired today.

• interpolate_in_space's returned Spectrum.depth feeds real physics if consumed directly (roguewavespectrum's wavenumber/groupspeed/relative_depth/shallow-water calcs; roguewave's own wavephysics/balance/* source-term functions all take spectrum.depth.values).
• But the one real production consumer, esm-python-scripts's preprocess/rewrite_restart_file.py, discards .depth entirely — it extracts only the spectral array before calling write_restart_file, which recomputes the energy↔action Jacobian from the target grid's own bathymetry, never touching the interpolated depth. So today, the biased value is computed and discarded, not written into a restart file.
• Note: while tracing this, found that rewrite_restart_file.py's interpolation branch is currently broken for an unrelated reason (calls .spectral_values, an attribute that doesn't exist on the roguewavespectrum.Spectrum class interpolate_in_space has returned since 2024-03-26 — the script wasn't updated for that switch). Filing that separately; noting here only because it means this depth bug can't be observed end-to-end via that script until it's fixed.

Why this was never caught

No existing unit test asserts interpolated depth values near a masked/land boundary — test_interpolate_in_space in tests/restart_files/test_restart_file.py only checks spectra-derived hm0() at two open-ocean points.

Suggested fix

Add output[~mask] = numpy.nan to _get_depth, mirroring _get_data and RestartFileTimeStack._get_depth_data. Add a test exercising a point whose stencil includes at least one masked corner, asserting the interpolated depth isn't biased toward zero.

Related

Found while scoping the NaN-aware nearest-neighbor fallback for restart-file spatial interpolation (see Notion: Roguewave — Safe Restart-File Interpolation). Same function family (interpolate_in_space), same root-cause class (masked corners not excluded from the bilinear average) as the spectra NaN issue that plan addresses, but a distinct code path and not itself part of that plan's scope.

[pieterbartsmit]

lgtm - one liner.