Getting Started with NEMO Cookbook

PROMOTE Sprint - Day 1¶

Tutorial: Getting Started with NEMO Cookbook¶

13th July 2026¶

Contact:¶

Ollie Tooth (oliver.tooth@noc.ac.uk), National Oceanography Centre

Background:¶

NEMO Cookbook extends the familiar xarray data model with grid-aware data structures designed for performing reproducible analyses of the Nucleus for European Modelling of the Ocean (NEMO) ocean general circulation model outputs.

In these hands-on tutorials, we will use the NEMODataTree object to analyse ocean & sea-ice outputs of the UKESM1-2-LL TIPMIP idealised emissions scenarios and the CANARI Large-Ensemble.

We will cover:

Creating NEMODataTree objects directly from netCDF files, including CMORISED outputs from coupled climate models, and from virtual Icechunk repositories.
Exploring NEMO ocean & sea-ice variables stored in a NEMODataTree.
Geographical plotting of NEMODataArray.
Performing grid-aware operations using NEMODataArray.

1.1 Creating NEMODataTrees¶

Let's begin by importing the NEMO Cookbook Python library and the NEMODataTree class:

In [ ]:

Copied!





import os

import cartopy.crs as ccrs
import icechunk
import matplotlib.pyplot as plt
import nemo_cookbook as ncb
import pandas as pd
import xarray as xr
from nemo_cookbook import NEMODataTree
from xarray.coders import CFDatetimeCoder
import os

import cartopy.crs as ccrs
import icechunk
import matplotlib.pyplot as plt
import nemo_cookbook as ncb
import pandas as pd
import xarray as xr
from nemo_cookbook import NEMODataTree
from xarray.coders import CFDatetimeCoder

1.1.1 Creating a NEMODataTree from Local Files¶

We'll begin by creating a simple NEMODataTree using local NEMO model output files from the AGRIF_DEMO reference configuration.

AGRIF_DEMO

AGRIF_DEMO is based on the ORCA2_ICE_PISCES global reference configuration with the inclusion of 3 online nested domains.
Here, we will only consider the 2° global parent domain.
NEMO Cookbook includes a subset of the AGRIF_DEMO model outputs made accessible via cloud object storage.
Further information on this reference configuration can be found here.

Let's download and collect the filepaths for the AGRIF_DEMO configuration using the nemo_cookbook.examples.get_filepaths() convenience function, which caches and generates local filepaths for available NEMO reference configurations:

In [ ]:

Copied!

filepaths = ncb.examples.get_filepaths("AGRIF_DEMO")
filepaths
filepaths = ncb.examples.get_filepaths("AGRIF_DEMO")
filepaths

We can create a NEMODataTree from a dictionary of paths to local netCDF files using the .from_paths() constructor:

from_paths()

Creates a NEMODataTree from a dictionary of paths to NEMO model outputs files organised into a hierarchy (e.g., parent, child grandchild) of NEMO model grids (e.g., gridT, gridU, etc.):

    <xarray.DataTree 'nemo'>
    Group: /
    ├── Group: /gridT
    ├── Group: /gridU
    ├── Group: /gridV
    ├── Group: /gridW
    └── Group: /gridF

where the gridT node contains time series of scalar variables stored in the ...grid_T.nc files in a single xarray.Dataset and so on.

Let's start by defining the paths dictionary, which contains the filepaths corresponding to our global parent domain.
We populate the parent dictionary with the filepaths to the domain_cfg and gridT/U/V/W netCDF files produced for the AGRIF_DEMO (ORCA2) parent domain.

In [ ]:

Copied!





paths = {"parent": {"domain": filepaths["domain_cfg.nc"],
                    "gridT": filepaths["ORCA2_5d_00010101_00010110_grid_T.nc"],
                    "gridU": filepaths["ORCA2_5d_00010101_00010110_grid_U.nc"],
                    "gridV": filepaths["ORCA2_5d_00010101_00010110_grid_V.nc"],
                    "gridW": filepaths["ORCA2_5d_00010101_00010110_grid_W.nc"],
                    "icemod": filepaths["ORCA2_5d_00010101_00010110_icemod.nc"]
                    },
         }
