MOM6/CESM Ocean Model Analysis

This notebook shows how to load and analyze ocean data from an out-of-the-box MOM6/CESM G-case simulation (coupled ocean ocean/sea ice).

NOTE: MOM6/CESM is not ready to be used for research.

In [ ]:
%matplotlib inline

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import holoviews as hv
import datashader
from holoviews.operation.datashader import regrid, shade, datashade

hv.extension('bokeh', width=100)

Create and Connect to Dask Distributed Cluster

This will launch a cluster of virtual machines in the cloud.

In [ ]:
from dask.distributed import Client, progress
from dask_gateway import Gateway

gateway = Gateway()
cluster = gateway.new_cluster()

👆 Don't forget to click this link to get the cluster dashboard

In [ ]:
client = Client(cluster)

Load MOM6/CESM Data

This data is stored in xarray-zarr format in Google Cloud Storage. This format is optimized for parallel distributed reads from within the cloud environment.

It may take up to a minute to initialize the dataset when you run this cell.

In [ ]:
from intake import open_catalog

cat = open_catalog("")
In [ ]:
ds = cat.ocean.cesm_mom6_example.to_dask()

Visualize SST Data with Holoviews and Datashader

The cells below show how to interactively explore the dataset.

In [ ]:
sst_ds = hv.Dataset(ds['SST'], kdims=['time', 'geolon', 'geolat'])
sst =, kdims=["geolon", "geolat"], dynamic=True)
%opts RGB [width=900 height=600] 
datashade(sst, precompute=True,

Visualize SSS Data with Holoviews and Datashader

In [ ]:
sss_ds = hv.Dataset(ds['SSS'], kdims=['time', 'geolon', 'geolat'])
sss =, kdims=["geolon", "geolat"], dynamic=True)
%opts RGB [width=900 height=600] 
datashade(sss, precompute=True,

Data reduction

Here we make a data reduction by taking the time of SST and SSS. This demonstrates how the cluster distributes the reads from storage.

In [ ]:
SST_mean = ds.SST.mean(dim=('time'))
In [ ]:
SSS_mean = ds.SSS.mean(dim=('time'))
In [ ]:
%time SST_mean.load()
In [ ]:
# plot mean SST
qm = hv.QuadMesh((ds.geolon.values, ds.geolat.values, SST_mean))
datashade(qm, precompute=True,
In [ ]:
%time SSS_mean.load()
In [ ]:
# plot mean SSS
qm = hv.QuadMesh((ds.geolon.values, ds.geolat.values, SSS_mean))
datashade(qm, precompute=True,