NEMODataTree API

nemo_cookbook.NEMODataTree

Bases: DataTree

Hierarchical data structure containing collections of NEMO ocean model outputs.

This class extends xarray.DataTree to provide methods for processing and analysing NEMO output xarray objects defining one or more model domains.

It supports NEMO discrete scalar and vector operators such as computing gradients, divergence, curl, weighted averages, integrals, cumulative integrals, and transforming variables between grids.

Methods:

Name	Description
__getitem__	Access child nodes, variables, or coordinates stored in this NEMODataTree.
__init__	Create a single node of a NEMODataTree.
__setitem__	Set a child node or variable in this NEMODataTree.
add_geoindex	Add geographical index variables to a given NEMO model grid.
binned_statistic	Calculate binned statistic of a variable defined on a NEMO model grid.
cell_area	Calculate grid cell areas orthogonal to a given dimension of a NEMO model grid.
cell_volume	Calculate grid cell volumes for a given NEMO model grid.
clip_domain	Clip a NEMO model domain to specified longitude and latitude range.
clip_grid	Clip a NEMO model grid to specified longitude and latitude range.
curl	Calculate the vertical (k) curl component of a vector field on a NEMO model grid.
depth_integral	Integrate a variable in depth coordinates between two limits.
divergence	Calculate the horizontal divergence of a vector field defined on a NEMO model grid.
extract_mask_boundary	Extract the boundary of a masked region defined on a NEMO model grid.
extract_section	Extract hydrographic section from a NEMO model domain.
from_datasets	Create a NEMODataTree from a dictionary of xarray.Dataset objects created from NEMO model output files,
from_icechunk	Create a NEMODataTree from an Icechunk repository storing NEMO model outputs
from_paths	Create a NEMODataTree from a dictionary of paths to NEMO model output files,
from_zarr	Create a NEMODataTree from Zarr store groups storing NEMO model outputs
gradient	Calculate the gradient of a scalar variable along one dimension (e.g., 'i', 'j', 'k') of a NEMO model grid.
integral	Integrate a variable along one or more dimensions of a NEMO model grid.
mask_with_polygon	Create mask of NEMO model grid points contained within a polygon.
masked_statistic	Masked statistic of a variable defined on a NEMO model grid.
transform_to	Transform variable defined on a NEMO model grid to a neighbouring
transform_vertical_grid	Transform variable defined on a NEMO model grid to a new vertical grid using conservative interpolation.

Source code in nemo_cookbook/nemodatatree.py