paths = {"parent": {"domain": filepaths["domain_cfg.nc"],
                    "gridT": filepaths["ORCA2_5d_00010101_00010110_grid_T.nc"],
                    "gridU": filepaths["ORCA2_5d_00010101_00010110_grid_U.nc"],
                    "gridV": filepaths["ORCA2_5d_00010101_00010110_grid_V.nc"],
                    "gridW": filepaths["ORCA2_5d_00010101_00010110_grid_W.nc"],
                    "icemod": filepaths["ORCA2_5d_00010101_00010110_icemod.nc"]
                    },
         }

We can now construct a new NEMODataTree using the .from_paths() constructor.
Note that we also need to specify that our global parent domain is zonally periodic (iperio=True) and north folding on T-points (nftype = "T") rather than a closed (regional) domain.

In [ ]:

Copied!

NEMODataTree.from_paths(paths, iperio=True, nftype="T", chunks={})
NEMODataTree.from_paths(paths, iperio=True, nftype="T", chunks={})

1.1.2 Creating a NEMODataTree from CMORISED Datasets¶

In addition to reading netCDF files directly, we can also create a NEMODataTree from a dictionary of xarray.Datasets using the .from_datasets() constructor:

from_datasets()

Creates a NEMODataTree from a dictionary of xarray.Datasets organised into a hierarchy (e.g., parent, child grandchild) of NEMO model grids (e.g., gridT, gridU, etc.):

    <xarray.DataTree 'nemo'>
    Group: /
    ├── Group: /gridT
    ├── Group: /gridU
    ├── Group: /gridV
    ├── Group: /gridW
    └── Group: /gridF

where the gridT node contains time series of scalar variables stored in the input xarray.Dataset and so on.

Here, we show how to construct a NEMODataTree from the 1° global ocean sea-ice component of UKESM1-2-LL using outputs accessible via the terrafirma group workspace.
Since CMIP6 outputs are processed and formatted according to the CMIP Community Climate Model Output Rewriter (CMOR) software, we will need to include a few additional pre-processing steps to reformat our NEMO model outputs in order to construct a NEMODataTree

Important: only CMIP model outputs variables stored on their original NEMO ocean model grid (i.e, gn) can be used to construct a NEMODataTree

In [ ]:

Copied!

# -- Open UKESM1-2-LL domain_cfg -- #
ds_domain = xr.open_dataset("/gws/ssde/j25b/terrafirma/otooth/Ofx/domain_cfg_UKESM1.nc")

ds_domain
# -- Open UKESM1-2-LL domain_cfg -- #
ds_domain = xr.open_dataset("/gws/ssde/j25b/terrafirma/otooth/Ofx/domain_cfg_UKESM1.nc")

ds_domain

In [ ]:

Copied!





# -- Create NEMO gridT dataset from UKESM1-2-LL -- #
# Define time decoder to handle CMIP6 time units:
time_decoder = xr.coders.CFDatetimeCoder(use_cftime=True)

# Open UKESM1-2-LL thetao dataset, including potential temperature (°C):
base_filepath = "/gws/ssde/j25b/terrafirma/TerraFIRMA/MOHC/UKESM1-2/esm-hist/r1i1p1f1/Omon"
ds_ukesm_gridT = xr.open_mfdataset(f"{base_filepath}/thetao/gn/v20240619/thetao_Omon_UKESM1-2-LL_esm-hist_r1i1p1f1_gn_*.nc",
                                   data_vars='all',
                                   decode_times=time_decoder
                                   )

# Adding UKESM1-2-LL thkcello dataset, including vertical cell thickness (m):
ds_ukesm_gridT['thkcello'] = xr.open_mfdataset(f"{base_filepath}/thkcello/gn/v20240619/thkcello_Omon_UKESM1-2-LL_esm-hist_r1i1p1f1_gn_*.nc",
                                               data_vars='all',
                                               decode_times=time_decoder
                                               )['thkcello']

