Wavelet Analysis

Continuous wavelet transforms in 1-D (Torrence–Compo with Liu-rectified power) and 2-D (Morlet), plus scalograms, the cone of influence, dominant- period maps, and significance testing. The CWT of a signal \(x(t)\) at scale \(s\) and translation \(\tau\) is

\[ W_x(s, \tau) = \int x(t)\,\psi^{*}\!\left(\frac{t-\tau}{s}\right)\frac{dt}{\sqrt{s}}. \]

Operators

WaveletPowerSpectrum

WaveletPowerSpectrum(variable: str, scales: Sequence[float] | DataArray, *, dim: tuple[str, str] = ('y', 'x'), x0: float = 50000.0, ntheta: int = 16, k0: float = 1.0, isotropic: bool = True, output_var: str | None = None)

Bases: Operator

2-D Morlet wavelet power spectrum of one variable.

Parameters:

Name	Type	Description	Default
variable	str	Name of the variable to transform.	required
scales	Sequence[float] | DataArray	Wavelet scales (a sequence or DataArray).	required
dim	tuple[str, str]	The two spatial dims (y, x).	('y', 'x')
x0	float	Reference length scale (m).	50000.0
ntheta	int	Number of orientation angles.	16
k0	float	Morlet central wavenumber.	1.0
isotropic	bool	Average over orientation when True.	True
output_var	str | None	Name for the output variable (default "<variable>_wpsd").	None

Returns:

Type	Description
	The input dataset with the wavelet power-spectrum variable added.

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

def __init__(
    self,
    variable: str,
    scales: Sequence[float] | xr.DataArray,
    *,
    dim: tuple[str, str] = ("y", "x"),
    x0: float = 50e3,
    ntheta: int = 16,
    k0: float = 1.0,
    isotropic: bool = True,
    output_var: str | None = None,
) -> None:
    self.variable = variable
    self.scales = scales
    self.dim = tuple(dim)
    self.x0 = float(x0)
    self.ntheta = int(ntheta)
    self.k0 = float(k0)
    self.isotropic = bool(isotropic)
    self.output_var = output_var

WaveletScalogram

WaveletScalogram(variable: str, *, dim: str = 'time', mother: str = 'morlet', param: float | None = None, s0: float | None = None, dj: float = 0.25, j_max: int | None = None, rectify: bool = True, output_prefix: str | None = None)

Bases: Operator

Compute a 1-D wavelet scalogram for one variable.

Adds six derived variables under a common prefix: <prefix>_wave, <prefix>_power, <prefix>_power_rect, <prefix>_scalogram, <prefix>_coi, <prefix>_coi_mask. The prefix defaults to variable and is overridable via output_prefix.

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

def __init__(
    self,
    variable: str,
    *,
    dim: str = "time",
    mother: str = "morlet",
    param: float | None = None,
    s0: float | None = None,
    dj: float = 0.25,
    j_max: int | None = None,
    rectify: bool = True,
    output_prefix: str | None = None,
) -> None:
    self.variable = variable
    self.dim = dim
    self.mother = mother
    self.param = param
    self.s0 = s0
    self.dj = float(dj)
    self.j_max = j_max
    self.rectify = bool(rectify)
    self.output_prefix = output_prefix

WaveletSignificance

WaveletSignificance(variable: str, *, dim_time: str = 'time', dim_scale: str = 'scale', null: Literal['red', 'white'] = 'red', alpha: float | None = None, confidence: float = 0.95, mother: str = 'morlet', param: float | None = None, output_var: str | None = None)

Bases: Operator

Torrence-Compo significance mask for a wavelet-power variable.

Parameters:

Name	Type	Description	Default
variable	str	Name of the wavelet-power variable to test.	required
dim_time	str	Time dimension name.	'time'
dim_scale	str	Scale dimension name.	'scale'
null	Literal['red', 'white']	Null spectrum — "red" (AR(1)) or "white" noise.	'red'
alpha	float | None	Lag-1 autocorrelation; None uses the stored value.	None
confidence	float	Chi-square confidence level.	0.95
mother	str	Mother wavelet name.	'morlet'
param	float | None	Optional mother parameter.	None
output_var	str | None	Name for the mask variable (default "<variable>_signif_mask").	None

Returns:

Type	Description
	The input dataset with a boolean significance-mask variable added.

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

def __init__(
    self,
    variable: str,
    *,
    dim_time: str = "time",
    dim_scale: str = "scale",
    null: Literal["red", "white"] = "red",
    alpha: float | None = None,
    confidence: float = 0.95,
    mother: str = "morlet",
    param: float | None = None,
    output_var: str | None = None,
) -> None:
    self.variable = variable
    self.dim_time = dim_time
    self.dim_scale = dim_scale
    self.null = null
    self.alpha = alpha
    self.confidence = float(confidence)
    self.mother = mother
    self.param = param
    self.output_var = output_var

BandpassWavelength

