Interpolation

Value resampling between grids and point clouds: regrid, coarsen / refine, gap-fill, bin scattered observations onto a grid, smooth, and sample a grid at arbitrary points (e.g. along a satellite track). Grid, SpaceTimeGrid, and Period are the lightweight target-grid carriers.

Regridding & resolution

Three ways to move a field between grids, by what they preserve:

• RegridLike — xr.interp onto another grid's coordinates. Smooth point values; does not preserve integrals.
• Coarsen(conservative=True) — integer-factor block averaging with cosine-latitude weights.
• RegridConservative — first-order conservative regrid onto an arbitrary rectilinear grid (mode="mean" for intensive fields such as concentrations, mode="sum" for extensive per-cell totals such as emissions). Cell bounds are inferred from the centres (midpoints, end cells symmetric) unless CF bounds are present; on the default spherical geometry longitude is periodic with period 360°, so a 0–360 source regrids onto a −180–180 target directly. For an aligned integer factor it reproduces Coarsen(conservative=True).

import numpy as np
from xrops.interpolate import RegridConservative

# 0.1° inventory -> analysis grid; basin totals are preserved.
regrid = RegridConservative(
    {"lat": np.arange(40.05, 45.0, 0.1), "lon": np.arange(-5.95, 0.0, 0.1)},
    mode="sum",
)
prior = regrid(inventory)

RegridLike

RegridLike(target: Dataset | DataArray, *, dims: tuple[str, ...] = ('lat', 'lon'), method: str = 'linear')

Bases: Operator

Wrap :func:xrops.interpolate.regrid_like — bilinear resample of the input onto another Dataset's coordinate grid along dims.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    target: xr.Dataset | xr.DataArray,
    *,
    dims: tuple[str, ...] = ("lat", "lon"),
    method: str = "linear",
):
    self.target = target
    self.dims = tuple(dims)
    self.method = method

RegridConservative

RegridConservative(target: GridLike, *, dims: tuple[str, str] = ('lat', 'lon'), geometry: Geometry = 'spherical', target_bounds: Mapping[str, BoundsLike] | None = None, mode: RegridMode = 'mean', normalize: Normalize = 'fracarea', skipna: bool = True)

Bases: Operator

Wrap :func:xrops.interpolate.regrid_conservative — first-order conservative (area-weighted) regrid onto another rectilinear grid.

mode="mean" preserves the area integral of intensive fields; mode="sum" preserves the total of extensive per-cell fields. See the primitive for the bounds inference, longitude wrap and NaN rules.

The effective target cell bounds are resolved at construction (explicit target_bounds > CF bounds attribute > midpoints between centres) and serialised by :meth:get_config as {dim: [[lo, hi], ...]}, so an operator rebuilt from its config regrids onto exactly the same cells.

Parameters:

Name	Type	Description	Default
target	GridLike	Target grid — a dataset / data array whose coordinates include dims (a CF bounds attribute is honoured), or a mapping {dim: centres}.	required
dims	tuple[str, str]	The two dims to regrid, (lat, lon) order on a spherical geometry.	('lat', 'lon')
geometry	Geometry	"spherical" or "planar".	'spherical'
target_bounds	Mapping[str, BoundsLike] | None	Explicit target cell bounds per dim, as (n + 1,) edges or CF (n, 2) bounds; overrides the target's own.	None
mode	RegridMode	"mean" or "sum".	'mean'
normalize	Normalize	"fracarea" or "destarea".	'fracarea'
skipna	bool	Whether NaN source cells are skipped.	True

Example

>>> import numpy as np, xarray as xr
>>> from xrops.interpolate import RegridConservative
>>> lat, lon = np.arange(-9.5, 10.0, 1.0), np.arange(0.5, 20.0, 1.0)
>>> ds = xr.Dataset(
...     {"flux": (("lat", "lon"), np.ones((lat.size, lon.size)))},
...     coords={"lat": lat, "lon": lon},
... )
>>> target = {"lat": np.arange(-8.0, 10.0, 4.0), "lon": [5.0, 15.0]}
>>> op = RegridConservative(target)
>>> op(ds)["flux"].shape
(5, 2)

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    target: _grid_to_grid.GridLike,
    *,
    dims: tuple[str, str] = ("lat", "lon"),
    geometry: _grid_to_grid.Geometry = "spherical",
    target_bounds: Mapping[str, _grid_to_grid.BoundsLike] | None = None,
    mode: _grid_to_grid.RegridMode = "mean",
    normalize: _grid_to_grid.Normalize = "fracarea",
    skipna: bool = True,
):
    self.dims = tuple(dims)
    if len(self.dims) != 2:
        raise ValueError(f"dims must name exactly two dimensions, got {dims!r}")
    self.target = target
    self._target_coords = {
        d: _grid_to_grid._centres(target, d, what="target") for d in self.dims
    }
    self._target_bounds = {
        d: _grid_to_grid._resolve_intervals(target, d, target_bounds, what="target")
        for d in self.dims
    }
    self.geometry = geometry
    self.mode = mode
    self.normalize = normalize
    self.skipna = skipna

Coarsen

Coarsen(factor: dict[str, int], method: str = 'mean', boundary: str = 'trim', conservative: bool = False, lat: str = 'lat')

Bases: Operator

Wrap :func:xrops.interpolate.coarsen.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    factor: dict[str, int],
    method: str = "mean",
    boundary: str = "trim",
    conservative: bool = False,
    lat: str = "lat",
):
    if boundary not in self._VALID_BOUNDARY:
        raise ValueError(
            f"Coarsen boundary must be one of {self._VALID_BOUNDARY!r}, "
            f"got {boundary!r}."
        )
    if conservative and method != "mean":
        raise ValueError(
            f"conservative coarsen only supports method='mean', got {method!r}."
        )
    # Reuse the layer-0 validator so int-likes (np.int64) are accepted and
    # negative / zero / non-integer factors fail at construction time.
    self.factor = _grid_to_grid._validate_coarsen_factor(factor)
    self.method = method
    self.boundary = boundary
    self.conservative = conservative
    self.lat = lat

Refine

Refine(factor: dict[str, int | float], method: str = 'linear', *, order: int | None = None, lat: str = 'lat', lon: str = 'lon', anti_aliasing: bool | None = None, mode: ResizeMode = 'reflect', cval: float = 0.0)

Bases: Operator

Wrap :func:xrops.interpolate.refine.

If order is set, dispatches to the scikit-image-backed 2-D resize path and requires factor to include both lat and lon.

Parameters:

Name	Type	Description	Default
factor	dict[str, int | float]	Per-dimension refinement factors.	required
method	str	Interpolation method for the default xr.interp path.	'linear'
order	int | None	Optional scikit-image spline order (0..5). When set, uses the 2-D resize path.	None
lat	str	Latitude-like dimension name for the 2-D resize path.	'lat'
lon	str	Longitude-like dimension name for the 2-D resize path.	'lon'
anti_aliasing	bool | None	Anti-aliasing setting passed to scikit-image.	None
mode	ResizeMode	Boundary mode for scikit-image ("reflect", "constant", "edge", "symmetric", or "wrap").	'reflect'
cval	float	Fill value used when mode="constant".	0.0

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    factor: dict[str, int | float],
    method: str = "linear",
    *,
    order: int | None = None,
    lat: str = "lat",
    lon: str = "lon",
    anti_aliasing: bool | None = None,
    mode: ResizeMode = "reflect",
    cval: float = 0.0,
):
    if order is not None:
        extra = set(factor) - {lat, lon}
        if extra:
            raise ValueError(
                "Refine(order=...) only resizes the (lat, lon) plane; "
                f"got extra factor dims {sorted(extra)!r}. Drop them or "
                "use order=None."
            )
    self.factor = dict(factor)
    self.method = method
    self.order = order
    self.lat = lat
    self.lon = lon
    self.anti_aliasing = anti_aliasing
    self.mode = mode
    self.cval = cval

Upscale module-attribute

Upscale = _downscale.Upscale

Downscale module-attribute

Downscale = _downscale.Downscale

ResampleTime

ResampleTime(freq: str = '1D', method: str = 'mean', time: str = 'time', interp_method: Literal['linear', 'nearest', 'cubic'] = 'linear')

Bases: Operator

Wrap :func:xrops.interpolate.resample_time.

Parameters:

Name	Type	Description	Default
freq	str	Pandas-style frequency string (e.g. "1D", "6H").	'1D'
method	str	Zero-argument aggregation method ("mean", "sum", "max", …) or "interpolate" for upsampling.	'mean'
time	str	Name of the time dimension.	'time'
interp_method	Literal['linear', 'nearest', 'cubic']	Interpolation method forwarded to xarray when method="interpolate".	'linear'

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    freq: str = "1D",
    method: str = "mean",
    time: str = "time",
    interp_method: Literal["linear", "nearest", "cubic"] = "linear",
):
    self.freq = freq
    self.method = method
    self.time = time
    self.interp_method = interp_method

Gap-filling

FillNaNSpatial

FillNaNSpatial(method: str = 'linear', lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Wrap :func:xrops.interpolate.fillnan_spatial.

Parameters:

Name	Type	Description	Default
method	str	"linear", "nearest", or "cubic".	'linear'
lon	str	Name of the longitude coordinate.	'lon'
lat	str	Name of the latitude coordinate.	'lat'

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(self, method: str = "linear", lon: str = "lon", lat: str = "lat"):
    self.method = method
    self.lon = lon
    self.lat = lat

FillNaNTemporal

FillNaNTemporal(method: str = 'linear', time: str = 'time', max_gap: Any = None)

Bases: Operator

Wrap :func:xrops.interpolate.fillnan_temporal.

Parameters:

Name	Type	Description	Default
method	str	Any method accepted by xarray's interpolate_na ("linear", "nearest", "cubic", etc.).	'linear'
time	str	Name of the time dimension.	'time'
max_gap	Any	Maximum time delta to interpolate across; gaps wider than this stay NaN.	None

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    method: str = "linear",
    time: str = "time",
    max_gap: Any = None,
):
    self.method = method
    self.time = time
    self.max_gap = max_gap

FillNaNLaplacian

FillNaNLaplacian(*, max_iter: int = 1000, tol: float = 0.0001, relaxation: float = 1.0, boundary: str = 'reflect', lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Wrap :func:xrops.interpolate.fillnan_laplacian.

Parameters:

Name	Type	Description	Default
max_iter	int	Maximum red/black Gauss-Seidel sweeps.	1000
tol	float	Absolute convergence tolerance for missing-cell updates.	0.0001
relaxation	float	SOR relaxation factor. 1.0 is Gauss-Seidel.	1.0
boundary	str	"reflect" for Neumann edges, or "wrap" to wrap longitude edges while reflecting latitude edges.	'reflect'
lon	str	Name of the longitude dimension.	'lon'
lat	str	Name of the latitude dimension.	'lat'

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    *,
    max_iter: int = 1000,
    tol: float = 1e-4,
    relaxation: float = 1.0,
    boundary: str = "reflect",
    lon: str = "lon",
    lat: str = "lat",
):
    # Validate eagerly so misconfigured operators fail at construction
    # time rather than deep inside _apply (mirrors Coarsen).
    _gap_fill._validate_laplacian_args(max_iter, tol, relaxation, boundary)
    self.max_iter = max_iter
    self.tol = tol
    self.relaxation = relaxation
    self.boundary = boundary
    self.lon = lon
    self.lat = lat

FillNaNBiharmonic

FillNaNBiharmonic(*, lon: str = 'lon', lat: str = 'lat', split_into_regions: bool = True)

Bases: Operator

Wrap :func:xrops.interpolate.fillnan_biharmonic.

Parameters:

Name	Type	Description	Default
lon	str	Name of the longitude dimension.	'lon'
lat	str	Name of the latitude dimension.	'lat'
split_into_regions	bool	Forwarded to scikit-image. When True, each connected masked region is solved independently.	True

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    *,
    lon: str = "lon",
    lat: str = "lat",
    split_into_regions: bool = True,
):
    self.lon = lon
    self.lat = lat
    self.split_into_regions = split_into_regions

FillNaNRBF

FillNaNRBF(kernel: str = 'thin_plate_spline', neighbors: int | None = 32, lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Wrap :func:xrops.interpolate.fillnan_rbf.

Parameters:

Name	Type	Description	Default
kernel	str	RBF kernel passed to :class:scipy.interpolate.RBFInterpolator.	'thin_plate_spline'
neighbors	int | None	Number of nearest neighbors used per query point; None uses all points.	32
lon	str	Name of the longitude coordinate.	'lon'
lat	str	Name of the latitude coordinate.	'lat'

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    kernel: str = "thin_plate_spline",
    neighbors: int | None = 32,
    lon: str = "lon",
    lat: str = "lat",
):
    self.kernel = kernel
    self.neighbors = neighbors
    self.lon = lon
    self.lat = lat

FillNaNIDW

FillNaNIDW(*, lon: str = 'lon', lat: str = 'lat', k: int = 8, power: float = 2.0, metric: Metric = 'euclidean', max_distance: float | None = None, eps: float = 1e-12)

Bases: Operator

Wrap :func:xrops.interpolate.fillnan_idw.

Parameters:

Name	Type	Description	Default
lon	str	Name of the longitude coordinate and dimension.	'lon'
lat	str	Name of the latitude coordinate and dimension.	'lat'
k	int	Number of nearest finite neighbours to use.	8
power	float	Inverse-distance exponent. 0 gives an unweighted kNN mean.	2.0
metric	Metric	"euclidean" for array-index distances or "haversine" for great-circle distances from coordinate values.	'euclidean'
max_distance	float | None	Optional neighbour cutoff (radians for "haversine", index units for "euclidean").	None
eps	float	Small non-negative offset for non-exact distance weights.	1e-12

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    *,
    lon: str = "lon",
    lat: str = "lat",
    k: int = 8,
    power: float = 2.0,
    metric: _knn.Metric = "euclidean",
    max_distance: float | None = None,
    eps: float = 1e-12,
):
    _knn._validate_idw_args(k, power, metric, max_distance, eps)
    self.lon = lon
    self.lat = lat
    self.k = k
    self.power = power
    self.metric = metric
    self.max_distance = max_distance
    self.eps = eps

FillNaNClimatology

FillNaNClimatology(*, time: str = 'time', group: Literal['month', 'dayofyear', 'season'] = 'month', residual: Literal['zero', 'linear'] = 'linear', min_count: int = 1)

Bases: Operator

Wrap :func:xrops.interpolate.fillnan_climatology.

Parameters:

Name	Type	Description	Default
time	str	Name of the time dimension.	'time'
group	Literal['month', 'dayofyear', 'season']	Calendar grouping used to compute the climatology.	'month'
residual	Literal['zero', 'linear']	"zero" fills missing values with climatology only; "linear" linearly interpolates the anomaly and adds it back.	'linear'
min_count	int	Minimum finite observations required for each climatology group.	1

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    *,
    time: str = "time",
    group: Literal["month", "dayofyear", "season"] = "month",
    residual: Literal["zero", "linear"] = "linear",
    min_count: int = 1,
):
    self.time = time
    self.group = group
    self.residual = residual
    self.min_count = min_count

Binning & gridding scattered data

Bin2D

Bin2D(grid: Grid, statistic: str = 'mean', lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Wrap :func:xrops.interpolate.bin_2d.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    grid: _binning.Grid,
    statistic: str = "mean",
    lon: str = "lon",
    lat: str = "lat",
):
    self.grid = grid
    self.statistic = statistic
    self.lon = lon
    self.lat = lat

Histogram2D

Histogram2D(grid: Grid, lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Wrap :func:xrops.interpolate.histogram_2d.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(self, grid: _binning.Grid, lon: str = "lon", lat: str = "lat"):
    self.grid = grid
    self.lon = lon
    self.lat = lat

KDEToGrid

KDEToGrid(grid: Grid, *, bandwidth: float | str = 'scott', kernel: str = 'gaussian', metric: str = 'euclidean', algorithm: str = 'auto', output: str = 'density', rtol: float = 0.0001)

Bases: Operator

Wrap :func:xrops.interpolate.kde_to_grid.

Expects (lons, lats) or (lons, lats, weights) as input.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    grid: _binning.Grid,
    *,
    bandwidth: float | str = "scott",
    kernel: str = "gaussian",
    metric: str = "euclidean",
    algorithm: str = "auto",
    output: str = "density",
    rtol: float = 1e-4,
):
    self.grid = grid
    self.bandwidth = bandwidth
    self.kernel = kernel
    self.metric = metric
    self.algorithm = algorithm
    self.output = output
    self.rtol = float(rtol)

PointsToGrid

PointsToGrid(grid: Grid, statistic: str = 'mean')

Bases: Operator

Wrap :func:xrops.interpolate.points_to_grid.

Expects a 3-tuple (lons, lats, values) as input.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(self, grid: _binning.Grid, statistic: str = "mean"):
    self.grid = grid
    self.statistic = statistic

IDWToGrid

