Gap Filling

Filling missing values across the mesh, as opposed to along time.

fill_na_spatial takes each gap from the nearest valid element; laplace_interpolate solves Laplace's equation on the mesh adjacency graph, which reproduces a field linear in x/y exactly. Where UgridDataset.interpolate_na fills along the time axis, these fill across space.

Both are pure with respect to the mesh they are given: the adjacency is derived through graph.sparse_adjacency rather than built onto the mesh, which matters because UgridDataset.copy() shares its Mesh2d — a mutating fill on a copy would reach back into the original.

pyramids.netcdf.ugrid.fill

Mesh-space gap filling for UGRID unstructured meshes.

Pure functions that fill NaN gaps in a per-element data array defined on a :class:~pyramids.netcdf.ugrid.mesh.Mesh2d. Two strategies are provided:

• :func:fill_na_spatial — nearest-valid fill by geometric (centroid/node/edge) distance, using a :class:scipy.spatial.cKDTree over the valid elements.
• :func:laplace_interpolate — smooth fill by solving Laplace's equation on the mesh adjacency graph, with the valid elements as Dirichlet boundary conditions.

Naming distinction (mesh space vs the time axis): these functions fill across mesh space using connectivity/geometry. Filling along the time axis is a separate operation (interpolate_na(dim="time") on the NetCDF side); it is not what these functions do.

The element axis is always the last axis (-1); any leading time/layer axes are broadcast over and filled per slab (each slab uses its own null mask). The input array is never mutated — a filled copy is returned.

fill_na_spatial(mesh, values, location='face', *, max_distance=None)

Fill NaN gaps from the nearest valid element by geometric distance.

Each NaN element takes the value of the nearest non-NaN element of the same location, where "nearest" is Euclidean distance between element coordinates (face centroids, node coordinates, or edge centers). A :class:scipy.spatial.cKDTree is built over the valid elements of each slab.

This fills across mesh space; it is distinct from filling along the time axis (interpolate_na(dim="time")).

Parameters:

Name	Type	Description	Default
mesh	Mesh2d	Source mesh topology.	required
values	NDArray	Data array whose last axis is the element axis; NaN marks a gap. Any leading axes (time/layers) are broadcast over and filled per slab (each slab uses its own null mask).	required
location	str	Mesh location of values: "face", "node" or "edge".	'face'
max_distance	float | None	Maximum distance for a gap to be filled. None (default) is unbounded; when set, gaps with no valid element within this distance stay NaN.	None

Returns:

Type	Description
NDArray	A filled copy of values (float64); the input is not mutated.

Raises:

Type	Description
ValueError	If location is unknown, edge coordinates are missing, or the element axis length does not match the mesh.

Source code in src/pyramids/netcdf/ugrid/fill.py

def fill_na_spatial(
    mesh: Mesh2d,
    values: np.typing.NDArray,
    location: str = "face",
    *,
    max_distance: float | None = None,
) -> np.typing.NDArray:
    """Fill NaN gaps from the nearest valid element by geometric distance.

    Each NaN element takes the value of the nearest non-NaN element of the same
    location, where "nearest" is Euclidean distance between element coordinates
    (face centroids, node coordinates, or edge centers). A :class:`scipy.spatial.cKDTree`
    is built over the valid elements of each slab.

    This fills across mesh **space**; it is distinct from filling along the time
    axis (``interpolate_na(dim="time")``).

    Args:
        mesh: Source mesh topology.
        values: Data array whose **last** axis is the element axis; NaN marks a
            gap. Any leading axes (time/layers) are broadcast over and filled per
            slab (each slab uses its own null mask).
        location: Mesh location of ``values``: ``"face"``, ``"node"`` or ``"edge"``.
        max_distance: Maximum distance for a gap to be filled. ``None`` (default)
            is unbounded; when set, gaps with no valid element within this distance
            stay NaN.

    Returns:
        A filled copy of ``values`` (``float64``); the input is not mutated.

    Raises:
        ValueError: If ``location`` is unknown, edge coordinates are missing, or the
            element axis length does not match the mesh.
    """
    coords = _location_coords(mesh, location)
    data = np.asarray(values, dtype=np.float64)
    if data.shape[-1] != coords.shape[0]:
        raise ValueError(
            f"values last axis ({data.shape[-1]}) does not match the number of "
            f"{location} elements ({coords.shape[0]})."
        )
    out = data.copy()
    flat = out.reshape(-1, out.shape[-1])
    for slab in flat:
        _fill_slab_nearest(coords, slab, max_distance)
    return out

laplace_interpolate(mesh, values, location='face', *, use_weights=True, direct_solve=False, rtol=1e-05, atol=0.0, maxiter=500)

Fill NaN gaps by solving Laplace's equation on the mesh adjacency graph.