# Rename variables to NEMO standard names:
ds_ukesm_gridT = ds_ukesm_gridT.rename({'time': 'time_counter',
                                        'i': 'x',
                                        'j': 'y',
                                        'lev': 'deptht',
                                        'thkcello': 'e3t'}
                                        )

ds_ukesm_gridT
# -- Create NEMO gridT dataset from UKESM1-2-LL -- #
# Define time decoder to handle CMIP6 time units:
time_decoder = xr.coders.CFDatetimeCoder(use_cftime=True)

# Open UKESM1-2-LL thetao dataset, including potential temperature (°C):
base_filepath = "/gws/ssde/j25b/terrafirma/TerraFIRMA/MOHC/UKESM1-2/esm-hist/r1i1p1f1/Omon"
ds_ukesm_gridT = xr.open_mfdataset(f"{base_filepath}/thetao/gn/v20240619/thetao_Omon_UKESM1-2-LL_esm-hist_r1i1p1f1_gn_*.nc",
                                   data_vars='all',
                                   decode_times=time_decoder
                                   )

# Adding UKESM1-2-LL thkcello dataset, including vertical cell thickness (m):
ds_ukesm_gridT['thkcello'] = xr.open_mfdataset(f"{base_filepath}/thkcello/gn/v20240619/thkcello_Omon_UKESM1-2-LL_esm-hist_r1i1p1f1_gn_*.nc",
                                               data_vars='all',
                                               decode_times=time_decoder
                                               )['thkcello']

# Rename variables to NEMO standard names:
ds_ukesm_gridT = ds_ukesm_gridT.rename({'time': 'time_counter',
                                        'i': 'x',
                                        'j': 'y',
                                        'lev': 'deptht',
                                        'thkcello': 'e3t'}
                                        )

ds_ukesm_gridT

Now we have defined our domain and gridT datasets, let's define a datasets dictionary ensuring that we rename CMORISED dimensions to be consistent with standard NEMO model outputs.

Note that we must rename the CMORISED variable thkcello to standard NEMO name e3t

We can then define a NEMODataTree using the .from_datasets() constructor, specifying that our global parent domain is zonally periodic and north-folding on F-points.

Note that we will still obtain a NEMODataTree with a complete set of grid nodes (i.e., gridT, gridU, etc.)

In [ ]:

Copied!

datasets = {"parent": {"domain": ds_domain, "gridT": ds_ukesm_gridT}}

NEMODataTree.from_datasets(datasets=datasets, iperio=True, nftype="F", read_mask=True)
datasets = {"parent": {"domain": ds_domain, "gridT": ds_ukesm_gridT}}

NEMODataTree.from_datasets(datasets=datasets, iperio=True, nftype="F", read_mask=True)

1.1.3 Virtual NEMODataTrees¶

Rather that having to merge O(100) - O(1000s) of netCDF files to create a NEMODataTree at the begining of every script, we can alternatively build a virtual NEMODataTree:

Virtual NEMODataTree -> lightweight, persistent map describing how to organise our netCDF files into a NEMODataTree.

For the PROMOTE Analysis Sprint, we have created virtual NEMODataTrees for the following simulations:
- UKESM1-2-LL: /gws/ssde/j25b/terrafirma/otooth/repos/
- CANARI-LE: /gws/ssde/j25b/canari/shared/large-ensemble/ocean/repos/
To access virtual NEMODataTrees, we can use the .from_icechunk() constructor:

from_icechunk()

Creates a NEMODataTree from an Icechunk repository:

    <xarray.DataTree 'nemo'>
    Group: /
    ├── Group: /gridT
    ├── Group: /gridU
    ├── Group: /gridV
    ├── Group: /gridW
    └── Group: /gridF

where each grid node is stored as a virtual dataset in separate groups of a hierarchical Icechunk store.

Here, we show how to construct a NEMODataTree from the 1° global ocean sea-ice component of UKESM1-2-LL using the available Icechunk repositories accessible via the terrafirma group workspace.
Let's start by exploring the available Icechunk repositories...

