Calculus

Finite-difference vector calculus on xarray datasets. Each operator works on cartesian, rectilinear, or spherical geometries by deriving the grid metrics from the coordinates. For a scalar \(\phi\) and vector field \(\mathbf{u}=(u,v)\):

\[ \nabla\phi=\left(\frac{\partial\phi}{\partial x},\frac{\partial\phi}{\partial y}\right),\quad \nabla\cdot\mathbf{u}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y},\quad \nabla\times\mathbf{u}=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y},\quad \Delta\phi=\nabla\cdot\nabla\phi. \]

The stencils are geoops's diagnostic difference kernel (geoops.geom.difference, geostack OP-O7), so a field gives the same numbers in xrops.grad and geoops.geom. On a lon/lat grid the metric is a sphere of radius EARTH_RADIUS by default; ellipsoid="WGS84" on partial, gradient, Gradient and Difference uses that ellipsoid's radii, as geoops does on a geographic CRS.

Operators

Gradient

Gradient(dims: Sequence[str], *, geometry: Geometry = 'cartesian', accuracy: int = 1, method: str = 'central', periodic: str | Sequence[str] = (), nan_policy: NanPolicy = 'propagate', lon: str = 'lon', lat: str = 'lat', radius: float = EARTH_RADIUS, ellipsoid: str | None = None)

Bases: Operator

The first partial derivatives along dims, stacked on a component dim.

The result has a new leading component dim labelled "d/d<dim>" in dims order; a Dataset gets one per variable.

Parameters:

Name	Type	Description	Default
dims	Sequence[str]	The dims to differentiate along.	required
geometry	Geometry	As for Difference.	'cartesian'
accuracy	int	The accuracy order of the one-sided stencils.	1
method	str	"central", "forward" or "backward".	'central'
periodic	str | Sequence[str]	The dim, or dims, that wrap (a global longitude).	()
nan_policy	NanPolicy	"propagate" or "adaptive".	'propagate'
lon	str	Longitude coordinate name (spherical).	'lon'
lat	str	Latitude coordinate name (spherical).	'lat'
radius	float	Sphere radius in metres (spherical).	EARTH_RADIUS
ellipsoid	str | None	A pyproj ellipsoid name for an ellipsoidal metric.	None

Example

import numpy as np import xarray as xr y, x = np.arange(3.0), np.arange(4.0) da = xr.DataArray(3.0 * x + 2.0 * y[:, None], dims=("y", "x"), ... coords={"y": y, "x": x}) g = Gradient(dims=("y", "x"))(da) g.dims, g["component"].values.tolist() (('component', 'y', 'x'), ['d/dy', 'd/dx']) float(g.sel(component="d/dx").mean()), float(g.sel(component="d/dy").mean()) (3.0, 2.0)

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

def __init__(
    self,
    dims: Sequence[str],
    *,
    geometry: Geometry = "cartesian",
    accuracy: int = 1,
    method: str = "central",
    periodic: str | Sequence[str] = (),
    nan_policy: NanPolicy = "propagate",
    lon: str = "lon",
    lat: str = "lat",
    radius: float = EARTH_RADIUS,
    ellipsoid: str | None = None,
) -> None:
    if isinstance(dims, str) or not dims:
        raise ValueError("Gradient needs a sequence of one or more dims.")
    self.dims = tuple(dims)
    self.geometry = geometry
    self.accuracy = accuracy
    self.method = method
    self.periodic = (periodic,) if isinstance(periodic, str) else tuple(periodic)
    self.nan_policy = nan_policy
    self.lon = lon
    self.lat = lat
    self.radius = radius
    self.ellipsoid = ellipsoid
    _geom_kw(geometry, lon, lat, radius, ellipsoid)

Difference

Difference(dim: str, *, derivative: int = 1, geometry: Geometry = 'cartesian', accuracy: int = 1, method: str = 'central', periodic: bool = False, nan_policy: NanPolicy = 'propagate', lon: str = 'lon', lat: str = 'lat', radius: float = EARTH_RADIUS, ellipsoid: str | None = None)

Bases: Operator

The derivative-th partial derivative along one dim.

Parameters:

Name	Type	Description	Default
dim	str	The dim to differentiate along.	required
derivative	int	The derivative order.	1
geometry	Geometry	"cartesian" (uniform spacing), "rectilinear" (non-uniform 1-D coordinate) or "spherical" (lon/lat in degrees, derivatives per metre).	'cartesian'
accuracy	int	The accuracy order of the one-sided stencils.	1
method	str	"central", "forward" or "backward".	'central'
periodic	bool	The dim wraps (a global longitude).	False
nan_policy	NanPolicy	"propagate" (a NaN costs a halo) or "adaptive" (per point, the widest stencil on finite data).	'propagate'
lon	str	Longitude coordinate name (spherical).	'lon'
lat	str	Latitude coordinate name (spherical).	'lat'
radius	float	Sphere radius in metres (spherical, the default metric).	EARTH_RADIUS
ellipsoid	str | None	A pyproj ellipsoid name for an ellipsoidal metric (spherical), as geoops's on a geographic CRS.	None

Example

import numpy as np import xarray as xr x = np.linspace(0.0, 3.0, 4) da = xr.DataArray(x**2, dims="x", coords={"x": x}) Difference("x")(da).values array([1., 2., 4., 5.])

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

def __init__(
    self,
    dim: str,
    *,
    derivative: int = 1,
    geometry: Geometry = "cartesian",
    accuracy: int = 1,
    method: str = "central",
    periodic: bool = False,
    nan_policy: NanPolicy = "propagate",
    lon: str = "lon",
    lat: str = "lat",
    radius: float = EARTH_RADIUS,
    ellipsoid: str | None = None,
) -> None:
    self.dim = dim
    self.derivative = derivative
    self.geometry = geometry
    self.accuracy = accuracy
    self.method = method
    self.periodic = periodic
    self.nan_policy = nan_policy
    self.lon = lon
    self.lat = lat
    self.radius = radius
    self.ellipsoid = ellipsoid
    _geom_kw(geometry, lon, lat, radius, ellipsoid)

Differential operators

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

partial

partial(da: DataArray, dim: str, *, order: int = 1, geometry: Geometry = 'cartesian', accuracy: int = 1, method: str = 'central', periodic: bool = False, nan_policy: NanPolicy = 'propagate', **geom_kw: Any) -> DataArray

