Atmosphere

Atmospheric physics on CF datasets. The design rule still holds: only true physics lives here; anything domain-agnostic belongs in cf. Trace-gas (methane) physics sits under xrops.atm.gas.ch4.

WRF

open_wrfout turns a WRF-ARW wrfout file into a CF Dataset on (time, level, y, x): the Times char rows (or XTIME minutes) become a datetime64[s] axis, U / V / W are destaggered to mass points (U / V rotated from grid- to earth-relative when the file carries COSALPHA / SINALPHA), and the split base + perturbation fields are recombined,

\[ p = P + PB,\qquad T = (T' + 300\,\mathrm{K})\left(\frac{p}{10^5\,\mathrm{Pa}}\right)^{R_d/c_p},\qquad z_\text{agl} = \tfrac12\left(z_\text{stag}[k] + z_\text{stag}[k+1]\right) - HGT, \quad z_\text{stag} = \frac{PH + PHB}{g}. \]

Everything stays lazy, so chunks={"time": 1} keeps the file dask backed.

atm

Atmospheric physics operators.

Layer-0 primitives are implemented in :mod:xrops.atm._src and re-exported here.

Content:

• open_wrfout: WRF-ARW wrfout → CF Dataset on (time, level, y, x) (time decode, destagger, physical temperature / pressure, height above ground), with the destagger / wrf_time / wrf_wind / wrf_height_agl / wrf_temperature helpers it is built from.
• Wind diagnostics: wind_speed, wind_direction (meteorological "from" bearing by default) and their inverse wind_components.
• Vertical-column diagnostics: column_integral (trapezoid or cell-width sum), hypsometric_height (pressure levels → geopotential height) and pbl_height_bulk_richardson.
• Layer-1 Operator wrappers (:mod:xrops.atm.operators), all re-exported here.

Trace-gas physics lives under :mod:xrops.atm.gas. Potential temperature is still planned (see docs/design/).

open_wrfout

open_wrfout(path: str | Path, *, times: slice | None = None, variables: str | Sequence[str] | None = None, chunks: Mapping[str, int] | None = None) -> Dataset

Open a wrfout file as a CF-conformant Dataset.

Returns dims (time, level, y, x) with:

• coords: time (datetime64[s]), x / y (metres from the south-west mass point, from DX / DY), level (int index), z_agl(time, level, y, x) (m above ground), and 2-D lon / lat from XLONG / XLAT when present.
• data_vars (all on mass points): u, v (m s⁻¹, destaggered and rotated to earth-relative when COSALPHA / SINALPHA are in the file — see :func:wrf_wind), w (m s⁻¹, destaggered), temperature (K), pressure (Pa), pbl_height (m, from PBLH), terrain (m, from HGT), plus any raw variable named in variables (dims renamed, otherwise untouched).
• attrs: DX, DY, MAP_PROJ, CEN_LAT, CEN_LON, TRUELAT1 / 2, STAND_LON copied through for CRS work.

HGT, XLONG and XLAT lose their time axis only when they are the same at every step; on a moving nest they keep it. CF standard_name / units are stamped via :func:xrproducts.apply_cf_attrs. Nothing is materialised, so chunks={"time": 1} keeps every variable dask backed. The returned Dataset owns the file handle: closing it (or leaving a with block) closes the file.

Parameters:

Name	Type	Description	Default
path	str | Path	wrfout NetCDF file.	required
times	slice | None	Optional slice over the file's time axis.	None
variables	str | Sequence[str] | None	Raw WRF variable name(s) to pass through (e.g. "QVAPOR" or ["QVAPOR", "QCLOUD"]).	None
chunks	Mapping[str, int] | None	Dask chunks keyed by the output dimension names.	None

Returns:

Type	Description
Dataset	The CF Dataset.

Raises:

Type	Description
ValueError	If no time axis can be built (see :func:wrf_time).
KeyError	If a name in variables is not in the file.

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/wrf.py

def open_wrfout(
    path: str | Path,
    *,
    times: slice | None = None,
    variables: str | Sequence[str] | None = None,
    chunks: Mapping[str, int] | None = None,
) -> xr.Dataset:
    """Open a ``wrfout`` file as a CF-conformant Dataset.

    Returns dims ``(time, level, y, x)`` with:

    - coords: ``time`` (``datetime64[s]``), ``x`` / ``y`` (metres from the
      south-west mass point, from ``DX`` / ``DY``), ``level`` (int index),
      ``z_agl(time, level, y, x)`` (m above ground), and 2-D ``lon`` /
      ``lat`` from ``XLONG`` / ``XLAT`` when present.
    - data_vars (all on mass points): ``u``, ``v`` (m s⁻¹, destaggered and
      rotated to earth-relative when ``COSALPHA`` / ``SINALPHA`` are in the
      file — see :func:`wrf_wind`), ``w`` (m s⁻¹, destaggered),
      ``temperature`` (K), ``pressure`` (Pa), ``pbl_height`` (m, from
      ``PBLH``), ``terrain`` (m, from ``HGT``), plus any raw variable named
      in ``variables`` (dims renamed, otherwise untouched).
    - attrs: ``DX``, ``DY``, ``MAP_PROJ``, ``CEN_LAT``, ``CEN_LON``,
      ``TRUELAT1`` / ``2``, ``STAND_LON`` copied through for CRS work.

    ``HGT``, ``XLONG`` and ``XLAT`` lose their ``time`` axis only when they
    are the same at every step; on a moving nest they keep it. CF
    ``standard_name`` / ``units`` are stamped via
    :func:`xrproducts.apply_cf_attrs`. Nothing is materialised, so
    ``chunks={"time": 1}`` keeps every variable dask backed. The returned
    Dataset owns the file handle: closing it (or leaving a ``with`` block)
    closes the file.

    Args:
        path: ``wrfout`` NetCDF file.
        times: Optional ``slice`` over the file's time axis.
        variables: Raw WRF variable name(s) to pass through (e.g.
            ``"QVAPOR"`` or ``["QVAPOR", "QCLOUD"]``).
        chunks: Dask chunks keyed by the *output* dimension names.

    Returns:
        The CF Dataset.

    Raises:
        ValueError: If no time axis can be built (see :func:`wrf_time`).
        KeyError: If a name in ``variables`` is not in the file.
    """
    if isinstance(variables, str):
        variables = [variables]
    raw = xr.open_dataset(path, decode_times=False)
    close = raw._close  # rename/isel/chunk drop it; re-attached on the output
    raw = raw.rename({k: v for k, v in _DIMS.items() if k in raw.dims})
    if times is not None:
        raw = raw.isel(time=times)
    if chunks:
        raw = raw.chunk(chunks)

    u, v = wrf_wind(raw)
    data_vars: dict[str, xr.DataArray] = {
        "u": u,
        "v": v,
        "w": apply_cf_attrs(
            destagger(_clean(raw["W"]), "level_stag"), _W, overwrite=True
        ),
        "temperature": apply_cf_attrs(
            _clean(wrf_temperature(raw)), "t", overwrite=True
        ),
        "pressure": apply_cf_attrs(
            _clean(raw["P"] + raw["PB"]), _PRESSURE, overwrite=True
        ),
        "pbl_height": apply_cf_attrs(_clean(raw["PBLH"]), "blh", overwrite=True),
        "terrain": apply_cf_attrs(
            _static(_clean(raw["HGT"])), _TERRAIN, overwrite=True
        ),
    }
    for name in variables or ():
        data_vars[name] = raw[name]

    dx, dy = float(raw.attrs["DX"]), float(raw.attrs["DY"])
    z_agl = apply_cf_attrs(_clean(wrf_height_agl(raw)), _Z_AGL, overwrite=True)
    z_agl.attrs["positive"] = "up"
    coords: dict[str, xr.DataArray] = {
        "time": wrf_time(raw),
        "level": xr.DataArray(
            np.arange(raw.sizes["level"]),
            dims="level",
            attrs={"long_name": "Mass level index"},
        ),
        "y": apply_cf_attrs(xr.DataArray(dy * np.arange(raw.sizes["y"]), dims="y"), _Y),
        "x": apply_cf_attrs(xr.DataArray(dx * np.arange(raw.sizes["x"]), dims="x"), _X),
        "z_agl": z_agl,
    }
    for name, key, variable in (("lon", "XLONG", _LON), ("lat", "XLAT", _LAT)):
        if key in raw:
            coords[name] = apply_cf_attrs(
                _static(_clean(raw[key])), variable, overwrite=True
            )

    attrs = {k: raw.attrs[k] for k in _PROJECTION_ATTRS if k in raw.attrs}
    out = xr.Dataset(data_vars, coords=coords, attrs=attrs)
    out.set_close(close)
    return out

destagger

destagger(da: DataArray, dim: str, *, new_dim: str | None = None) -> DataArray

Average a face-staggered variable onto cell centres along dim.

out[..., i, ...] = 0.5 * (da[..., i, ...] + da[..., i + 1, ...]), so the output has da.sizes[dim] - 1 entries along new_dim. Numeric coordinates carrying dim are averaged the same way (non-numeric ones are dropped); other coordinates and attrs are preserved. Exact for fields linear in the staggered coordinate.

Parameters:

Name	Type	Description	Default
da	DataArray	Staggered variable (eager or dask backed; stays lazy).	required
dim	str	Staggered dimension to average over.	required
new_dim	str | None	Name of the output dimension; defaults to dim with a trailing _stag stripped.	None

Returns:

Type	Description
DataArray	The cell-centred variable.

Raises:

Type	Description
ValueError	If dim is not a dimension of da.

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/wrf.py

def destagger(
    da: xr.DataArray, dim: str, *, new_dim: str | None = None
) -> xr.DataArray:
    """Average a face-staggered variable onto cell centres along ``dim``.

    ``out[..., i, ...] = 0.5 * (da[..., i, ...] + da[..., i + 1, ...])``, so
    the output has ``da.sizes[dim] - 1`` entries along ``new_dim``. Numeric
    coordinates carrying ``dim`` are averaged the same way (non-numeric ones
    are dropped); other coordinates and ``attrs`` are preserved. Exact for
    fields linear in the staggered coordinate.

    Args:
        da: Staggered variable (eager or dask backed; stays lazy).
        dim: Staggered dimension to average over.
        new_dim: Name of the output dimension; defaults to ``dim`` with a
            trailing ``_stag`` stripped.

    Returns:
        The cell-centred variable.

    Raises:
        ValueError: If ``dim`` is not a dimension of ``da``.
    """
    if dim not in da.dims:
        raise ValueError(f"dim {dim!r} not in DataArray dims {tuple(da.dims)}")
    if new_dim is None:
        new_dim = dim.removesuffix("_stag")
    lo, hi = {dim: slice(None, -1)}, {dim: slice(1, None)}

    def _centre(var: xr.Variable) -> xr.Variable:
        return _renamed(0.5 * (var.isel(lo) + var.isel(hi)), dim, new_dim)

    coords: dict[str, xr.Variable] = {}
    for name, coord in da.coords.items():
        if dim not in coord.dims:
            coords[str(name)] = coord.variable
        elif np.issubdtype(coord.dtype, np.number):
            coords[str(name) if name != dim else new_dim] = _centre(coord.variable)
    return xr.DataArray(
        _centre(da.variable), coords=coords, name=da.name, attrs=da.attrs
    )

wrf_time

wrf_time(ds: Dataset) -> DataArray

Time axis of a wrfout Dataset as datetime64[s].

Decodes the Times char rows (YYYY-MM-DD_HH:MM:SS) when present; otherwise uses XTIME — raw minutes since the SIMULATION_START_DATE global attribute, or datetimes if a caller already decoded it.

Parameters:

Name	Type	Description	Default
ds	Dataset	Raw wrfout Dataset (dimension names may be renamed already).	required

Returns:

Type	Description
DataArray	1-D datetime64[s] DataArray named "time" on the source's
DataArray	time dimension.

Raises:

Type	Description
ValueError	If neither Times nor XTIME is present, or XTIME holds minutes and SIMULATION_START_DATE is missing.

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/wrf.py

def wrf_time(ds: xr.Dataset) -> xr.DataArray:
    """Time axis of a ``wrfout`` Dataset as ``datetime64[s]``.

    Decodes the ``Times`` char rows (``YYYY-MM-DD_HH:MM:SS``) when present;
    otherwise uses ``XTIME`` — raw minutes since the ``SIMULATION_START_DATE``
    global attribute, or datetimes if a caller already decoded it.

    Args:
        ds: Raw ``wrfout`` Dataset (dimension names may be renamed already).

    Returns:
        1-D ``datetime64[s]`` DataArray named ``"time"`` on the source's
        time dimension.

    Raises:
        ValueError: If neither ``Times`` nor ``XTIME`` is present, or
            ``XTIME`` holds minutes and ``SIMULATION_START_DATE`` is missing.
    """
    if "Times" in ds:
        source = ds["Times"]
        stamps = _decode_times(np.asarray(source.values))
    elif "XTIME" in ds:
        source = ds["XTIME"]
        xtime = np.asarray(source.values)
        if np.issubdtype(xtime.dtype, np.datetime64):
            stamps = xtime.astype("datetime64[s]")
        else:
            start = ds.attrs.get("SIMULATION_START_DATE")
            if start is None:
                raise ValueError(
                    "`XTIME` holds minutes but the `SIMULATION_START_DATE` "
                    "global attribute is missing"
                )
            origin = np.datetime64(str(start).strip().replace("_", "T"), "s")
            seconds = np.round(np.asarray(xtime, dtype=np.float64) * 60.0)
            stamps = origin + seconds.astype("timedelta64[s]")
    else:
        raise ValueError(
            "neither `Times` nor `XTIME` present; cannot build a time axis"
        )
    return xr.DataArray(stamps, dims=(source.dims[0],), name="time")

wrf_wind

wrf_wind(ds: Dataset) -> tuple[DataArray, DataArray]

Horizontal wind on mass points, earth-relative where possible.

U / V are destaggered onto mass points. When the file carries the map-rotation fields COSALPHA / SINALPHA the components are rotated from grid-relative to earth-relative, u_e = u cos α − v sin α, v_e = v cos α + u sin α, and stamped with the CF eastward_wind / northward_wind names. Otherwise the grid-relative components are returned with no standard_name and the attr grid_relative = "true".

Parameters:

Name	Type	Description	Default
ds	Dataset	Raw wrfout Dataset (dims renamed) with U and V.	required

Returns:

Type	Description
tuple[DataArray, DataArray]	(u, v) on mass points [m s⁻¹].

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/wrf.py

def wrf_wind(ds: xr.Dataset) -> tuple[xr.DataArray, xr.DataArray]:
    """Horizontal wind on mass points, earth-relative where possible.

    ``U`` / ``V`` are destaggered onto mass points. When the file carries
    the map-rotation fields ``COSALPHA`` / ``SINALPHA`` the components are
    rotated from grid-relative to earth-relative,
    ``u_e = u cos α − v sin α``, ``v_e = v cos α + u sin α``, and stamped
    with the CF ``eastward_wind`` / ``northward_wind`` names. Otherwise the
    grid-relative components are returned with no ``standard_name`` and the
    attr ``grid_relative = "true"``.

    Args:
        ds: Raw ``wrfout`` Dataset (dims renamed) with ``U`` and ``V``.

    Returns:
        ``(u, v)`` on mass points [m s⁻¹].
    """
    u = destagger(_clean(ds["U"]), "x_stag")
    v = destagger(_clean(ds["V"]), "y_stag")
    if "COSALPHA" in ds and "SINALPHA" in ds:
        cos_a = _static(_clean(ds["COSALPHA"]))
        sin_a = _static(_clean(ds["SINALPHA"]))
        u, v = u * cos_a - v * sin_a, v * cos_a + u * sin_a
        return (
            apply_cf_attrs(u, "u", overwrite=True),
            apply_cf_attrs(v, "v", overwrite=True),
        )
    u = apply_cf_attrs(u, _U_GRID, overwrite=True)
    v = apply_cf_attrs(v, _V_GRID, overwrite=True)
    u.attrs["grid_relative"] = "true"
    v.attrs["grid_relative"] = "true"
    return u, v

wrf_height_agl

wrf_height_agl(ds: Dataset) -> DataArray

Height above ground on mass levels.

z_stag = (PH + PHB) / GRAVITY on the staggered levels, averaged to mass levels and reduced by the terrain height: z_agl = 0.5 * (z_stag[k] + z_stag[k + 1]) - HGT. HGT may carry a time axis or not.

Parameters:

Name	Type	Description	Default
ds	Dataset	Raw wrfout Dataset with PH, PHB and HGT.	required

Returns:

Type	Description
DataArray	(time, level, y, x) height above ground [m].

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/wrf.py

def wrf_height_agl(ds: xr.Dataset) -> xr.DataArray:
    """Height above ground on mass levels.

    ``z_stag = (PH + PHB) / GRAVITY`` on the staggered levels, averaged to
    mass levels and reduced by the terrain height:
    ``z_agl = 0.5 * (z_stag[k] + z_stag[k + 1]) - HGT``. ``HGT`` may carry a
    time axis or not.

    Args:
        ds: Raw ``wrfout`` Dataset with ``PH``, ``PHB`` and ``HGT``.

    Returns:
        ``(time, level, y, x)`` height above ground [m].
    """
    z_stag = (ds["PH"] + ds["PHB"]) / GRAVITY
    return destagger(z_stag, _stag_dim(z_stag)) - ds["HGT"]

wrf_temperature

wrf_temperature(ds: Dataset) -> DataArray

Air temperature from WRF's perturbation potential temperature.

T = (T' + THETA_BASE) * ((P + PB) / P_REFERENCE) ** KAPPA [K].

Parameters:

Name	Type	Description	Default
ds	Dataset	Raw wrfout Dataset with T, P and PB.	required

Returns:

Type	Description
DataArray	Air temperature on mass points [K].

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/wrf.py

def wrf_temperature(ds: xr.Dataset) -> xr.DataArray:
    """Air temperature from WRF's perturbation potential temperature.

    ``T = (T' + THETA_BASE) * ((P + PB) / P_REFERENCE) ** KAPPA`` [K].

    Args:
        ds: Raw ``wrfout`` Dataset with ``T``, ``P`` and ``PB``.

    Returns:
        Air temperature on mass points [K].
    """
    pressure = ds["P"] + ds["PB"]
    return (ds["T"] + THETA_BASE) * (pressure / P_REFERENCE) ** KAPPA

Met diagnostics

Wind in the meteorological convention (bearing the wind blows from, degrees clockwise from north),

\[ |V|=\sqrt{u^2+v^2},\qquad \theta_\text{from}=\left(270^\circ-\tfrac{180}{\pi}\operatorname{atan2}(v,u)\right)\bmod 360^\circ,\qquad u=-|V|\sin\theta,\; v=-|V|\cos\theta, \]

the column integral (trapezoid on the coordinate, or \(\sum_k f_k\,\Delta z_k\) with cell widths from the coordinate), the hypsometric equation \(z_{k+1}=z_k+\frac{R_d\,\bar T_{v,k}}{g}\ln\frac{p_k}{p_{k+1}}\), and the bulk-Richardson boundary-layer height \(Ri_b(z)=\frac{g\,(\theta_v-\theta_{v,s})(z-z_s)}{\theta_{v,s}\,(u^2+v^2)}\).

atm

Atmospheric physics operators.

Layer-0 primitives are implemented in :mod:xrops.atm._src and re-exported here.

Content:

• open_wrfout: WRF-ARW wrfout → CF Dataset on (time, level, y, x) (time decode, destagger, physical temperature / pressure, height above ground), with the destagger / wrf_time / wrf_wind / wrf_height_agl / wrf_temperature helpers it is built from.
• Wind diagnostics: wind_speed, wind_direction (meteorological "from" bearing by default) and their inverse wind_components.
• Vertical-column diagnostics: column_integral (trapezoid or cell-width sum), hypsometric_height (pressure levels → geopotential height) and pbl_height_bulk_richardson.
• Layer-1 Operator wrappers (:mod:xrops.atm.operators), all re-exported here.

Trace-gas physics lives under :mod:xrops.atm.gas. Potential temperature is still planned (see docs/design/).

WindSpeed

WindSpeed(u: str = 'u', v: str = 'v', name: str = 'wind_speed')

Bases: Operator

Horizontal wind speed sqrt(u² + v²).

Parameters:

Name	Type	Description	Default
u	str	Name of the eastward wind component.	'u'
v	str	Name of the northward wind component.	'v'
name	str	Name of the output variable.	'wind_speed'

Returns:

Type	Description
	The input dataset with name added.

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

def __init__(self, u: str = "u", v: str = "v", name: str = "wind_speed"):
    self.u = u
    self.v = v
    self.name = name

WindDirection

WindDirection(u: str = 'u', v: str = 'v', name: str = 'wind_direction', convention: Literal['from', 'to'] = 'from')

Bases: Operator

Wind direction in degrees clockwise from north.

Meteorological "from" bearing by default, θ = (270° − atan2(v, u)) mod 360°; convention="to" adds 180°.

Parameters:

Name	Type	Description	Default
u	str	Name of the eastward wind component.	'u'
v	str	Name of the northward wind component.	'v'
name	str	Name of the output variable.	'wind_direction'
convention	Literal['from', 'to']	"from" or "to".	'from'

Returns:

Type	Description
	The input dataset with name added.

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

def __init__(
    self,
    u: str = "u",
    v: str = "v",
    name: str = "wind_direction",
    convention: Literal["from", "to"] = "from",
):
    self.u = u
    self.v = v
    self.name = name
    self.convention = convention

WindComponents

WindComponents(speed: str = 'wind_speed', direction: str = 'wind_direction', convention: Literal['from', 'to'] = 'from', u: str = 'u', v: str = 'v')

Bases: Operator

Eastward / northward components from speed and direction.

Inverse of :class:WindSpeed + :class:WindDirection: u = −|V| sin θ, v = −|V| cos θ for a "from" bearing.

Parameters:

Name	Type	Description	Default
speed	str	Name of the wind-speed variable.	'wind_speed'
direction	str	Name of the direction variable (degrees from north).	'wind_direction'
convention	Literal['from', 'to']	Whether direction is a "from" or "to" bearing.	'from'
u	str	Name of the eastward output component.	'u'
v	str	Name of the northward output component.	'v'

Returns:

Type	Description
	The input dataset with u and v added.

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

def __init__(
    self,
    speed: str = "wind_speed",
    direction: str = "wind_direction",
    convention: Literal["from", "to"] = "from",
    u: str = "u",
    v: str = "v",
):
    self.speed = speed
    self.direction = direction
    self.convention = convention
    self.u = u
    self.v = v

ColumnIntegral

ColumnIntegral(variable: str, *, dim: str = 'z', coord: str | None = None, method: Literal['trapezoid', 'sum'] = 'trapezoid', name: str | None = None)

Bases: Operator

Integrate a variable along its vertical coordinate.

"trapezoid" applies the trapezoid rule on the coordinate values; "sum" is Σ f_k Δz_k with cell widths from the coordinate. See :func:xrops.atm.column_integral.

Parameters:

Name	Type	Description	Default
variable	str	Name of the variable to integrate.	required
dim	str	Vertical dimension.	'z'
coord	str | None	Name of the coordinate locating the levels, or None for the index coordinate on dim.	None
method	Literal['trapezoid', 'sum']	"trapezoid" or "sum".	'trapezoid'
name	str | None	Name of the output variable; defaults to f"{variable}_column".	None

Returns:

Type	Description
	The input dataset with name added (dim reduced away).

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

def __init__(
    self,
    variable: str,
    *,
    dim: str = "z",
    coord: str | None = None,
    method: Literal["trapezoid", "sum"] = "trapezoid",
    name: str | None = None,
):
    self.variable = variable
    self.dim = dim
    self.coord = coord
    self.method = method
    self.name = name

HypsometricHeight

HypsometricHeight(pressure: str = 'pressure', temperature: str = 'temperature', specific_humidity: str | None = None, surface_pressure: str = 'sp', surface_height: str | None = None, level: str = 'level', name: str = 'height')

Bases: Operator

Geopotential height of pressure levels (hypsometric equation).

z_{k+1} = z_k + (R_d T̄_v / g) ln(p_k / p_{k+1}) integrated upward from the surface. See :func:xrops.atm.hypsometric_height.

Parameters:

Name	Type	Description	Default
pressure	str	Pressure per level.	'pressure'
temperature	str	Air temperature per level [K].	'temperature'
specific_humidity	str | None	Optional specific humidity for T_v.	None
surface_pressure	str	Surface pressure.	'sp'
surface_height	str | None	Optional terrain height added to every level.	None
level	str	Vertical dimension.	'level'
name	str	Name of the output variable.	'height'

Returns:

Type	Description
	The input dataset with name added.

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

def __init__(
    self,
    pressure: str = "pressure",
    temperature: str = "temperature",
    specific_humidity: str | None = None,
    surface_pressure: str = "sp",
    surface_height: str | None = None,
    level: str = "level",
    name: str = "height",
):
    self.pressure = pressure
    self.temperature = temperature
    self.specific_humidity = specific_humidity
    self.surface_pressure = surface_pressure
    self.surface_height = surface_height
    self.level = level
    self.name = name

PBLHeightBulkRichardson

PBLHeightBulkRichardson(theta_v: str, u: str = 'u', v: str = 'v', height: str = 'z_agl', ri_crit: float = 0.25, level: str = 'level', name: str = 'pbl_height')

Bases: Operator

Boundary-layer height from the bulk Richardson number.

First height where Ri_b = g (θ_v − θ_{v,s}) (z − z_s) / (θ_{v,s} (u² + v²)) reaches ri_crit, interpolated between levels; NaN when no level crosses. See :func:xrops.atm.pbl_height_bulk_richardson.

Parameters:

Name	Type	Description	Default
theta_v	str	Virtual potential temperature per level [K].	required
u	str	Eastward wind per level.	'u'
v	str	Northward wind per level.	'v'
height	str	Height above ground per level [m].	'z_agl'
ri_crit	float	Critical bulk Richardson number.	0.25
level	str	Vertical dimension.	'level'
name	str	Name of the output variable.	'pbl_height'

Returns:

Type	Description
	The input dataset with name added (level reduced away).

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

def __init__(
    self,
    theta_v: str,
    u: str = "u",
    v: str = "v",
    height: str = "z_agl",
    ri_crit: float = 0.25,
    level: str = "level",
    name: str = "pbl_height",
):
    self.theta_v = theta_v
    self.u = u
    self.v = v
    self.height = height
    self.ri_crit = ri_crit
    self.level = level
    self.name = name

wind_speed

wind_speed(ds: Dataset, *, u: str = 'u', v: str = 'v', name: str = 'wind_speed') -> Dataset

Horizontal wind speed |V| = sqrt(u² + v²).

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset holding the eastward and northward wind components.	required
u	str	Name of the eastward component.	'u'
v	str	Name of the northward component.	'v'
name	str	Name of the output variable.	'wind_speed'

Returns:

Type	Description
Dataset	ds with name added (CF wind_speed, same units as the
Dataset	components).

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/wind.py

def wind_speed(
    ds: xr.Dataset,
    *,
    u: str = "u",
    v: str = "v",
    name: str = "wind_speed",
) -> xr.Dataset:
    """Horizontal wind speed ``|V| = sqrt(u² + v²)``.

    Args:
        ds: Dataset holding the eastward and northward wind components.
        u: Name of the eastward component.
        v: Name of the northward component.
        name: Name of the output variable.

    Returns:
        ``ds`` with ``name`` added (CF ``wind_speed``, same units as the
        components).
    """
    speed = (ds[u] ** 2 + ds[v] ** 2) ** 0.5
    out = ds.copy()
    out[name] = apply_cf_attrs(speed, _WIND_SPEED, overwrite=True)
    return out

wind_direction

wind_direction(ds: Dataset, *, u: str = 'u', v: str = 'v', name: str = 'wind_direction', convention: Literal['from', 'to'] = 'from') -> Dataset

Wind direction in degrees clockwise from north.

θ_from = (270° − atan2(v, u)·180/π) mod 360° is the bearing the wind blows from (meteorological convention); convention="to" returns θ_to = (θ_from + 180°) mod 360°. Calm air (u = v = 0) maps to 270° under "from" (atan2(0, 0) = 0).

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset holding the wind components.	required
u	str	Name of the eastward component.	'u'
v	str	Name of the northward component.	'v'
name	str	Name of the output variable.	'wind_direction'
convention	Literal['from', 'to']	"from" (CF wind_from_direction) or "to" (CF wind_to_direction).	'from'

Returns:

Type	Description
Dataset	ds with name added, in degrees on [0, 360).

Raises:

Type	Description
ValueError	If convention is not "from" or "to".

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/wind.py

def wind_direction(
    ds: xr.Dataset,
    *,
    u: str = "u",
    v: str = "v",
    name: str = "wind_direction",
    convention: Literal["from", "to"] = "from",
) -> xr.Dataset:
    """Wind direction in degrees clockwise from north.

    ``θ_from = (270° − atan2(v, u)·180/π) mod 360°`` is the bearing the wind
    blows *from* (meteorological convention); ``convention="to"`` returns
    ``θ_to = (θ_from + 180°) mod 360°``. Calm air (``u = v = 0``) maps to
    ``270°`` under ``"from"`` (``atan2(0, 0) = 0``).

    Args:
        ds: Dataset holding the wind components.
        u: Name of the eastward component.
        v: Name of the northward component.
        name: Name of the output variable.
        convention: ``"from"`` (CF ``wind_from_direction``) or ``"to"``
            (CF ``wind_to_direction``).

    Returns:
        ``ds`` with ``name`` added, in degrees on ``[0, 360)``.

    Raises:
        ValueError: If ``convention`` is not ``"from"`` or ``"to"``.
    """
    if convention not in ("from", "to"):
        raise ValueError(f"convention must be 'from' or 'to', got {convention!r}")
    theta = np.mod(270.0 - np.degrees(np.arctan2(ds[v], ds[u])), 360.0)
    if convention == "to":
        theta = np.mod(theta + 180.0, 360.0)
    out = ds.copy()
    variable = _WIND_FROM if convention == "from" else _WIND_TO
    out[name] = apply_cf_attrs(theta, variable, overwrite=True)
    return out

wind_components

wind_components(ds: Dataset, *, speed: str = 'wind_speed', direction: str = 'wind_direction', convention: Literal['from', 'to'] = 'from', u: str = 'u', v: str = 'v') -> Dataset

Eastward / northward components from speed and direction.

Inverse of :func:wind_speed + :func:wind_direction. With the meteorological "from" bearing θ: u = −|V| sin θ, v = −|V| cos θ; under "to" the signs flip.

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset holding the speed and direction.	required
speed	str	Name of the wind-speed variable.	'wind_speed'
direction	str	Name of the direction variable (degrees from north).	'wind_direction'
convention	Literal['from', 'to']	Whether direction is a "from" or "to" bearing.	'from'
u	str	Name of the eastward output component.	'u'
v	str	Name of the northward output component.	'v'

Returns:

Type	Description
Dataset	ds with u and v added (CF eastward_wind /
Dataset	northward_wind).

Raises:

Type	Description
ValueError	If convention is not "from" or "to".

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/wind.py

def wind_components(
    ds: xr.Dataset,
    *,
    speed: str = "wind_speed",
    direction: str = "wind_direction",
    convention: Literal["from", "to"] = "from",
    u: str = "u",
    v: str = "v",
) -> xr.Dataset:
    """Eastward / northward components from speed and direction.

    Inverse of :func:`wind_speed` + :func:`wind_direction`. With the
    meteorological ``"from"`` bearing ``θ``:
    ``u = −|V| sin θ``, ``v = −|V| cos θ``; under ``"to"`` the signs flip.

    Args:
        ds: Dataset holding the speed and direction.
        speed: Name of the wind-speed variable.
        direction: Name of the direction variable (degrees from north).
        convention: Whether ``direction`` is a ``"from"`` or ``"to"`` bearing.
        u: Name of the eastward output component.
        v: Name of the northward output component.

    Returns:
        ``ds`` with ``u`` and ``v`` added (CF ``eastward_wind`` /
        ``northward_wind``).

    Raises:
        ValueError: If ``convention`` is not ``"from"`` or ``"to"``.
    """
    if convention not in ("from", "to"):
        raise ValueError(f"convention must be 'from' or 'to', got {convention!r}")
    theta = np.radians(ds[direction])
    sign = -1.0 if convention == "from" else 1.0
    out = ds.copy()
    out[u] = apply_cf_attrs(sign * ds[speed] * np.sin(theta), "u", overwrite=True)
    out[v] = apply_cf_attrs(sign * ds[speed] * np.cos(theta), "v", overwrite=True)
    return out

column_integral

column_integral(da: DataArray, *, dim: str = 'z', coord: str | DataArray | None = None, method: Literal['trapezoid', 'sum'] = 'trapezoid') -> DataArray

Integrate da along the vertical coordinate.

method="trapezoid" applies the trapezoid rule on the coordinate values (exact for fields linear in z, the right choice for endpoint-inclusive level grids):

∫ f dz ≈ Σ_k ½ (f_k + f_{k+1}) (z_{k+1} − z_k)

method="sum" is the cell-centred Riemann sum Σ_k f_k Δz_k with Δz_k the cell width around level k (half-way to each neighbour, half a spacing at the ends). Both methods give c (z_1 − z_0) for a uniform field c on [z_0, z_1] — unlike sum · dz with a constant dz, which over-counts an endpoint-inclusive grid by n_z / (n_z − 1). The sign follows the coordinate's orientation.

Parameters:

Name	Type	Description	Default
da	DataArray	Field to integrate (any dims, must include dim).	required
dim	str	Vertical dimension.	'z'
coord	str | DataArray | None	Coordinate giving the position of each level: None uses the index coordinate on dim; a str names a coordinate of da (may be multi-dimensional, e.g. a terrain-following z_agl); a DataArray is used directly and must carry dim.	None
method	Literal['trapezoid', 'sum']	"trapezoid" or "sum".	'trapezoid'

Returns:

Type	Description
DataArray	da reduced over dim; units are da units × coordinate
DataArray	units. A NaN at any level makes the whole column NaN.

Raises:

Type	Description
ValueError	If dim is missing, has fewer than two levels, no coordinate is available, or method is unknown.

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/vertical.py

def column_integral(
    da: xr.DataArray,
    *,
    dim: str = "z",
    coord: str | xr.DataArray | None = None,
    method: Literal["trapezoid", "sum"] = "trapezoid",
) -> xr.DataArray:
    """Integrate ``da`` along the vertical coordinate.

    ``method="trapezoid"`` applies the trapezoid rule on the coordinate
    values (exact for fields linear in ``z``, the right choice for
    endpoint-inclusive level grids):

    ``∫ f dz ≈ Σ_k ½ (f_k + f_{k+1}) (z_{k+1} − z_k)``

    ``method="sum"`` is the cell-centred Riemann sum ``Σ_k f_k Δz_k`` with
    ``Δz_k`` the cell width around level ``k`` (half-way to each neighbour,
    half a spacing at the ends). Both methods give ``c (z_1 − z_0)`` for a
    uniform field ``c`` on ``[z_0, z_1]`` — unlike ``sum · dz`` with a
    constant ``dz``, which over-counts an endpoint-inclusive grid by
    ``n_z / (n_z − 1)``. The sign follows the coordinate's orientation.

    Args:
        da: Field to integrate (any dims, must include ``dim``).
        dim: Vertical dimension.
        coord: Coordinate giving the position of each level: ``None`` uses
            the index coordinate on ``dim``; a ``str`` names a coordinate
            of ``da`` (may be multi-dimensional, e.g. a terrain-following
            ``z_agl``); a ``DataArray`` is used directly and must carry
            ``dim``.
        method: ``"trapezoid"`` or ``"sum"``.

    Returns:
        ``da`` reduced over ``dim``; units are ``da`` units × coordinate
        units. A NaN at any level makes the whole column NaN.

    Raises:
        ValueError: If ``dim`` is missing, has fewer than two levels, no
            coordinate is available, or ``method`` is unknown.
    """
    if dim not in da.dims:
        raise ValueError(f"dim {dim!r} not in DataArray dims {tuple(da.dims)}")
    if da.sizes[dim] < 2:
        raise ValueError(f"column_integral needs at least two levels along {dim!r}")
    if method not in ("trapezoid", "sum"):
        raise ValueError(f"method must be 'trapezoid' or 'sum', got {method!r}")
    if coord is None:
        if dim not in da.coords:
            raise ValueError(
                f"no coordinate on {dim!r}; pass `coord=` to locate the levels"
            )
        z = da[dim]
    elif isinstance(coord, str):
        z = da[coord]
    else:
        z = coord
    if dim not in z.dims:
        raise ValueError(f"coordinate {z.name!r} does not carry dim {dim!r}")

    f = _positional(da, dim)
    z = _positional(z, dim)
    if method == "trapezoid":
        f_lo, f_hi = _lo_hi(f, dim)
        z_lo, z_hi = _lo_hi(z, dim)
        out = (0.5 * (f_lo + f_hi) * (z_hi - z_lo)).sum(dim, skipna=False)
    else:
        out = (f * _cell_widths(z, dim)).sum(dim, skipna=False)
    return out.rename(da.name)

hypsometric_height

hypsometric_height(ds: Dataset, *, pressure: str = 'pressure', temperature: str = 'temperature', specific_humidity: str | None = None, surface_pressure: str = 'sp', surface_height: str | None = None, level: str = 'level', name: str = 'height') -> Dataset

Geopotential height of pressure levels from the hypsometric equation.

Integrates layer by layer upward from the surface with the virtual temperature T_v = T (1 + 0.61 q):

z_{k+1} = z_k + (R_d T̄_{v,k} / g) ln(p_k / p_{k+1})

where T̄_{v,k} is the mean of the two bounding levels. The first level is anchored to the surface, z_0 = z_s + (R_d T_{v,0} / g) ln(p_s / p_0), using that level's virtual temperature for the surface layer. For an isothermal column this reduces to z = (R_d T / g) ln(p_s / p) exactly. Levels with p > p_s get negative heights (extrapolated below ground). A NaN temperature at a level invalidates that level and every level above it.

pressure and surface_pressure must share units (only ratios enter). A 1-D pressure coordinate may be ordered either way; a multi-dimensional pressure field must be ordered surface first.

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset with the pressure, temperature and surface-pressure fields (pressure may be the level coordinate itself).	required
pressure	str	Pressure per level.	'pressure'
temperature	str	Air temperature per level [K].	'temperature'
specific_humidity	str | None	Optional specific humidity [kg kg⁻¹] for the virtual-temperature correction; None treats the air as dry.	None
surface_pressure	str	Surface pressure (no level dim).	'sp'
surface_height	str | None	Optional terrain height [m] added to every level; None measures height above the surface.	None
level	str	Vertical dimension.	'level'
name	str	Name of the output variable.	'height'

Returns:

Type	Description
Dataset	ds with name added (CF geopotential_height, m).

Raises:

Type	Description
ValueError	If pressure or temperature lack level.

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/vertical.py

def hypsometric_height(
    ds: xr.Dataset,
    *,
    pressure: str = "pressure",
    temperature: str = "temperature",
    specific_humidity: str | None = None,
    surface_pressure: str = "sp",
    surface_height: str | None = None,
    level: str = "level",
    name: str = "height",
) -> xr.Dataset:
    """Geopotential height of pressure levels from the hypsometric equation.

    Integrates layer by layer upward from the surface with the virtual
    temperature ``T_v = T (1 + 0.61 q)``:

    ``z_{k+1} = z_k + (R_d T̄_{v,k} / g) ln(p_k / p_{k+1})``

    where ``T̄_{v,k}`` is the mean of the two bounding levels. The first
    level is anchored to the surface, ``z_0 = z_s + (R_d T_{v,0} / g)
    ln(p_s / p_0)``, using that level's virtual temperature for the
    surface layer. For an isothermal column this reduces to
    ``z = (R_d T / g) ln(p_s / p)`` exactly. Levels with ``p > p_s`` get
    negative heights (extrapolated below ground). A NaN temperature at a
    level invalidates that level and every level above it.

    ``pressure`` and ``surface_pressure`` must share units (only ratios
    enter). A 1-D pressure coordinate may be ordered either way; a
    multi-dimensional pressure field must be ordered surface first.

    Args:
        ds: Dataset with the pressure, temperature and surface-pressure
            fields (``pressure`` may be the ``level`` coordinate itself).
        pressure: Pressure per level.
        temperature: Air temperature per level [K].
        specific_humidity: Optional specific humidity [kg kg⁻¹] for the
            virtual-temperature correction; ``None`` treats the air as dry.
        surface_pressure: Surface pressure (no ``level`` dim).
        surface_height: Optional terrain height [m] added to every level;
            ``None`` measures height above the surface.
        level: Vertical dimension.
        name: Name of the output variable.

    Returns:
        ``ds`` with ``name`` added (CF ``geopotential_height``, m).

    Raises:
        ValueError: If ``pressure`` or ``temperature`` lack ``level``.
    """
    p = ds[pressure]
    t = ds[temperature]
    for label, da in ((pressure, p), (temperature, t)):
        if level not in da.dims:
            raise ValueError(f"{label!r} does not carry the level dim {level!r}")
    tv = t * (1.0 + EPSILON_TV * ds[specific_humidity]) if specific_humidity else t
    ps = ds[surface_pressure]
    # 1-D pressure axes may run top-down (CDS/ERA5 style); integrate
    # surface-first and flip the result back into the dataset's order.
    flip = p.ndim == 1 and bool(p.values[0] < p.values[-1])
    if flip:
        p = p.isel({level: slice(None, None, -1)})
        tv = tv.isel({level: slice(None, None, -1)})

    p_pos = _positional(p, level)
    tv_pos = _positional(tv, level)
    p_lo, p_hi = _lo_hi(p_pos, level)
    tv_lo, tv_hi = _lo_hi(tv_pos, level)
    dz_layers = (R_D / G) * 0.5 * (tv_lo + tv_hi) * np.log(p_lo / p_hi)
    first = {level: slice(0, 1)}
    dz_first = (R_D / G) * tv_pos.isel(first) * np.log(ps / p_pos.isel(first))
    z = xr.concat([dz_first, dz_layers], dim=level).cumsum(level, skipna=False)
    if surface_height is not None:
        z = z + ds[surface_height]
    if flip:
        z = z.isel({level: slice(None, None, -1)})
    z = z.transpose(*[d for d in t.dims if d in z.dims], ...)
    if level in ds.coords:
        z = z.assign_coords({level: ds[level]})
    out = ds.copy()
    out[name] = apply_cf_attrs(z, _GEOPOTENTIAL_HEIGHT, overwrite=True)
    return out

pbl_height_bulk_richardson

pbl_height_bulk_richardson(ds: Dataset, *, theta_v: str, u: str = 'u', v: str = 'v', height: str = 'z_agl', ri_crit: float = 0.25, level: str = 'level', name: str = 'pbl_height') -> Dataset

Boundary-layer height from the bulk Richardson number.

Vogelezang & Holtslag (1996) form without the surface-friction term, with the surface values taken from the lowest level s:

Ri_b(z) = g (θ_v(z) − θ_{v,s}) (z − z_s) / (θ_{v,s} ((u(z) − u_s)² + (v(z) − v_s)²))

h is the first height (scanning upward from the level above the surface) where Ri_b ≥ Ri_crit, linearly interpolated between the bracketing levels. Columns that never reach Ri_crit get NaN rather than the top level. The surface level itself is assigned Ri_b = 0; a level aloft with no shear relative to the surface (u = u_s, v = v_s) takes the sign of the numerator: +∞ when stably stratified (so h stops at the level below it), −∞ when unstable (never a crossing) and 0 when θ_v = θ_{v,s}.

The vertical ordering is read from height: profiles may be surface first or top first (they are flipped internally), but every column must share one orientation.

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset with virtual potential temperature, wind components and height per level.	required
theta_v	str	Virtual potential temperature [K].	required
u	str	Eastward wind [m s⁻¹].	'u'
v	str	Northward wind [m s⁻¹].	'v'
height	str	Height above ground per level [m] (may be a coordinate, 1-D or per column).	'z_agl'
ri_crit	float	Critical bulk Richardson number.	0.25
level	str	Vertical dimension.	'level'
name	str	Name of the output variable.	'pbl_height'

Returns:

Type	Description
Dataset	ds with name added (CF
Dataset	atmosphere_boundary_layer_thickness, m); the level dim is
Dataset	reduced away.

Raises:

Type	Description
ValueError	If theta_v or height lack level, or if the vertical orientation of height differs between columns.

Source code in packages/xrtoolz-ops/src/xrops/atm/_src/vertical.py

def pbl_height_bulk_richardson(
    ds: xr.Dataset,
    *,
    theta_v: str,
    u: str = "u",
    v: str = "v",
    height: str = "z_agl",
    ri_crit: float = 0.25,
    level: str = "level",
    name: str = "pbl_height",
) -> xr.Dataset:
    """Boundary-layer height from the bulk Richardson number.

    Vogelezang & Holtslag (1996) form without the surface-friction term,
    with the surface values taken from the lowest level ``s``:

    ``Ri_b(z) = g (θ_v(z) − θ_{v,s}) (z − z_s)
    / (θ_{v,s} ((u(z) − u_s)² + (v(z) − v_s)²))``

    ``h`` is the first height (scanning upward from the level above the
    surface) where ``Ri_b ≥ Ri_crit``, linearly interpolated between the
    bracketing levels. Columns that never reach ``Ri_crit`` get NaN rather
    than the top level. The surface level itself is assigned ``Ri_b = 0``;
    a level aloft with no shear relative to the surface (``u = u_s``,
    ``v = v_s``) takes the sign of the numerator: ``+∞`` when stably
    stratified (so ``h`` stops at the level below it), ``−∞`` when
    unstable (never a crossing) and ``0`` when ``θ_v = θ_{v,s}``.

    The vertical ordering is read from ``height``: profiles may be surface
    first or top first (they are flipped internally), but every column must
    share one orientation.

    Args:
        ds: Dataset with virtual potential temperature, wind components and
            height per level.
        theta_v: Virtual potential temperature [K].
        u: Eastward wind [m s⁻¹].
        v: Northward wind [m s⁻¹].
        height: Height above ground per level [m] (may be a coordinate,
            1-D or per column).
        ri_crit: Critical bulk Richardson number.
        level: Vertical dimension.
        name: Name of the output variable.

    Returns:
        ``ds`` with ``name`` added (CF
        ``atmosphere_boundary_layer_thickness``, m); the ``level`` dim is
        reduced away.

    Raises:
        ValueError: If ``theta_v`` or ``height`` lack ``level``, or if the
            vertical orientation of ``height`` differs between columns.
    """
    thv = ds[theta_v]
    z = ds[height]
    for label, da in ((theta_v, thv), (height, z)):
        if level not in da.dims:
            raise ValueError(f"{label!r} does not carry the level dim {level!r}")
    # Orientation from the height field: flip top-first profiles so the
    # surface is level 0. Columns whose end levels are NaN are ignored.
    z_first, z_last = z.isel({level: 0}), z.isel({level: -1})
    n_top_first = int((z_first > z_last).sum())
    n_surface_first = int((z_first < z_last).sum())
    if n_top_first and n_surface_first:
        raise ValueError(
            f"{height!r} is surface-first in some columns and top-first in "
            "others; the vertical orientation must be consistent"
        )
    if n_top_first:
        ds = ds.isel({level: slice(None, None, -1)})
        thv, z = ds[theta_v], ds[height]
    surface = {level: 0}
    thv_s = thv.isel(surface)
    z_s = z.isel(surface)
    shear2 = (ds[u] - ds[u].isel(surface)) ** 2 + (ds[v] - ds[v].isel(surface)) ** 2
    numerator = G * (thv - thv_s) * (z - z_s)
    # Calm levels (zero shear relative to the surface) have infinite Ri_b
    # with the sign of the numerator (0 where the numerator is zero too):
    # mask the denominator so the division never warns, then restore the
    # limit where the shear was zero (NaN inputs stay NaN). The surface
    # level is zeroed in the kernel.
    calm = shear2 == 0.0
    calm_ri = xr.where(numerator > 0.0, np.inf, 0.0).where(numerator >= 0.0, -np.inf)
    calm_ri = calm_ri.where(numerator.notnull())
    ri = (numerator / (thv_s * shear2.where(~calm))).where(~calm, calm_ri)
    ri, z_b = xr.broadcast(_positional(ri, level), _positional(z, level))
    if ri.chunks is not None or z_b.chunks is not None:
        # The kernel needs the whole column: one chunk along ``level``.
        ri, z_b = ri.chunk({level: -1}), z_b.chunk({level: -1})
    h = xr.apply_ufunc(
        _first_crossing,
        ri,
        z_b,
        input_core_dims=[[level], [level]],
        kwargs={"ri_crit": ri_crit},
        dask="parallelized",
        output_dtypes=[np.float64],
    )
    out = ds.copy()
    out[name] = apply_cf_attrs(h, "blh", overwrite=True)
    return out

Methane columns

Column averaging kernel in the TROPOMI / OCO convention, \(\hat y=\mathbf h^\top\mathbf x_a+\sum_l h_l A_l (x_l-x_{a,l})\), the hydrostatic dry-air column \(N_\text{dry}=\frac{N_A}{g\,M_\text{dry}}\int_0^{p_s}(1-q)\,dp\), the mixing-ratio column \(\Omega=\frac{N_A}{g\,M_\text{dry}}\int\chi\,(1-q)\,dp\), and the single-layer column mass to \(\Delta\chi\) conversion \(\Delta\chi=\frac{m\,N_A/M_\text{gas}}{n_\text{air}\,L}\) with \(n_\text{air}=p/(k_B T)\).

ch4

Methane-specific physics operators.

Layer-0 primitives are implemented in :mod:xrops.atm.gas.ch4._src and re-exported here.

Content:

• apply_column_averaging_kernel: TROPOMI / OCO-style smoothing of a model profile with a column averaging kernel and pressure weights.
• dry_air_column: dry-air column number density from surface pressure (and optional total column water vapour).
• mixing_ratio_to_column: dry-air mole-fraction profile → column number density; column_mass_to_delta_vmr: column mass enhancement → mixing-ratio enhancement over a well-mixed layer.
• Layer-1 Operator wrappers (:mod:xrops.atm.gas.ch4.operators), all re-exported here.

ApplyColumnAveragingKernel

ApplyColumnAveragingKernel(profile: str, prior: str, averaging_kernel: str, pressure_weights: str, level: str = 'layer', name: str = 'xch4_smoothed')

Bases: Operator

Smooth a model profile with a satellite column averaging kernel.

ŷ = hᵀ x_a + Σ_l h_l A_l (x_l − x_{a,l}). See :func:xrops.atm.gas.ch4.apply_column_averaging_kernel.

Parameters:

Name	Type	Description	Default
profile	str	Model profile x per layer.	required
prior	str	Retrieval prior x_a per layer.	required
averaging_kernel	str	Column averaging kernel A per layer.	required
pressure_weights	str	Pressure weights h per layer.	required
level	str	Layer dimension.	'layer'
name	str	Name of the output variable.	'xch4_smoothed'

Returns:

Type	Description
	The input dataset with name added (level reduced away).

Source code in packages/xrtoolz-ops/src/xrops/atm/gas/ch4/operators.py

def __init__(
    self,
    profile: str,
    prior: str,
    averaging_kernel: str,
    pressure_weights: str,
    level: str = "layer",
    name: str = "xch4_smoothed",
):
    self.profile = profile
    self.prior = prior
    self.averaging_kernel = averaging_kernel
    self.pressure_weights = pressure_weights
    self.level = level
    self.name = name

DryAirColumn

DryAirColumn(surface_pressure: str = 'sp', water_vapour_column: str | None = None, name: str = 'dry_air_column')

Bases: Operator

Dry-air column number density N_A (p_s / g − W) / M_dry.

See :func:xrops.atm.gas.ch4.dry_air_column.

Parameters:

Name	Type	Description	Default
surface_pressure	str	Surface pressure [Pa].	'sp'
water_vapour_column	str | None	Optional total column water vapour [kg m⁻²].	None
name	str	Name of the output variable.	'dry_air_column'

Returns:

Type	Description
	The input dataset with name added [molecules m⁻²].

Source code in packages/xrtoolz-ops/src/xrops/atm/gas/ch4/operators.py

def __init__(
    self,
    surface_pressure: str = "sp",
    water_vapour_column: str | None = None,
    name: str = "dry_air_column",
):
    self.surface_pressure = surface_pressure
    self.water_vapour_column = water_vapour_column
    self.name = name

MixingRatioToColumn

MixingRatioToColumn(vmr: str, pressure: str = 'pressure', specific_humidity: str | None = None, level: str = 'level', name: str = 'column')

Bases: Operator

Column number density Ω = N_A / (g M_dry) ∫ χ (1 − q) dp.

See :func:xrops.atm.gas.ch4.mixing_ratio_to_column.

Parameters:

Name	Type	Description	Default
vmr	str	Dry-air mole fraction per level.	required
pressure	str	Pressure per level [Pa].	'pressure'
specific_humidity	str | None	Optional specific humidity per level.	None
level	str	Vertical dimension.	'level'
name	str	Name of the output variable.	'column'

Returns:

Type	Description
	The input dataset with name added (level reduced away).

Source code in packages/xrtoolz-ops/src/xrops/atm/gas/ch4/operators.py

def __init__(
    self,
    vmr: str,
    pressure: str = "pressure",
    specific_humidity: str | None = None,
    level: str = "level",
    name: str = "column",
):
    self.vmr = vmr
    self.pressure = pressure
    self.specific_humidity = specific_humidity
    self.level = level
    self.name = name

apply_column_averaging_kernel

apply_column_averaging_kernel(ds: Dataset, *, profile: str, prior: str, averaging_kernel: str, pressure_weights: str, level: str = 'layer', name: str = 'xch4_smoothed') -> Dataset

Smooth a model profile with a satellite column averaging kernel.

TROPOMI / OCO convention, with pressure weights h summing to one, column averaging kernel A per layer and retrieval prior x_a:

ŷ = hᵀ x_a + Σ_l h_l A_l (x_l − x_{a,l})

A ≡ 1 reduces to the column average hᵀ x; x = x_a returns the prior column hᵀ x_a.

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset with the four per-layer fields.	required
profile	str	Model profile x per layer (mole fraction).	required
prior	str	Retrieval prior x_a per layer, same units as profile.	required
averaging_kernel	str	Column averaging kernel A per layer.	required
pressure_weights	str	Pressure weights h per layer (sum to one).	required
level	str	Layer dimension.	'layer'
name	str	Name of the output variable.	'xch4_smoothed'

Returns:

Type	Description
Dataset	ds with name added (CF
Dataset	dry_atmosphere_mole_fraction_of_methane, units copied from
Dataset	profile); the level dim is reduced away. A NaN in any layer
Dataset	(e.g. a QA-masked pixel) makes the result NaN.

Raises:

Type	Description
ValueError	If profile lacks level.

Source code in packages/xrtoolz-ops/src/xrops/atm/gas/ch4/_src/column.py

def apply_column_averaging_kernel(
    ds: xr.Dataset,
    *,
    profile: str,
    prior: str,
    averaging_kernel: str,
    pressure_weights: str,
    level: str = "layer",
    name: str = "xch4_smoothed",
) -> xr.Dataset:
    """Smooth a model profile with a satellite column averaging kernel.

    TROPOMI / OCO convention, with pressure weights ``h`` summing to one,
    column averaging kernel ``A`` per layer and retrieval prior ``x_a``:

    ``ŷ = hᵀ x_a + Σ_l h_l A_l (x_l − x_{a,l})``

    ``A ≡ 1`` reduces to the column average ``hᵀ x``; ``x = x_a`` returns
    the prior column ``hᵀ x_a``.

    Args:
        ds: Dataset with the four per-layer fields.
        profile: Model profile ``x`` per layer (mole fraction).
        prior: Retrieval prior ``x_a`` per layer, same units as ``profile``.
        averaging_kernel: Column averaging kernel ``A`` per layer.
        pressure_weights: Pressure weights ``h`` per layer (sum to one).
        level: Layer dimension.
        name: Name of the output variable.

    Returns:
        ``ds`` with ``name`` added (CF
        ``dry_atmosphere_mole_fraction_of_methane``, units copied from
        ``profile``); the ``level`` dim is reduced away. A NaN in any layer
        (e.g. a QA-masked pixel) makes the result NaN.

    Raises:
        ValueError: If ``profile`` lacks ``level``.
    """
    x = ds[profile]
    if level not in x.dims:
        raise ValueError(f"{profile!r} does not carry the layer dim {level!r}")
    x_a = ds[prior]
    h = ds[pressure_weights]
    a = ds[averaging_kernel]
    smoothed = (h * x_a).sum(level, skipna=False) + (h * a * (x - x_a)).sum(
        level, skipna=False
    )
    variable = Variable(
        name=name,
        standard_name="dry_atmosphere_mole_fraction_of_methane",
        long_name="Averaging-kernel-smoothed XCH4",
        units=x.attrs.get("units"),
    )
    out = ds.copy()
    out[name] = apply_cf_attrs(smoothed, variable, overwrite=True)
    return out

dry_air_column

dry_air_column(ds: Dataset, *, surface_pressure: str = 'sp', water_vapour_column: str | None = None, name: str = 'dry_air_column') -> Dataset

Dry-air column number density from surface pressure.

Hydrostatic balance gives the total-air column as p_s / g kg m⁻²; removing the water-vapour mass column W (total column water vapour, kg m⁻²) and converting to molecules with N_A / M_dry:

N_dry = N_A / (g M_dry) ∫₀^{p_s} (1 − q) dp = N_A (p_s / g − W) / M_dry

With W = 0 and p_s = 10⁵ Pa this is ≈ 2.12 × 10²⁹ m⁻².

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset with the surface pressure (and optional water column).	required
surface_pressure	str	Surface pressure [Pa].	'sp'
water_vapour_column	str | None	Optional total column water vapour [kg m⁻²] (CF atmosphere_mass_content_of_water_vapor, ERA5 tcwv); None treats the column as dry.	None
name	str	Name of the output variable.	'dry_air_column'

Returns:

Type	Description
Dataset	ds with name added [molecules m⁻²].

Source code in packages/xrtoolz-ops/src/xrops/atm/gas/ch4/_src/column.py

def dry_air_column(
    ds: xr.Dataset,
    *,
    surface_pressure: str = "sp",
    water_vapour_column: str | None = None,
    name: str = "dry_air_column",
) -> xr.Dataset:
    """Dry-air column number density from surface pressure.

    Hydrostatic balance gives the total-air column as ``p_s / g`` kg m⁻²;
    removing the water-vapour mass column ``W`` (total column water vapour,
    kg m⁻²) and converting to molecules with ``N_A / M_dry``:

    ``N_dry = N_A / (g M_dry) ∫₀^{p_s} (1 − q) dp = N_A (p_s / g − W) / M_dry``

    With ``W = 0`` and ``p_s = 10⁵ Pa`` this is ``≈ 2.12 × 10²⁹ m⁻²``.

    Args:
        ds: Dataset with the surface pressure (and optional water column).
        surface_pressure: Surface pressure [Pa].
        water_vapour_column: Optional total column water vapour [kg m⁻²]
            (CF ``atmosphere_mass_content_of_water_vapor``, ERA5 ``tcwv``);
            ``None`` treats the column as dry.
        name: Name of the output variable.

    Returns:
        ``ds`` with ``name`` added [molecules m⁻²].
    """
    mass = ds[surface_pressure] / G
    if water_vapour_column is not None:
        mass = mass - ds[water_vapour_column]
    n_dry = mass * (N_A / M_DRY)
    out = ds.copy()
    out[name] = apply_cf_attrs(n_dry, _DRY_AIR_COLUMN, overwrite=True)
    return out

mixing_ratio_to_column

mixing_ratio_to_column(ds: Dataset, *, vmr: str, pressure: str = 'pressure', specific_humidity: str | None = None, level: str = 'level', name: str = 'column') -> Dataset

Column number density of a gas from its dry-air mole-fraction profile.

Trapezoid integral over the pressure levels provided (hydrostatic):

Ω = N_A / (g M_dry) ∫ χ (1 − q) dp

so a uniform profile χ ≡ c over the full column gives Ω = c N_dry. The integral spans the levels present in ds only — include p = 0 and p = p_s levels to cover the whole column. The result is orientation-independent (absolute value).

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset with the mole-fraction profile and pressure per level.	required
vmr	str	Dry-air mole fraction χ per level (dimensionless; ppb in gives ppb-scaled molecules out).	required
pressure	str	Pressure per level [Pa]; may be the level coordinate.	'pressure'
specific_humidity	str | None	Optional specific humidity [kg kg⁻¹] per level; None treats the air as dry.	None
level	str	Vertical dimension.	'level'
name	str	Name of the output variable.	'column'

Returns:

Type	Description
Dataset	ds with name added (molecules m⁻², units="m-2"; no CF
Dataset	standard_name since the CF mole-content name is in mol m⁻²);
Dataset	the level dim is reduced away.

Source code in packages/xrtoolz-ops/src/xrops/atm/gas/ch4/_src/column.py

def mixing_ratio_to_column(
    ds: xr.Dataset,
    *,
    vmr: str,
    pressure: str = "pressure",
    specific_humidity: str | None = None,
    level: str = "level",
    name: str = "column",
) -> xr.Dataset:
    """Column number density of a gas from its dry-air mole-fraction profile.

    Trapezoid integral over the pressure levels provided (hydrostatic):

    ``Ω = N_A / (g M_dry) ∫ χ (1 − q) dp``

    so a uniform profile ``χ ≡ c`` over the full column gives
    ``Ω = c N_dry``. The integral spans the levels present in ``ds`` only —
    include ``p = 0`` and ``p = p_s`` levels to cover the whole column. The
    result is orientation-independent (absolute value).

    Args:
        ds: Dataset with the mole-fraction profile and pressure per level.
        vmr: Dry-air mole fraction ``χ`` per level (dimensionless; ppb in
            gives ppb-scaled molecules out).
        pressure: Pressure per level [Pa]; may be the ``level`` coordinate.
        specific_humidity: Optional specific humidity [kg kg⁻¹] per level;
            ``None`` treats the air as dry.
        level: Vertical dimension.
        name: Name of the output variable.

    Returns:
        ``ds`` with ``name`` added (molecules m⁻², ``units="m-2"``; no CF
        ``standard_name`` since the CF mole-content name is in mol m⁻²);
        the ``level`` dim is reduced away.
    """
    chi = ds[vmr]
    if specific_humidity is not None:
        chi = chi * (1.0 - ds[specific_humidity])
    integral = column_integral(chi, dim=level, coord=ds[pressure], method="trapezoid")
    column = abs(integral) * _MOLECULES_PER_PA
    out = ds.copy()
    out[name] = apply_cf_attrs(column, _CH4_COLUMN, overwrite=True)
    return out

column_mass_to_delta_vmr

column_mass_to_delta_vmr(column_mass: DataArray, *, pressure_pa: float | DataArray, temperature_k: float | DataArray, path_length_m: float | DataArray, molar_mass_kg_mol: float = M_CH4) -> DataArray

Column mass enhancement → volume-mixing-ratio enhancement.

Spreads the column mass m [kg m⁻²] over a single well-mixed layer of length L at the ideal-gas number density n_air = p / (k_B T):

Δχ = (m N_A / M_gas) / (n_air L)

Parameters:

Name	Type	Description	Default
column_mass	DataArray	Column mass enhancement [kg m⁻²].	required
pressure_pa	float | DataArray	Layer pressure [Pa].	required
temperature_k	float | DataArray	Layer temperature [K].	required
path_length_m	float | DataArray	Vertical path length L [m] over which the enhancement is distributed.	required
molar_mass_kg_mol	float	Molar mass of the gas [kg mol⁻¹] (methane).	M_CH4

Returns:

Type	Description
DataArray	Dimensionless mixing-ratio enhancement, same shape as
DataArray	column_mass broadcast against the layer parameters.

Source code in packages/xrtoolz-ops/src/xrops/atm/gas/ch4/_src/column.py

def column_mass_to_delta_vmr(
    column_mass: xr.DataArray,
    *,
    pressure_pa: float | xr.DataArray,
    temperature_k: float | xr.DataArray,
    path_length_m: float | xr.DataArray,
    molar_mass_kg_mol: float = M_CH4,
) -> xr.DataArray:
    """Column mass enhancement → volume-mixing-ratio enhancement.

    Spreads the column mass ``m`` [kg m⁻²] over a single well-mixed layer
    of length ``L`` at the ideal-gas number density ``n_air = p / (k_B T)``:

    ``Δχ = (m N_A / M_gas) / (n_air L)``

    Args:
        column_mass: Column mass enhancement [kg m⁻²].
        pressure_pa: Layer pressure [Pa].
        temperature_k: Layer temperature [K].
        path_length_m: Vertical path length ``L`` [m] over which the
            enhancement is distributed.
        molar_mass_kg_mol: Molar mass of the gas [kg mol⁻¹] (methane).

    Returns:
        Dimensionless mixing-ratio enhancement, same shape as
        ``column_mass`` broadcast against the layer parameters.
    """
    n_air = pressure_pa / (K_B * temperature_k)
    molecules_per_m2 = column_mass / molar_mass_kg_mol * N_A
    delta = molecules_per_m2 / (n_air * path_length_m)
    return delta.assign_attrs(long_name="Mixing-ratio enhancement", units="1")