In [ ]:

Copied!

os.listdir("/gws/ssde/j25b/terrafirma/otooth/repos/")
os.listdir("/gws/ssde/j25b/terrafirma/otooth/repos/")

In [ ]:

Copied!

os.listdir("/gws/ssde/j25b/canari/shared/large-ensemble/ocean/repos/")
os.listdir("/gws/ssde/j25b/canari/shared/large-ensemble/ocean/repos/")

Next, let's define a reusable function open_virtual_NEMODataTree() to open a virtual NEMODataTree given the experiment and member of a UKESM1-2-LL simulation...

In [ ]:

Copied!





def open_virtual_NEMODataTree(
        experiment: str,
        member: str = "r1i1p1f1"
    ) -> NEMODataTree:

    # Authorise Icechunk to fetch virtual chunks from the original netCDF files on JASMIN terrafirma group workspace:
    store_fpath = "/gws/ssde/j25b/terrafirma"
    credentials = icechunk.credentials.containers_credentials({f"file://{store_fpath}/": icechunk.credentials.LocalFileSystemAccess})

    # Open Icechunk repository stored on JASMIN terrafirma group workspace:
    repo_fpath = f"/gws/ssde/j25b/terrafirma/otooth/repos/UKESM1-2-LL_{experiment}_{member}"
    repo = icechunk.Repository.open(storage=icechunk.local_filesystem_storage(repo_fpath),
                                    authorize_virtual_chunk_access=credentials
                                    )
    
    # Create virtual NEMODataTree from Icechunk repository:
    nemo = NEMODataTree.from_icechunk(repo=repo,
                                      branch="main",
                                      iperio=True,
                                      nftype="F",
                                      name=f"UKESM1-2-LL_{experiment}_{member}",
                                      open_kwargs={"decode_times": CFDatetimeCoder(time_unit="s"),
                                                   "chunks": {},
                                                   }
                                      )

    return nemo
def open_virtual_NEMODataTree(
        experiment: str,
        member: str = "r1i1p1f1"
    ) -> NEMODataTree:

    # Authorise Icechunk to fetch virtual chunks from the original netCDF files on JASMIN terrafirma group workspace:
    store_fpath = "/gws/ssde/j25b/terrafirma"
    credentials = icechunk.credentials.containers_credentials({f"file://{store_fpath}/": icechunk.credentials.LocalFileSystemAccess})

    # Open Icechunk repository stored on JASMIN terrafirma group workspace:
    repo_fpath = f"/gws/ssde/j25b/terrafirma/otooth/repos/UKESM1-2-LL_{experiment}_{member}"
    repo = icechunk.Repository.open(storage=icechunk.local_filesystem_storage(repo_fpath),
                                    authorize_virtual_chunk_access=credentials
                                    )
    
    # Create virtual NEMODataTree from Icechunk repository:
    nemo = NEMODataTree.from_icechunk(repo=repo,
                                      branch="main",
                                      iperio=True,
                                      nftype="F",
                                      name=f"UKESM1-2-LL_{experiment}_{member}",
                                      open_kwargs={"decode_times": CFDatetimeCoder(time_unit="s"),
                                                   "chunks": {},
                                                   }
                                      )

    return nemo

In [ ]:

Copied!

# Open virtual NEMODataTree from UKESM1-2-LL esm-up2p0 Icechunk repository:
nemo = open_virtual_NEMODataTree(experiment="esm-up2p0", member="r1i1p1f1")

nemo
# Open virtual NEMODataTree from UKESM1-2-LL esm-up2p0 Icechunk repository:
nemo = open_virtual_NEMODataTree(experiment="esm-up2p0", member="r1i1p1f1")

nemo

1.2 Navigating a NEMODataTree¶

Now we have created an example NEMODataTree, let's take a closer look at its contents and how to access NEMO model output variables.
We'll start by looking at the nodes in our NEMODataTree:

In [ ]:

Copied!

nemo.groups
nemo.groups

Some Important Points:

A NEMODataTree does not have a domain node since grid scale factors and masks associated with each model domain are assigned to their respective grid nodes during pre-processing (e.g., horizontal grid scale factors e1t and e2t are stored in gridT etc.).
During the construction of a NEMODataTree, the standard (depth{p}, y, x) dimensions of NEMO model outputs (where p is the grid point type) are transformed to NEMO grid indices (i, j, k). This has two important implications:
- xarray.Datasets stored in each grid node share the same coordinate dimension names (i, j, k), but are staggered according to the location of variables on the NEMO model grid.
- All grid indices use Fortran (1-based) indexing to be consistent with the original NEMO model code.
| Grid Type | Grid Indices | | ----------- | -------------------------- | | T | $(i, j, k)$ | | U | $(i + \frac{1}{2}, j, k)$ | | V | $(i, j + \frac{1}{2}, k)$ | | W | $(i, j, k + \frac{1}{2})$ | | F | $(i + \frac{1}{2}, j + \frac{1}{2}, k)$ |

To see what this means in practice, we will consider the gridT and gridU nodes of the 1° UKESM1-2-LL NEMODataTree below:

In [ ]:

Copied!

nemo['gridT']
nemo['gridT']

In [ ]:

Copied!

nemo['gridU']
nemo['gridU']

Above, we have seen how to access the nodes in an NEMODataTree using a dictionary-like syntax.
We can also extend this to include output variables stored in given grid node as follows...
To access an unmasked variable (i.e., unchanged from the original model output files) as an xarray.DataArray:

In [ ]:

Copied!

nemo['gridT']['thetao']
nemo['gridT']['thetao']

To access a unmasked variable as a NEMODataArray (i.e., NEMO grid-aware extension of xarray.DataArray), we can instead provide a direct path to the variable:

In [ ]:

Copied!

nemo["gridT/thetao"]
nemo["gridT/thetao"]

To access a masked (i.e., the appropriate land-sea mask is applied automatically) variable as an NEMODataArray, we can instead provide a direct path to the variable:

In [ ]:

Copied!

nemo["gridT/thetao"].masked
nemo["gridT/thetao"].masked

1.3 NEMODataTree Grid Properties¶

We often need access to variables describing the NEMO ocean model domain to calculate derived quantities, such as volume, heat and salt transports.
For convenience, the NEMODataTree class includes methods for calculating grid cell areas and volumes.
For example, to compute the horizontal area of grid cells centered on T-points in the parent domain, we can use the .cell_area() method:

cell_area()

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 |

In [ ]:

Copied!

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

Next, let's compute the volume of each grid cell centered on a T grid point in the model parent domain using the .cell_volume() method:

In [ ]:

Copied!

nemo.cell_volume(grid="gridT")
nemo.cell_volume(grid="gridT")

1.4 NEMODataArray Properties¶

In addition to simply storing ocean model output variable, a NEMODataArray also gives each variable knowledge of its:
- NEMO model grid location (e.g., T, U, V, W, F)
- parent NEMODataTree
- associated NEMO grid metrics (grid scale factors)

Each NEMODataArray has several key properties to support grid-aware computation...
- .data Underlying xarray.DataArray of the NEMO output variable.
- .grid Path to NEMO model grid node where variable is stored.
- .grid_type Type of NEMO model grid where variable is defined.
- .metrics Dictionary of NEMO model grid scale factors (e.g., e1t, e2t etc.) associated with the variable.
- .mask Variable land-sea mask (xarray.DataArray).
- .masked Returns variable NEMODataArray with land-sea mask applied.

In [ ]:

Copied!

# Collect grid scale factors of sea water potential temperature (thetao) variable on T-grid:
nemo['gridT/thetao'].metrics
# Collect grid scale factors of sea water potential temperature (thetao) variable on T-grid:
nemo['gridT/thetao'].metrics

In [ ]:

Copied!

# Collect land-sea mask of sea water potential temperature (thetao) variable on T-grid:
nemo['gridT/thetao'].mask
# Collect land-sea mask of sea water potential temperature (thetao) variable on T-grid:
nemo['gridT/thetao'].mask