Partial derivative ∂ᵒʳᵈᵉʳda/∂<dim>ᵒʳᵈᵉʳ under the given geometry.

Parameters:

Name	Type	Description	Default
da	DataArray	Input field.	required
dim	str	Dimension along which to differentiate.	required
order	int	Derivative order. 1 (default) is the first derivative. order > 1 is supported for geometry="cartesian" (and for "rectilinear" when the coordinate turns out to be uniformly spaced, which delegates to the cartesian backend); see Raises for the other geometries.	1
geometry	Geometry	One of "cartesian" (uniform grid), "rectilinear" (non-uniform 1-D coords), "spherical" (lon/lat in degrees).	'cartesian'
accuracy	int	Accuracy order forwarded to :func:geoops.geom.difference.	1
method	str	"central" (default), "forward", or "backward".	'central'
periodic	bool	Treat dim as wrapping, so the first and last points are differentiated against each other rather than falling back to one-sided stencils. For geometry="spherical" this is valid on the longitude axis of a full 360° grid only.	False
nan_policy	NanPolicy	How to treat NaN (typically a land mask). "propagate" (default) is plain stencil arithmetic, so a masked cell costs a halo of roughly accuracy valid cells around it. "adaptive" instead picks, per point, the widest stencil that lies wholly on finite data, recovering those cells at one order lower accuracy; masked cells stay NaN. Raising accuracy improves the recovered cells.	'propagate'
**geom_kw	Any	Geometry-specific keyword arguments (e.g. radius and coord names for spherical, uniform_rtol for cartesian).	{}

Returns:

Type	Description
DataArray	DataArray with the same dims/coords as da.

Raises:

Type	Description
ValueError	for an unknown geometry, a non-positive order, or periodic on a spherical grid that is not global.
NotImplementedError	for order > 1 under geometry="spherical", or under geometry="rectilinear" with a genuinely non-uniform coordinate; likewise for periodic on a latitude axis or a non-uniform coordinate.

Example

d(x²)/dx = 2x — exact on interior points with the central stencil; the boundary points fall back to one-sided differences of the requested accuracy:

>>> import numpy as np
>>> import xarray as xr
>>> x = np.linspace(0.0, 3.0, 4)
>>> da = xr.DataArray(x**2, dims="x", coords={"x": x})
>>> partial(da, "x").values
array([1., 2., 4., 5.])

The second derivative d²(x²)/dx² = 2 comes from a single order-2 stencil rather than two composed passes:

>>> partial(da, "x", order=2).values
array([2., 2., 2., 2.])

A masked cell costs its neighbour too, because the centred stencil there reads the NaN. nan_policy="adaptive" gives that neighbour a one-sided estimate instead, while the masked cell itself stays NaN:

>>> v = np.linspace(0.0, 5.0, 6) ** 2
>>> v[0] = np.nan
>>> masked = xr.DataArray(v, dims="x", coords={"x": np.linspace(0.0, 5.0, 6)})
>>> partial(masked, "x").values
array([nan, nan,  4.,  6.,  8.,  9.])
>>> partial(masked, "x", nan_policy="adaptive").values
array([nan,  3.,  4.,  6.,  8.,  9.])

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

def partial(
    da: xr.DataArray,
    dim: str,
    *,
    order: int = 1,
    geometry: Geometry = "cartesian",
    accuracy: int = 1,
    method: str = "central",
    periodic: bool = False,
    nan_policy: NanPolicy = "propagate",
    **geom_kw: Any,
) -> xr.DataArray:
    """Partial derivative ``∂ᵒʳᵈᵉʳda/∂<dim>ᵒʳᵈᵉʳ`` under the given geometry.

    Args:
        da: Input field.
        dim: Dimension along which to differentiate.
        order: Derivative order. ``1`` (default) is the first derivative.
            ``order > 1`` is supported for ``geometry="cartesian"`` (and
            for ``"rectilinear"`` when the coordinate turns out to be
            uniformly spaced, which delegates to the cartesian backend);
            see Raises for the other geometries.
        geometry: One of ``"cartesian"`` (uniform grid), ``"rectilinear"``
            (non-uniform 1-D coords), ``"spherical"`` (lon/lat in degrees).
        accuracy: Accuracy order forwarded to :func:`geoops.geom.difference`.
        method: ``"central"`` (default), ``"forward"``, or ``"backward"``.
        periodic: Treat ``dim`` as wrapping, so the first and last points
            are differentiated against each other rather than falling
            back to one-sided stencils. For ``geometry="spherical"`` this
            is valid on the longitude axis of a full 360° grid only.
        nan_policy: How to treat NaN (typically a land mask).
            ``"propagate"`` (default) is plain stencil arithmetic, so a
            masked cell costs a halo of roughly ``accuracy`` valid cells
            around it. ``"adaptive"`` instead picks, per point, the
            widest stencil that lies wholly on finite data, recovering
            those cells at one order lower accuracy; masked cells stay
            NaN. Raising ``accuracy`` improves the recovered cells.
        **geom_kw: Geometry-specific keyword arguments (e.g. ``radius``
            and coord names for spherical, ``uniform_rtol`` for cartesian).

    Returns:
        DataArray with the same dims/coords as ``da``.

    Raises:
        ValueError: for an unknown ``geometry``, a non-positive ``order``,
            or ``periodic`` on a spherical grid that is not global.
        NotImplementedError: for ``order > 1`` under
            ``geometry="spherical"``, or under ``geometry="rectilinear"``
            with a genuinely non-uniform coordinate; likewise for
            ``periodic`` on a latitude axis or a non-uniform coordinate.

    Example:
        ``d(x²)/dx = 2x`` — exact on interior points with the central
        stencil; the boundary points fall back to one-sided differences
        of the requested accuracy:

        ```pycon
        >>> import numpy as np
        >>> import xarray as xr
        >>> x = np.linspace(0.0, 3.0, 4)
        >>> da = xr.DataArray(x**2, dims="x", coords={"x": x})
        >>> partial(da, "x").values
        array([1., 2., 4., 5.])

        ```

        The second derivative ``d²(x²)/dx² = 2`` comes from a single
        order-2 stencil rather than two composed passes:

        ```pycon
        >>> partial(da, "x", order=2).values
        array([2., 2., 2., 2.])

        ```

        A masked cell costs its neighbour too, because the centred
        stencil there reads the NaN. ``nan_policy="adaptive"`` gives that
        neighbour a one-sided estimate instead, while the masked cell
        itself stays NaN:

        ```pycon
        >>> v = np.linspace(0.0, 5.0, 6) ** 2
        >>> v[0] = np.nan
        >>> masked = xr.DataArray(v, dims="x", coords={"x": np.linspace(0.0, 5.0, 6)})
        >>> partial(masked, "x").values
        array([nan, nan,  4.,  6.,  8.,  9.])
        >>> partial(masked, "x", nan_policy="adaptive").values
        array([nan,  3.,  4.,  6.,  8.,  9.])

        ```
    """
    kw: dict[str, Any] = dict(
        order=order,
        accuracy=accuracy,
        method=method,
        periodic=periodic,
        nan_policy=nan_policy,
        **geom_kw,
    )
    if geometry == "cartesian":
        return cartesian.cartesian_partial(da, dim, **kw)
    if geometry == "rectilinear":
        return rectilinear.rectilinear_partial(da, dim, **kw)
    if geometry == "spherical":
        return spherical.spherical_partial(da, dim, **kw)
    raise ValueError(
        f"Unknown geometry {geometry!r}; expected one of "
        "'cartesian', 'rectilinear', 'spherical'."
    )

