Using THREDDS/Opendap to subset CCI Sea State L4 dataset#

A THREDDS Data Server (TDS, https://www.unidata.ucar.edu/software/tds) provides catalog, metadata, and data access services for scientific data. Every TDS publishes THREDDS catalogs that advertise the datasets and services it makes available. The available remote data access protocols include OPeNDAP, OGC WCS, OGC WMS, and HTTP.

The TDS also supports several dataset collection services including some sophisticated dataset aggregation capabilities. This allows the TDS to aggregate a collection of datasets into a single virtual dataset, greatly simplifying user access to that data collection.

In CCI Sea State, we provide access through a THREDDS server to the Altimeter L4 dataset (gridded monthly statistics of SWH). All monthly files are nicely aggregated by THREDDS into a 3-dimensional data array (time, lat, lon) that can be remotely subsetted and read through OpenDAP API.

The end-point to read remotely the CCI Sea State L4 dataset is: https://tds0.ifremer.fr/thredds/dodsC/ESACCI-SEASTATE-L4-SWH-MULTI_1M

In this tutorial, we give you an example on how to use it from python (but libraries are available in other languages for OpenDAP).

Let’s first import xarray, this is the only package you need to read remotely data through OPENDAP:

import xarray as xr

Open the OpenDAP end-point as you would open a file on local storage. We see the result as a 3-dimensional Dataset object:

THREDDS convert datetime information to strings.We can convert them back to numpy datetime64 data type:

import numpy as np
ds['time'] = ds['time'].astype(np.datetime64)

/home/jfpiolle/mambaforge/envs/cciseastate/lib/python3.10/site-packages/xarray/core/duck_array_ops.py:226: UserWarning: no explicit representation of timezones available for np.datetime64
  return data.astype(dtype, **kwargs)
/tmp/ipykernel_12784/2712952193.py:2: UserWarning: Converting non-nanosecond precision datetime values to nanosecond precision. This behavior can eventually be relaxed in xarray, as it is an artifact from pandas which is now beginning to support non-nanosecond precision values. This warning is caused by passing non-nanosecond np.datetime64 or np.timedelta64 values to the DataArray or Variable constructor; it can be silenced by converting the values to nanosecond precision ahead of time.
  ds['time'] = ds['time'].astype(np.datetime64)

Let’s now read a full time series of the mean significant wave height (swh_mean variable) at a given point (longitude=-4.5, latitude=46.5). Note that you must provide latitude and longitude that exist in the dataset coordinates (it won’t approximate to the closest grid point but raise an error instead).

This is a lazy operation so compute() must be executed to actully donwload the data.

We can now display the returned time series:

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(10,3))
swh.plot()

[<matplotlib.lines.Line2D at 0x7fa0d2ff8eb0>]

../_images/4be258e279dc5fc1ab55f8692252a1c4e512355c36c13df4613cb58a572f16e7.png

To read a subset over an area instead of at a single point, just use slices in the selection call:

%%time
swh_area = ds.sel(lon=slice(-6.5, -3.5), lat=slice(44.5, 48.5)).swh_mean.compute()
swh_area.shape

CPU times: user 11 ms, sys: 310 μs, total: 11.3 ms
Wall time: 3.28 s

(385, 5, 4)

We can no display a time series of the mean value of swh_mean over the selected area (swh.mean(dim=['lon','lat'])):

fig = plt.figure(figsize=(10,3))

swh_area.mean(dim=['lon','lat']).plot()

[<matplotlib.lines.Line2D at 0x7fa0c88ec820>]

../_images/b207f9ca2425714a607c9b4446480c3ab0638b4b7b56563f3324784cadd23398.png