The valid elements act as Dirichlet boundary conditions and the NaN elements are the unknowns of a discrete Laplace problem (∇² = 0) on the face-adjacency graph (face variables) or node-adjacency graph (node variables). A sparse graph Laplacian L = diag(degree) - W is built with :mod:scipy.sparse, split by connected component, and the unknown block is solved with :func:scipy.sparse.linalg.spsolve (direct) or :func:scipy.sparse.linalg.cg.

This fills across mesh space; it is distinct from filling along the time axis (interpolate_na(dim="time")).

Parameters:

Name	Type	Description	Default
mesh	Mesh2d	Source mesh topology.	required
values	NDArray	Data array whose last axis is the element axis; NaN marks a gap. Any leading axes (time/layers) are broadcast over and filled per slab (each slab uses its own null mask).	required
location	str	Mesh location of values: "face" (face adjacency) or "node" (node adjacency). "edge" is not supported.	'face'
use_weights	bool	If True (default) use inverse-distance edge weights normalized to a mean of ~1.0; if False use uniform weights.	True
direct_solve	bool	If True solve with spsolve; otherwise use CG.	False
rtol	float	Relative tolerance for the CG solver.	1e-05
atol	float	Absolute tolerance for the CG solver.	0.0
maxiter	int	Maximum CG iterations.	500

Returns:

Type	Description
NDArray	A filled copy of values (float64); the input is not mutated.
NDArray	Elements in connected components that are entirely NaN in a slab stay NaN.

Raises:

Type	Description
ValueError	If location is "edge" or unknown, the element axis length does not match the mesh, or every value is NaN.

Source code in src/pyramids/netcdf/ugrid/fill.py

def laplace_interpolate(
    mesh: Mesh2d,
    values: np.typing.NDArray,
    location: str = "face",
    *,
    use_weights: bool = True,
    direct_solve: bool = False,
    rtol: float = 1e-5,
    atol: float = 0.0,
    maxiter: int = 500,
) -> np.typing.NDArray:
    """Fill NaN gaps by solving Laplace's equation on the mesh adjacency graph.

    The valid elements act as Dirichlet boundary conditions and the NaN elements are
    the unknowns of a discrete Laplace problem (``∇² = 0``) on the face-adjacency
    graph (face variables) or node-adjacency graph (node variables). A sparse graph
    Laplacian ``L = diag(degree) - W`` is built with :mod:`scipy.sparse`, split by
    connected component, and the unknown block is solved with
    :func:`scipy.sparse.linalg.spsolve` (direct) or :func:`scipy.sparse.linalg.cg`.

    This fills across mesh **space**; it is distinct from filling along the time
    axis (``interpolate_na(dim="time")``).

    Args:
        mesh: Source mesh topology.
        values: Data array whose **last** axis is the element axis; NaN marks a
            gap. Any leading axes (time/layers) are broadcast over and filled per
            slab (each slab uses its own null mask).
        location: Mesh location of ``values``: ``"face"`` (face adjacency) or
            ``"node"`` (node adjacency). ``"edge"`` is not supported.
        use_weights: If ``True`` (default) use inverse-distance edge weights
            normalized to a mean of ~1.0; if ``False`` use uniform weights.
        direct_solve: If ``True`` solve with ``spsolve``; otherwise use CG.
        rtol: Relative tolerance for the CG solver.
        atol: Absolute tolerance for the CG solver.
        maxiter: Maximum CG iterations.

    Returns:
        A filled copy of ``values`` (``float64``); the input is not mutated.
        Elements in connected components that are entirely NaN in a slab stay NaN.

    Raises:
        ValueError: If ``location`` is ``"edge"`` or unknown, the element axis length
            does not match the mesh, or every value is NaN.
    """
    if location == "edge":
        raise ValueError(
            "Laplace interpolation is not supported for edge-located variables."
        )
    if location == "face":
        adjacency = _face_adjacency(mesh)
    elif location == "node":
        adjacency = _node_adjacency(mesh)
    else:
        raise ValueError(f"Unknown location: {location!r}")

    coords = _location_coords(mesh, location)
    data = np.asarray(values, dtype=np.float64)
    n_elements = adjacency.shape[0]
    if data.shape[-1] != n_elements:
        raise ValueError(
            f"values last axis ({data.shape[-1]}) does not match the number of "
            f"{location} elements ({n_elements})."
        )
    if np.isnan(data).all():
        raise ValueError("Cannot interpolate: every value is NaN.")

    laplacian = _laplacian(adjacency, coords, use_weights)
    n_components, labels = connected_components(adjacency, directed=False)

    out = data.copy()
    flat = out.reshape(-1, out.shape[-1])
    for slab in flat:
        _solve_slab(
            laplacian, labels, n_components, slab, direct_solve, rtol, atol, maxiter
        )
    return out