gradient

gradient(da: DataArray, *, dims: tuple[str, ...] | None = None, geometry: Geometry = 'cartesian', accuracy: int | tuple[int, ...] = 1, method: str = 'central', periodic: Periodic = None, nan_policy: NanPolicy = 'propagate', **geom_kw: Any) -> Dataset

Gradient ∇da under the given geometry.

Parameters:

Name	Type	Description	Default
da	DataArray	Input scalar field.	required
dims	tuple[str, ...] | None	Dimensions to differentiate against (defaults to da.dims).	None
geometry	Geometry	"cartesian" | "rectilinear" | "spherical".	'cartesian'
accuracy	int | tuple[int, ...]	Scalar or per-dim tuple.	1
method	str	Forwarded to :func:geoops.geom.difference.	'central'
periodic	Periodic	Dimension name, or sequence of names, to treat as wrapping. Must be among the differentiated dims.	None
nan_policy	NanPolicy	"propagate" or "adaptive"; see :func:partial.	'propagate'
**geom_kw	Any	Geometry-specific keyword arguments.	{}

Returns:

Type	Description
Dataset	Dataset with one DataArray per dim in dims, named
Dataset	d<name>_d<axis> (f when the input is unnamed).

Example

>>> import numpy as np
>>> import xarray as xr
>>> x = np.linspace(0.0, 3.0, 4)
>>> X, Y = np.meshgrid(x, x, indexing="xy")
>>> f = xr.DataArray(
...     X**2 + Y**2, dims=("y", "x"), coords={"y": x, "x": x},
... )
>>> grads = gradient(f, dims=("x", "y"))
>>> sorted(grads.data_vars)
['df_dx', 'df_dy']

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

def gradient(
    da: xr.DataArray,
    *,
    dims: tuple[str, ...] | None = None,
    geometry: Geometry = "cartesian",
    accuracy: int | tuple[int, ...] = 1,
    method: str = "central",
    periodic: Periodic = None,
    nan_policy: NanPolicy = "propagate",
    **geom_kw: Any,
) -> xr.Dataset:
    """Gradient ``∇da`` under the given geometry.

    Args:
        da: Input scalar field.
        dims: Dimensions to differentiate against (defaults to ``da.dims``).
        geometry: ``"cartesian"`` | ``"rectilinear"`` | ``"spherical"``.
        accuracy: Scalar or per-dim tuple.
        method: Forwarded to :func:`geoops.geom.difference`.
        periodic: Dimension name, or sequence of names, to treat as
            wrapping. Must be among the differentiated dims.
        nan_policy: ``"propagate"`` or ``"adaptive"``; see
            :func:`partial`.
        **geom_kw: Geometry-specific keyword arguments.

    Returns:
        Dataset with one DataArray per dim in ``dims``, named
        ``d<name>_d<axis>`` (``f`` when the input is unnamed).

    Example:
        ```pycon
        >>> import numpy as np
        >>> import xarray as xr
        >>> x = np.linspace(0.0, 3.0, 4)
        >>> X, Y = np.meshgrid(x, x, indexing="xy")
        >>> f = xr.DataArray(
        ...     X**2 + Y**2, dims=("y", "x"), coords={"y": x, "x": x},
        ... )
        >>> grads = gradient(f, dims=("x", "y"))
        >>> sorted(grads.data_vars)
        ['df_dx', 'df_dy']

        ```
    """
    if geometry == "spherical":
        target = (
            (geom_kw.get("lon", "lon"), geom_kw.get("lat", "lat"))
            if dims is None
            else tuple(dims)
        )
    else:
        target = tuple(da.dims) if dims is None else tuple(dims)
    wrap = _periodic_set(periodic, target)

    if geometry == "cartesian":
        return cartesian.cartesian_gradient(
            da,
            dims=dims,
            accuracy=accuracy,
            method=method,
            periodic=wrap,
            nan_policy=nan_policy,
            **geom_kw,
        )
    if geometry == "rectilinear":
        return rectilinear.rectilinear_gradient(
            da,
            dims=dims,
            accuracy=accuracy,
            method=method,
            periodic=wrap,
            nan_policy=nan_policy,
            **geom_kw,
        )
    if geometry == "spherical":
        return spherical.spherical_gradient(
            da,
            dims=dims,
            accuracy=accuracy,
            method=method,
            periodic=wrap,
            nan_policy=nan_policy,
            **geom_kw,
        )
    raise ValueError(
        f"Unknown geometry {geometry!r}; expected one of "
        "'cartesian', 'rectilinear', 'spherical'."
    )

divergence

divergence(ds: Dataset, components: tuple[str, str] | tuple[str, str, str], *, dims: tuple[str, ...], geometry: Geometry | Sequence[Geometry] = 'cartesian', accuracy: int | tuple[int, ...] = 1, method: str = 'central', periodic: Periodic = None, nan_policy: NanPolicy = 'propagate', **geom_kw: Any) -> DataArray