IDWToGrid(grid: Grid, *, k: int = 8, power: float = 2.0, metric: Metric = 'euclidean', max_distance: float | None = None, eps: float = 1e-12)

Bases: Operator

Wrap :func:xrops.interpolate.idw_to_grid.

Expects a 3-tuple (lons, lats, values) as input.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    grid: _binning.Grid,
    *,
    k: int = 8,
    power: float = 2.0,
    metric: _knn.Metric = "euclidean",
    max_distance: float | None = None,
    eps: float = 1e-12,
):
    _knn._validate_idw_args(k, power, metric, max_distance, eps)
    self.grid = grid
    self.k = k
    self.power = power
    self.metric = metric
    self.max_distance = max_distance
    self.eps = eps

IDWToPoints

IDWToPoints(dst_lons: ndarray, dst_lats: ndarray, *, k: int = 8, power: float = 2.0, metric: Metric = 'euclidean', max_distance: float | None = None, eps: float = 1e-12)

Bases: Operator

Wrap :func:xrops.interpolate.idw_to_points.

Expects source (lons, lats, values) as input.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    dst_lons: np.ndarray,
    dst_lats: np.ndarray,
    *,
    k: int = 8,
    power: float = 2.0,
    metric: _knn.Metric = "euclidean",
    max_distance: float | None = None,
    eps: float = 1e-12,
):
    _knn._validate_idw_args(k, power, metric, max_distance, eps)
    self.dst_lons = np.asarray(dst_lons)
    self.dst_lats = np.asarray(dst_lats)
    self.k = k
    self.power = power
    self.metric = metric
    self.max_distance = max_distance
    self.eps = eps

Smoothing & filtering

GaussianSmooth

GaussianSmooth(dim: str, sigma: float, *, truncate: float = 4.0)

Bases: Operator

Wrap :func:xrops.interpolate._src.smooth.gaussian_smooth.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(self, dim: str, sigma: float, *, truncate: float = 4.0):
    if sigma <= 0:
        raise ValueError(f"sigma must be > 0, got {sigma}")
    self.dim = dim
    self.sigma = float(sigma)
    self.truncate = float(truncate)

GaussianSmoothMasked

GaussianSmoothMasked(*, dim: str | Sequence[str], sigma: float | Mapping[str, float], truncate: float = 4.0, mode: str = 'reflect', nan_aware: bool = True, min_weight: float = 1e-06)

Bases: Operator

Wrap :func:xrops.interpolate._src.smooth.gaussian_smooth_masked.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    *,
    dim: str | Sequence[str],
    sigma: float | Mapping[str, float],
    truncate: float = 4.0,
    mode: str = "reflect",
    nan_aware: bool = True,
    min_weight: float = 1e-6,
):
    dims = _smooth._normalize_dims(dim)
    sigmas = _smooth._normalize_sigmas(sigma, dims)
    if truncate <= 0:
        raise ValueError(f"truncate must be > 0, got {truncate}")
    if min_weight < 0:
        raise ValueError(f"min_weight must be >= 0, got {min_weight}")

    self.dim: str | tuple[str, ...] = dim if isinstance(dim, str) else dims
    self.sigma: float | dict[str, float]
    if isinstance(sigma, Mapping):
        self.sigma = {d: sigmas[i] for i, d in enumerate(dims)}
    else:
        self.sigma = sigmas[0]
    self.truncate = float(truncate)
    self.mode = mode
    self.nan_aware = bool(nan_aware)
    self.min_weight = float(min_weight)

MovingAverage

MovingAverage(dim: str, window: int, *, center: bool = True, min_periods: int | None = None)

Bases: Operator

Wrap :func:xrops.interpolate._src.smooth.moving_average.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    dim: str,
    window: int,
    *,
    center: bool = True,
    min_periods: int | None = None,
):
    if not isinstance(window, int) or isinstance(window, bool):
        raise TypeError(f"window must be an int, got {type(window).__name__}")
    if window < 1:
        raise ValueError(f"window must be >= 1, got {window}")
    if min_periods is not None and (
        not isinstance(min_periods, int) or min_periods < 0
    ):
        raise ValueError(
            f"min_periods must be a non-negative int or None, got {min_periods!r}"
        )
    self.dim = dim
    self.window = window
    self.center = bool(center)
    self.min_periods = min_periods

LowpassFilter

LowpassFilter(dim: str, cutoff: Any, *, order: int = 4, btype: str = 'low')

Bases: Operator

Wrap :func:xrops.interpolate._src.smooth.lowpass_filter.

For btype in {"low", "high", "lowpass", "highpass"} cutoff is a scalar in (0, 1). For btype in {"bandpass", "bandstop"} it is a length-2 sequence. Validation is delegated to the private kernel.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    dim: str,
    cutoff: Any,
    *,
    order: int = 4,
    btype: str = "low",
):
    if not isinstance(order, int) or isinstance(order, bool):
        raise TypeError(f"order must be an int, got {type(order).__name__}")
    self.dim = dim
    if np.isscalar(cutoff):
        self.cutoff: Any = float(cutoff)
    else:
        pair = tuple(float(v) for v in cutoff)
        if len(pair) != 2:
            raise ValueError(f"cutoff sequence must have length 2, got {len(pair)}")
        self.cutoff = pair
    self.order = order
    self.btype = btype

Mask cleanup

CleanMask

CleanMask(*, fill_holes_area: int | None = 4, drop_objects_area: int | None = None, closing_footprint: Footprint | None = None, opening_footprint: Footprint | None = None, lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Convenience operator for the common mask-cleanup pipeline.

The fixed order is remove small holes, remove small objects, binary closing, then binary opening; each step is opt-in through its keyword.

Parameters:

Name	Type	Description	Default
fill_holes_area	int | None	Optional hole area threshold.	4
drop_objects_area	int | None	Optional object area threshold.	None
closing_footprint	Footprint | None	Optional binary-closing footprint.	None
opening_footprint	Footprint | None	Optional binary-opening footprint.	None
lon	str	Longitude dimension name.	'lon'
lat	str	Latitude dimension name.	'lat'

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    *,
    fill_holes_area: int | None = 4,
    drop_objects_area: int | None = None,
    closing_footprint: _mask_ops.Footprint | None = None,
    opening_footprint: _mask_ops.Footprint | None = None,
    lon: str = "lon",
    lat: str = "lat",
):
    self.fill_holes_area = fill_holes_area
    self.drop_objects_area = drop_objects_area
    self.closing_footprint = _decode_footprint_config(closing_footprint)
    self.opening_footprint = _decode_footprint_config(opening_footprint)
    self.lon = lon
    self.lat = lat

MaskBinaryOpening

MaskBinaryOpening(*, footprint: Footprint = 1, lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Operator wrapper for binary opening inside Sequential pipelines.

Parameters:

Name	Type	Description	Default
footprint	Footprint	Structuring element specification.	1
lon	str	Longitude dimension name.	'lon'
lat	str	Latitude dimension name.	'lat'

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    *,
    footprint: _mask_ops.Footprint = 1,
    lon: str = "lon",
    lat: str = "lat",
):
    self.footprint = _decode_footprint_config(footprint)
    self.lon = lon
    self.lat = lat

MaskBinaryClosing

MaskBinaryClosing(*, footprint: Footprint = 1, lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Operator wrapper for binary closing inside Sequential pipelines.

Parameters:

Name	Type	Description	Default
footprint	Footprint	Structuring element specification.	1
lon	str	Longitude dimension name.	'lon'
lat	str	Latitude dimension name.	'lat'

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    *,
    footprint: _mask_ops.Footprint = 1,
    lon: str = "lon",
    lat: str = "lat",
):
    self.footprint = _decode_footprint_config(footprint)
    self.lon = lon
    self.lat = lat

MaskRemoveSmallHoles

MaskRemoveSmallHoles(*, area: int = 4, lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Operator wrapper for filling small unmasked holes in mask pipelines.

Parameters:

Name	Type	Description	Default
area	int	False regions smaller than this pixel count are filled.	4
lon	str	Longitude dimension name.	'lon'
lat	str	Latitude dimension name.	'lat'

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(self, *, area: int = 4, lon: str = "lon", lat: str = "lat"):
    self.area = area
    self.lon = lon
    self.lat = lat

MaskRemoveSmallObjects

MaskRemoveSmallObjects(*, area: int = 4, lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Operator wrapper for dropping small masked specks in mask pipelines.

Parameters:

Name	Type	Description	Default
area	int	True regions smaller than this pixel count are dropped.	4
lon	str	Longitude dimension name.	'lon'
lat	str	Latitude dimension name.	'lat'

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(self, *, area: int = 4, lon: str = "lon", lat: str = "lat"):
    self.area = area
    self.lon = lon
    self.lat = lat

Point sampling

SampleAtPoints

SampleAtPoints(points: Dataset, *, coords: tuple[str, ...] = ('lat', 'lon'), method: Method = 'linear', fill_value: float | None = nan, bounds_error: bool = False, point_dim: str = 'points')

Bases: Operator

Wrap :func:xrops.interpolate.sample_at_points.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    points: xr.Dataset,
    *,
    coords: tuple[str, ...] = ("lat", "lon"),
    method: _grid_to_points.Method = "linear",
    fill_value: float | None = np.nan,
    bounds_error: bool = False,
    point_dim: str = "points",
):
    self.points = points
    self.coords = tuple(coords)
    self.method = method
    self.fill_value = fill_value
    self.bounds_error = bounds_error
    self.point_dim = point_dim

AlongTrack

AlongTrack(track: Dataset, *, coords: tuple[str, ...] = ('time', 'lat', 'lon'), method: Method = 'linear', fill_value: float | None = nan, point_dim: str = 'points')

Bases: Operator

Wrap :func:xrops.interpolate.along_track.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    track: xr.Dataset,
    *,
    coords: tuple[str, ...] = ("time", "lat", "lon"),
    method: _grid_to_points.Method = "linear",
    fill_value: float | None = np.nan,
    point_dim: str = "points",
):
    self.track = track
    self.coords = tuple(coords)
    self.method = method
    self.fill_value = fill_value
    self.point_dim = point_dim

Vertical & axis remapping

RemapAxis replaces one dimension of every numeric variable by a new axis of target values; the vertical presets pin conventional source / target names. Two options cover terrain-following (WRF η) and hybrid (ERA5 model-level) grids, whose source levels vary per column:

• source_coords= names a coordinate on the input — e.g. z_agl(time, level, y, x) — that supplies the source axis values for every column. Each column must be strictly monotone along the source dim (ascending or descending, and the orientation may differ between columns).
• extrapolate="nearest" holds each column's end value beyond its own range instead of the default "nan". An analysis level that sits below the lowest model level of a column thus gets the lowest-level value.

import numpy as np
from xrops.interpolate import ToHeight

# ds: (time, level, y, x) model-level fields plus a 4-D ``z_agl`` coordinate.
to_height = ToHeight(
    np.array([10.0, 50.0, 100.0, 250.0]),
    source_coords="z_agl",
    extrapolate="nearest",
)
out = to_height(ds)  # ``level`` -> ``height``; variables without ``level`` pass through

Chunked inputs stay lazy along every dim except the source dim, which must sit in a single chunk. With a 1-D dimension coordinate and the default extrapolate="nan" the operator is unchanged from earlier releases.

RemapAxis

RemapAxis(source_axis: str, target_axis: DataArray | ndarray | list, *, target_name: str | None = None, method: str = 'linear', source_coords: str | None = None, extrapolate: Extrapolate = 'nan')

Bases: Operator

Generic axis remapping (D12).

Replaces the source_axis dimension in the input Dataset with a new dimension whose coordinate values are target_axis. Every numeric variable that carries source_axis is interpolated onto the target axis.

Parameters:

Name	Type	Description	Default
source_axis	str	Name of the existing dimension to remap.	required
target_axis	DataArray | ndarray | list	Target coordinate values. If an :class:xr.DataArray, its .name becomes the new dim name; otherwise the new dim name defaults to target_name or source_axis.	required
target_name	str | None	Optional explicit new dim name.	None
method	str	"linear" or "nearest".	'linear'
source_coords	str | None	Name of a coordinate on the input supplying per-column source axis values (e.g. "z_agl" on a terrain-following grid). It must carry source_axis and may vary along the other dims. None uses the 1-D dimension coordinate.	None
extrapolate	Extrapolate	"nan" (default) — targets outside a column's range become NaN; "nearest" — they hold the column's end value.	'nan'

Example

Model-level temperature with a 4-D height-above-ground coordinate onto fixed heights, holding the end values beyond each column:

>>> import numpy as np, xarray as xr
>>> from xrops.interpolate import ToHeight
>>> level = np.arange(4)[None, :, None, None]
>>> z_agl = 20.0 * (level + 1) + np.zeros((2, 4, 3, 3))
>>> ds = xr.Dataset(
...     {"t": (("time", "level", "y", "x"), 300.0 - 0.01 * z_agl)},
...     coords={"z_agl": (("time", "level", "y", "x"), z_agl)},
... )
>>> heights = np.array([10.0, 50.0, 90.0])
>>> op = ToHeight(heights, source_coords="z_agl", extrapolate="nearest")
>>> op(ds)["t"].dims
('time', 'height', 'y', 'x')

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    source_axis: str,
    target_axis: xr.DataArray | np.ndarray | list,
    *,
    target_name: str | None = None,
    method: str = "linear",
    source_coords: str | None = None,
    extrapolate: Extrapolate = "nan",
):
    self.source_axis = source_axis
    if isinstance(target_axis, xr.DataArray):
        self._target_da = target_axis
        self._target_values: np.ndarray = np.asarray(
            target_axis.values, dtype=float
        )
        self._inferred_name = target_axis.name
    else:
        self._target_da = None
        self._target_values = np.asarray(target_axis, dtype=float)
        self._inferred_name = None
    self.target_name = target_name
    self.method = method
    self.source_coords = source_coords
    self.extrapolate = extrapolate

ToHeight

ToHeight(target_axis: DataArray | ndarray | list, *, source_axis: str | None = None, target_name: str | None = None, method: str = 'linear', source_coords: str | None = None, extrapolate: Extrapolate = 'nan')

Bases: _VerticalPreset

Remap a pressure or hybrid axis to geometric height.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    target_axis: xr.DataArray | np.ndarray | list,
    *,
    source_axis: str | None = None,
    target_name: str | None = None,
    method: str = "linear",
    source_coords: str | None = None,
    extrapolate: Extrapolate = "nan",
):
    super().__init__(
        source_axis=source_axis or self._DEFAULT_SOURCE,
        target_axis=target_axis,
        target_name=target_name or self._DEFAULT_TARGET_NAME,
        method=method,
        source_coords=source_coords,
        extrapolate=extrapolate,
    )

ToPressureLevels

ToPressureLevels(target_axis: DataArray | ndarray | list, *, source_axis: str | None = None, target_name: str | None = None, method: str = 'linear', source_coords: str | None = None, extrapolate: Extrapolate = 'nan')

Bases: _VerticalPreset

Remap a height/depth axis to standard pressure levels.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    target_axis: xr.DataArray | np.ndarray | list,
    *,
    source_axis: str | None = None,
    target_name: str | None = None,
    method: str = "linear",
    source_coords: str | None = None,
    extrapolate: Extrapolate = "nan",
):
    super().__init__(
        source_axis=source_axis or self._DEFAULT_SOURCE,
        target_axis=target_axis,
        target_name=target_name or self._DEFAULT_TARGET_NAME,
        method=method,
        source_coords=source_coords,
        extrapolate=extrapolate,
    )

ToSigma

ToSigma(target_axis: DataArray | ndarray | list, *, source_axis: str | None = None, target_name: str | None = None, method: str = 'linear', source_coords: str | None = None, extrapolate: Extrapolate = 'nan')

Bases: _VerticalPreset

Remap a depth axis to terrain-following sigma values.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    target_axis: xr.DataArray | np.ndarray | list,
    *,
    source_axis: str | None = None,
    target_name: str | None = None,
    method: str = "linear",
    source_coords: str | None = None,
    extrapolate: Extrapolate = "nan",
):
    super().__init__(
        source_axis=source_axis or self._DEFAULT_SOURCE,
        target_axis=target_axis,
        target_name=target_name or self._DEFAULT_TARGET_NAME,
        method=method,
        source_coords=source_coords,
        extrapolate=extrapolate,
    )

FromSigma

FromSigma(target_axis: DataArray | ndarray | list, *, source_axis: str | None = None, target_name: str | None = None, method: str = 'linear', source_coords: str | None = None, extrapolate: Extrapolate = 'nan')

Bases: _VerticalPreset

Remap a sigma axis back to a fixed depth grid.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    target_axis: xr.DataArray | np.ndarray | list,
    *,
    source_axis: str | None = None,
    target_name: str | None = None,
    method: str = "linear",
    source_coords: str | None = None,
    extrapolate: Extrapolate = "nan",
):
    super().__init__(
        source_axis=source_axis or self._DEFAULT_SOURCE,
        target_axis=target_axis,
        target_name=target_name or self._DEFAULT_TARGET_NAME,
        method=method,
        source_coords=source_coords,
        extrapolate=extrapolate,
    )