class NEMODataTree(xr.DataTree):
    """
    Hierarchical data structure containing collections of NEMO ocean model outputs.

    This class extends `xarray.DataTree` to provide methods for processing
    and analysing NEMO output xarray objects defining one or more model domains.

    It supports NEMO discrete scalar and vector operators such as computing gradients,
    divergence, curl, weighted averages, integrals, cumulative integrals, and
    transforming variables between grids.
    """

    def __init__(self, *args, **kwargs):
        """
        Create a single node of a NEMODataTree.

        The node may optionally contain data in the form of data
        and coordinate variables, stored in the same way as data
        is stored in an `xarray.Dataset`.

        Parameters
        ----------
        *args : tuple
            Positional arguments to pass to the parent class.
        **kwargs : dict
            Keyword arguments to pass to the parent class.
        """
        super().__init__(*args, **kwargs)

    @classmethod
    def from_paths(
        cls,
        paths: dict[str, str],
        nests: dict[str, str] | None = None,
        name : str = "NEMO model",
        iperio: bool = False,
        nftype: str | None = None,
        read_mask: bool = False,
        key_linssh: bool = False,
        nbghost_child: int = 4,
        **open_kwargs: dict[str, any],
    ) -> Self:
        """
        Create a NEMODataTree from a dictionary of paths to NEMO model output files,
        organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

        Parameters
        ----------
        paths : dict[str, str]
            Dictionary containing paths to NEMO grid files, structured as:
            {
                'parent': {'domain': 'path/to/domain.nc',
                           'gridT': 'path/to/gridT.nc',
                            , ... ,
                            'icemod': 'path/to/icemod.nc',
                            },
                'child': {'1': {'domain': 'path/to/child_domain.nc',
                                'gridT': 'path/to/child_gridT.nc',
                                , ... ,
                                'icemod': 'path/to/child_icemod.nc',
                                },
                          },
                'grandchild': {'2': {'domain': 'path/to/grandchild_domain.nc',
                                     'gridT': 'path/to/grandchild_gridT.nc',
                                     , ...,
                                     'icemod': 'path/to/grandchild_icemod.nc',
                                     },
                               }
            }

        nests : dict[str, str], optional
            Dictionary describing the properties of nested domains, structured as:
            {
                "1": {
                    "parent": "/",
                    "rx": rx,
                    "ry": ry,
                    "imin": imin,
                    "imax": imax,
                    "jmin": jmin,
                    "jmax": jmax,
                    "iperio": iperio,
                    },
            }
            where `rx` and `ry` are the horizontal refinement factors, and `imin`, `imax`, `jmin`, `jmax`
            define the indices of the child (grandchild) domain within the parent (child) domain. Zonally
            periodic nested domains should be specified with `iperio=True`.

        name: str, optional
            Name of the NEMODataTree. Default is "NEMO model".

        iperio: bool = False
            Zonal periodicity of the parent domain. Default is False.

        nftype: str, optional
            Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point
            pivot. By default, no north fold lateral boundary condition is applied (None).

        read_mask: bool = False
            If True, read NEMO model land/sea mask from domain files. Default is False, meaning masks are computed from top_level and bottom_level domain variables.

        key_linssh: bool = False
            Linear free-surface approximation. If True, vertical coordinates are time-independent and given by (e3t_0, e3u_0, e3v_0, e3w_0) in domain_cfg.
            If False, vertical coordinates are time-dependent and must be specified in NEMO model grid datasets. Default is False.

        nbghost_child : int = 4
            Number of ghost cells to remove from the western/southern boundaries of the (grand)child domains. Default is 4.

        **open_kwargs : dict, optional
            Additional keyword arguments to pass to `xarray.open_dataset` or `xr.open_mfdataset` when opening NEMO model output files.
        Returns
        -------
        NEMODataTree
            A hierarchical DataTree storing NEMO model outputs.

        Examples
        --------
        Create a zonally periodic `NEMODataTree` with north folding on T-points from a dictionary of paths to local netCDF files:

        >>> from nemo_cookbook import NEMODataTree
        >>> paths = {"parent": {
        ...          "domain": "/path/to/domain_cfg.nc",
        ...          "gridT": "path/to/*_gridT.nc",
        ...          "gridU": "path/to/*_gridV.nc",
        ...          "gridV": "path/to/*_gridV.nc",
        ...          "gridW": "path/to/*_gridW.nc",
        ...          "icemod": "path/to/*_icemod.nc",
        ...          }}
        >>> nemo = NEMODataTree.from_paths(paths, name="My NEMO model", iperio=True, nftype="T")

        Create a regional `NEMODataTree` using a linear free-surface approximation from a dictionary of paths to remote netCDF files:

        >>> nemo = NEMODataTree.from_paths(paths, name="My NEMO model", iperio=False, nftype=None, key_linssh=True)

        See Also
        --------
        from_datasets
        """
        # -- Validate Inputs -- #
        if not isinstance(paths, dict):
            raise TypeError("`paths` must be a dictionary or nested dictionary.")
        if not isinstance(nests, (dict, type(None))):
            raise TypeError("`nests` must be a dictionary or None.")
        if not isinstance(name, str):
            raise TypeError("`name` must be a string.")
        if not isinstance(iperio, bool):
            raise TypeError("zonal periodicity (`iperio`) of parent domain must be a boolean.")
        if nftype is not None and nftype not in ("T", "F"):
            raise ValueError(
                "north fold type (`nftype`) of parent domain must be 'T' (T-pivot fold), 'F' (F-pivot fold), or None."
            )
        if not isinstance(read_mask, bool):
            raise TypeError("`read_mask` must be a boolean.")
        if not isinstance(key_linssh, bool):
            raise TypeError("linear free-surface approximation (`key_linssh`) must be a boolean.")
        if not isinstance(nbghost_child, int):
            raise TypeError(
                "number of ghost cells along the western/southern boundaries (`nbghost_child`) must be an integer."
            )
        if not isinstance(open_kwargs, dict):
            raise TypeError("`open_kwargs` must be a dictionary.")

        # -- Define parent, child, grandchild filepath collections -- #
        d_child, d_grandchild = None, None
        if "parent" in paths.keys() and isinstance(paths["parent"], dict):
            for key in paths.keys():
                if key not in ("parent", "child", "grandchild"):
                    raise KeyError(f"Unexpected key '{key}' in `paths` dictionary.")
                if key == "parent":
                    d_parent = paths["parent"]
                elif key == "child":
                    if nests is None:
                        raise ValueError(
                            "`nests` dictionary must be provided when defining NEMO child domains."
                        )
                    else:
                        d_child = paths["child"]
                elif key == "grandchild":
                    d_grandchild = paths["grandchild"]
        else:
            raise ValueError(
                "Invalid `paths` structure. Expected a nested dictionary defining NEMO 'parent', 'child' and 'grandchild' domains."
            )

        # -- Construct DataTree from parent / child / grandchild domains -- #
        d_tree = create_datatree_dict(
            d_parent=d_parent,
            d_child=d_child,
            d_grandchild=d_grandchild,
            nests=nests,
            iperio=iperio,
            nftype=nftype,
            read_mask=read_mask,
            nbghost_child=nbghost_child,
            key_linssh=key_linssh,
            open_kwargs=dict(**open_kwargs),
        )

        nemo = super().from_dict(d_tree)
        nemo.name = name

        return nemo

    @classmethod
    def from_datasets(
        cls,
        datasets: dict[str, dict[str, xr.Dataset]],
        nests: dict[str, dict[str, str]] | None = None,
        name : str = "NEMO model",
        iperio: bool = False,
        nftype: str | None = None,
        read_mask: bool = False,
        key_linssh: bool = False,
        nbghost_child: int = 4,
    ) -> Self:
        """
        Create a NEMODataTree from a dictionary of `xarray.Dataset` objects created from NEMO model output files,
        organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

        Parameters
        ----------
        datasets : dict[str, dict[str, xr.Dataset]]
            Dictionary containing `xarray.Datasets` created from NEMO grid files, structured as:
            {
                'parent': {'domain': ds_domain, 'gridT': ds_gridT, ... , 'icemod': ds_icemod.nc},
                'child': {'1': {'domain': ds_domain_1, 'gridT': d_gridT_1, ...}},
                'grandchild': {'2': {'domain': ds_domain_2, 'gridT': ds_gridT_2, ...}}
            }

        nests : dict[str, dict[str, str]], optional
            Dictionary describing the properties of nested domains, structured as:
            {
                "1": {
                    "parent": "/",
                    "rx": rx,
                    "ry": ry,
                    "imin": imin,
                    "imax": imax,
                    "jmin": jmin,
                    "jmax": jmax,
                    },
            }
            where `rx` and `ry` are the horizontal refinement factors, and `imin`, `imax`, `jmin`, `jmax`
            define the indices of the child (grandchild) domain within the parent (child) domain.

        name: str, optional
            Name of the NEMODataTree. Default is "NEMO model".

        iperio: bool = False
            Zonal periodicity of the parent domain.

        nftype: str, optional
            Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point
            pivot. By default, no north fold lateral boundary condition is applied (None).

        read_mask: bool = False
            If True, read NEMO model land/sea mask from domain files. Default is False, meaning masks are computed from top_level and bottom_level domain variables.

        key_linssh: bool = False
            Linear free-surface approximation. If True, vertical coordinates are time-independent and given by (e3t_0, e3u_0, e3v_0, e3w_0) in domain_cfg.
            If False, vertical coordinates are time-dependent and must be specified in NEMO model grid datasets. Default is False.

        nbghost_child : int = 4
            Number of ghost cells to remove from the western/southern boundaries of the (grand)child domains. Default is 4.

        Returns
        -------
        NEMODataTree
            Hierarchical DataTree of NEMO model outputs.

        Examples
        --------
        Create a zonally periodic `NEMODataTree` with north folding on T-points from a dictionary of xarray.Dataset objects:

        >>> import xarray as xr
        >>> from nemo_cookbook import NEMODataTree
        >>> ds_domain = xr.open_zarr("https://some_remote_data/domain_cfg.zarr")
        >>> ds_gridT = xr.open_zarr("https://some_remote_data/my_model_gridT.zarr")
        >>> datasets = {"parent": {"domain": ds_domain, "gridT": ds_gridT}}
        >>> nemo = NEMODataTree.from_datasets(datasets=datasets, name="My NEMO Model", iperio=True, nftype="T")

        Create a regional `NEMODataTree` using a linear free-surface approximation from a dictionary of xarray.Dataset objects:

        >>> nemo = NEMODataTree.from_datasets(datasets=datasets, name="My NEMO Model", iperio=False, nftype=None, key_linssh=True)

        See Also
        --------
        from_paths
        """
        # -- Validate Inputs -- #
        if not isinstance(datasets, dict):
            raise TypeError("`datasets` must be a dictionary or nested dictionary.")
        if not isinstance(nests, (dict, type(None))):
            raise TypeError("`nests` must be a dictionary or None.")
        if not isinstance(name, str):
            raise TypeError("`name` must be a string.")
        if not isinstance(iperio, bool):
            raise TypeError("zonal periodicity (`iperio`) of parent domain must be a boolean.")
        if nftype is not None and nftype not in ("T", "F"):
            raise ValueError(
                "north fold type (`nftype`) of parent domain must be 'T' (T-pivot fold), 'F' (F-pivot fold), or None."
            )
        if not isinstance(read_mask, bool):
            raise TypeError("`read_mask` must be a boolean.")
        if not isinstance(key_linssh, bool):
            raise TypeError("linear free-surface approximation (`key_linssh`) must be a boolean.")
        if not isinstance(nbghost_child, int):
            raise TypeError(
                "number of ghost cells along the western/southern boundaries (`nbghost_child`) must be an integer."
            )

        # -- Define parent, child, grandchild dataset collections -- #
        d_child, d_grandchild = None, None
        if "parent" in datasets.keys() and isinstance(datasets["parent"], dict):
            for key in datasets.keys():
                if key not in ("parent", "child", "grandchild"):
                    raise KeyError(f"Unexpected key '{key}' in `datasets` dictionary.")
                if key == "parent":
                    d_parent = datasets["parent"]
                elif key == "child":
                    if nests is None:
                        raise ValueError(
                            "`nests` dictionary must be provided when defining NEMO child domains."
                        )
                    else:
                        d_child = datasets["child"]
                elif key == "grandchild":
                    d_grandchild = datasets["grandchild"]
        else:
            raise ValueError(
                "Invalid `datasets` structure. Expected a nested dictionary defining NEMO 'parent', 'child' and 'grandchild' domains."
            )

        # -- Construct DataTree from parent / child / grandchild domains -- #
        d_tree = create_datatree_dict(
            d_parent=d_parent,
            d_child=d_child,
            d_grandchild=d_grandchild,
            nests=nests,
            iperio=iperio,
            nftype=nftype,
            read_mask=read_mask,
            key_linssh=key_linssh,
            nbghost_child=nbghost_child,
        )

        nemo = super().from_dict(d_tree)
        nemo.name = name

        return nemo

    @classmethod
    def from_icechunk(
        cls,
        repo: icechunk.repository.Repository,
        name: str = "NEMO model",
        iperio: bool = False,
        nftype: str | None = None,
        open_kwargs: dict[str, any] | None = None,
        **session_kwargs: dict[str, any],
    ) -> Self:
        """
        Create a NEMODataTree from an Icechunk repository storing NEMO model outputs
        organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

        Parameters
        ----------
        repo : icechunk.repository.Repository
            Icechunk repository containing NEMO model outputs organised into a
            hierarchy of domains.

        name : str, optional
            Name of the NEMODataTree. Default is "NEMO model".

        iperio: bool = False
            Zonal periodicity of the parent domain. Default is False.

        nftype: str, optional
            Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point

        open_kwargs : dict[str, Any], optional
            Additional keyword arguments to pass to `xarray.open_datatree`. Default is None.

        **session_kwargs : dict[str, Any], optional
            Additional keyword arguments to pass to Icechunk `repo.readonly_session`.

        Returns
        -------
        NEMODataTree
            Hierarchical DataTree of NEMO model outputs.

        Examples
        --------
        Create a zonally periodic `NEMODataTree` with north folding on F-points from the main branch of an Icechunk repository:

        >>> from nemo_cookbook import NEMODataTree
        >>> nemo = NEMODataTree.from_icechunk(repo=repo, branch="main", iperio=True, nftype="F")

        See Also
        --------
        from_zarr
        """
        # -- Validate Inputs -- #
        if not hasattr(repo, "readonly_session"):
            raise TypeError("`repo` must implement readonly_session().")
        if not isinstance(name, str):
            raise TypeError("`name` must be a string.")
        if not isinstance(iperio, bool):
            raise TypeError("zonal periodicity (`iperio`) of parent domain must be a boolean.")
        if nftype is not None and nftype not in ("T", "F"):
            raise ValueError(
                "north fold type (`nftype`) of parent domain must be 'T' (T-pivot fold), 'F' (F-pivot fold), or None."
            )

        # -- Create NEMODataTree from Icechunk repository -- #:
        session = repo.readonly_session(**session_kwargs)
        datatree = xr.open_datatree(session.store, engine="zarr", **(open_kwargs or {}))
        nemo = super().from_dict(datatree.to_dict())

        # -- Update NEMODataTree properties -- #
        nemo["/"].attrs.update({"nftype": nftype, "iperio": iperio})
        nemo.name = name

        # -- Validate NEMO grid node Datasets -- #
        for key in [grid for grid in nemo.groups if grid.startswith("grid")]:
            validate_nemo_grid_node(key=key, value=nemo[key])

        return nemo

    @classmethod
    def from_zarr(
        cls,
        store: str,
        name: str = "NEMO model",
        iperio: bool = False,
        nftype: str | None = None,
        **open_kwargs: dict[str, any],
    ) -> Self:
        """
        Create a NEMODataTree from Zarr store groups storing NEMO model outputs
        organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

        Parameters
        ----------
        store : str
            Path to the Zarr store containing NEMO model outputs in hierarchical groups.

        name : str, optional
            Name of the NEMODataTree. Default is "NEMO model".

        iperio: bool = False
            Zonal periodicity of the parent domain. Default is False.

        nftype: str, optional
            Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point

        **open_kwargs : dict[str, Any], optional
            Additional keyword arguments to pass to `xarray.open_datatree`.

        Returns
        -------
        NEMODataTree
            Hierarchical DataTree of NEMO model outputs.

        Examples
        --------
        Create a zonally periodic `NEMODataTree` with north folding on T-points from a hierarchical Zarr store:

        >>> from nemo_cookbook import NEMODataTree
        >>> nemo = NEMODataTree.from_zarr(store="path/to/zarr/store", iperio=True, nftype="T")

        See Also
        --------
        from_icechunk
        """
        # -- Validate Inputs -- #
        if not isinstance(store, str):
            raise TypeError("`store` must be a string.")
        if not isinstance(name, str):
            raise TypeError("`name` must be a string.")
        if not isinstance(iperio, bool):
            raise TypeError("zonal periodicity (`iperio`) of parent domain must be a boolean.")
        if nftype is not None and nftype not in ("T", "F"):
            raise ValueError(
                "north fold type (`nftype`) of parent domain must be 'T' (T-pivot fold), 'F' (F-pivot fold), or None."
            )

        # -- Create NEMODataTree from Zarr store -- #:
        datatree = xr.open_datatree(store, engine="zarr", **(open_kwargs or {}))
        nemo = super().from_dict(datatree.to_dict())

        # -- Update NEMODataTree properties -- #
        nemo["/"].attrs.update({"nftype": nftype, "iperio": iperio})
        nemo.name = name

        # -- Validate NEMO grid node Datasets -- #
        for key in [grid for grid in nemo.groups if grid.startswith("grid")]:
            validate_nemo_grid_node(key=key, value=nemo[key])

        return nemo

    def __setitem__(
            self,
            key: str,
            value: NEMODataArray | xr.DataArray | xr.Dataset,
            strict: bool = True
        ) -> None:
        """
        Set a child node or variable in this NEMODataTree.

        Overloads the __setitem__() method of xarray.DataTree to allow
        setting NEMODataArrays via variable paths (i.e, /grid/var).

        Optionally set strict=False to bypass validation of child grid nodes.

        Parameters
        ----------
        key : str
            Name of variable or child node, or unix-like path to variable
            within a child node.

        value : NEMODataArray | xarray.DataArray | xarray.Dataset
            Object to set at the specified key. If a NEMODataArray is provided,
            the underlying xarray.DataArray will be set at the specified key.

        strict : bool, optional
            Validate Datasets assigned to NEMO grid nodes to ensure they contain
            the required dimensions and coordinates. Default is True.

        Returns
        -------
        None
        """
        # -- Access DataArray from NEMODataArray -- #
        if isinstance(value, NEMODataArray):
            value = value.data

        if strict and isinstance(value, xr.Dataset):
            # -- Validate NEMO grid node Dataset -- #
            validate_nemo_grid_node(key=key, value=value)

        return super().__setitem__(key, value)

    def __getitem__(self, key: str) -> Self | NEMODataArray:
        """
        Access child nodes, variables, or coordinates stored in this NEMODataTree.

        Returned object will be either a NEMODataTree or NEMODataArray object depending
        on whether the key points to a child node or variable.

        Overloads the __getitem__() method of xarray.DataTree to return NEMODataArrays
        accessed via variable paths (i.e, /grid/var).

        Parameters
        ----------
        key : str
            Name of variable / child within this node, or unix-like path to variable
            / child within another node.

        Returns
        -------
        NEMODataTree | NEMODataArray
        """
        # -- Access child node or variable -- #
        item = super().__getitem__(key=key)
        is_gridpath = key.startswith("/grid") or key.startswith("grid")

        # -- Return NEMODataArray -- #
        if isinstance(item, xr.DataArray) and is_gridpath:
            var_name = key.split("/")[-1]
            grid = key.replace(f"/{var_name}", "")
            item = NEMODataArray(da=item,
                                 tree=self,
                                 grid=grid
                                 )

        return item

    def _get_properties(
        self, dom: str | None = None, grid: str | None = None, infer_dom: bool = False
    ) -> str | tuple[str, ...]:
        """
        Get NEMO model domain and grid properties.

        The domain prefix & suffix (e.g., '1_', '1') are returned
        if only the NEMO model domain (`dom`) is specified.

        The grid suffix (e.g., 't', 'u', 'v', 'w') is returned if
        only the NEMO model grid (`grid`) is specified.

        The domain number, domain prefix & suffix, and grid suffix
        are returned if both the NEMO model grid (`grid`) and
        `infer_dom = True` are specified.

        Parameters
        ----------
        dom : str, optional
            Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
        grid : str, optional
            Path to NEMO model grid (e.g., 'gridT').
        infer_dom : bool, optional
            Whether to infer the domain number & domain name from only the
            grid path. Default is False.

        Returns
        -------
        str | tuple[str, ...]
            NEMO model domain and grid properties.
        """
        if (grid is None) and (dom is not None):
            dom_prefix = "" if dom == "." else f"{dom}_"
            dom_suffix = "" if dom == "." else f"{dom}"
            return dom_prefix, dom_suffix
        else:
            grid_keys = list(dict(self.subtree_with_keys).keys())
            if grid not in grid_keys:
                raise KeyError(
                    f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
                )
            grid_suffix = f"{grid.lower()[-1]}"

            if infer_dom:
                dom_inds = [char for char in grid if char.isdigit()]
                dom_prefix = f"{dom_inds[-1]}_" if len(dom_inds) != 0 else ""
                dom = dom_prefix[:-1] if dom_prefix != "" else "."
                dom_suffix = dom if dom != "." else ""
                return dom, dom_prefix, dom_suffix, grid_suffix
            else:
                return grid_suffix

    def _get_grid_paths(self, dom: str) -> dict[str, str]:
        """
        Get paths to NEMO model grids in a given domain.

        Parameters
        ----------
        dom : str, optional
            Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).

        Returns
        -------
        dict[str, str]
            Dictionary of NEMO model grid paths.
        """
        # Collect paths to all NEMO model grids:
        grid_paths = list(dict(self.subtree_with_keys).keys())

        if dom == ".":
            grid_paths = [
                path for path in grid_paths if ("_" not in path) and ("grid" in path)
            ]
        else:
            grid_paths = [path for path in grid_paths if dom in path]

        d_paths = {path.split("/")[0]: path for path in grid_paths}

        return d_paths

    def _get_ijk_names(self, dom: str | None = None, grid: str | None = None) -> dict[str, str]:
        """
        Get (i, j, k) grid index names for given NEMO model domain.

        If path to NEMO model grid is provided, domain is inferred.

        Parameters
        ----------
        dom : str, optional
            Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
        grid : str, optional
            Path to NEMO model grid (e.g., 'gridT').

        Returns
        -------
        dict[str, str]
            NEMO model grid index names.
        """
        if grid is not None:
            dom, _, dom_suffix, _ = self._get_properties(grid=grid, infer_dom=True)
        else:
            _, dom_suffix = self._get_properties(dom=dom)

        indexes = ["i", "j", "k"]
        if dom == ".":
            d_ijk = {index: index for index in indexes}
        else:
            d_ijk = {index: f"{index}{dom_suffix}" for index in indexes}

        return d_ijk

    def _get_weights(self, grid: str, dims: list, fillna: bool = True) -> xr.DataArray:
        """
        Get the weights (scale factors) for specified dimensions of a NEMO model grid.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid where weights are stored (e.g., 'gridT').
        dims : list
            Dimensions to collect weights for.
        fillna : bool, optional
            Fill NaN values in weights with zeros. Default is True.

        Returns
        -------
        xr.DataArray
            Weights (scale factors) for the specified dimensions of the NEMO model grid.
        """
        if any(dim not in ["i", "j", "k"] for dim in dims):
            raise ValueError(
                "dims must be a list containing one or more of the following dimensions: ['i', 'j', 'k']."
            )

        grid_suffix = self._get_properties(grid=grid)

        weights_dict = {
            "i": f"e1{grid_suffix}",
            "j": f"e2{grid_suffix}",
            "k": f"e3{grid_suffix}",
        }
        try:
            weights_list = [self[f"{grid}"][weights_dict[dim]] for dim in dims]
        except KeyError as e:
            raise KeyError(
                f"weights missing for dimensions {dims} of NEMO model grid {grid}"
            ) from e

        if len(weights_list) == 1:
            weights = weights_list[0]
        elif len(weights_list) == 2:
            weights = weights_list[0] * weights_list[1]
        elif len(weights_list) == 3:
            weights = weights_list[0] * weights_list[1] * weights_list[2]

        if fillna:
            weights = weights.fillna(value=0)

        return weights

    def add_geoindex(
        self,
        grid: str,
    ) -> Self:
        """
        Add geographical index variables to a given NEMO model grid.

        This enables users to index grid variables using geographical
        coordinates (e.g., glamt, gphit) in addition to their (i, j, k)
        dimensions.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid to add geographical indexes (e.g., 'gridT').

        Returns
        -------
        NEMODataTree
            NEMO DataTree with geographical indexes added to specified model grid.

        Examples
        --------
        Add glamt, gphit as geographical indexes to the T-grid of the NEMO parent domain:

        >>> nemo.add_geoindex(grid="gridT")

        """
        # -- Set geographical indexes -- #
        _, dom_prefix, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
        lon_name = f"{dom_prefix}glam{grid_suffix}"
        lat_name = f"{dom_prefix}gphi{grid_suffix}"
        self_copy = self.copy()
        self_copy[grid] = (
            self_copy[grid]
            .dataset.assign_coords(
                {lat_name: self_copy[grid][lat_name], lon_name: self_copy[grid][lon_name]}
            )
            .set_xindex(
                (lat_name, lon_name),
                NDPointIndex,
                tree_adapter_cls=SklearnGeoBallTreeAdapter,
            )
        )

        return self_copy

    def cell_area(
        self,
        grid: str,
        dim: str,
    ) -> xr.DataArray:
        """
        Calculate grid cell areas orthogonal to a given dimension of a NEMO model grid.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid from which to calculate grid cell areas
            (e.g., 'gridT').
        dim : str
            Dimension orthogonal to grid cell area to
            calculate (e.g., 'k' returns e1 * e2).

        Returns
        -------
        xr.DataArray
            Grid cell areas (m^2) for the specified NEMO model grid.

        Examples
        --------
        Compute the horizontal area of each grid cell centered on a V-grid point
        in the NEMO parent domain:

        >>> nemo.cell_area(grid="gridT", dim="k")

        Note, `dim` represents the dimension orthogonal to the grid cell
        area to be computed.

        See Also
        --------
        cell_volume
        """
        grid_suffix = self._get_properties(grid=grid)

        if dim not in ["i", "j", "k"]:
            raise ValueError(f"dim {dim} must be one of ['i', 'j', 'k'].")

        match dim:
            case "i":
                cell_area = (
                    self[f"{grid}/e3{grid_suffix}"].masked.data * self[f"{grid}/e2{grid_suffix}"].masked.data
                )
            case "j":
                cell_area = (
                    self[f"{grid}/e3{grid_suffix}"].masked.data * self[f"{grid}/e1{grid_suffix}"].masked.data
                )
            case "k":
                cell_area = (
                    self[f"{grid}/e1{grid_suffix}"].masked.data * self[f"{grid}/e2{grid_suffix}"].masked.data
                )
        cell_area.name = "areacello"

        return cell_area

    def cell_volume(self, grid: str) -> xr.DataArray:
        """
        Calculate grid cell volumes for a given NEMO model grid.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid to calculate grid cell volumes
            (e.g., 'gridT').

        Returns
        -------
        xr.DataArray
            Grid cell volumes for the specified NEMO model grid.

        Examples
        --------
        Compute the volume of each grid cell centered on a V-grid point
        in the NEMO parent domain:

        >>> nemo.cell_volumes(grid="gridV")

        See Also
        --------
        cell_area
        """
        grid_suffix = self._get_properties(grid=grid)

        cell_volume = (
            self[f"{grid}/e3{grid_suffix}"].masked.data
            * self[f"{grid}/e1{grid_suffix}"].masked.data
            * self[f"{grid}/e2{grid_suffix}"].masked.data
        )
        cell_volume.name = "volcello"

        return cell_volume

    @deprecated(version_since="2026.03.b1",
                version_removed="2026.05",
                alternative="NEMODataArray.derivative or NEMOVectorField.gradient from v2026.04 onwards"
                )
    def gradient(
        self,
        var: str,
        dim: str,
        dom: str = ".",
    ) -> xr.DataArray:
        """
        Calculate the gradient of a scalar variable along one dimension (e.g., 'i', 'j', 'k') of a NEMO model grid.

        Parameters
        ----------
        var : str
            Name of the scalar variable.
        dim : str
            Dimension along which to calculate gradient (e.g., 'i', 'j', 'k').
        dom : str, optional
            Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
            Default is '.' for the parent domain.

        Returns
        -------
        xr.DataArray
            Gradient of scalar variable defined on a NEMO model grid.

        Examples
        --------
        Compute the 'meridional' gradient of sea surface temperature `tos_con`
        along the NEMO parent domain `j` dimension:

        >>> nemo.gradient(dom='.', var="tos_con", dim="j")

        Compute the vertical gradient of absolute salinity in the first NEMO
        nested child domain:

        >>> nemo.gradient(dom="1", var="so_abs", dim="k")

        See Also
        --------
        integral
        """
        # -- Validate input -- #
        if not isinstance(var, str):
            raise ValueError(
                "var must be a string specifying name of the scalar variable."
            )
        if not isinstance(dim, str):
            raise ValueError(
                "dim must be a string specifying dimension along which to calculate the gradient (e.g., 'i', 'j', 'k')."
            )
        if not isinstance(dom, str):
            raise ValueError(
                "dom must be a string specifying prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
            )

        # -- Get NEMO model grid properties -- #
        dom_prefix, dom_suffix = self._get_properties(dom=dom)
        grid_paths = self._get_grid_paths(dom=dom)
        gridT, gridU, gridV, gridW = (
            grid_paths["gridT"],
            grid_paths["gridU"],
            grid_paths["gridV"],
            grid_paths["gridW"],
        )

        if var not in self[gridT].data_vars:
            raise KeyError(f"variable '{var}' not found in grid '{gridT}'.")

        da = self[f"{gridT}/{var}"].masked.data
        dim_name = f"{dim}{dom_suffix}"
        if dim_name not in da.dims:
            raise KeyError(
                f"dimension '{dim_name}' not found in variable '{var}'. Dimensions available: {da.dims}."
            )

        match dim:
            case "i":
                if f"{dom_prefix}deptht" in da.coords:
                    # 3-dimensional umask:
                    umask = self[gridU]["umask"]
                else:
                    # 2-dimensional umask:
                    umask = self[gridU]["umaskutil"]

                # Zonally Periodic Domain:
                if self[gridT].attrs.get("iperio", False):
                    da_end = da.isel({dim_name: 0})
                    da_end[dim_name] = da[dim_name].max() + 1
                    da = xr.concat([da, da_end], dim=dim_name)
                    dvar = da.diff(dim=dim_name, label="lower")
                else:
                    # Non-Periodic: pad with NaN values after differencing:
                    dvar = da.diff(dim=dim_name, label="lower").pad({dim_name: (0, 1)})
                # Apply u-mask & transform coords -> calculate gradient:
                dvar.coords[dim_name] = dvar.coords[dim_name] + 0.5
                gradient = dvar.where(umask) / self[f"{gridU}/e1u"].masked.data

                # Remove redundant depth coordinates:
                if f"{dom_prefix}deptht" in gradient.coords:
                    gradient = gradient.drop_vars(
                        [f"{dom_prefix}deptht"]
                    ).assign_coords(
                        {f"{dom_prefix}depthu": self[gridU][f"{dom_prefix}depthu"]}
                    )
            case "j":
                # 3-dimensional vmask:
                if f"{dom_prefix}deptht" in da.coords:
                    vmask = self[gridV]["vmask"]
                else:
                    # 2-dimensional vmask (unique points):
                    vmask = self[gridV]["vmaskutil"]

                # Pad with zeros after differencing (zero gradient at jmaxdom):
                dvar = da.diff(dim=dim_name, label="lower").pad(
                    {dim_name: (0, 1)}, constant_values=0
                )
                # Apply vmask & transform coords -> calculate gradient:
                dvar.coords[dim_name] = dvar.coords[dim_name] + 0.5
                gradient = dvar.where(vmask) / self[f"{gridV}/e2v"].masked.data

                if f"{dom_prefix}deptht" in gradient.coords:
                    gradient = gradient.drop_vars(
                        [f"{dom_prefix}deptht"]
                    ).assign_coords(
                        {f"{dom_prefix}depthv": self[gridV][f"{dom_prefix}depthv"]}
                    )

            case "k":
                dvar = da.diff(dim=dim_name, label="lower")
                # Transform coords & apply w-mask -> calculate gradient:
                dvar.coords[dim_name] = dvar.coords[dim_name] + 0.5
                dvar = dvar.where(self[gridW]["wmask"].isel({dim_name: slice(1, None)}))
                try:
                    gradient = -dvar / self[f"{gridW}/e3w"].masked.data.isel(
                        {dim_name: slice(1, None)}
                    )
                    gradient = gradient.drop_vars([f"{dom_prefix}deptht"])
                except KeyError as e:
                    raise KeyError(
                        f"NEMO model grid: '{gridW}' does not contain vertical scale factor 'e3w' required to calculate gradients along the k-dimension."
                    ) from e

        # Update DataArray properties:
        gradient.name = f"grad_{dim_name}({var})"
        gradient = gradient.drop_vars([f"{dom_prefix}glamt", f"{dom_prefix}gphit"])

        return gradient

    @deprecated(version_since="2026.03.b1",
                version_removed="2026.05",
                alternative="NEMOVectorField.divergence from v2026.04 onwards"
                )
    def divergence(
        self,
        vars: list[str],
        dom: str = ".",
    ) -> xr.DataArray:
        """
        Calculate the horizontal divergence of a vector field defined on a NEMO model grid.

        Parameters
        ----------
        vars : list[str]
            Name of vector variables, structured as: ['u', 'v'], where
            'u' and 'v' are the i and j components of the vector field,
            respectively.
        dom : str, optional
            Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
            Default is '.' for the parent domain.

        Returns
        -------
        xr.DataArray
            Horizontal divergence of vector field defined on a NEMO model grid.

        Examples
        --------
        Compute the horizontal divergence of the seawater velocity field in the
        NEMO parent domain:

        >>> nemo.divergence(dom=".", vars=["uo", "vo"])

        Note, `vars` expects a list of the `i` and `j` components of the vector
        field, respectively.

        See Also
        --------
        divergence
        """
        # -- Validate input -- #
        if not isinstance(vars, list) or len(vars) != 2:
            raise ValueError(
                "vars must be a list of two elements structured as ['u', 'v']."
            )
        if not isinstance(dom, str):
            raise ValueError(
                "dom must be a string specifying the prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
            )

        # -- Get NEMO model grid properties -- #
        dom_prefix, _ = self._get_properties(dom=dom)
        grid_paths = self._get_grid_paths(dom=dom)
        gridT, gridU, gridV = (
            grid_paths["gridT"],
            grid_paths["gridU"],
            grid_paths["gridV"],
        )
        ijk_names = self._get_ijk_names(dom=dom)
        i_name, j_name = ijk_names["i"], ijk_names["j"]

        # -- Define i,j vector components -- #
        var_i, var_j = vars[0], vars[1]
        if var_i not in self[gridU].data_vars:
            raise KeyError(f"variable '{var_i}' not found in grid '{gridU}'.")
        if var_j not in self[gridV].data_vars:
            raise KeyError(f"variable '{var_j}' not found in grid '{gridV}'.")

        da_i = self[f"{gridU}/{var_i}"].masked.data
        da_j = self[f"{gridV}/{var_j}"].masked.data

        # -- Collect mask -- #
        if (f"{dom_prefix}depthu" in da_i.coords) and (
            f"{dom_prefix}depthv" in da_j.coords
        ):
            # 3-dimensional tmask:
            tmask = self[gridT]["tmask"]
        else:
            # 2-dimensional tmask (unique points):
            tmask = self[gridT]["tmaskutil"]

        # -- Neglecting the first T-grid points along i, j dimensions -- #
        e1t = self[f"{gridT}/e1t"].masked.data.isel({i_name: slice(1, None), j_name: slice(1, None)})
        e2t = self[f"{gridT}/e2t"].masked.data.isel({i_name: slice(1, None), j_name: slice(1, None)})
        e3t = self[f"{gridT}/e3t"].masked.data.isel({i_name: slice(1, None), j_name: slice(1, None)})

        e2u, e3u = self[f"{gridU}/e2u"].masked.data, self[f"{gridU}/e3u"].masked.data
        e1v, e3v = self[f"{gridV}/e1v"].masked.data, self[f"{gridV}/e3v"].masked.data  

        # -- Calculate divergence on T-points -- #
        dvar_i = (e2u * e3u * da_i).diff(dim=i_name, label="lower")
        dvar_i.coords[i_name] = dvar_i.coords[i_name] + 0.5

        dvar_j = (e1v * e3v * da_j).diff(dim=j_name, label="lower")
        dvar_j.coords[j_name] = dvar_j.coords[j_name] + 0.5

        divergence = (1 / (e1t * e2t * e3t)) * (dvar_i + dvar_j).where(tmask)

        # -- Update DataArray properties -- #
        divergence.name = f"div({var_i}, {var_j})"
        divergence = divergence.drop_vars(
            [
                f"{dom_prefix}glamu",
                f"{dom_prefix}gphiu",
                f"{dom_prefix}glamv",
                f"{dom_prefix}gphiv",
                f"{dom_prefix}depthu",
                f"{dom_prefix}depthv",
            ]
        )

        return divergence

    @deprecated(version_since="2026.03.b1",
                version_removed="2026.05",
                alternative="NEMOVectorField.curl from v2026.04 onwards"
                )
    def curl(
        self,
        vars: list[str],
        dom: str = ".",
    ) -> xr.DataArray:
        """
        Calculate the vertical (k) curl component of a vector field on a NEMO model grid.

        Parameters
        ----------
        vars : list[str]
            Name of the vector variables, structured as: ['u', 'v'], where 'u' and 'v' are
            the i and j components of the vector field, respectively.
        dom : str, optional
            Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
            Default is '.' for the parent domain.

        Returns
        -------
        xr.DataArray
            Vertical curl component of vector field defined on a NEMO model grid.

        Examples
        --------
        Compute the vertical component of the curl of the seawater velocity field in
        the second NEMO nested child domain:

        >>> nemo.curl(dom="2", vars=["uo", "vo"])

        Note, `vars` expects a list of the `i` and `j` components of the vector field,
        respectively.

        See Also
        --------
        divergence
        """
        # -- Validate input -- #
        if not isinstance(vars, list) or len(vars) != 2:
            raise ValueError(
                "vars must be a list of two elements structured as ['u', 'v']."
            )
        if not isinstance(dom, str):
            raise ValueError(
                "dom must be a string specifying the prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
            )

        # -- Get NEMO model grid properties -- #
        dom_prefix, _ = self._get_properties(dom=dom)
        grid_paths = self._get_grid_paths(dom=dom)
        gridU, gridV, gridF = (
            grid_paths["gridU"],
            grid_paths["gridV"],
            grid_paths["gridF"],
        )
        ijk_names = self._get_ijk_names(dom=dom)
        i_name, j_name = ijk_names["i"], ijk_names["j"]

        # -- Define i,j vector components -- #
        var_i, var_j = vars[0], vars[1]
        if var_i not in self[gridU].data_vars:
            raise KeyError(f"variable '{var_i}' not found in grid '{gridU}'.")
        if var_j not in self[gridV].data_vars:
            raise KeyError(f"variable '{var_j}' not found in grid '{gridV}'.")

        da_i = self[f"{gridU}/{var_i}"].masked.data
        da_j = self[f"{gridV}/{var_j}"].masked.data

        # -- Collect mask -- #
        if (f"{dom_prefix}depthu" in da_i.coords) and (
            f"{dom_prefix}depthv" in da_j.coords
        ):
            # 3-dimensional fmask
            fmask = self[gridF]["fmask"]
        else:
            # 2-dimensional fmask (unique points):
            fmask = self[gridF]["fmaskutil"]

        # -- Neglecting the final F-grid points along i, j dimensions -- #
        e1f = self[f"{gridF}/e1f"].masked.data.isel(
            {i_name: slice(None, -1), j_name: slice(None, -1)}
        )
        e2f = self[f"{gridF}/e2f"].masked.data.isel(
            {i_name: slice(None, -1), j_name: slice(None, -1)}
        )

        e1u = self[f"{gridU}/e1u"].masked.data
        e2v = self[f"{gridV}/e2v"].masked.data
        # -- Calculate vertical curl component on F-points -- #
        dvar_i = (e2v * da_j).diff(dim=i_name, label="lower")
        dvar_i.coords[i_name] = dvar_i.coords[i_name] + 0.5

        dvar_j = (e1u * da_i).diff(dim=j_name, label="lower")
        dvar_j.coords[j_name] = dvar_j.coords[j_name] + 0.5

        curl = (1 / (e1f * e2f)) * (dvar_i - dvar_j).where(fmask)

        # -- Update DataArray properties -- #
        curl.name = f"curl({var_i}, {var_j})"
        curl = curl.drop_vars(
            [
                f"{dom_prefix}glamu",
                f"{dom_prefix}gphiu",
                f"{dom_prefix}glamv",
                f"{dom_prefix}gphiv",
            ]
        )

        return curl

    @deprecated(version_since="2026.03.b1",
                version_removed="2026.05",
                alternative="NEMODataArray.integral"
                )
    def integral(
        self,
        grid: str,
        var: str,
        dims: list,
        cum_dims: list | None = None,
        dir: str | None = None,
        mask: xr.DataArray | None = None,
    ) -> xr.DataArray:
        """
        Integrate a variable along one or more dimensions of a NEMO model grid.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid where variable is stored (e.g., 'gridT').
        var : str
            Name of variable to integrate.
        dims : list
            Dimensions over which to integrate (e.g., ['i', 'k']).
        cum_dims : list, optional
            Dimensions over which to cumulatively integrate (e.g., ['k']).
            Specified dimensions must also be included in `dims`.
        dir : str, optional
            Direction of cumulative integration. Options are '+1' (along
            increasing cum_dims) or '-1' (along decreasing cum_dims).
        mask: xr.DataArray, optional
            Boolean mask identifying NEMO model grid points to be included (1)
            or neglected (0) from integration.

        Returns
        -------
        xr.DataArray
            Variable integrated along specified dimensions of the NEMO model grid.


        Examples
        --------
        Compute the integral of conservative temperature `thetao_con` along the vertical
        `k` dimension in the NEMO parent domain:


        >>> nemo.integral(grid="gridT",
        ...               var="thetao_con",
        ...               dims=["k"]
        ...               )

        Compute the vertical meridional overturning stream function from the meridional
        velocity `vo` (zonally integrated meridional velocity accumulated with increasing
        depth):

        >>> nemo.integral(grid="gridV",
        ...               var="vo",
        ...               dims=["i", "k"],
        ...               cum_dims=["k"],
        ...               dir="+1",
        ...               )

        See Also
        --------
        gradient
        """
        # -- Validate input -- #
        grid_keys = list(dict(self.subtree_with_keys).keys())
        if grid not in grid_keys:
            raise KeyError(
                f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
            )
        if var not in self[grid].data_vars:
            raise KeyError(f"variable '{var}' not found in grid '{grid}'.")
        if cum_dims is not None:
            for dim in cum_dims:
                if dim not in dims:
                    raise ValueError(
                        f"cumulative integration dimension '{dim}' not included in `dims`."
                    )
            if dir not in ["+1", "-1"]:
                raise ValueError(
                    f"invalid direction of cumulative integration '{dir}'. Expected '+1' or '-1'."
                )
        if mask is not None:
            if not isinstance(mask, xr.DataArray):
                raise ValueError("mask must be an xarray.DataArray.")
            if any(dim not in self[grid].dims for dim in mask.dims):
                raise ValueError(
                    f"mask must have dimensions subset from {self[grid].dims}."
                )

        # -- Collect variable, weights & mask -- #
        da = (
            self[f"{grid}/{var}"].masked.data.where(mask)
            if mask is not None
            else self[f"{grid}/{var}"].masked.data
        )
        weights = self._get_weights(grid=grid, dims=dims)

        # -- Perform integration -- #
        if cum_dims is not None:
            sum_dims = [dim for dim in dims if dim not in cum_dims]
            if dir == "+1":
                # Cumulative integration along ordered dimension:
                result = (
                    da.weighted(weights)
                    .sum(dim=sum_dims, skipna=True)
                    .cumsum(dim=cum_dims, skipna=True)
                )
            elif dir == "-1":
                # Cumulative integration along reversed dimension:
                result = (
                    da.weighted(weights)
                    .sum(dim=sum_dims, skipna=True)
                    .reindex({dim: self[grid][dim][::-1] for dim in cum_dims})
                    .cumsum(dim=cum_dims, skipna=True)
                )
        else:
            # Integration only:
            result = da.weighted(weights).sum(dim=dims, skipna=True)

        # -- Update DataArray properties -- #
        result.name = f"integral_{', '.join(dims)}({var})"

        return result

    @deprecated(version_since="2026.03.b1",
                version_removed="2026.05",
                alternative="NEMODataArray.depth_integral"
                )
    def depth_integral(
        self, grid: str, var: str, limits: tuple[int | float]
    ) -> xr.Dataset:
        """
        Integrate a variable in depth coordinates between two limits.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid where variable is stored (e.g., 'gridT').
        var : str
            Name of the variable to vertically integrate.
        limits : tuple[int | float]
            Limits of depth integration given as a tuple of the form
            (depth_lower, depth_upper) where depth_lower and depth_upper are
            the lower and upper limits of vertical integration, respectively.

        Returns
        -------
        xr.DataArray
            Vertical integral of chosen variable between depth surfaces (depth_lower, depth_upper).

        Examples
        --------
        Vertically integrate the conservative temperature variable `thetao_con` defined in a
        NEMO model parent domain from the sea surface to 100 m depth:

        >>> nemo.depth_integral(grid='gridT',
        ...                     var='thetao_con',
        ...                     limits=(0, 100)
        ...                              )

        See Also
        --------
        integral
        """
        # -- Validate input -- #
        grid_keys = list(dict(self.subtree_with_keys).keys())
        if grid not in grid_keys:
            raise KeyError(
                f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
            )
        if var not in self[grid].data_vars:
            raise KeyError(f"Variable '{var}' not found in grid '{grid}'.")
        if (not isinstance(limits, tuple)) | (len(limits) != 2):
            raise TypeError(
                "depth limits of integration should be given by a tuple of the form (depth_lower, depth_upper)"
            )
        if (limits[0] < 0) | (limits[1] < 0):
            raise ValueError("depth limits of integration must be non-negative.")
        if limits[0] >= limits[1]:
            raise ValueError(
                "lower depth limit must be less than upper depth limit."
            )

        # -- Get NEMO model grid properties -- #
        dom, _, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
        ijk_names = self._get_ijk_names(dom=dom)
        i_name, j_name, k_name = ijk_names["i"], ijk_names["j"], ijk_names["k"]

        # -- Define input variables -- #
        var_in = self[f"{grid}/{var}"].masked.data
        e3_in = self[f"{grid}/e3{grid_suffix}"].masked.data

        # -- Vertically integrate w.r.t depth -- #
        result = xr.apply_ufunc(
            compute_depth_integral,
            e3_in,
            var_in,
            np.array([limits[1]]),
            np.array([limits[0]]),
            input_core_dims=[[k_name], [k_name], [None], [None]],
            output_core_dims=[["k_new"]],
            dask="parallelized",
            output_dtypes=[var_in.dtype],
            dask_gufunc_kwargs={"output_sizes": {"k_new": 1}},
        )

        # -- Create variable integral DataArray -- #
        t_name = self[f"{grid}/{var}"].t_name
        if t_name is not None:
            result = result.transpose(t_name, "k_new", j_name, i_name).squeeze()
        else:
            result = result.transpose("k_new", j_name, i_name).squeeze()
        result.name = f"integral_z({var})"

        return result

    def clip_grid(
        self,
        grid: str,
        bbox: tuple,
    ) -> Self:
        """
        Clip a NEMO model grid to specified longitude and latitude range.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid to clip (e.g., 'gridT').
        bbox : tuple
            Bounding box to clip to (lon_min, lon_max, lat_min, lat_max).

        Returns
        -------
        NEMODataTree
            NEMO DataTree with specified model grid clipped to bounding box.

        Examples
        --------
        Clip T-grid in a NEMO parent domain in the bounding box (-80°E, 0°E, 40°N, 80°N):

        >>> bbox = (-80, 0, 40, 80)

        >>> nemo.clip_grid(grid="gridT", bbox=bbox)

        See Also
        --------
        clip_domain
        """
        # -- Validate input -- #
        if not isinstance(bbox, tuple) or len(bbox) != 4:
            raise ValueError(
                "bounding box must be a tuple (lon_min, lon_max, lat_min, lat_max)."
            )

        # -- Get NEMO model grid properties -- #
        _, dom_prefix, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)

        # -- Clip the grid to given bounding box -- #
        # Indexing with a mask requires loading coords into memory:
        glam = self[grid][f"{dom_prefix}glam{grid_suffix}"].load()
        gphi = self[grid][f"{dom_prefix}gphi{grid_suffix}"].load()

        grid_clipped = self[grid].dataset.where(
            (glam >= bbox[0])
            & (glam <= bbox[1])
            & (gphi >= bbox[2])
            & (gphi <= bbox[3]),
            drop=True,
        )

        d_dtypes = {var: self[grid][var].dtype for var in self[grid].dataset.data_vars}
        for var, dtype in d_dtypes.items():
            if dtype in [np.int32, np.int64, bool]:
                grid_clipped[var] = grid_clipped[var].fillna(0).astype(dtype)

        if bbox != (-180, 180, -90, 90):
            # Update zonal periodicity of grid node:
            grid_clipped = grid_clipped.assign_attrs({"iperio": False})

        # Update shallow copy of NEMODataTree:
        self_copy = self.copy()
        self_copy[grid].dataset = grid_clipped

        return self_copy

    def clip_domain(
        self,
        dom: str,
        bbox: tuple,
    ) -> Self:
        """
        Clip a NEMO model domain to specified longitude and latitude range.

        Parameters
        ----------
        dom : str
            Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
            Default is '.' for the parent domain.
        bbox : tuple
            Bounding box to clip to (lon_min, lon_max, lat_min, lat_max).

        Returns
        -------
        NEMODataTree
            NEMO DataTree with specified model domain clipped to bounding box.

        Examples
        --------
        Clip all model grids in a NEMO parent domain in the bounding box
        (-80°E, 0°E, 40°N, 80°N):

        >>> bbox = (-80, 0, 40, 80)

        >>> nemo.clip_domain(dom=".", bbox=bbox)

        See Also
        --------
        clip_grid
        """
        # -- Validate input -- #
        if not isinstance(dom, str):
            raise ValueError(
                "dom must be a string specifying the prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
            )
        if not isinstance(bbox, tuple) or len(bbox) != 4:
            raise ValueError(
                "bounding box must be a tuple: (lon_min, lon_max, lat_min, lat_max)."
            )

        # -- Get NEMO model grid properties -- #
        grid_paths = self._get_grid_paths(dom=dom)
        ijk_names = self._get_ijk_names(dom=dom)
        i_name, j_name = ijk_names["i"], ijk_names["j"]

        # -- Clip grids to given bounding box -- #
        if not grid_paths:
            raise ValueError(f"NEMO model domain '{dom}' not found in the DataTree.")
        else:
            for grid in grid_paths.values():
                # Identify grid type:
                grid_suffix = self._get_properties(grid=grid)

                if grid_suffix == "t":
                    # Clip shallow copy of NEMODataTree T-grid using lon/lat bbox:
                    self_copy = self.copy().clip_grid(grid=grid, bbox=bbox)
                    # Store (i, j) coords of bbox on T-grid:
                    i_bbox = self_copy[grid][i_name]
                    j_bbox = self_copy[grid][j_name]

                else:
                    # Clip adjacent horizontal grid using (i, j) coords of clipped T-grid:
                    match grid_suffix:
                        case "u":
                            grid_clipped = self[grid].dataset.sel(i=i_bbox + 0.5, j=j_bbox)
                        case "v":
                            grid_clipped = self[grid].dataset.sel(i=i_bbox, j=j_bbox + 0.5)
                        case "w":
                            grid_clipped = self[grid].dataset.sel(i=i_bbox, j=j_bbox)
                        case "f":
                            grid_clipped = self[grid].dataset.sel(i=i_bbox + 0.5, j=j_bbox + 0.5)

                    if bbox != (-180, 180, -90, 90):
                        # Update zonal periodicity of NEMO model grid:
                        grid_clipped = grid_clipped.assign_attrs({"iperio": False})
                    # Update shallow copy of NEMODataTree:
                    self_copy[grid].dataset = grid_clipped

        return self_copy

    def mask_with_polygon(
        self,
        grid: str,
        lon_poly: list | np.ndarray,
        lat_poly: list | np.ndarray,
    ):
        """
        Create mask of NEMO model grid points contained within a polygon.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid where longitude and latitude coordinates
            are stored (e.g., 'gridT').
        lon_poly : list | ndarray
            Longitudes of closed polygon.
        lat_poly : list | ndarray
            Latitudes of closed polygon.

        Returns
        -------
        xr.DataArray
            Boolean mask identifying NEMO model grid points which are inside
            the polygon.

        Examples
        --------
        Create a regional boolean mask using the geographical coordinates of a closed
        polygon `lon_poly` and `lat_poly` in a NEMO parent domain:


        >>> nemo.mask_with_polygon(grid="gridT",
        ...                        lon_poly=lon_poly,
        ...                        lat_poly=lat_poly,
        ...                        )

        See Also
        --------
        masked_statistic
        """
        # -- Validate input -- #
        if not isinstance(lon_poly, (np.ndarray, list)) or not isinstance(
            lat_poly, (np.ndarray, list)
        ):
            raise TypeError(
                "longitude & latitude coordinates of polygon must be numpy arrays or lists."
            )
        if (lon_poly[0] != lon_poly[-1]) or (lat_poly[0] != lat_poly[-1]):
            raise ValueError(
                "longitude & latitude coordinates must form a closed polygon."
            )

        # -- Get NEMO model grid properties -- #
        dom, dom_prefix, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
        ijk_names = self._get_ijk_names(grid=grid)
        i_name, j_name = ijk_names["i"], ijk_names["j"]

        if dom == ".":
            lon_name = f"glam{grid_suffix}"
            lat_name = f"gphi{grid_suffix}"
        else:
            lon_name = f"{dom_prefix}glam{grid_suffix}"
            lat_name = f"{dom_prefix}gphi{grid_suffix}"

        # -- Create mask using polygon coordinates -- #
        mask = create_polygon_mask(
            lon_grid=self[grid][lon_name],
            lat_grid=self[grid][lat_name],
            lon_poly=lon_poly,
            lat_poly=lat_poly,
            dims=(j_name, i_name),
        )

        return mask

    @deprecated(version_since="2026.03.b1",
                version_removed="2026.05",
                alternative="NEMODataArray.masked_statistic"
                )
    def masked_statistic(
        self,
        grid: str,
        var: str,
        lon_poly: list | np.ndarray,
        lat_poly: list | np.ndarray,
        statistic: str,
        dims: list,
    ) -> xr.DataArray:
        """
        Masked statistic of a variable defined on a NEMO model grid.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid where variable is stored (e.g., 'gridT').
        var : str
            Name of the variable to compute statistic.
        lon_poly : list | np.ndarray
            Longitudes of closed polygon.
        lat_poly : list | np.ndarray
            Latitudes of closed polygon.
        statistic : str
            Name of the statistic to calculate (e.g., 'mean', 'weighted_mean' 'sum').
        dims : list
            Dimensions over which to apply statistic (e.g., ['i', 'j']).

        Returns
        -------
        xr.DataArray
            Masked statistic of specified variable.

        Examples
        --------
        Compute the grid cell area-weighted mean sea surface temperature `tos_con` for a
        region enclosed in a polygon defined by `lon_poly` and `lat_poly` in a NEMO nested
        child domain:

        >>> nemo.masked_statistic(grid="gridT/1_gridT",
        ...                       var="tos_con",
        ...                       lon_poly=lon_poly,
        ...                       lat_poly=lat_poly,
        ...                       statistic="weighted_mean",
        ...                       dims=["i", "j"]
        ...                       )

        See Also
        --------
        binned_statistic
        """
        # -- Validate input -- #
        grid_keys = list(dict(self.subtree_with_keys).keys())
        if grid not in grid_keys:
            raise KeyError(
                f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
            )
        if var not in self[grid].data_vars:
            raise KeyError(f"variable '{var}' not found in grid '{grid}'.")

        # -- Create polygon mask using coordinates -- #
        mask_poly = self.mask_with_polygon(
            lon_poly=lon_poly, lat_poly=lat_poly, grid=grid
        )

        # -- Get NEMO model grid properties -- #
        _, _, dom_suffix, _ = self._get_properties(grid=grid, infer_dom=True)

        # -- Apply masks & calculate statistic -- #
        da = self[f"{grid}/{var}"].masked.data.where(mask_poly)

        match statistic:
            case "mean":
                result = da.mean(dim=dims, skipna=True)

            case "weighted_mean":
                weight_dims = [dim.replace(dom_suffix, "") for dim in dims]
                weights = self._get_weights(grid=grid, dims=weight_dims)
                result = da.weighted(weights).mean(dim=dims, skipna=True)

            case "min":
                result = da.min(dim=dims, skipna=True)

            case "max":
                result = da.max(dim=dims, skipna=True)

            case "sum":
                result = da.sum(dim=dims, skipna=True)

            case _:
                raise ValueError(
                    f"Unsupported statistic '{statistic}'. Supported statistics are: 'mean', 'weighted_mean', 'min', 'max', 'sum'."
                )

        return result

    def extract_mask_boundary(
        self,
        mask: xr.DataArray,
        uv_vars: list | None = None,
        vars: list | None = None,
        dom: str = ".",
    ) -> xr.Dataset:
        """
        Extract the boundary of a masked region defined on a NEMO model grid.

        Parameters
        ----------
        mask : xr.DataArray
            Boolean mask identifying NEMO model grid points which
            are inside the region of interest.
        uv_vars : list, optional
            Names of velocity variables to extract along the boundary.
            Default is ['uo', 'vo'].
        vars : list, optional
            Names of scalar variables to extract along the boundary.
        dom : str, optional
            Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
            Default is '.' for the parent domain.

        Returns
        -------
        xr.Dataset
            Dataset containing variables and NEMO model coordinates
            extracted along the boundary of the mask.

        Examples
        --------
        Extract normal velocities and absolute salinity along the boundary of
        a masked region in the NEMO parent domain:

        >>> nemo.extract_mask_boundary(mask=mask_osnap,
        ...                            uv_vars=["uo", "vo"],
        ...                            vars=["so_abs"],
        ...                            dom=".",
        ...                            )

        See Also
        --------
        extract_section
        """
        # -- Validate Input -- #
        if not isinstance(mask, xr.DataArray):
            raise TypeError("mask must be an xarray DataArray")
        if not isinstance(dom, str):
            raise TypeError(
                "dom must be a string specifying prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
            )
        if uv_vars is None:
            uv_vars = ["uo", "vo"]
        else:
            if not isinstance(uv_vars, list) or len(uv_vars) != 2:
                raise TypeError(
                    "uv_vars must be a list of velocity variables to extract (e.g., ['uo', 'vo'])."
                )

        # -- Get NEMO model grid properties -- #
        _, dom_suffix = self._get_properties(dom=dom)

        # -- Extract mask boundary -- #
        if f"i{dom_suffix}" not in mask.dims or f"j{dom_suffix}" not in mask.dims:
            raise ValueError(
                f"mask must have dimensions 'i{dom_suffix}' and 'j{dom_suffix}'"
            )
        i_bdy, j_bdy, flux_type, flux_dir = get_mask_boundary(mask)

        # -- Construct boundary dataset -- #
        # Neglecting final indices -> duplicate of the first indices:
        ds_bdy = create_boundary_dataset(
            nemo=self,
            dom=dom,
            i_bdy=i_bdy[:-1],
            j_bdy=j_bdy[:-1],
            flux_type=flux_type[:-1],
            flux_dir=flux_dir[:-1],
        )

        # -- Update boundary dataset with extracted section data -- #
        ds_bdy = update_boundary_dataset(
            ds_bdy=ds_bdy,
            nemo=self,
            dom=dom,
            sec_indexes=None,
            uv_vars=uv_vars,
            vars=vars,
        )

        return ds_bdy

    def extract_section(
        self,
        lon_section: np.ndarray,
        lat_section: np.ndarray,
        uv_vars: list | None = None,
        vars: list | None = None,
        dom: str = ".",
    ) -> xr.Dataset:
        """
        Extract hydrographic section from a NEMO model domain.

        Parameters
        ----------
        lon_section : np.ndarray
            Longitudes defining the section polygon.
        lat_section : np.ndarray
            Latitudes defining the section polygon.
        uv_vars : list, optional
            Names of velocity variables to extract along the boundary.
            Default is ['uo', 'vo'].
        vars : list, optional
            Names of scalar variables to extract along the boundary.
        dom : str
            Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
            Default is '.' for the parent domain.

        Returns
        -------
        xr.Dataset
            Dataset containing hydrographic section extracted from NEMO model grid.

        Examples
        --------
        Extract normal velocities and potential density along the Overturning in the Subpolar
        North Atlantic (OSNAP) array defined by `lon_osnap` and `lat_osnap` coordinates in the
        NEMO parent domain:

        >>> nemo.extract_section(lon_section=lon_osnap,
        ...                      lat_section=lat_osnap,
        ...                      uv_vars=["uo", "vo"],
        ...                      vars=["sigma0"],
        ...                      dom=".",
        ...                      )

        See Also
        --------
        extract_mask_boundary
        """
        # -- Validate Input -- #
        if not isinstance(lon_section, np.ndarray):
            raise TypeError("lon_section must be a numpy array.")
        if not isinstance(lat_section, np.ndarray):
            raise TypeError("lat_section must be a numpy array.")
        if not isinstance(dom, str):
            raise TypeError(
                "dom must be a string specifying prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
            )
        if uv_vars is None:
            uv_vars = ["uo", "vo"]
        else:
            if not isinstance(uv_vars, list) or len(uv_vars) != 2:
                raise TypeError(
                    "uv_vars must be a list of velocity variables to extract (e.g., ['uo', 'vo'])."
                )

        # -- Get NEMO model grid properties -- #
        grid_paths = self._get_grid_paths(dom=dom)

        # -- Define hydrographic section using polygon -- #
        lon_poly, lat_poly = create_section_polygon(
            lon_sec=lon_section,
            lat_sec=lat_section,
        )

        mask = self.mask_with_polygon(
            grid=grid_paths["gridT"], lon_poly=lon_poly, lat_poly=lat_poly
        )

        i_bdy, j_bdy, flux_type, flux_dir = get_mask_boundary(mask)

        # -- Create mask boundary dataset -- #
        ds_bdy = create_boundary_dataset(
            nemo=self,
            dom=dom,
            i_bdy=i_bdy,
            j_bdy=j_bdy,
            flux_type=flux_type,
            flux_dir=flux_dir,
        )

        # -- Get indexes of hydrographic section along mask boundary -- #
        dom_prefix, _ = self._get_properties(dom=dom)
        sec_indexes = get_section_indexes(
            lon_section=lon_section,
            lat_section=lat_section,
            gphib=ds_bdy[f"{dom_prefix}gphib"].values,
            glamb=ds_bdy[f"{dom_prefix}glamb"].values,
            bdy=ds_bdy["bdy"].values,
        )

        # -- Update boundary dataset with extracted section data -- #
        ds_bdy = update_boundary_dataset(
            ds_bdy=ds_bdy,
            nemo=self,
            dom=dom,
            sec_indexes=sec_indexes,
            uv_vars=uv_vars,
            vars=vars,
        )

        return ds_bdy

    def binned_statistic(
        self,
        grid: str,
        vars: list[str],
        values: str,
        keep_dims: list[str] | None,
        bins: list[list | np.ndarray],
        statistic: str,
        mask: xr.DataArray | None,
    ) -> xr.DataArray:
        """
        Calculate binned statistic of a variable defined on a NEMO model grid.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid where variables and values are stored
            (e.g., 'gridT').
        vars : list[str]
            Names of variable(s) to be grouped in discrete bins.
        values : str
            Name of the values with which to calculate binned statistic.
        keep_dims : list[str] | None
            Names of dimensions in values to keep as labels in binned statistic.
        bins : list[list | np.ndarray]
            Bin edges used to group each of the variables in `vars`.
        statistic : str
            Statistic to calculate (e.g., 'count', 'sum', 'nansum', 'mean', 'nanmean',
            'max', 'nanmax', 'min', 'nanmin'). See flox.xarray.xarray_reduce for a
            complete list of aggregation statistics.
        mask : xr.DataArray | None
            Boolean mask identifying NEMO model grid points to be included (1)
            or neglected (0) from calculation.

        Returns
        -------
        xr.DataArray
            Values of the selected statistic in each bin.

        Examples
        --------
        Compute the mean depth associated with each isopycnal in discrete potential
        density `sigma0` coordinates:

        >>> sigma0_bins = np.arange(22, 29.05, 0.05)

        >>> nemo.binned_statistic(grid="gridT",
        ...                       vars=["sigma0"],
        ...                       values="deptht",
        ...                       keep_dims=["time_counter"],
        ...                       bins=[sigma0_bins],
        ...                       statistic="nanmean",
        ...                       )

        See Also
        --------
        masked_statistic
        """
        # -- Validate input -- #
        grid_keys = list(dict(self.subtree_with_keys).keys())
        if grid not in grid_keys:
            raise KeyError(
                f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
            )
        if any(var not in self[grid].data_vars for var in vars):
            raise KeyError(f"one or more variables {vars} not found in grid '{grid}'.")
        if values not in self[grid].data_vars:
            raise KeyError(f"values '{values}' not found in grid '{grid}'.")
        if keep_dims is not None:
            if any(dim not in self[grid][values].dims for dim in keep_dims):
                raise KeyError(
                    f"one or more dimensions {keep_dims} not found in values '{values}'."
                )
        if not all(isinstance(bin, (list, np.ndarray)) for bin in bins):
            raise ValueError("bins must be a list of lists or numpy arrays.")
        if statistic not in [
            "all",
            "any",
            "count",
            "sum",
            "nansum",
            "mean",
            "nanmean",
            "max",
            "nanmax",
            "min",
            "nanmin",
            "argmax",
            "nanargmax",
            "argmin",
            "nanargmin",
            "quantile",
            "nanquantile",
            "median",
            "nanmedian",
            "mode",
            "nanmode",
            "first",
            "nanfirst",
            "last",
            "nanlast",
        ]:
            raise ValueError(f"statistic '{statistic}' is not supported.")
        if mask is not None:
            if not isinstance(mask, xr.DataArray):
                raise ValueError("mask must be an xarray.DataArray.")
            if mask.dtype != bool:
                raise TypeError("mask dtype must be boolean.")
            if any(dim not in self[grid].dims for dim in mask.dims):
                raise ValueError(
                    f"mask must have dimensions subset from {self[grid].dims}."
                )

        # -- Calculate binned statistics -- #
        values_data = self[f"{grid}/{values}"].masked.data
        var_data = [self[f"{grid}/{var}"].masked.data for var in vars]

        result = compute_binned_statistic(
            vars=var_data,
            values=values_data,
            keep_dims=keep_dims,
            bins=bins,
            statistic=statistic,
            mask=mask,
        )

        return result

    @deprecated(version_since="2026.03.b1",
                version_removed="2026.05",
                alternative="NEMODataArray.transform_vertical_grid"
                )
    def transform_vertical_grid(
        self, grid: str, var: str, e3_new: xr.DataArray
    ) -> xr.Dataset:
        """
        Transform variable defined on a NEMO model grid to a new vertical grid using conservative interpolation.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid where variable is stored
            (e.g., 'gridT').
        var : str
            Name of the variable to transform.
        e3_new : xarray.DataArray
            Grid cell thicknesses of the new vertical grid.
            Must be a 1-dimensional xarray.DataArray with
            dimension 'k_new'.

        Returns
        -------
        xr.Dataset
            Variable defined at the centre of each vertical
            grid cell on the new grid, and vertical grid cell
            thicknesses adjusted for model bathymetry.

        Examples
        --------
        Transform the conservative temperature variable `thetao_con` defined in a
        NEMO model parent domain from it's native 75 unevenly-spaced z-levels to
        regularly spaced z-levels at 200 m intervals:

        >>> e3t_target = xr.DataArray(np.repeat(200.0, 30), dims=['k_new'])

        >>> nemo.transform_vertical_grid(grid='gridT',
        ...                              var='thetao_con',
        ...                              e3_new=e3t_target
        ...                              )

        See Also
        --------
        transform_to
        """
        # -- Validate input -- #
        grid_keys = list(dict(self.subtree_with_keys).keys())
        if grid not in grid_keys:
            raise KeyError(
                f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
            )
        if var not in self[grid].data_vars:
            raise KeyError(f"Variable '{var}' not found in grid '{grid}'.")
        if e3_new.dims != ("k_new",) or (e3_new.ndim != 1):
            raise ValueError(
                "e3_new must be a 1-dimensional xarray.DataArray with dimension 'k_new'."
            )

        # -- Get NEMO model grid properties -- #
        dom, _, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
        ijk_names = self._get_ijk_names(dom=dom)
        i_name, j_name, k_name = ijk_names["i"], ijk_names["j"], ijk_names["k"]
        t_name = [dim for dim in self[grid].dims if "time" in dim][0]

        # -- Define input variables -- #
        var_in = self[f"{grid}/{var}"].masked.data
        e3_in = self[f"{grid}/e3{grid_suffix}"].masked.data
        if e3_new.sum(dim="k_new") < self[grid][f"depth{grid_suffix}"].max(dim=k_name):
            raise ValueError(
                f"e3_new must sum to at least the maximum depth ({self[grid][f'depth{grid_suffix}'].max(dim=k_name).item()} m) of the original vertical grid."
            )

        # -- Transform variable to target vertical grid -- #
        var_out, e3_out = xr.apply_ufunc(
            transform_vertical_coords,
            e3_in,
            var_in,
            e3_new.astype(e3_in.dtype),
            input_core_dims=[[k_name], [k_name], ["k_new"]],
            output_core_dims=[["k_new"], ["k_new"]],
            dask="parallelized",
            output_dtypes=[var_in.dtype, e3_in.dtype],
            dask_gufunc_kwargs={"output_sizes": {"k_new": e3_new.sizes["k_new"]}},
        )

        # -- Construct transformed variable Dataset -- #
        var_out = var_out.transpose(t_name, "k_new", j_name, i_name)

        ds_out = xr.Dataset(
            data_vars={var: var_out, f"e3{grid_suffix}_new": e3_out},
            coords={
                f"depth{grid_suffix}_new": ("k_new", e3_new.cumsum(dim="k_new").data)
            },
        )

        return ds_out

    @deprecated(version_since="2026.03.b1",
                version_removed="2026.05",
                alternative="NEMODataArray.interp_to"
                )
    def transform_to(
        self,
        grid: str,
        var: str,
        to: str,
    ) -> xr.DataArray:
        """
        Transform variable defined on a NEMO model grid to a neighbouring
        horizontal grid using linear interpolation.

        For flux variables defined at U- or V-points, the specified variable
        is first weighted by grid cell face areas prior to linear interpolation,
        and is then normalised by the target grid cell face areas following
        interpolation.

        Parameters
        ----------
        grid : str
            Path to NEMO model grid where variable is stored (e.g., 'gridT').
        var : str
            Name of the variable to transform.
        to : str
            Suffix of the neighbouring horizontal NEMO model grid to
            transform variable to. Options are 'T', 'U', 'V', 'F'.

        Returns
        -------
        xr.DataArray
            Values of variable linearly interpolated onto a neighbouring
            horizontal grid.

        Examples
        --------
        Transform conservative temperature `thetao_con` defined on scalar T-points
        to neighbouring V-points in a NEMO model parent domain:

        >>> nemo.transform_to(grid='gridT', var='thetao_con', to='V')

        See Also
        --------
        transform_vertical_grid
        """
        # -- Validate input -- #
        grid_keys = list(dict(self.subtree_with_keys).keys())
        if grid not in grid_keys:
            raise KeyError(
                f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
            )
        if var not in self[grid].data_vars:
            raise KeyError(f"variable '{var}' not found in grid '{grid}'.")
        if not isinstance(to, str):
            raise TypeError(f"'to' must be a string, got {type(to)}.")
        if to not in ["T", "U", "V", "F"]:
            raise ValueError(f"'to' must be one of ['T', 'U', 'V', 'F'], got {to}.")

        # -- Get NEMO model grid properties -- #
        _, dom_prefix, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
        ijk_names = self._get_ijk_names(grid=grid)
        i_name, j_name, k_name = ijk_names["i"], ijk_names["j"], ijk_names["k"]
        iperio = self[grid].attrs.get("iperio", False)
        target_grid = f"{grid.replace(grid[-1], to)}"

        # -- Prepare variable for linear interpolation -- #
        if grid_suffix.upper() in ["U", "V"]:
            weight_dims = (
                [k_name, j_name] if grid_suffix.upper() == "U" else [k_name, i_name]
            )
            if f"{dom_prefix}depth{grid_suffix}" in self[grid][var].coords:
                # 3-D variables - weight by grid cell face area:
                weights = self._get_weights(grid=grid, dims=weight_dims, fillna=False)
                target_weights = self._get_weights(
                    grid=target_grid, dims=weight_dims, fillna=False
                )
            else:
                # 2-D variables - weight by grid cell width:
                weights = self._get_weights(grid=grid, dims=weight_dims[1], fillna=False)
                target_weights = self._get_weights(
                    grid=target_grid, dims=weight_dims[1], fillna=False
                )
            da = self[f"{grid}/{var}"].masked.data * weights
        else:
            # Scalar variables:
            da = self[f"{grid}/{var}"].masked.data

        # -- Linearly interpolate variable -- #
        # Define source grid mask:
        if f"{dom_prefix}depth{grid_suffix}" in da.coords:
            mask = self[grid][f"{grid_suffix}mask"]
        else:
            mask = self[grid][f"{grid_suffix}maskutil"]

        result = interpolate_grid(
            da=da,
            mask=mask,
            source_grid=grid[-1],
            target_grid=target_grid[-1],
            iperio=iperio,
            ijk_names=ijk_names,
        )

        # -- Update interpolated variable -- #
        # Retain input variable name:
        result.name = da.name
        # Reorder dimensions (time_counter, [k], j, i):
        new_dims = (result.dims[-1], *result.dims[:-1])
        result = result.transpose(*new_dims)

        # Update NEMO grid coords:
        result[i_name] = self[target_grid][i_name]
        result[j_name] = self[target_grid][j_name]
        if k_name in result.dims:
            result[k_name] = self[target_grid][k_name]

        # Drop NEMO source grid coords:
        drop_vars = [f"{dom_prefix}glam{grid_suffix}", f"{dom_prefix}gphi{grid_suffix}"]
        if f"{dom_prefix}depth{grid_suffix}" in da.coords:
            drop_vars.append(f"{dom_prefix}depth{grid_suffix}")
        result = result.drop_vars(drop_vars)

        # Normalise by target grid cell weights for flux variables:
        if grid_suffix.upper() in ["U", "V"]:
            result = result / target_weights

        # Apply target grid mask:
        if f"{dom_prefix}depth{grid_suffix}" in da.coords:
            target_mask = f"{to.lower()}mask"
        else:
            target_mask = f"{to.lower()}maskutil"
        result = result.where(self[target_grid][target_mask])

        return result