Divergence ∇·F of a vector field stored as Dataset variables.

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset containing the vector components.	required
components	tuple[str, str] | tuple[str, str, str]	Variable names for (F_x, F_y, …). Must match the length of dims.	required
dims	tuple[str, ...]	Spatial dimensions that pair with each component (e.g. ("x", "y") for cartesian, ("lon", "lat") for spherical).	required
geometry	Geometry | Sequence[Geometry]	A single geometry applied to every dim, or a per-dim sequence pairing with dims. The per-dim form supports mixed grids — the common case being a spherical horizontal pair plus a rectilinear vertical axis, e.g. ("spherical", "spherical", "rectilinear") with dims=("lon", "lat", "depth").	'cartesian'
accuracy	int | tuple[int, ...]	Scalar (applied to every dim) or per-dim tuple matching the length of dims.	1
method	str	"central" | "forward" | "backward".	'central'
periodic	Periodic	Dimension name, or sequence of names, to treat as wrapping. Must be among dims.	None
nan_policy	NanPolicy	"propagate" or "adaptive"; see :func:partial. Applies to each component's stencil; the spherical curvature term is pointwise and inherits the mask.	'propagate'
**geom_kw	Any	Geometry-specific kwargs (radius, lon, lat, uniform_rtol). Each is routed only to the backends that accept it, so radius may accompany a mixed geometry without reaching the cartesian or rectilinear axes.	{}

Returns:

Type	Description
DataArray	Scalar DataArray with the same dims/coords as the components.

Raises:

Type	Description
ValueError	if components and dims lengths differ, if a geometry sequence does not match dims, or if exactly one axis (or more than two) is marked spherical — the metric needs the lon/lat pair together.

Notes

The curvature correction -(v tan φ) / R is added whenever both the longitude and latitude axes are spherical, so the result matches :func:metpy.calc.divergence on lon/lat fields:

.. math::

\nabla \cdot \mathbf{F} =
\frac{1}{R\cos\varphi} \frac{\partial F_x}{\partial \lambda}
+ \frac{1}{R} \frac{\partial F_y}{\partial \varphi}
- \frac{F_y \tan\varphi}{R}

Example

The radial field F = (x, y) has constant divergence 2:

>>> import numpy as np
>>> import xarray as xr
>>> x = np.linspace(0.0, 3.0, 4)
>>> X, Y = np.meshgrid(x, x, indexing="xy")
>>> coords = {"y": x, "x": x}
>>> ds = xr.Dataset(
...     {
...         "u": xr.DataArray(X, dims=("y", "x"), coords=coords),
...         "v": xr.DataArray(Y, dims=("y", "x"), coords=coords),
...     }
... )
>>> np.unique(divergence(ds, ("u", "v"), dims=("x", "y")).values)
array([2.])

A 3-D ocean budget mixes a spherical horizontal pair with a rectilinear (stretched) vertical axis in one call:

>>> depth = np.array([0.0, 10.0, 30.0, 70.0])
>>> lon, lat = np.array([0.0, 1.0, 2.0]), np.array([10.0, 11.0, 12.0])
>>> shape = (depth.size, lat.size, lon.size)
>>> dims3 = ("depth", "lat", "lon")
>>> coords3 = {"depth": depth, "lat": lat, "lon": lon}
>>> flow = xr.Dataset(
...     {
...         name: xr.DataArray(np.zeros(shape), dims=dims3, coords=coords3)
...         for name in ("u", "v", "w")
...     }
... )
>>> div = divergence(
...     flow,
...     ("u", "v", "w"),
...     dims=("lon", "lat", "depth"),
...     geometry=("spherical", "spherical", "rectilinear"),
... )
>>> div.dims
('depth', 'lat', 'lon')

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