ToIsopycnal

ToIsopycnal(target_axis: DataArray | ndarray | list, *, source_axis: str | None = None, target_name: str | None = None, method: str = 'linear', source_coords: str | None = None, extrapolate: Extrapolate = 'nan')

Bases: _VerticalPreset

Remap a depth axis to potential-density (isopycnal) levels.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    target_axis: xr.DataArray | np.ndarray | list,
    *,
    source_axis: str | None = None,
    target_name: str | None = None,
    method: str = "linear",
    source_coords: str | None = None,
    extrapolate: Extrapolate = "nan",
):
    super().__init__(
        source_axis=source_axis or self._DEFAULT_SOURCE,
        target_axis=target_axis,
        target_name=target_name or self._DEFAULT_TARGET_NAME,
        method=method,
        source_coords=source_coords,
        extrapolate=extrapolate,
    )

ToPhase

ToPhase(time_dim: str, period: float, n_bins: int, *, epoch: float = 0.0)

Bases: Operator

Fold a time axis onto a phase axis by binning + averaging.

Phase is computed as ((t - epoch) / period) mod 1 and binned into n_bins evenly-spaced bins on [0, 1).

Parameters

time_dim Name of the time dimension. period Length of one cycle, in the same units as the time coordinate. n_bins Number of phase bins. epoch Reference time at which phase = 0.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/operators.py

def __init__(
    self,
    time_dim: str,
    period: float,
    n_bins: int,
    *,
    epoch: float = 0.0,
):
    if period <= 0:
        raise ValueError(f"period must be > 0, got {period}")
    if n_bins < 1:
        raise ValueError(f"n_bins must be >= 1, got {n_bins}")
    self.time_dim = time_dim
    self.period = float(period)
    self.n_bins = int(n_bins)
    self.epoch = float(epoch)

Functional primitives (Layer 0)

These pure functions back the operators above; each takes xr.DataArray/xr.Dataset and a dim: argument.

regrid_like

regrid_like(ds: Dataset | DataArray, target: Dataset | DataArray, *, dims: Iterable[str] = ('lat', 'lon'), method: str = 'linear') -> Dataset | DataArray

Resample ds onto target's coordinates along dims.

Thin :meth:xr.Dataset.interp wrapper for the common "regrid model output to observation grid" step. Coordinates listed in dims must exist on both ds and target; values along other dims pass through.

Assumes ds and target share a coordinate space — if they are in different CRSs, use :func:xrops.geom.reproject_match instead, which is CRS-aware and matches target's pixel grid exactly. Interpolation does not preserve integrals; for fluxes, inventories and other fields whose area integral must survive the regrid use :func:regrid_conservative.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/grid_to_grid.py

def regrid_like(
    ds: xr.Dataset | xr.DataArray,
    target: xr.Dataset | xr.DataArray,
    *,
    dims: Iterable[str] = ("lat", "lon"),
    method: str = "linear",
) -> xr.Dataset | xr.DataArray:
    """Resample ``ds`` onto ``target``'s coordinates along ``dims``.

    Thin :meth:`xr.Dataset.interp` wrapper for the common
    "regrid model output to observation grid" step. Coordinates listed
    in ``dims`` must exist on both ``ds`` and ``target``; values along
    other dims pass through.

    Assumes ``ds`` and ``target`` share a coordinate space — if they are
    in different CRSs, use :func:`xrops.geom.reproject_match` instead,
    which is CRS-aware and matches ``target``'s pixel grid exactly.
    Interpolation does not preserve integrals; for fluxes, inventories
    and other fields whose area integral must survive the regrid use
    :func:`regrid_conservative`.
    """
    dim_list = list(dims)
    missing_target = [d for d in dim_list if d not in target.coords]
    if missing_target:
        raise ValueError(
            f"target is missing requested dims {missing_target!r} as coords; "
            f"got coords {tuple(target.coords)}. Pass `dims=` explicitly to "
            "regrid only the dims that are actually shared."
        )
    missing_source = [d for d in dim_list if d not in ds.coords]
    if missing_source:
        raise ValueError(
            f"input is missing requested dims {missing_source!r} as coords; "
            f"got coords {tuple(ds.coords)}."
        )
    target_coords = {d: target[d] for d in dim_list}
    return ds.interp(target_coords, method=method)

regrid_conservative

regrid_conservative(ds: Dataset | DataArray, target: GridLike, *, dims: tuple[str, str] = ('lat', 'lon'), geometry: Geometry = 'spherical', bounds: Mapping[str, BoundsLike] | None = None, target_bounds: Mapping[str, BoundsLike] | None = None, mode: RegridMode = 'mean', normalize: Normalize = 'fracarea', skipna: bool = True) -> Dataset | DataArray

First-order conservative regrid along two rectilinear dims.

Source cells \(i\) (area \(a_i\), value \(f_i\)) and target cells \(j\) overlap on areas \(w_{ij} = |c_i \cap c_j|\). Because both grids are rectilinear the overlap factorises into one 1-D weight matrix per axis (:func:overlap_weights_1d), so the whole regrid is two matrix products per field:

\[\text{mean: } \bar f_j = \frac{\sum_i w_{ij} f_i}{\sum_i w_{ij}\,[f_i \text{ valid}]},\qquad \text{sum: } F_j = \sum_i \frac{w_{ij}}{a_i} f_i\]

mode="mean" preserves intensive fields (concentration, mixing ratio, flux density). With normalize="destarea" the area integral \(\sum_j A_j \bar f_j = \sum_i a_i f_i\) is conserved whenever the target covers the source; with normalize="fracarea" (default) it is conserved only over target cells that are fully covered by unmasked source cells — a partially covered or partially masked target cell holds the mean of its valid part, which over-counts that cell's contribution to the integral. mode="sum" preserves extensive per-cell fields (mass, emission per cell): each source cell's value is split among the target cells in proportion to the overlapped fraction, so \(\sum_j F_j = \sum_i f_i\).

On a "spherical" geometry dims[0] is latitude and dims[1] longitude, both in degrees, and cell areas are \(R^2 \Delta\lambda (\sin\varphi_+ - \sin\varphi_-)\) (\(R\) cancels). Longitude is periodic with period 360 degrees: the target bounds are unwrapped into the source's range, so a source on 0–360 regrids onto a target on −180–180 (and vice versa) without pre-rolling. On a "planar" geometry areas are \(\Delta x\,\Delta y\) in coordinate units and nothing is periodic.

Cell bounds default to the midpoints between centres, with the end cells extrapolated symmetrically. A CF bounds attribute on a coordinate (naming an (n, 2) variable) is honoured, and bounds / target_bounds override both, as (n + 1,) edges or CF (n, 2) bounds per dim. CF bounds need not be contiguous: cells may leave gaps (or overlap), and only the overlap with each cell's own (lo, hi) interval counts.

Missing values: with skipna=True (default) NaN source cells carry no weight. normalize="fracarea" divides by the valid overlap, so a partially masked (or partially covered) target cell keeps the mean of its valid part; "destarea" divides by the full target cell area (ESMF / xESMF naming), scaling the mean down by the valid fraction. A target cell with no valid overlap is NaN in every mode. With skipna=False any NaN source cell that touches a target cell makes it NaN. Complex fields are regridded component-wise (a cell is valid when both parts are finite) and stay complex.

Extra dims (time, level) broadcast; dask inputs stay lazy along them, while the two regridded dims must each sit in a single chunk. Datasets are regridded per variable: variables without either dim pass through, CF bounds variables of the source grid are dropped (those of the target grid, when it is a dataset / data array carrying them, are attached), and a variable carrying only one of the two dims raises. A data array output has no room for a bounds variable, so its coordinates carry no dangling bounds attribute.

Parameters:

Name	Type	Description	Default
ds	Dataset | DataArray	Source dataset or data array with 1-D, strictly monotone coordinates for both dims.	required
target	GridLike	Target grid — a dataset / data array whose coordinates include dims, or a mapping {dim: centres}.	required
dims	tuple[str, str]	The two dims to regrid, (lat, lon) order on a spherical geometry.	('lat', 'lon')
geometry	Geometry	"spherical" (degrees, sin-latitude weighting, periodic longitude) or "planar".	'spherical'
bounds	Mapping[str, BoundsLike] | None	Explicit source cell bounds per dim.	None
target_bounds	Mapping[str, BoundsLike] | None	Explicit target cell bounds per dim.	None
mode	RegridMode	"mean" (intensive) or "sum" (extensive).	'mean'
normalize	Normalize	"fracarea" or "destarea"; only affects mode="mean".	'fracarea'
skipna	bool	Whether NaN source cells are skipped (default) or propagate.	True

Returns:

Type	Description
Dataset | DataArray	ds on the target grid, with the target centres as coordinates
Dataset | DataArray	along dims and every other dim unchanged. Values are
Dataset | DataArray	float64 (complex128 for complex input).

Raises:

Type	Description
ValueError	If an option is unknown, a dim or its coordinate is missing, a coordinate is not strictly monotone, bounds have the wrong shape, longitude bounds span more than one period, a dask input is chunked along a regridded dim, or a Dataset variable carries only one of the two dims.

Example

Move a 1-degree per-cell emission inventory (0–360 longitudes) onto a 2-degree grid on −180–180; the global total is preserved to round-off:

>>> import numpy as np, xarray as xr
>>> from xrops.interpolate import regrid_conservative
>>> lat, lon = np.arange(-89.5, 90.0, 1.0), np.arange(0.5, 360.0, 1.0)
>>> rng = np.random.default_rng(0)
>>> emis = xr.DataArray(
...     rng.uniform(0.0, 1.0, size=(lat.size, lon.size)),
...     dims=("lat", "lon"),
...     coords={"lat": lat, "lon": lon},
... )
>>> target = {
...     "lat": np.arange(-89.0, 90.0, 2.0),
...     "lon": np.arange(-179.0, 180.0, 2.0),
... }
>>> out = regrid_conservative(emis, target, mode="sum")
>>> out.sizes
Frozen({'lat': 90, 'lon': 180})
>>> bool(np.isclose(float(out.sum()), float(emis.sum()), rtol=1e-12))
True

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/grid_to_grid.py

def regrid_conservative(
    ds: xr.Dataset | xr.DataArray,
    target: GridLike,
    *,
    dims: tuple[str, str] = ("lat", "lon"),
    geometry: Geometry = "spherical",
    bounds: Mapping[str, BoundsLike] | None = None,
    target_bounds: Mapping[str, BoundsLike] | None = None,
    mode: RegridMode = "mean",
    normalize: Normalize = "fracarea",
    skipna: bool = True,
) -> xr.Dataset | xr.DataArray:
    r"""First-order conservative regrid along two rectilinear dims.

    Source cells $i$ (area $a_i$, value $f_i$) and target cells $j$ overlap
    on areas $w_{ij} = |c_i \cap c_j|$. Because both grids are rectilinear
    the overlap factorises into one 1-D weight matrix per axis
    (:func:`overlap_weights_1d`), so the whole regrid is two matrix
    products per field:

    $$\text{mean: } \bar f_j = \frac{\sum_i w_{ij} f_i}{\sum_i w_{ij}\,[f_i
    \text{ valid}]},\qquad \text{sum: } F_j = \sum_i \frac{w_{ij}}{a_i} f_i$$

    ``mode="mean"`` preserves *intensive* fields (concentration, mixing
    ratio, flux density). With ``normalize="destarea"`` the area integral
    $\sum_j A_j \bar f_j = \sum_i a_i f_i$ is conserved whenever the target
    covers the source; with ``normalize="fracarea"`` (default) it is
    conserved only over target cells that are fully covered by unmasked
    source cells — a partially covered or partially masked target cell
    holds the mean of its valid part, which over-counts that cell's
    contribution to the integral. ``mode="sum"`` preserves *extensive*
    per-cell fields (mass, emission per cell): each source cell's value
    is split among the target cells in proportion to the overlapped
    fraction, so $\sum_j F_j = \sum_i f_i$.

    On a ``"spherical"`` geometry ``dims[0]`` is latitude and ``dims[1]``
    longitude, both in degrees, and cell areas are $R^2 \Delta\lambda
    (\sin\varphi_+ - \sin\varphi_-)$ ($R$ cancels). Longitude is periodic
    with period 360 degrees: the target bounds are unwrapped into the
    source's range, so a source on 0–360 regrids onto a target on
    −180–180 (and vice versa) without pre-rolling. On a ``"planar"``
    geometry areas are $\Delta x\,\Delta y$ in coordinate units and nothing
    is periodic.

    Cell bounds default to the midpoints between centres, with the end
    cells extrapolated symmetrically. A CF ``bounds`` attribute on a
    coordinate (naming an ``(n, 2)`` variable) is honoured, and ``bounds``
    / ``target_bounds`` override both, as ``(n + 1,)`` edges or CF
    ``(n, 2)`` bounds per dim. CF bounds need not be contiguous: cells
    may leave gaps (or overlap), and only the overlap with each cell's
    own ``(lo, hi)`` interval counts.

    Missing values: with ``skipna=True`` (default) NaN source cells carry
    no weight. ``normalize="fracarea"`` divides by the *valid* overlap, so
    a partially masked (or partially covered) target cell keeps the mean
    of its valid part; ``"destarea"`` divides by the full target cell area
    (ESMF / xESMF naming), scaling the mean down by the valid fraction.
    A target cell with no valid overlap is NaN in every mode. With
    ``skipna=False`` any NaN source cell that touches a target cell makes
    it NaN. Complex fields are regridded component-wise (a cell is valid
    when both parts are finite) and stay complex.

    Extra dims (``time``, ``level``) broadcast; dask inputs stay lazy
    along them, while the two regridded dims must each sit in a single
    chunk. Datasets are regridded per variable: variables without either
    dim pass through, CF bounds variables of the source grid are dropped
    (those of the target grid, when it is a dataset / data array carrying
    them, are attached), and a variable carrying only one of the two dims
    raises. A data array output has no room for a bounds variable, so its
    coordinates carry no dangling ``bounds`` attribute.

    Args:
        ds: Source dataset or data array with 1-D, strictly monotone
            coordinates for both ``dims``.
        target: Target grid — a dataset / data array whose coordinates
            include ``dims``, or a mapping ``{dim: centres}``.
        dims: The two dims to regrid, ``(lat, lon)`` order on a spherical
            geometry.
        geometry: ``"spherical"`` (degrees, sin-latitude weighting,
            periodic longitude) or ``"planar"``.
        bounds: Explicit source cell bounds per dim.
        target_bounds: Explicit target cell bounds per dim.
        mode: ``"mean"`` (intensive) or ``"sum"`` (extensive).
        normalize: ``"fracarea"`` or ``"destarea"``; only affects
            ``mode="mean"``.
        skipna: Whether NaN source cells are skipped (default) or
            propagate.

    Returns:
        ``ds`` on the target grid, with the target centres as coordinates
        along ``dims`` and every other dim unchanged. Values are
        ``float64`` (``complex128`` for complex input).

    Raises:
        ValueError: If an option is unknown, a dim or its coordinate is
            missing, a coordinate is not strictly monotone, bounds have the
            wrong shape, longitude bounds span more than one period, a dask
            input is chunked along a regridded dim, or a Dataset variable
            carries only one of the two dims.

    Example:
        Move a 1-degree per-cell emission inventory (0–360 longitudes) onto
        a 2-degree grid on −180–180; the global total is preserved to
        round-off:

        ```pycon
        >>> import numpy as np, xarray as xr
        >>> from xrops.interpolate import regrid_conservative
        >>> lat, lon = np.arange(-89.5, 90.0, 1.0), np.arange(0.5, 360.0, 1.0)
        >>> rng = np.random.default_rng(0)
        >>> emis = xr.DataArray(
        ...     rng.uniform(0.0, 1.0, size=(lat.size, lon.size)),
        ...     dims=("lat", "lon"),
        ...     coords={"lat": lat, "lon": lon},
        ... )
        >>> target = {
        ...     "lat": np.arange(-89.0, 90.0, 2.0),
        ...     "lon": np.arange(-179.0, 180.0, 2.0),
        ... }
        >>> out = regrid_conservative(emis, target, mode="sum")
        >>> out.sizes
        Frozen({'lat': 90, 'lon': 180})
        >>> bool(np.isclose(float(out.sum()), float(emis.sum()), rtol=1e-12))
        True

        ```
    """
    if geometry not in ("spherical", "planar"):
        raise ValueError(f"geometry must be 'spherical' or 'planar', got {geometry!r}")
    if mode not in ("mean", "sum"):
        raise ValueError(f"mode must be 'mean' or 'sum', got {mode!r}")
    if normalize not in ("destarea", "fracarea"):
        raise ValueError(
            f"normalize must be 'destarea' or 'fracarea', got {normalize!r}"
        )
    dims = tuple(dims)
    if len(dims) != 2:
        raise ValueError(f"dims must name exactly two dimensions, got {dims!r}")
    missing = [d for d in dims if d not in ds.dims]
    if missing:
        raise ValueError(f"input is missing dims {missing!r}; got {tuple(ds.dims)}")

    src_iv = {d: _resolve_intervals(ds, d, bounds, what="input") for d in dims}
    tgt_iv = {
        d: _resolve_intervals(target, d, target_bounds, what="target") for d in dims
    }
    tgt_centres = {d: _centres(target, d, what="target") for d in dims}
    d0, d1 = dims
    if geometry == "spherical":
        src_len = {d0: _sin_lat(src_iv[d0]), d1: np.deg2rad(src_iv[d1])}
        tgt_len = {d0: _sin_lat(tgt_iv[d0]), d1: np.deg2rad(tgt_iv[d1])}
        w0 = _overlap_intervals(src_len[d0], tgt_len[d0])
        w1 = _lon_overlap_weights(src_iv[d1], tgt_iv[d1])
    else:
        src_len, tgt_len = src_iv, tgt_iv
        w0 = _overlap_intervals(src_len[d0], tgt_len[d0])
        w1 = _overlap_intervals(src_len[d1], tgt_len[d1])
    src_area = np.outer(_cell_lengths(src_len[d0]), _cell_lengths(src_len[d1]))
    tgt_area = np.outer(_cell_lengths(tgt_len[d0]), _cell_lengths(tgt_len[d1]))

    # A Dataset output carries the target's CF bounds variables that were
    # actually used (not overridden by ``target_bounds``) alongside the
    # ``bounds`` attr; otherwise the attr would dangle, so it is dropped.
    keep_bounds = isinstance(ds, xr.Dataset)
    tgt_bounds_vars: dict[str, xr.DataArray] = {}
    new_coords = {}
    for d in dims:
        attrs = (
            dict(target[d].attrs)
            if isinstance(target, xr.Dataset | xr.DataArray)
            else {}
        )
        cf_bounds = _cf_bounds_variable(target, d)
        overridden = target_bounds is not None and d in target_bounds
        if keep_bounds and cf_bounds is not None and not overridden:
            tgt_bounds_vars[str(cf_bounds.name)] = cf_bounds
        else:
            attrs.pop("bounds", None)
        new_coords[d] = xr.DataArray(tgt_centres[d], dims=(d,), attrs=attrs)

    def _regrid(da: xr.DataArray) -> xr.DataArray:
        return _regrid_conservative_dataarray(
            da,
            dims=dims,
            w0=w0,
            w1=w1,
            src_area=src_area,
            tgt_area=tgt_area,
            new_coords=new_coords,
            mode=mode,
            normalize=normalize,
            skipna=skipna,
        )

    if isinstance(ds, xr.DataArray):
        return _regrid(ds)

    bounds_vars = {ds[d].attrs.get("bounds") for d in dims} - {None}
    out_vars: dict[str, xr.DataArray] = {}
    for name, da in ds.data_vars.items():
        if name in bounds_vars:
            continue
        present = [d for d in dims if d in da.dims]
        if not present:
            out_vars[str(name)] = da
        elif len(present) == 2:
            out_vars[str(name)] = _regrid(da)
        else:
            raise ValueError(
                f"variable {name!r} carries only {present[0]!r} of the regridded "
                f"dims {dims!r}; drop it or regrid it separately"
            )
    out = xr.Dataset(out_vars, attrs=dict(ds.attrs))
    for name, bnds in tgt_bounds_vars.items():
        out[name] = bnds.variable
    extra_coords = {
        name: coord
        for name, coord in ds.coords.items()
        if name not in out.coords and not set(dims) & set(coord.dims)
    }
    return out.assign_coords(extra_coords) if extra_coords else out