__getitem__

__getitem__(key)

Access child nodes, variables, or coordinates stored in this NEMODataTree.

Returned object will be either a NEMODataTree or NEMODataArray object depending on whether the key points to a child node or variable.

Overloads the getitem() method of xarray.DataTree to return NEMODataArrays accessed via variable paths (i.e, /grid/var).

Parameters:

Name	Type	Description	Default
key	str	Name of variable / child within this node, or unix-like path to variable / child within another node.	required

Returns:

Type	Description
NEMODataTree | NEMODataArray

Source code in nemo_cookbook/nemodatatree.py

def __getitem__(self, key: str) -> Self | NEMODataArray:
    """
    Access child nodes, variables, or coordinates stored in this NEMODataTree.

    Returned object will be either a NEMODataTree or NEMODataArray object depending
    on whether the key points to a child node or variable.

    Overloads the __getitem__() method of xarray.DataTree to return NEMODataArrays
    accessed via variable paths (i.e, /grid/var).

    Parameters
    ----------
    key : str
        Name of variable / child within this node, or unix-like path to variable
        / child within another node.

    Returns
    -------
    NEMODataTree | NEMODataArray
    """
    # -- Access child node or variable -- #
    item = super().__getitem__(key=key)
    is_gridpath = key.startswith("/grid") or key.startswith("grid")

    # -- Return NEMODataArray -- #
    if isinstance(item, xr.DataArray) and is_gridpath:
        var_name = key.split("/")[-1]
        grid = key.replace(f"/{var_name}", "")
        item = NEMODataArray(da=item,
                             tree=self,
                             grid=grid
                             )

    return item

__init__

__init__(*args, **kwargs)

Create a single node of a NEMODataTree.

The node may optionally contain data in the form of data and coordinate variables, stored in the same way as data is stored in an xarray.Dataset.

Parameters:

Name	Type	Description	Default
*args	tuple	Positional arguments to pass to the parent class.	()
**kwargs	dict	Keyword arguments to pass to the parent class.	{}

Source code in nemo_cookbook/nemodatatree.py

def __init__(self, *args, **kwargs):
    """
    Create a single node of a NEMODataTree.

    The node may optionally contain data in the form of data
    and coordinate variables, stored in the same way as data
    is stored in an `xarray.Dataset`.

    Parameters
    ----------
    *args : tuple
        Positional arguments to pass to the parent class.
    **kwargs : dict
        Keyword arguments to pass to the parent class.
    """
    super().__init__(*args, **kwargs)

__setitem__

__setitem__(key, value, strict=True)

Set a child node or variable in this NEMODataTree.

Overloads the setitem() method of xarray.DataTree to allow setting NEMODataArrays via variable paths (i.e, /grid/var).

Optionally set strict=False to bypass validation of child grid nodes.

Parameters:

Name	Type	Description	Default
key	str	Name of variable or child node, or unix-like path to variable within a child node.	required
value	NEMODataArray | DataArray | Dataset	Object to set at the specified key. If a NEMODataArray is provided, the underlying xarray.DataArray will be set at the specified key.	required
strict	bool	Validate Datasets assigned to NEMO grid nodes to ensure they contain the required dimensions and coordinates. Default is True.	True

Returns:

Type	Description
None

Source code in nemo_cookbook/nemodatatree.py

def __setitem__(
        self,
        key: str,
        value: NEMODataArray | xr.DataArray | xr.Dataset,
        strict: bool = True
    ) -> None:
    """
    Set a child node or variable in this NEMODataTree.

    Overloads the __setitem__() method of xarray.DataTree to allow
    setting NEMODataArrays via variable paths (i.e, /grid/var).

    Optionally set strict=False to bypass validation of child grid nodes.

    Parameters
    ----------
    key : str
        Name of variable or child node, or unix-like path to variable
        within a child node.

    value : NEMODataArray | xarray.DataArray | xarray.Dataset
        Object to set at the specified key. If a NEMODataArray is provided,
        the underlying xarray.DataArray will be set at the specified key.

    strict : bool, optional
        Validate Datasets assigned to NEMO grid nodes to ensure they contain
        the required dimensions and coordinates. Default is True.

    Returns
    -------
    None
    """
    # -- Access DataArray from NEMODataArray -- #
    if isinstance(value, NEMODataArray):
        value = value.data

    if strict and isinstance(value, xr.Dataset):
        # -- Validate NEMO grid node Dataset -- #
        validate_nemo_grid_node(key=key, value=value)

    return super().__setitem__(key, value)

add_geoindex

add_geoindex(grid)

Add geographical index variables to a given NEMO model grid.

This enables users to index grid variables using geographical coordinates (e.g., glamt, gphit) in addition to their (i, j, k) dimensions.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid to add geographical indexes (e.g., 'gridT').	required

Returns:

Type	Description
NEMODataTree	NEMO DataTree with geographical indexes added to specified model grid.

Examples:

Add glamt, gphit as geographical indexes to the T-grid of the NEMO parent domain:

>>> nemo.add_geoindex(grid="gridT")

Source code in nemo_cookbook/nemodatatree.py

def add_geoindex(
    self,
    grid: str,
) -> Self:
    """
    Add geographical index variables to a given NEMO model grid.

    This enables users to index grid variables using geographical
    coordinates (e.g., glamt, gphit) in addition to their (i, j, k)
    dimensions.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid to add geographical indexes (e.g., 'gridT').

    Returns
    -------
    NEMODataTree
        NEMO DataTree with geographical indexes added to specified model grid.

    Examples
    --------
    Add glamt, gphit as geographical indexes to the T-grid of the NEMO parent domain:

    >>> nemo.add_geoindex(grid="gridT")

    """
    # -- Set geographical indexes -- #
    _, dom_prefix, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
    lon_name = f"{dom_prefix}glam{grid_suffix}"
    lat_name = f"{dom_prefix}gphi{grid_suffix}"
    self_copy = self.copy()
    self_copy[grid] = (
        self_copy[grid]
        .dataset.assign_coords(
            {lat_name: self_copy[grid][lat_name], lon_name: self_copy[grid][lon_name]}
        )
        .set_xindex(
            (lat_name, lon_name),
            NDPointIndex,
            tree_adapter_cls=SklearnGeoBallTreeAdapter,
        )
    )

    return self_copy

binned_statistic

binned_statistic(grid, vars, values, keep_dims, bins, statistic, mask)

Calculate binned statistic of a variable defined on a NEMO model grid.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid where variables and values are stored (e.g., 'gridT').	required
vars	list[str]	Names of variable(s) to be grouped in discrete bins.	required
values	str	Name of the values with which to calculate binned statistic.	required
keep_dims	list[str] | None	Names of dimensions in values to keep as labels in binned statistic.	required
bins	list[list | ndarray]	Bin edges used to group each of the variables in vars.	required
statistic	str	Statistic to calculate (e.g., 'count', 'sum', 'nansum', 'mean', 'nanmean', 'max', 'nanmax', 'min', 'nanmin'). See flox.xarray.xarray_reduce for a complete list of aggregation statistics.	required
mask	DataArray | None	Boolean mask identifying NEMO model grid points to be included (1) or neglected (0) from calculation.	required

Returns:

Type	Description
DataArray	Values of the selected statistic in each bin.

Examples:

Compute the mean depth associated with each isopycnal in discrete potential density sigma0 coordinates:

>>> sigma0_bins = np.arange(22, 29.05, 0.05)

>>> nemo.binned_statistic(grid="gridT",
...                       vars=["sigma0"],
...                       values="deptht",
...                       keep_dims=["time_counter"],
...                       bins=[sigma0_bins],
...                       statistic="nanmean",
...                       )

See Also

masked_statistic

Source code in nemo_cookbook/nemodatatree.py

def binned_statistic(
    self,
    grid: str,
    vars: list[str],
    values: str,
    keep_dims: list[str] | None,
    bins: list[list | np.ndarray],
    statistic: str,
    mask: xr.DataArray | None,
) -> xr.DataArray:
    """
    Calculate binned statistic of a variable defined on a NEMO model grid.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid where variables and values are stored
        (e.g., 'gridT').
    vars : list[str]
        Names of variable(s) to be grouped in discrete bins.
    values : str
        Name of the values with which to calculate binned statistic.
    keep_dims : list[str] | None
        Names of dimensions in values to keep as labels in binned statistic.
    bins : list[list | np.ndarray]
        Bin edges used to group each of the variables in `vars`.
    statistic : str
        Statistic to calculate (e.g., 'count', 'sum', 'nansum', 'mean', 'nanmean',
        'max', 'nanmax', 'min', 'nanmin'). See flox.xarray.xarray_reduce for a
        complete list of aggregation statistics.
    mask : xr.DataArray | None
        Boolean mask identifying NEMO model grid points to be included (1)
        or neglected (0) from calculation.

    Returns
    -------
    xr.DataArray
        Values of the selected statistic in each bin.

    Examples
    --------
    Compute the mean depth associated with each isopycnal in discrete potential
    density `sigma0` coordinates:

    >>> sigma0_bins = np.arange(22, 29.05, 0.05)

    >>> nemo.binned_statistic(grid="gridT",
    ...                       vars=["sigma0"],
    ...                       values="deptht",
    ...                       keep_dims=["time_counter"],
    ...                       bins=[sigma0_bins],
    ...                       statistic="nanmean",
    ...                       )

    See Also
    --------
    masked_statistic
    """
    # -- Validate input -- #
    grid_keys = list(dict(self.subtree_with_keys).keys())
    if grid not in grid_keys:
        raise KeyError(
            f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
        )
    if any(var not in self[grid].data_vars for var in vars):
        raise KeyError(f"one or more variables {vars} not found in grid '{grid}'.")
    if values not in self[grid].data_vars:
        raise KeyError(f"values '{values}' not found in grid '{grid}'.")
    if keep_dims is not None:
        if any(dim not in self[grid][values].dims for dim in keep_dims):
            raise KeyError(
                f"one or more dimensions {keep_dims} not found in values '{values}'."
            )
    if not all(isinstance(bin, (list, np.ndarray)) for bin in bins):
        raise ValueError("bins must be a list of lists or numpy arrays.")
    if statistic not in [
        "all",
        "any",
        "count",
        "sum",
        "nansum",
        "mean",
        "nanmean",
        "max",
        "nanmax",
        "min",
        "nanmin",
        "argmax",
        "nanargmax",
        "argmin",
        "nanargmin",
        "quantile",
        "nanquantile",
        "median",
        "nanmedian",
        "mode",
        "nanmode",
        "first",
        "nanfirst",
        "last",
        "nanlast",
    ]:
        raise ValueError(f"statistic '{statistic}' is not supported.")
    if mask is not None:
        if not isinstance(mask, xr.DataArray):
            raise ValueError("mask must be an xarray.DataArray.")
        if mask.dtype != bool:
            raise TypeError("mask dtype must be boolean.")
        if any(dim not in self[grid].dims for dim in mask.dims):
            raise ValueError(
                f"mask must have dimensions subset from {self[grid].dims}."
            )

    # -- Calculate binned statistics -- #
    values_data = self[f"{grid}/{values}"].masked.data
    var_data = [self[f"{grid}/{var}"].masked.data for var in vars]

    result = compute_binned_statistic(
        vars=var_data,
        values=values_data,
        keep_dims=keep_dims,
        bins=bins,
        statistic=statistic,
        mask=mask,
    )

    return result

cell_area

cell_area(grid, dim)

Calculate grid cell areas orthogonal to a given dimension of a NEMO model grid.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid from which to calculate grid cell areas (e.g., 'gridT').	required
dim	str	Dimension orthogonal to grid cell area to calculate (e.g., 'k' returns e1 * e2).	required

Returns:

Type	Description
DataArray	Grid cell areas (m^2) for the specified NEMO model grid.

Examples:

Compute the horizontal area of each grid cell centered on a V-grid point in the NEMO parent domain:

>>> nemo.cell_area(grid="gridT", dim="k")

Note, dim represents the dimension orthogonal to the grid cell area to be computed.

See Also

cell_volume

Source code in nemo_cookbook/nemodatatree.py

def cell_area(
    self,
    grid: str,
    dim: str,
) -> xr.DataArray:
    """
    Calculate grid cell areas orthogonal to a given dimension of a NEMO model grid.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid from which to calculate grid cell areas
        (e.g., 'gridT').
    dim : str
        Dimension orthogonal to grid cell area to
        calculate (e.g., 'k' returns e1 * e2).

    Returns
    -------
    xr.DataArray
        Grid cell areas (m^2) for the specified NEMO model grid.

    Examples
    --------
    Compute the horizontal area of each grid cell centered on a V-grid point
    in the NEMO parent domain:

    >>> nemo.cell_area(grid="gridT", dim="k")

    Note, `dim` represents the dimension orthogonal to the grid cell
    area to be computed.

    See Also
    --------
    cell_volume
    """
    grid_suffix = self._get_properties(grid=grid)

    if dim not in ["i", "j", "k"]:
        raise ValueError(f"dim {dim} must be one of ['i', 'j', 'k'].")

    match dim:
        case "i":
            cell_area = (
                self[f"{grid}/e3{grid_suffix}"].masked.data * self[f"{grid}/e2{grid_suffix}"].masked.data
            )
        case "j":
            cell_area = (
                self[f"{grid}/e3{grid_suffix}"].masked.data * self[f"{grid}/e1{grid_suffix}"].masked.data
            )
        case "k":
            cell_area = (
                self[f"{grid}/e1{grid_suffix}"].masked.data * self[f"{grid}/e2{grid_suffix}"].masked.data
            )
    cell_area.name = "areacello"

    return cell_area

cell_volume

cell_volume(grid)

Calculate grid cell volumes for a given NEMO model grid.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid to calculate grid cell volumes (e.g., 'gridT').	required

Returns:

Type	Description
DataArray	Grid cell volumes for the specified NEMO model grid.

Examples:

Compute the volume of each grid cell centered on a V-grid point in the NEMO parent domain:

>>> nemo.cell_volumes(grid="gridV")

See Also

cell_area

Source code in nemo_cookbook/nemodatatree.py

def cell_volume(self, grid: str) -> xr.DataArray:
    """
    Calculate grid cell volumes for a given NEMO model grid.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid to calculate grid cell volumes
        (e.g., 'gridT').

    Returns
    -------
    xr.DataArray
        Grid cell volumes for the specified NEMO model grid.

    Examples
    --------
    Compute the volume of each grid cell centered on a V-grid point
    in the NEMO parent domain:

    >>> nemo.cell_volumes(grid="gridV")

    See Also
    --------
    cell_area
    """
    grid_suffix = self._get_properties(grid=grid)

    cell_volume = (
        self[f"{grid}/e3{grid_suffix}"].masked.data
        * self[f"{grid}/e1{grid_suffix}"].masked.data
        * self[f"{grid}/e2{grid_suffix}"].masked.data
    )
    cell_volume.name = "volcello"

    return cell_volume

clip_domain

clip_domain(dom, bbox)

Clip a NEMO model domain to specified longitude and latitude range.

Parameters:

Name	Type	Description	Default
dom	str	Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.). Default is '.' for the parent domain.	required
bbox	tuple	Bounding box to clip to (lon_min, lon_max, lat_min, lat_max).	required

Returns:

Type	Description
NEMODataTree	NEMO DataTree with specified model domain clipped to bounding box.

Examples:

Clip all model grids in a NEMO parent domain in the bounding box (-80°E, 0°E, 40°N, 80°N):

>>> bbox = (-80, 0, 40, 80)

>>> nemo.clip_domain(dom=".", bbox=bbox)

See Also

clip_grid

Source code in nemo_cookbook/nemodatatree.py

def clip_domain(
    self,
    dom: str,
    bbox: tuple,
) -> Self:
    """
    Clip a NEMO model domain to specified longitude and latitude range.

    Parameters
    ----------
    dom : str
        Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
        Default is '.' for the parent domain.
    bbox : tuple
        Bounding box to clip to (lon_min, lon_max, lat_min, lat_max).

    Returns
    -------
    NEMODataTree
        NEMO DataTree with specified model domain clipped to bounding box.

    Examples
    --------
    Clip all model grids in a NEMO parent domain in the bounding box
    (-80°E, 0°E, 40°N, 80°N):

    >>> bbox = (-80, 0, 40, 80)

    >>> nemo.clip_domain(dom=".", bbox=bbox)

    See Also
    --------
    clip_grid
    """
    # -- Validate input -- #
    if not isinstance(dom, str):
        raise ValueError(
            "dom must be a string specifying the prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
        )
    if not isinstance(bbox, tuple) or len(bbox) != 4:
        raise ValueError(
            "bounding box must be a tuple: (lon_min, lon_max, lat_min, lat_max)."
        )

    # -- Get NEMO model grid properties -- #
    grid_paths = self._get_grid_paths(dom=dom)
    ijk_names = self._get_ijk_names(dom=dom)
    i_name, j_name = ijk_names["i"], ijk_names["j"]

    # -- Clip grids to given bounding box -- #
    if not grid_paths:
        raise ValueError(f"NEMO model domain '{dom}' not found in the DataTree.")
    else:
        for grid in grid_paths.values():
            # Identify grid type:
            grid_suffix = self._get_properties(grid=grid)

            if grid_suffix == "t":
                # Clip shallow copy of NEMODataTree T-grid using lon/lat bbox:
                self_copy = self.copy().clip_grid(grid=grid, bbox=bbox)
                # Store (i, j) coords of bbox on T-grid:
                i_bbox = self_copy[grid][i_name]
                j_bbox = self_copy[grid][j_name]

            else:
                # Clip adjacent horizontal grid using (i, j) coords of clipped T-grid:
                match grid_suffix:
                    case "u":
                        grid_clipped = self[grid].dataset.sel(i=i_bbox + 0.5, j=j_bbox)
                    case "v":
                        grid_clipped = self[grid].dataset.sel(i=i_bbox, j=j_bbox + 0.5)
                    case "w":
                        grid_clipped = self[grid].dataset.sel(i=i_bbox, j=j_bbox)
                    case "f":
                        grid_clipped = self[grid].dataset.sel(i=i_bbox + 0.5, j=j_bbox + 0.5)

                if bbox != (-180, 180, -90, 90):
                    # Update zonal periodicity of NEMO model grid:
                    grid_clipped = grid_clipped.assign_attrs({"iperio": False})
                # Update shallow copy of NEMODataTree:
                self_copy[grid].dataset = grid_clipped

    return self_copy

clip_grid

clip_grid(grid, bbox)

Clip a NEMO model grid to specified longitude and latitude range.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid to clip (e.g., 'gridT').	required
bbox	tuple	Bounding box to clip to (lon_min, lon_max, lat_min, lat_max).	required

Returns:

Type	Description
NEMODataTree	NEMO DataTree with specified model grid clipped to bounding box.

Examples:

Clip T-grid in a NEMO parent domain in the bounding box (-80°E, 0°E, 40°N, 80°N):

>>> bbox = (-80, 0, 40, 80)

>>> nemo.clip_grid(grid="gridT", bbox=bbox)

See Also

clip_domain

Source code in nemo_cookbook/nemodatatree.py

def clip_grid(
    self,
    grid: str,
    bbox: tuple,
) -> Self:
    """
    Clip a NEMO model grid to specified longitude and latitude range.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid to clip (e.g., 'gridT').
    bbox : tuple
        Bounding box to clip to (lon_min, lon_max, lat_min, lat_max).

    Returns
    -------
    NEMODataTree
        NEMO DataTree with specified model grid clipped to bounding box.

    Examples
    --------
    Clip T-grid in a NEMO parent domain in the bounding box (-80°E, 0°E, 40°N, 80°N):

    >>> bbox = (-80, 0, 40, 80)

    >>> nemo.clip_grid(grid="gridT", bbox=bbox)

    See Also
    --------
    clip_domain
    """
    # -- Validate input -- #
    if not isinstance(bbox, tuple) or len(bbox) != 4:
        raise ValueError(
            "bounding box must be a tuple (lon_min, lon_max, lat_min, lat_max)."
        )

    # -- Get NEMO model grid properties -- #
    _, dom_prefix, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)

    # -- Clip the grid to given bounding box -- #
    # Indexing with a mask requires loading coords into memory:
    glam = self[grid][f"{dom_prefix}glam{grid_suffix}"].load()
    gphi = self[grid][f"{dom_prefix}gphi{grid_suffix}"].load()

    grid_clipped = self[grid].dataset.where(
        (glam >= bbox[0])
        & (glam <= bbox[1])
        & (gphi >= bbox[2])
        & (gphi <= bbox[3]),
        drop=True,
    )

    d_dtypes = {var: self[grid][var].dtype for var in self[grid].dataset.data_vars}
    for var, dtype in d_dtypes.items():
        if dtype in [np.int32, np.int64, bool]:
            grid_clipped[var] = grid_clipped[var].fillna(0).astype(dtype)

    if bbox != (-180, 180, -90, 90):
        # Update zonal periodicity of grid node:
        grid_clipped = grid_clipped.assign_attrs({"iperio": False})

    # Update shallow copy of NEMODataTree:
    self_copy = self.copy()
    self_copy[grid].dataset = grid_clipped

    return self_copy

curl

curl(vars, dom='.')

Calculate the vertical (k) curl component of a vector field on a NEMO model grid.

Parameters:

Name	Type	Description	Default
vars	list[str]	Name of the vector variables, structured as: ['u', 'v'], where 'u' and 'v' are the i and j components of the vector field, respectively.	required
dom	str	Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.). Default is '.' for the parent domain.	'.'

Returns:

Type	Description
DataArray	Vertical curl component of vector field defined on a NEMO model grid.

Examples:

Compute the vertical component of the curl of the seawater velocity field in the second NEMO nested child domain:

>>> nemo.curl(dom="2", vars=["uo", "vo"])

Note, vars expects a list of the i and j components of the vector field, respectively.

See Also

divergence

Source code in nemo_cookbook/nemodatatree.py

@deprecated(version_since="2026.03.b1",
            version_removed="2026.05",
            alternative="NEMOVectorField.curl from v2026.04 onwards"
            )
