Interpolation

Interpolation utilities for structured and unstructured grids.

class nos_utils.interp.StructuredGridInterpolator(lon_2d: numpy.ndarray, lat_2d: numpy.ndarray, mask: numpy.ndarray | None = None)[source]

Bases: object

Bilinear interpolation on a 2D curvilinear structured grid.

The grid has shape (ny, nx) with 2D lon/lat arrays. Grid cells are quadrilaterals defined by 4 corner nodes.

interpolate(target_lon: numpy.ndarray, target_lat: numpy.ndarray, data: numpy.ndarray) → numpy.ndarray[source]

Interpolate data field to target points using bilinear interpolation.

Parameters:

target_lon – Target longitudes (-180/180), shape (n_targets,)
target_lat – Target latitudes, shape (n_targets,)
data – Field values on the grid, shape (ny, nx). NaN/fill for land.

Returns:

Interpolated values at target points, shape (n_targets,)

Overview

The nos_utils.interp subpackage provides grid-to-grid interpolation utilities used by the forcing processors. Two flavors are available:

Structured bilinear (nos_utils.interp.structured_interp) — on-the-fly bilinear weights built from a source lon/lat mesh. Useful for ad-hoc interpolation from a regular RTOFS/GFS grid to scattered target points.
Precomputed REMESH weights (nos_utils.interp.precomputed_weights) — weight replay from one-time Fortran INTERP_REMESH exports. This is the hot path used every cycle for SSH, 3D T/S, and T/S nudging on SECOFS and STOFS-3D-ATL. By reusing the weights that the original Fortran prep computed, the Python path is numerically equivalent to the Fortran baseline and roughly 20x faster than rebuilding a Delaunay triangulation every cycle.

Precomputed weight workflow

There are two steps:

One-time build (offline, after any grid change): run the Fortran prep once with NOS_EXPORT_WEIGHTS=YES to dump *_remesh_export.txt text files, then convert them to NPZ using build_ssh_npz(), build_nudge_npz(), or build_3d_npz(). The NPZ is committed to the FIX directory and keyed to the RTOFS grid via an MD5 checksum.
Per-cycle apply: each operational run loads the NPZ with numpy.load(), calls apply_precomputed_ssh() or apply_precomputed_nudge() on the current RTOFS field, and writes the result to elev2D.th.nc, TEM_3D.th.nc, SAL_3D.th.nc, TEM_nu.nc, or SAL_nu.nc.

Example

import numpy as np
from pathlib import Path
from nos_utils.interp.precomputed_weights import (
    build_ssh_npz,
    apply_precomputed_ssh,
)

# --- One-time build (offline) ---
build_ssh_npz(
    export_txt=Path("obc_ssh_remesh_export.txt"),
    rtofs_lon_2d=rtofs_lon,  # (ny, nx) from rtofs_glo_2ds
    rtofs_lat_2d=rtofs_lat,
    out_npz=Path("secofs.obc_ssh_weights.npz"),
)

# --- Per-cycle apply ---
npz = dict(np.load("secofs.obc_ssh_weights.npz"))
ssh_at_obc = apply_precomputed_ssh(npz, ssh_2d=current_rtofs_ssh)
# ssh_at_obc has shape (n_boundary_nodes,)

Structured bilinear interpolator

Structured grid bilinear interpolation for curvilinear grids (e.g., RTOFS).

Replaces Delaunay-based interpolation with cell-search + bilinear on the native grid structure. This eliminates triangulation algorithm differences and convex-hull edge effects that cause ~3cm SSH bias.

Algorithm: 1. Build KDTree from grid cell centers for fast lookup 2. For each target point, find the nearest cell center 3. Check if the target is inside that cell (or its neighbors) 4. Compute bilinear weights within the enclosing quadrilateral 5. Fall back to nearest-neighbor for points outside the grid

class nos_utils.interp.structured_interp.StructuredGridInterpolator(lon_2d: numpy.ndarray, lat_2d: numpy.ndarray, mask: numpy.ndarray | None = None)[source]

Bases: object

Bilinear interpolation on a 2D curvilinear structured grid.

The grid has shape (ny, nx) with 2D lon/lat arrays. Grid cells are quadrilaterals defined by 4 corner nodes.

interpolate(target_lon: numpy.ndarray, target_lat: numpy.ndarray, data: numpy.ndarray) → numpy.ndarray[source]

Interpolate data field to target points using bilinear interpolation.

Parameters:

target_lon – Target longitudes (-180/180), shape (n_targets,)
target_lat – Target latitudes, shape (n_targets,)
data – Field values on the grid, shape (ny, nx). NaN/fill for land.

Returns:

Interpolated values at target points, shape (n_targets,)

Precomputed REMESH weights

Precomputed INTERP_REMESH weight replay for SSH, 3D T/S, and nudging interpolation.

Uses weights exported from the Fortran INTERP_REMESH subroutine to produce numerically equivalent interpolation without Delaunay triangulation. The weights are computed once and stored in .npz files in the FIX directory.