overlap_weights_1d

overlap_weights_1d(src_bounds: Float[ndarray, 'n+1'], tgt_bounds: Float[ndarray, 'm+1']) -> Float[ndarray, 'm n']

Overlap length of each source interval with each target interval.

Entry (j, i) is max(0, min(hi_i, hi_j) - max(lo_i, lo_j)) for source interval i and target interval j. Bounds may be ascending or descending (each interval is normalised to lo < hi), and in any unit — pass sin(lat) bounds to obtain spherical latitude weights.

Parameters:

Name	Type	Description	Default
src_bounds	Float[ndarray, 'n+1']	n + 1 cell edges of the source axis.	required
tgt_bounds	Float[ndarray, 'm+1']	m + 1 cell edges of the target axis.	required

Returns:

Type	Description
Float[ndarray, 'm n']	(m, n) array of non-negative overlap lengths.

Raises:

Type	Description
ValueError	If either bounds array is not 1-D with at least two edges.

Example

>>> import numpy as np
>>> from xrops.interpolate import overlap_weights_1d
>>> src, tgt = np.array([0.0, 1.0, 2.0, 3.0]), np.array([0.0, 1.5, 3.0])
>>> overlap_weights_1d(src, tgt)
array([[1. , 0.5, 0. ],
       [0. , 0.5, 1. ]])

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/grid_to_grid.py

def overlap_weights_1d(
    src_bounds: Float[np.ndarray, "n+1"], tgt_bounds: Float[np.ndarray, "m+1"]
) -> Float[np.ndarray, "m n"]:
    """Overlap length of each source interval with each target interval.

    Entry ``(j, i)`` is ``max(0, min(hi_i, hi_j) - max(lo_i, lo_j))`` for
    source interval ``i`` and target interval ``j``. Bounds may be
    ascending or descending (each interval is normalised to ``lo < hi``),
    and in any unit — pass ``sin(lat)`` bounds to obtain spherical
    latitude weights.

    Args:
        src_bounds: ``n + 1`` cell edges of the source axis.
        tgt_bounds: ``m + 1`` cell edges of the target axis.

    Returns:
        ``(m, n)`` array of non-negative overlap lengths.

    Raises:
        ValueError: If either bounds array is not 1-D with at least two
            edges.

    Example:
        ```pycon
        >>> import numpy as np
        >>> from xrops.interpolate import overlap_weights_1d
        >>> src, tgt = np.array([0.0, 1.0, 2.0, 3.0]), np.array([0.0, 1.5, 3.0])
        >>> overlap_weights_1d(src, tgt)
        array([[1. , 0.5, 0. ],
               [0. , 0.5, 1. ]])

        ```
    """
    src = np.asarray(src_bounds, dtype=float)
    tgt = np.asarray(tgt_bounds, dtype=float)
    if src.ndim != 1 or tgt.ndim != 1 or src.size < 2 or tgt.size < 2:
        raise ValueError(
            "src_bounds and tgt_bounds must be 1-D arrays of at least two edges; "
            f"got shapes {src.shape} and {tgt.shape}"
        )
    return _overlap_intervals(
        _edges_to_intervals(src), _edges_to_intervals(tgt)
    ).toarray()

coarsen

coarsen(ds: Dataset | DataArray, factor: dict[str, int], method: str = 'mean', boundary: str = 'trim') -> Dataset | DataArray

Coarsen ds along one or more dimensions by integer factors.

Thin wrapper around xr.Dataset.coarsen.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/grid_to_grid.py

def coarsen(
    ds: xr.Dataset | xr.DataArray,
    factor: dict[str, int],
    method: str = "mean",
    boundary: str = "trim",
) -> xr.Dataset | xr.DataArray:
    """Coarsen ``ds`` along one or more dimensions by integer factors.

    Thin wrapper around ``xr.Dataset.coarsen``.
    """
    coarsened = ds.coarsen(dim=factor, boundary=boundary)
    return getattr(coarsened, method)()

coarsen_conservative

coarsen_conservative(ds: Dataset | DataArray, factor: Mapping[str, int], *, lat: str = 'lat', boundary: Literal['trim', 'exact', 'pad'] = 'trim') -> Dataset | DataArray

Area-weighted coarsen using cosine-of-latitude weights.

This preserves cosine-latitude-weighted integrals for aligned, integer coarsening of regular latitude/longitude grids. Non-latitude dimensions use uniform weights, and missing values are skipped with weights renormalized within each block. For non-integer factors or an arbitrary rectilinear target grid use :func:regrid_conservative, which reduces to this function for aligned integer factors.

Parameters:

Name	Type	Description	Default
ds	Dataset | DataArray	Dataset or data array to coarsen.	required
factor	Mapping[str, int]	Mapping from dimension name to integer coarsen factor.	required
lat	str	Latitude dimension name. Values are expected to be cell centers in degrees within the usual [-90, 90] latitude range.	'lat'
boundary	Literal['trim', 'exact', 'pad']	Boundary mode forwarded to :meth:xarray.DataArray.coarsen.	'trim'

Returns:

Type	Description
Dataset | DataArray	Coarsened dataset or data array with the same dimension names.

Examples:

>>> coarsen_conservative(da, {"lat": 4, "lon": 4})
>>> coarsen_conservative(ds, {"latitude": 2}, lat="latitude")

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/grid_to_grid.py

def coarsen_conservative(
    ds: xr.Dataset | xr.DataArray,
    factor: Mapping[str, int],
    *,
    lat: str = "lat",
    boundary: Literal["trim", "exact", "pad"] = "trim",
) -> xr.Dataset | xr.DataArray:
    """Area-weighted coarsen using cosine-of-latitude weights.

    This preserves cosine-latitude-weighted integrals for aligned, integer
    coarsening of regular latitude/longitude grids. Non-latitude dimensions
    use uniform weights, and missing values are skipped with weights
    renormalized within each block. For non-integer factors or an
    arbitrary rectilinear target grid use :func:`regrid_conservative`,
    which reduces to this function for aligned integer factors.

    Args:
        ds: Dataset or data array to coarsen.
        factor: Mapping from dimension name to integer coarsen factor.
        lat: Latitude dimension name. Values are expected to be cell centers in
            degrees within the usual [-90, 90] latitude range.
        boundary: Boundary mode forwarded to :meth:`xarray.DataArray.coarsen`.

    Returns:
        Coarsened dataset or data array with the same dimension names.

    Examples:
        ```pycon
        >>> coarsen_conservative(da, {"lat": 4, "lon": 4})
        >>> coarsen_conservative(ds, {"latitude": 2}, lat="latitude")
        ```
    """
    factor_dict = _validate_coarsen_factor(factor)
    if isinstance(ds, xr.Dataset):
        return ds.map(
            lambda da: _coarsen_conservative_dataset_variable(
                da, factor_dict, lat=lat, boundary=boundary
            )
        )
    return _coarsen_conservative_dataarray(ds, factor_dict, lat=lat, boundary=boundary)

refine

refine(ds: Dataset | DataArray, factor: dict[str, int], method: str = 'linear') -> Dataset | DataArray

Refine ds along one or more dimensions by integer factors.

Produces a factor[dim]-times-denser grid along each dimension via :meth:xr.Dataset.interp.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/grid_to_grid.py

def refine(
    ds: xr.Dataset | xr.DataArray,
    factor: dict[str, int],
    method: str = "linear",
) -> xr.Dataset | xr.DataArray:
    """Refine ``ds`` along one or more dimensions by integer factors.

    Produces a ``factor[dim]``-times-denser grid along each dimension
    via :meth:`xr.Dataset.interp`.
    """
    new_coords: dict[str, Sequence[float]] = {}
    for dim, f in factor.items():
        old = ds[dim].values
        if f <= 0:
            raise ValueError(f"refinement factor for {dim!r} must be positive.")
        new_coords[dim] = np.linspace(old.min(), old.max(), (len(old) - 1) * f + 1)
    return ds.interp(new_coords, method=method)

resample_time

resample_time(da: DataArray, freq: str = '1D', method: str = 'mean', time: str = 'time', *, interp_method: Literal['linear', 'nearest', 'cubic'] = 'linear') -> DataArray

Resample a DataArray along the time axis via xarray's built-in resampler.

Per the PR β primitive-flip, this primitive is DataArray-only. The Layer-1 ResampleTime operator handles ds.map-style dispatch when a Dataset is supplied.

Parameters:

Name	Type	Description	Default
da	DataArray	Input DataArray.	required
freq	str	Pandas-style frequency string (e.g. "1D", "6H", "1M").	'1D'
method	str	Zero-argument aggregation method; one of "mean", "sum", "max", "min", "median", "std", "var", "first", "last", "count". Methods that require additional arguments (e.g. reduce, quantile) are intentionally not exposed here. Use "interpolate" for upsampling via xarray's resampler interpolation.	'mean'
time	str	Name of the time dimension.	'time'
interp_method	Literal['linear', 'nearest', 'cubic']	Interpolation method forwarded to xarray when method="interpolate".	'linear'

Returns:

Type	Description
DataArray	Resampled DataArray.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/resample.py

def resample_time(
    da: xr.DataArray,
    freq: str = "1D",
    method: str = "mean",
    time: str = "time",
    *,
    interp_method: Literal["linear", "nearest", "cubic"] = "linear",
) -> xr.DataArray:
    """Resample a DataArray along the time axis via xarray's built-in resampler.

    Per the PR β primitive-flip, this primitive is DataArray-only.
    The Layer-1 ``ResampleTime`` operator handles ``ds.map``-style
    dispatch when a Dataset is supplied.

    Args:
        da: Input DataArray.
        freq: Pandas-style frequency string (e.g. ``"1D"``, ``"6H"``,
            ``"1M"``).
        method: Zero-argument aggregation method; one of ``"mean"``,
            ``"sum"``, ``"max"``, ``"min"``, ``"median"``, ``"std"``,
            ``"var"``, ``"first"``, ``"last"``, ``"count"``. Methods
            that require additional arguments (e.g. ``reduce``, ``quantile``)
            are intentionally not exposed here. Use ``"interpolate"`` for
            upsampling via xarray's resampler interpolation.
        time: Name of the time dimension.
        interp_method: Interpolation method forwarded to xarray when
            ``method="interpolate"``.

    Returns:
        Resampled DataArray.
    """
    if method not in _ALLOWED_RESAMPLE_METHODS | {"interpolate"}:
        raise ValueError(
            f"Unknown resample method {method!r}; expected one of "
            f"{sorted(_ALLOWED_RESAMPLE_METHODS | {'interpolate'})}."
        )

    if method == "interpolate":
        _ensure_datetime_like_time(da, time=time)
        if _target_is_coarser_than_source(da, freq=freq, time=time):
            raise ValueError(
                "resample_time(method='interpolate') only supports upsampling; "
                f"got target frequency {freq!r} coarser than the input."
            )
        return da.resample({time: freq}).interpolate(interp_method)

    resampler = da.resample({time: freq})
    return getattr(resampler, method)()

fillnan_spatial

fillnan_spatial(da: DataArray, method: str = 'linear', lon: str = 'lon', lat: str = 'lat') -> DataArray

Fill NaNs in a 2-D lon/lat field by scattered interpolation.

Operates slice-by-slice along any leading dimensions. Uses :func:scipy.interpolate.griddata over the non-NaN support.

Parameters:

Name	Type	Description	Default
da	DataArray	Input DataArray with at least lon and lat dims.	required
method	str	"linear", "nearest", or "cubic".	'linear'
lon	str	Name of the longitude coordinate.	'lon'
lat	str	Name of the latitude coordinate.	'lat'

Returns:

Type	Description
DataArray	Same-shaped DataArray with NaNs filled where interpolation is
DataArray	possible; points outside the convex hull of non-NaN samples stay
DataArray	NaN (except for method="nearest", which extrapolates).

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/gap_fill.py

def fillnan_spatial(
    da: xr.DataArray,
    method: str = "linear",
    lon: str = "lon",
    lat: str = "lat",
) -> xr.DataArray:
    """Fill NaNs in a 2-D lon/lat field by scattered interpolation.

    Operates slice-by-slice along any leading dimensions. Uses
    :func:`scipy.interpolate.griddata` over the non-NaN support.

    Args:
        da: Input DataArray with at least ``lon`` and ``lat`` dims.
        method: ``"linear"``, ``"nearest"``, or ``"cubic"``.
        lon: Name of the longitude coordinate.
        lat: Name of the latitude coordinate.

    Returns:
        Same-shaped DataArray with NaNs filled where interpolation is
        possible; points outside the convex hull of non-NaN samples stay
        NaN (except for ``method="nearest"``, which extrapolates).
    """
    lon_vals = da[lon].values
    lat_vals = da[lat].values
    lon_grid, lat_grid = np.meshgrid(lon_vals, lat_vals, indexing="xy")
    targets = np.column_stack([lon_grid.ravel(), lat_grid.ravel()])

    def _fill_slice(arr: np.ndarray) -> np.ndarray:
        finite = _finite_mask(arr)
        if finite.all() or not finite.any():
            return arr
        samples = np.column_stack([lon_grid[finite].ravel(), lat_grid[finite].ravel()])
        values = arr[finite].ravel()
        filled = griddata(samples, values, targets, method=method)
        out = arr.copy()
        flat = filled.reshape(arr.shape)
        nan_positions = np.isnan(arr)
        out[nan_positions] = np.where(
            np.isnan(flat[nan_positions]),
            arr[nan_positions],
            flat[nan_positions],
        )
        return out

    return xr.apply_ufunc(
        _fill_slice,
        da,
        input_core_dims=[[lat, lon]],
        output_core_dims=[[lat, lon]],
        vectorize=True,
        dask="parallelized",
        output_dtypes=[_floating_output_dtype(da)],
        dask_gufunc_kwargs={"allow_rechunk": False},
    )