def curl(
    self,
    vars: list[str],
    dom: str = ".",
) -> xr.DataArray:
    """
    Calculate the vertical (k) curl component of a vector field on a NEMO model grid.

    Parameters
    ----------
    vars : list[str]
        Name of the vector variables, structured as: ['u', 'v'], where 'u' and 'v' are
        the i and j components of the vector field, respectively.
    dom : str, optional
        Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
        Default is '.' for the parent domain.

    Returns
    -------
    xr.DataArray
        Vertical curl component of vector field defined on a NEMO model grid.

    Examples
    --------
    Compute the vertical component of the curl of the seawater velocity field in
    the second NEMO nested child domain:

    >>> nemo.curl(dom="2", vars=["uo", "vo"])

    Note, `vars` expects a list of the `i` and `j` components of the vector field,
    respectively.

    See Also
    --------
    divergence
    """
    # -- Validate input -- #
    if not isinstance(vars, list) or len(vars) != 2:
        raise ValueError(
            "vars must be a list of two elements structured as ['u', 'v']."
        )
    if not isinstance(dom, str):
        raise ValueError(
            "dom must be a string specifying the prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
        )

    # -- Get NEMO model grid properties -- #
    dom_prefix, _ = self._get_properties(dom=dom)
    grid_paths = self._get_grid_paths(dom=dom)
    gridU, gridV, gridF = (
        grid_paths["gridU"],
        grid_paths["gridV"],
        grid_paths["gridF"],
    )
    ijk_names = self._get_ijk_names(dom=dom)
    i_name, j_name = ijk_names["i"], ijk_names["j"]

    # -- Define i,j vector components -- #
    var_i, var_j = vars[0], vars[1]
    if var_i not in self[gridU].data_vars:
        raise KeyError(f"variable '{var_i}' not found in grid '{gridU}'.")
    if var_j not in self[gridV].data_vars:
        raise KeyError(f"variable '{var_j}' not found in grid '{gridV}'.")

    da_i = self[f"{gridU}/{var_i}"].masked.data
    da_j = self[f"{gridV}/{var_j}"].masked.data

    # -- Collect mask -- #
    if (f"{dom_prefix}depthu" in da_i.coords) and (
        f"{dom_prefix}depthv" in da_j.coords
    ):
        # 3-dimensional fmask
        fmask = self[gridF]["fmask"]
    else:
        # 2-dimensional fmask (unique points):
        fmask = self[gridF]["fmaskutil"]

    # -- Neglecting the final F-grid points along i, j dimensions -- #
    e1f = self[f"{gridF}/e1f"].masked.data.isel(
        {i_name: slice(None, -1), j_name: slice(None, -1)}
    )
    e2f = self[f"{gridF}/e2f"].masked.data.isel(
        {i_name: slice(None, -1), j_name: slice(None, -1)}
    )

    e1u = self[f"{gridU}/e1u"].masked.data
    e2v = self[f"{gridV}/e2v"].masked.data
    # -- Calculate vertical curl component on F-points -- #
    dvar_i = (e2v * da_j).diff(dim=i_name, label="lower")
    dvar_i.coords[i_name] = dvar_i.coords[i_name] + 0.5

    dvar_j = (e1u * da_i).diff(dim=j_name, label="lower")
    dvar_j.coords[j_name] = dvar_j.coords[j_name] + 0.5

    curl = (1 / (e1f * e2f)) * (dvar_i - dvar_j).where(fmask)

    # -- Update DataArray properties -- #
    curl.name = f"curl({var_i}, {var_j})"
    curl = curl.drop_vars(
        [
            f"{dom_prefix}glamu",
            f"{dom_prefix}gphiu",
            f"{dom_prefix}glamv",
            f"{dom_prefix}gphiv",
        ]
    )

    return curl

depth_integral

depth_integral(grid, var, limits)

Integrate a variable in depth coordinates between two limits.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid where variable is stored (e.g., 'gridT').	required
var	str	Name of the variable to vertically integrate.	required
limits	tuple[int | float]	Limits of depth integration given as a tuple of the form (depth_lower, depth_upper) where depth_lower and depth_upper are the lower and upper limits of vertical integration, respectively.	required

Returns:

Type	Description
DataArray	Vertical integral of chosen variable between depth surfaces (depth_lower, depth_upper).

Examples:

Vertically integrate the conservative temperature variable thetao_con defined in a NEMO model parent domain from the sea surface to 100 m depth:

>>> nemo.depth_integral(grid='gridT',
...                     var='thetao_con',
...                     limits=(0, 100)
...                              )

See Also

integral

Source code in nemo_cookbook/nemodatatree.py

@deprecated(version_since="2026.03.b1",
            version_removed="2026.05",
            alternative="NEMODataArray.depth_integral"
            )
def depth_integral(
    self, grid: str, var: str, limits: tuple[int | float]
) -> xr.Dataset:
    """
    Integrate a variable in depth coordinates between two limits.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid where variable is stored (e.g., 'gridT').
    var : str
        Name of the variable to vertically integrate.
    limits : tuple[int | float]
        Limits of depth integration given as a tuple of the form
        (depth_lower, depth_upper) where depth_lower and depth_upper are
        the lower and upper limits of vertical integration, respectively.

    Returns
    -------
    xr.DataArray
        Vertical integral of chosen variable between depth surfaces (depth_lower, depth_upper).

    Examples
    --------
    Vertically integrate the conservative temperature variable `thetao_con` defined in a
    NEMO model parent domain from the sea surface to 100 m depth:

    >>> nemo.depth_integral(grid='gridT',
    ...                     var='thetao_con',
    ...                     limits=(0, 100)
    ...                              )

    See Also
    --------
    integral
    """
    # -- Validate input -- #
    grid_keys = list(dict(self.subtree_with_keys).keys())
    if grid not in grid_keys:
        raise KeyError(
            f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
        )
    if var not in self[grid].data_vars:
        raise KeyError(f"Variable '{var}' not found in grid '{grid}'.")
    if (not isinstance(limits, tuple)) | (len(limits) != 2):
        raise TypeError(
            "depth limits of integration should be given by a tuple of the form (depth_lower, depth_upper)"
        )
    if (limits[0] < 0) | (limits[1] < 0):
        raise ValueError("depth limits of integration must be non-negative.")
    if limits[0] >= limits[1]:
        raise ValueError(
            "lower depth limit must be less than upper depth limit."
        )

    # -- Get NEMO model grid properties -- #
    dom, _, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
    ijk_names = self._get_ijk_names(dom=dom)
    i_name, j_name, k_name = ijk_names["i"], ijk_names["j"], ijk_names["k"]

    # -- Define input variables -- #
    var_in = self[f"{grid}/{var}"].masked.data
    e3_in = self[f"{grid}/e3{grid_suffix}"].masked.data

    # -- Vertically integrate w.r.t depth -- #
    result = xr.apply_ufunc(
        compute_depth_integral,
        e3_in,
        var_in,
        np.array([limits[1]]),
        np.array([limits[0]]),
        input_core_dims=[[k_name], [k_name], [None], [None]],
        output_core_dims=[["k_new"]],
        dask="parallelized",
        output_dtypes=[var_in.dtype],
        dask_gufunc_kwargs={"output_sizes": {"k_new": 1}},
    )

    # -- Create variable integral DataArray -- #
    t_name = self[f"{grid}/{var}"].t_name
    if t_name is not None:
        result = result.transpose(t_name, "k_new", j_name, i_name).squeeze()
    else:
        result = result.transpose("k_new", j_name, i_name).squeeze()
    result.name = f"integral_z({var})"

    return result

divergence

divergence(vars, dom='.')

Calculate the horizontal divergence of a vector field defined on a NEMO model grid.

Parameters:

Name	Type	Description	Default
vars	list[str]	Name of vector variables, structured as: ['u', 'v'], where 'u' and 'v' are the i and j components of the vector field, respectively.	required
dom	str	Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.). Default is '.' for the parent domain.	'.'

Returns:

Type	Description
DataArray	Horizontal divergence of vector field defined on a NEMO model grid.

Examples:

Compute the horizontal divergence of the seawater velocity field in the NEMO parent domain:

>>> nemo.divergence(dom=".", vars=["uo", "vo"])

Note, vars expects a list of the i and j components of the vector field, respectively.

See Also

divergence

Source code in nemo_cookbook/nemodatatree.py

@deprecated(version_since="2026.03.b1",
            version_removed="2026.05",
            alternative="NEMOVectorField.divergence from v2026.04 onwards"
            )
def divergence(
    self,
    vars: list[str],
    dom: str = ".",
) -> xr.DataArray:
    """
    Calculate the horizontal divergence of a vector field defined on a NEMO model grid.

    Parameters
    ----------
    vars : list[str]
        Name of vector variables, structured as: ['u', 'v'], where
        'u' and 'v' are the i and j components of the vector field,
        respectively.
    dom : str, optional
        Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
        Default is '.' for the parent domain.

    Returns
    -------
    xr.DataArray
        Horizontal divergence of vector field defined on a NEMO model grid.

    Examples
    --------
    Compute the horizontal divergence of the seawater velocity field in the
    NEMO parent domain:

    >>> nemo.divergence(dom=".", vars=["uo", "vo"])

    Note, `vars` expects a list of the `i` and `j` components of the vector
    field, respectively.

    See Also
    --------
    divergence
    """
    # -- Validate input -- #
    if not isinstance(vars, list) or len(vars) != 2:
        raise ValueError(
            "vars must be a list of two elements structured as ['u', 'v']."
        )
    if not isinstance(dom, str):
        raise ValueError(
            "dom must be a string specifying the prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
        )

    # -- Get NEMO model grid properties -- #
    dom_prefix, _ = self._get_properties(dom=dom)
    grid_paths = self._get_grid_paths(dom=dom)
    gridT, gridU, gridV = (
        grid_paths["gridT"],
        grid_paths["gridU"],
        grid_paths["gridV"],
    )
    ijk_names = self._get_ijk_names(dom=dom)
    i_name, j_name = ijk_names["i"], ijk_names["j"]

    # -- Define i,j vector components -- #
    var_i, var_j = vars[0], vars[1]
    if var_i not in self[gridU].data_vars:
        raise KeyError(f"variable '{var_i}' not found in grid '{gridU}'.")
    if var_j not in self[gridV].data_vars:
        raise KeyError(f"variable '{var_j}' not found in grid '{gridV}'.")

    da_i = self[f"{gridU}/{var_i}"].masked.data
    da_j = self[f"{gridV}/{var_j}"].masked.data

    # -- Collect mask -- #
    if (f"{dom_prefix}depthu" in da_i.coords) and (
        f"{dom_prefix}depthv" in da_j.coords
    ):
        # 3-dimensional tmask:
        tmask = self[gridT]["tmask"]
    else:
        # 2-dimensional tmask (unique points):
        tmask = self[gridT]["tmaskutil"]

    # -- Neglecting the first T-grid points along i, j dimensions -- #
    e1t = self[f"{gridT}/e1t"].masked.data.isel({i_name: slice(1, None), j_name: slice(1, None)})
    e2t = self[f"{gridT}/e2t"].masked.data.isel({i_name: slice(1, None), j_name: slice(1, None)})
    e3t = self[f"{gridT}/e3t"].masked.data.isel({i_name: slice(1, None), j_name: slice(1, None)})

    e2u, e3u = self[f"{gridU}/e2u"].masked.data, self[f"{gridU}/e3u"].masked.data
    e1v, e3v = self[f"{gridV}/e1v"].masked.data, self[f"{gridV}/e3v"].masked.data  

    # -- Calculate divergence on T-points -- #
    dvar_i = (e2u * e3u * da_i).diff(dim=i_name, label="lower")
    dvar_i.coords[i_name] = dvar_i.coords[i_name] + 0.5

    dvar_j = (e1v * e3v * da_j).diff(dim=j_name, label="lower")
    dvar_j.coords[j_name] = dvar_j.coords[j_name] + 0.5

    divergence = (1 / (e1t * e2t * e3t)) * (dvar_i + dvar_j).where(tmask)

    # -- Update DataArray properties -- #
    divergence.name = f"div({var_i}, {var_j})"
    divergence = divergence.drop_vars(
        [
            f"{dom_prefix}glamu",
            f"{dom_prefix}gphiu",
            f"{dom_prefix}glamv",
            f"{dom_prefix}gphiv",
            f"{dom_prefix}depthu",
            f"{dom_prefix}depthv",
        ]
    )

    return divergence

extract_mask_boundary

extract_mask_boundary(mask, uv_vars=None, vars=None, dom='.')

Extract the boundary of a masked region defined on a NEMO model grid.

Parameters:

Name	Type	Description	Default
mask	DataArray	Boolean mask identifying NEMO model grid points which are inside the region of interest.	required
uv_vars	list	Names of velocity variables to extract along the boundary. Default is ['uo', 'vo'].	None
vars	list	Names of scalar variables to extract along the boundary.	None
dom	str	Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.). Default is '.' for the parent domain.	'.'

Returns:

Type	Description
Dataset	Dataset containing variables and NEMO model coordinates extracted along the boundary of the mask.

Examples:

Extract normal velocities and absolute salinity along the boundary of a masked region in the NEMO parent domain:

>>> nemo.extract_mask_boundary(mask=mask_osnap,
...                            uv_vars=["uo", "vo"],
...                            vars=["so_abs"],
...                            dom=".",
...                            )

See Also

extract_section

Source code in nemo_cookbook/nemodatatree.py

def extract_mask_boundary(
    self,
    mask: xr.DataArray,
    uv_vars: list | None = None,
    vars: list | None = None,
    dom: str = ".",
) -> xr.Dataset:
    """
    Extract the boundary of a masked region defined on a NEMO model grid.

    Parameters
    ----------
    mask : xr.DataArray
        Boolean mask identifying NEMO model grid points which
        are inside the region of interest.
    uv_vars : list, optional
        Names of velocity variables to extract along the boundary.
        Default is ['uo', 'vo'].
    vars : list, optional
        Names of scalar variables to extract along the boundary.
    dom : str, optional
        Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
        Default is '.' for the parent domain.

    Returns
    -------
    xr.Dataset
        Dataset containing variables and NEMO model coordinates
        extracted along the boundary of the mask.

    Examples
    --------
    Extract normal velocities and absolute salinity along the boundary of
    a masked region in the NEMO parent domain:

    >>> nemo.extract_mask_boundary(mask=mask_osnap,
    ...                            uv_vars=["uo", "vo"],
    ...                            vars=["so_abs"],
    ...                            dom=".",
    ...                            )

    See Also
    --------
    extract_section
    """
    # -- Validate Input -- #
    if not isinstance(mask, xr.DataArray):
        raise TypeError("mask must be an xarray DataArray")
    if not isinstance(dom, str):
        raise TypeError(
            "dom must be a string specifying prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
        )
    if uv_vars is None:
        uv_vars = ["uo", "vo"]
    else:
        if not isinstance(uv_vars, list) or len(uv_vars) != 2:
            raise TypeError(
                "uv_vars must be a list of velocity variables to extract (e.g., ['uo', 'vo'])."
            )

    # -- Get NEMO model grid properties -- #
    _, dom_suffix = self._get_properties(dom=dom)

    # -- Extract mask boundary -- #
    if f"i{dom_suffix}" not in mask.dims or f"j{dom_suffix}" not in mask.dims:
        raise ValueError(
            f"mask must have dimensions 'i{dom_suffix}' and 'j{dom_suffix}'"
        )
    i_bdy, j_bdy, flux_type, flux_dir = get_mask_boundary(mask)

    # -- Construct boundary dataset -- #
    # Neglecting final indices -> duplicate of the first indices:
    ds_bdy = create_boundary_dataset(
        nemo=self,
        dom=dom,
        i_bdy=i_bdy[:-1],
        j_bdy=j_bdy[:-1],
        flux_type=flux_type[:-1],
        flux_dir=flux_dir[:-1],
    )

    # -- Update boundary dataset with extracted section data -- #
    ds_bdy = update_boundary_dataset(
        ds_bdy=ds_bdy,
        nemo=self,
        dom=dom,
        sec_indexes=None,
        uv_vars=uv_vars,
        vars=vars,
    )

    return ds_bdy

extract_section

extract_section(lon_section, lat_section, uv_vars=None, vars=None, dom='.')

Extract hydrographic section from a NEMO model domain.

Parameters:

Name	Type	Description	Default
lon_section	ndarray	Longitudes defining the section polygon.	required
lat_section	ndarray	Latitudes defining the section polygon.	required
uv_vars	list	Names of velocity variables to extract along the boundary. Default is ['uo', 'vo'].	None
vars	list	Names of scalar variables to extract along the boundary.	None
dom	str	Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.). Default is '.' for the parent domain.	'.'

Returns:

Type	Description
Dataset	Dataset containing hydrographic section extracted from NEMO model grid.

Examples:

Extract normal velocities and potential density along the Overturning in the Subpolar North Atlantic (OSNAP) array defined by lon_osnap and lat_osnap coordinates in the NEMO parent domain:

>>> nemo.extract_section(lon_section=lon_osnap,
...                      lat_section=lat_osnap,
...                      uv_vars=["uo", "vo"],
...                      vars=["sigma0"],
...                      dom=".",
...                      )

See Also

extract_mask_boundary

Source code in nemo_cookbook/nemodatatree.py

def extract_section(
    self,
    lon_section: np.ndarray,
    lat_section: np.ndarray,
    uv_vars: list | None = None,
    vars: list | None = None,
    dom: str = ".",
) -> xr.Dataset:
    """
    Extract hydrographic section from a NEMO model domain.

    Parameters
    ----------
    lon_section : np.ndarray
        Longitudes defining the section polygon.
    lat_section : np.ndarray
        Latitudes defining the section polygon.
    uv_vars : list, optional
        Names of velocity variables to extract along the boundary.
        Default is ['uo', 'vo'].
    vars : list, optional
        Names of scalar variables to extract along the boundary.
    dom : str
        Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
        Default is '.' for the parent domain.

    Returns
    -------
    xr.Dataset
        Dataset containing hydrographic section extracted from NEMO model grid.

    Examples
    --------
    Extract normal velocities and potential density along the Overturning in the Subpolar
    North Atlantic (OSNAP) array defined by `lon_osnap` and `lat_osnap` coordinates in the
    NEMO parent domain:

    >>> nemo.extract_section(lon_section=lon_osnap,
    ...                      lat_section=lat_osnap,
    ...                      uv_vars=["uo", "vo"],
    ...                      vars=["sigma0"],
    ...                      dom=".",
    ...                      )

    See Also
    --------
    extract_mask_boundary
    """
    # -- Validate Input -- #
    if not isinstance(lon_section, np.ndarray):
        raise TypeError("lon_section must be a numpy array.")
    if not isinstance(lat_section, np.ndarray):
        raise TypeError("lat_section must be a numpy array.")
    if not isinstance(dom, str):
        raise TypeError(
            "dom must be a string specifying prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
        )
    if uv_vars is None:
        uv_vars = ["uo", "vo"]
    else:
        if not isinstance(uv_vars, list) or len(uv_vars) != 2:
            raise TypeError(
                "uv_vars must be a list of velocity variables to extract (e.g., ['uo', 'vo'])."
            )

    # -- Get NEMO model grid properties -- #
    grid_paths = self._get_grid_paths(dom=dom)

    # -- Define hydrographic section using polygon -- #
    lon_poly, lat_poly = create_section_polygon(
        lon_sec=lon_section,
        lat_sec=lat_section,
    )

    mask = self.mask_with_polygon(
        grid=grid_paths["gridT"], lon_poly=lon_poly, lat_poly=lat_poly
    )

    i_bdy, j_bdy, flux_type, flux_dir = get_mask_boundary(mask)

    # -- Create mask boundary dataset -- #
    ds_bdy = create_boundary_dataset(
        nemo=self,
        dom=dom,
        i_bdy=i_bdy,
        j_bdy=j_bdy,
        flux_type=flux_type,
        flux_dir=flux_dir,
    )

    # -- Get indexes of hydrographic section along mask boundary -- #
    dom_prefix, _ = self._get_properties(dom=dom)
    sec_indexes = get_section_indexes(
        lon_section=lon_section,
        lat_section=lat_section,
        gphib=ds_bdy[f"{dom_prefix}gphib"].values,
        glamb=ds_bdy[f"{dom_prefix}glamb"].values,
        bdy=ds_bdy["bdy"].values,
    )

    # -- Update boundary dataset with extracted section data -- #
    ds_bdy = update_boundary_dataset(
        ds_bdy=ds_bdy,
        nemo=self,
        dom=dom,
        sec_indexes=sec_indexes,
        uv_vars=uv_vars,
        vars=vars,
    )

    return ds_bdy

from_datasets classmethod

from_datasets(datasets, nests=None, name='NEMO model', iperio=False, nftype=None, read_mask=False, key_linssh=False, nbghost_child=4)

Create a NEMODataTree from a dictionary of xarray.Dataset objects created from NEMO model output files, organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

Parameters:

Name	Type	Description	Default
datasets	dict[str, dict[str, Dataset]]	Dictionary containing xarray.Datasets created from NEMO grid files, structured as: { 'parent': {'domain': ds_domain, 'gridT': ds_gridT, ... , 'icemod': ds_icemod.nc}, 'child': {'1': {'domain': ds_domain_1, 'gridT': d_gridT_1, ...}}, 'grandchild': {'2': {'domain': ds_domain_2, 'gridT': ds_gridT_2, ...}} }	required
nests	dict[str, dict[str, str]]	Dictionary describing the properties of nested domains, structured as: { "1": { "parent": "/", "rx": rx, "ry": ry, "imin": imin, "imax": imax, "jmin": jmin, "jmax": jmax, }, } where rx and ry are the horizontal refinement factors, and imin, imax, jmin, jmax define the indices of the child (grandchild) domain within the parent (child) domain.	None
name	str	Name of the NEMODataTree. Default is "NEMO model".	'NEMO model'
iperio	bool	Zonal periodicity of the parent domain.	False
nftype	str | None	Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point pivot. By default, no north fold lateral boundary condition is applied (None).	None
read_mask	bool	If True, read NEMO model land/sea mask from domain files. Default is False, meaning masks are computed from top_level and bottom_level domain variables.	False
key_linssh	bool	Linear free-surface approximation. If True, vertical coordinates are time-independent and given by (e3t_0, e3u_0, e3v_0, e3w_0) in domain_cfg. If False, vertical coordinates are time-dependent and must be specified in NEMO model grid datasets. Default is False.	False
nbghost_child	int = 4	Number of ghost cells to remove from the western/southern boundaries of the (grand)child domains. Default is 4.	4

Returns:

Type	Description
NEMODataTree	Hierarchical DataTree of NEMO model outputs.

Examples:

Create a zonally periodic NEMODataTree with north folding on T-points from a dictionary of xarray.Dataset objects:

>>> import xarray as xr
>>> from nemo_cookbook import NEMODataTree
>>> ds_domain = xr.open_zarr("https://some_remote_data/domain_cfg.zarr")
>>> ds_gridT = xr.open_zarr("https://some_remote_data/my_model_gridT.zarr")
>>> datasets = {"parent": {"domain": ds_domain, "gridT": ds_gridT}}
>>> nemo = NEMODataTree.from_datasets(datasets=datasets, name="My NEMO Model", iperio=True, nftype="T")

Create a regional NEMODataTree using a linear free-surface approximation from a dictionary of xarray.Dataset objects:

>>> nemo = NEMODataTree.from_datasets(datasets=datasets, name="My NEMO Model", iperio=False, nftype=None, key_linssh=True)

See Also

from_paths

Source code in nemo_cookbook/nemodatatree.py

@classmethod
def from_datasets(
    cls,
    datasets: dict[str, dict[str, xr.Dataset]],
    nests: dict[str, dict[str, str]] | None = None,
    name : str = "NEMO model",
    iperio: bool = False,
    nftype: str | None = None,
    read_mask: bool = False,
    key_linssh: bool = False,
    nbghost_child: int = 4,
) -> Self:
    """
    Create a NEMODataTree from a dictionary of `xarray.Dataset` objects created from NEMO model output files,
    organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

    Parameters
    ----------
    datasets : dict[str, dict[str, xr.Dataset]]
        Dictionary containing `xarray.Datasets` created from NEMO grid files, structured as:
        {
            'parent': {'domain': ds_domain, 'gridT': ds_gridT, ... , 'icemod': ds_icemod.nc},
            'child': {'1': {'domain': ds_domain_1, 'gridT': d_gridT_1, ...}},
            'grandchild': {'2': {'domain': ds_domain_2, 'gridT': ds_gridT_2, ...}}
        }

    nests : dict[str, dict[str, str]], optional
        Dictionary describing the properties of nested domains, structured as:
        {
            "1": {
                "parent": "/",
                "rx": rx,
                "ry": ry,
                "imin": imin,
                "imax": imax,
                "jmin": jmin,
                "jmax": jmax,
                },
        }
        where `rx` and `ry` are the horizontal refinement factors, and `imin`, `imax`, `jmin`, `jmax`
        define the indices of the child (grandchild) domain within the parent (child) domain.

    name: str, optional
        Name of the NEMODataTree. Default is "NEMO model".

    iperio: bool = False
        Zonal periodicity of the parent domain.

    nftype: str, optional
        Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point
        pivot. By default, no north fold lateral boundary condition is applied (None).

    read_mask: bool = False
        If True, read NEMO model land/sea mask from domain files. Default is False, meaning masks are computed from top_level and bottom_level domain variables.

    key_linssh: bool = False
        Linear free-surface approximation. If True, vertical coordinates are time-independent and given by (e3t_0, e3u_0, e3v_0, e3w_0) in domain_cfg.
        If False, vertical coordinates are time-dependent and must be specified in NEMO model grid datasets. Default is False.

    nbghost_child : int = 4
        Number of ghost cells to remove from the western/southern boundaries of the (grand)child domains. Default is 4.

    Returns
    -------
    NEMODataTree
        Hierarchical DataTree of NEMO model outputs.

    Examples
    --------
    Create a zonally periodic `NEMODataTree` with north folding on T-points from a dictionary of xarray.Dataset objects:

    >>> import xarray as xr
    >>> from nemo_cookbook import NEMODataTree
    >>> ds_domain = xr.open_zarr("https://some_remote_data/domain_cfg.zarr")
    >>> ds_gridT = xr.open_zarr("https://some_remote_data/my_model_gridT.zarr")
    >>> datasets = {"parent": {"domain": ds_domain, "gridT": ds_gridT}}
    >>> nemo = NEMODataTree.from_datasets(datasets=datasets, name="My NEMO Model", iperio=True, nftype="T")

    Create a regional `NEMODataTree` using a linear free-surface approximation from a dictionary of xarray.Dataset objects:

    >>> nemo = NEMODataTree.from_datasets(datasets=datasets, name="My NEMO Model", iperio=False, nftype=None, key_linssh=True)

    See Also
    --------
    from_paths
    """
    # -- Validate Inputs -- #
    if not isinstance(datasets, dict):
        raise TypeError("`datasets` must be a dictionary or nested dictionary.")
    if not isinstance(nests, (dict, type(None))):
        raise TypeError("`nests` must be a dictionary or None.")
    if not isinstance(name, str):
        raise TypeError("`name` must be a string.")
    if not isinstance(iperio, bool):
        raise TypeError("zonal periodicity (`iperio`) of parent domain must be a boolean.")
    if nftype is not None and nftype not in ("T", "F"):
        raise ValueError(
            "north fold type (`nftype`) of parent domain must be 'T' (T-pivot fold), 'F' (F-pivot fold), or None."
        )
    if not isinstance(read_mask, bool):
        raise TypeError("`read_mask` must be a boolean.")
    if not isinstance(key_linssh, bool):
        raise TypeError("linear free-surface approximation (`key_linssh`) must be a boolean.")
    if not isinstance(nbghost_child, int):
        raise TypeError(
            "number of ghost cells along the western/southern boundaries (`nbghost_child`) must be an integer."
        )

    # -- Define parent, child, grandchild dataset collections -- #
    d_child, d_grandchild = None, None
    if "parent" in datasets.keys() and isinstance(datasets["parent"], dict):
        for key in datasets.keys():
            if key not in ("parent", "child", "grandchild"):
                raise KeyError(f"Unexpected key '{key}' in `datasets` dictionary.")
            if key == "parent":
                d_parent = datasets["parent"]
            elif key == "child":
                if nests is None:
                    raise ValueError(
                        "`nests` dictionary must be provided when defining NEMO child domains."
                    )
                else:
                    d_child = datasets["child"]
            elif key == "grandchild":
                d_grandchild = datasets["grandchild"]
    else:
        raise ValueError(
            "Invalid `datasets` structure. Expected a nested dictionary defining NEMO 'parent', 'child' and 'grandchild' domains."
        )

    # -- Construct DataTree from parent / child / grandchild domains -- #
    d_tree = create_datatree_dict(
        d_parent=d_parent,
        d_child=d_child,
        d_grandchild=d_grandchild,
        nests=nests,
        iperio=iperio,
        nftype=nftype,
        read_mask=read_mask,
        key_linssh=key_linssh,
        nbghost_child=nbghost_child,
    )

    nemo = super().from_dict(d_tree)
    nemo.name = name

    return nemo

from_icechunk classmethod

from_icechunk(repo, name='NEMO model', iperio=False, nftype=None, open_kwargs=None, **session_kwargs)

Create a NEMODataTree from an Icechunk repository storing NEMO model outputs organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

Parameters:

Name	Type	Description	Default
repo	Repository	Icechunk repository containing NEMO model outputs organised into a hierarchy of domains.	required
name	str	Name of the NEMODataTree. Default is "NEMO model".	'NEMO model'
iperio	bool	Zonal periodicity of the parent domain. Default is False.	False
nftype	str | None	Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point	None
open_kwargs	dict[str, Any]	Additional keyword arguments to pass to xarray.open_datatree. Default is None.	None
**session_kwargs	dict[str, Any]	Additional keyword arguments to pass to Icechunk repo.readonly_session.	{}

Returns:

Type	Description
NEMODataTree	Hierarchical DataTree of NEMO model outputs.

Examples:

Create a zonally periodic NEMODataTree with north folding on F-points from the main branch of an Icechunk repository:

>>> from nemo_cookbook import NEMODataTree
>>> nemo = NEMODataTree.from_icechunk(repo=repo, branch="main", iperio=True, nftype="F")

See Also

from_zarr

Source code in nemo_cookbook/nemodatatree.py

@classmethod
def from_icechunk(
    cls,
    repo: icechunk.repository.Repository,
    name: str = "NEMO model",
    iperio: bool = False,
    nftype: str | None = None,
    open_kwargs: dict[str, any] | None = None,
    **session_kwargs: dict[str, any],
) -> Self:
    """
    Create a NEMODataTree from an Icechunk repository storing NEMO model outputs
    organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

    Parameters
    ----------
    repo : icechunk.repository.Repository
        Icechunk repository containing NEMO model outputs organised into a
        hierarchy of domains.

    name : str, optional
        Name of the NEMODataTree. Default is "NEMO model".

    iperio: bool = False
        Zonal periodicity of the parent domain. Default is False.

    nftype: str, optional
        Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point

    open_kwargs : dict[str, Any], optional
        Additional keyword arguments to pass to `xarray.open_datatree`. Default is None.

    **session_kwargs : dict[str, Any], optional
        Additional keyword arguments to pass to Icechunk `repo.readonly_session`.

    Returns
    -------
    NEMODataTree
        Hierarchical DataTree of NEMO model outputs.

    Examples
    --------
    Create a zonally periodic `NEMODataTree` with north folding on F-points from the main branch of an Icechunk repository:

    >>> from nemo_cookbook import NEMODataTree
    >>> nemo = NEMODataTree.from_icechunk(repo=repo, branch="main", iperio=True, nftype="F")

    See Also
    --------
    from_zarr
    """
    # -- Validate Inputs -- #
    if not hasattr(repo, "readonly_session"):
        raise TypeError("`repo` must implement readonly_session().")
    if not isinstance(name, str):
        raise TypeError("`name` must be a string.")
    if not isinstance(iperio, bool):
        raise TypeError("zonal periodicity (`iperio`) of parent domain must be a boolean.")
    if nftype is not None and nftype not in ("T", "F"):
        raise ValueError(
            "north fold type (`nftype`) of parent domain must be 'T' (T-pivot fold), 'F' (F-pivot fold), or None."
        )

    # -- Create NEMODataTree from Icechunk repository -- #:
    session = repo.readonly_session(**session_kwargs)
    datatree = xr.open_datatree(session.store, engine="zarr", **(open_kwargs or {}))
    nemo = super().from_dict(datatree.to_dict())

    # -- Update NEMODataTree properties -- #
    nemo["/"].attrs.update({"nftype": nftype, "iperio": iperio})
    nemo.name = name

    # -- Validate NEMO grid node Datasets -- #
    for key in [grid for grid in nemo.groups if grid.startswith("grid")]:
        validate_nemo_grid_node(key=key, value=nemo[key])

    return nemo

from_paths classmethod

from_paths(paths, nests=None, name='NEMO model', iperio=False, nftype=None, read_mask=False, key_linssh=False, nbghost_child=4, **open_kwargs)

Create a NEMODataTree from a dictionary of paths to NEMO model output files, organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

Parameters:

Name	Type	Description	Default
paths	dict[str, str]	Dictionary containing paths to NEMO grid files, structured as: { 'parent': {'domain': 'path/to/domain.nc', 'gridT': 'path/to/gridT.nc', , ... , 'icemod': 'path/to/icemod.nc', }, 'child': {'1': {'domain': 'path/to/child_domain.nc', 'gridT': 'path/to/child_gridT.nc', , ... , 'icemod': 'path/to/child_icemod.nc', }, }, 'grandchild': {'2': {'domain': 'path/to/grandchild_domain.nc', 'gridT': 'path/to/grandchild_gridT.nc', , ..., 'icemod': 'path/to/grandchild_icemod.nc', }, } }	required
nests	dict[str, str]	Dictionary describing the properties of nested domains, structured as: { "1": { "parent": "/", "rx": rx, "ry": ry, "imin": imin, "imax": imax, "jmin": jmin, "jmax": jmax, "iperio": iperio, }, } where rx and ry are the horizontal refinement factors, and imin, imax, jmin, jmax define the indices of the child (grandchild) domain within the parent (child) domain. Zonally periodic nested domains should be specified with iperio=True.	None
name	str	Name of the NEMODataTree. Default is "NEMO model".	'NEMO model'
iperio	bool	Zonal periodicity of the parent domain. Default is False.	False
nftype	str | None	Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point pivot. By default, no north fold lateral boundary condition is applied (None).	None
read_mask	bool	If True, read NEMO model land/sea mask from domain files. Default is False, meaning masks are computed from top_level and bottom_level domain variables.	False
key_linssh	bool	Linear free-surface approximation. If True, vertical coordinates are time-independent and given by (e3t_0, e3u_0, e3v_0, e3w_0) in domain_cfg. If False, vertical coordinates are time-dependent and must be specified in NEMO model grid datasets. Default is False.	False
nbghost_child	int = 4	Number of ghost cells to remove from the western/southern boundaries of the (grand)child domains. Default is 4.	4
**open_kwargs	dict	Additional keyword arguments to pass to xarray.open_dataset or xr.open_mfdataset when opening NEMO model output files.	{}

Returns:

Type	Description
NEMODataTree	A hierarchical DataTree storing NEMO model outputs.

Examples:

Create a zonally periodic NEMODataTree with north folding on T-points from a dictionary of paths to local netCDF files:

>>> from nemo_cookbook import NEMODataTree
>>> paths = {"parent": {
...          "domain": "/path/to/domain_cfg.nc",
...          "gridT": "path/to/*_gridT.nc",
...          "gridU": "path/to/*_gridV.nc",
...          "gridV": "path/to/*_gridV.nc",
...          "gridW": "path/to/*_gridW.nc",
...          "icemod": "path/to/*_icemod.nc",
...          }}
>>> nemo = NEMODataTree.from_paths(paths, name="My NEMO model", iperio=True, nftype="T")

Create a regional NEMODataTree using a linear free-surface approximation from a dictionary of paths to remote netCDF files:

>>> nemo = NEMODataTree.from_paths(paths, name="My NEMO model", iperio=False, nftype=None, key_linssh=True)

See Also

from_datasets

Source code in nemo_cookbook/nemodatatree.py

@classmethod
def from_paths(
    cls,
    paths: dict[str, str],
    nests: dict[str, str] | None = None,
    name : str = "NEMO model",
    iperio: bool = False,
    nftype: str | None = None,
    read_mask: bool = False,
    key_linssh: bool = False,
    nbghost_child: int = 4,
    **open_kwargs: dict[str, any],
) -> Self:
    """
    Create a NEMODataTree from a dictionary of paths to NEMO model output files,
    organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

    Parameters
    ----------
    paths : dict[str, str]
        Dictionary containing paths to NEMO grid files, structured as:
        {
            'parent': {'domain': 'path/to/domain.nc',
                       'gridT': 'path/to/gridT.nc',
                        , ... ,
                        'icemod': 'path/to/icemod.nc',
                        },
            'child': {'1': {'domain': 'path/to/child_domain.nc',
                            'gridT': 'path/to/child_gridT.nc',
                            , ... ,
                            'icemod': 'path/to/child_icemod.nc',
                            },
                      },
            'grandchild': {'2': {'domain': 'path/to/grandchild_domain.nc',
                                 'gridT': 'path/to/grandchild_gridT.nc',
                                 , ...,
                                 'icemod': 'path/to/grandchild_icemod.nc',
                                 },
                           }
        }

    nests : dict[str, str], optional
        Dictionary describing the properties of nested domains, structured as:
        {
            "1": {
                "parent": "/",
                "rx": rx,
                "ry": ry,
                "imin": imin,
                "imax": imax,
                "jmin": jmin,
                "jmax": jmax,
                "iperio": iperio,
                },
        }
        where `rx` and `ry` are the horizontal refinement factors, and `imin`, `imax`, `jmin`, `jmax`
        define the indices of the child (grandchild) domain within the parent (child) domain. Zonally
        periodic nested domains should be specified with `iperio=True`.

    name: str, optional
        Name of the NEMODataTree. Default is "NEMO model".

    iperio: bool = False
        Zonal periodicity of the parent domain. Default is False.

    nftype: str, optional
        Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point
        pivot. By default, no north fold lateral boundary condition is applied (None).

    read_mask: bool = False
        If True, read NEMO model land/sea mask from domain files. Default is False, meaning masks are computed from top_level and bottom_level domain variables.

    key_linssh: bool = False
        Linear free-surface approximation. If True, vertical coordinates are time-independent and given by (e3t_0, e3u_0, e3v_0, e3w_0) in domain_cfg.
        If False, vertical coordinates are time-dependent and must be specified in NEMO model grid datasets. Default is False.

    nbghost_child : int = 4
        Number of ghost cells to remove from the western/southern boundaries of the (grand)child domains. Default is 4.

    **open_kwargs : dict, optional
        Additional keyword arguments to pass to `xarray.open_dataset` or `xr.open_mfdataset` when opening NEMO model output files.
    Returns
    -------
    NEMODataTree
        A hierarchical DataTree storing NEMO model outputs.

    Examples
    --------
    Create a zonally periodic `NEMODataTree` with north folding on T-points from a dictionary of paths to local netCDF files:

    >>> from nemo_cookbook import NEMODataTree
    >>> paths = {"parent": {
    ...          "domain": "/path/to/domain_cfg.nc",
    ...          "gridT": "path/to/*_gridT.nc",
    ...          "gridU": "path/to/*_gridV.nc",
    ...          "gridV": "path/to/*_gridV.nc",
    ...          "gridW": "path/to/*_gridW.nc",
    ...          "icemod": "path/to/*_icemod.nc",
    ...          }}
    >>> nemo = NEMODataTree.from_paths(paths, name="My NEMO model", iperio=True, nftype="T")

    Create a regional `NEMODataTree` using a linear free-surface approximation from a dictionary of paths to remote netCDF files:

    >>> nemo = NEMODataTree.from_paths(paths, name="My NEMO model", iperio=False, nftype=None, key_linssh=True)

    See Also
    --------
    from_datasets
    """
    # -- Validate Inputs -- #
    if not isinstance(paths, dict):
        raise TypeError("`paths` must be a dictionary or nested dictionary.")
    if not isinstance(nests, (dict, type(None))):
        raise TypeError("`nests` must be a dictionary or None.")
    if not isinstance(name, str):
        raise TypeError("`name` must be a string.")
    if not isinstance(iperio, bool):
        raise TypeError("zonal periodicity (`iperio`) of parent domain must be a boolean.")
    if nftype is not None and nftype not in ("T", "F"):
        raise ValueError(
            "north fold type (`nftype`) of parent domain must be 'T' (T-pivot fold), 'F' (F-pivot fold), or None."
        )
    if not isinstance(read_mask, bool):
        raise TypeError("`read_mask` must be a boolean.")
    if not isinstance(key_linssh, bool):
        raise TypeError("linear free-surface approximation (`key_linssh`) must be a boolean.")
    if not isinstance(nbghost_child, int):
        raise TypeError(
            "number of ghost cells along the western/southern boundaries (`nbghost_child`) must be an integer."
        )
    if not isinstance(open_kwargs, dict):
        raise TypeError("`open_kwargs` must be a dictionary.")

    # -- Define parent, child, grandchild filepath collections -- #
    d_child, d_grandchild = None, None
    if "parent" in paths.keys() and isinstance(paths["parent"], dict):
        for key in paths.keys():
            if key not in ("parent", "child", "grandchild"):
                raise KeyError(f"Unexpected key '{key}' in `paths` dictionary.")
            if key == "parent":
                d_parent = paths["parent"]
            elif key == "child":
                if nests is None:
                    raise ValueError(
                        "`nests` dictionary must be provided when defining NEMO child domains."
                    )
                else:
                    d_child = paths["child"]
            elif key == "grandchild":
                d_grandchild = paths["grandchild"]
    else:
        raise ValueError(
            "Invalid `paths` structure. Expected a nested dictionary defining NEMO 'parent', 'child' and 'grandchild' domains."
        )

    # -- Construct DataTree from parent / child / grandchild domains -- #
    d_tree = create_datatree_dict(
        d_parent=d_parent,
        d_child=d_child,
        d_grandchild=d_grandchild,
        nests=nests,
        iperio=iperio,
        nftype=nftype,
        read_mask=read_mask,
        nbghost_child=nbghost_child,
        key_linssh=key_linssh,
        open_kwargs=dict(**open_kwargs),
    )

    nemo = super().from_dict(d_tree)
    nemo.name = name

    return nemo

from_zarr classmethod

from_zarr(store, name='NEMO model', iperio=False, nftype=None, **open_kwargs)

Create a NEMODataTree from Zarr store groups storing NEMO model outputs organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

Parameters:

Name	Type	Description	Default
store	str	Path to the Zarr store containing NEMO model outputs in hierarchical groups.	required
name	str	Name of the NEMODataTree. Default is "NEMO model".	'NEMO model'
iperio	bool	Zonal periodicity of the parent domain. Default is False.	False
nftype	str | None	Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point	None
**open_kwargs	dict[str, Any]	Additional keyword arguments to pass to xarray.open_datatree.	{}

Returns:

Type	Description
NEMODataTree	Hierarchical DataTree of NEMO model outputs.

Examples:

Create a zonally periodic NEMODataTree with north folding on T-points from a hierarchical Zarr store:

>>> from nemo_cookbook import NEMODataTree
>>> nemo = NEMODataTree.from_zarr(store="path/to/zarr/store", iperio=True, nftype="T")

See Also

from_icechunk

Source code in nemo_cookbook/nemodatatree.py

@classmethod
def from_zarr(
    cls,
    store: str,
    name: str = "NEMO model",
    iperio: bool = False,
    nftype: str | None = None,
    **open_kwargs: dict[str, any],
) -> Self:
    """
    Create a NEMODataTree from Zarr store groups storing NEMO model outputs
    organised into a hierarchy of domains (i.e., 'parent', 'child', 'grandchild').

    Parameters
    ----------
    store : str
        Path to the Zarr store containing NEMO model outputs in hierarchical groups.

    name : str, optional
        Name of the NEMODataTree. Default is "NEMO model".

    iperio: bool = False
        Zonal periodicity of the parent domain. Default is False.

    nftype: str, optional
        Type of north fold lateral boundary condition to apply. Options are 'T' for T-point pivot or 'F' for F-point

    **open_kwargs : dict[str, Any], optional
        Additional keyword arguments to pass to `xarray.open_datatree`.

    Returns
    -------
    NEMODataTree
        Hierarchical DataTree of NEMO model outputs.

    Examples
    --------
    Create a zonally periodic `NEMODataTree` with north folding on T-points from a hierarchical Zarr store:

    >>> from nemo_cookbook import NEMODataTree
    >>> nemo = NEMODataTree.from_zarr(store="path/to/zarr/store", iperio=True, nftype="T")

    See Also
    --------
    from_icechunk
    """
    # -- Validate Inputs -- #
    if not isinstance(store, str):
        raise TypeError("`store` must be a string.")
    if not isinstance(name, str):
        raise TypeError("`name` must be a string.")
    if not isinstance(iperio, bool):
        raise TypeError("zonal periodicity (`iperio`) of parent domain must be a boolean.")
    if nftype is not None and nftype not in ("T", "F"):
        raise ValueError(
            "north fold type (`nftype`) of parent domain must be 'T' (T-pivot fold), 'F' (F-pivot fold), or None."
        )

    # -- Create NEMODataTree from Zarr store -- #:
    datatree = xr.open_datatree(store, engine="zarr", **(open_kwargs or {}))
    nemo = super().from_dict(datatree.to_dict())

    # -- Update NEMODataTree properties -- #
    nemo["/"].attrs.update({"nftype": nftype, "iperio": iperio})
    nemo.name = name

    # -- Validate NEMO grid node Datasets -- #
    for key in [grid for grid in nemo.groups if grid.startswith("grid")]:
        validate_nemo_grid_node(key=key, value=nemo[key])

    return nemo

gradient

gradient(var, dim, dom='.')

Calculate the gradient of a scalar variable along one dimension (e.g., 'i', 'j', 'k') of a NEMO model grid.

Parameters:

Name	Type	Description	Default
var	str	Name of the scalar variable.	required
dim	str	Dimension along which to calculate gradient (e.g., 'i', 'j', 'k').	required
dom	str	Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.). Default is '.' for the parent domain.	'.'

Returns:

Type	Description
DataArray	Gradient of scalar variable defined on a NEMO model grid.

Examples:

Compute the 'meridional' gradient of sea surface temperature tos_con along the NEMO parent domain j dimension:

>>> nemo.gradient(dom='.', var="tos_con", dim="j")

Compute the vertical gradient of absolute salinity in the first NEMO nested child domain:

>>> nemo.gradient(dom="1", var="so_abs", dim="k")

See Also

integral

Source code in nemo_cookbook/nemodatatree.py

@deprecated(version_since="2026.03.b1",
            version_removed="2026.05",
            alternative="NEMODataArray.derivative or NEMOVectorField.gradient from v2026.04 onwards"
            )
def gradient(
    self,
    var: str,
    dim: str,
    dom: str = ".",
) -> xr.DataArray:
    """
    Calculate the gradient of a scalar variable along one dimension (e.g., 'i', 'j', 'k') of a NEMO model grid.

    Parameters
    ----------
    var : str
        Name of the scalar variable.
    dim : str
        Dimension along which to calculate gradient (e.g., 'i', 'j', 'k').
    dom : str, optional
        Prefix of NEMO domain in the DataTree (e.g., '1', '2', '3', etc.).
        Default is '.' for the parent domain.

    Returns
    -------
    xr.DataArray
        Gradient of scalar variable defined on a NEMO model grid.

    Examples
    --------
    Compute the 'meridional' gradient of sea surface temperature `tos_con`
    along the NEMO parent domain `j` dimension:

    >>> nemo.gradient(dom='.', var="tos_con", dim="j")

    Compute the vertical gradient of absolute salinity in the first NEMO
    nested child domain:

    >>> nemo.gradient(dom="1", var="so_abs", dim="k")

    See Also
    --------
    integral
    """
    # -- Validate input -- #
    if not isinstance(var, str):
        raise ValueError(
            "var must be a string specifying name of the scalar variable."
        )
    if not isinstance(dim, str):
        raise ValueError(
            "dim must be a string specifying dimension along which to calculate the gradient (e.g., 'i', 'j', 'k')."
        )
    if not isinstance(dom, str):
        raise ValueError(
            "dom must be a string specifying prefix of a NEMO domain (e.g., '.', '1', '2', etc.)."
        )

    # -- Get NEMO model grid properties -- #
    dom_prefix, dom_suffix = self._get_properties(dom=dom)
    grid_paths = self._get_grid_paths(dom=dom)
    gridT, gridU, gridV, gridW = (
        grid_paths["gridT"],
        grid_paths["gridU"],
        grid_paths["gridV"],
        grid_paths["gridW"],
    )

    if var not in self[gridT].data_vars:
        raise KeyError(f"variable '{var}' not found in grid '{gridT}'.")

    da = self[f"{gridT}/{var}"].masked.data
    dim_name = f"{dim}{dom_suffix}"
    if dim_name not in da.dims:
        raise KeyError(
            f"dimension '{dim_name}' not found in variable '{var}'. Dimensions available: {da.dims}."
        )

    match dim:
        case "i":
            if f"{dom_prefix}deptht" in da.coords:
                # 3-dimensional umask:
                umask = self[gridU]["umask"]
            else:
                # 2-dimensional umask:
                umask = self[gridU]["umaskutil"]

            # Zonally Periodic Domain:
            if self[gridT].attrs.get("iperio", False):
                da_end = da.isel({dim_name: 0})
                da_end[dim_name] = da[dim_name].max() + 1
                da = xr.concat([da, da_end], dim=dim_name)
                dvar = da.diff(dim=dim_name, label="lower")
            else:
                # Non-Periodic: pad with NaN values after differencing:
                dvar = da.diff(dim=dim_name, label="lower").pad({dim_name: (0, 1)})
            # Apply u-mask & transform coords -> calculate gradient:
            dvar.coords[dim_name] = dvar.coords[dim_name] + 0.5
            gradient = dvar.where(umask) / self[f"{gridU}/e1u"].masked.data

            # Remove redundant depth coordinates:
            if f"{dom_prefix}deptht" in gradient.coords:
                gradient = gradient.drop_vars(
                    [f"{dom_prefix}deptht"]
                ).assign_coords(
                    {f"{dom_prefix}depthu": self[gridU][f"{dom_prefix}depthu"]}
                )
        case "j":
            # 3-dimensional vmask:
            if f"{dom_prefix}deptht" in da.coords:
                vmask = self[gridV]["vmask"]
            else:
                # 2-dimensional vmask (unique points):
                vmask = self[gridV]["vmaskutil"]

            # Pad with zeros after differencing (zero gradient at jmaxdom):
            dvar = da.diff(dim=dim_name, label="lower").pad(
                {dim_name: (0, 1)}, constant_values=0
            )
            # Apply vmask & transform coords -> calculate gradient:
            dvar.coords[dim_name] = dvar.coords[dim_name] + 0.5
            gradient = dvar.where(vmask) / self[f"{gridV}/e2v"].masked.data

            if f"{dom_prefix}deptht" in gradient.coords:
                gradient = gradient.drop_vars(
                    [f"{dom_prefix}deptht"]
                ).assign_coords(
                    {f"{dom_prefix}depthv": self[gridV][f"{dom_prefix}depthv"]}
                )

        case "k":
            dvar = da.diff(dim=dim_name, label="lower")
            # Transform coords & apply w-mask -> calculate gradient:
            dvar.coords[dim_name] = dvar.coords[dim_name] + 0.5
            dvar = dvar.where(self[gridW]["wmask"].isel({dim_name: slice(1, None)}))
            try:
                gradient = -dvar / self[f"{gridW}/e3w"].masked.data.isel(
                    {dim_name: slice(1, None)}
                )
                gradient = gradient.drop_vars([f"{dom_prefix}deptht"])
            except KeyError as e:
                raise KeyError(
                    f"NEMO model grid: '{gridW}' does not contain vertical scale factor 'e3w' required to calculate gradients along the k-dimension."
                ) from e

    # Update DataArray properties:
    gradient.name = f"grad_{dim_name}({var})"
    gradient = gradient.drop_vars([f"{dom_prefix}glamt", f"{dom_prefix}gphit"])

    return gradient

integral

integral(grid, var, dims, cum_dims=None, dir=None, mask=None)

Integrate a variable along one or more dimensions of a NEMO model grid.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid where variable is stored (e.g., 'gridT').	required
var	str	Name of variable to integrate.	required
dims	list	Dimensions over which to integrate (e.g., ['i', 'k']).	required
cum_dims	list	Dimensions over which to cumulatively integrate (e.g., ['k']). Specified dimensions must also be included in dims.	None
dir	str	Direction of cumulative integration. Options are '+1' (along increasing cum_dims) or '-1' (along decreasing cum_dims).	None
mask	DataArray | None	Boolean mask identifying NEMO model grid points to be included (1) or neglected (0) from integration.	None

Returns:

Type	Description
DataArray	Variable integrated along specified dimensions of the NEMO model grid.

Examples:

Compute the integral of conservative temperature thetao_con along the vertical k dimension in the NEMO parent domain:

>>> nemo.integral(grid="gridT",
...               var="thetao_con",
...               dims=["k"]
...               )

Compute the vertical meridional overturning stream function from the meridional velocity vo (zonally integrated meridional velocity accumulated with increasing depth):

>>> nemo.integral(grid="gridV",
...               var="vo",
...               dims=["i", "k"],
...               cum_dims=["k"],
...               dir="+1",
...               )

See Also

gradient

Source code in nemo_cookbook/nemodatatree.py

@deprecated(version_since="2026.03.b1",
            version_removed="2026.05",
            alternative="NEMODataArray.integral"
            )
def integral(
    self,
    grid: str,
    var: str,
    dims: list,
    cum_dims: list | None = None,
    dir: str | None = None,
    mask: xr.DataArray | None = None,
) -> xr.DataArray:
    """
    Integrate a variable along one or more dimensions of a NEMO model grid.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid where variable is stored (e.g., 'gridT').
    var : str
        Name of variable to integrate.
    dims : list
        Dimensions over which to integrate (e.g., ['i', 'k']).
    cum_dims : list, optional
        Dimensions over which to cumulatively integrate (e.g., ['k']).
        Specified dimensions must also be included in `dims`.
    dir : str, optional
        Direction of cumulative integration. Options are '+1' (along
        increasing cum_dims) or '-1' (along decreasing cum_dims).
    mask: xr.DataArray, optional
        Boolean mask identifying NEMO model grid points to be included (1)
        or neglected (0) from integration.

    Returns
    -------
    xr.DataArray
        Variable integrated along specified dimensions of the NEMO model grid.


    Examples
    --------
    Compute the integral of conservative temperature `thetao_con` along the vertical
    `k` dimension in the NEMO parent domain:


    >>> nemo.integral(grid="gridT",
    ...               var="thetao_con",
    ...               dims=["k"]
    ...               )

    Compute the vertical meridional overturning stream function from the meridional
    velocity `vo` (zonally integrated meridional velocity accumulated with increasing
    depth):

    >>> nemo.integral(grid="gridV",
    ...               var="vo",
    ...               dims=["i", "k"],
    ...               cum_dims=["k"],
    ...               dir="+1",
    ...               )

    See Also
    --------
    gradient
    """
    # -- Validate input -- #
    grid_keys = list(dict(self.subtree_with_keys).keys())
    if grid not in grid_keys:
        raise KeyError(
            f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
        )
    if var not in self[grid].data_vars:
        raise KeyError(f"variable '{var}' not found in grid '{grid}'.")
    if cum_dims is not None:
        for dim in cum_dims:
            if dim not in dims:
                raise ValueError(
                    f"cumulative integration dimension '{dim}' not included in `dims`."
                )
        if dir not in ["+1", "-1"]:
            raise ValueError(
                f"invalid direction of cumulative integration '{dir}'. Expected '+1' or '-1'."
            )
    if mask is not None:
        if not isinstance(mask, xr.DataArray):
            raise ValueError("mask must be an xarray.DataArray.")
        if any(dim not in self[grid].dims for dim in mask.dims):
            raise ValueError(
                f"mask must have dimensions subset from {self[grid].dims}."
            )

    # -- Collect variable, weights & mask -- #
    da = (
        self[f"{grid}/{var}"].masked.data.where(mask)
        if mask is not None
        else self[f"{grid}/{var}"].masked.data
    )
    weights = self._get_weights(grid=grid, dims=dims)

    # -- Perform integration -- #
    if cum_dims is not None:
        sum_dims = [dim for dim in dims if dim not in cum_dims]
        if dir == "+1":
            # Cumulative integration along ordered dimension:
            result = (
                da.weighted(weights)
                .sum(dim=sum_dims, skipna=True)
                .cumsum(dim=cum_dims, skipna=True)
            )
        elif dir == "-1":
            # Cumulative integration along reversed dimension:
            result = (
                da.weighted(weights)
                .sum(dim=sum_dims, skipna=True)
                .reindex({dim: self[grid][dim][::-1] for dim in cum_dims})
                .cumsum(dim=cum_dims, skipna=True)
            )
    else:
        # Integration only:
        result = da.weighted(weights).sum(dim=dims, skipna=True)

    # -- Update DataArray properties -- #
    result.name = f"integral_{', '.join(dims)}({var})"

    return result

mask_with_polygon

mask_with_polygon(grid, lon_poly, lat_poly)

Create mask of NEMO model grid points contained within a polygon.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid where longitude and latitude coordinates are stored (e.g., 'gridT').	required
lon_poly	list | ndarray	Longitudes of closed polygon.	required
lat_poly	list | ndarray	Latitudes of closed polygon.	required

Returns:

Type	Description
DataArray	Boolean mask identifying NEMO model grid points which are inside the polygon.

Examples:

Create a regional boolean mask using the geographical coordinates of a closed polygon lon_poly and lat_poly in a NEMO parent domain:

>>> nemo.mask_with_polygon(grid="gridT",
...                        lon_poly=lon_poly,
...                        lat_poly=lat_poly,
...                        )

See Also

masked_statistic

Source code in nemo_cookbook/nemodatatree.py

def mask_with_polygon(
    self,
    grid: str,
    lon_poly: list | np.ndarray,
    lat_poly: list | np.ndarray,
):
    """
    Create mask of NEMO model grid points contained within a polygon.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid where longitude and latitude coordinates
        are stored (e.g., 'gridT').
    lon_poly : list | ndarray
        Longitudes of closed polygon.
    lat_poly : list | ndarray
        Latitudes of closed polygon.

    Returns
    -------
    xr.DataArray
        Boolean mask identifying NEMO model grid points which are inside
        the polygon.

    Examples
    --------
    Create a regional boolean mask using the geographical coordinates of a closed
    polygon `lon_poly` and `lat_poly` in a NEMO parent domain:


    >>> nemo.mask_with_polygon(grid="gridT",
    ...                        lon_poly=lon_poly,
    ...                        lat_poly=lat_poly,
    ...                        )

    See Also
    --------
    masked_statistic
    """
    # -- Validate input -- #
    if not isinstance(lon_poly, (np.ndarray, list)) or not isinstance(
        lat_poly, (np.ndarray, list)
    ):
        raise TypeError(
            "longitude & latitude coordinates of polygon must be numpy arrays or lists."
        )
    if (lon_poly[0] != lon_poly[-1]) or (lat_poly[0] != lat_poly[-1]):
        raise ValueError(
            "longitude & latitude coordinates must form a closed polygon."
        )

    # -- Get NEMO model grid properties -- #
    dom, dom_prefix, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
    ijk_names = self._get_ijk_names(grid=grid)
    i_name, j_name = ijk_names["i"], ijk_names["j"]

    if dom == ".":
        lon_name = f"glam{grid_suffix}"
        lat_name = f"gphi{grid_suffix}"
    else:
        lon_name = f"{dom_prefix}glam{grid_suffix}"
        lat_name = f"{dom_prefix}gphi{grid_suffix}"

    # -- Create mask using polygon coordinates -- #
    mask = create_polygon_mask(
        lon_grid=self[grid][lon_name],
        lat_grid=self[grid][lat_name],
        lon_poly=lon_poly,
        lat_poly=lat_poly,
        dims=(j_name, i_name),
    )

    return mask

masked_statistic

masked_statistic(grid, var, lon_poly, lat_poly, statistic, dims)

Masked statistic of a variable defined on a NEMO model grid.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid where variable is stored (e.g., 'gridT').	required
var	str	Name of the variable to compute statistic.	required
lon_poly	list | ndarray	Longitudes of closed polygon.	required
lat_poly	list | ndarray	Latitudes of closed polygon.	required
statistic	str	Name of the statistic to calculate (e.g., 'mean', 'weighted_mean' 'sum').	required
dims	list	Dimensions over which to apply statistic (e.g., ['i', 'j']).	required

Returns:

Type	Description
DataArray	Masked statistic of specified variable.

Examples:

Compute the grid cell area-weighted mean sea surface temperature tos_con for a region enclosed in a polygon defined by lon_poly and lat_poly in a NEMO nested child domain:

>>> nemo.masked_statistic(grid="gridT/1_gridT",
...                       var="tos_con",
...                       lon_poly=lon_poly,
...                       lat_poly=lat_poly,
...                       statistic="weighted_mean",
...                       dims=["i", "j"]
...                       )

See Also

binned_statistic

Source code in nemo_cookbook/nemodatatree.py

@deprecated(version_since="2026.03.b1",
            version_removed="2026.05",
            alternative="NEMODataArray.masked_statistic"
            )
def masked_statistic(
    self,
    grid: str,
    var: str,
    lon_poly: list | np.ndarray,
    lat_poly: list | np.ndarray,
    statistic: str,
    dims: list,
) -> xr.DataArray:
    """
    Masked statistic of a variable defined on a NEMO model grid.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid where variable is stored (e.g., 'gridT').
    var : str
        Name of the variable to compute statistic.
    lon_poly : list | np.ndarray
        Longitudes of closed polygon.
    lat_poly : list | np.ndarray
        Latitudes of closed polygon.
    statistic : str
        Name of the statistic to calculate (e.g., 'mean', 'weighted_mean' 'sum').
    dims : list
        Dimensions over which to apply statistic (e.g., ['i', 'j']).

    Returns
    -------
    xr.DataArray
        Masked statistic of specified variable.

    Examples
    --------
    Compute the grid cell area-weighted mean sea surface temperature `tos_con` for a
    region enclosed in a polygon defined by `lon_poly` and `lat_poly` in a NEMO nested
    child domain:

    >>> nemo.masked_statistic(grid="gridT/1_gridT",
    ...                       var="tos_con",
    ...                       lon_poly=lon_poly,
    ...                       lat_poly=lat_poly,
    ...                       statistic="weighted_mean",
    ...                       dims=["i", "j"]
    ...                       )

    See Also
    --------
    binned_statistic
    """
    # -- Validate input -- #
    grid_keys = list(dict(self.subtree_with_keys).keys())
    if grid not in grid_keys:
        raise KeyError(
            f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
        )
    if var not in self[grid].data_vars:
        raise KeyError(f"variable '{var}' not found in grid '{grid}'.")

    # -- Create polygon mask using coordinates -- #
    mask_poly = self.mask_with_polygon(
        lon_poly=lon_poly, lat_poly=lat_poly, grid=grid
    )

    # -- Get NEMO model grid properties -- #
    _, _, dom_suffix, _ = self._get_properties(grid=grid, infer_dom=True)

    # -- Apply masks & calculate statistic -- #
    da = self[f"{grid}/{var}"].masked.data.where(mask_poly)

    match statistic:
        case "mean":
            result = da.mean(dim=dims, skipna=True)

        case "weighted_mean":
            weight_dims = [dim.replace(dom_suffix, "") for dim in dims]
            weights = self._get_weights(grid=grid, dims=weight_dims)
            result = da.weighted(weights).mean(dim=dims, skipna=True)

        case "min":
            result = da.min(dim=dims, skipna=True)

        case "max":
            result = da.max(dim=dims, skipna=True)

        case "sum":
            result = da.sum(dim=dims, skipna=True)

        case _:
            raise ValueError(
                f"Unsupported statistic '{statistic}'. Supported statistics are: 'mean', 'weighted_mean', 'min', 'max', 'sum'."
            )

    return result

transform_to

transform_to(grid, var, to)

Transform variable defined on a NEMO model grid to a neighbouring horizontal grid using linear interpolation.

For flux variables defined at U- or V-points, the specified variable is first weighted by grid cell face areas prior to linear interpolation, and is then normalised by the target grid cell face areas following interpolation.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid where variable is stored (e.g., 'gridT').	required
var	str	Name of the variable to transform.	required
to	str	Suffix of the neighbouring horizontal NEMO model grid to transform variable to. Options are 'T', 'U', 'V', 'F'.	required

Returns:

Type	Description
DataArray	Values of variable linearly interpolated onto a neighbouring horizontal grid.

Examples:

Transform conservative temperature thetao_con defined on scalar T-points to neighbouring V-points in a NEMO model parent domain:

>>> nemo.transform_to(grid='gridT', var='thetao_con', to='V')

See Also

transform_vertical_grid

Source code in nemo_cookbook/nemodatatree.py

@deprecated(version_since="2026.03.b1",
            version_removed="2026.05",
            alternative="NEMODataArray.interp_to"
            )
