A Quickstart Guide to Common Operations with NEMO Cookbook
In this section, we describe some of the most common NEMODataTree and NEMODataArray operations in a concise how-to guide (inspired by the excellent documentation of Icechunk).
For more detailed documentation on NEMODataTree and NEMODataArray structures, users should visit the API Reference.
Create a NEMODataTree from Local Files
We can create a NEMODataTree from a dictionary of paths to local netCDF files using the .from_paths() constructor:
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",
},
}
NEMODataTree.from_paths(paths, iperio=True, nftype="T")
In the example above, we consider only a global parent domain, which is zonally periodic (iperio=True) and north-folding on T grid points (nftype="T"). Note, that we are only required to specify paths for one or more NEMO model grid types (e.g., *_gridT.nc).
For NEMO models using a linear free surface approximation (i.e., vertical scale factors are time-independent), we should also specify linssh=True to read these directly from the domain_cfg file included in the paths dictionary.
Additionally, for NEMO models configured with more complex vertical coordinates (e.g., MEs or sigma-coordinates), such that vertical reference variables (e.g., depth) vary spatially, we should also specify vco="3d" to include all vertical reference variables as 3-dimensional arrays analogously to using key_vco_3d within NEMO itself. By default, a NEMODataTree is constructed using 1-dimensional vertical reference variables.
Create a NEMODataTree from xarray.Datasets
Alternatively, we can create a NEMODataTree from a dictionary of single or multi-file xarray.Datasets. This is particularly valuable when working with remote NEMO model data or Coupled Model Intercomparison Project (CMIP) outputs which require us to reformat coordinate dimensions (see Example NEMODataTrees).
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, read_masks=True)
This example would be applicable to the outputs of a regional NEMO model configuration which is neither zonally periodic nor north-folding (by default, iperio=False & nftype=None).
If all land-sea masks expected by NEMODataTree are included in ds_domain, we can also specify read_mask=True to read rather than calculate land-sea masks when constructing our NEMODataTree. This is recommended for larger NEMO model domains (e.g., eORCA12, ORCA36 etc.)
The following core dimensions are required to construct a new NEMODataTree from xarray.Datasets:
-
domain: (
nav_lev,y,x) -
grid{T/U/V/W}: (
time_counter,depth{p},y,x), where p is the grid point type.
Merging NEMODataTrees
To combine the variables in multiple NEMODataTrees into one larger NEMODataTree, we can use the merge() method, which takes a list of NEMODataTrees to merge:
We can also pass additional keyword arguments to xarray.merge() alongside our list of NEMODataTree objects. In the example above, we pass compat="no_conflicts" to compare variables of the same name such that only values which are not null in both NEMODataTrees must be equal.
Concatenating NEMODataTrees
To combine NEMODataTrees along an existing or new dimension into a larger NEMODataTree, we can use the concat() method, which takes a list of NEMODataTrees and a dimension name or xarray.DataArray defining the dimension along which to concatenate variables.
Note, the merge() and concat() methods are simply convienence wrappers to xarray.merge() and xarray.concat() functions, which handle the process of assembling the resulting xarray.DataTree object back into a compliant NEMODataTree.
Access NEMO Variables
To access a variable stored in a given grid node of a NEMODataTree as an xarray.DataArray, we can use the following syntax:
However, if we want to access the chosen variable as a grid-aware NEMODataArray (recommended), we can instead provide the direct path to the variable:
Adding a NEMO Variable to a NEMODataTree
To add a new variable stored as a NEMODataArray to the grid node of an existing NEMODataTree, we can use the following syntax:
Note, this simply adds the underlying xarray.DataArray to the NEMO model T-grid node xarray.Dataset in our NEMODataTree with the name my_var.
Adding a New Grid Node to a NEMODataTree
To create a new grid node from an xarray.Dataset containing NEMO model grid variables in an existing NEMODataTree, we can use the following syntax:
Before the new grid node is added to an existing NEMODataTree, the xarray.Dataset will be validated to ensure it contains:
-
NEMO grid dimension coordinates (
i,j). -
Longitude and Latitude coordinates named
gphi{x}andglam{x}. -
Depth coordinate named
depth{x}if NEMO grid dimensionkexists.
where x is final character of the specified grid node. In this example, gphip, glamp and depthp would be expected since we are introducing a new grid node named "gridP".
Access NEMO Variable Properties
Each variable stored as a NEMODataArray has a collection of useful properties to support grid-aware calculations:
-
Path to NEMO model grid node where variable is stored:
-
Type of NEMO model grid where variable is defined (i.e., t, u, v, etc.):
-
Dictionary of variable NEMO grid scale factors (e.g., e1t, e2t, e3t):
-
Variable land-sea mask:
Apply Land-Sea Mask to NEMO Variable
To mask a given NEMODataArray variable with its associated land-sea mask, we can use the .masked property:
In the example above, the 4-dimensional conservative temperature variable thetao_con is masked using the 3-dimensional tmask before being returned as a NEMODataArray.
Apply Custom Mask to NEMO Variable
To apply a custom mask to a given NEMODataArray variable, we can use the .apply_mask() method:
Here, the 4-dimensional absolute salinity variable so_abs is masked using the my_mask boolean mask before being returned as a NEMODataArray.
To drop the absolute salinity values where my_mask is False, we can optionally use .apply_mask(mask=my_mask, drop=True).
Index a NEMODataArray to Match the Dimension Labels of Another (NEMO)DataArray
In addition to the more familiar .sel() and .isel() label based selection methods, the .sel_like() method can be used to index a NEMODataArray according to the dimension index labels of another NEMODataArray or xarray.DataArray as follows:
nda = nemo["gridT/so_abs"].sel(time_counter=slice('2000-01', '2025-01'), k=1)
nemo["gridT/thetao_con"].sel_like(nda)
Calculate Grid Cell Areas
To calculate the area of a model grid cell face, we can use the .cell_area() method.
For example, to compute the horizontal area of cells centered on T grid points in the parent domain:
Importantly, the dim argument represents the dimensional orthogonal to the grid cell area to be computed. For T grid points, this results in the following grid cell areas:
dim |
Grid Cell Area |
|---|---|
i |
e2t * e3t |
j |
e1t * e3t |
k |
e1t * e2t |
Calculate Grid Cell Volumes
To calculate the volume of model grid cells, we can use the .cell_volume() method.
For example, to compute the volume of each grid cell centered on a V grid point in the model parent domain:
Plotting a NEMODataArray using Geographical Coordinates
To plot a 2-dimensional slice of a NEMODataArray as a Cartopy GeoQuadMesh using its longitude & latitude coordinates (i.e., glam{t/u/v/w}(j, i) & gphi{t/u/v/w}(j, i)), we can use the .geoplot() method.
For example, to create a geographical plot of sea surface temperature:
Similarly to xarray.DataArray.plot, the .geoplot() method allows for significant customisation of geographical plots.
For example, to create a geographical plot of sea ice concentration using a North Polar Stereographic projection with a bounding box (-180°E, 180°E, 50°N, 90°N):
(nemo["gridT/siconc"]
.isel(time_counter=0)
.geoplot(projection=ccrs.NorthPolarStereo(),
vmin=0, vmax=1,
extent=(-180, 180, 50, 90),
cmap='Blues_r',
clabel_kwargs={'label': 'Sea Ice Concentration [fraction of unity]'}
)
)
Indexing with Geographical Coordinates
To subset variables of a given model grid using their longitude & latitude coordinates (i.e., glam{t/u/v/w}(j, i) & gphi{t/u/v/w}(j, i)), we can add these geographical variables as indexes using the .add_geoindex() method.
For example, to enable geographical indexing of the parent T grid points & select the values of this dataset nearest to (-30°E, 60°N):
nemo_geo = nemo.add_geoindex(grid="gridT")
nemo_geo["gridT"].dataset.sel(gphit=60, glamt=-30, method="nearest")
Clip a NEMODataArray
To clip a NEMODataArray using a geographical bounding box defined by a tuple of the form (lon_min, lon_max, lat_min, lat_max), we can use the .clip() method.
For example, to clip the sea surface temperature variable tos_con defined on T-points in the bounding box (-40°E, 10°E, 35°N, 60°N):
Clip a NEMO Model Grid
To clip a given model grid using a geographical bounding box defined by a tuple of the form (lon_min, lon_max, lat_min, lat_max), we can use the .clip_grid() method.
For example, to clip the parent T-grid in the bounding box (-80°E, 0°E, 40°N, 80°N):
Clip a NEMO Model Domain
To clip all of the model grids of a given NEMO model domain:
where dom is the prefix of the chosen NEMO model domain. Note dom="." for the parent domain.
Calculate Horizontal Derivatives
To calculate the derivative of a scalar variable var along one of the horizontal dimensions (e.g., i, j) of a given NEMO model grid, we can use the .derivative() method.
For example, to compute the 'meridional' derivative of sea surface temperature tos_con along the NEMO model parent domain j dimension:
Alternatively, to compute the derivative of sea surface temperature tos_con along a regional subset of a global, zonally periodic domain NEMO model parent domain i dimension:
Note, using iperio=False overrides the zonal periodicity inherited from the NEMO model grid since the selected subset of the global domain is no longer zonally periodic.
Both the .diff() and .derivative() methods provide flexibility on how to handle NaN values using the fillna argument. By default, fillna=False, meaning that NaN values are not filled with zeros prior to performing finite difference operations. However, in the case of velocity variables this may be inappropriate in the vicinity of coastlines and users may prefer to use fillna=True to impose zero-magnitude velocity components along land-sea boundaries prior to calculating derivatives.
Calculate Vertical Derivative
To calculate the vertical derivative of a scalar variable var along the k dimension of a given NEMO model grid, we can also use the .derivative() method.
For example, to compute the vertical gradient of absolute salinity in our first NEMO model nested child domain:
Calculate Divergence
To calculate the horizontal divergence from the i and j components of a vector field, we can use the .divergence() method.
For example, to compute the horizontal divergence from the seawater velocity field in the NEMO model parent domain:
where uv_vars is a list specifying the names of i and j vector components, respectively.
Calculate Curl
To calculate the vertical k component of the curl of a horizontal vector field, we can use the .curl() method.
For example, to compute the vertical component of the curl of the seawater velocity field in the second NEMO nested child domain:
where, as in the case of .divergence(), the uv_vars argument expects a list of the i and j components of the vector field, respectively.
Calculate Integrals
To integrate a variable along one or more dimensions of a given NEMO model grid, we can use the .integral() method.
For example, to compute the integral of conservative temperature thetao_con along the vertical k dimension in the NEMO model parent domain:
which will return an NEMODataArray with one less dimension than thetao_con, in this case k since we have integrated vertically.
Calculate Cumulative Integrals
We can also use the .integral() method to calculate cumulative integrals along one or more dimensions of a given NEMO model grid.
For example, to calculate the vertical meridional overturning stream function from the meridional velocity vo (zonally integrated meridional velocity accumulated with increasing depth):
dims is a list of grid dimension names along which integration will be performed, and cum_dims specifies which of the dimensions in dims should be cumulatively integrated.
The dir argument is used to define the direction of cumulative integration, where dir = "+1" means accumulating along the chosen dimension, such that grid indices are increasing. Conversely, dir = "-1" means that cumulative integration is performed after reversing the chosen dimension, such that grid dimensions are decreasing.
Note, we can also pass the mask argument to .integral() to mask the variable var prior to performing the integration.
Calculate Depth Integrals
To integrate a variable of a given NEMO model grid in depth coordinates between two limits, we can use the .depth_integral() method.
For example, to compute the vertical integral of conservative temperature thetao_con in the upper 100 m in the NEMO model parent domain:
where limits is 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.
Calculate Weighted Average
To calculate a grid-aware weighted mean of a variable defined on a NEMO model grid, we can use the .weighted_mean() method.
For example, to compute the grid cell area-weighted mean sea surface temperature tos_con in a NEMO model nested child domain:
where dims represents the dimensions of the NEMO model grid to average over. In this example, dims=["i", "j"] is equivalent to computing the mean of variable tos_con using the horizontal cell area of T grid points (i.e., e1t * e2t) as weights.
Create Regional Masks using Polygons
To define a regional mask using the geographical coordinates of a closed polygon, we can use the .mask_with_polygon() method:
where lon_poly and lat_poly are lists or ndarrays containing the longitude and latitude coordinates defining the closed polygon.
Calculate Statistics for a Region Masked using a Polygon
To calculate an aggregated statistic from only the model grid cells contained inside a geographical polygon, we can use the .masked_statistic() method.
For example, to 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 model nested child domain:
nemo["gridT/1_gridT/tos_con"].masked_statistic(lon_poly,
lat_poly,
statistic="weighted_mean",
dims=["i", "j"]
)
where dims represent the dimensions of the NEMO model grid used for aggregation. In this example, combining statistic="weighted_mean" and dims=["i", "j"] is equivalent to computing the mean of variable tos_con using the horizontal cell area of T grid points (i.e., e1t * e2t) as weights.
Calculate Binned Statistics
To calculate aggregated statistics of a variable binned according to the values of one or more variables, we can use the .binned_statistic() method.
This is a generalization of a histogram function, enabling the computation of the sum, mean, median, or other statistic of the values assigned to each bin.
For example, to 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",
)
where vars is a list of the names of variables to be binned using the bin edges passed to bins, and values is the name of the variable over which the statistic will be performed once values have been grouped into each bin.
We can use keep_dims to specify the dimensions of the xarray.DataArray named values to retain. In the example above, using keep_dims="time_counter" will return the average depths of water in each potential density bin for each time-slice of available NEMO model output.
Interpolate Variable to a Neighbouring Horizontal Grid
To linearly interpolate a variable defined on a given NEMO horizontal grid to a neighbouring grid, we can use the .interp_to() method.
For example, to interpolate conservative temperature thetao_con defined on scalar T-points to neighbouring V-points in a NEMO model parent domain:
We can also interpolate variables defined on U- and V-points to either scalar or vector grid points. Unlike interpolating scalar variables defined on T-points, this is achieved by linearly interpolating the grid cell face area-weighted flux onto the target grid, before then normalising by the grid cell face area defined on the target horizontal grid.
For example, to interpolate the zonal wind stress defined on U-points to neighbouring V-points in a NEMO model parent domain and store this in the V-grid node of our NEMODataTree:
Transform a Vertical Grid
To transform a variable defined on a given NEMO model vertical grid to a new vertical grid using conservative interpolation, we can use the .transform_vertical_grid() method.
For example, if we wanted to 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)
where e3_new represents the time-invariant vertical grid cell thicknesses defing the vertical grid onto which the variable var will be conservatively interpolated.
There are some important points to remember when transforming variables onto new vertical grids with NEMODataTree:
-
New vertical grid cell thicknesses
e3_newmust sum to at least the maximum depth of the original vertical grid cell thicknesses (e.g., e3t). -
Currently,
e3_newmust be a 1-dimensionalxarray.DataArraywith dimension 'k_new'. -
The output
xarray.Datasetwill contain multi-dimensionalxarray.DataArraysfor both the vertically remapped variablevar(time_counter, k_new, j, i)and the vertical grid cell thicknessese3t_new(time_counter, k_new, j, i)(updated to explicitly account for partial grid cells above the seafloor).
Extract a Zonal Section
To extract an approximately zonal hydrographic section, we can use the .extract_zonal_section() method.
For example, to extract the RAPID-MOCHA array located at 26.5°N:
nemo.extract_zonal_section(lat=26.5,
lon_min=-82.0,
lon_max=-10.0,
u_vars=["tauuo"],
scalar_vars=['thetao_con', 'so_abs'],
dom='.'
)
Here, we have also provided the names of the zonal vector variables (tauuo) and any scalar variables (thetao_con and so_abs) to be linearly interpolated onto the V-grid.
The resulting xarray.Dataset will be defined on the NEMO model V-grid at the constant j-coordinate whose average latitude is closest to 26.5°N and is structured analogously to a NEMO model grid node within a NEMODataTree, including:
- Inherited NEMO Grid Dimension:
k - New Along-Boundary Grid Dimension:
bdy, such thati_bdy(bdy),j_bdy(bdy) - Geographical Coordinates:
glamb(longitude),gphib(latitude). - Horizontal Grid Scale Factors:
e1b,e2b - Vertical Grid Scale Factors:
e3b - Land-Sea Masks:
bmask,bmaskutil
Extract a Meridional Section
To extract an approximately meridional hydrographic section, we can use the .extract_meridional_section() method.
For example, to extract a section across Drake Passage located at -67°E in the Southern Ocean:
nemo.extract_meridional_section(lon=-67,
lat_min=-62,
lat_max=-55,
v_vars=["tauvo"],
scalar_vars=['thetao_con', 'so_abs'],
dom='.'
)
Here, we have also provided the names of the meridional vector variable (tauvo) and any scalar variables (thetao_con and so_abs) to be linearly interpolated onto the U-grid.
The resulting xarray.Dataset will be defined on the NEMO model U-grid at the constant i-coordinate whose average longitude is closest to -67°E and is structured analogously to a NEMO model grid node within a NEMODataTree as shown above in the extract_zonal_section() example above.
Extract a Piecewise Section
To extract a piecewise hydrographic section, we can use the .extract_section() method.
For example, to extract the Overturning in the Subpolar North Atlantic (OSNAP) array as a continous hydrographic section comprised of connected U and V-grid cell faces in the NEMO model parent domain:
nemo.extract_section(lon_section=lon_osnap,
lat_section=lat_osnap,
uv_vars=["uo", "vo"],
vars=["thetao_con"],
dom=".",
)
where lon_osnap and lat_osnap are the longitude and latitude coordinates defining the OSNAP section, respectively.
Here, we have also provided the names of the zonal (uo) and meridional (vo) sea water velocity variables, and the scalar variable thetao_con to be linearly interpolated onto the U and V-grid points comprising the section.
The resulting xarray.Dataset is structured as outlined above, but also includes a variable flux_type which defines the type of vector point the grid cell faces is and flux_dir which specifies the sign of each normal velocity component.
The velocity variable included in the output xarray.Dataset therefore represents a combination of zonal (U) and meridional (V) velocities to which the appropriate flux_dir has already been applied.
Extract the Variables along the Boundary of a Masked Region
To extract scalar and vector variables along the boundary of a masked region, we can use the .extract_mask_boundary() method.
For example, to extract the velocities and sea water temperature along the boundary of the Barents Sea, we first define a regional mask using the geographical coordinates (lon_BSea, lat_BSea) of the closed polygon enclosing the Barents Sea.
Next, we pass our regional mask to the .extract_mask_boundary() method:
mask = nemo.mask_with_polygon(grid='gridT',
lon_poly=lon_BSea,
lat_poly=lat_BSea
)
nemo.extract_mask_boundary(mask=mask,
uv_vars=["uo", "vo"],
vars=["thetao_con"],
dom="."
)
Here, we define the names of the zonal and meridional velocity components as a list uv_vars and the names of any scalar variables defined on T-grid points (vars) that we would like to extract along the boundary.
The resulting xarray.Dataset is structured analogously to that returned from .extract_section(), which is outlined above.
Export Variable to an xESMF Dataset
To access a variable stored in a NEMODataArray as an xESMF compatible xarray.Dataset to perform regridding, we can use the .to_xesmf() method.
For example, to access the sea surface temperature variable tos_con defined on scalar T-points in a NEMO model parent domain:
The resulting xarray.Dataset includes the following coordinate variables as expected by xESMF:
lon- longitude of T-grid cell centers .lat- latitude of T-grid cell centers.lon_b- longitude of T-grid cell corners (lower left F-grid points).lat_b- longitude of T-grid cell corners (lower left F-grid points).
By default, to_xesmf() also returns a 2-dimensional land-sea mask to support masked regridding to ensure source land grid cells do not contribute to target ocean grid cells. Following ESMF convention, the mask varible is defined ocean = 1 and land = 0 consistent with NEMO land-sea masks.
Note, the to_xesmf() accessor only supports variables defined on scalar T-points currently.
For more information on regridding NEMO outputs using xESMF, visit the Regridding using xESMF recipe on the Recipes page.
Create a NEMODataArray using an xESMF Dataset
To create a NEMODataArray from a variable stored in an xESMF compatible xarray.Dataset following regridding, we can use the .from_xesmf() method.
For example, to create a NEMODataArray called tos_obs from sea surface temperature observations regridded to scalar T-points in a NEMO model parent domain using xESMF:
The resulting NEMODataArray can now be used for grid-aware computation and for performing model validation.