BandpassWavelength(*, dim: str, lambda_min_km: float | None = None, lambda_max_km: float | None = None, spacing_km: float | None = None, method: str = 'lanczos', num_taps: int | None = None, attenuation_db: float | None = None, lon: str = 'lon', lat: str = 'lat')

Bases: Operator

Band-pass a field by wavelength with a windowed-sinc FIR filter.

Keeps along-dim structure with wavelengths between lambda_min_km and lambda_max_km. When spacing_km is omitted it is derived from the lon/lat ride-along coordinates.

Parameters:

Name	Type	Description	Default
dim	str	Dimension to filter along (e.g. an along-track index).	required
lambda_min_km	float | None	Shortest wavelength to keep (km); None for a low-pass.	None
lambda_max_km	float | None	Longest wavelength to keep (km); None for a high-pass.	None
spacing_km	float | None	Sample spacing (km); None derives it from lon/lat.	None
method	str	FIR window — "lanczos" or "kaiser".	'lanczos'
num_taps	int | None	Filter length; None chooses a default.	None
attenuation_db	float | None	Stop-band attenuation for the Kaiser window.	None
lon	str	Longitude coordinate name (used when spacing_km is None).	'lon'
lat	str	Latitude coordinate name (used when spacing_km is None).	'lat'

Returns:

Type	Description
	The band-pass-filtered dataset (each numeric variable along dim).

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

def __init__(
    self,
    *,
    dim: str,
    lambda_min_km: float | None = None,
    lambda_max_km: float | None = None,
    spacing_km: float | None = None,
    method: str = "lanczos",
    num_taps: int | None = None,
    attenuation_db: float | None = None,
    lon: str = "lon",
    lat: str = "lat",
) -> None:
    self.dim = dim
    self.lambda_min_km = lambda_min_km
    self.lambda_max_km = lambda_max_km
    self.spacing_km = spacing_km
    self.method = method
    self.num_taps = num_taps
    self.attenuation_db = attenuation_db
    self.lon = lon
    self.lat = lat

Functional primitives (Layer 0)

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

cwt1d

cwt1d(da: DataArray, *, dim: str = 'time', mother: str = 'morlet', param: float | None = None, s0: float | None = None, dj: float = 0.25, j_max: int | None = None, pad: bool = True) -> Dataset

Compute a 1-D continuous wavelet transform along one dimension.

The implementation follows Torrence and Compo (1998): zero-pad to the next power of two, multiply the FFT of the signal by the analytic Fourier transform of the selected mother wavelet, and strip the padding after the inverse FFT. The returned power_rect variable is Liu et al. (2007) rectified power, |W|² / scale.

Parameters:

Name	Type	Description	Default
da	DataArray	Input time series or array with one transform dimension.	required
dim	str	Dimension to transform.	'time'
mother	str	Mother wavelet: "morlet", "paul", or "dog".	'morlet'
param	float | None	Mother parameter. Defaults are Morlet k0=6, Paul m=4, and DOG m=2.	None
s0	float | None	Smallest scale. Defaults to 2 * dt.	None
dj	float	Scale spacing in octaves. 0.25 gives four voices per octave.	0.25
j_max	int | None	Largest scale index. Defaults to log2(n * dt / s0) / dj.	None
pad	bool	If True, zero-pad to the next power of two before the FFT.	True

Returns:

Type	Description
Dataset	Dataset with wave, power, power_rect, coi, and
Dataset	coi_mask variables plus scale and period coordinates.

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/wavelet1d.py

