Volume Census in T-S Space
Description¶
This recipe shows how to calculate the volume census in discrete temperature-salinity coordinates using annual-mean outputs from the National Oceanography Centre Near-Present-Day global eORCA1 configuration of NEMO forced using JRA55-do from 1976-2024.
For more details on this model configuration and the available outputs, users can explore the Near-Present-Day documentation here.
In [1]:
Copied!
# -- Import required packages -- #
import numpy as np
import xarray as xr
from nemo_cookbook import NEMODataTree
xr.set_options(display_style="text")
# -- Import required packages -- #
import numpy as np
import xarray as xr
from nemo_cookbook import NEMODataTree
xr.set_options(display_style="text")
Out[1]:
<xarray.core.options.set_options at 0x30df91010>
Using Dask¶
Optional: Connect Client to Dask Local Cluster to run analysis in parallel.
Note: Although using Dask is not strictly necessary for this simple example using eORCA1, if we wanted to generalise this recipe to eORCA025 or eORCA12 outputs, using Dask would be essential to avoid unnecessary slow calculations using only a single process.
In [ ]:
Copied!
# -- Initialise Dask Local Cluster -- #
import os
import dask
from dask.distributed import Client, LocalCluster
# Update temporary directory for Dask workers:
dask.config.set({'temporary_directory': f"{os.getcwd()}/dask_tmp",
'local_directory': f"{os.getcwd()}/dask_tmp"
})
# Create Local Cluster:
cluster = LocalCluster(n_workers=4, threads_per_worker=3, memory_limit='5GB')
client = Client(cluster)
client
# -- Initialise Dask Local Cluster -- #
import os
import dask
from dask.distributed import Client, LocalCluster
# Update temporary directory for Dask workers:
dask.config.set({'temporary_directory': f"{os.getcwd()}/dask_tmp",
'local_directory': f"{os.getcwd()}/dask_tmp"
})
# Create Local Cluster:
cluster = LocalCluster(n_workers=4, threads_per_worker=3, memory_limit='5GB')
client = Client(cluster)
client
Accessing NEMO Model Data¶
Let's begin by loading the grid variables for our eORCA1 NEMO model from the JASMIN Object Store.
Alternatively, you can replace the domain_filepath below with a file path to your domain_cfg.nc file and read this with xarray's open_dataset() function.
In [2]:
Copied!
# Define directory path to ancillary files:
domain_filepath = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain"
# Open eORCA1 NEMO model domain_cfg:
ds_domain = xr.open_zarr(f"{domain_filepath}/domain_cfg", consolidated=True, chunks={})
ds_subbasins = xr.open_zarr(f"{domain_filepath}/subbasins", consolidated=True, chunks={})
ds_domain
# Define directory path to ancillary files:
domain_filepath = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain"
# Open eORCA1 NEMO model domain_cfg:
ds_domain = xr.open_zarr(f"{domain_filepath}/domain_cfg", consolidated=True, chunks={})
ds_subbasins = xr.open_zarr(f"{domain_filepath}/subbasins", consolidated=True, chunks={})
ds_domain
Out[2]:
<xarray.Dataset> Size: 667MB
Dimensions: (t: 1, y: 331, x: 360, z: 75)
Dimensions without coordinates: t, y, x, z
Data variables: (12/54)
e1t (t, y, x) float64 953kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
e2v (t, y, x) float64 953kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
bottom_level (t, y, x) int32 477kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
e2t (t, y, x) float64 953kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
e2u (t, y, x) float64 953kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
bathy_metry (t, y, x) float32 477kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
... ...
nav_lev (z) float32 300B dask.array<chunksize=(75,), meta=np.ndarray>
mask_csundef (t, y, x) int32 477kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
nav_lat (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
nav_lon (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
time_counter (t) float64 8B dask.array<chunksize=(1,), meta=np.ndarray>
top_level (t, y, x) int32 477kB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>
Attributes:
DOMAIN_number_total: 1
DOMAIN_number: 0
DOMAIN_dimensions_ids: [1, 2]
DOMAIN_size_global: [362, 332]
DOMAIN_size_local: [362, 332]
DOMAIN_position_first: [1, 1]
DOMAIN_position_last: [362, 332]
DOMAIN_halo_size_start: [0, 0]
DOMAIN_halo_size_end: [0, 0]
DOMAIN_type: BOX
history: Mon Jun 5 12:41:32 2023: ncks -A mask.nc ORCA1_...
NCO: 4.4.7
Next, we need to import the sea water conservative temperature and absolute salinity stored at T-points in a single dataset.
In [3]:
Copied!
# Define directory path to model output files:
output_dir = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/T1y"
# Construct NEMO model grid dataset, including seawater conservative temperature (C) and absolute salinity (g/kg):
ds_gridT = xr.merge([xr.open_zarr(f"{output_dir}/{var}", consolidated=True, chunks={})[var] for var in ['thetao_con', 'so_abs', 'e3t']], compat="override")
ds_gridT
# Define directory path to model output files:
output_dir = "https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/T1y"
# Construct NEMO model grid dataset, including seawater conservative temperature (C) and absolute salinity (g/kg):
ds_gridT = xr.merge([xr.open_zarr(f"{output_dir}/{var}", consolidated=True, chunks={})[var] for var in ['thetao_con', 'so_abs', 'e3t']], compat="override")
ds_gridT
Out[3]:
<xarray.Dataset> Size: 5GB
Dimensions: (deptht: 75, y: 331, x: 360, time_counter: 48)
Coordinates:
* deptht (deptht) float32 300B 0.5058 1.556 ... 5.698e+03 5.902e+03
nav_lat (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
nav_lon (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
time_centered (time_counter) datetime64[ns] 384B dask.array<chunksize=(1,), meta=np.ndarray>
* time_counter (time_counter) datetime64[ns] 384B 1976-07-02 ... 2023-07-...
Dimensions without coordinates: y, x
Data variables:
thetao_con (time_counter, deptht, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
so_abs (time_counter, deptht, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
e3t (time_counter, deptht, y, x) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
Attributes:
cell_methods: time: mean
interval_operation: 1 yr
interval_write: 1 yr
long_name: sea_water_conservative_temperature
online_operation: average
standard_name: sea_water_conservative_temperature
units: degC
Creating a NEMODataTree¶
Next, let's create a NEMODataTree to store our domain and T- & V-grid variables for the eORCA1 model.
In [4]:
Copied!
# Define dictionary of grid datasets defining eORCA1 parent model domain with no child/grand-child nests:
# Note: domain_cfg z-dimension is expected to be named 'nav_lev'.
datasets = {"parent": {"domain": ds_domain.rename({"z": "nav_lev"}), "gridT": ds_gridT}}
# Initialise a new NEMODataTree whose parent domain is zonally periodic & north-folding on F-points:
nemo = NEMODataTree.from_datasets(datasets=datasets, iperio=True, nftype="F")
nemo
# Define dictionary of grid datasets defining eORCA1 parent model domain with no child/grand-child nests:
# Note: domain_cfg z-dimension is expected to be named 'nav_lev'.
datasets = {"parent": {"domain": ds_domain.rename({"z": "nav_lev"}), "gridT": ds_gridT}}
# Initialise a new NEMODataTree whose parent domain is zonally periodic & north-folding on F-points:
nemo = NEMODataTree.from_datasets(datasets=datasets, iperio=True, nftype="F")
nemo
Out[4]:
<xarray.DataTree>
Group: /
│ Dimensions: (time_counter: 48)
│ Coordinates:
│ time_centered (time_counter) datetime64[ns] 384B dask.array<chunksize=(1,), meta=np.ndarray>
│ * time_counter (time_counter) datetime64[ns] 384B 1976-07-02 ... 2023-07-...
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridT
│ Dimensions: (time_counter: 48, k: 75, j: 331, i: 360)
│ Coordinates:
│ * deptht (k) float32 300B 0.5058 1.556 2.668 ... 5.698e+03 5.902e+03
│ time_centered (time_counter) datetime64[ns] 384B dask.array<chunksize=(1,), meta=np.ndarray>
│ gphit (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamt (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331
│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
│ Data variables:
│ thetao_con (time_counter, k, j, i) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
│ so_abs (time_counter, k, j, i) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
│ e3t (time_counter, k, j, i) float32 2GB dask.array<chunksize=(1, 25, 331, 360), meta=np.ndarray>
│ e1t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ top_level (j, i) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ bottom_level (j, i) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ tmask (k, j, i) bool 9MB False False False ... False False False
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridU
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphiu (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamu (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331
│ * i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5
│ Data variables:
│ e1u (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2u (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ umask (k, j, i) bool 9MB False False False ... False False False
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridV
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphiv (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamv (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
│ * j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
│ Data variables:
│ e1v (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2v (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ vmask (k, j, i) bool 9MB False False False ... False False False
│ Attributes:
│ nftype: F
│ iperio: True
├── Group: /gridW
│ Dimensions: (j: 331, i: 360, k: 75)
│ Coordinates:
│ gphit (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ glamt (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ * k (k) float64 600B 0.5 1.5 2.5 3.5 4.5 ... 71.5 72.5 73.5 74.5
│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331
│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360
│ Data variables:
│ e1t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ e2t (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
│ wmask (k, j, i) bool 9MB False False False ... False False False
│ Attributes:
│ nftype: F
│ iperio: True
└── Group: /gridF
Dimensions: (j: 331, i: 360, k: 75)
Coordinates:
gphif (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
glamf (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
* k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75
* j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5
* i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5
Data variables:
e1f (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
e2f (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>
fmask (k, j, i) bool 9MB False False False ... False False False
Attributes:
nftype: F
iperio: True
Calculating Volume Census¶
Now we have constructed our NEMODataTree, let's calculate the volume census in in T-S coordinates using the .binned_statistic() method:
In [5]:
Copied!
# Define Atlantic Ocean basin mask:
atlmask = ds_subbasins['atlmsk'].rename({"x":"i", "y":"j"}).astype(bool)
# Define Atlantic Ocean basin mask:
atlmask = ds_subbasins['atlmsk'].rename({"x":"i", "y":"j"}).astype(bool)
In [6]:
Copied!
# Define discrete conservative temperature [C] and absolute salinity [g/kg] bins:
thetao_bins = np.arange(-2, 35, 0.5)
so_bins = np.arange(20, 38, 0.1)
# Compute volume of each T-grid cell [m3].
nemo['gridT']['volcello'] = nemo.cell_volume(grid='/gridT')
# Compute total volume in discrete conservative temperature - absolute salinity coords:
vol_thetao_so_atl = nemo.binned_statistic(grid="/gridT",
vars=["thetao_con", "so_abs"],
values="volcello",
keep_dims=["time_counter"],
bins=[thetao_bins, so_bins],
statistic="nansum",
mask=atlmask
)
vol_thetao_so_atl
# Define discrete conservative temperature [C] and absolute salinity [g/kg] bins:
thetao_bins = np.arange(-2, 35, 0.5)
so_bins = np.arange(20, 38, 0.1)
# Compute volume of each T-grid cell [m3].
nemo['gridT']['volcello'] = nemo.cell_volume(grid='/gridT')
# Compute total volume in discrete conservative temperature - absolute salinity coords:
vol_thetao_so_atl = nemo.binned_statistic(grid="/gridT",
vars=["thetao_con", "so_abs"],
values="volcello",
keep_dims=["time_counter"],
bins=[thetao_bins, so_bins],
statistic="nansum",
mask=atlmask
)
vol_thetao_so_atl
Out[6]:
<xarray.DataArray 'volcello' (time_counter: 48, thetao_con_bins: 73,
so_abs_bins: 179)> Size: 5MB
dask.array<reshape, shape=(48, 73, 179), dtype=float64, chunksize=(48, 73, 179), chunktype=numpy.ndarray>
Coordinates:
* time_counter (time_counter) datetime64[ns] 384B 1976-07-02 ... 2023-0...
* thetao_con_bins (thetao_con_bins) float64 584B -1.75 -1.25 ... 33.75 34.25
* so_abs_bins (so_abs_bins) float64 1kB 20.05 20.15 20.25 ... 37.75 37.85
Notice that the output above contains dask arrays, so we haven't actually computed the volume census yet. To do this, we need to call the .compute() method:
In [7]:
Copied!
vol_thetao_so_atl = vol_thetao_so_atl.compute()
vol_thetao_so_atl = vol_thetao_so_atl.compute()
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
Visualising the time-mean volume census in T-S coordinates¶
Finally, let's visualise the results by plotting the time-mean volume census in conservative temperature - absolute salinity space:
In [8]:
Copied!
vol_thetao_so_atl.mean(dim='time_counter').plot()
vol_thetao_so_atl.mean(dim='time_counter').plot()
Out[8]:
<matplotlib.collections.QuadMesh at 0x33a884050>