1.5 NEMODataArray Core Operations¶

Summary of Core Operations

Plotting: geoplot()
Indexing: sel_like()
Clipping & Masking: clip(), apply_mask()
Grid Operators: diff(), derivative(), integral(), depth_integral().
Statistics: weighted_mean(), masked_statistic().
Grid Transformations: interp_to(), transform_vertical_grid().
For more details on the above core operations see the follow section of the How To... Guide.

1.5.0 Plotting¶

Since NEMODataArray is a wrapper around the familiar xarray.DataArray, we can still access xarray.DataArray.plot for exploratory plotting.
However, for geographical plots of 2-dimensional NEMODataArrays with dimensions (j, i) only, we can use the geoplot() method to produce a Cartopy GeoQuadMesh as follows:

In [ ]:

Copied!

# -- Plot example sea surface temperature (°C) -- #
nemo["gridT/tos"].isel(time_counter=0).geoplot()
# -- Plot example sea surface temperature (°C) -- #
nemo["gridT/tos"].isel(time_counter=0).geoplot()

In [ ]:

Copied!





# -- Plot example sea ice concentration in the Northern Hemisphere -- #
proj = ccrs.NorthPolarStereo()

(nemo["gridT/siconc"]
 .isel(time_counter=0)
 .geoplot(projection=proj,
          vmin=0, vmax=1,
          extent=(-180, 180, 50, 90),
          cmap='Blues_r',
          clabel_kwargs={'label': 'Sea Ice Concentration [fraction of unity]'}
          )
 )
# -- Plot example sea ice concentration in the Northern Hemisphere -- #
proj = ccrs.NorthPolarStereo()

(nemo["gridT/siconc"]
 .isel(time_counter=0)
 .geoplot(projection=proj,
          vmin=0, vmax=1,
          extent=(-180, 180, 50, 90),
          cmap='Blues_r',
          clabel_kwargs={'label': 'Sea Ice Concentration [fraction of unity]'}
          )
 )

1.5.1 Indexing¶

Alongside 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:

In [ ]:

Copied!

nda = nemo["gridT/thetao"].sel(i=slice(1, 90), k=20)

nemo["gridT/thetao"].sel_like(nda)
nda = nemo["gridT/thetao"].sel(i=slice(1, 90), k=20)

nemo["gridT/thetao"].sel_like(nda)

We can also use the add_geoindex() method of our NEMODataTree to enable indexing of a NEMO model variables with geographical coordinates:

Note that the .dataset property must be used to access the xarray.Dataset view of the NEMO model grid node

In [ ]:

Copied!

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

nemo_geo["gridT"].dataset.sel(gphit=60, glamt=-40, method="nearest")
nemo_geo = nemo.add_geoindex(grid="gridT")

nemo_geo["gridT"].dataset.sel(gphit=60, glamt=-40, method="nearest")

1.5.2 Clipping & Masking¶

We often want to perform a calculation using only NEMO grid points within a specific region.
To clip a NEMODataArray to include only data within a geographical bounding box, we can use the clip() method.

In [ ]:

Copied!





# Define subpolar North Atlantic boundary box for analysis:
bbox_spna = (-68, -2, 45, 67)  # (-68°E, -2°E, 45°N, 67°N)

# Plot the sea surface temperature (°C) in the subpolar North Atlantic region:
(nemo["gridT/tos"]
 .clip(bbox_spna)
 .isel(time_counter=0)
 .geoplot(projection=ccrs.Orthographic(central_longitude=-35, central_latitude=50),
          extent=(-75, 0, 42, 70), # (-75°E, -0°E, 42°N, 70°N)
          cmap='RdBu_r',
          clabel_kwargs={'label': 'Sea Surface Temperature (°C)'}
          )
 )
# Define subpolar North Atlantic boundary box for analysis:
bbox_spna = (-68, -2, 45, 67)  # (-68°E, -2°E, 45°N, 67°N)