def cwt1d(
    da: xr.DataArray,
    *,
    dim: str = "time",
    mother: str = "morlet",
    param: float | None = None,
    s0: float | None = None,
    dj: float = 0.25,
    j_max: int | None = None,
    pad: bool = True,
) -> xr.Dataset:
    """Compute a 1-D continuous wavelet transform along one dimension.

    The implementation follows Torrence and Compo (1998): zero-pad to the
    next power of two, multiply the FFT of the signal by the analytic Fourier
    transform of the selected mother wavelet, and strip the padding after the
    inverse FFT. The returned ``power_rect`` variable is Liu et al. (2007)
    rectified power, ``|W|² / scale``.

    Args:
        da: Input time series or array with one transform dimension.
        dim: Dimension to transform.
        mother: Mother wavelet: ``"morlet"``, ``"paul"``, or ``"dog"``.
        param: Mother parameter. Defaults are Morlet ``k0=6``, Paul ``m=4``,
            and DOG ``m=2``.
        s0: Smallest scale. Defaults to ``2 * dt``.
        dj: Scale spacing in octaves. ``0.25`` gives four voices per octave.
        j_max: Largest scale index. Defaults to ``log2(n * dt / s0) / dj``.
        pad: If ``True``, zero-pad to the next power of two before the FFT.

    Returns:
        Dataset with ``wave``, ``power``, ``power_rect``, ``coi``, and
        ``coi_mask`` variables plus ``scale`` and ``period`` coordinates.
    """
    _require_dim(da, dim)
    _raise_if_chunked_along_dim(da, dim, "cwt1d")
    dt = coord_spacing(da, dim, require_coord=False)
    scale = _scale_grid(da.sizes[dim], dt=dt, s0=s0, dj=dj, j_max=j_max)
    mother_name = _normalize_mother(mother)
    param_value = _default_param(mother_name, param)
    fourier_factor, coi_factor, _ = _wavelet_factors(mother_name, param_value)
    period = scale * fourier_factor

    wave = xr.apply_ufunc(
        _cwt1d_numpy,
        da,
        input_core_dims=[[dim]],
        output_core_dims=[["scale", dim]],
        kwargs={
            "dt": dt,
            "scale": scale,
            "mother": mother_name,
            "param": param_value,
            "pad": pad,
        },
        vectorize=True,
        dask="parallelized",
        output_dtypes=[np.complex128],
        dask_gufunc_kwargs={"output_sizes": {"scale": scale.size}},
    )
    outer_dims = tuple(d for d in da.dims if d != dim)
    wave = wave.transpose("scale", dim, *outer_dims)
    coords = {
        "scale": ("scale", scale),
        "period": ("scale", period),
        dim: da[dim],
        "fourier_factor": fourier_factor,
    }
    wave = wave.assign_coords(coords)
    wave.name = "wave"
    wave.attrs.update(
        {
            "mother": mother_name,
            "param": param_value,
            "dt": dt,
            "dj": dj,
            "lag1_autocorrelation": _lag1_autocorrelation(da, dim),
        }
    )
    # Per-series variance preserves outer-dim broadcasting for significance
    # thresholds — a scalar `source_variance` would over/under-flag stacks of
    # series with different amplitudes.
    source_variance = da.var(dim, skipna=True)
    wave = wave.assign_coords(source_variance=source_variance)
    power = (xr.apply_ufunc(np.abs, wave) ** 2).rename("power")
    power.attrs.update(wave.attrs)
    power_rect = (
        power / xr.DataArray(scale, dims=("scale",), coords={"scale": scale})
    ).rename("power_rect")
    power_rect.attrs.update(power.attrs)
    power_rect.attrs["bias_corrected"] = "liu_2007"
    coi = _coi(da[dim], dt=dt, fourier_factor=fourier_factor, coi_factor=coi_factor)
    coi_mask = _coi_mask(scale=scale, period=period, coi=coi, dim=dim)
    return xr.Dataset(
        {
            "wave": wave,
            "power": power,
            "power_rect": power_rect,
            "coi": coi,
            "coi_mask": coi_mask,
        }
    )

icwt1d

icwt1d(wave: DataArray, *, band: tuple[float, float] | None = None, mother: str = 'morlet', param: float | None = None, dj: float | None = None, dt: float | None = None) -> DataArray

Reconstruct a signal from 1-D CWT coefficients.

Parameters:

Name	Type	Description	Default
wave	DataArray	Complex coefficients from :func:cwt1d.	required
band	tuple[float, float] | None	Optional (min, max) band to reconstruct. Interpreted in period units when wave carries a period coordinate (the usual case after :func:cwt1d); otherwise interpreted in scale units.	None
mother	str	Mother wavelet name.	'morlet'
param	float | None	Optional mother parameter. Reconstruction constants are available for the Torrence-Compo defaults.	None
dj	float | None	Scale spacing used by :func:cwt1d. Defaults to the value stored in wave.attrs["dj"] (falling back to 0.25 if absent). Reconstruction is only correct when dj matches the value used during decomposition.	None
dt	float | None	Sample spacing. Defaults to the coefficient dt attribute.	None

Returns:

Type	Description
DataArray	Reconstructed DataArray without the scale dimension.

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/wavelet1d.py

def icwt1d(
    wave: xr.DataArray,
    *,
    band: tuple[float, float] | None = None,
    mother: str = "morlet",
    param: float | None = None,
    dj: float | None = None,
    dt: float | None = None,
) -> xr.DataArray:
    """Reconstruct a signal from 1-D CWT coefficients.

    Args:
        wave: Complex coefficients from :func:`cwt1d`.
        band: Optional ``(min, max)`` band to reconstruct. Interpreted in
            **period** units when ``wave`` carries a ``period`` coordinate
            (the usual case after :func:`cwt1d`); otherwise interpreted in
            **scale** units.
        mother: Mother wavelet name.
        param: Optional mother parameter. Reconstruction constants are
            available for the Torrence-Compo defaults.
        dj: Scale spacing used by :func:`cwt1d`. Defaults to the value
            stored in ``wave.attrs["dj"]`` (falling back to ``0.25`` if
            absent). Reconstruction is only correct when ``dj`` matches
            the value used during decomposition.
        dt: Sample spacing. Defaults to the coefficient ``dt`` attribute.

    Returns:
        Reconstructed DataArray without the scale dimension.
    """
    _require_dim(wave, "scale")
    mother_name = _normalize_mother(mother)
    param_value = _default_param(mother_name, param)
    _, _, cdelta = _wavelet_factors(mother_name, param_value)
    psi0 = _psi0(mother_name, param_value)
    dt_value = float(wave.attrs.get("dt", 1.0) if dt is None else dt)
    dj_value = float(wave.attrs.get("dj", 0.25) if dj is None else dj)
    coeffs = wave
    if band is not None:
        band_coord = coeffs["period"] if "period" in coeffs.coords else coeffs["scale"]
        lo, hi = sorted((float(band[0]), float(band[1])))
        coeffs = coeffs.where((band_coord >= lo) & (band_coord <= hi), other=0.0)
    scale = xr.DataArray(
        np.asarray(coeffs["scale"].values, dtype=float),
        dims=("scale",),
        coords={"scale": coeffs["scale"]},
    )
    reconstructed = (coeffs.real / np.sqrt(scale)).sum("scale") * (
        dj_value * np.sqrt(dt_value) / (cdelta * psi0)
    )
    reconstructed.name = (
        _strip_suffix(str(wave.name) if wave.name is not None else None, "_wave")
        or "icwt"
    )
    return reconstructed