fillnan_temporal

fillnan_temporal(ds: Dataset | DataArray, method: str = 'linear', time: str = 'time', max_gap: Any = None) -> Dataset | DataArray

Interpolate NaNs along the time dimension.

Parameters:

Name	Type	Description	Default
ds	Dataset | DataArray	Input Dataset or DataArray.	required
method	str	Any method accepted by xarray's interpolate_na ("linear", "nearest", "quadratic", "cubic", "spline", etc.).	'linear'
time	str	Name of the time dimension.	'time'
max_gap	Any	Maximum time delta to interpolate across; gaps wider than this stay NaN.	None

Returns:

Type	Description
Dataset | DataArray	Same-shaped container with temporal NaNs interpolated.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/gap_fill.py

def fillnan_temporal(
    ds: xr.Dataset | xr.DataArray,
    method: str = "linear",
    time: str = "time",
    max_gap: Any = None,
) -> xr.Dataset | xr.DataArray:
    """Interpolate NaNs along the time dimension.

    Args:
        ds: Input Dataset or DataArray.
        method: Any method accepted by xarray's ``interpolate_na``
            (``"linear"``, ``"nearest"``, ``"quadratic"``, ``"cubic"``,
            ``"spline"``, etc.).
        time: Name of the time dimension.
        max_gap: Maximum time delta to interpolate across; gaps wider
            than this stay NaN.

    Returns:
        Same-shaped container with temporal NaNs interpolated.
    """
    return ds.interpolate_na(dim=time, method=method, max_gap=max_gap)

fillnan_laplacian

fillnan_laplacian(da: DataArray, *, max_iter: int = 1000, tol: float = 0.0001, relaxation: float = 1.0, boundary: str = 'reflect', lon: str = 'lon', lat: str = 'lat') -> DataArray

Fill NaN cells via iterative harmonic relaxation.

Solves ∇²u = 0 on the masked region with Dirichlet boundary values from finite cells. Iteration stops when the maximum absolute update over missing cells is below tol or after max_iter iterations. All-NaN and all-finite slices pass through unchanged.

Parameters:

Name	Type	Description	Default
da	DataArray	Input DataArray with at least lat and lon dims.	required
max_iter	int	Maximum red/black Gauss-Seidel sweeps.	1000
tol	float	Absolute convergence tolerance for missing-cell updates.	0.0001
relaxation	float	SOR relaxation factor. 1.0 is Gauss-Seidel.	1.0
boundary	str	"reflect" for Neumann edges, or "wrap" to wrap longitude edges while reflecting latitude edges.	'reflect'
lon	str	Name of the longitude dimension.	'lon'
lat	str	Name of the latitude dimension.	'lat'

Returns:

Type	Description
DataArray	Same-shaped DataArray with NaNs filled slice-by-slice along any
DataArray	leading dimensions.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/gap_fill.py

def fillnan_laplacian(
    da: xr.DataArray,
    *,
    max_iter: int = 1000,
    tol: float = 1e-4,
    relaxation: float = 1.0,
    boundary: str = "reflect",
    lon: str = "lon",
    lat: str = "lat",
) -> xr.DataArray:
    """Fill NaN cells via iterative harmonic relaxation.

    Solves ∇²u = 0 on the masked region with Dirichlet boundary values
    from finite cells. Iteration stops when the maximum absolute update
    over missing cells is below ``tol`` or after ``max_iter`` iterations.
    All-NaN and all-finite slices pass through unchanged.

    Args:
        da: Input DataArray with at least ``lat`` and ``lon`` dims.
        max_iter: Maximum red/black Gauss-Seidel sweeps.
        tol: Absolute convergence tolerance for missing-cell updates.
        relaxation: SOR relaxation factor. ``1.0`` is Gauss-Seidel.
        boundary: ``"reflect"`` for Neumann edges, or ``"wrap"`` to wrap
            longitude edges while reflecting latitude edges.
        lon: Name of the longitude dimension.
        lat: Name of the latitude dimension.

    Returns:
        Same-shaped DataArray with NaNs filled slice-by-slice along any
        leading dimensions.
    """
    _validate_laplacian_args(max_iter, tol, relaxation, boundary)

    def _fill_slice(arr: np.ndarray) -> np.ndarray:
        filled, _ = _fillnan_laplacian_slice(
            arr,
            max_iter=max_iter,
            tol=tol,
            relaxation=relaxation,
            boundary=boundary,
        )
        return filled

    return xr.apply_ufunc(
        _fill_slice,
        da,
        input_core_dims=[[lat, lon]],
        output_core_dims=[[lat, lon]],
        vectorize=True,
        dask="parallelized",
        output_dtypes=[_floating_output_dtype(da)],
        dask_gufunc_kwargs={"allow_rechunk": False},
    )

fillnan_biharmonic

fillnan_biharmonic(da: DataArray, *, lon: str = 'lon', lat: str = 'lat', mask: DataArray | None = None, split_into_regions: bool = True) -> DataArray

Fill NaNs in a 2-D lon/lat field by biharmonic inpainting.

Operates slice-by-slice along any leading dimensions and wraps :func:skimage.restoration.inpaint_biharmonic. The optional mask follows scikit-image semantics: True marks pixels to fill.

Parameters:

Name	Type	Description	Default
da	DataArray	Input DataArray with at least lon and lat dims.	required
lon	str	Name of the longitude dimension.	'lon'
lat	str	Name of the latitude dimension.	'lat'
mask	DataArray | None	Optional explicit boolean mask. If None, NaNs in da are filled.	None
split_into_regions	bool	Forwarded to scikit-image. When True, each connected masked region is solved independently.	True

Returns:

Type	Description
DataArray	Same-shaped DataArray with masked pixels inpainted. Fully masked
DataArray	slices pass through unchanged.

Raises:

Type	Description
ImportError	If scikit-image is not installed.
ValueError	If da or mask is missing the spatial dimensions.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/gap_fill.py

def fillnan_biharmonic(
    da: xr.DataArray,
    *,
    lon: str = "lon",
    lat: str = "lat",
    mask: xr.DataArray | None = None,
    split_into_regions: bool = True,
) -> xr.DataArray:
    """Fill NaNs in a 2-D lon/lat field by biharmonic inpainting.

    Operates slice-by-slice along any leading dimensions and wraps
    :func:`skimage.restoration.inpaint_biharmonic`. The optional ``mask``
    follows scikit-image semantics: ``True`` marks pixels to fill.

    Args:
        da: Input DataArray with at least ``lon`` and ``lat`` dims.
        lon: Name of the longitude dimension.
        lat: Name of the latitude dimension.
        mask: Optional explicit boolean mask. If ``None``, NaNs in ``da``
            are filled.
        split_into_regions: Forwarded to scikit-image. When ``True``, each
            connected masked region is solved independently.

    Returns:
        Same-shaped DataArray with masked pixels inpainted. Fully masked
        slices pass through unchanged.

    Raises:
        ImportError: If scikit-image is not installed.
        ValueError: If ``da`` or ``mask`` is missing the spatial dimensions.
    """
    _require_dims(da, lon, lat, name="da")

    if mask is None:
        mask_da = da.isnull()
    else:
        _require_dims(mask, lon, lat, name="mask")
        mask_da = mask.astype(bool)

    output_dtype = _floating_output_dtype(da)

    def _fill_slice(arr: np.ndarray, mask_arr: np.ndarray) -> np.ndarray:
        mask_bool = np.array(mask_arr, dtype=bool, copy=True)
        if not mask_bool.any() or mask_bool.all():
            return arr.astype(output_dtype, copy=True)

        # Only import scikit-image when we actually inpaint, so users without
        # the [image] extra still get fully-finite / empty-mask passthrough.
        inpaint_biharmonic = _require_inpaint_biharmonic()

        # Extend the solver mask to cover any non-finite pixels outside the
        # caller's mask — using NaNs as zero-valued boundary conditions would
        # bias the solve. After the solve we restore those positions back to
        # their original (NaN) value so the caller's mask remains the only
        # set of cells whose values we modify.
        nonfinite_outside_mask = ~_finite_mask(arr) & ~mask_bool
        solver_mask = mask_bool | nonfinite_outside_mask
        arr_filled = arr.astype(np.float64, copy=True)
        # Inside the solver mask the value is unused but must be finite.
        arr_filled[solver_mask] = 0.0
        out = inpaint_biharmonic(
            arr_filled,
            solver_mask,
            split_into_regions=split_into_regions,
            channel_axis=None,
        )
        out = np.asarray(out, dtype=output_dtype)
        # Restore everything outside the caller's mask (including unrelated
        # NaNs and inf values) to its original input value.
        out[~mask_bool] = arr[~mask_bool]
        return out

    return xr.apply_ufunc(
        _fill_slice,
        da,
        mask_da,
        input_core_dims=[[lat, lon], [lat, lon]],
        output_core_dims=[[lat, lon]],
        vectorize=True,
        dask="parallelized",
        output_dtypes=[output_dtype],
        dask_gufunc_kwargs={"allow_rechunk": False},
    )

fillnan_rbf

fillnan_rbf(da: DataArray, kernel: str = 'thin_plate_spline', neighbors: int | None = 32, lon: str = 'lon', lat: str = 'lat') -> DataArray

Fill NaNs using a radial-basis-function interpolator.

Uses :class:scipy.interpolate.RBFInterpolator. More expensive than fillnan_spatial but extrapolates smoothly and respects the global shape of the signal.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/gap_fill.py

def fillnan_rbf(
    da: xr.DataArray,
    kernel: str = "thin_plate_spline",
    neighbors: int | None = 32,
    lon: str = "lon",
    lat: str = "lat",
) -> xr.DataArray:
    """Fill NaNs using a radial-basis-function interpolator.

    Uses :class:`scipy.interpolate.RBFInterpolator`. More expensive than
    ``fillnan_spatial`` but extrapolates smoothly and respects the
    global shape of the signal.
    """
    lon_vals = da[lon].values
    lat_vals = da[lat].values
    lon_grid, lat_grid = np.meshgrid(lon_vals, lat_vals, indexing="xy")

    def _fill_slice(arr: np.ndarray) -> np.ndarray:
        finite = _finite_mask(arr)
        if finite.all() or not finite.any():
            return arr
        samples = np.column_stack([lon_grid[finite].ravel(), lat_grid[finite].ravel()])
        values = arr[finite].ravel()
        rbf = RBFInterpolator(samples, values, kernel=kernel, neighbors=neighbors)
        # Only patch missing positions; leave observed values untouched.
        missing = ~finite
        missing_points = np.column_stack(
            [lon_grid[missing].ravel(), lat_grid[missing].ravel()]
        )
        out = arr.copy()
        out[missing] = rbf(missing_points)
        return out

    return xr.apply_ufunc(
        _fill_slice,
        da,
        input_core_dims=[[lat, lon]],
        output_core_dims=[[lat, lon]],
        vectorize=True,
        dask="parallelized",
        output_dtypes=[_floating_output_dtype(da)],
        dask_gufunc_kwargs={"allow_rechunk": False},
    )

fillnan_idw

fillnan_idw(da: DataArray, *, lon: str = 'lon', lat: str = 'lat', k: int = 8, power: float = 2.0, metric: Metric = 'euclidean', max_distance: float | None = None, eps: float = 1e-12) -> DataArray

Fill NaNs in a 2-D lon/lat field by IDW from finite neighbours.

Operates slice-by-slice along any leading dimensions. For metric="haversine", coordinate values are interpreted as lon/lat degrees and converted to radians internally; max_distance is in radians. For metric="euclidean", distances are measured in array-index space (column/row positions), not coordinate-value space. This differs from :func:idw_to_points and :func:idw_to_grid, whose Euclidean mode uses source and target lon/lat coordinate values directly.

Parameters:

Name	Type	Description	Default
da	DataArray	Input DataArray with lat and lon dimensions.	required
lon	str	Name of the longitude coordinate and dimension.	'lon'
lat	str	Name of the latitude coordinate and dimension.	'lat'
k	int	Number of nearest finite neighbours to use.	8
power	float	Inverse-distance exponent. 0 gives an unweighted kNN mean.	2.0
metric	Metric	"euclidean" for array-index distances or "haversine" for great-circle distances from coordinate values.	'euclidean'
max_distance	float | None	Optional neighbour cutoff. For metric="haversine", this is interpreted in radians; for "euclidean", in index units.	None
eps	float	Small non-negative offset for non-exact distance weights.	1e-12

Returns:

Type	Description
DataArray	Same-shaped DataArray with NaNs filled where finite neighbours are
DataArray	available.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/knn.py

def fillnan_idw(
    da: xr.DataArray,
    *,
    lon: str = "lon",
    lat: str = "lat",
    k: int = 8,
    power: float = 2.0,
    metric: Metric = "euclidean",
    max_distance: float | None = None,
    eps: float = 1e-12,
) -> xr.DataArray:
    """Fill NaNs in a 2-D lon/lat field by IDW from finite neighbours.

    Operates slice-by-slice along any leading dimensions. For
    ``metric="haversine"``, coordinate values are interpreted as lon/lat
    degrees and converted to radians internally; ``max_distance`` is in
    radians. For ``metric="euclidean"``, distances are measured in array-index
    space (column/row positions), not coordinate-value space.
    This differs from :func:`idw_to_points` and :func:`idw_to_grid`, whose
    Euclidean mode uses source and target lon/lat coordinate values directly.

    Args:
        da: Input DataArray with ``lat`` and ``lon`` dimensions.
        lon: Name of the longitude coordinate and dimension.
        lat: Name of the latitude coordinate and dimension.
        k: Number of nearest finite neighbours to use.
        power: Inverse-distance exponent. ``0`` gives an unweighted kNN mean.
        metric: ``"euclidean"`` for array-index distances or ``"haversine"``
            for great-circle distances from coordinate values.
        max_distance: Optional neighbour cutoff. For ``metric="haversine"``,
            this is interpreted in radians; for ``"euclidean"``, in index
            units.
        eps: Small non-negative offset for non-exact distance weights.

    Returns:
        Same-shaped DataArray with NaNs filled where finite neighbours are
        available.
    """
    _validate_idw_args(k, power, metric, max_distance, eps)
    lon_coord = np.asarray(da.coords[lon].values, dtype=float)
    lat_coord = np.asarray(da.coords[lat].values, dtype=float)

    if metric == "haversine":
        lon_grid, lat_grid = np.meshgrid(lon_coord, lat_coord, indexing="xy")
    else:
        rows, cols = np.indices((lat_coord.size, lon_coord.size))

    def _fill_slice(arr: np.ndarray) -> np.ndarray:
        finite = _finite_mask(arr)
        if finite.all() or not finite.any():
            return arr

        if metric == "haversine":
            src_xy = _to_metric_xy(lon_grid[finite], lat_grid[finite], "haversine")
            query_xy = _to_metric_xy(
                lon_grid[~finite],
                lat_grid[~finite],
                "haversine",
            )
            tree = _build_tree(src_xy, "haversine")
        else:
            src_xy = np.column_stack([cols[finite], rows[finite]])
            query_xy = np.column_stack([cols[~finite], rows[~finite]])
            tree = _build_tree(src_xy, "euclidean")

        filled = _idw_kernel(
            tree,
            arr[finite].astype(float, copy=False),
            query_xy,
            k,
            power,
            max_distance,
            eps,
        )
        out = arr.copy()
        out[~finite] = filled
        return out

    return xr.apply_ufunc(
        _fill_slice,
        da,
        input_core_dims=[[lat, lon]],
        output_core_dims=[[lat, lon]],
        vectorize=True,
        dask="parallelized",
        output_dtypes=[_floating_output_dtype(da)],
        dask_gufunc_kwargs={"allow_rechunk": False},
    )

fillnan_climatology

fillnan_climatology(da: DataArray, *, time: str = 'time', group: Literal['month', 'dayofyear', 'season'] = 'month', residual: Literal['zero', 'linear'] = 'linear', min_count: int = 1) -> DataArray

Fill temporal NaNs from a climatological mean.

Parameters:

Name	Type	Description	Default
da	DataArray	Input DataArray with a datetime-like time dimension.	required
time	str	Name of the time dimension.	'time'
group	Literal['month', 'dayofyear', 'season']	Calendar grouping used to compute the climatology.	'month'
residual	Literal['zero', 'linear']	"zero" fills missing values with climatology only; "linear" linearly interpolates the anomaly and adds it back.	'linear'
min_count	int	Minimum finite observations required for each climatology group.	1

Returns:

Type	Description
DataArray	Same-shaped DataArray with fillable temporal NaNs replaced.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/gap_fill.py

