Compute MOC along the approximate OSNAP array from ECCO

Here we compute volumetric transport along the OSNAP lines in depth space, which can be compared to recent observations.

Here we show:

  • how to get masks denoting the great circle arc between two points in space

  • how to compute the transport or streamfunction across this section

  • a comparison to observations

Datasets to download

If you don’t have any of the following datasets already, you will need to download them to complete the tutorial. Aside from the grid geometry file (which has no time dimension), you will need 2 monthly datasets for the full time span of ECCOv4r4 output (1992 through 2017). The ShortNames of the datasets are:




If you haven’t yet been through the download tutorial or used the ecco_download module, it may help you to review that information before downloading the datasets.

import warnings
import os
import sys
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import glob
import cartopy as cart
import as ccrs
notebook_path = os.getcwd()
# setting up a dask LocalCluster
from dask.distributed import Client
from dask.distributed import LocalCluster
cluster = LocalCluster()
client = Client(cluster)



Connection method: Cluster object Cluster type: distributed.LocalCluster

Cluster Info

## Import the ecco_v4_py library into Python
## =========================================
##    If ecco_v4_py is not installed in your local Python library,
##    tell Python where to find it.  The example below adds
##    ecco_v4_py to the user's path if it is stored in the folder
##    ECCOv4-py under the user's home directory

from os.path import join,expanduser,exists,split
user_home_dir = expanduser('~')


import ecco_v4_py as ecco
## Set top-level file directory for the ECCO NetCDF files
## =================================================================

## currently set to ~/Downloads/ECCO_V4r4_PODAAC,
## the default if ecco_podaac_download was used to download dataset granules
ECCO_dir = join(user_home_dir,'Downloads','ECCO_V4r4_PODAAC')
## Load the model grid
ecco_grid = xr.open_dataset(glob.glob(join(ECCO_dir,'*GEOMETRY*','*.nc'))[0])