SSH weights (obc_ssh_weights.npz):: 1488 boundary nodes for elev2D.th.nc Source grid: global RTOFS 2D (3298x4500)
3D weights (obc_3d_weights.npz):: 1488 boundary nodes for TEM_3D.th.nc / SAL_3D.th.nc Source grid: regional RTOFS 3D (e.g., US_east 1710x742)
Nudge weights (obc_nudge_weights.npz):: ~32K interior nodes for TEM_nu.nc / SAL_nu.nc Source grid: global RTOFS 2D (3298x4500)

Runtime: 1. First call: KDTree maps source coordinates to RTOFS grid cells (cached) 2. Per timestep: gather data at source cells, apply weights (pure indexed ops)

nos_utils.interp.precomputed_weights.load_remesh_export(filepath: Path) → dict[source]

Load the one-time Fortran INTERP_REMESH export text file.

Returns dict with source coordinates, target mapping, weights, modes, donors.

nos_utils.interp.precomputed_weights.build_npz(export_txt: Path, rtofs_lon_2d: numpy.ndarray, rtofs_lat_2d: numpy.ndarray, out_npz: Path, tol: float = 0.01) → None[source]

One-time conversion: export text → production .npz with grid indices.

Parameters:

export_txt – Path to obc_ssh_remesh_export.txt
rtofs_lon_2d – RTOFS 2D longitude array (ny, nx)
rtofs_lat_2d – RTOFS 2D latitude array (ny, nx)
out_npz – Output NPZ path
tol – Max allowed distance (degrees) for source-to-grid matching

nos_utils.interp.precomputed_weights.validate_grid(npz: dict, rtofs_lon_2d: numpy.ndarray, rtofs_lat_2d: numpy.ndarray) → None[source]: Verify the RTOFS grid matches the stored weight mapping.

nos_utils.interp.precomputed_weights.apply_precomputed_ssh(npz: dict, ssh_2d: numpy.ndarray) → numpy.ndarray[source]

Apply precomputed Fortran REMESH weights to RTOFS SSH field.

Numerically equivalent to Fortran, with residual differences on the order of floating-point roundoff (REAL*4 weights, float64 arithmetic).

Parameters:

npz – Dict from np.load(‘secofs.obc_ssh_weights.npz’)
ssh_2d – Current cycle RTOFS SSH field, shape matching grid_shape

Returns:

SSH at boundary nodes, shape (n_target,). Does NOT include SSH offset.

nos_utils.interp.precomputed_weights.build_nudge_npz(export_txt: Path, rtofs_lon_2d: numpy.ndarray, rtofs_lat_2d: numpy.ndarray, out_npz: Path, tol: float = 0.01) → None[source]

One-time conversion: nudge export text -> production .npz with grid indices.

The export text has the same format as the SSH export (produced by INTERP_REMESH with NUDGE_REMESH_ACTIVE), but for the ~32K interior nudging nodes instead of ~1488 boundary nodes.

Parameters:

export_txt – Path to obc_nudge_remesh_export.txt
rtofs_lon_2d – RTOFS 2D longitude array (ny, nx)
rtofs_lat_2d – RTOFS 2D latitude array (ny, nx)
out_npz – Output NPZ path
tol – Max allowed distance (degrees) for source-to-grid matching

nos_utils.interp.precomputed_weights.build_3d_npz(export_txt: Path, rtofs_lon_2d: numpy.ndarray, rtofs_lat_2d: numpy.ndarray, out_npz: Path, tol: float = 0.01) → None[source]

One-time conversion: 3D T/S export text -> production .npz with grid indices.

The export text has the same format as the SSH and nudge exports (produced by INTERP_REMESH with TS3D_REMESH_ACTIVE), but uses the regional RTOFS 3D grid (e.g., US_east 1710x742) instead of the global 2D grid (3298x4500). Target nodes are the same 1488 boundary nodes as SSH.

Parameters:

export_txt – Path to obc_3d_remesh_export.txt
rtofs_lon_2d – RTOFS 3D regional longitude array (ny, nx)
rtofs_lat_2d – RTOFS 3D regional latitude array (ny, nx)
out_npz – Output NPZ path
tol – Max allowed distance (degrees) for source-to-grid matching

nos_utils.interp.precomputed_weights.apply_precomputed_nudge(npz: dict, field_2d: numpy.ndarray, fill_value: float = numpy.nan) → numpy.ndarray[source]

Apply precomputed Fortran REMESH weights to a 2D RTOFS field for nudging.

Works identically to apply_precomputed_ssh but for the ~32K interior nudging target nodes. Each call interpolates ONE 2D slice (one depth level, one timestep) from the RTOFS grid to the nudge target nodes.

Parameters:

npz – Dict from np.load(‘secofs.obc_nudge_weights.npz’)
field_2d – Single 2D RTOFS field, shape matching grid_shape. Land/fill values should be set to NaN or abs >= 99.
fill_value – Value to use for corner-mean when all data is NaN. Defaults to NaN (let caller handle fill).

Returns:

Interpolated values at nudge nodes, shape (n_target,), float32.