def transform_to(
    self,
    grid: str,
    var: str,
    to: str,
) -> xr.DataArray:
    """
    Transform variable defined on a NEMO model grid to a neighbouring
    horizontal grid using linear interpolation.

    For flux variables defined at U- or V-points, the specified variable
    is first weighted by grid cell face areas prior to linear interpolation,
    and is then normalised by the target grid cell face areas following
    interpolation.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid where variable is stored (e.g., 'gridT').
    var : str
        Name of the variable to transform.
    to : str
        Suffix of the neighbouring horizontal NEMO model grid to
        transform variable to. Options are 'T', 'U', 'V', 'F'.

    Returns
    -------
    xr.DataArray
        Values of variable linearly interpolated onto a neighbouring
        horizontal grid.

    Examples
    --------
    Transform conservative temperature `thetao_con` defined on scalar T-points
    to neighbouring V-points in a NEMO model parent domain:

    >>> nemo.transform_to(grid='gridT', var='thetao_con', to='V')

    See Also
    --------
    transform_vertical_grid
    """
    # -- Validate input -- #
    grid_keys = list(dict(self.subtree_with_keys).keys())
    if grid not in grid_keys:
        raise KeyError(
            f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
        )
    if var not in self[grid].data_vars:
        raise KeyError(f"variable '{var}' not found in grid '{grid}'.")
    if not isinstance(to, str):
        raise TypeError(f"'to' must be a string, got {type(to)}.")
    if to not in ["T", "U", "V", "F"]:
        raise ValueError(f"'to' must be one of ['T', 'U', 'V', 'F'], got {to}.")

    # -- Get NEMO model grid properties -- #
    _, dom_prefix, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
    ijk_names = self._get_ijk_names(grid=grid)
    i_name, j_name, k_name = ijk_names["i"], ijk_names["j"], ijk_names["k"]
    iperio = self[grid].attrs.get("iperio", False)
    target_grid = f"{grid.replace(grid[-1], to)}"

    # -- Prepare variable for linear interpolation -- #
    if grid_suffix.upper() in ["U", "V"]:
        weight_dims = (
            [k_name, j_name] if grid_suffix.upper() == "U" else [k_name, i_name]
        )
        if f"{dom_prefix}depth{grid_suffix}" in self[grid][var].coords:
            # 3-D variables - weight by grid cell face area:
            weights = self._get_weights(grid=grid, dims=weight_dims, fillna=False)
            target_weights = self._get_weights(
                grid=target_grid, dims=weight_dims, fillna=False
            )
        else:
            # 2-D variables - weight by grid cell width:
            weights = self._get_weights(grid=grid, dims=weight_dims[1], fillna=False)
            target_weights = self._get_weights(
                grid=target_grid, dims=weight_dims[1], fillna=False
            )
        da = self[f"{grid}/{var}"].masked.data * weights
    else:
        # Scalar variables:
        da = self[f"{grid}/{var}"].masked.data

    # -- Linearly interpolate variable -- #
    # Define source grid mask:
    if f"{dom_prefix}depth{grid_suffix}" in da.coords:
        mask = self[grid][f"{grid_suffix}mask"]
    else:
        mask = self[grid][f"{grid_suffix}maskutil"]

    result = interpolate_grid(
        da=da,
        mask=mask,
        source_grid=grid[-1],
        target_grid=target_grid[-1],
        iperio=iperio,
        ijk_names=ijk_names,
    )

    # -- Update interpolated variable -- #
    # Retain input variable name:
    result.name = da.name
    # Reorder dimensions (time_counter, [k], j, i):
    new_dims = (result.dims[-1], *result.dims[:-1])
    result = result.transpose(*new_dims)

    # Update NEMO grid coords:
    result[i_name] = self[target_grid][i_name]
    result[j_name] = self[target_grid][j_name]
    if k_name in result.dims:
        result[k_name] = self[target_grid][k_name]

    # Drop NEMO source grid coords:
    drop_vars = [f"{dom_prefix}glam{grid_suffix}", f"{dom_prefix}gphi{grid_suffix}"]
    if f"{dom_prefix}depth{grid_suffix}" in da.coords:
        drop_vars.append(f"{dom_prefix}depth{grid_suffix}")
    result = result.drop_vars(drop_vars)

    # Normalise by target grid cell weights for flux variables:
    if grid_suffix.upper() in ["U", "V"]:
        result = result / target_weights

    # Apply target grid mask:
    if f"{dom_prefix}depth{grid_suffix}" in da.coords:
        target_mask = f"{to.lower()}mask"
    else:
        target_mask = f"{to.lower()}maskutil"
    result = result.where(self[target_grid][target_mask])

    return result

transform_vertical_grid

transform_vertical_grid(grid, var, e3_new)

Transform variable defined on a NEMO model grid to a new vertical grid using conservative interpolation.

Parameters:

Name	Type	Description	Default
grid	str	Path to NEMO model grid where variable is stored (e.g., 'gridT').	required
var	str	Name of the variable to transform.	required
e3_new	DataArray	Grid cell thicknesses of the new vertical grid. Must be a 1-dimensional xarray.DataArray with dimension 'k_new'.	required

Returns:

Type	Description
Dataset	Variable defined at the centre of each vertical grid cell on the new grid, and vertical grid cell thicknesses adjusted for model bathymetry.

Examples:

Transform the conservative temperature variable thetao_con defined in a NEMO model parent domain from it's native 75 unevenly-spaced z-levels to regularly spaced z-levels at 200 m intervals:

>>> e3t_target = xr.DataArray(np.repeat(200.0, 30), dims=['k_new'])

>>> nemo.transform_vertical_grid(grid='gridT',
...                              var='thetao_con',
...                              e3_new=e3t_target
...                              )

See Also

transform_to

Source code in nemo_cookbook/nemodatatree.py

@deprecated(version_since="2026.03.b1",
            version_removed="2026.05",
            alternative="NEMODataArray.transform_vertical_grid"
            )
def transform_vertical_grid(
    self, grid: str, var: str, e3_new: xr.DataArray
) -> xr.Dataset:
    """
    Transform variable defined on a NEMO model grid to a new vertical grid using conservative interpolation.

    Parameters
    ----------
    grid : str
        Path to NEMO model grid where variable is stored
        (e.g., 'gridT').
    var : str
        Name of the variable to transform.
    e3_new : xarray.DataArray
        Grid cell thicknesses of the new vertical grid.
        Must be a 1-dimensional xarray.DataArray with
        dimension 'k_new'.

    Returns
    -------
    xr.Dataset
        Variable defined at the centre of each vertical
        grid cell on the new grid, and vertical grid cell
        thicknesses adjusted for model bathymetry.

    Examples
    --------
    Transform the conservative temperature variable `thetao_con` defined in a
    NEMO model parent domain from it's native 75 unevenly-spaced z-levels to
    regularly spaced z-levels at 200 m intervals:

    >>> e3t_target = xr.DataArray(np.repeat(200.0, 30), dims=['k_new'])

    >>> nemo.transform_vertical_grid(grid='gridT',
    ...                              var='thetao_con',
    ...                              e3_new=e3t_target
    ...                              )

    See Also
    --------
    transform_to
    """
    # -- Validate input -- #
    grid_keys = list(dict(self.subtree_with_keys).keys())
    if grid not in grid_keys:
        raise KeyError(
            f"grid '{grid}' not found in available NEMODataTree grids {grid_keys}."
        )
    if var not in self[grid].data_vars:
        raise KeyError(f"Variable '{var}' not found in grid '{grid}'.")
    if e3_new.dims != ("k_new",) or (e3_new.ndim != 1):
        raise ValueError(
            "e3_new must be a 1-dimensional xarray.DataArray with dimension 'k_new'."
        )

    # -- Get NEMO model grid properties -- #
    dom, _, _, grid_suffix = self._get_properties(grid=grid, infer_dom=True)
    ijk_names = self._get_ijk_names(dom=dom)
    i_name, j_name, k_name = ijk_names["i"], ijk_names["j"], ijk_names["k"]
    t_name = [dim for dim in self[grid].dims if "time" in dim][0]

    # -- Define input variables -- #
    var_in = self[f"{grid}/{var}"].masked.data
    e3_in = self[f"{grid}/e3{grid_suffix}"].masked.data
    if e3_new.sum(dim="k_new") < self[grid][f"depth{grid_suffix}"].max(dim=k_name):
        raise ValueError(
            f"e3_new must sum to at least the maximum depth ({self[grid][f'depth{grid_suffix}'].max(dim=k_name).item()} m) of the original vertical grid."
        )

    # -- Transform variable to target vertical grid -- #
    var_out, e3_out = xr.apply_ufunc(
        transform_vertical_coords,
        e3_in,
        var_in,
        e3_new.astype(e3_in.dtype),
        input_core_dims=[[k_name], [k_name], ["k_new"]],
        output_core_dims=[["k_new"], ["k_new"]],
        dask="parallelized",
        output_dtypes=[var_in.dtype, e3_in.dtype],
        dask_gufunc_kwargs={"output_sizes": {"k_new": e3_new.sizes["k_new"]}},
    )

    # -- Construct transformed variable Dataset -- #
    var_out = var_out.transpose(t_name, "k_new", j_name, i_name)

    ds_out = xr.Dataset(
        data_vars={var: var_out, f"e3{grid_suffix}_new": e3_out},
        coords={
            f"depth{grid_suffix}_new": ("k_new", e3_new.cumsum(dim="k_new").data)
        },
    )

    return ds_out

nemo_cookbook.NEMODataArray

Grid-aware xarray.DataArray wrapper for NEMO ocean model variables.

This class interfaces with the NEMODataTree to provide discrete operators such as differences, derivatives, and integrals alongside transformation between NEMO model grids.

Methods:

Name	Description
__getattr__	Delegate attributes to xarray.DataArray and attempt to return
__getitem__	Delegate indexing to xarray.DataArray and attempt to return
__init__	Create a NEMODataArray from an xarray.DataArray and NEMODataTree.
__setitem__	Delegate item assignment to xarray.DataArray.
apply_mask	Apply NEMO parent grid land-sea mask or combined land-sea & custom mask
depth_integral	Integrate a variable in depth coordinates between two limits.
derivative	Calculate the derivative of a variable along a given dimension
diff	Calculate the discrete difference of a variable along a given dimension
integral	Integrate variable along one or more dimensions of a NEMO model grid.
interp_to	Linearly interpolate variable to a neighbouring horizontal grid.
masked_statistic	Compute masked statistic of a variable defined on a NEMO model grid.
sel_like	Return a new NEMODataArray whose data is given by matching the dimension
to_xesmf	Create an xESMF-compatible dataset from a variable stored in a NEMODataArray
transform_vertical_grid	Transform variable defined on a NEMO model grid to a new vertical grid using conservative interpolation.
weighted_mean	Calculate grid-aware weighted mean of a variable defined on a NEMO model grid.

Attributes:

Name	Type	Description
data	DataArray	Access underlying xarray.DataArray.
grid	str	Access path to parent NEMO model grid.
grid_type	str	Access type of parent NEMO model grid (e.g. 't', 'u', etc.).
iperio		Zonal periodicity of variable defined on NEMO model grid.
mask	DataArray	Access variable land-sea mask for the parent NEMO model grid.
masked		Apply land-sea mask to variable defined on NEMO model grid.
metrics	dict[str, 'NEMODataArray']	Access grid scale factors for the parent NEMO model grid.

Source code in nemo_cookbook/nemodataarray.py

class NEMODataArray:
    """
    Grid-aware xarray.DataArray wrapper for NEMO ocean model variables.

    This class interfaces with the NEMODataTree to provide discrete
    operators such as differences, derivatives, and integrals alongside
    transformation between NEMO model grids.
    """
    def __init__(self, da: xr.DataArray, tree: NEMODataTree, grid: str):
        """
        Create a NEMODataArray from an xarray.DataArray and NEMODataTree.

        Parameters
        ----------
        da : xr.DataArray
            Variable defined on a NEMO model grid.
        tree: NEMODataTree
            NEMODataTree to which the variable belongs.
        grid: str
            Path to NEMO model grid where variable is defined (e.g., 'gridT').

        Returns
        -------
        NEMODataArray
            Grid-aware xarray.DataArray storing variable defined on NEMO model grid.
        """
        # -- Validate Inputs -- #
        if not isinstance(da, xr.DataArray):
            raise TypeError("da must be specified as an xarray.DataArray.")
        else:
            self._da = da

        if not isinstance(tree, xr.DataTree):
            raise TypeError("tree must be specified as a NEMODataTree.")
        else:
            self._tree = tree

        if not isinstance(grid, str):
            raise TypeError("grid must be specified as a string.")
        else:
            self._grid = grid

        # -- Validate DataArray compatibility -- #
        validate_nemo_dataarray(da=self._da, grid=self._grid, tree=self._tree)

        # -- Assign NEMO domain number, grid path and grid type -- #
        dom, dom_prefix, dom_suffix, grid_suffix = self._tree._get_properties(grid=self._grid, infer_dom=True)
        self._dom = dom
        self._dom_prefix = dom_prefix
        self._dom_suffix = dom_suffix
        self._grid_suffix = grid_suffix

        # -- Assign NEMO dimension names -- #
        # Spatial dimensions:
        ijk_names = self._tree._get_ijk_names(grid=self._grid)
        self.i_name = ijk_names["i"]
        self.j_name = ijk_names["j"]
        self.k_name = ijk_names["k"]
        # Temporal dimension except for time-independent variables:
        # Use coords to handle single time slices of time-dependent variables.
        t_list = [dim for dim in self._da.coords if "time" in dim]
        self.t_name = t_list[0] if len(t_list) != 0 else None


    # -----------
    # Properties 
    # -----------
    @property
    def data(self) -> xr.DataArray:
        """
        Access underlying xarray.DataArray.
        """
        return self._da

    @property
    def grid(self) -> str:
        """
        Access path to parent NEMO model grid.
        """
        return self._grid

    @property
    def grid_type(self) -> str:
        """
        Access type of parent NEMO model grid (e.g. 't', 'u', etc.).
        """
        return self._grid_suffix

    @property
    def metrics(self) -> dict[str, "NEMODataArray"]:
        """
        Access grid scale factors for the parent NEMO model grid.
        """
        # 2-dimensional: Horizontal grid scale factors:
        d_metrics = {"e1": self._tree[f"{self._grid}/e1{self._grid_suffix}"].sel_like(self._da),
                     "e2": self._tree[f"{self._grid}/e2{self._grid_suffix}"].sel_like(self._da),
                     }

        if f"{self._dom_prefix}depth{self._grid_suffix}" in self._da.coords:
            # 3-dimensional: Vertical grid scale factor:
            d_metrics["e3"] = self._tree[f"{self._grid}/e3{self._grid_suffix}"].sel_like(self._da)

        return d_metrics

    @property
    def mask(self) -> xr.DataArray:
        """
        Access variable land-sea mask for the parent NEMO model grid.
        """
        if (self._da.ndim == 1) and (self.t_name in self._da.dims):
            raise ValueError("land-sea mask does not exist for variables without spatial dimensions.")
        else:
            if (f"{self._dom_prefix}depth{self._grid_suffix}" in self._da.coords) and (self.k_name in self._da.dims):
                # 3-dimensional land-sea mask:
                mask_name = f"{self._grid_suffix}mask"
            else:
                # 2-dimensional land-sea mask:
                mask_name = f"{self._grid_suffix}maskutil"

        # Select land-sea mask values according to variable dimension labels:
        result = self._tree[f"{self._grid}/{mask_name}"].sel_like(self._da).data

        return result

    @property
    def masked(self):
        """
        Apply land-sea mask to variable defined on NEMO model grid.
        """
        return self.apply_mask(mask=None, drop=False)

    @property
    def iperio(self):
        """
        Zonal periodicity of variable defined on NEMO model grid.
        """
        # Determine zonal periodicity of parent NEMO model grid:
        iperio = self._tree[self._grid].attrs.get("iperio", False)

        # Update zonal periodicity if variable is defined on a subset
        # of i-dimension labels of the parent NEMO model grid:
        if iperio:
            if self._da.coords[self.i_name].size != self._tree[self._grid].coords[self.i_name].size:
                iperio = False

        return iperio

    # ---------------
    # Public Methods 
    # ---------------
    def apply_mask(
        self,
        mask: xr.DataArray | None = None,
        drop: bool = False
    ) -> Self:
        """
        Apply NEMO parent grid land-sea mask or combined land-sea & custom mask
        to variable defined on a NEMO model grid.

        Parameters
        ----------
        mask : xr.DataArray | None
            Boolean mask to apply to variable defined on NEMO model grid. Default is None
            meaning only land-sea mask of NEMO parent grid is applied.
        drop : bool, optional
            If True, coordinate labels that only correspond to False values of the condition
            are dropped from the result. Default is False.

        Returns
        -------
        NEMODataArray

        Examples
        --------
        Apply a custom boolean mask `my_mask` to sea surface temperature `tos_con` defined on
        scalar T-points in a NEMO model parent domain:

        >>> nemo['gridT/tos_con'].apply_mask(mask=my_mask)

        Apply custom mask `my_mask` to sea surface salinity `sos_abs` and drop masked data values:

        >>> nemo['gridT/sos_abs'].apply_mask(mask=my_mask, drop=True)

        See Also
        --------
        masked
        """
        # -- Validate Inputs -- #
        if mask is not None:
            if not isinstance(mask, xr.DataArray):
                raise TypeError("mask must be specified as an xarray.DataArray or None.")
        if mask is None:
            mask = self.mask
        else:
            mask = (self.mask & mask)

        # -- Apply mask & return NEMODataArray -- #
        if drop and isinstance(mask.data, dask.array.Array):
            warnings.warn(message="Indexing with a boolean dask array is not allowed. Mask will be materialised first using .load(). This may result in high memory usage for large masks.",
                          category=RuntimeWarning,
                          stacklevel=2
                          )
            mask = mask.load()

        result = self._da.where(mask, drop=drop)

        return self._wrap(result)

    def sel_like(
        self,
        other: Self | xr.DataArray,
    ) -> Self:
        """
        Return a new NEMODataArray whose data is given by matching the dimension
        index labels present in another NEMODataArray or xarray.DataArray.

        Parameters
        ----------
        other : NEMODataArray | xr.DataArray
            NEMODataArray or xarray.DataArray used to select dimension index labels.

        Returns
        -------
        NEMODataArray
             A new NEMODataArray with data selected according to dimension index labels
             of input object.

        Examples
        --------
        Indexing conservative temperature `thetao_con` defined on scalar T-points in a NEMO
        model parent domain to match a subset of the absolute salinity `so_abs`:

        >>> nda = nemo['gridT/so_abs'].sel(time_counter=slice("2020-01", "2024-01"), k=1)

        >>> nemo['gridT/thetao_con'].sel_like(nda)

        See Also
        --------
        sel
        """
        # -- Validate Inputs -- #
        if not (isinstance(other, NEMODataArray) or isinstance(other, xr.DataArray)):
            raise TypeError("other must be specified as a NEMODataArray or xarray.DataArray.")

        if isinstance(other, NEMODataArray):
            other = other.data

        # -- Index data using dimension labels of other object & return NEMODataArray -- #
        coord_dims = [self.t_name, self.k_name, self.j_name, self.i_name]
        d_dims = {}
        for dim in coord_dims:
            # Select only dimensions present in both objects:
            if (dim in other.coords) and (dim in self.data.coords):
                # Select only dimensions whose index labels do not already match:
                if other.coords[dim].size != self.data.coords[dim].size:
                    d_dims[dim] = other.coords[dim].data

        # Indexing only required if dimensions modified:
        if len(d_dims) > 0:
            result = self.data.sel(d_dims)
            return self._wrap(result)
        # Otherwise return original NEMODataArray:
        else:
            return self

    def weighted_mean(
        self,
        dims : list,
        mask: xr.DataArray | None = None,
        skipna : bool | None = None
    ) -> Self:
        """
        Calculate grid-aware weighted mean of a variable defined on a NEMO model grid.

        Parameters
        ----------
        dims : list
            Dimensions over which to apply weighted mean (e.g., ['i', 'j']).
        mask: xr.DataArray, optional
            Boolean mask identifying NEMO model grid points to be included (1)
            or neglected (0) from integration.
        skipna : bool | None
            If True, skip missing values (as marked by NaN).
            By default, only skips missing values for float dtypes.

        Returns
        -------
        NEMODataArray
            Grid-aware weighted mean of variable defined on a NEMO model grid.

        Examples
        --------
        Area-weighted mean of the sea surface temperature `tos_con` defined on scalar T-points
        in a NEMO model nested child domain:

        >>> nemo["gridT/1_gridT/tos_con"].weighted_mean(dims=["i", "j"], skipna=True)

        See Also
        --------
        masked_statistic
        """
        # -- Validate Input -- #
        if not isinstance(dims, list):
            raise TypeError("dims must be specified as a list.")
        if mask is not None:
            if not isinstance(mask, xr.DataArray):
                raise TypeError("mask must be an xarray.DataArray.")
            if any(dim not in self.dims for dim in mask.dims):
                raise ValueError(
                    f"mask must have dimensions subset from {self.dims}."
                )
        if skipna is not None:
            if not isinstance(skipna, bool):
                raise TypeError("skipna must be specified as a boolean or None.")

        # -- Calculate weighted mean & return NEMODataArray -- #
        da = self.apply_mask(mask=mask).data
        weight_dims = [dim.replace(self._dom_suffix, "") for dim in dims]
        weights = self._tree._get_weights(grid=self._grid, dims=weight_dims)
        result = da.weighted(weights).mean(dim=dims, skipna=skipna)
        result.name = f"wmean_{'_'.join(dims)}({self.name})"

        return self._wrap(result)

    def diff(
        self,
        dim: str,
        fillna: bool = False,
    ) -> Self:
        """
        Calculate the discrete difference of a variable along a given dimension
        (e.g., 'i', 'j', 'k') of a NEMO model grid.

        Parameters
        ----------
        dim : str
            Dimension over which to calculate the finite difference (e.g., 'i', 'j', 'k').
        fillna : bool, optional
            Fill NaN values in NEMODataArray with zeros prior to finite differencing. Default is False.
        Returns
        -------
        NEMODataArray
            Discrete difference of variable defined on a new NEMO model grid. For example, the
            discrete difference along the i-dimension of a scalar variable defined on a T-grid
            returns a NEMODataArray defined on the U-grid. 

        Examples
        --------
        Compute the difference of sea surface temperature `tos_con` values
        along the NEMO parent domain `j` dimension:

        >>> nemo['gridT/tos_con'].diff(dim="j")

        Compute the difference of sea surface temperature `tos_con` values along a regional subset of a global,
        zonally periodic domain NEMO parent domain `i` dimension:

        >>> nemo['gridT/tos_con'].sel(i=slice(10, 80)).diff(dim="i")

        Here, the zonal periodicity inherited from the NEMO model grid is automatically set to `False` since a subset
        of the global domain is no longer zonally periodic.

        Compute the discrete difference of absolute salinity `so_abs` values along the first NEMO
        nested child domain `k` dimension:

        >>> nemo['gridT/1_gridT/so_abs'].diff(dim="k")

        See Also
        --------
        derivative
        """
        # -- Validate Inputs -- #
        if not isinstance(dim, str):
            raise ValueError(
                "dim must be a string specifying dimension along which to calculate the gradient (e.g., 'i', 'i1', 'j', 'j1', 'k', 'k1')."
            )
        if dim not in self.dims:
            raise KeyError(
                f"dimension '{dim}' not found in {self.name or 'unnamed'} dimensions {self.dims}."
            )
        if not isinstance(fillna, bool):
            raise TypeError(
                "`fillna` must be specified as a boolean. Default is False."
            )

        # -- Get NEMO model grid properties -- #
        grid_paths = self._tree._get_grid_paths(dom=self._dom)

        # -- Calculate 1st-order discrete difference -- #
        if fillna:
            da = self.masked.data.fillna(value=0)
        else:
            da = self.masked.data

        if dim == self.k_name:
            match self._grid_suffix:
                case "w":
                    # W-grid located at k = 0.5, 1.5 ... Nk-0.5
                    result = da.diff(dim=self.k_name, n=1, label="lower")
                    # Fill final T-point [k=Nk] -> NaN:
                    result = result.pad({self.k_name: (0, 1)}, constant_values=np.nan)
                    result.coords[self.k_name] = (result.coords[self.k_name] + 0.5).astype(int)
                case _:
                    # T/U/V/F-grids located at k = 1, 2 ... Nk
                    result = da.diff(dim=self.k_name, n=1, label="lower")
                    # Fill initial W-point [k=0.5] -> NaN:
                    result = result.pad({self.k_name: (1, 0)}, constant_values=np.nan)
                    result.coords[self.k_name] = (result.coords[self.k_name].fillna(0) + 0.5)

        elif dim == self.i_name:
            match self._grid_suffix:
                case "t" | "v" | "w":
                    # T/V/W-grids located at i = 1, 2 ... Ni
                    if self.iperio:
                        result = da.roll({self.i_name: -1}) - da
                    else:
                        result = da.shift({self.i_name: -1}) - da
                    result.coords[self.i_name] = result.coords[self.i_name] + 0.5
                case "u" | "f":
                    # U/F-grids located at i = 1.5, 2.5 ... Ni+0.5
                    if self.iperio:
                        result = da - da.roll({self.i_name: 1})
                        result.coords[self.i_name] = result.coords[self.i_name] - 0.5
                    else:
                        # Fill initial U/F-point [i=0.5] -> NaN or 0:
                        result = da.pad({self.i_name: (1, 0)}, constant_values=0)
                        result = result.diff(dim=self.i_name, n=1, label="lower")
                        result.coords[self.i_name] = (result.coords[self.i_name].fillna(0.5) + 0.5).astype(int)

        elif dim == self.j_name:
            match self._grid_suffix:
                case "t" | "u" | "w":
                    # T/U/W-grids located at j = 1, 2 ... Nj
                    result = da.shift({self.j_name: -1}) - da
                    result.coords[self.j_name] = result.coords[self.j_name] + 0.5
                case "v" | "f":
                    # V/F-grids located at j = 1.5, 2.5 ... Nj+0.5
                    # Fill initial V/F-point [j=0.5] -> NaN or 0:
                    result = da.pad({self.j_name: (1, 0)}, constant_values=0)
                    result = result.diff(dim=self.j_name, n=1, label="lower")
                    result.coords[self.j_name] = (result.coords[self.j_name].fillna(0.5) + 0.5).astype(int)
        else:
            raise ValueError(f"Invalid dimension {dim}. Dimension must be one of (i{self._dom}, j{self._dom}, k{self._dom}).")

        # -- Updating NEMO model grid coordinates -- #
        geo_coords = [coord for coord in da.coords if coord not in (self.t_name, self.k_name, self.i_name, self.j_name)]

        # Determine new NEMO model grid type & dims:
        new_ijk_dims = [dim for dim in result.dims if dim != 'time_counter']
        new_grid_suffix, new_hgrid_suffix, _ = _NEMO_DIFF_MAP[(self._grid_suffix.upper(), dim.replace(self._dom_suffix, ""))]
        if len(new_grid_suffix) > 1:
            new_vgrid_suffix = new_grid_suffix[1]
        else:
            new_vgrid_suffix = new_grid_suffix

        # Define new geographical coordinates (glam, gphi, depth) based on new grid type:
        # NOTE: Subsets of NEMODataTree geographical coordinates are supported implicitly since the
        # (i, j, k) coordinates of the grid dimensions (i, j, k) can be used to align incoming (i.e., _tree) coordinates.
        new_coords = {
            f"{self._dom_prefix}glam{new_hgrid_suffix}": self._tree[grid_paths[f"grid{new_hgrid_suffix.upper()}"]][f"{self._dom_prefix}glam{new_hgrid_suffix}"],
            f"{self._dom_prefix}gphi{new_hgrid_suffix}": self._tree[grid_paths[f"grid{new_hgrid_suffix.upper()}"]][f"{self._dom_prefix}gphi{new_hgrid_suffix}"],
        }
        if self.k_name in new_ijk_dims:
            new_coords[f"{self._dom_prefix}depth{new_vgrid_suffix}"] = self._tree[grid_paths[f"grid{new_vgrid_suffix.upper()}"]][f"{self._dom_prefix}depth{new_vgrid_suffix}"]

        # -- Update DataArray properties -- #
        result = result.drop_vars(geo_coords)
        result = result.assign_coords(new_coords)
        result.name = f"diff_{dim}({self.name})"

        # -- Apply land-sea mask & return NEMODataArray -- #
        new_grid = f"{self._grid.replace(self._grid[-1], new_grid_suffix.upper())}"
        result = NEMODataArray(da=result, tree=self._tree, grid=new_grid)
        result = result.masked

        return result

    def derivative(
        self,
        dim: str,
        fillna: bool = False,
    ) -> Self:
        """
        Calculate the derivative of a variable along a given dimension
        (e.g., 'i', 'j', 'k') of a NEMO model grid.

        Parameters
        ----------
        dim : str
            Dimension over which to calculate the derivative (e.g., 'i', 'j', 'k').
        fillna : bool, optional
            Fill NaN values in NEMODataArray with zeros prior to finite differencing. Default is False.

        Returns
        -------
        NEMODataArray
            Derivative of variable defined on a NEMO model grid.

        Examples
        --------
        Compute the derivative of sea surface temperature `tos_con` values
        along the NEMO parent domain `j` dimension:

        >>> nemo['gridT/tos_con'].derivative(dim="j")

        Compute the derivative of sea surface temperature `tos_con` values along a regional subset of a global,
        zonally periodic domain NEMO parent domain `i` dimension:

        >>> nemo['gridT/tos_con'].sel(i=slice(10, 80)).derivative(dim="i")

        Here, the zonal periodicity inherited from the NEMO model grid is automatically set to `False` since a subset
        of the global domain is no longer zonally periodic.

        Compute the discrete difference of absolute salinity `so_abs` values along the first NEMO
        nested child domain `k` dimension:

        >>> nemo['gridT/1_gridT/so_abs'].derivative(dim="k")

        See Also
        --------
        diff
        """
        # -- Validate Inputs -- #
        if not isinstance(dim, str):
            raise ValueError(
                "dim must be a string specifying dimension along which to calculate the gradient (e.g., 'i', 'i1', 'j', 'j1', 'k', 'k1')."
            )
        if dim not in self.dims:
            raise KeyError(
                f"dimension '{dim}' not found in {self.name or 'unnamed'} dimensions {self.dims}."
            )
        if not isinstance(fillna, bool):
            raise TypeError(
                "`fillna` must be specified as a boolean. Default is False."
            )


        # -- Calculate 1st-discrete derivative along dimension -- #
        # Determine NEMO model grid type and scale factors for derivative:
        new_grid_suffix, _, new_grid_weights = _NEMO_DIFF_MAP[(self._grid_suffix.upper(),
                                                               dim.replace(self._dom_suffix, "")
                                                               )]
        # Define path to derivative NEMO model grid:
        new_grid = f"{self._grid.replace(self._grid[-1], new_grid_suffix.upper())}"

        # Collect derivative grid scale factors (e.g., e1u, e3t etc.):
        try:
            weights = self._tree[new_grid][new_grid_weights]
        except KeyError as e:
            raise KeyError(
                f"NEMO model grid: '{new_grid}' does not contain grid scale factor '{new_grid_weights}' required to calculate derivatives along the {dim}-dimension."
            ) from e

        # Calculate 1st-finite difference along dimension:
        da = self.diff(dim=dim, fillna=fillna)

        # Calculate derivative (i.e., diff(var) / e{1/2/3}{t/u/v/w}):
        if dim in [self.k_name]:
            # Vertical derivative [k increasing downward]:
            result = - da.data / weights
        else:
            # Horizontal derivative [i/j increasing eastward/northward]:
            result = da.data / weights

        # -- Update DataArray properties & return NEMODataArray -- #
        result.name = f"d({self.name})/d{dim}"
        result = NEMODataArray(da=result, tree=self._tree, grid=new_grid)

        return result


    def integral(
        self,
        dims: list,
        cum_dims: list | None = None,
        dir: str | None = None,
        mask: xr.DataArray | None = None,
    ) -> Self:
        """
        Integrate variable along one or more dimensions of a NEMO model grid.

        Parameters
        ----------
        dims : list
            Dimensions over which to integrate (e.g., ['i', 'k']).
        cum_dims : list, optional
            Dimensions over which to cumulatively integrate (e.g., ['k']).
            Specified dimensions must also be included in `dims`.
        dir : str, optional
            Direction of cumulative integration. Options are '+1' (along
            increasing cum_dims) or '-1' (along decreasing cum_dims).
        mask: xr.DataArray, optional
            Boolean mask identifying NEMO model grid points to be included (1)
            or neglected (0) from integration.

        Returns
        -------
        NEMODataArray
            Variable integrated along specified dimensions of the NEMO model grid.


        Examples
        --------
        Compute the integral of conservative temperature `thetao_con` along the vertical
        `k` dimension in the NEMO parent domain:

        >>> nemo['gridT/thetao_con'].integral(dims=["k"])

        Compute the vertical meridional overturning stream function from the meridional
        velocity `vo` (zonally integrated meridional velocity accumulated with increasing
        depth):

        >>> nemo['gridV/vo'].integral(dims=["i", "k"],
        ...                           cum_dims=["k"],
        ...                           dir="+1",
        ...                           )

        See Also
        --------
        depth_integral
        """
        # -- Validate Input -- #
        if cum_dims is not None:
            for dim in cum_dims:
                if dim not in dims:
                    raise ValueError(
                        f"cumulative integration dimension '{dim}' not included in `dims`."
                    )
            if dir not in ["+1", "-1"]:
                raise ValueError(
                    f"invalid direction of cumulative integration '{dir}'. Expected '+1' or '-1'."
                )
        if mask is not None:
            if not isinstance(mask, xr.DataArray):
                raise ValueError("mask must be an xarray.DataArray.")
            if any(dim not in self.dims for dim in mask.dims):
                raise ValueError(
                    f"mask must have dimensions subset from {self.dims}."
                )

        # -- Collect variable, weights & mask -- #
        da = self.apply_mask(mask=mask).data
        weights = self._tree._get_weights(grid=self._grid, dims=dims)

        # -- Perform integration -- #
        if cum_dims is not None:
            sum_dims = [dim for dim in dims if dim not in cum_dims]
            if dir == "+1":
                # Cumulative integration along ordered dimension:
                result = (
                    da.weighted(weights)
                    .sum(dim=sum_dims, skipna=True)
                    .cumsum(dim=cum_dims, skipna=True)
                )
            elif dir == "-1":
                # Cumulative integration along reversed dimension:
                result = (
                    da.weighted(weights)
                    .sum(dim=sum_dims, skipna=True)
                    .reindex({dim: self.data[dim][::-1] for dim in cum_dims})
                    .cumsum(dim=cum_dims, skipna=True)
                )
        else:
            result = da.weighted(weights).sum(dim=dims, skipna=True)

        # -- Apply land-sea mask & return NEMODataArray -- #
        result.name = f"integral_{''.join(dims)}({self.name})"
        result = self._wrap(result)

        return result

    def depth_integral(
        self,
        limits: tuple[int | float]
    ) -> Self:
        """
        Integrate a variable in depth coordinates between two limits.

        Parameters
        ----------
        limits : tuple[int | float]
            Limits of depth integration given as a tuple of the form
            (depth_min, depth_max) where depth_min and depth_max are
            the lower and upper limits of vertical integration, respectively.

        Returns
        -------
        NEMODataArray
            Vertical integral of chosen variable between two depth surfaces (depth_min, depth_max).

        Examples
        --------
        Vertically integrate the conservative temperature variable `thetao_con` defined in a
        NEMO model parent domain from the sea surface to 100 m depth:

        >>> nemo["gridT/thetao_con"].depth_integral(limits=(0, 100))

        See Also
        --------
        integral
        """
        # -- Validate Input -- #
        if (not isinstance(limits, tuple)) | (len(limits) != 2):
            raise TypeError(
                "depth limits of integration should be given by a tuple of the form (depth_min, depth_max)"
            )
        if (limits[0] < 0) | (limits[1] < 0):
            raise ValueError("depth limits of integration must be non-negative.")
        if limits[0] >= limits[1]:
            raise ValueError(
                "lower depth limit must be less than upper depth limit."
            )

        # -- Define input variables -- #
        var_in = self.masked.data
        e3_in = self.metrics["e3"].masked.data

        # -- Vertically integrate w.r.t depth -- #
        result = xr.apply_ufunc(
            compute_depth_integral,
            e3_in,
            var_in,
            np.array([limits[1]]),
            np.array([limits[0]]),
            input_core_dims=[[self.k_name], [self.k_name], [], []],
            output_core_dims=[["k_new"]],
            dask="parallelized",
            output_dtypes=[var_in.dtype],
            dask_gufunc_kwargs={"output_sizes": {"k_new": 1}},
        )

        # -- Define integral variable DataArray -- #
        dim_list = [dim for dim in [self.t_name, "k_new", self.j_name, self.i_name] if (dim is not None) and (dim in result.dims)]
        result = result.transpose(*dim_list).squeeze(dim="k_new")
        result.name = f"integral_z({self.name})"

        # -- Apply land-sea mask & return NEMODataArray -- #
        result = self._wrap(result)

        return result

    def masked_statistic(
        self,
        lon_poly: list | np.ndarray,
        lat_poly: list | np.ndarray,
        statistic: str,
        dims: list,
        skipna: bool | None = None,
    ) -> Self:
        """
        Compute masked statistic of a variable defined on a NEMO model grid.

        Parameters
        ----------
        lon_poly : list | np.ndarray
            Longitudes of closed polygon.
        lat_poly : list | np.ndarray
            Latitudes of closed polygon.
        statistic : str
            Name of the statistic to calculate (e.g., 'mean', 'weighted_mean' 'sum').
        dims : list
            Dimensions over which to apply statistic (e.g., ['i', 'j']).
        skipna : bool | None
            If True, skip missing values (as marked by NaN). By default, only skips missing values for float dtypes.

        Returns
        -------
        NEMODataArray
            Masked statistic of variable defined on a NEMO model grid.

        Examples
        --------
        Compute the grid cell area-weighted mean sea surface temperature `tos_con` for a
        region enclosed in a polygon defined by `lon_poly` and `lat_poly` in a NEMO nested
        child domain:

        >>> nemo["gridT/1_gridT/tos_con"].masked_statistic(lon_poly=lon_poly,
        ...                                                lat_poly=lat_poly,
        ...                                                statistic="weighted_mean",
        ...                                                dims=["i", "j"]
        ...                                                )

        See Also
        --------
        NEMODataTree.binned_statistic
        """
        # -- Validate input -- #
        if not isinstance(statistic, str):
            raise TypeError("statistic must be specified as a string.")
        if not isinstance(dims, list):
            raise TypeError("dims must be specified as a list.")
        if skipna is not None:
            if not isinstance(skipna, bool):
                raise TypeError("skipna must be specified as a boolean or None.")

        # -- Create polygon mask using coordinates -- #
        mask_poly = self._tree.mask_with_polygon(
            lon_poly=lon_poly, lat_poly=lat_poly, grid=self._grid
        )

        # -- Apply masks & calculate statistic -- #
        da = self.apply_mask(mask=mask_poly).data

        match statistic:
            case "mean":
                result = da.mean(dim=dims, skipna=skipna)

            case "weighted_mean":
                weight_dims = [dim.replace(self._dom_suffix, "") for dim in dims]
                weights = self._tree._get_weights(grid=self._grid, dims=weight_dims)
                result = da.weighted(weights).mean(dim=dims, skipna=skipna)

            case "min":
                result = da.min(dim=dims, skipna=skipna)

            case "max":
                result = da.max(dim=dims, skipna=skipna)

            case "sum":
                result = da.sum(dim=dims, skipna=skipna)

            case _:
                raise ValueError(
                    f"Unsupported statistic '{statistic}'. Supported statistics are: 'mean', 'weighted_mean', 'min', 'max', 'sum'."
                )

        # -- Update DataArray properties -- #
        result.name = f"masked_{statistic}({self.name})"
        result = self._wrap(result)

        return result

    def interp_to(
        self,
        to: str,
    ) -> Self:
        """
        Linearly interpolate variable to a neighbouring horizontal grid.

        For flux variables defined at U/V-points, the specified variable
        is first weighted by grid cell face area prior to linear interpolation,
        and is then normalised by the area of the target grid cell face following
        interpolation.

        Parameters
        ----------
        to : str
            Suffix of neighbouring horizontal NEMO model grid to linear interpolate
            variable to. Options are 'T', 'U', 'V', 'F'.

        Returns
        -------
        NEMODataArray
            Variable linearly interpolated onto a neighbouring horizontal grid.

        Examples
        --------
        Linearly interpolate conservative temperature `thetao_con` defined on scalar T-points
        to neighbouring V-points in a NEMO model parent domain:

        >>> nemo['gridT/thetao_con'].interp_to(to='V')

        See Also
        --------
        transform_vertical_grid
        """
        # -- Validate input -- #
        if not isinstance(to, str):
            raise TypeError(f"'to' must be a string, got {type(to)}.")
        if to not in ["T", "U", "V", "F"]:
            raise ValueError(f"'to' must be one of ['T', 'U', 'V', 'F'], got {to}.")

        # -- Get NEMO model grid properties -- #
        ijk_names = self._tree._get_ijk_names(grid=self._grid)
        target_grid = f"{self._grid.replace(self._grid[-1], to)}"

        # -- Collect variable grid scale factors -- #
        if self._grid_suffix.upper() in ["U", "V"]:
            weight_dims = (
                [self.k_name, self.j_name] if self._grid_suffix.upper() == "U" else [self.k_name, self.i_name]
            )
            if f"{self._dom_prefix}depth{self._grid_suffix}" in self.coords:
                # 3-D variables - weight by grid cell face area:
                weights = self._tree._get_weights(grid=self._grid, dims=weight_dims, fillna=False)
                target_weights = self._tree._get_weights(
                    grid=target_grid, dims=weight_dims, fillna=False
                )
            else:
                # 2-D variables - weight by grid cell width:
                weights = self._tree._get_weights(grid=self._grid, dims=weight_dims[1], fillna=False)
                target_weights = self._tree._get_weights(
                    grid=target_grid, dims=weight_dims[1], fillna=False
                )
            da = self.masked.data * weights
        else:
            # Scalar variables:
            da = self.masked.data

        # -- Linearly interpolate variable -- #
        result = interpolate_grid(
            da=da,
            mask=self.mask,
            source_grid=self._grid_suffix.upper(),
            target_grid=to,
            iperio=self.iperio,
            ijk_names=ijk_names,
        )

        # -- Updating NEMO model grid coordinates -- #
        geo_coords = [coord for coord in da.coords if coord not in (self.t_name, self.k_name, self.i_name, self.j_name)]

        # Define new geographical coordinates (glam, gphi, depth) based on new grid type:
        # NOTE: Subsets of NEMODataTree geographical coordinates are supported implicitly since the
        # (i, j, k) coordinates of the grid dimensions (i, j, k) can be used to align incoming (i.e., _tree) coordinates.
        new_coords = {
            f"{self._dom_prefix}glam{to.lower()}": self._tree[target_grid][f"{self._dom_prefix}glam{to.lower()}"],
            f"{self._dom_prefix}gphi{to.lower()}": self._tree[target_grid][f"{self._dom_prefix}gphi{to.lower()}"],
        }
        if self.k_name in result.coords:
            new_coords[f"{self._dom_prefix}depth{to.lower()}"] = self._tree[target_grid][f"{self._dom_prefix}depth{to.lower()}"]

        # -- Update DataArray properties -- #
        result = result.drop_vars(geo_coords)
        result = result.assign_coords(new_coords)
        # Retain original variable name:
        result.name = da.name

        # Reorder dimensions (time_counter, [k], j, i):
        var_dims = [dim for dim in [self.t_name, self.k_name, self.j_name, self.i_name] if (dim is not None) and (dim in result.dims)]
        result = result.transpose(*var_dims)

        # Normalise by target grid cell weights for flux variables:
        if self._grid_suffix.upper() in ["U", "V"]:
            result = result / target_weights

        # -- Apply land-sea mask & return NEMODataArray -- #
        result = NEMODataArray(da=result, tree=self._tree, grid=target_grid)
        result = result.masked

        return result

    def transform_vertical_grid(
        self, e3_new: xr.DataArray
    ) -> xr.Dataset:
        """
        Transform variable defined on a NEMO model grid to a new vertical grid using conservative interpolation.

        Parameters
        ----------
        e3_new : xarray.DataArray
            Grid cell thicknesses of the new vertical grid.
            Must be a 1-dimensional xarray.DataArray with
            dimension 'k_new'.

        Returns
        -------
        xr.Dataset
            Variable defined at the centre of each vertical
            grid cell on the new grid, and vertical grid cell
            thicknesses adjusted for model bathymetry.

        Examples
        --------
        Transform the conservative temperature variable `thetao_con` defined in a
        NEMO model parent domain from it's native 75 unevenly-spaced z-levels to
        regularly spaced z-levels at 200 m intervals:

        >>> e3t_target = xr.DataArray(np.repeat(200.0, 30), dims=['k_new'])

        >>> nemo['gridT/thetao_con'].transform_vertical_grid(e3_new=e3t_target)

        See Also
        --------
        transform_to
        """
        # -- Validate Input -- #
        if e3_new.dims != ("k_new",) or (e3_new.ndim != 1):
            raise ValueError(
                "e3_new must be a 1-dimensional xarray.DataArray with dimension 'k_new'."
            )

        # -- Define input variables -- #
        var_in = self.masked.data
        e3_in = self.metrics["e3"].masked.data

        # Ensure total depth of new vertical grid >= total depth of NEMO model vertical grid:
        depth_max = self._tree[self._grid][f"{self._dom_prefix}depth{self._grid_suffix}"].max(self.k_name)
        if e3_new.sum(dim="k_new") < depth_max:
            raise ValueError(
                f"e3_new must sum to at least the maximum depth ({depth_max.item()} m) of the original vertical grid."
            )

        # -- Transform variable to target vertical grid -- #
        var_out, e3_out = xr.apply_ufunc(
            transform_vertical_coords,
            e3_in,
            var_in,
            e3_new.astype(e3_in.dtype),
            input_core_dims=[[self.k_name], [self.k_name], ["k_new"]],
            output_core_dims=[["k_new"], ["k_new"]],
            dask="parallelized",
            output_dtypes=[var_in.dtype, e3_in.dtype],
            dask_gufunc_kwargs={"output_sizes": {"k_new": e3_new.sizes["k_new"]}},
        )

        # -- Construct transformed variable Dataset -- #
        var_dims = [dim for dim in [self.t_name, "k_new", self.j_name, self.i_name] if (dim is not None) and (dim in var_out.dims)]
        var_out = var_out.transpose(*var_dims)

        e3_dims = [dim for dim in [self.t_name, "k_new", self.j_name, self.i_name] if (dim is not None) and (dim in e3_out.dims)]
        e3_out = e3_out.transpose(*e3_dims)

        result = xr.Dataset(
            data_vars={self.name: var_out, f"e3{self._grid_suffix}_new": e3_out},
            coords={
                f"depth{self._grid_suffix}_new": ("k_new", e3_new.cumsum(dim="k_new").data)
            },
        )

        return result

    def to_xesmf(
            self, mask: bool = True
        ) -> xr.Dataset:
        """
        Create an xESMF-compatible dataset from a variable stored in a NEMODataArray
        with the following properties:

        - lon (longitudes of grid cell centers)
        - lat (latitudes of grid cell centers)
        - lon_b (longitudes of grid cell boundaries)
        - lat_b (latitudes of grid cell boundaries)
        - mask (boolean land-sea mask identifying valid grid cells) [ Optional ]

        Only NEMODataArrays defined on NEMO model T-grids are currently supported
        for transformation to xESMF-compatible datasets.

        Users can use interp_to(to='T') to linearly interpolate NEMO variables defined
        on neighbouring U/V/F-grids to T-grids prior to creating an xESMF-compatible
        dataset.

        Parameters
        ----------
        mask : bool, optional
            Whether to include the land-sea mask in the resulting dataset.
            Default is True.

        Returns
        -------
        xr.Dataset
            Dataset containing variable and associated geographical coordinates
            formatted for use with xESMF.

        Examples
        --------
        Create an xESMF-compatible dataset from the conservative temperature variable
        `thetao_con` defined on a NEMO model parent domain T-grid:

        >>> nemo['gridT/thetao_con'].to_xesmf(mask=True)

        See Also
        --------
        interp_to
        """
        # -- Validate Input -- #
        if not isinstance(mask, bool):
            raise ValueError("mask must be specified as a boolean. Default is True.")
        if self._grid_suffix.upper() != "T":
            raise ValueError("`to_xesmf()` accessor only supports variables defined on a NEMO model T-grid." \
            "Use `interp_to(to='T')` to linearly interpolate variables onto the T-grid first."
            )

        # -- Create xESMF-compatible dataset -- #
        ds = create_xesmf_dataset(nda=self, mask=mask)

        return ds

    # ----------------
    # Binary Operators 
    # ----------------
    def _binary_op(self, other: Self | xr.DataArray | int | float, op) -> Self:
        if isinstance(other, NEMODataArray):
            if self.grid != other.grid:
                raise ValueError(
                    f"Cannot perform binary operation between NEMODataArrays defined on different NEMO grids ({self.name} -> {self.grid}, {other.name} -> {other.grid})."
                )
        result = op(self._da, self._unwrap(other))
        return self._wrap(result)

    def _rbinary_op(self, other: Self | xr.DataArray | int | float, op) -> Self:
        if isinstance(other, NEMODataArray):
            if self.grid != other.grid:
                raise ValueError(
                    f"Cannot perform binary operation between NEMODataArrays defined on different NEMO grids ({self.name} -> {self.grid}, {other.name} -> {other.grid})."
                )
        result = op(self._unwrap(other), self._da)
        return self._wrap(result)

    def __add__(self, other: Self | xr.DataArray | int | float) -> Self:
        return self._binary_op(other, operator.add)
    def __radd__(self, other: Self | xr.DataArray | int | float) -> Self:
        return self._rbinary_op(other, operator.add)

    def __sub__(self, other: Self | xr.DataArray | int | float) -> Self:
        return self._binary_op(other, operator.sub)
    def __rsub__(self, other: Self | xr.DataArray | int | float) -> Self:
        return self._rbinary_op(other, operator.sub)

    def __mul__(self, other: Self | xr.DataArray | int | float) -> Self:
        return self._binary_op(other, operator.mul)
    def __rmul__(self, other: Self | xr.DataArray | int | float) -> Self:
        return self._rbinary_op(other, operator.mul)

    def __truediv__(self, other: Self | xr.DataArray | int | float) -> Self:
        return self._binary_op(other, operator.truediv)
    def __rtruediv__(self, other: Self | xr.DataArray | int | float) -> Self:
        return self._rbinary_op(other, operator.truediv)

    # ----------------
    # Utility Methods 
    # ----------------
    def _wrap(self, da: xr.DataArray) -> Self:
        """
        Wrap xarray.DataArray to preserve existing NEMO domain & grid attributes.
        """
        return NEMODataArray(da=da, tree=self._tree, grid=self._grid)

    def _unwrap(self, other: Any) -> Any:
        """
        Unwrap NEMODataArray to access underlying xarray.DataArray.
        """
        if isinstance(other, NEMODataArray):
            return other.data
        return other

    def _conditional_wrap(self, result: Any) -> Any:
        """
        Conditionally wrap result of methods when returning xarray.DataArray.
        """
        if isinstance(result, xr.DataArray):
            try:
                # Attempt to return NEMODataArray:
                return self._wrap(result)
            except Exception:
                # Otherwise, return unwrapped xarray.DataArray:
                return result 
        return result

    def __setitem__(self, key, value):
        """
        Delegate item assignment to xarray.DataArray.
        """
        self._da[key] = self._unwrap(value)

    def __getitem__(self, key):
        """
        Delegate indexing to xarray.DataArray and attempt to return
        result as NEMODataArray.
        """
        result = self._da[key]
        return self._conditional_wrap(result)

    def __getattr__(self, name: str) -> Any:
        """
        Delegate attributes to xarray.DataArray and attempt to return
        result as NEMODataArray.
        """
        attr = getattr(self._da, name)

        if callable(attr):
            @wraps(attr)
            def wrapper(*args: Any, **kwargs: Any) -> Any:
                result = attr(*args, **kwargs)
                return self._conditional_wrap(result)
            return wrapper

        return attr

    def __repr__(self) -> str:
        return (
            f"<NEMODataTree '{self._tree.name or 'unnamed'}'>\n"
            f"  <NEMODataArray '{self.name or 'unnamed'}' (Domain: '{self._dom}', "
            f"Grid: '{self._grid}', Grid Type: '{self._grid_suffix.upper()}')>\n\n"
            f"{repr(self._da)}"
        )

    def _repr_html_(self) -> str:
        banner = f"""
        <div style="margin-bottom: 10px;">
            <div">
                <span style="color:#6a737d;">NEMODataTree</span> '{self._tree.name or "unnamed"}'
            </div>

            <div style="margin-left: 18px; margin-top:4px;">
                <span style="color:#6a737d;">NEMODataArray</span> '{self.name or "unnamed"}'
                <div>├─ <strong>Domain:</strong> '{self._dom}'</div>
                <div>├─ <strong>Grid Path:</strong> '{self._grid}'</div>
                <div>└─ <strong>Grid Type:</strong> '{self._grid_suffix.upper()}'</div>
            </div>
        </div>
        """

        return f"""
        <div style="border:1px solid #ddd; padding:12px; border-radius:8px;">
            {banner}
            {self.data._repr_html_()}
        </div>
        """