def fillnan_climatology(
    da: xr.DataArray,
    *,
    time: str = "time",
    group: Literal["month", "dayofyear", "season"] = "month",
    residual: Literal["zero", "linear"] = "linear",
    min_count: int = 1,
) -> xr.DataArray:
    """Fill temporal NaNs from a climatological mean.

    Args:
        da: Input DataArray with a datetime-like time dimension.
        time: Name of the time dimension.
        group: Calendar grouping used to compute the climatology.
        residual: ``"zero"`` fills missing values with climatology only;
            ``"linear"`` linearly interpolates the anomaly and adds it back.
        min_count: Minimum finite observations required for each climatology
            group.

    Returns:
        Same-shaped DataArray with fillable temporal NaNs replaced.
    """
    if time not in da.dims:
        raise ValueError(f"da must have dim {time!r}.")
    if group not in {"month", "dayofyear", "season"}:
        raise ValueError(
            f"group must be one of {{'month', 'dayofyear', 'season'}}, got {group!r}."
        )
    if residual not in {"zero", "linear"}:
        raise ValueError(f"residual must be 'zero' or 'linear', got {residual!r}.")
    if min_count < 1:
        raise ValueError(f"min_count must be >= 1, got {min_count}.")

    try:
        grouper = getattr(da[time].dt, group)
    except AttributeError as exc:
        raise ValueError(
            f"{time!r} coordinate must be datetime-like for climatology grouping."
        ) from exc
    # Name of the synthetic group coordinate created by xarray (e.g. "month").
    group_name = grouper.name
    grouped = da.groupby(grouper)
    climatology = grouped.mean(time, skipna=True)
    counts = grouped.count(time)
    climatology = climatology.where(counts >= min_count)
    climatology_broadcast = climatology.sel({group_name: grouper})

    if residual == "zero":
        out = da.fillna(climatology_broadcast)
    else:
        anomaly = (da - climatology_broadcast).interpolate_na(dim=time, method="linear")
        out = da.fillna(anomaly + climatology_broadcast)

    # Drop the synthetic groupby coordinate if it was not present on input.
    if group_name not in da.coords and group_name in out.coords:
        out = out.drop_vars(group_name)
    return out

bin_2d

bin_2d(da: DataArray, grid: Grid, statistic: str = 'mean', lon: str = 'lon', lat: str = 'lat') -> DataArray

Bin a scattered DataArray onto a regular lon/lat grid.

Parameters:

Name	Type	Description	Default
da	DataArray	Input DataArray whose non-NaN samples will be binned. Must carry 1-D lon and lat coordinates (typically as sibling arrays of the same length as the data).	required
grid	Grid	Target :class:Grid.	required
statistic	str	"mean", "median", "sum", "count", "min", "max", "std" — anything accepted by :func:scipy.stats.binned_statistic_2d.	'mean'
lon	str	Name of the longitude coordinate on da.	'lon'
lat	str	Name of the latitude coordinate on da.	'lat'

Returns:

Type	Description
DataArray	DataArray on the grid with (lat, lon) dims.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/binning.py

def bin_2d(
    da: xr.DataArray,
    grid: Grid,
    statistic: str = "mean",
    lon: str = "lon",
    lat: str = "lat",
) -> xr.DataArray:
    """Bin a scattered DataArray onto a regular lon/lat grid.

    Args:
        da: Input DataArray whose non-NaN samples will be binned. Must
            carry 1-D ``lon`` and ``lat`` coordinates (typically as
            sibling arrays of the same length as the data).
        grid: Target :class:`Grid`.
        statistic: ``"mean"``, ``"median"``, ``"sum"``, ``"count"``,
            ``"min"``, ``"max"``, ``"std"`` — anything accepted by
            :func:`scipy.stats.binned_statistic_2d`.
        lon: Name of the longitude coordinate on ``da``.
        lat: Name of the latitude coordinate on ``da``.

    Returns:
        DataArray on the grid with ``(lat, lon)`` dims.
    """
    values = np.ravel(da.values)
    lons = np.ravel(np.asarray(da[lon].values))
    lats = np.ravel(np.asarray(da[lat].values))
    finite = _finite_mask(values)

    lon_edges, lat_edges = grid.bin_edges()
    stat, _, _, _ = binned_statistic_2d(
        lons[finite],
        lats[finite],
        values[finite],
        statistic=statistic,
        bins=[lon_edges, lat_edges],
    )
    return xr.DataArray(
        data=stat.T,
        dims=("lat", "lon"),
        coords={"lon": grid.lon, "lat": grid.lat},
        name=da.name,
        attrs=dict(da.attrs),
    )

histogram_2d

histogram_2d(da: DataArray, grid: Grid, lon: str = 'lon', lat: str = 'lat') -> DataArray

2-D histogram count on grid.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/binning.py

def histogram_2d(
    da: xr.DataArray,
    grid: Grid,
    lon: str = "lon",
    lat: str = "lat",
) -> xr.DataArray:
    """2-D histogram count on ``grid``."""
    return bin_2d(da, grid, statistic="count", lon=lon, lat=lat)

kde_to_grid

kde_to_grid(lons: Float[ndarray, 'n'], lats: Float[ndarray, 'n'], grid: Grid, *, weights: Float[ndarray, 'n'] | None = None, bandwidth: float | BandwidthRule = 'scott', kernel: KernelName = 'gaussian', metric: MetricName = 'euclidean', algorithm: AlgorithmName = 'auto', output: OutputMode = 'density', rtol: float = 0.0001) -> DataArray

Evaluate a KDE of scattered lon/lat points on grid.

Parameters:

Name	Type	Description	Default
lons	Float[ndarray, 'n']	1-D point longitudes. Non-finite entries are dropped.	required
lats	Float[ndarray, 'n']	1-D point latitudes. Non-finite entries are dropped.	required
grid	Grid	Target :class:Grid.	required
weights	Float[ndarray, 'n'] | None	Optional per-point weights passed to :class:sklearn.neighbors.KernelDensity. Non-finite weights are dropped with their corresponding lon/lat points.	None
bandwidth	float | BandwidthRule	Positive float, or "scott" / "silverman". For metric="haversine", float bandwidths are radians.	'scott'
kernel	KernelName	sklearn KDE kernel name.	'gaussian'
metric	MetricName	"euclidean" treats lon/lat as planar; "haversine" converts lon/lat to radians and uses great-circle distances.	'euclidean'
algorithm	AlgorithmName	sklearn tree backend. "haversine" uses "ball_tree".	'auto'
output	OutputMode	Output normalization. See semantics below.	'density'
rtol	float	Relative tolerance for sklearn tree pruning.	0.0001

Output modes

"density" is a probability density that integrates to 1 over the grid (per unit area). For metric="haversine" the output is renormalized to integrate to 1 on the sphere using cos(lat) * dlat * dlon (steradian) cell areas, since sklearn's KernelDensity only normalizes correctly for the Euclidean metric.

"counts_per_area" is the expected per-unit-area count (density * n_eff). Integrating over the grid recovers the input point count (or weight sum). Comparable across grids of different resolutions.

"counts" is the expected per-cell count (density * n_eff * cell_area). Summing over the grid recovers the input point count (or weight sum). Scales with cell size, so it is only comparable across grids of identical resolution.

Returns:

Type	Description
DataArray	DataArray on grid with ("lat", "lon") dimensions.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/points_to_grid.py

def kde_to_grid(
    lons: Float[np.ndarray, "n"],
    lats: Float[np.ndarray, "n"],
    grid: Grid,
    *,
    weights: Float[np.ndarray, "n"] | None = None,
    bandwidth: float | BandwidthRule = "scott",
    kernel: KernelName = "gaussian",
    metric: MetricName = "euclidean",
    algorithm: AlgorithmName = "auto",
    output: OutputMode = "density",
    rtol: float = 1e-4,
) -> xr.DataArray:
    """Evaluate a KDE of scattered lon/lat points on ``grid``.

    Args:
        lons: 1-D point longitudes. Non-finite entries are dropped.
        lats: 1-D point latitudes. Non-finite entries are dropped.
        grid: Target :class:`Grid`.
        weights: Optional per-point weights passed to
            :class:`sklearn.neighbors.KernelDensity`. Non-finite weights are
            dropped with their corresponding lon/lat points.
        bandwidth: Positive float, or ``"scott"`` / ``"silverman"``. For
            ``metric="haversine"``, float bandwidths are radians.
        kernel: sklearn KDE kernel name.
        metric: ``"euclidean"`` treats lon/lat as planar; ``"haversine"``
            converts lon/lat to radians and uses great-circle distances.
        algorithm: sklearn tree backend. ``"haversine"`` uses ``"ball_tree"``.
        output: Output normalization. See semantics below.
        rtol: Relative tolerance for sklearn tree pruning.

    Output modes:
        ``"density"`` is a probability density that integrates to 1 over the
        grid (per unit area). For ``metric="haversine"`` the output is
        renormalized to integrate to 1 on the sphere using
        ``cos(lat) * dlat * dlon`` (steradian) cell areas, since sklearn's
        ``KernelDensity`` only normalizes correctly for the Euclidean metric.

        ``"counts_per_area"`` is the expected per-unit-area count
        (``density * n_eff``). Integrating over the grid recovers the input
        point count (or weight sum). Comparable across grids of different
        resolutions.

        ``"counts"`` is the expected per-cell count
        (``density * n_eff * cell_area``). Summing over the grid recovers
        the input point count (or weight sum). Scales with cell size, so it
        is only comparable across grids of identical resolution.

    Returns:
        DataArray on ``grid`` with ``("lat", "lon")`` dimensions.
    """
    if output not in ("density", "counts", "counts_per_area"):
        raise ValueError(f"unknown output mode {output!r}")
    if metric not in ("euclidean", "haversine"):
        raise ValueError(f"unknown metric {metric!r}")
    if grid.lon.size < 2 or grid.lat.size < 2:
        raise ValueError(
            "kde_to_grid requires grid.lon and grid.lat to each have at least "
            f"2 points; got {grid.lon.size} lon, {grid.lat.size} lat"
        )

    lon_values = np.ravel(np.asarray(lons, dtype=float))
    lat_values = np.ravel(np.asarray(lats, dtype=float))
    if lon_values.shape != lat_values.shape:
        raise ValueError("lons and lats must have the same shape")

    finite = _finite_mask(lon_values, lat_values)
    if weights is None:
        sample_weight = None
    else:
        weight_values = np.ravel(np.asarray(weights, dtype=float))
        if weight_values.shape != lon_values.shape:
            raise ValueError("weights must have the same shape as lons and lats")
        finite &= _finite_mask(weight_values)
        sample_weight = weight_values[finite]

    pts = np.column_stack([lon_values[finite], lat_values[finite]])
    if pts.shape[0] < 2:
        raise ValueError("KDE requires at least 2 finite points")

    if metric == "haversine":
        fit_pts = _to_metric_xy(pts[:, 0], pts[:, 1], "haversine")
        algorithm = "ball_tree"
    else:
        fit_pts = pts

    if isinstance(bandwidth, str):
        bandwidth_value = _bandwidth_rule(fit_pts, bandwidth)
    else:
        bandwidth_value = float(bandwidth)
        if bandwidth_value <= 0 or not np.isfinite(bandwidth_value):
            raise ValueError(f"bandwidth must be positive, got {bandwidth!r}")

    neighbors = _require_optional(
        "sklearn.neighbors",
        extra="sklearn",
        feature="kernel density gridding",
        package="scikit-learn",
    )
    kde = neighbors.KernelDensity(
        bandwidth=bandwidth_value,
        kernel=kernel,
        metric=metric,
        algorithm=algorithm,
        rtol=rtol,
    ).fit(fit_pts, sample_weight=sample_weight)

    lon_grid, lat_grid = np.meshgrid(grid.lon, grid.lat, indexing="xy")
    if metric == "haversine":
        queries = _to_metric_xy(lon_grid.ravel(), lat_grid.ravel(), "haversine")
    else:
        queries = np.column_stack([lon_grid.ravel(), lat_grid.ravel()])

    density = np.exp(kde.score_samples(queries)).reshape(len(grid.lat), len(grid.lon))

    # Cell area: deg² for Euclidean, steradian (cos(lat)*dlat*dlon, radians²)
    # for haversine. sklearn normalizes assuming flat (lat, lon) integration in
    # haversine mode, so we re-normalize on the sphere to recover ∫p dA = 1.
    if metric == "haversine":
        dlat = float(np.mean(np.abs(np.diff(np.deg2rad(np.asarray(grid.lat))))))
        dlon = float(np.mean(np.abs(np.diff(np.deg2rad(np.asarray(grid.lon))))))
        cell_area = (
            np.cos(np.deg2rad(np.asarray(grid.lat, dtype=float)))[:, None] * dlat * dlon
        )
        total = float((density * cell_area).sum())
        if total > 0:
            density = density / total
    else:
        dlat = float(np.mean(np.abs(np.diff(np.asarray(grid.lat, dtype=float)))))
        dlon = float(np.mean(np.abs(np.diff(np.asarray(grid.lon, dtype=float)))))
        cell_area = np.full(density.shape, dlat * dlon)

    n_eff = pts.shape[0] if sample_weight is None else float(sample_weight.sum())
    if output == "density":
        data = density
    elif output == "counts_per_area":
        data = density * n_eff
    else:  # "counts" — per-cell expected count
        data = density * n_eff * cell_area

    return xr.DataArray(
        data=data,
        dims=("lat", "lon"),
        coords={"lon": grid.lon, "lat": grid.lat},
    )

points_to_grid

points_to_grid(lons: Float[ndarray, 'n'], lats: Float[ndarray, 'n'], values: Float[ndarray, 'n'], grid: Grid, statistic: str = 'mean') -> DataArray

Bin raw (lon, lat, value) tuples onto grid.

Thin wrapper around :func:scipy.stats.binned_statistic_2d that doesn't require constructing a scattered DataArray first.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/points_to_grid.py

def points_to_grid(
    lons: Float[np.ndarray, "n"],
    lats: Float[np.ndarray, "n"],
    values: Float[np.ndarray, "n"],
    grid: Grid,
    statistic: str = "mean",
) -> xr.DataArray:
    """Bin raw (lon, lat, value) tuples onto ``grid``.

    Thin wrapper around :func:`scipy.stats.binned_statistic_2d` that
    doesn't require constructing a scattered DataArray first.
    """
    finite = np.isfinite(values)
    lon_edges, lat_edges = grid.bin_edges()
    stat, _, _, _ = binned_statistic_2d(
        np.ravel(lons)[finite],
        np.ravel(lats)[finite],
        np.ravel(values)[finite],
        statistic=statistic,
        bins=[lon_edges, lat_edges],
    )
    return xr.DataArray(
        data=stat.T,
        dims=("lat", "lon"),
        coords={"lon": grid.lon, "lat": grid.lat},
    )

idw_to_grid

idw_to_grid(lons: Float[ndarray, '...'], lats: Float[ndarray, '...'], values: Float[ndarray, '...'], grid: Grid, *, k: int = 8, power: float = 2.0, metric: Metric = 'euclidean', max_distance: float | None = None, eps: float = 1e-12) -> DataArray

Inverse-distance interpolate scattered points onto grid.

For metric="haversine", max_distance is interpreted in radians (for example, 50 / 6371 for approximately 50 km).

Parameters:

Name	Type	Description	Default
lons	Float[ndarray, '...']	Source longitudes.	required
lats	Float[ndarray, '...']	Source latitudes.	required
values	Float[ndarray, '...']	Source values.	required
grid	Grid	Target regular lon/lat grid.	required
k	int	Number of nearest finite source neighbours to use.	8
power	float	Inverse-distance exponent. 0 gives an unweighted kNN mean.	2.0
metric	Metric	"euclidean" for degree-space distances or "haversine" for great-circle distances.	'euclidean'
max_distance	float | None	Optional neighbour cutoff. For metric="haversine", this is interpreted in radians.	None
eps	float	Small non-negative offset for non-exact distance weights.	1e-12

Returns:

Type	Description
DataArray	DataArray with ("lat", "lon") dimensions on grid.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/knn.py

def idw_to_grid(
    lons: Float[np.ndarray, "..."],
    lats: Float[np.ndarray, "..."],
    values: Float[np.ndarray, "..."],
    grid: Grid,
    *,
    k: int = 8,
    power: float = 2.0,
    metric: Metric = "euclidean",
    max_distance: float | None = None,
    eps: float = 1e-12,
) -> xr.DataArray:
    """Inverse-distance interpolate scattered points onto ``grid``.

    For ``metric="haversine"``, ``max_distance`` is interpreted in radians
    (for example, ``50 / 6371`` for approximately 50 km).

    Args:
        lons: Source longitudes.
        lats: Source latitudes.
        values: Source values.
        grid: Target regular lon/lat grid.
        k: Number of nearest finite source neighbours to use.
        power: Inverse-distance exponent. ``0`` gives an unweighted kNN mean.
        metric: ``"euclidean"`` for degree-space distances or ``"haversine"``
            for great-circle distances.
        max_distance: Optional neighbour cutoff. For ``metric="haversine"``,
            this is interpreted in radians.
        eps: Small non-negative offset for non-exact distance weights.

    Returns:
        DataArray with ``("lat", "lon")`` dimensions on ``grid``.
    """
    lon_grid, lat_grid = np.meshgrid(grid.lon, grid.lat, indexing="xy")
    out = idw_to_points(
        lons,
        lats,
        values,
        lon_grid,
        lat_grid,
        k=k,
        power=power,
        metric=metric,
        max_distance=max_distance,
        eps=eps,
    )
    return xr.DataArray(
        out,
        dims=("lat", "lon"),
        coords={"lat": grid.lat, "lon": grid.lon},
    )