cwt2

cwt2(da: DataArray, s: DataArray, dim: tuple[str, str] = ('y', 'x'), *, x0: float = 50000.0, ntheta: int = 16, k0: float = 1.0, a: float = 1.0) -> DataArray

Compute a 2-D directional Morlet CWT on a Cartesian grid.

Parameters:

Name	Type	Description	Default
da	DataArray	Two-dimensional field.	required
s	DataArray	One-dimensional positive scale coordinate.	required
dim	tuple[str, str]	Spatial dimensions as (y, x).	('y', 'x')
x0	float	Reference length scale in the same units as dim coords.	50000.0
ntheta	int	Number of evenly spaced angles in [0, 2π).	16
k0	float	Dimensionless Morlet central wavenumber.	1.0
a	float	Fourier-domain amplitude multiplier.	1.0

Returns:

Type	Description
DataArray	Complex coefficients with dims ("scale", "angle", y, x) and a
DataArray	boolean coi_mask coordinate with dims ("scale", y, x).

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/wavelet.py

def cwt2(
    da: xr.DataArray,
    s: xr.DataArray,
    dim: tuple[str, str] = ("y", "x"),
    *,
    x0: float = 50e3,
    ntheta: int = 16,
    k0: float = 1.0,
    a: float = 1.0,
) -> xr.DataArray:
    """Compute a 2-D directional Morlet CWT on a Cartesian grid.

    Args:
        da: Two-dimensional field.
        s: One-dimensional positive scale coordinate.
        dim: Spatial dimensions as ``(y, x)``.
        x0: Reference length scale in the same units as ``dim`` coords.
        ntheta: Number of evenly spaced angles in ``[0, 2π)``.
        k0: Dimensionless Morlet central wavenumber.
        a: Fourier-domain amplitude multiplier.

    Returns:
        Complex coefficients with dims ``("scale", "angle", y, x)`` and a
        boolean ``coi_mask`` coordinate with dims ``("scale", y, x)``.
    """
    ydim, xdim = dim
    s = _as_scale_dataarray(s)
    _validate_cwt_input(da, s, dim=dim, ntheta=ntheta)
    dy = coord_spacing(da, ydim)
    dx = coord_spacing(da, xdim)
    arr = np.asarray(da.transpose(ydim, xdim).values)
    nan_mask = np.isnan(arr)
    if nan_mask.any():
        fill = float(np.nanmean(arr)) if not np.all(nan_mask) else 0.0
        arr = np.where(nan_mask, fill, arr)

    ky = fftfreq(arr.shape[0], d=dy)
    kx = fftfreq(arr.shape[1], d=dx)
    kx_grid, ky_grid = np.meshgrid(kx, ky)
    field_ft = fft2(arr)

    scale_values = _scale_values(s)
    angles = np.linspace(0.0, 2.0 * np.pi, int(ntheta), endpoint=False)
    coeffs = np.empty((scale_values.size, angles.size, *arr.shape), dtype=np.complex128)
    for i, scale in enumerate(scale_values):
        for j, angle in enumerate(angles):
            kernel = morlet2_ft(
                kx_grid,
                ky_grid,
                float(scale),
                float(angle),
                x0=x0,
                k0=k0,
                a=a,
            )
            coeffs[i, j] = ifft2(field_ft * np.conj(kernel))

    scale_dim = s.dims[0]
    out = xr.DataArray(
        coeffs,
        dims=(scale_dim, "angle", ydim, xdim),
        coords={
            scale_dim: s,
            "angle": angles,
            ydim: da[ydim],
            xdim: da[xdim],
        },
        name=f"{da.name or 'field'}_cwt2",
        attrs=dict(da.attrs),
    )
    out = out.assign_coords(
        coi_mask=build_coi_mask(da, s, dim=dim, x0=x0),
        wavenumber=scale_to_wavenumber(s, x0=x0, k0=k0),
    )
    out["wavenumber"].attrs["units"] = f"cycles per {xdim}/{ydim} coordinate unit"
    return out

