Sea Surface Altimetry Data Analysis

For this example we will use gridded sea-surface altimetry data from The Copernicus Marine Environment:

http://marine.copernicus.eu/services-portfolio/access-to-products/?option=com_csw&view=details&product_id=SEALEVEL_GLO_PHY_L4_REP_OBSERVATIONS_008_047

This is a widely used dataset in physical oceanography and climate.

globe image

The dataset has already been extracted from copernicus and stored in google cloud storage in xarray-zarr format.

In [ ]:
import numpy as np
import xarray as xr
import holoviews as hv
import hvplot.xarray
import hvplot.pandas

Initialize Dataset

Here we load the dataset from the zarr store. Note that this very large dataset initializes nearly instantly, and we can see the full list of variables and coordinates.

In [ ]:
from intake import open_catalog

cat = open_catalog("https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/master.yaml")
In [ ]:
ds = cat.ocean.sea_surface_height.to_dask()
ds

Examine Metadata

For those unfamiliar with this dataset, the variable metadata is very helpful for understanding what the variables actually represent

In [ ]:
for v in ds.data_vars:
    print('{:>10}: {}'.format(v, ds[v].attrs['long_name']))

Create and Connect to Dask Distributed Cluster

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

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

☝️ Don't forget to click the link above to view the scheduler dashboard!

In [ ]:
client = Client(cluster)
client

Visually Examine Some of the Data

Let's do a sanity check that the data looks reasonable:

In [ ]:
ds.sla.sel(time='1982-08-07', method='nearest').hvplot(colormap='RdBu_r', width=900, height=550, rasterize=True)

Timeseries of Global Mean Sea Level

Here we make a simple yet fundamental calculation: the rate of increase of global mean sea level over the observational period.

In [ ]:
# the number of GB involved in the reduction
ds.sla.nbytes/1e9
In [ ]:
# the computationally intensive step
sla_timeseries = ds.sla.mean(dim=('latitude', 'longitude')).load()
In [ ]:
sla_full = sla_timeseries.hvplot(label='full data', grid=True,
                          title='Global Sea Level Rise', 
                          width=800, height=400)
sla_filt = sla_timeseries.rolling(time=365, center=True).mean().hvplot(label='rolling annual mean')
hv.Overlay([sla_full, sla_filt]).options(legend_position='top_left')

In order to understand how the sea level rise is distributed in latitude, we can make a sort of Hovmöller diagram.

In [ ]:
sla_hov = ds.sla.mean(dim='longitude').load()
In [ ]:
sla_hov.transpose().hvplot(rasterize=True, colormap='RdBu_r', width=900, height=400, clim=(-.2,.2))

We can see that most sea level rise is actually in the Southern Hemisphere.

Sea Level Variability

We can examine the natural variability in sea level by looking at its standard deviation in time.

In [ ]:
sla_std = ds.sla.std(dim='time').load()
sla_std.name = 'Sea Level Variability [m]'
In [ ]:
sla_std.hvplot(colormap='viridis', width=900, height=550, rasterize=True)