idw_to_points

idw_to_points(src_lons: Float[ndarray, '...'], src_lats: Float[ndarray, '...'], src_values: Float[ndarray, '...'], dst_lons: Float[ndarray, '...'], dst_lats: Float[ndarray, '...'], *, k: int = 8, power: float = 2.0, metric: Metric = 'euclidean', max_distance: float | None = None, eps: float = 1e-12) -> Float[ndarray, '...']

Inverse-distance interpolate scattered samples onto target points.

Parameters:

Name	Type	Description	Default
src_lons	Float[ndarray, '...']	Source longitudes.	required
src_lats	Float[ndarray, '...']	Source latitudes.	required
src_values	Float[ndarray, '...']	Source values.	required
dst_lons	Float[ndarray, '...']	Target longitudes.	required
dst_lats	Float[ndarray, '...']	Target latitudes.	required
k	int	Number of nearest finite source neighbours to use.	8
power	float	Inverse-distance exponent. 0 gives an unweighted kNN mean.	2.0
metric	Metric	"euclidean" for degree-space distances or "haversine" for great-circle distances.	'euclidean'
max_distance	float | None	Optional neighbour cutoff. For metric="haversine", this is interpreted in radians (approximately km / 6371).	None
eps	float	Small non-negative offset for non-exact distance weights.	1e-12

Returns:

Type	Description
Float[ndarray, '...']	NumPy array matching the broadcast shape of dst_lons and
Float[ndarray, '...']	dst_lats.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/knn.py

def idw_to_points(
    src_lons: Float[np.ndarray, "..."],
    src_lats: Float[np.ndarray, "..."],
    src_values: Float[np.ndarray, "..."],
    dst_lons: Float[np.ndarray, "..."],
    dst_lats: Float[np.ndarray, "..."],
    *,
    k: int = 8,
    power: float = 2.0,
    metric: Metric = "euclidean",
    max_distance: float | None = None,
    eps: float = 1e-12,
) -> Float[np.ndarray, "..."]:
    """Inverse-distance interpolate scattered samples onto target points.

    Args:
        src_lons: Source longitudes.
        src_lats: Source latitudes.
        src_values: Source values.
        dst_lons: Target longitudes.
        dst_lats: Target latitudes.
        k: Number of nearest finite source neighbours to use.
        power: Inverse-distance exponent. ``0`` gives an unweighted kNN mean.
        metric: ``"euclidean"`` for degree-space distances or ``"haversine"``
            for great-circle distances.
        max_distance: Optional neighbour cutoff. For ``metric="haversine"``,
            this is interpreted in radians (approximately ``km / 6371``).
        eps: Small non-negative offset for non-exact distance weights.

    Returns:
        NumPy array matching the broadcast shape of ``dst_lons`` and
        ``dst_lats``.
    """
    _validate_idw_args(k, power, metric, max_distance, eps)
    src_lon, src_lat, src_value = _prepare_sources(src_lons, src_lats, src_values)
    dst_lon, dst_lat = np.broadcast_arrays(
        np.asarray(dst_lons, dtype=float),
        np.asarray(dst_lats, dtype=float),
    )

    out = np.full(dst_lon.shape, np.nan, dtype=float)
    finite_queries = _finite_mask(dst_lon, dst_lat)
    if not finite_queries.any():
        return out

    src_xy = _to_metric_xy(src_lon, src_lat, metric)
    query_xy = _to_metric_xy(dst_lon[finite_queries], dst_lat[finite_queries], metric)
    tree = _build_tree(src_xy, metric)
    out[finite_queries] = _idw_kernel(
        tree,
        src_value,
        query_xy,
        k,
        power,
        max_distance,
        eps,
    )
    return out

gaussian_smooth

gaussian_smooth(da: DataArray, *, dim: str, sigma: float, truncate: float = 4.0) -> DataArray

Gaussian smoothing along dim with standard deviation sigma.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/smooth.py

def gaussian_smooth(
    da: xr.DataArray,
    *,
    dim: str,
    sigma: float,
    truncate: float = 4.0,
) -> xr.DataArray:
    """Gaussian smoothing along ``dim`` with standard deviation ``sigma``."""

    def _fn(arr: np.ndarray, *, axis: int) -> np.ndarray:
        return _kernels.gaussian_smooth(arr, axis=axis, sigma=sigma, truncate=truncate)

    return _apply_kernel_to_dataarray(da, dim, _fn)

gaussian_smooth_masked

gaussian_smooth_masked(da: DataArray, *, dim: str | Sequence[str], sigma: float | Mapping[str, float], truncate: float = 4.0, mode: str = 'reflect', nan_aware: bool = True, min_weight: float = 1e-06) -> DataArray

NaN-aware N-D Gaussian smoothing along one or more dimensions.

Operates on a single :class:xr.DataArray. The Layer-1 GaussianSmoothMasked operator handles the Dataset loop.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/smooth.py

def gaussian_smooth_masked(
    da: xr.DataArray,
    *,
    dim: str | Sequence[str],
    sigma: float | Mapping[str, float],
    truncate: float = 4.0,
    mode: str = "reflect",
    nan_aware: bool = True,
    min_weight: float = 1e-6,
) -> xr.DataArray:
    """NaN-aware N-D Gaussian smoothing along one or more dimensions.

    Operates on a single :class:`xr.DataArray`. The Layer-1
    ``GaussianSmoothMasked`` operator handles the Dataset loop.
    """
    if truncate <= 0:
        raise ValueError(f"truncate must be > 0, got {truncate}")
    if min_weight < 0:
        raise ValueError(f"min_weight must be >= 0, got {min_weight}")

    dims = _normalize_dims(dim)
    sigmas = _normalize_sigmas(sigma, dims)

    missing = tuple(d for d in dims if d not in da.dims)
    if missing:
        raise ValueError(f"dim entries {missing} not in DataArray dims {da.dims}")
    return _gaussian_smooth_masked_dataarray(
        da,
        dims=dims,
        sigmas=sigmas,
        truncate=truncate,
        mode=mode,
        nan_aware=nan_aware,
        min_weight=min_weight,
    )

moving_average

moving_average(da: DataArray, *, dim: str, window: int, center: bool = True, min_periods: int | None = None) -> DataArray

Sliding-window mean along dim.

Computes a centred (center=True) or trailing sliding-window mean of length window along dim. min_periods controls the minimum number of valid (non-NaN) samples required inside the window for the output to be non-NaN; if unset it defaults to window. Inputs with fewer than min_periods valid samples in a window produce NaN at that position.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/smooth.py

def moving_average(
    da: xr.DataArray,
    *,
    dim: str,
    window: int,
    center: bool = True,
    min_periods: int | None = None,
) -> xr.DataArray:
    """Sliding-window mean along ``dim``.

    Computes a centred (``center=True``) or trailing sliding-window mean
    of length ``window`` along ``dim``. ``min_periods`` controls the
    minimum number of valid (non-NaN) samples required inside the window
    for the output to be non-NaN; if unset it defaults to ``window``.
    Inputs with fewer than ``min_periods`` valid samples in a window
    produce NaN at that position.
    """

    def _fn(arr: np.ndarray, *, axis: int) -> np.ndarray:
        return _kernels.moving_average(
            arr,
            axis=axis,
            window=window,
            center=center,
            min_periods=min_periods,
        )

    return _apply_kernel_to_dataarray(da, dim, _fn)

lowpass_filter

lowpass_filter(da: DataArray, *, dim: str, cutoff: float | tuple[float, float] | list[float] | ndarray, order: int = 4, btype: str = 'low') -> DataArray

Zero-phase Butterworth filter along dim.

cutoff is the normalized critical frequency (fraction of the Nyquist rate). For btype in {"bandpass", "bandstop"} pass a length-2 (low, high) sequence. The filter is applied with scipy.signal.sosfiltfilt (forward-backward) for zero phase distortion; order is the per-direction order of the SOS sections.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/smooth.py

def lowpass_filter(
    da: xr.DataArray,
    *,
    dim: str,
    cutoff: float | tuple[float, float] | list[float] | np.ndarray,
    order: int = 4,
    btype: str = "low",
) -> xr.DataArray:
    """Zero-phase Butterworth filter along ``dim``.

    ``cutoff`` is the normalized critical frequency (fraction of the
    Nyquist rate). For ``btype`` in ``{"bandpass", "bandstop"}`` pass a
    length-2 ``(low, high)`` sequence. The filter is applied with
    ``scipy.signal.sosfiltfilt`` (forward-backward) for zero phase
    distortion; ``order`` is the per-direction order of the SOS sections.
    """

    def _fn(arr: np.ndarray, *, axis: int) -> np.ndarray:
        return _kernels.lowpass_filter(
            arr, axis=axis, cutoff=cutoff, order=order, btype=btype
        )

    return _apply_kernel_to_dataarray(da, dim, _fn)

fir_filter

fir_filter(da: DataArray, *, dim: str, cutoff: float | tuple[float, float] | list[float] | ndarray, method: str = 'lanczos', btype: str = 'low', num_taps: int | None = None, attenuation_db: float | None = None) -> DataArray

Zero-phase FIR filter along dim.

cutoff is the normalized critical frequency (fraction of the Nyquist rate); for btype in {"bandpass", "bandstop"} pass a length-2 (low, high) sequence. method selects the window family used to taper the ideal sinc response — "lanczos" or "kaiser". num_taps is an odd FIR tap count; if omitted, a conservative default is chosen (Lanczos picks from cutoff; Kaiser estimates from attenuation_db). attenuation_db is the Kaiser stop-band attenuation target in decibels (default 60.0) and is ignored when method="lanczos". The filter is applied with scipy.signal.filtfilt for zero phase.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/smooth.py

def fir_filter(
    da: xr.DataArray,
    *,
    dim: str,
    cutoff: float | tuple[float, float] | list[float] | np.ndarray,
    method: str = "lanczos",
    btype: str = "low",
    num_taps: int | None = None,
    attenuation_db: float | None = None,
) -> xr.DataArray:
    """Zero-phase FIR filter along ``dim``.

    ``cutoff`` is the normalized critical frequency (fraction of the
    Nyquist rate); for ``btype`` in ``{"bandpass", "bandstop"}`` pass a
    length-2 ``(low, high)`` sequence. ``method`` selects the window
    family used to taper the ideal sinc response — ``"lanczos"`` or
    ``"kaiser"``. ``num_taps`` is an odd FIR tap count; if omitted, a
    conservative default is chosen (Lanczos picks from ``cutoff``;
    Kaiser estimates from ``attenuation_db``). ``attenuation_db`` is the
    Kaiser stop-band attenuation target in decibels (default ``60.0``)
    and is ignored when ``method="lanczos"``. The filter is applied with
    ``scipy.signal.filtfilt`` for zero phase.
    """

    def _fn(arr: np.ndarray, *, axis: int) -> np.ndarray:
        return _kernels.fir_filter(
            arr,
            axis=axis,
            cutoff=cutoff,
            method=method,
            btype=btype,
            num_taps=num_taps,
            attenuation_db=attenuation_db,
        )

    return _apply_kernel_to_dataarray(da, dim, _fn)

sample_at_points

sample_at_points(da: DataArray, points: Dataset | Mapping[str, object], *, coords: tuple[str, ...] = ('lat', 'lon'), method: Method = 'linear', fill_value: float | None = nan, bounds_error: bool = False, point_dim: str = 'points') -> DataArray

Sample a gridded field at scattered target points using RGI.

Parameters:

Name	Type	Description	Default
da	DataArray	Source data with 1-D, strictly monotone coordinates for every interpolation dimension in coords.	required
points	Dataset | Mapping[str, object]	Target locations as a dataset or mapping. Each coordinate in coords must be a 1-D variable on point_dim.	required
coords	tuple[str, ...]	Source/target coordinate names to interpolate over.	('lat', 'lon')
method	Method	Interpolation method supported by :class:scipy.interpolate.RegularGridInterpolator.	'linear'
fill_value	float | None	Value used for finite out-of-bounds points when bounds_error is false. None enables scipy extrapolation.	nan
bounds_error	bool	Whether finite out-of-bounds points raise an error.	False
point_dim	str	Name of the target-point dimension.	'points'

Returns:

Type	Description
DataArray	Values sampled at points with all non-interpolated source
DataArray	dimensions retained and point_dim appended.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/grid_to_points.py

def sample_at_points(
    da: xr.DataArray,
    points: xr.Dataset | Mapping[str, object],
    *,
    coords: tuple[str, ...] = ("lat", "lon"),
    method: Method = "linear",
    fill_value: float | None = np.nan,
    bounds_error: bool = False,
    point_dim: str = "points",
) -> xr.DataArray:
    """Sample a gridded field at scattered target points using RGI.

    Args:
        da: Source data with 1-D, strictly monotone coordinates for every
            interpolation dimension in ``coords``.
        points: Target locations as a dataset or mapping. Each coordinate in
            ``coords`` must be a 1-D variable on ``point_dim``.
        coords: Source/target coordinate names to interpolate over.
        method: Interpolation method supported by
            :class:`scipy.interpolate.RegularGridInterpolator`.
        fill_value: Value used for finite out-of-bounds points when
            ``bounds_error`` is false. ``None`` enables scipy extrapolation.
        bounds_error: Whether finite out-of-bounds points raise an error.
        point_dim: Name of the target-point dimension.

    Returns:
        Values sampled at ``points`` with all non-interpolated source
        dimensions retained and ``point_dim`` appended.
    """
    if not coords:
        raise ValueError("coords must contain at least one interpolation dimension")

    coords = tuple(coords)
    points_ds = _normalize_points(points, coords=coords, point_dim=point_dim)
    _check_dask_chunks(da, coords)

    grid_axes: list[np.ndarray] = []
    normalized = da
    for name in coords:
        normalized, axis = _normalize_axis(normalized, name)
        grid_axes.append(axis)

    point_columns: list[np.ndarray] = []
    valid = np.ones(points_ds.sizes[point_dim], dtype=bool)
    for name in coords:
        numeric, finite = _as_numeric_points(points_ds[name].values, name=name)
        point_columns.append(numeric)
        valid &= finite
    interp_pts = np.column_stack(point_columns)

    if np.issubdtype(normalized.dtype, np.floating) or np.issubdtype(
        normalized.dtype,
        np.complexfloating,
    ):
        result_dtype = normalized.dtype
    else:
        result_dtype = np.dtype(np.float64)

    def _interp_block(arr: np.ndarray) -> np.ndarray:
        out = np.full(interp_pts.shape[0], np.nan, dtype=result_dtype)
        if not np.any(valid):
            return out
        rgi = RegularGridInterpolator(
            tuple(grid_axes),
            arr,
            method=method,
            bounds_error=bounds_error,
            fill_value=fill_value,
        )
        out[valid] = rgi(interp_pts[valid])
        return out

    out = xr.apply_ufunc(
        _interp_block,
        normalized,
        input_core_dims=[list(coords)],
        output_core_dims=[[point_dim]],
        exclude_dims=set(coords),
        vectorize=True,
        dask="parallelized",
        output_dtypes=[result_dtype],
        dask_gufunc_kwargs={
            "output_sizes": {point_dim: interp_pts.shape[0]},
            "allow_rechunk": False,
        },
    )
    # `out.dims` keeps the source dims that are NOT being interpolated over;
    # skip those names when collecting per-point coords so we don't try to
    # attach a 1-D point coord on top of an existing output dim coord.
    retained = tuple(str(d) for d in out.dims if d != point_dim)
    return out.assign_coords(
        _collect_point_coordinates(
            points_ds,
            coords,
            point_dim,
            retained_dims=retained,
        )
    )

along_track

along_track(da: DataArray, track: Dataset, *, coords: tuple[str, ...] = ('time', 'lat', 'lon'), method: Method = 'linear', fill_value: float | None = nan, point_dim: str = 'points') -> DataArray

Collocate a gridded field onto an observation track.

This is a convenience wrapper around :func:sample_at_points for time-collocated tracks such as satellite altimetry, drifters, or buoys. It validates that datetime-like coordinates are datetime-like in both the source and target before delegating to the generic sampler.

Parameters:

Name	Type	Description	Default
da	DataArray	Source gridded field with 1-D coordinates for every name in coords.	required
track	Dataset	Track dataset containing each coordinate in coords as a 1-D variable on point_dim.	required
coords	tuple[str, ...]	Coordinate names to interpolate over. Defaults to ("time", "lat", "lon").	('time', 'lat', 'lon')
method	Method	Interpolation method supported by :class:scipy.interpolate.RegularGridInterpolator.	'linear'
fill_value	float | None	Value used for finite out-of-bounds track points.	nan
point_dim	str	Name of the target-track dimension.	'points'