def divergence(
    ds: xr.Dataset,
    components: tuple[str, str] | tuple[str, str, str],
    *,
    dims: tuple[str, ...],
    geometry: Geometry | Sequence[Geometry] = "cartesian",
    accuracy: int | tuple[int, ...] = 1,
    method: str = "central",
    periodic: Periodic = None,
    nan_policy: NanPolicy = "propagate",
    **geom_kw: Any,
) -> xr.DataArray:
    r"""Divergence ``∇·F`` of a vector field stored as Dataset variables.

    Args:
        ds: Dataset containing the vector components.
        components: Variable names for ``(F_x, F_y, …)``. Must match the
            length of ``dims``.
        dims: Spatial dimensions that pair with each component (e.g.
            ``("x", "y")`` for cartesian, ``("lon", "lat")`` for spherical).
        geometry: A single geometry applied to every dim, or a per-dim
            sequence pairing with ``dims``. The per-dim form supports
            mixed grids — the common case being a spherical horizontal
            pair plus a rectilinear vertical axis, e.g.
            ``("spherical", "spherical", "rectilinear")`` with
            ``dims=("lon", "lat", "depth")``.
        accuracy: Scalar (applied to every dim) or per-dim tuple matching
            the length of ``dims``.
        method: ``"central"`` | ``"forward"`` | ``"backward"``.
        periodic: Dimension name, or sequence of names, to treat as
            wrapping. Must be among ``dims``.
        nan_policy: ``"propagate"`` or ``"adaptive"``; see
            :func:`partial`. Applies to each component's stencil; the
            spherical curvature term is pointwise and inherits the mask.
        **geom_kw: Geometry-specific kwargs (``radius``, ``lon``, ``lat``,
            ``uniform_rtol``). Each is routed only to the backends that
            accept it, so ``radius`` may accompany a mixed geometry
            without reaching the cartesian or rectilinear axes.

    Returns:
        Scalar DataArray with the same dims/coords as the components.

    Raises:
        ValueError: if ``components`` and ``dims`` lengths differ, if a
            geometry sequence does not match ``dims``, or if exactly one
            axis (or more than two) is marked spherical — the metric
            needs the lon/lat pair together.

    Notes:
        The curvature correction ``-(v tan φ) / R`` is added whenever both
        the longitude and latitude axes are spherical, so the result
        matches :func:`metpy.calc.divergence` on lon/lat fields:

        .. math::

            \nabla \cdot \mathbf{F} =
            \frac{1}{R\cos\varphi} \frac{\partial F_x}{\partial \lambda}
            + \frac{1}{R} \frac{\partial F_y}{\partial \varphi}
            - \frac{F_y \tan\varphi}{R}

    Example:
        The radial field ``F = (x, y)`` has constant divergence 2:

        ```pycon
        >>> import numpy as np
        >>> import xarray as xr
        >>> x = np.linspace(0.0, 3.0, 4)
        >>> X, Y = np.meshgrid(x, x, indexing="xy")
        >>> coords = {"y": x, "x": x}
        >>> ds = xr.Dataset(
        ...     {
        ...         "u": xr.DataArray(X, dims=("y", "x"), coords=coords),
        ...         "v": xr.DataArray(Y, dims=("y", "x"), coords=coords),
        ...     }
        ... )
        >>> np.unique(divergence(ds, ("u", "v"), dims=("x", "y")).values)
        array([2.])

        ```

        A 3-D ocean budget mixes a spherical horizontal pair with a
        rectilinear (stretched) vertical axis in one call:

        ```pycon
        >>> depth = np.array([0.0, 10.0, 30.0, 70.0])
        >>> lon, lat = np.array([0.0, 1.0, 2.0]), np.array([10.0, 11.0, 12.0])
        >>> shape = (depth.size, lat.size, lon.size)
        >>> dims3 = ("depth", "lat", "lon")
        >>> coords3 = {"depth": depth, "lat": lat, "lon": lon}
        >>> flow = xr.Dataset(
        ...     {
        ...         name: xr.DataArray(np.zeros(shape), dims=dims3, coords=coords3)
        ...         for name in ("u", "v", "w")
        ...     }
        ... )
        >>> div = divergence(
        ...     flow,
        ...     ("u", "v", "w"),
        ...     dims=("lon", "lat", "depth"),
        ...     geometry=("spherical", "spherical", "rectilinear"),
        ... )
        >>> div.dims
        ('depth', 'lat', 'lon')

        ```
    """
    _sphere_only(geom_kw, "divergence")
    if len(components) != len(dims):
        raise ValueError(
            f"components ({len(components)}) and dims ({len(dims)}) "
            "must have the same length."
        )
    geometries = _resolve_geometries(geometry, dims)
    _validate_geom_kw(geometries, geom_kw)
    per_dim = cartesian._per_dim_accuracy(accuracy, dims)
    wrap = _periodic_set(periodic, dims)

    lon = geom_kw.get("lon", "lon")
    lat = geom_kw.get("lat", "lat")
    spherical_dims = [
        d for d, g in zip(dims, geometries, strict=True) if g == "spherical"
    ]
    # Validated up front so the geometry-level message wins over whichever
    # per-axis partial would otherwise fail first.
    if spherical_dims and sorted(spherical_dims) != sorted({lon, lat}):
        if all(g == "spherical" for g in geometries):
            raise ValueError(
                "spherical divergence is defined for 2-D (lon, lat) "
                f"fields; got dims={tuple(dims)!r}. For a 3-D field pass "
                "a per-dim geometry sequence marking the vertical axis "
                "'rectilinear' or 'cartesian', e.g. geometry="
                "('spherical', 'spherical', 'rectilinear')."
            )
        raise ValueError(
            f"spherical axes {spherical_dims!r} must be exactly the "
            f"longitude/latitude pair ({lon!r}, {lat!r}): the metric "
            "term -(v tan φ)/R couples the two, so neither can be "
            "spherical on its own. Mark the remaining axes "
            "'rectilinear' or 'cartesian'."
        )

    total: xr.DataArray | None = None
    for comp_name, dim, geom, acc in zip(
        components, dims, geometries, per_dim, strict=True
    ):
        d = partial(
            ds[comp_name],
            dim,
            geometry=geom,
            accuracy=acc,
            method=method,
            periodic=dim in wrap,
            nan_policy=nan_policy,
            **_geom_kw_for(geom, geom_kw),
        )
        total = d if total is None else total + d
    assert total is not None  # narrows Optional; only empty components trips this

    if spherical_dims:
        radius = geom_kw.get("radius", EARTH_RADIUS)
        v = ds[components[dims.index(lat)]]
        phi = np.deg2rad(np.asarray(v[lat].values))
        tan_phi = _broadcast_lat_factor(np.tan(phi), ref=v, lat=lat)
        curvature = -(v.data * tan_phi) / radius  # lazy on dask
        total = total + xr.DataArray(
            curvature,
            dims=v.dims,
            coords={k: v.coords[k] for k in v.coords},
        )

    return total.rename(None)

curl

curl(ds: Dataset, components: tuple[str, str], *, dims: tuple[str, str], geometry: Geometry = 'cartesian', accuracy: int | tuple[int, ...] = 1, method: str = 'central', periodic: Periodic = None, nan_policy: NanPolicy = 'propagate', **geom_kw: Any) -> DataArray

2-D scalar curl ∂v/∂x − ∂u/∂y (vertical component of ∇×F).

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset with the two horizontal components.	required
components	tuple[str, str]	(u_name, v_name) — eastward then northward.	required
dims	tuple[str, str]	(x_dim, y_dim) paired with the components.	required
geometry	Geometry	"cartesian" | "rectilinear" | "spherical".	'cartesian'
accuracy	int | tuple[int, ...]	Scalar (applied to both dims) or a 2-tuple pairing with dims, forwarded to :func:geoops.geom.difference.	1
periodic	Periodic	Dimension name, or sequence of names, to treat as wrapping. Must be among dims.	None
nan_policy	NanPolicy	"propagate" or "adaptive"; see :func:partial.	'propagate'
method	str	Stencil family (e.g. "central"), forwarded to :func:geoops.geom.difference.	'central'
**geom_kw	Any	Forwarded to :func:partial.	{}

Returns:

Type	Description
DataArray	Scalar DataArray with the same dims/coords as the components.

Notes

For spherical geometry the curvature correction +(u tan φ) / R is added so the result matches :func:metpy.calc.vorticity on lon/lat fields:

.. math::

\zeta =
\frac{1}{R\cos\varphi} \frac{\partial v}{\partial \lambda}
- \frac{1}{R} \frac{\partial u}{\partial \varphi}
+ \frac{u \tan\varphi}{R}

Example

Rigid rotation F = (-y, x) has constant vorticity 2:

>>> import numpy as np
>>> import xarray as xr
>>> x = np.linspace(0.0, 3.0, 4)
>>> X, Y = np.meshgrid(x, x, indexing="xy")
>>> coords = {"y": x, "x": x}
>>> ds = xr.Dataset(
...     {
...         "u": xr.DataArray(-Y, dims=("y", "x"), coords=coords),
...         "v": xr.DataArray(X, dims=("y", "x"), coords=coords),
...     }
... )
>>> np.unique(curl(ds, ("u", "v"), dims=("x", "y")).values)
array([2.])

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