# Plot the sea surface temperature (°C) in the subpolar North Atlantic region:
(nemo["gridT/tos"]
 .clip(bbox_spna)
 .isel(time_counter=0)
 .geoplot(projection=ccrs.Orthographic(central_longitude=-35, central_latitude=50),
          extent=(-75, 0, 42, 70), # (-75°E, -0°E, 42°N, 70°N)
          cmap='RdBu_r',
          clabel_kwargs={'label': 'Sea Surface Temperature (°C)'}
          )
 )

To define a regional mask using the geographical coordinates of a closed polygon, we can use the mask_with_polygon() method.
In this example, we use the IHO World Seas v3 Labrador Sea polygon to generate the longitude and latitude coordinates of the closed polygon to mask.

In [ ]:

Copied!

# Open IHO World Seas v3 polygons from JASMIN Object Store as Pandas DataFrame:
filepaths = ncb.examples.get_filepaths("IHO")

df_IHO_World_Seas = pd.read_parquet(filepaths['IHO_World_Seas_v3_polygons.parquet'])

# Define IHO World Seas Labrador Sea polygon:
lon_poly = df_IHO_World_Seas[df_IHO_World_Seas['Name'] == 'Labrador Sea']['Longitudes'].item()[0]
lat_poly = df_IHO_World_Seas[df_IHO_World_Seas['Name'] == 'Labrador Sea']['Latitudes'].item()[0]

# Define boolean mask for Labrador Sea:
LSea_mask = nemo.mask_with_polygon(grid='gridT', lon_poly=lon_poly, lat_poly=lat_poly)

LSea_mask
# Open IHO World Seas v3 polygons from JASMIN Object Store as Pandas DataFrame:
filepaths = ncb.examples.get_filepaths("IHO")

df_IHO_World_Seas = pd.read_parquet(filepaths['IHO_World_Seas_v3_polygons.parquet'])

# Define IHO World Seas Labrador Sea polygon:
lon_poly = df_IHO_World_Seas[df_IHO_World_Seas['Name'] == 'Labrador Sea']['Longitudes'].item()[0]
lat_poly = df_IHO_World_Seas[df_IHO_World_Seas['Name'] == 'Labrador Sea']['Latitudes'].item()[0]

# Define boolean mask for Labrador Sea:
LSea_mask = nemo.mask_with_polygon(grid='gridT', lon_poly=lon_poly, lat_poly=lat_poly)

LSea_mask

Next, let's create a plot applying our regional mask to a given NEMODataArray variable using the apply_mask() method.

Note that using apply_mask() is preferred over where() since this automatically applies the appropriate land-sea mask in addition to our custom mask

In [ ]:

Copied!





# Plotting Labrador Sea sub-domain sea surface temperature:
(nemo['gridT/tos']
 .isel(time_counter=0)
 .apply_mask(LSea_mask)
 .geoplot(projection=ccrs.Orthographic(central_longitude=-35, central_latitude=50),
          extent=(-75, 0, 42, 70), # (-75°E, -0°E, 42°N, 70°N)
          cmap='RdBu_r',
          clabel_kwargs={'label': 'Sea Surface Temperature (°C)'}
          )
 )

# Add IHO World Seas Labrador Sea polygon:
plt.plot(lon_poly, lat_poly, transform=ccrs.PlateCarree(), color='0.1', linewidth=1.5)
# Plotting Labrador Sea sub-domain sea surface temperature:
(nemo['gridT/tos']
 .isel(time_counter=0)
 .apply_mask(LSea_mask)
 .geoplot(projection=ccrs.Orthographic(central_longitude=-35, central_latitude=50),
          extent=(-75, 0, 42, 70), # (-75°E, -0°E, 42°N, 70°N)
          cmap='RdBu_r',
          clabel_kwargs={'label': 'Sea Surface Temperature (°C)'}
          )
 )

# Add IHO World Seas Labrador Sea polygon:
plt.plot(lon_poly, lat_poly, transform=ccrs.PlateCarree(), color='0.1', linewidth=1.5)