Returns:

Type	Description
DataArray	Values from da collocated onto track along point_dim.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/grid_to_points.py

def along_track(
    da: xr.DataArray,
    track: xr.Dataset,
    *,
    coords: tuple[str, ...] = ("time", "lat", "lon"),
    method: Method = "linear",
    fill_value: float | None = np.nan,
    point_dim: str = "points",
) -> xr.DataArray:
    """Collocate a gridded field onto an observation track.

    This is a convenience wrapper around :func:`sample_at_points` for
    time-collocated tracks such as satellite altimetry, drifters, or buoys.
    It validates that datetime-like coordinates are datetime-like in both the
    source and target before delegating to the generic sampler.

    Args:
        da: Source gridded field with 1-D coordinates for every name in
            ``coords``.
        track: Track dataset containing each coordinate in ``coords`` as a
            1-D variable on ``point_dim``.
        coords: Coordinate names to interpolate over. Defaults to
            ``("time", "lat", "lon")``.
        method: Interpolation method supported by
            :class:`scipy.interpolate.RegularGridInterpolator`.
        fill_value: Value used for finite out-of-bounds track points.
        point_dim: Name of the target-track dimension.

    Returns:
        Values from ``da`` collocated onto ``track`` along ``point_dim``.
    """
    for name in coords:
        # Validate presence first so the user gets a clear ValueError instead
        # of a bare KeyError from da[name] / track[name] indexing.
        if name not in da.coords and name not in da.dims:
            raise ValueError(f"da is missing coord {name!r}")
        if name not in track:
            raise ValueError(f"track is missing variable {name!r}")
        source_is_time = np.issubdtype(da[name].values.dtype, np.datetime64)
        target_is_time = np.issubdtype(
            track[name].values.dtype,
            np.datetime64,
        )
        if source_is_time != target_is_time:
            raise ValueError(
                f"coord {name!r} must be datetime-like in both da and track, "
                "or numeric in both"
            )
    return sample_at_points(
        da,
        track,
        coords=coords,
        method=method,
        fill_value=fill_value,
        bounds_error=False,
        point_dim=point_dim,
    )

remap_axis

remap_axis(da: DataArray, *, source_dim: str, target_coords: DataArray | ndarray, target_name: str | None = None, method: str = 'linear', source_coords: str | DataArray | None = None, extrapolate: Extrapolate = 'nan') -> DataArray

Remap da along source_dim onto target_coords.

If target_coords is a :class:xr.DataArray and target_name is None, the new dim name is taken from target_coords.name; otherwise target_name (or source_dim) is used.

Every coordinate that depends on source_dim is dropped from the output; the new target coordinate is attached under the new name. The Layer-1 RemapAxis operator handles Dataset selection / per-variable looping.

Parameters:

Name	Type	Description	Default
da	DataArray	Numeric DataArray carrying source_dim.	required
source_dim	str	Name of the dimension being replaced.	required
target_coords	DataArray | ndarray	1-D target axis values.	required
target_name	str | None	Explicit name for the new dimension.	None
method	str	"linear" or "nearest".	'linear'
source_coords	str | DataArray | None	Source axis values. None (default) uses the 1-D dimension coordinate da[source_dim]. A coordinate name on da or a DataArray broadcastable against da supplies per-column levels — it must carry source_dim and may vary along any of the other dims (e.g. z_agl(time, level, y, x) on a terrain-following grid). Every column must be strictly monotone along source_dim.	None
extrapolate	Extrapolate	"nan" (default) — targets outside a column's range become NaN; "nearest" — they hold that column's end value.	'nan'

Returns:

Type	Description
DataArray	da with source_dim replaced by the target axis. When
DataArray	source_coords is given, chunked (dask) inputs stay lazy as
DataArray	long as source_dim sits in a single chunk.

Raises:

Type	Description
ValueError	If source_dim is missing, source_coords does not carry it, a column is not monotone, or extrapolate / method is unknown.
TypeError	If da is not numeric.

Example

Remap a model-level field with a 4-D height-above-ground coordinate onto fixed heights, holding the lowest / highest model level beyond each column's range:

>>> import numpy as np, xarray as xr
>>> from xrops.transforms import remap_axis
>>> level = np.arange(4)
>>> z_agl = 20.0 * (level[None, :, None, None] + 1) + np.zeros((2, 4, 3, 3))
>>> z_agl[..., 0, 0] += 15.0  # one column starts higher up
>>> t = xr.DataArray(
...     300.0 - 0.01 * z_agl,
...     dims=("time", "level", "y", "x"),
...     coords={"z_agl": (("time", "level", "y", "x"), z_agl)},
... )
>>> out = remap_axis(
...     t,
...     source_dim="level",
...     target_coords=np.array([10.0, 50.0, 90.0]),
...     target_name="height",
...     source_coords="z_agl",
...     extrapolate="nearest",
... )
>>> out.dims, out.sizes["height"]
(('time', 'height', 'y', 'x'), 3)
>>> float(out.isel(time=0, height=1, y=1, x=1))  # in range: linear
299.5

Source code in packages/xrtoolz-ops/src/xrops/transforms/_src/coord_remap.py

def remap_axis(
    da: xr.DataArray,
    *,
    source_dim: str,
    target_coords: xr.DataArray | np.ndarray,
    target_name: str | None = None,
    method: str = "linear",
    source_coords: str | xr.DataArray | None = None,
    extrapolate: Extrapolate = "nan",
) -> xr.DataArray:
    """Remap ``da`` along ``source_dim`` onto ``target_coords``.

    If ``target_coords`` is a :class:`xr.DataArray` and ``target_name``
    is None, the new dim name is taken from ``target_coords.name``;
    otherwise ``target_name`` (or ``source_dim``) is used.

    Every coordinate that depends on ``source_dim`` is dropped from the
    output; the new target coordinate is attached under the new name.
    The Layer-1 ``RemapAxis`` operator handles Dataset selection /
    per-variable looping.

    Args:
        da: Numeric DataArray carrying ``source_dim``.
        source_dim: Name of the dimension being replaced.
        target_coords: 1-D target axis values.
        target_name: Explicit name for the new dimension.
        method: ``"linear"`` or ``"nearest"``.
        source_coords: Source axis *values*. ``None`` (default) uses the
            1-D dimension coordinate ``da[source_dim]``. A coordinate name
            on ``da`` or a DataArray broadcastable against ``da`` supplies
            per-column levels — it must carry ``source_dim`` and may vary
            along any of the other dims (e.g. ``z_agl(time, level, y, x)``
            on a terrain-following grid). Every column must be strictly
            monotone along ``source_dim``.
        extrapolate: ``"nan"`` (default) — targets outside a column's
            range become ``NaN``; ``"nearest"`` — they hold that column's
            end value.

    Returns:
        ``da`` with ``source_dim`` replaced by the target axis. When
        ``source_coords`` is given, chunked (dask) inputs stay lazy as
        long as ``source_dim`` sits in a single chunk.

    Raises:
        ValueError: If ``source_dim`` is missing, ``source_coords`` does
            not carry it, a column is not monotone, or ``extrapolate`` /
            ``method`` is unknown.
        TypeError: If ``da`` is not numeric.

    Example:
        Remap a model-level field with a 4-D height-above-ground
        coordinate onto fixed heights, holding the lowest / highest model
        level beyond each column's range:

        ```pycon
        >>> import numpy as np, xarray as xr
        >>> from xrops.transforms import remap_axis
        >>> level = np.arange(4)
        >>> z_agl = 20.0 * (level[None, :, None, None] + 1) + np.zeros((2, 4, 3, 3))
        >>> z_agl[..., 0, 0] += 15.0  # one column starts higher up
        >>> t = xr.DataArray(
        ...     300.0 - 0.01 * z_agl,
        ...     dims=("time", "level", "y", "x"),
        ...     coords={"z_agl": (("time", "level", "y", "x"), z_agl)},
        ... )
        >>> out = remap_axis(
        ...     t,
        ...     source_dim="level",
        ...     target_coords=np.array([10.0, 50.0, 90.0]),
        ...     target_name="height",
        ...     source_coords="z_agl",
        ...     extrapolate="nearest",
        ... )
        >>> out.dims, out.sizes["height"]
        (('time', 'height', 'y', 'x'), 3)
        >>> float(out.isel(time=0, height=1, y=1, x=1))  # in range: linear
        299.5

        ```
    """
    if source_dim not in da.dims:
        raise ValueError(
            f"source_dim {source_dim!r} not in DataArray dims {tuple(da.dims)}"
        )
    if extrapolate not in _kernels._EXTRAPOLATE_MODES:
        raise ValueError(
            f"unknown extrapolate {extrapolate!r}; expected one of "
            f"{_kernels._EXTRAPOLATE_MODES!r}"
        )
    if source_coords is None and source_dim not in da.coords:
        raise ValueError(
            f"DataArray must carry a coordinate named {source_dim!r} "
            "for the source axis values"
        )
    if not np.issubdtype(da.dtype, np.number):
        raise TypeError(
            f"remap_axis requires numeric data; got dtype {da.dtype}. "
            "Drop or convert the variable before calling remap_axis."
        )

    if isinstance(target_coords, xr.DataArray):
        new_name = target_name or target_coords.name or source_dim
        tgt = np.asarray(target_coords.values, dtype=float)
    else:
        new_name = target_name or source_dim
        tgt = np.asarray(target_coords, dtype=float)

    if source_coords is not None or extrapolate != "nan":
        return _remap_axis_columns(
            da,
            source_dim=source_dim,
            tgt=tgt,
            new_name=str(new_name),
            method=method,
            source_coords=source_coords,
            extrapolate=extrapolate,
        )

    src = np.asarray(da[source_dim].values, dtype=float)

    axis = da.get_axis_num(source_dim)
    new_values = _kernels.remap_axis(
        da.values,
        axis=axis,
        source_coords=src,
        target_coords=tgt,
        method=method,
    )
    new_dims = tuple(new_name if d == source_dim else d for d in da.dims)
    new_coords = {
        cname: c
        for cname, c in da.coords.items()
        if source_dim not in c.dims and cname != source_dim
    }
    new_coords[new_name] = xr.DataArray(tgt, dims=(new_name,), name=new_name)
    return xr.DataArray(
        new_values,
        dims=new_dims,
        coords=new_coords,
        attrs=dict(da.attrs),
        name=da.name,
    )

to_phase

to_phase(da: DataArray, *, time_dim: str, period: float, n_bins: int, epoch: float = 0.0) -> DataArray

Fold time_dim onto a phase axis by binning + averaging.

The time coordinate must be numeric in the same units as period. Phase is computed as ((t - epoch) / period) mod 1 and binned into n_bins evenly-spaced bins on [0, 1). Output carries a "phase" dim with coordinate values at bin centers.

The Layer-1 ToPhase operator handles Dataset selection / per-variable looping.

Source code in packages/xrtoolz-ops/src/xrops/transforms/_src/coord_remap.py

def to_phase(
    da: xr.DataArray,
    *,
    time_dim: str,
    period: float,
    n_bins: int,
    epoch: float = 0.0,
) -> xr.DataArray:
    """Fold ``time_dim`` onto a phase axis by binning + averaging.

    The time coordinate must be numeric in the same units as ``period``.
    Phase is computed as ``((t - epoch) / period) mod 1`` and binned
    into ``n_bins`` evenly-spaced bins on ``[0, 1)``. Output carries
    a ``"phase"`` dim with coordinate values at bin centers.

    The Layer-1 ``ToPhase`` operator handles Dataset selection /
    per-variable looping.
    """
    if time_dim not in da.dims:
        raise ValueError(
            f"time_dim {time_dim!r} not in DataArray dims {tuple(da.dims)}"
        )
    if time_dim not in da.coords:
        raise ValueError(
            f"DataArray must carry a coordinate named {time_dim!r} to compute phase"
        )
    if period <= 0:
        raise ValueError(f"period must be > 0, got {period}")
    if n_bins < 1:
        raise ValueError(f"n_bins must be >= 1, got {n_bins}")
    if not np.issubdtype(da.dtype, np.number):
        raise TypeError(
            f"to_phase requires numeric data; got dtype {da.dtype}. "
            "Drop or convert the variable before calling to_phase."
        )

    t = np.asarray(da[time_dim].values, dtype=float)
    finite_t = _finite_mask(t)
    phase = np.where(finite_t, ((t - epoch) / period) % 1.0, 0.0)
    edges = np.linspace(0.0, 1.0, n_bins + 1)
    centers = 0.5 * (edges[:-1] + edges[1:])
    bin_idx = np.clip(np.searchsorted(edges, phase, side="right") - 1, 0, n_bins - 1)

    axis = da.get_axis_num(time_dim)
    moved = np.moveaxis(da.values, axis, 0)
    flat = moved.reshape(moved.shape[0], -1)
    # Use a complex accumulator if the data is complex so the imaginary
    # part isn't dropped (P1 review).
    is_complex = np.iscomplexobj(flat)
    acc_dtype = np.complex128 if is_complex else float
    # cast: with dtype resolved at runtime (complex128 or float), the
    # specialised ndarray type numpy's stubs infer would reject the row
    # assignments below.
    sums = cast(np.ndarray, np.zeros((n_bins, flat.shape[1]), dtype=acc_dtype))
    counts = np.zeros((n_bins, flat.shape[1]), dtype=float)
    # A sample is valid only if its time coord is finite AND every data
    # value is finite. Excluding NaN-time rows keeps stale rows out of
    # the phase means (P2 review). ``np.isfinite`` handles integer
    # dtypes (treated as all-finite) where ``np.isnan`` would raise.
    valid_value = (
        np.isfinite(flat)
        if not is_complex
        else (np.isfinite(flat.real) & np.isfinite(flat.imag))
    )
    valid = finite_t[:, None] & valid_value
    for b in range(n_bins):
        m = (bin_idx == b) & finite_t
        if not m.any():
            continue
        sub = flat[m]
        sub_valid = valid[m]
        zero = acc_dtype(0)
        sums[b] = np.where(sub_valid, sub, zero).sum(axis=0)
        counts[b] = sub_valid.sum(axis=0)
    with np.errstate(invalid="ignore"):
        nan_fill = (np.nan + 0j) if is_complex else np.nan
        mean = np.where(counts > 0, sums / np.where(counts > 0, counts, 1.0), nan_fill)
    out_shape = (n_bins, *moved.shape[1:])
    new_values = mean.reshape(out_shape)
    new_values = np.moveaxis(new_values, 0, axis)
    new_dims = tuple("phase" if d == time_dim else d for d in da.dims)
    new_coords = {
        cname: c
        for cname, c in da.coords.items()
        if time_dim not in c.dims and cname != time_dim
    }
    new_coords["phase"] = xr.DataArray(centers, dims=("phase",), name="phase")
    return xr.DataArray(
        new_values,
        dims=new_dims,
        coords=new_coords,
        attrs=dict(da.attrs),
        name=da.name,
    )

Grid carriers

Grid dataclass

Grid(lon: Float[ndarray, 'lon'], lat: Float[ndarray, 'lat'])

Regular lon/lat grid defined by 1-D coordinate arrays.

from_bounds classmethod

from_bounds(lon_bnds: tuple[float, float], lat_bnds: tuple[float, float], resolution: float) -> Grid

Build a grid whose cells have width resolution degrees.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/binning.py

@classmethod
def from_bounds(
    cls,
    lon_bnds: tuple[float, float],
    lat_bnds: tuple[float, float],
    resolution: float,
) -> Grid:
    """Build a grid whose cells have width ``resolution`` degrees."""
    lon = np.arange(lon_bnds[0], lon_bnds[1] + resolution / 2, resolution)
    lat = np.arange(lat_bnds[0], lat_bnds[1] + resolution / 2, resolution)
    return cls(lon=lon, lat=lat)

bin_edges

bin_edges() -> tuple[Float[ndarray, 'lon_edge'], Float[ndarray, 'lat_edge']]

Return (lon_edges, lat_edges) for use by np.histogram2d.

Source code in packages/xrtoolz-ops/src/xrops/interpolate/_src/binning.py

def bin_edges(
    self,
) -> tuple[Float[np.ndarray, "lon_edge"], Float[np.ndarray, "lat_edge"]]:
    """Return ``(lon_edges, lat_edges)`` for use by ``np.histogram2d``."""
    return _cell_edges(self.lon), _cell_edges(self.lat)

SpaceTimeGrid dataclass

SpaceTimeGrid(lon: Float[ndarray, 'lon'], lat: Float[ndarray, 'lat'], time: DatetimeIndex = (lambda: DatetimeIndex([]))())

Lon/lat/time grid. Use :meth:from_bounds for the common case.

Period dataclass

Period(time_min: str, time_max: str, freq: str = '1D')

A pandas-compatible time window.