def curl(
    ds: xr.Dataset,
    components: tuple[str, str],
    *,
    dims: tuple[str, str],
    geometry: Geometry = "cartesian",
    accuracy: int | tuple[int, ...] = 1,
    method: str = "central",
    periodic: Periodic = None,
    nan_policy: NanPolicy = "propagate",
    **geom_kw: Any,
) -> xr.DataArray:
    r"""2-D scalar curl ``∂v/∂x − ∂u/∂y`` (vertical component of ``∇×F``).

    Args:
        ds: Dataset with the two horizontal components.
        components: ``(u_name, v_name)`` — eastward then northward.
        dims: ``(x_dim, y_dim)`` paired with the components.
        geometry: ``"cartesian"`` | ``"rectilinear"`` | ``"spherical"``.
        accuracy: Scalar (applied to both dims) or a 2-tuple pairing with
            ``dims``, forwarded to :func:`geoops.geom.difference`.
        periodic: Dimension name, or sequence of names, to treat as
            wrapping. Must be among ``dims``.
        nan_policy: ``"propagate"`` or ``"adaptive"``; see
            :func:`partial`.
        method: Stencil family (e.g. ``"central"``), forwarded to
            :func:`geoops.geom.difference`.
        **geom_kw: Forwarded to :func:`partial`.

    Returns:
        Scalar DataArray with the same dims/coords as the components.

    Notes:
        For spherical geometry the curvature correction
        ``+(u tan φ) / R`` is added so the result matches
        :func:`metpy.calc.vorticity` on lon/lat fields:

        .. math::

            \zeta =
            \frac{1}{R\cos\varphi} \frac{\partial v}{\partial \lambda}
            - \frac{1}{R} \frac{\partial u}{\partial \varphi}
            + \frac{u \tan\varphi}{R}

    Example:
        Rigid rotation ``F = (-y, x)`` has constant vorticity 2:

        ```pycon
        >>> import numpy as np
        >>> import xarray as xr
        >>> x = np.linspace(0.0, 3.0, 4)
        >>> X, Y = np.meshgrid(x, x, indexing="xy")
        >>> coords = {"y": x, "x": x}
        >>> ds = xr.Dataset(
        ...     {
        ...         "u": xr.DataArray(-Y, dims=("y", "x"), coords=coords),
        ...         "v": xr.DataArray(X, dims=("y", "x"), coords=coords),
        ...     }
        ... )
        >>> np.unique(curl(ds, ("u", "v"), dims=("x", "y")).values)
        array([2.])

        ```
    """
    _sphere_only(geom_kw, "curl")
    if len(components) != 2 or len(dims) != 2:
        raise ValueError("2-D curl needs exactly two components and two dims.")
    u_name, v_name = components
    x_dim, y_dim = dims
    acc_x, acc_y = cartesian._per_dim_accuracy(accuracy, dims)
    wrap = _periodic_set(periodic, dims)
    dvdx = partial(
        ds[v_name],
        x_dim,
        geometry=geometry,
        accuracy=acc_x,
        method=method,
        periodic=x_dim in wrap,
        nan_policy=nan_policy,
        **geom_kw,
    )
    dudy = partial(
        ds[u_name],
        y_dim,
        geometry=geometry,
        accuracy=acc_y,
        method=method,
        periodic=y_dim in wrap,
        nan_policy=nan_policy,
        **geom_kw,
    )
    out = dvdx - dudy

    if geometry == "spherical":
        lat = geom_kw.get("lat", "lat")
        radius = geom_kw.get("radius", EARTH_RADIUS)
        u = ds[u_name]
        phi = np.deg2rad(np.asarray(u[lat].values))
        tan_phi = _broadcast_lat_factor(np.tan(phi), ref=u, lat=lat)
        curvature = (u.data * tan_phi) / radius  # lazy on dask
        out = out + xr.DataArray(
            curvature,
            dims=u.dims,
            coords={k: u.coords[k] for k in u.coords},
        )

    return out.rename(None)

laplacian

laplacian(da: DataArray, *, dims: tuple[str, ...] | None = None, geometry: Geometry = 'cartesian', accuracy: int | tuple[int, ...] = 1, method: str = 'central', periodic: Periodic = None, nan_policy: NanPolicy = 'propagate', **geom_kw: Any) -> DataArray

Laplacian Δf = ∇·∇f.

For geometry="cartesian" this sums direct second-derivative stencils (Σ ∂²f/∂dim²), which is exact for quadratic fields at the default accuracy=1 and contaminates only a single boundary ring. An axis too short to host a second-derivative stencil falls back to the composed form below, so short axes keep working.

The other geometries compose gradient followed by divergence so the spherical curvature correction is inherited automatically. That applies to all rectilinear grids, including those whose coordinate is uniform enough to delegate each partial to the cartesian backend. On those paths the second derivative is two stacked first-derivative stencils, so pass accuracy=2 (or higher) when you need the interior to reproduce polynomial fields exactly.

Parameters:

Name	Type	Description	Default
da	DataArray	Input scalar field.	required
dims	tuple[str, ...] | None	Dimensions to differentiate against. Defaults to da.dims ((lon, lat) for spherical geometry).	None
geometry	Geometry	"cartesian" | "rectilinear" | "spherical".	'cartesian'
accuracy	int | tuple[int, ...]	Scalar (applied to every dim) or per-dim tuple matching the length of dims.	1
method	str	Forwarded to :func:geoops.geom.difference.	'central'
periodic	Periodic	Dimension name, or sequence of names, to treat as wrapping. Must be among the differentiated dims.	None
nan_policy	NanPolicy	"propagate" or "adaptive"; see :func:partial. On the composed (non-cartesian) path it is applied to both the gradient and the divergence stage, so a masked cell costs one halo rather than two.	'propagate'
**geom_kw	Any	Geometry-specific keyword arguments.	{}

Returns:

Type	Description
DataArray	Unnamed scalar DataArray with the same dims/coords as da.

Raises:

Type	Description
ValueError	if dims resolves to an empty tuple, or for an unknown geometry.

Example

Δ(x² + y²) = 4, exact on interior points:

>>> import numpy as np
>>> import xarray as xr
>>> x = np.linspace(0.0, 3.0, 4)
>>> X, Y = np.meshgrid(x, x, indexing="xy")
>>> f = xr.DataArray(
...     X**2 + Y**2, dims=("y", "x"), coords={"y": x, "x": x},
... )
>>> lap = laplacian(f, dims=("x", "y"))
>>> np.unique(lap.values[1:-1, 1:-1])
array([4.])

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

def laplacian(
    da: xr.DataArray,
    *,
    dims: tuple[str, ...] | None = None,
    geometry: Geometry = "cartesian",
    accuracy: int | tuple[int, ...] = 1,
    method: str = "central",
    periodic: Periodic = None,
    nan_policy: NanPolicy = "propagate",
    **geom_kw: Any,
) -> xr.DataArray:
    """Laplacian ``Δf = ∇·∇f``.

    For ``geometry="cartesian"`` this sums direct second-derivative
    stencils (``Σ ∂²f/∂dim²``), which is exact for quadratic fields at
    the default ``accuracy=1`` and contaminates only a single boundary
    ring. An axis too short to host a second-derivative stencil falls
    back to the composed form below, so short axes keep working.

    The other geometries compose gradient followed by divergence so the
    spherical curvature correction is inherited automatically. That
    applies to *all* rectilinear grids, including those whose coordinate
    is uniform enough to delegate each partial to the cartesian backend.
    On those paths the second derivative is two stacked first-derivative
    stencils, so pass ``accuracy=2`` (or higher) when you need the
    interior to reproduce polynomial fields exactly.

    Args:
        da: Input scalar field.
        dims: Dimensions to differentiate against. Defaults to ``da.dims``
            (``(lon, lat)`` for spherical geometry).
        geometry: ``"cartesian"`` | ``"rectilinear"`` | ``"spherical"``.
        accuracy: Scalar (applied to every dim) or per-dim tuple matching
            the length of ``dims``.
        method: Forwarded to :func:`geoops.geom.difference`.
        periodic: Dimension name, or sequence of names, to treat as
            wrapping. Must be among the differentiated dims.
        nan_policy: ``"propagate"`` or ``"adaptive"``; see
            :func:`partial`. On the composed (non-cartesian) path it is
            applied to both the gradient and the divergence stage, so a
            masked cell costs one halo rather than two.
        **geom_kw: Geometry-specific keyword arguments.

    Returns:
        Unnamed scalar DataArray with the same dims/coords as ``da``.

    Raises:
        ValueError: if ``dims`` resolves to an empty tuple, or for an
            unknown ``geometry``.

    Example:
        ``Δ(x² + y²) = 4``, exact on interior points:

        ```pycon
        >>> import numpy as np
        >>> import xarray as xr
        >>> x = np.linspace(0.0, 3.0, 4)
        >>> X, Y = np.meshgrid(x, x, indexing="xy")
        >>> f = xr.DataArray(
        ...     X**2 + Y**2, dims=("y", "x"), coords={"y": x, "x": x},
        ... )
        >>> lap = laplacian(f, dims=("x", "y"))
        >>> np.unique(lap.values[1:-1, 1:-1])
        array([4.])

        ```
    """
    _sphere_only(geom_kw, "laplacian")
    if dims is None:
        if geometry == "spherical":
            target_dims: tuple[str, ...] = (
                geom_kw.get("lon", "lon"),
                geom_kw.get("lat", "lat"),
            )
        else:
            target_dims = tuple(da.dims)
    else:
        target_dims = tuple(dims)

    if not target_dims:
        raise ValueError("laplacian needs at least one dimension to sum over.")

    per_dim = cartesian._per_dim_accuracy(accuracy, target_dims)
    wrap = _periodic_set(periodic, target_dims)

    if geometry == "cartesian":
        total: xr.DataArray | None = None
        for dim, acc in zip(target_dims, per_dim, strict=True):
            wrapped = dim in wrap
            try:
                second = cartesian.cartesian_partial(
                    da,
                    dim,
                    order=2,
                    accuracy=acc,
                    method=method,
                    periodic=wrapped,
                    nan_policy=nan_policy,
                    **geom_kw,
                )
            except ValueError:
                # The stencil kernel rejects a direct second-derivative stencil when
                # the axis is too short to host one (e.g. two samples). Fall
                # back to the composed two-pass form, which is what this
                # function used before the order=2 fast path existed, so short
                # axes keep working. A ValueError raised for any other reason
                # — missing, non-1-D, or non-uniform coordinate — is raised
                # again by these same calls.
                first = cartesian.cartesian_partial(
                    da,
                    dim,
                    accuracy=acc,
                    method=method,
                    periodic=wrapped,
                    nan_policy=nan_policy,
                    **geom_kw,
                )
                second = cartesian.cartesian_partial(
                    first,
                    dim,
                    accuracy=acc,
                    method=method,
                    periodic=wrapped,
                    nan_policy=nan_policy,
                    **geom_kw,
                )
            total = second if total is None else total + second
        assert total is not None  # narrowed by the empty-dims guard above
        return total.rename(None)

    grad = gradient(
        da,
        dims=target_dims,
        geometry=geometry,
        accuracy=per_dim,
        method=method,
        periodic=sorted(wrap),  # type: ignore[arg-type]
        nan_policy=nan_policy,
        **geom_kw,
    )
    component_names = tuple(grad.data_vars)
    return divergence(
        grad,
        components=component_names,  # type: ignore[arg-type]
        dims=target_dims,
        geometry=geometry,
        accuracy=per_dim,
        method=method,
        periodic=sorted(wrap),  # type: ignore[arg-type]
        nan_policy=nan_policy,
        **geom_kw,
    )

grid_metrics_from_coords

grid_metrics_from_coords(ds: Dataset, *, lat: str = 'lat', lon: str = 'lon', depth: str | None = None, sphere: bool = True, radius: float = EARTH_RADIUS) -> tuple[Dataset, Dataset]

Derive (volume_metrics, face_metrics) from a Dataset's coords.

Parameters:

Name	Type	Description	Default
ds	Dataset	Dataset whose coordinates define the grid.	required
lat	str	Name of the latitude coordinate. Spherical mode assumes degrees.	'lat'
lon	str	Name of the longitude coordinate. Spherical mode assumes degrees.	'lon'
depth	str | None	Optional vertical coordinate name. If None, dz is omitted from volume_metrics and cell_volume is populated with the 2-D cell_area (in m², not m³) so downstream code that always reads cell_volume keeps working in 2-D mode.	None
sphere	bool	If True (default), apply spherical metric dx = R cos(φ) Δλ, dy = R Δφ with degrees-to-radians conversion. If False, treat lon/lat as Cartesian with units already in metres.	True
radius	float	Earth radius (m), used in spherical mode.	EARTH_RADIUS