For more details on clipping grids and domains within a NEMODataTree see the following section of the How To... Guide.

1.5.3 Grid Operators¶

Grid operators are the foundational set of finite difference operators defined on the NEMO model grid.
For example, the difference operator $\delta_{i} [q]$ defined by:

$$\delta_{i} [q] = q(i + 1/2) - q(i - 1/2)$$

For a complete list of grid operators implemented in NEMO Cookbook, see the following sections of the How To... Guide.
Here, we will focus on the integral() method, which can be used to integrate a variable var along one or more dimensions dims of a given NEMO model grid grid.
As an example, let's compute the total March mixed layer volume in the subpolar North Atlantic in esm-up2p0...

In [ ]:

Copied!





(nemo["gridT/mlotst"]
 # Select only March months:
 .sel(time_counter=nemo["gridT/time_counter"].dt.month == 3)
 # Clip to the SPNA bbox (-68°E, -2°E, 45°N, 67°N):
 .clip(bbox_spna)
 # Calculate area integral of mixed layer depth (m -> m³) in the SPNA:
 .integral(dims=["i", "j"])
 # Plot March mixed layer volume in the SPNA (m³):
 .plot(linewidth=2, color='0.1')
 )
(nemo["gridT/mlotst"]
 # Select only March months:
 .sel(time_counter=nemo["gridT/time_counter"].dt.month == 3)
 # Clip to the SPNA bbox (-68°E, -2°E, 45°N, 67°N):
 .clip(bbox_spna)
 # Calculate area integral of mixed layer depth (m -> m³) in the SPNA:
 .integral(dims=["i", "j"])
 # Plot March mixed layer volume in the SPNA (m³):
 .plot(linewidth=2, color='0.1')
 )

We can also use the .integral() method to calculate cumulative integrals along one or more dimensions of a given NEMO model grid.
For more details on cumulative integration see the following section of the How To... Guide.

1.5.4 Grid Transformations¶

Given that scalar and vector NEMO output variables are not co-located on the model grid, we often need to transform a variable defined on a given NEMO horizontal grid to a neighbouring grid via linear interpolation.
For example, we can transform the seawater potential temperature thetao defined on scalar T-points to neighbouring V-points in a NEMO model parent domain using the .interp_to() method:

In [ ]:

Copied!

nemo['gridT/thetao'].interp_to(to='V')
nemo['gridT/thetao'].interp_to(to='V')

Note that the original T-grid coordinate variables have been replaced by those of the V-grid onto which the conservative temperature field has been linearly interpolated.

We can also transform a variable defined on a given NEMO model vertical grid to a new vertical grid using conservative interpolation via the .transform_vertical_grid() method.
For more details on vertical grid transformations see the following section of the How To... Guide.

1.5.5 Statistics¶

Finally, we often want to calculate grid-aware statistics, such as area-weighted averages, using our NEMO model outputs.
Here, we leverage method chaining and the weighted_mean() method to calculate the horizontal grid cell area-weighted average sea surface temperature in the Northern Hemisphere:

In [ ]:

Copied!





# Calculate the area weighted-mean sea surface temperature in the Labrador Sea:
tos_LSea_wmean = nemo["gridT/tos"].apply_mask(LSea_mask).weighted_mean(dims=["i", "j"])

(tos_LSea_wmean
 .rolling(time_counter=120, center=True) # -> 120-month rolling average.
 .mean()
 .plot(linewidth=2, color='0.1')
)
# Calculate the area weighted-mean sea surface temperature in the Labrador Sea:
tos_LSea_wmean = nemo["gridT/tos"].apply_mask(LSea_mask).weighted_mean(dims=["i", "j"])

(tos_LSea_wmean
 .rolling(time_counter=120, center=True) # -> 120-month rolling average.
 .mean()
 .plot(linewidth=2, color='0.1')
)

We can also calculate aggregated statistics of a variable binned according to the values of one or more variables using the NEMODataTree.binned_statistic() method.
For more details on statistical methods see the following section of the How To... Guide.