morlet2_ft

morlet2_ft(kx: Float[ndarray, 'y x'], ky: Float[ndarray, 'y x'], s: float, theta: float, *, x0: float, k0: float = 1.0, a: float = 1.0) -> Complex[ndarray, 'y x']

Analytic Fourier transform of a rotated 2-D Morlet wavelet.

Parameters:

Name	Type	Description	Default
kx	Float[ndarray, 'y x']	Zonal wavenumber grid in cycles per coordinate unit.	required
ky	Float[ndarray, 'y x']	Meridional wavenumber grid in cycles per coordinate unit.	required
s	float	Positive dimensionless scale.	required
theta	float	Wavelet orientation in radians.	required
x0	float	Reference length scale in coordinate units.	required
k0	float	Dimensionless central wavenumber.	1.0
a	float	Amplitude multiplier.	1.0

Returns:

Type	Description
Complex[ndarray, 'y x']	Complex Fourier-domain wavelet sampled on kx/ky.

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/wavelet.py

def morlet2_ft(
    kx: Float[np.ndarray, "y x"],
    ky: Float[np.ndarray, "y x"],
    s: float,
    theta: float,
    *,
    x0: float,
    k0: float = 1.0,
    a: float = 1.0,
) -> Complex[np.ndarray, "y x"]:
    """Analytic Fourier transform of a rotated 2-D Morlet wavelet.

    Args:
        kx: Zonal wavenumber grid in cycles per coordinate unit.
        ky: Meridional wavenumber grid in cycles per coordinate unit.
        s: Positive dimensionless scale.
        theta: Wavelet orientation in radians.
        x0: Reference length scale in coordinate units.
        k0: Dimensionless central wavenumber.
        a: Amplitude multiplier.

    Returns:
        Complex Fourier-domain wavelet sampled on ``kx``/``ky``.
    """
    if s <= 0:
        raise ValueError("s must be strictly positive")
    if x0 <= 0:
        raise ValueError("x0 must be strictly positive")
    if k0 <= 0:
        raise ValueError("k0 must be strictly positive")
    sigma = s * x0
    kc = k0 / sigma
    kcx = kc * np.cos(theta)
    kcy = kc * np.sin(theta)
    envelope = np.exp(-2.0 * np.pi**2 * sigma**2 * ((kx - kcx) ** 2 + (ky - kcy) ** 2))
    return (a * 2.0 * np.pi * sigma**2 * envelope).astype(np.complex128)

build_coi_mask

build_coi_mask(da: DataArray, scales: DataArray, *, dim: tuple[str, str] = ('y', 'x'), x0: float) -> DataArray

Build a cone-of-influence mask for boundaries and NaN land cells.

A sample is trustworthy when its distance to the nearest boundary or invalid cell exceeds scale * x0 in coordinate units.

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/wavelet_utils.py

def build_coi_mask(
    da: xr.DataArray,
    scales: xr.DataArray,
    *,
    dim: tuple[str, str] = ("y", "x"),
    x0: float,
) -> xr.DataArray:
    """Build a cone-of-influence mask for boundaries and NaN land cells.

    A sample is trustworthy when its distance to the nearest boundary or
    invalid cell exceeds ``scale * x0`` in coordinate units.
    """
    ydim, xdim = dim
    _require_dims(da, *dim)
    dy = coord_spacing(da, ydim)
    dx = coord_spacing(da, xdim)
    arr = np.asarray(da.transpose(ydim, xdim).values, dtype=float)
    valid = np.isfinite(arr)

    yy, xx = np.indices(valid.shape, dtype=float)
    boundary_distance = np.minimum.reduce(
        [
            yy * dy,
            (valid.shape[0] - 1 - yy) * dy,
            xx * dx,
            (valid.shape[1] - 1 - xx) * dx,
        ]
    )

    if valid.all():
        invalid_distance = np.full(valid.shape, np.inf, dtype=float)
    else:
        invalid_distance = distance_transform_edt(valid, sampling=(dy, dx))
    distance = np.minimum(boundary_distance, invalid_distance)

    scale_values = _scale_values(scales)
    mask = valid[None, :, :] & (distance[None, :, :] > scale_values[:, None, None] * x0)
    scale_dim = scales.dims[0]
    return xr.DataArray(
        mask,
        dims=(scale_dim, ydim, xdim),
        coords={
            scale_dim: scales,
            ydim: da[ydim],
            xdim: da[xdim],
        },
        name="coi_mask",
    )

dominant_period_map

dominant_period_map(power_rect: DataArray, *, dim_time: str = 'time', dim_scale: str = 'scale', coi_mask: DataArray | None = None, signif_mask: DataArray | None = None, reduce: DominantReduce = 'argmax') -> DataArray