## Load ECCO variables
ecco_vars_adv_th = xr.open_mfdataset(join(ECCO_dir,'*_OCEAN_3D_TEMPERATURE_FLUX*MONTHLY*','*.nc'),\
ecco_vars_uvw = xr.open_mfdataset(join(ECCO_dir,'*_OCEAN_3D_VOLUME_FLUX*MONTHLY*','*.nc'),\

ecco_vars = xr.merge((ecco_vars_adv_th[['ADVx_TH','ADVy_TH']],ecco_vars_uvw[['UVELMASS','VVELMASS']]))

## Merge the ecco_grid with the ecco_vars to make the ecco_ds
ds = xr.merge((ecco_grid , ecco_vars))

Load O-SNAP data

These data are available at and are published in

Lozier, M. S., Li, F., Bacon, S., Bahr, F., Bower, A. S., Cunningham, S. A., … Zhao, J. (2019). A sea change in our view of overturning in the subpolar North Atlantic. Science, 363(6426), 516–521.

Download these files from the o-snap data page: : :

and put them in a path:

# change this path as needed based on where you have stored the O-SNAP data files above
osnap_data_dir = join(user_home_dir,'Downloads','OSNAP')
obs1 = xr.open_dataset(join(osnap_data_dir,''))
obs2 = xr.open_dataset(join(osnap_data_dir,''))
obs = xr.merge((obs1,obs2))
Dimensions:       (TIME: 47)
  * TIME          (TIME) datetime64[ns] 2014-07-31 2014-08-30 ... 2018-05-11
Data variables: (12/24)
    MOC_ALL       (TIME) float64 ...
    MOC_ALL_ERR   (TIME) float64 ...
    MOC_WEST      (TIME) float64 ...
    MOC_WEST_ERR  (TIME) float64 ...
    MOC_EAST      (TIME) float64 ...
    MOC_EAST_ERR  (TIME) float64 ...
    ...            ...
    MHT_EAST      (TIME) float64 ...
    MHT_EAST_ERR  (TIME) float64 ...
    MST_EAST      (TIME) float64 ...
    MST_EAST_ERR  (TIME) float64 ...
    FW_EAST       (TIME) float64 ...
    FW_EAST_ERR   (TIME) float64 ...
Attributes: (12/17)
    title:                     OSNAP MOC time series (2014-2018)
    project:                   OSNAP
    contributor_name:          F. Li; M.S. Lozier; S. Bacon; A. Bower; S.A. C...
    contributor_role:          data design, collection and/or processing
    contributor_institution:   Georgia Institute of Technology, USA; National...
    publisher_name:            M. Susan Lozier; Feili Li
    ...                        ...
    date_created:              2021-11-19 18:09:07
    time_coverage_start:       2014-07-31 00:00:00
    time_coverage_end:         2018-06-09 00:00:00
    time_coverage_duratinon:   P630D
    time_coverage_resolution:  PT30D
    netcdf_version:            4.3

Define the OSNAP lines

We define the OSNAP lines roughly by point pairs in [longitude, latitude] space. The lines are then computed via the function ecco_v4_py.get_section_line_masks as the great circle arc between these two points.

See below for similar MATLAB functions

pt1_east = [-44, 60]
pt2_east = [-5, 56]

pt1_west = [-56, 51]
pt2_west = [-45, 60]
maskC_east, maskW_east, maskS_east = ecco.get_section_line_masks(pt1_east,pt2_east,ds)
maskC_west, maskW_west, maskS_west = ecco.get_section_line_masks(pt1_west,pt2_west,ds)
maskC_tot = (maskC_east+maskC_west).where(maskC_east+maskC_west==1,0)
maskW_tot = (maskW_east+maskW_west).where(np.abs(maskW_east)+np.abs(maskW_west)==1,0)
maskS_tot = (maskS_east+maskS_west).where(np.abs(maskS_east)+np.abs(maskS_west)==1,0)

-179.875 179.875
-180.0 180.0
-89.875 89.875
-90.0 90.0
# use dx=.1, dy=.1 so that plot shows the osnap array as a thin
# line.  remember, plot_proj_to_latlon_grid first interpolates
# the model grid to lat-lon with grid spacing as dx, dy
P = ecco.plot_proj_to_latlon_grid(ds.XC, ds.YC, maskC_tot, \
P[1].set_extent([-80, 20, 40, 80], crs=ccrs.PlateCarree())
-179.95 179.95000000000002
-180.0 180.0
-89.95 89.95
-90.0 90.0

We have defined many commonly used sections in oceanography so that users can access these easily. The available sections are shown below. This allows one to, e.g. compute volumetric transport across the Drake Passage as follows:

drake_vol = ecco_v4_py.calc_section_vol_trsp(ds,section_name='Drake Passage')

Similarly, we can do the same with calc_section_heat_trsp and calc_section_salt_trsp.

One can also get these pre-defined section masks as follows:

pt1,pt2 = ecco_v4_py.get_section_endpoints('Drake Passage')
maskC, maskW, maskS = ecco_v4_py.get_section_line_masks(ds,pt1,pt2)

Finally, one can see similar functions in MATLAB:

['Bering Strait',
 'Florida Strait',
 'Florida Strait W1',
 'Florida Strait S1',
 'Florida Strait E1',
 'Florida Strait E2',
 'Florida Strait E3',
 'Florida Strait E4',
 'Davis Strait',
 'Denmark Strait',
 'Iceland Faroe',
 'Faroe Scotland',
 'Scotland Norway',
 'Fram Strait',
 'Barents Sea',
 'Labrador Greenland',
 'Hudson Strait',
 'English Channel',
 'Newfoundland Iberia',
 'Drake Passage',
 'Indonesia W1',
 'Indonesia W2',
 'Indonesia W3',
 'Indonesia W4',
 'Australia Antarctica',
 'Madagascar Channel',
 'Madagascar Antarctica',
 'South Africa Antarctica']

Compute the overturning streamfunction in depth space

The function calc_section_stf computes the overturning streamfunction across the plane normal to the section denoted by the west and south masks. It is also possible to compute the overturning streamfunction at a particular latitude band, as is done to compare to the RAPID array, for instance. Please see the function calc_meridional_stf to do this, which is also in ecco_v4_py.calc_stf.

Note that we can also compute the volumetric, heat, or salt transport across these sections using the first three functions defined in ecco_v4_py.calc_section_trsp.

In MATLAB, we can compute meridional overturning streamfunctions with gcmfaces_calc/calc_overturn.m. Section transports can be computed with gcmfaces_calc/calc_transports.m once the masks are defined.

osnap_z_stf_east = ecco.calc_section_stf(ds,\
                                         pt1=pt1_east, \
                                         section_name='OSNAP East Overturning Streamfunction').compute()

osnap_z_stf_west = ecco.calc_section_stf(ds, \
                                         pt1=pt1_west, \
                                         section_name='OSNAP West Overturning Streamfunction').compute()

osnap_z_stf_tot = ecco.calc_section_stf(ds,\
                                        maskW=maskW_tot, \
                                        section_name='OSNAP Total Overturning Streamfunction').compute()
Wall time: 2min 51s
def osnap_depth_stf_vs_time(stf_ds,label):
    fig = plt.figure(figsize=(18,6))

    # Time evolving
    time_edge_extrap = np.hstack((stf_ds['time'].values[0] - (0.5*np.diff(stf_ds['time'].values[0:2])),\
                                  stf_ds['time'].values[:-1] + (0.5*np.diff(stf_ds['time'].values)),\
                                  stf_ds['time'].values[-1] + (0.5*np.diff(stf_ds['time'].values[-2:]))))
    Z_edge_extrap = np.hstack((np.array([0]),\
                               stf_ds['Z'].values[:-1] + (0.5*np.diff(stf_ds['Z'].values)),\
    plt.title('ECCOv4r4\nOverturning streamfunction across OSNAP %s [Sv]' % label)
    plt.ylabel('Depth [m]')
    cb = plt.colorbar()

    plt.title('ECCOv4r4\nTime mean streamfunction, OSNAP %s' % label)
    plt.ylabel('Depth [m]')
osnap_z_ov_east = osnap_z_stf_east['psi_moc'].max(dim='k')
osnap_z_ov_west = osnap_z_stf_west['psi_moc'].max(dim='k')
osnap_z_ov_tot = osnap_z_stf_tot['psi_moc'].max(dim='k')
fig = plt.figure(figsize=(12,6))
plt.title('ECCOv4r4\nOverturning strength across OSNAP lines [Sv]')
plt.legend((('East, avg. = %.2f $\pm$ %.2f [Sv]' %
            ('West, avg. = %.2f $\pm$ %.2f [Sv]' %
            ('Total, avg. = %.2f $\pm$ %.2f [Sv]' %

Compare to Observations


var_list = ['MOC_EAST','MOC_WEST','MOC_ALL']
obs_handles = []
for var, err in zip(var_list,err_list):
    if 'MOC_ALL' in var:
        clr = 'k'
    elif 'MOC_EAST' in var:
        clr = 'r'
    elif 'MOC_WEST' in var:
        clr = 'b'

    curr_hand = plt.plot(obs['TIME'],obs[var], color=clr)


# plot ECCO v4r4 equivalent
ecco_hand0 = plt.plot(osnap_z_ov_east['time'],osnap_z_ov_east, 'r--')
ecco_hand1 = plt.plot(osnap_z_ov_west['time'],osnap_z_ov_west, 'b--',)
ecco_hand2 = plt.plot(osnap_z_ov_tot['time'],osnap_z_ov_tot, 'k--')
ecco_handles = [ecco_hand0[0],ecco_hand1[0],ecco_hand2[0]]

plt.ylabel('[%s]' % obs['MOC_ALL'].units)
plt.title('OSNAP MOC Derived from observations [%s]' % obs['MOC_ALL'].units)
plt.legend(handles=(obs_handles + ecco_handles),labels=(\
    ('East MOC, (Lozier et al, 2019)\nTime mean: %.2f $\pm$ %.2f %s' %
    ('West MOC, (Lozier et al, 2019)\nTime mean: %.2f $\pm$ %.2f %s' %
    ('Total MOC, (Lozier et al, 2019)\nTime mean: %.2f $\pm$ %.2f %s' %
    ('East MOC, ECCOv4r4\nTime mean: %.2f $\pm$ %.2f %s' %
            (osnap_z_ov_east.mean('time'),osnap_z_ov_east.std('time'), \
    ('West MOC, ECCOv4r4\nTime mean: %.2f $\pm$ %.2f %s' %
            (osnap_z_ov_west.mean('time'),osnap_z_ov_west.std('time'), \
    ('Total MOC, ECCOv4r4\nTime mean: %.2f $\pm$ %.2f %s' %
            (osnap_z_ov_tot.mean('time'),osnap_z_ov_tot.std('time'), \
    loc='center left', bbox_to_anchor=(1, 0.5))
2023-05-02 21:25:35,339 - distributed.scheduler - WARNING - Worker failed to heartbeat within 300 seconds. Closing: <WorkerState 'tcp://', name: 2, status: running, memory: 0, processing: 0>
2023-05-02 21:25:35,379 - distributed.scheduler - WARNING - Worker failed to heartbeat within 300 seconds. Closing: <WorkerState 'tcp://', name: 3, status: running, memory: 0, processing: 0>
2023-05-02 21:25:35,389 - distributed.scheduler - WARNING - Received heartbeat from unregistered worker 'tcp://'.
2023-05-02 21:25:35,398 - distributed.scheduler - WARNING - Received heartbeat from unregistered worker 'tcp://'.
2023-05-02 21:25:39,424 - distributed.nanny - WARNING - Restarting worker
2023-05-02 21:25:39,459 - distributed.nanny - WARNING - Restarting worker
2023-05-03 19:30:38,689 - distributed.scheduler - WARNING - Worker failed to heartbeat within 300 seconds. Closing: <WorkerState 'tcp://', name: 0, status: running, memory: 0, processing: 0>
2023-05-03 19:30:38,770 - distributed.scheduler - WARNING - Worker failed to heartbeat within 300 seconds. Closing: <WorkerState 'tcp://', name: 3, status: running, memory: 0, processing: 0>
2023-05-03 19:30:38,786 - distributed.scheduler - WARNING - Received heartbeat from unregistered worker 'tcp://'.
2023-05-03 19:30:38,791 - distributed.scheduler - WARNING - Received heartbeat from unregistered worker 'tcp://'.
2023-05-03 19:30:43,501 - distributed.nanny - WARNING - Restarting worker
2023-05-03 19:30:47,019 - distributed.nanny - WARNING - Restarting worker