Returns:

Type	Description
Dataset	(volume_metrics, face_metrics) — two Datasets with the
Dataset	keys documented in the module docstring.

Notes

Edge cells extrapolate the nearest interior spacing; this is adequate for budget closure tests and the demo notebook but will bias real-world budgets near the domain boundary.

Example

A 1°×1° grid near the equator has ≈111 km cell edges:

>>> import xarray as xr
>>> ds = xr.Dataset(coords={"lat": [0.0, 1.0, 2.0], "lon": [10.0, 11.0, 12.0]})
>>> vol, face = grid_metrics_from_coords(ds)
>>> sorted(vol.data_vars)
['cell_area', 'cell_volume', 'dx', 'dy']
>>> sorted(face.data_vars)
['area_e', 'area_n', 'dx_e', 'dy_n']
>>> round(float(vol["dy"].isel(lat=0)) / 1000)  # km per degree latitude
111

Source code in packages/xrtoolz-ops/src/xrops/grad/_src/grid_metrics.py

def grid_metrics_from_coords(
    ds: xr.Dataset,
    *,
    lat: str = "lat",
    lon: str = "lon",
    depth: str | None = None,
    sphere: bool = True,
    radius: float = EARTH_RADIUS,
) -> tuple[xr.Dataset, xr.Dataset]:
    """Derive ``(volume_metrics, face_metrics)`` from a Dataset's coords.

    Args:
        ds: Dataset whose coordinates define the grid.
        lat: Name of the latitude coordinate. Spherical mode assumes
            degrees.
        lon: Name of the longitude coordinate. Spherical mode assumes
            degrees.
        depth: Optional vertical coordinate name. If ``None``, ``dz``
            is omitted from ``volume_metrics`` and ``cell_volume`` is
            populated with the 2-D ``cell_area`` (in m², not m³) so
            downstream code that always reads ``cell_volume`` keeps
            working in 2-D mode.
        sphere: If ``True`` (default), apply spherical metric
            ``dx = R cos(φ) Δλ``, ``dy = R Δφ`` with degrees-to-radians
            conversion. If ``False``, treat ``lon``/``lat`` as Cartesian
            with units already in metres.
        radius: Earth radius (m), used in spherical mode.

    Returns:
        ``(volume_metrics, face_metrics)`` — two Datasets with the
        keys documented in the module docstring.

    Notes:
        Edge cells extrapolate the nearest interior spacing; this is
        adequate for budget closure tests and the demo notebook but
        will bias real-world budgets near the domain boundary.

    Example:
        A 1°×1° grid near the equator has ≈111 km cell edges:

        ```pycon
        >>> import xarray as xr
        >>> ds = xr.Dataset(coords={"lat": [0.0, 1.0, 2.0], "lon": [10.0, 11.0, 12.0]})
        >>> vol, face = grid_metrics_from_coords(ds)
        >>> sorted(vol.data_vars)
        ['cell_area', 'cell_volume', 'dx', 'dy']
        >>> sorted(face.data_vars)
        ['area_e', 'area_n', 'dx_e', 'dy_n']
        >>> round(float(vol["dy"].isel(lat=0)) / 1000)  # km per degree latitude
        111

        ```
    """
    if lat not in ds.coords:
        raise ValueError(f"Dataset is missing latitude coord {lat!r}.")
    if lon not in ds.coords:
        raise ValueError(f"Dataset is missing longitude coord {lon!r}.")
    lat_c = ds[lat]
    lon_c = ds[lon]
    if sphere:
        rad = float(np.pi / 180.0)
        dlat_rad = _cell_widths_from_coord(lat_c) * rad
        dlon_rad = _cell_widths_from_coord(lon_c) * rad
        dlat_face_rad = _face_widths_from_centres(lat_c) * rad
        dlon_face_rad = _face_widths_from_centres(lon_c) * rad
        cos_phi = np.cos(lat_c * rad)
        dx = (radius * cos_phi * dlon_rad).rename("dx")
        dy = (radius * dlat_rad).rename("dy")
        dx_e = (radius * cos_phi * dlon_face_rad).rename("dx_e")
        dy_n = (radius * dlat_face_rad).rename("dy_n")
    else:
        dx = _cell_widths_from_coord(lon_c).rename("dx")
        dy = _cell_widths_from_coord(lat_c).rename("dy")
        dx_e = _face_widths_from_centres(lon_c).rename("dx_e")
        dy_n = _face_widths_from_centres(lat_c).rename("dy_n")

    cell_area = (dx * dy).rename("cell_area")
    area_e = (dx_e * dy).rename("area_e")
    area_n = (dx * dy_n).rename("area_n")

    vol_data: dict[str, xr.DataArray] = {
        "dx": dx,
        "dy": dy,
        "cell_area": cell_area,
    }
    face_data: dict[str, xr.DataArray] = {
        "dx_e": dx_e,
        "dy_n": dy_n,
        "area_e": area_e,
        "area_n": area_n,
    }

    if depth is not None:
        if depth not in ds.coords:
            raise ValueError(f"Dataset is missing depth coord {depth!r}.")
        dz = _cell_widths_from_coord(ds[depth]).rename("dz")
        dz_top = _face_widths_from_centres(ds[depth]).rename("dz_top")
        cell_volume = (cell_area * dz).rename("cell_volume")
        area_top = cell_area.rename("area_top")
        vol_data["dz"] = dz
        vol_data["cell_volume"] = cell_volume
        face_data["dz_top"] = dz_top
        face_data["area_top"] = area_top
    else:
        vol_data["cell_volume"] = cell_area.rename("cell_volume")

    return xr.Dataset(vol_data), xr.Dataset(face_data)

Physical constants

EARTH_RADIUS module-attribute

EARTH_RADIUS: float = 6371000.0

Mean Earth radius in metres.

GRAVITY module-attribute

GRAVITY: float = 9.80665

Standard gravitational acceleration in m/s².

OMEGA module-attribute

OMEGA: float = 7.2921159e-05

Earth's rotation rate in rad/s.