Return the dominant Fourier period after time-averaging rectified power.

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/wavelet1d.py

def dominant_period_map(
    power_rect: xr.DataArray,
    *,
    dim_time: str = "time",
    dim_scale: str = "scale",
    coi_mask: xr.DataArray | None = None,
    signif_mask: xr.DataArray | None = None,
    reduce: DominantReduce = "argmax",
) -> xr.DataArray:
    """Return the dominant Fourier period after time-averaging rectified power."""
    _require_dim(power_rect, dim_time)
    _require_dim(power_rect, dim_scale)
    if reduce not in {"argmax", "median_argmax"}:
        raise ValueError("reduce must be 'argmax' or 'median_argmax'")
    masked = power_rect
    if coi_mask is not None:
        masked = masked.where(coi_mask)
    if signif_mask is not None:
        masked = masked.where(signif_mask)
    if reduce == "median_argmax":
        summary = masked.median(dim_time, skipna=True)
    else:
        summary = masked.mean(dim_time, skipna=True)
    peak_scale = summary.idxmax(dim_scale)
    if "period" in power_rect.coords and power_rect["period"].dims == (dim_scale,):
        period_lookup = xr.DataArray(
            np.asarray(power_rect["period"].values, dtype=float),
            dims=(dim_scale,),
            coords={dim_scale: power_rect[dim_scale]},
        )
        out = period_lookup.sel({dim_scale: peak_scale})
    else:
        out = peak_scale
    out = out.rename("dominant_period")
    out.attrs["long_name"] = "dominant Fourier period"
    return out

wvlt_power_spectrum

wvlt_power_spectrum(da: DataArray, s: DataArray, *, dim: tuple[str, str] = ('y', 'x'), x0: float = 50000.0, ntheta: int = 16, k0: float = 1.0, isotropic: bool = True) -> DataArray

Return the 2-D Morlet wavelet power spectrum of da.

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/wavelet.py

def wvlt_power_spectrum(
    da: xr.DataArray,
    s: xr.DataArray,
    *,
    dim: tuple[str, str] = ("y", "x"),
    x0: float = 50e3,
    ntheta: int = 16,
    k0: float = 1.0,
    isotropic: bool = True,
) -> xr.DataArray:
    """Return the 2-D Morlet wavelet power spectrum of ``da``."""
    s = _as_scale_dataarray(s)
    coeffs = cwt2(da, s, dim=dim, x0=x0, ntheta=ntheta, k0=k0)
    power = xr.apply_ufunc(np.abs, coeffs) ** 2
    power = _normalize_power(power, da)
    power.name = f"{da.name or 'field'}_wpsd"
    power.attrs.update({"long_name": "2-D Morlet wavelet power spectrum"})
    if isotropic:
        power = power.mean("angle", keep_attrs=True)
    return power

wvlt_cross_spectrum

wvlt_cross_spectrum(da1: DataArray, da2: DataArray, s: DataArray, *, dim: tuple[str, str] = ('y', 'x'), x0: float = 50000.0, ntheta: int = 16, k0: float = 1.0, isotropic: bool = True) -> DataArray

Return the complex 2-D Morlet cross-spectrum of two fields.

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/wavelet.py

def wvlt_cross_spectrum(
    da1: xr.DataArray,
    da2: xr.DataArray,
    s: xr.DataArray,
    *,
    dim: tuple[str, str] = ("y", "x"),
    x0: float = 50e3,
    ntheta: int = 16,
    k0: float = 1.0,
    isotropic: bool = True,
) -> xr.DataArray:
    """Return the complex 2-D Morlet cross-spectrum of two fields."""
    s = _as_scale_dataarray(s)
    w1 = cwt2(da1, s, dim=dim, x0=x0, ntheta=ntheta, k0=k0)
    w2 = cwt2(da2, s, dim=dim, x0=x0, ntheta=ntheta, k0=k0)
    cross = w1 * np.conj(w2)
    if isotropic:
        cross = cross.mean("angle", keep_attrs=True)
    cross.name = f"{da1.name or 'field'}_{da2.name or 'field'}_wcsd"
    return cross

wavelet_significance

wavelet_significance(power: DataArray, *, dim_time: str = 'time', dim_scale: str = 'scale', null: NullModel = 'red', alpha: float | None = None, confidence: float = 0.95, mother: str = 'morlet', param: float | None = None) -> DataArray

Return a mask where wavelet power exceeds a red/white-noise threshold.

Parameters:

Name	Type	Description	Default
power	DataArray	Wavelet power with scale and time dimensions.	required
dim_time	str	Time dimension name.	'time'
dim_scale	str	Scale dimension name.	'scale'
null	NullModel	Null spectrum, either AR(1) red noise or white noise.	'red'
alpha	float | None	Lag-1 autocorrelation. When omitted, the value saved by :func:cwt1d is used, falling back to white noise.	None
confidence	float	Chi-square confidence level.	0.95
mother	str	Mother wavelet name, used for the scale-period mapping if no period coordinate is attached.	'morlet'
param	float | None	Optional mother parameter.	None