data property

data

Access underlying xarray.DataArray.

grid property

grid

Access path to parent NEMO model grid.

grid_type property

grid_type

Access type of parent NEMO model grid (e.g. 't', 'u', etc.).

iperio property

iperio

Zonal periodicity of variable defined on NEMO model grid.

mask property

mask

Access variable land-sea mask for the parent NEMO model grid.

masked property

masked

Apply land-sea mask to variable defined on NEMO model grid.

metrics property

metrics

Access grid scale factors for the parent NEMO model grid.

__getattr__

__getattr__(name)

Delegate attributes to xarray.DataArray and attempt to return result as NEMODataArray.

Source code in nemo_cookbook/nemodataarray.py

def __getattr__(self, name: str) -> Any:
    """
    Delegate attributes to xarray.DataArray and attempt to return
    result as NEMODataArray.
    """
    attr = getattr(self._da, name)

    if callable(attr):
        @wraps(attr)
        def wrapper(*args: Any, **kwargs: Any) -> Any:
            result = attr(*args, **kwargs)
            return self._conditional_wrap(result)
        return wrapper

    return attr

__getitem__

__getitem__(key)

Delegate indexing to xarray.DataArray and attempt to return result as NEMODataArray.

Source code in nemo_cookbook/nemodataarray.py

def __getitem__(self, key):
    """
    Delegate indexing to xarray.DataArray and attempt to return
    result as NEMODataArray.
    """
    result = self._da[key]
    return self._conditional_wrap(result)

__init__

__init__(da, tree, grid)

Create a NEMODataArray from an xarray.DataArray and NEMODataTree.

Parameters:

Name	Type	Description	Default
da	DataArray	Variable defined on a NEMO model grid.	required
tree	NEMODataTree	NEMODataTree to which the variable belongs.	required
grid	str	Path to NEMO model grid where variable is defined (e.g., 'gridT').	required

Returns:

Type	Description
NEMODataArray	Grid-aware xarray.DataArray storing variable defined on NEMO model grid.

Source code in nemo_cookbook/nemodataarray.py

def __init__(self, da: xr.DataArray, tree: NEMODataTree, grid: str):
    """
    Create a NEMODataArray from an xarray.DataArray and NEMODataTree.

    Parameters
    ----------
    da : xr.DataArray
        Variable defined on a NEMO model grid.
    tree: NEMODataTree
        NEMODataTree to which the variable belongs.
    grid: str
        Path to NEMO model grid where variable is defined (e.g., 'gridT').

    Returns
    -------
    NEMODataArray
        Grid-aware xarray.DataArray storing variable defined on NEMO model grid.
    """
    # -- Validate Inputs -- #
    if not isinstance(da, xr.DataArray):
        raise TypeError("da must be specified as an xarray.DataArray.")
    else:
        self._da = da

    if not isinstance(tree, xr.DataTree):
        raise TypeError("tree must be specified as a NEMODataTree.")
    else:
        self._tree = tree

    if not isinstance(grid, str):
        raise TypeError("grid must be specified as a string.")
    else:
        self._grid = grid

    # -- Validate DataArray compatibility -- #
    validate_nemo_dataarray(da=self._da, grid=self._grid, tree=self._tree)

    # -- Assign NEMO domain number, grid path and grid type -- #
    dom, dom_prefix, dom_suffix, grid_suffix = self._tree._get_properties(grid=self._grid, infer_dom=True)
    self._dom = dom
    self._dom_prefix = dom_prefix
    self._dom_suffix = dom_suffix
    self._grid_suffix = grid_suffix

    # -- Assign NEMO dimension names -- #
    # Spatial dimensions:
    ijk_names = self._tree._get_ijk_names(grid=self._grid)
    self.i_name = ijk_names["i"]
    self.j_name = ijk_names["j"]
    self.k_name = ijk_names["k"]
    # Temporal dimension except for time-independent variables:
    # Use coords to handle single time slices of time-dependent variables.
    t_list = [dim for dim in self._da.coords if "time" in dim]
    self.t_name = t_list[0] if len(t_list) != 0 else None

__setitem__

__setitem__(key, value)

Delegate item assignment to xarray.DataArray.

Source code in nemo_cookbook/nemodataarray.py

def __setitem__(self, key, value):
    """
    Delegate item assignment to xarray.DataArray.
    """
    self._da[key] = self._unwrap(value)

apply_mask

apply_mask(mask=None, drop=False)

Apply NEMO parent grid land-sea mask or combined land-sea & custom mask to variable defined on a NEMO model grid.

Parameters:

Name	Type	Description	Default
mask	DataArray | None	Boolean mask to apply to variable defined on NEMO model grid. Default is None meaning only land-sea mask of NEMO parent grid is applied.	None
drop	bool	If True, coordinate labels that only correspond to False values of the condition are dropped from the result. Default is False.	False

Returns:

Type	Description
NEMODataArray

Examples:

Apply a custom boolean mask my_mask to sea surface temperature tos_con defined on scalar T-points in a NEMO model parent domain:

>>> nemo['gridT/tos_con'].apply_mask(mask=my_mask)

Apply custom mask my_mask to sea surface salinity sos_abs and drop masked data values:

>>> nemo['gridT/sos_abs'].apply_mask(mask=my_mask, drop=True)

See Also

masked

Source code in nemo_cookbook/nemodataarray.py

def apply_mask(
    self,
    mask: xr.DataArray | None = None,
    drop: bool = False
) -> Self:
    """
    Apply NEMO parent grid land-sea mask or combined land-sea & custom mask
    to variable defined on a NEMO model grid.

    Parameters
    ----------
    mask : xr.DataArray | None
        Boolean mask to apply to variable defined on NEMO model grid. Default is None
        meaning only land-sea mask of NEMO parent grid is applied.
    drop : bool, optional
        If True, coordinate labels that only correspond to False values of the condition
        are dropped from the result. Default is False.

    Returns
    -------
    NEMODataArray

    Examples
    --------
    Apply a custom boolean mask `my_mask` to sea surface temperature `tos_con` defined on
    scalar T-points in a NEMO model parent domain:

    >>> nemo['gridT/tos_con'].apply_mask(mask=my_mask)

    Apply custom mask `my_mask` to sea surface salinity `sos_abs` and drop masked data values:

    >>> nemo['gridT/sos_abs'].apply_mask(mask=my_mask, drop=True)

    See Also
    --------
    masked
    """
    # -- Validate Inputs -- #
    if mask is not None:
        if not isinstance(mask, xr.DataArray):
            raise TypeError("mask must be specified as an xarray.DataArray or None.")
    if mask is None:
        mask = self.mask
    else:
        mask = (self.mask & mask)

    # -- Apply mask & return NEMODataArray -- #
    if drop and isinstance(mask.data, dask.array.Array):
        warnings.warn(message="Indexing with a boolean dask array is not allowed. Mask will be materialised first using .load(). This may result in high memory usage for large masks.",
                      category=RuntimeWarning,
                      stacklevel=2
                      )
        mask = mask.load()

    result = self._da.where(mask, drop=drop)

    return self._wrap(result)

depth_integral

depth_integral(limits)

Integrate a variable in depth coordinates between two limits.

Parameters:

Name	Type	Description	Default
limits	tuple[int | float]	Limits of depth integration given as a tuple of the form (depth_min, depth_max) where depth_min and depth_max are the lower and upper limits of vertical integration, respectively.	required

Returns:

Type	Description
NEMODataArray	Vertical integral of chosen variable between two depth surfaces (depth_min, depth_max).

Examples:

Vertically integrate the conservative temperature variable thetao_con defined in a NEMO model parent domain from the sea surface to 100 m depth:

>>> nemo["gridT/thetao_con"].depth_integral(limits=(0, 100))

See Also

integral

Source code in nemo_cookbook/nemodataarray.py

def depth_integral(
    self,
    limits: tuple[int | float]
) -> Self:
    """
    Integrate a variable in depth coordinates between two limits.

    Parameters
    ----------
    limits : tuple[int | float]
        Limits of depth integration given as a tuple of the form
        (depth_min, depth_max) where depth_min and depth_max are
        the lower and upper limits of vertical integration, respectively.

    Returns
    -------
    NEMODataArray
        Vertical integral of chosen variable between two depth surfaces (depth_min, depth_max).

    Examples
    --------
    Vertically integrate the conservative temperature variable `thetao_con` defined in a
    NEMO model parent domain from the sea surface to 100 m depth:

    >>> nemo["gridT/thetao_con"].depth_integral(limits=(0, 100))

    See Also
    --------
    integral
    """
    # -- Validate Input -- #
    if (not isinstance(limits, tuple)) | (len(limits) != 2):
        raise TypeError(
            "depth limits of integration should be given by a tuple of the form (depth_min, depth_max)"
        )
    if (limits[0] < 0) | (limits[1] < 0):
        raise ValueError("depth limits of integration must be non-negative.")
    if limits[0] >= limits[1]:
        raise ValueError(
            "lower depth limit must be less than upper depth limit."
        )

    # -- Define input variables -- #
    var_in = self.masked.data
    e3_in = self.metrics["e3"].masked.data

    # -- Vertically integrate w.r.t depth -- #
    result = xr.apply_ufunc(
        compute_depth_integral,
        e3_in,
        var_in,
        np.array([limits[1]]),
        np.array([limits[0]]),
        input_core_dims=[[self.k_name], [self.k_name], [], []],
        output_core_dims=[["k_new"]],
        dask="parallelized",
        output_dtypes=[var_in.dtype],
        dask_gufunc_kwargs={"output_sizes": {"k_new": 1}},
    )

    # -- Define integral variable DataArray -- #
    dim_list = [dim for dim in [self.t_name, "k_new", self.j_name, self.i_name] if (dim is not None) and (dim in result.dims)]
    result = result.transpose(*dim_list).squeeze(dim="k_new")
    result.name = f"integral_z({self.name})"

    # -- Apply land-sea mask & return NEMODataArray -- #
    result = self._wrap(result)

    return result

derivative

derivative(dim, fillna=False)

Calculate the derivative of a variable along a given dimension (e.g., 'i', 'j', 'k') of a NEMO model grid.

Parameters:

Name	Type	Description	Default
dim	str	Dimension over which to calculate the derivative (e.g., 'i', 'j', 'k').	required
fillna	bool	Fill NaN values in NEMODataArray with zeros prior to finite differencing. Default is False.	False

Returns:

Type	Description
NEMODataArray	Derivative of variable defined on a NEMO model grid.

Examples:

Compute the derivative of sea surface temperature tos_con values along the NEMO parent domain j dimension:

>>> nemo['gridT/tos_con'].derivative(dim="j")

Compute the derivative of sea surface temperature tos_con values along a regional subset of a global, zonally periodic domain NEMO parent domain i dimension:

>>> nemo['gridT/tos_con'].sel(i=slice(10, 80)).derivative(dim="i")

Here, the zonal periodicity inherited from the NEMO model grid is automatically set to False since a subset of the global domain is no longer zonally periodic.

Compute the discrete difference of absolute salinity so_abs values along the first NEMO nested child domain k dimension:

>>> nemo['gridT/1_gridT/so_abs'].derivative(dim="k")

See Also

diff

Source code in nemo_cookbook/nemodataarray.py

def derivative(
    self,
    dim: str,
    fillna: bool = False,
) -> Self:
    """
    Calculate the derivative of a variable along a given dimension
    (e.g., 'i', 'j', 'k') of a NEMO model grid.

    Parameters
    ----------
    dim : str
        Dimension over which to calculate the derivative (e.g., 'i', 'j', 'k').
    fillna : bool, optional
        Fill NaN values in NEMODataArray with zeros prior to finite differencing. Default is False.

    Returns
    -------
    NEMODataArray
        Derivative of variable defined on a NEMO model grid.

    Examples
    --------
    Compute the derivative of sea surface temperature `tos_con` values
    along the NEMO parent domain `j` dimension:

    >>> nemo['gridT/tos_con'].derivative(dim="j")

    Compute the derivative of sea surface temperature `tos_con` values along a regional subset of a global,
    zonally periodic domain NEMO parent domain `i` dimension:

    >>> nemo['gridT/tos_con'].sel(i=slice(10, 80)).derivative(dim="i")

    Here, the zonal periodicity inherited from the NEMO model grid is automatically set to `False` since a subset
    of the global domain is no longer zonally periodic.

    Compute the discrete difference of absolute salinity `so_abs` values along the first NEMO
    nested child domain `k` dimension:

    >>> nemo['gridT/1_gridT/so_abs'].derivative(dim="k")

    See Also
    --------
    diff
    """
    # -- Validate Inputs -- #
    if not isinstance(dim, str):
        raise ValueError(
            "dim must be a string specifying dimension along which to calculate the gradient (e.g., 'i', 'i1', 'j', 'j1', 'k', 'k1')."
        )
    if dim not in self.dims:
        raise KeyError(
            f"dimension '{dim}' not found in {self.name or 'unnamed'} dimensions {self.dims}."
        )
    if not isinstance(fillna, bool):
        raise TypeError(
            "`fillna` must be specified as a boolean. Default is False."
        )


    # -- Calculate 1st-discrete derivative along dimension -- #
    # Determine NEMO model grid type and scale factors for derivative:
    new_grid_suffix, _, new_grid_weights = _NEMO_DIFF_MAP[(self._grid_suffix.upper(),
                                                           dim.replace(self._dom_suffix, "")
                                                           )]
    # Define path to derivative NEMO model grid:
    new_grid = f"{self._grid.replace(self._grid[-1], new_grid_suffix.upper())}"

    # Collect derivative grid scale factors (e.g., e1u, e3t etc.):
    try:
        weights = self._tree[new_grid][new_grid_weights]
    except KeyError as e:
        raise KeyError(
            f"NEMO model grid: '{new_grid}' does not contain grid scale factor '{new_grid_weights}' required to calculate derivatives along the {dim}-dimension."
        ) from e

    # Calculate 1st-finite difference along dimension:
    da = self.diff(dim=dim, fillna=fillna)

    # Calculate derivative (i.e., diff(var) / e{1/2/3}{t/u/v/w}):
    if dim in [self.k_name]:
        # Vertical derivative [k increasing downward]:
        result = - da.data / weights
    else:
        # Horizontal derivative [i/j increasing eastward/northward]:
        result = da.data / weights

    # -- Update DataArray properties & return NEMODataArray -- #
    result.name = f"d({self.name})/d{dim}"
    result = NEMODataArray(da=result, tree=self._tree, grid=new_grid)

    return result

diff

diff(dim, fillna=False)

Calculate the discrete difference of a variable along a given dimension (e.g., 'i', 'j', 'k') of a NEMO model grid.

Parameters:

Name	Type	Description	Default
dim	str	Dimension over which to calculate the finite difference (e.g., 'i', 'j', 'k').	required
fillna	bool	Fill NaN values in NEMODataArray with zeros prior to finite differencing. Default is False.	False

Returns:

Type	Description
NEMODataArray	Discrete difference of variable defined on a new NEMO model grid. For example, the discrete difference along the i-dimension of a scalar variable defined on a T-grid returns a NEMODataArray defined on the U-grid.

Examples:

Compute the difference of sea surface temperature tos_con values along the NEMO parent domain j dimension:

>>> nemo['gridT/tos_con'].diff(dim="j")

Compute the difference of sea surface temperature tos_con values along a regional subset of a global, zonally periodic domain NEMO parent domain i dimension:

>>> nemo['gridT/tos_con'].sel(i=slice(10, 80)).diff(dim="i")

Here, the zonal periodicity inherited from the NEMO model grid is automatically set to False since a subset of the global domain is no longer zonally periodic.

Compute the discrete difference of absolute salinity so_abs values along the first NEMO nested child domain k dimension:

>>> nemo['gridT/1_gridT/so_abs'].diff(dim="k")

See Also

derivative

Source code in nemo_cookbook/nemodataarray.py

def diff(
    self,
    dim: str,
    fillna: bool = False,
) -> Self:
    """
    Calculate the discrete difference of a variable along a given dimension
    (e.g., 'i', 'j', 'k') of a NEMO model grid.

    Parameters
    ----------
    dim : str
        Dimension over which to calculate the finite difference (e.g., 'i', 'j', 'k').
    fillna : bool, optional
        Fill NaN values in NEMODataArray with zeros prior to finite differencing. Default is False.
    Returns
    -------
    NEMODataArray
        Discrete difference of variable defined on a new NEMO model grid. For example, the
        discrete difference along the i-dimension of a scalar variable defined on a T-grid
        returns a NEMODataArray defined on the U-grid. 

    Examples
    --------
    Compute the difference of sea surface temperature `tos_con` values
    along the NEMO parent domain `j` dimension:

    >>> nemo['gridT/tos_con'].diff(dim="j")

    Compute the difference of sea surface temperature `tos_con` values along a regional subset of a global,
    zonally periodic domain NEMO parent domain `i` dimension:

    >>> nemo['gridT/tos_con'].sel(i=slice(10, 80)).diff(dim="i")

    Here, the zonal periodicity inherited from the NEMO model grid is automatically set to `False` since a subset
    of the global domain is no longer zonally periodic.

    Compute the discrete difference of absolute salinity `so_abs` values along the first NEMO
    nested child domain `k` dimension:

    >>> nemo['gridT/1_gridT/so_abs'].diff(dim="k")

    See Also
    --------
    derivative
    """
    # -- Validate Inputs -- #
    if not isinstance(dim, str):
        raise ValueError(
            "dim must be a string specifying dimension along which to calculate the gradient (e.g., 'i', 'i1', 'j', 'j1', 'k', 'k1')."
        )
    if dim not in self.dims:
        raise KeyError(
            f"dimension '{dim}' not found in {self.name or 'unnamed'} dimensions {self.dims}."
        )
    if not isinstance(fillna, bool):
        raise TypeError(
            "`fillna` must be specified as a boolean. Default is False."
        )

    # -- Get NEMO model grid properties -- #
    grid_paths = self._tree._get_grid_paths(dom=self._dom)

    # -- Calculate 1st-order discrete difference -- #
    if fillna:
        da = self.masked.data.fillna(value=0)
    else:
        da = self.masked.data

    if dim == self.k_name:
        match self._grid_suffix:
            case "w":
                # W-grid located at k = 0.5, 1.5 ... Nk-0.5
                result = da.diff(dim=self.k_name, n=1, label="lower")
                # Fill final T-point [k=Nk] -> NaN:
                result = result.pad({self.k_name: (0, 1)}, constant_values=np.nan)
                result.coords[self.k_name] = (result.coords[self.k_name] + 0.5).astype(int)
            case _:
                # T/U/V/F-grids located at k = 1, 2 ... Nk
                result = da.diff(dim=self.k_name, n=1, label="lower")
                # Fill initial W-point [k=0.5] -> NaN:
                result = result.pad({self.k_name: (1, 0)}, constant_values=np.nan)
                result.coords[self.k_name] = (result.coords[self.k_name].fillna(0) + 0.5)

    elif dim == self.i_name:
        match self._grid_suffix:
            case "t" | "v" | "w":
                # T/V/W-grids located at i = 1, 2 ... Ni
                if self.iperio:
                    result = da.roll({self.i_name: -1}) - da
                else:
                    result = da.shift({self.i_name: -1}) - da
                result.coords[self.i_name] = result.coords[self.i_name] + 0.5
            case "u" | "f":
                # U/F-grids located at i = 1.5, 2.5 ... Ni+0.5
                if self.iperio:
                    result = da - da.roll({self.i_name: 1})
                    result.coords[self.i_name] = result.coords[self.i_name] - 0.5
                else:
                    # Fill initial U/F-point [i=0.5] -> NaN or 0:
                    result = da.pad({self.i_name: (1, 0)}, constant_values=0)
                    result = result.diff(dim=self.i_name, n=1, label="lower")
                    result.coords[self.i_name] = (result.coords[self.i_name].fillna(0.5) + 0.5).astype(int)

    elif dim == self.j_name:
        match self._grid_suffix:
            case "t" | "u" | "w":
                # T/U/W-grids located at j = 1, 2 ... Nj
                result = da.shift({self.j_name: -1}) - da
                result.coords[self.j_name] = result.coords[self.j_name] + 0.5
            case "v" | "f":
                # V/F-grids located at j = 1.5, 2.5 ... Nj+0.5
                # Fill initial V/F-point [j=0.5] -> NaN or 0:
                result = da.pad({self.j_name: (1, 0)}, constant_values=0)
                result = result.diff(dim=self.j_name, n=1, label="lower")
                result.coords[self.j_name] = (result.coords[self.j_name].fillna(0.5) + 0.5).astype(int)
    else:
        raise ValueError(f"Invalid dimension {dim}. Dimension must be one of (i{self._dom}, j{self._dom}, k{self._dom}).")

    # -- Updating NEMO model grid coordinates -- #
    geo_coords = [coord for coord in da.coords if coord not in (self.t_name, self.k_name, self.i_name, self.j_name)]

    # Determine new NEMO model grid type & dims:
    new_ijk_dims = [dim for dim in result.dims if dim != 'time_counter']
    new_grid_suffix, new_hgrid_suffix, _ = _NEMO_DIFF_MAP[(self._grid_suffix.upper(), dim.replace(self._dom_suffix, ""))]
    if len(new_grid_suffix) > 1:
        new_vgrid_suffix = new_grid_suffix[1]
    else:
        new_vgrid_suffix = new_grid_suffix

    # Define new geographical coordinates (glam, gphi, depth) based on new grid type:
    # NOTE: Subsets of NEMODataTree geographical coordinates are supported implicitly since the
    # (i, j, k) coordinates of the grid dimensions (i, j, k) can be used to align incoming (i.e., _tree) coordinates.
    new_coords = {
        f"{self._dom_prefix}glam{new_hgrid_suffix}": self._tree[grid_paths[f"grid{new_hgrid_suffix.upper()}"]][f"{self._dom_prefix}glam{new_hgrid_suffix}"],
        f"{self._dom_prefix}gphi{new_hgrid_suffix}": self._tree[grid_paths[f"grid{new_hgrid_suffix.upper()}"]][f"{self._dom_prefix}gphi{new_hgrid_suffix}"],
    }
    if self.k_name in new_ijk_dims:
        new_coords[f"{self._dom_prefix}depth{new_vgrid_suffix}"] = self._tree[grid_paths[f"grid{new_vgrid_suffix.upper()}"]][f"{self._dom_prefix}depth{new_vgrid_suffix}"]

    # -- Update DataArray properties -- #
    result = result.drop_vars(geo_coords)
    result = result.assign_coords(new_coords)
    result.name = f"diff_{dim}({self.name})"

    # -- Apply land-sea mask & return NEMODataArray -- #
    new_grid = f"{self._grid.replace(self._grid[-1], new_grid_suffix.upper())}"
    result = NEMODataArray(da=result, tree=self._tree, grid=new_grid)
    result = result.masked

    return result

integral

integral(dims, cum_dims=None, dir=None, mask=None)

Integrate variable along one or more dimensions of a NEMO model grid.

Parameters:

Name	Type	Description	Default
dims	list	Dimensions over which to integrate (e.g., ['i', 'k']).	required
cum_dims	list	Dimensions over which to cumulatively integrate (e.g., ['k']). Specified dimensions must also be included in dims.	None
dir	str	Direction of cumulative integration. Options are '+1' (along increasing cum_dims) or '-1' (along decreasing cum_dims).	None
mask	DataArray | None	Boolean mask identifying NEMO model grid points to be included (1) or neglected (0) from integration.	None

Returns:

Type	Description
NEMODataArray	Variable integrated along specified dimensions of the NEMO model grid.

Examples:

Compute the integral of conservative temperature thetao_con along the vertical k dimension in the NEMO parent domain:

>>> nemo['gridT/thetao_con'].integral(dims=["k"])

Compute the vertical meridional overturning stream function from the meridional velocity vo (zonally integrated meridional velocity accumulated with increasing depth):

>>> nemo['gridV/vo'].integral(dims=["i", "k"],
...                           cum_dims=["k"],
...                           dir="+1",
...                           )

See Also

depth_integral

Source code in nemo_cookbook/nemodataarray.py

def integral(
    self,
    dims: list,
    cum_dims: list | None = None,
    dir: str | None = None,
    mask: xr.DataArray | None = None,
) -> Self:
    """
    Integrate variable along one or more dimensions of a NEMO model grid.

    Parameters
    ----------
    dims : list
        Dimensions over which to integrate (e.g., ['i', 'k']).
    cum_dims : list, optional
        Dimensions over which to cumulatively integrate (e.g., ['k']).
        Specified dimensions must also be included in `dims`.
    dir : str, optional
        Direction of cumulative integration. Options are '+1' (along
        increasing cum_dims) or '-1' (along decreasing cum_dims).
    mask: xr.DataArray, optional
        Boolean mask identifying NEMO model grid points to be included (1)
        or neglected (0) from integration.

    Returns
    -------
    NEMODataArray
        Variable integrated along specified dimensions of the NEMO model grid.


    Examples
    --------
    Compute the integral of conservative temperature `thetao_con` along the vertical
    `k` dimension in the NEMO parent domain:

    >>> nemo['gridT/thetao_con'].integral(dims=["k"])

    Compute the vertical meridional overturning stream function from the meridional
    velocity `vo` (zonally integrated meridional velocity accumulated with increasing
    depth):

    >>> nemo['gridV/vo'].integral(dims=["i", "k"],
    ...                           cum_dims=["k"],
    ...                           dir="+1",
    ...                           )

    See Also
    --------
    depth_integral
    """
    # -- Validate Input -- #
    if cum_dims is not None:
        for dim in cum_dims:
            if dim not in dims:
                raise ValueError(
                    f"cumulative integration dimension '{dim}' not included in `dims`."
                )
        if dir not in ["+1", "-1"]:
            raise ValueError(
                f"invalid direction of cumulative integration '{dir}'. Expected '+1' or '-1'."
            )
    if mask is not None:
        if not isinstance(mask, xr.DataArray):
            raise ValueError("mask must be an xarray.DataArray.")
        if any(dim not in self.dims for dim in mask.dims):
            raise ValueError(
                f"mask must have dimensions subset from {self.dims}."
            )

    # -- Collect variable, weights & mask -- #
    da = self.apply_mask(mask=mask).data
    weights = self._tree._get_weights(grid=self._grid, dims=dims)

    # -- Perform integration -- #
    if cum_dims is not None:
        sum_dims = [dim for dim in dims if dim not in cum_dims]
        if dir == "+1":
            # Cumulative integration along ordered dimension:
            result = (
                da.weighted(weights)
                .sum(dim=sum_dims, skipna=True)
                .cumsum(dim=cum_dims, skipna=True)
            )
        elif dir == "-1":
            # Cumulative integration along reversed dimension:
            result = (
                da.weighted(weights)
                .sum(dim=sum_dims, skipna=True)
                .reindex({dim: self.data[dim][::-1] for dim in cum_dims})
                .cumsum(dim=cum_dims, skipna=True)
            )
    else:
        result = da.weighted(weights).sum(dim=dims, skipna=True)

    # -- Apply land-sea mask & return NEMODataArray -- #
    result.name = f"integral_{''.join(dims)}({self.name})"
    result = self._wrap(result)

    return result

interp_to

interp_to(to)

Linearly interpolate variable to a neighbouring horizontal grid.

For flux variables defined at U/V-points, the specified variable is first weighted by grid cell face area prior to linear interpolation, and is then normalised by the area of the target grid cell face following interpolation.

Parameters:

Name	Type	Description	Default
to	str	Suffix of neighbouring horizontal NEMO model grid to linear interpolate variable to. Options are 'T', 'U', 'V', 'F'.	required

Returns:

Type	Description
NEMODataArray	Variable linearly interpolated onto a neighbouring horizontal grid.

Examples:

Linearly interpolate conservative temperature thetao_con defined on scalar T-points to neighbouring V-points in a NEMO model parent domain:

>>> nemo['gridT/thetao_con'].interp_to(to='V')

See Also

transform_vertical_grid

Source code in nemo_cookbook/nemodataarray.py

def interp_to(
    self,
    to: str,
) -> Self:
    """
    Linearly interpolate variable to a neighbouring horizontal grid.

    For flux variables defined at U/V-points, the specified variable
    is first weighted by grid cell face area prior to linear interpolation,
    and is then normalised by the area of the target grid cell face following
    interpolation.

    Parameters
    ----------
    to : str
        Suffix of neighbouring horizontal NEMO model grid to linear interpolate
        variable to. Options are 'T', 'U', 'V', 'F'.

    Returns
    -------
    NEMODataArray
        Variable linearly interpolated onto a neighbouring horizontal grid.

    Examples
    --------
    Linearly interpolate conservative temperature `thetao_con` defined on scalar T-points
    to neighbouring V-points in a NEMO model parent domain:

    >>> nemo['gridT/thetao_con'].interp_to(to='V')

    See Also
    --------
    transform_vertical_grid
    """
    # -- Validate input -- #
    if not isinstance(to, str):
        raise TypeError(f"'to' must be a string, got {type(to)}.")
    if to not in ["T", "U", "V", "F"]:
        raise ValueError(f"'to' must be one of ['T', 'U', 'V', 'F'], got {to}.")

    # -- Get NEMO model grid properties -- #
    ijk_names = self._tree._get_ijk_names(grid=self._grid)
    target_grid = f"{self._grid.replace(self._grid[-1], to)}"

    # -- Collect variable grid scale factors -- #
    if self._grid_suffix.upper() in ["U", "V"]:
        weight_dims = (
            [self.k_name, self.j_name] if self._grid_suffix.upper() == "U" else [self.k_name, self.i_name]
        )
        if f"{self._dom_prefix}depth{self._grid_suffix}" in self.coords:
            # 3-D variables - weight by grid cell face area:
            weights = self._tree._get_weights(grid=self._grid, dims=weight_dims, fillna=False)
            target_weights = self._tree._get_weights(
                grid=target_grid, dims=weight_dims, fillna=False
            )
        else:
            # 2-D variables - weight by grid cell width:
            weights = self._tree._get_weights(grid=self._grid, dims=weight_dims[1], fillna=False)
            target_weights = self._tree._get_weights(
                grid=target_grid, dims=weight_dims[1], fillna=False
            )
        da = self.masked.data * weights
    else:
        # Scalar variables:
        da = self.masked.data

    # -- Linearly interpolate variable -- #
    result = interpolate_grid(
        da=da,
        mask=self.mask,
        source_grid=self._grid_suffix.upper(),
        target_grid=to,
        iperio=self.iperio,
        ijk_names=ijk_names,
    )

    # -- Updating NEMO model grid coordinates -- #
    geo_coords = [coord for coord in da.coords if coord not in (self.t_name, self.k_name, self.i_name, self.j_name)]

    # Define new geographical coordinates (glam, gphi, depth) based on new grid type:
    # NOTE: Subsets of NEMODataTree geographical coordinates are supported implicitly since the
    # (i, j, k) coordinates of the grid dimensions (i, j, k) can be used to align incoming (i.e., _tree) coordinates.
    new_coords = {
        f"{self._dom_prefix}glam{to.lower()}": self._tree[target_grid][f"{self._dom_prefix}glam{to.lower()}"],
        f"{self._dom_prefix}gphi{to.lower()}": self._tree[target_grid][f"{self._dom_prefix}gphi{to.lower()}"],
    }
    if self.k_name in result.coords:
        new_coords[f"{self._dom_prefix}depth{to.lower()}"] = self._tree[target_grid][f"{self._dom_prefix}depth{to.lower()}"]

    # -- Update DataArray properties -- #
    result = result.drop_vars(geo_coords)
    result = result.assign_coords(new_coords)
    # Retain original variable name:
    result.name = da.name

    # Reorder dimensions (time_counter, [k], j, i):
    var_dims = [dim for dim in [self.t_name, self.k_name, self.j_name, self.i_name] if (dim is not None) and (dim in result.dims)]
    result = result.transpose(*var_dims)

    # Normalise by target grid cell weights for flux variables:
    if self._grid_suffix.upper() in ["U", "V"]:
        result = result / target_weights

    # -- Apply land-sea mask & return NEMODataArray -- #
    result = NEMODataArray(da=result, tree=self._tree, grid=target_grid)
    result = result.masked

    return result

masked_statistic

masked_statistic(lon_poly, lat_poly, statistic, dims, skipna=None)

Compute masked statistic of a variable defined on a NEMO model grid.

Parameters:

Name	Type	Description	Default
lon_poly	list | ndarray	Longitudes of closed polygon.	required
lat_poly	list | ndarray	Latitudes of closed polygon.	required
statistic	str	Name of the statistic to calculate (e.g., 'mean', 'weighted_mean' 'sum').	required
dims	list	Dimensions over which to apply statistic (e.g., ['i', 'j']).	required
skipna	bool | None	If True, skip missing values (as marked by NaN). By default, only skips missing values for float dtypes.	None

Returns:

Type	Description
NEMODataArray	Masked statistic of variable defined on a NEMO model grid.

Examples:

Compute the grid cell area-weighted mean sea surface temperature tos_con for a region enclosed in a polygon defined by lon_poly and lat_poly in a NEMO nested child domain:

>>> nemo["gridT/1_gridT/tos_con"].masked_statistic(lon_poly=lon_poly,
...                                                lat_poly=lat_poly,
...                                                statistic="weighted_mean",
...                                                dims=["i", "j"]
...                                                )

See Also

NEMODataTree.binned_statistic

Source code in nemo_cookbook/nemodataarray.py

def masked_statistic(
    self,
    lon_poly: list | np.ndarray,
    lat_poly: list | np.ndarray,
    statistic: str,
    dims: list,
    skipna: bool | None = None,
) -> Self:
    """
    Compute masked statistic of a variable defined on a NEMO model grid.

    Parameters
    ----------
    lon_poly : list | np.ndarray
        Longitudes of closed polygon.
    lat_poly : list | np.ndarray
        Latitudes of closed polygon.
    statistic : str
        Name of the statistic to calculate (e.g., 'mean', 'weighted_mean' 'sum').
    dims : list
        Dimensions over which to apply statistic (e.g., ['i', 'j']).
    skipna : bool | None
        If True, skip missing values (as marked by NaN). By default, only skips missing values for float dtypes.

    Returns
    -------
    NEMODataArray
        Masked statistic of variable defined on a NEMO model grid.

    Examples
    --------
    Compute the grid cell area-weighted mean sea surface temperature `tos_con` for a
    region enclosed in a polygon defined by `lon_poly` and `lat_poly` in a NEMO nested
    child domain:

    >>> nemo["gridT/1_gridT/tos_con"].masked_statistic(lon_poly=lon_poly,
    ...                                                lat_poly=lat_poly,
    ...                                                statistic="weighted_mean",
    ...                                                dims=["i", "j"]
    ...                                                )

    See Also
    --------
    NEMODataTree.binned_statistic
    """
    # -- Validate input -- #
    if not isinstance(statistic, str):
        raise TypeError("statistic must be specified as a string.")
    if not isinstance(dims, list):
        raise TypeError("dims must be specified as a list.")
    if skipna is not None:
        if not isinstance(skipna, bool):
            raise TypeError("skipna must be specified as a boolean or None.")

    # -- Create polygon mask using coordinates -- #
    mask_poly = self._tree.mask_with_polygon(
        lon_poly=lon_poly, lat_poly=lat_poly, grid=self._grid
    )

    # -- Apply masks & calculate statistic -- #
    da = self.apply_mask(mask=mask_poly).data

    match statistic:
        case "mean":
            result = da.mean(dim=dims, skipna=skipna)

        case "weighted_mean":
            weight_dims = [dim.replace(self._dom_suffix, "") for dim in dims]
            weights = self._tree._get_weights(grid=self._grid, dims=weight_dims)
            result = da.weighted(weights).mean(dim=dims, skipna=skipna)

        case "min":
            result = da.min(dim=dims, skipna=skipna)

        case "max":
            result = da.max(dim=dims, skipna=skipna)

        case "sum":
            result = da.sum(dim=dims, skipna=skipna)

        case _:
            raise ValueError(
                f"Unsupported statistic '{statistic}'. Supported statistics are: 'mean', 'weighted_mean', 'min', 'max', 'sum'."
            )

    # -- Update DataArray properties -- #
    result.name = f"masked_{statistic}({self.name})"
    result = self._wrap(result)

    return result

sel_like

sel_like(other)

Return a new NEMODataArray whose data is given by matching the dimension index labels present in another NEMODataArray or xarray.DataArray.

Parameters:

Name	Type	Description	Default
other	NEMODataArray | DataArray	NEMODataArray or xarray.DataArray used to select dimension index labels.	required

Returns:

Type	Description
NEMODataArray	A new NEMODataArray with data selected according to dimension index labels of input object.

Examples:

Indexing conservative temperature thetao_con defined on scalar T-points in a NEMO model parent domain to match a subset of the absolute salinity so_abs:

>>> nda = nemo['gridT/so_abs'].sel(time_counter=slice("2020-01", "2024-01"), k=1)

>>> nemo['gridT/thetao_con'].sel_like(nda)

See Also

sel

Source code in nemo_cookbook/nemodataarray.py

def sel_like(
    self,
    other: Self | xr.DataArray,
) -> Self:
    """
    Return a new NEMODataArray whose data is given by matching the dimension
    index labels present in another NEMODataArray or xarray.DataArray.

    Parameters
    ----------
    other : NEMODataArray | xr.DataArray
        NEMODataArray or xarray.DataArray used to select dimension index labels.

    Returns
    -------
    NEMODataArray
         A new NEMODataArray with data selected according to dimension index labels
         of input object.

    Examples
    --------
    Indexing conservative temperature `thetao_con` defined on scalar T-points in a NEMO
    model parent domain to match a subset of the absolute salinity `so_abs`:

    >>> nda = nemo['gridT/so_abs'].sel(time_counter=slice("2020-01", "2024-01"), k=1)

    >>> nemo['gridT/thetao_con'].sel_like(nda)

    See Also
    --------
    sel
    """
    # -- Validate Inputs -- #
    if not (isinstance(other, NEMODataArray) or isinstance(other, xr.DataArray)):
        raise TypeError("other must be specified as a NEMODataArray or xarray.DataArray.")

    if isinstance(other, NEMODataArray):
        other = other.data

    # -- Index data using dimension labels of other object & return NEMODataArray -- #
    coord_dims = [self.t_name, self.k_name, self.j_name, self.i_name]
    d_dims = {}
    for dim in coord_dims:
        # Select only dimensions present in both objects:
        if (dim in other.coords) and (dim in self.data.coords):
            # Select only dimensions whose index labels do not already match:
            if other.coords[dim].size != self.data.coords[dim].size:
                d_dims[dim] = other.coords[dim].data

    # Indexing only required if dimensions modified:
    if len(d_dims) > 0:
        result = self.data.sel(d_dims)
        return self._wrap(result)
    # Otherwise return original NEMODataArray:
    else:
        return self

to_xesmf

to_xesmf(mask=True)

Create an xESMF-compatible dataset from a variable stored in a NEMODataArray with the following properties:

• lon (longitudes of grid cell centers)
• lat (latitudes of grid cell centers)
• lon_b (longitudes of grid cell boundaries)
• lat_b (latitudes of grid cell boundaries)
• mask (boolean land-sea mask identifying valid grid cells) [ Optional ]

Only NEMODataArrays defined on NEMO model T-grids are currently supported for transformation to xESMF-compatible datasets.

Users can use interp_to(to='T') to linearly interpolate NEMO variables defined on neighbouring U/V/F-grids to T-grids prior to creating an xESMF-compatible dataset.

Parameters:

Name	Type	Description	Default
mask	bool	Whether to include the land-sea mask in the resulting dataset. Default is True.	True

Returns:

Type	Description
Dataset	Dataset containing variable and associated geographical coordinates formatted for use with xESMF.

Examples:

Create an xESMF-compatible dataset from the conservative temperature variable thetao_con defined on a NEMO model parent domain T-grid:

>>> nemo['gridT/thetao_con'].to_xesmf(mask=True)

See Also

interp_to

Source code in nemo_cookbook/nemodataarray.py

def to_xesmf(
        self, mask: bool = True
    ) -> xr.Dataset:
    """
    Create an xESMF-compatible dataset from a variable stored in a NEMODataArray
    with the following properties:

    - lon (longitudes of grid cell centers)
    - lat (latitudes of grid cell centers)
    - lon_b (longitudes of grid cell boundaries)
    - lat_b (latitudes of grid cell boundaries)
    - mask (boolean land-sea mask identifying valid grid cells) [ Optional ]

    Only NEMODataArrays defined on NEMO model T-grids are currently supported
    for transformation to xESMF-compatible datasets.

    Users can use interp_to(to='T') to linearly interpolate NEMO variables defined
    on neighbouring U/V/F-grids to T-grids prior to creating an xESMF-compatible
    dataset.

    Parameters
    ----------
    mask : bool, optional
        Whether to include the land-sea mask in the resulting dataset.
        Default is True.

    Returns
    -------
    xr.Dataset
        Dataset containing variable and associated geographical coordinates
        formatted for use with xESMF.

    Examples
    --------
    Create an xESMF-compatible dataset from the conservative temperature variable
    `thetao_con` defined on a NEMO model parent domain T-grid:

    >>> nemo['gridT/thetao_con'].to_xesmf(mask=True)

    See Also
    --------
    interp_to
    """
    # -- Validate Input -- #
    if not isinstance(mask, bool):
        raise ValueError("mask must be specified as a boolean. Default is True.")
    if self._grid_suffix.upper() != "T":
        raise ValueError("`to_xesmf()` accessor only supports variables defined on a NEMO model T-grid." \
        "Use `interp_to(to='T')` to linearly interpolate variables onto the T-grid first."
        )

    # -- Create xESMF-compatible dataset -- #
    ds = create_xesmf_dataset(nda=self, mask=mask)

    return ds

transform_vertical_grid

transform_vertical_grid(e3_new)

Transform variable defined on a NEMO model grid to a new vertical grid using conservative interpolation.

Parameters:

Name	Type	Description	Default
e3_new	DataArray	Grid cell thicknesses of the new vertical grid. Must be a 1-dimensional xarray.DataArray with dimension 'k_new'.	required

Returns:

Type	Description
Dataset	Variable defined at the centre of each vertical grid cell on the new grid, and vertical grid cell thicknesses adjusted for model bathymetry.

Examples:

Transform the conservative temperature variable thetao_con defined in a NEMO model parent domain from it's native 75 unevenly-spaced z-levels to regularly spaced z-levels at 200 m intervals:

>>> e3t_target = xr.DataArray(np.repeat(200.0, 30), dims=['k_new'])

>>> nemo['gridT/thetao_con'].transform_vertical_grid(e3_new=e3t_target)

See Also

transform_to

Source code in nemo_cookbook/nemodataarray.py

def transform_vertical_grid(
    self, e3_new: xr.DataArray
) -> xr.Dataset:
    """
    Transform variable defined on a NEMO model grid to a new vertical grid using conservative interpolation.

    Parameters
    ----------
    e3_new : xarray.DataArray
        Grid cell thicknesses of the new vertical grid.
        Must be a 1-dimensional xarray.DataArray with
        dimension 'k_new'.

    Returns
    -------
    xr.Dataset
        Variable defined at the centre of each vertical
        grid cell on the new grid, and vertical grid cell
        thicknesses adjusted for model bathymetry.

    Examples
    --------
    Transform the conservative temperature variable `thetao_con` defined in a
    NEMO model parent domain from it's native 75 unevenly-spaced z-levels to
    regularly spaced z-levels at 200 m intervals:

    >>> e3t_target = xr.DataArray(np.repeat(200.0, 30), dims=['k_new'])

    >>> nemo['gridT/thetao_con'].transform_vertical_grid(e3_new=e3t_target)

    See Also
    --------
    transform_to
    """
    # -- Validate Input -- #
    if e3_new.dims != ("k_new",) or (e3_new.ndim != 1):
        raise ValueError(
            "e3_new must be a 1-dimensional xarray.DataArray with dimension 'k_new'."
        )

    # -- Define input variables -- #
    var_in = self.masked.data
    e3_in = self.metrics["e3"].masked.data

    # Ensure total depth of new vertical grid >= total depth of NEMO model vertical grid:
    depth_max = self._tree[self._grid][f"{self._dom_prefix}depth{self._grid_suffix}"].max(self.k_name)
    if e3_new.sum(dim="k_new") < depth_max:
        raise ValueError(
            f"e3_new must sum to at least the maximum depth ({depth_max.item()} m) of the original vertical grid."
        )

    # -- Transform variable to target vertical grid -- #
    var_out, e3_out = xr.apply_ufunc(
        transform_vertical_coords,
        e3_in,
        var_in,
        e3_new.astype(e3_in.dtype),
        input_core_dims=[[self.k_name], [self.k_name], ["k_new"]],
        output_core_dims=[["k_new"], ["k_new"]],
        dask="parallelized",
        output_dtypes=[var_in.dtype, e3_in.dtype],
        dask_gufunc_kwargs={"output_sizes": {"k_new": e3_new.sizes["k_new"]}},
    )

    # -- Construct transformed variable Dataset -- #
    var_dims = [dim for dim in [self.t_name, "k_new", self.j_name, self.i_name] if (dim is not None) and (dim in var_out.dims)]
    var_out = var_out.transpose(*var_dims)

    e3_dims = [dim for dim in [self.t_name, "k_new", self.j_name, self.i_name] if (dim is not None) and (dim in e3_out.dims)]
    e3_out = e3_out.transpose(*e3_dims)

    result = xr.Dataset(
        data_vars={self.name: var_out, f"e3{self._grid_suffix}_new": e3_out},
        coords={
            f"depth{self._grid_suffix}_new": ("k_new", e3_new.cumsum(dim="k_new").data)
        },
    )

    return result

weighted_mean

weighted_mean(dims, mask=None, skipna=None)

Calculate grid-aware weighted mean of a variable defined on a NEMO model grid.

Parameters:

Name	Type	Description	Default
dims	list	Dimensions over which to apply weighted mean (e.g., ['i', 'j']).	required
mask	DataArray | None	Boolean mask identifying NEMO model grid points to be included (1) or neglected (0) from integration.	None
skipna	bool | None	If True, skip missing values (as marked by NaN). By default, only skips missing values for float dtypes.	None

Returns:

Type	Description
NEMODataArray	Grid-aware weighted mean of variable defined on a NEMO model grid.

Examples:

Area-weighted mean of the sea surface temperature tos_con defined on scalar T-points in a NEMO model nested child domain:

>>> nemo["gridT/1_gridT/tos_con"].weighted_mean(dims=["i", "j"], skipna=True)

See Also

masked_statistic

Source code in nemo_cookbook/nemodataarray.py

def weighted_mean(
    self,
    dims : list,
    mask: xr.DataArray | None = None,
    skipna : bool | None = None
) -> Self:
    """
    Calculate grid-aware weighted mean of a variable defined on a NEMO model grid.

    Parameters
    ----------
    dims : list
        Dimensions over which to apply weighted mean (e.g., ['i', 'j']).
    mask: xr.DataArray, optional
        Boolean mask identifying NEMO model grid points to be included (1)
        or neglected (0) from integration.
    skipna : bool | None
        If True, skip missing values (as marked by NaN).
        By default, only skips missing values for float dtypes.

    Returns
    -------
    NEMODataArray
        Grid-aware weighted mean of variable defined on a NEMO model grid.

    Examples
    --------
    Area-weighted mean of the sea surface temperature `tos_con` defined on scalar T-points
    in a NEMO model nested child domain:

    >>> nemo["gridT/1_gridT/tos_con"].weighted_mean(dims=["i", "j"], skipna=True)

    See Also
    --------
    masked_statistic
    """
    # -- Validate Input -- #
    if not isinstance(dims, list):
        raise TypeError("dims must be specified as a list.")
    if mask is not None:
        if not isinstance(mask, xr.DataArray):
            raise TypeError("mask must be an xarray.DataArray.")
        if any(dim not in self.dims for dim in mask.dims):
            raise ValueError(
                f"mask must have dimensions subset from {self.dims}."
            )
    if skipna is not None:
        if not isinstance(skipna, bool):
            raise TypeError("skipna must be specified as a boolean or None.")

    # -- Calculate weighted mean & return NEMODataArray -- #
    da = self.apply_mask(mask=mask).data
    weight_dims = [dim.replace(self._dom_suffix, "") for dim in dims]
    weights = self._tree._get_weights(grid=self._grid, dims=weight_dims)
    result = da.weighted(weights).mean(dim=dims, skipna=skipna)
    result.name = f"wmean_{'_'.join(dims)}({self.name})"

    return self._wrap(result)