Returns:

Type	Description
DataArray	Boolean DataArray with the same dimensions as power.

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/wavelet1d.py

def wavelet_significance(
    power: xr.DataArray,
    *,
    dim_time: str = "time",
    dim_scale: str = "scale",
    null: NullModel = "red",
    alpha: float | None = None,
    confidence: float = 0.95,
    mother: str = "morlet",
    param: float | None = None,
) -> xr.DataArray:
    """Return a mask where wavelet power exceeds a red/white-noise threshold.

    Args:
        power: Wavelet power with scale and time dimensions.
        dim_time: Time dimension name.
        dim_scale: Scale dimension name.
        null: Null spectrum, either AR(1) red noise or white noise.
        alpha: Lag-1 autocorrelation. When omitted, the value saved by
            :func:`cwt1d` is used, falling back to white noise.
        confidence: Chi-square confidence level.
        mother: Mother wavelet name, used for the scale-period mapping if no
            ``period`` coordinate is attached.
        param: Optional mother parameter.

    Returns:
        Boolean DataArray with the same dimensions as ``power``.
    """
    _require_dim(power, dim_time)
    _require_dim(power, dim_scale)
    if not 0.0 < confidence < 1.0:
        raise ValueError("confidence must be between 0 and 1")
    null_name = null.lower()
    if null_name not in {"red", "white"}:
        raise ValueError("null must be either 'red' or 'white'")
    null_model = cast(NullModel, null_name)
    alpha_value = 0.0 if alpha is None else float(alpha)
    if alpha is None:
        alpha_value = float(power.attrs.get("lag1_autocorrelation", 0.0))
    if not -1.0 < alpha_value < 1.0:
        raise ValueError("alpha must be between -1 and 1")

    scale = np.asarray(power[dim_scale].values, dtype=float)
    if "period" in power.coords and power["period"].dims == (dim_scale,):
        period = np.asarray(power["period"].values, dtype=float)
    else:
        mother_name = _normalize_mother(mother)
        param_value = _default_param(mother_name, param)
        fourier_factor, _, _ = _wavelet_factors(mother_name, param_value)
        period = scale * fourier_factor
    # source_variance is per-outer-index when written by cwt1d so the
    # threshold broadcasts correctly across stacks of series. Fall back to
    # the per-series wavelet-power mean if the coord is missing (e.g. when
    # power was constructed by hand).
    if "source_variance" in power.coords:
        variance: xr.DataArray | float = power.coords["source_variance"].astype(float)
    elif "source_variance" in power.attrs:
        variance = float(power.attrs["source_variance"])
    else:
        variance = power.mean(dim_time, skipna=True).mean(dim_scale, skipna=True)
    if null_model == "white":
        background = np.ones_like(period)
    else:
        dt = float(power.attrs.get("dt", 1.0))
        background = (1.0 - alpha_value**2) / (
            1.0 + alpha_value**2 - 2.0 * alpha_value * np.cos(2.0 * np.pi * dt / period)
        )
    background_da = xr.DataArray(
        background,
        dims=(dim_scale,),
        coords={dim_scale: power[dim_scale]},
    )
    threshold = variance * background_da * _chi2_ppf(confidence, df=2.0) / 2.0
    if power.attrs.get("bias_corrected") == "liu_2007":
        threshold = threshold / xr.DataArray(
            scale, dims=(dim_scale,), coords={dim_scale: power[dim_scale]}
        )
    out = (power > threshold).rename("signif_mask")
    out.attrs.update(
        {
            "null": null_model,
            "alpha": alpha_value,
            "confidence": confidence,
        }
    )
    return out

bandpass_wavelength

bandpass_wavelength(da: DataArray, *, dim: str, lambda_min_km: float | None = None, lambda_max_km: float | None = None, spacing_km: float | None = None, method: str = 'lanczos', num_taps: int | None = None, attenuation_db: float | None = None, lon: DataArray | None = None, lat: DataArray | None = None) -> DataArray

Filter an along-track DataArray using wavelength cutoffs in kilometres.

Per the PR γ primitive-flip (docs/design/xarray-native-primitives.md), this Layer-0 primitive operates on a single DataArray. The Layer-1 BandpassWavelength operator handles Dataset selection (including per-variable iteration and string lon/lat name lookup) and forwards the resulting DataArrays here.

Parameters:

Name	Type	Description	Default
da	DataArray	Input DataArray. Must carry dim and a numeric dtype.	required
dim	str	Along-track dimension to filter.	required
lambda_min_km	float | None	Short-wavelength edge in kilometres. With lambda_max_km this becomes the high-frequency band edge.	None
lambda_max_km	float | None	Long-wavelength edge in kilometres. With lambda_min_km this becomes the low-frequency band edge.	None
spacing_km	float | None	Optional along-track sample spacing. If omitted, it is derived as the median WGS-84 spacing from lon / lat.	None
method	str	FIR window family, "lanczos" or "kaiser".	'lanczos'
num_taps	int | None	Optional odd FIR tap count.	None
attenuation_db	float | None	Kaiser attenuation target in decibels.	None
lon	DataArray | None	1-D longitude DataArray along dim, used to derive spacing_km when it is not supplied explicitly.	None
lat	DataArray | None	1-D latitude DataArray along dim, used to derive spacing_km when it is not supplied explicitly.	None

Returns:

Type	Description
DataArray	DataArray filtered along dim, same shape as da.

Source code in packages/xrtoolz-ops/src/xrops/wavelets/_src/along_track.py

def bandpass_wavelength(
    da: xr.DataArray,
    *,
    dim: str,
    lambda_min_km: float | None = None,
    lambda_max_km: float | None = None,
    spacing_km: float | None = None,
    method: str = "lanczos",
    num_taps: int | None = None,
    attenuation_db: float | None = None,
    lon: xr.DataArray | None = None,
    lat: xr.DataArray | None = None,
) -> xr.DataArray:
    """Filter an along-track DataArray using wavelength cutoffs in kilometres.

    Per the PR γ primitive-flip (``docs/design/xarray-native-primitives.md``),
    this Layer-0 primitive operates on a single ``DataArray``. The
    Layer-1 ``BandpassWavelength`` operator handles Dataset selection
    (including per-variable iteration and string ``lon``/``lat`` name
    lookup) and forwards the resulting DataArrays here.

    Args:
        da: Input DataArray. Must carry ``dim`` and a numeric dtype.
        dim: Along-track dimension to filter.
        lambda_min_km: Short-wavelength edge in kilometres. With
            ``lambda_max_km`` this becomes the high-frequency band edge.
        lambda_max_km: Long-wavelength edge in kilometres. With
            ``lambda_min_km`` this becomes the low-frequency band edge.
        spacing_km: Optional along-track sample spacing. If omitted, it is
            derived as the median WGS-84 spacing from ``lon`` / ``lat``.
        method: FIR window family, ``"lanczos"`` or ``"kaiser"``.
        num_taps: Optional odd FIR tap count.
        attenuation_db: Kaiser attenuation target in decibels.
        lon: 1-D longitude DataArray along ``dim``, used to derive
            ``spacing_km`` when it is not supplied explicitly.
        lat: 1-D latitude DataArray along ``dim``, used to derive
            ``spacing_km`` when it is not supplied explicitly.

    Returns:
        DataArray filtered along ``dim``, same shape as ``da``.
    """
    if lambda_min_km is None and lambda_max_km is None:
        raise ValueError("at least one of lambda_min_km or lambda_max_km is required")
    if (
        lambda_min_km is not None
        and lambda_max_km is not None
        and lambda_min_km >= lambda_max_km
    ):
        raise ValueError(
            f"lambda_min_km must be < lambda_max_km; got {lambda_min_km} >= "
            f"{lambda_max_km}"
        )
    if dim not in da.dims:
        raise ValueError(f"dim {dim!r} not in DataArray dims {tuple(da.dims)}")
    if spacing_km is None:
        if lon is None or lat is None:
            raise ValueError(
                "spacing_km must be provided when lon / lat DataArrays are not "
                "supplied for spacing inference"
            )
        # Restrict spacing inference to the filter dim only; otherwise a
        # 2-D coord (multi-track grid) would mix unrelated endpoints into
        # the median.
        if lon.dims != (dim,) or lat.dims != (dim,):
            raise ValueError(
                "Cannot infer spacing: lon/lat are not 1-D along "
                f"{dim!r} (got dims {lon.dims} / {lat.dims}). Pass "
                "spacing_km explicitly."
            )
        spacing_km = median_dx_km(lon.values, lat.values)
    if spacing_km <= 0:
        raise ValueError(f"spacing_km must be > 0, got {spacing_km}")

    if lambda_min_km is None:
        # Only an upper wavelength bound is set — keep wavelengths shorter
        # than lambda_max_km, i.e. high-pass with cutoff at the long edge.
        btype = "high"
        cutoff: float | tuple[float, float] = _cutoff_from_wavelength(
            lambda_max_km, spacing_km
        )
    elif lambda_max_km is None:
        # Only a lower wavelength bound is set — keep wavelengths longer
        # than lambda_min_km, i.e. low-pass with cutoff at the short edge.
        btype = "low"
        cutoff = _cutoff_from_wavelength(lambda_min_km, spacing_km)
    else:
        btype = "bandpass"
        cutoff = (
            _cutoff_from_wavelength(lambda_max_km, spacing_km),
            _cutoff_from_wavelength(lambda_min_km, spacing_km),
        )

    return fir_filter(
        da,
        dim=dim,
        cutoff=cutoff,
        method=method,
        btype=btype,
        num_taps=num_taps,
        attenuation_db=attenuation_db,
    )