Global waves — significant wave height¶

Demonstrates the 3-hourly wave-field pattern: pull a few days of the global wave reanalysis (MFWAM) significant wave height VHM0 over the North Atlantic, then show both a snapshot map and a point time-series through a winter storm.

Dataset cmems_mod_glo_wav_my_0.2deg_PT3H-i is a multi-year reanalysis (stable historical coverage) on a 0.2° grid with a 3-hourly step, so a fixed date works fine here.

Reads credentials from COPERNICUSMARINE_SERVICE_USERNAME / COPERNICUSMARINE_SERVICE_PASSWORD.

Setup¶

Consolidate the imports up front: pyramids' NetCDF for reading the downloaded field, Dataset for the snapshot map, and the Basemap / Feature / Contour / ColorBar style objects for the map's coastlines, banding and legend; cleopatra's LineGlyph for the time-series; plus the unified EarthLens entry point and the CMEMS Catalog. matplotlib is imported for the two-panel layout only — every mark on both panels is drawn by pyramids or cleopatra. We also pick an output directory for the download.

In [ ]:

import os
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from cleopatra.glyphs.primitives.line_glyph import LineGlyph
from pyramids.dataset import Dataset, GeoReference
from pyramids.netcdf import NetCDF
from pyramids.plot import Basemap, ColorBar, Contour, Feature

from earthlens.cmems import Catalog
from earthlens.core import EarthLens

OUT_DIR = Path('data/cmems-wave')
OUT_DIR.mkdir(parents=True, exist_ok=True)

import os from pathlib import Path import matplotlib.pyplot as plt import numpy as np from cleopatra.glyphs.primitives.line_glyph import LineGlyph from pyramids.dataset import Dataset, GeoReference from pyramids.netcdf import NetCDF from pyramids.plot import Basemap, ColorBar, Contour, Feature from earthlens.cmems import Catalog from earthlens.core import EarthLens OUT_DIR = Path('data/cmems-wave') OUT_DIR.mkdir(parents=True, exist_ok=True)

Inspect the dataset in the catalog¶

Before downloading, look the dataset up in the CMEMS Catalog to confirm its domain, cadence, and the units of the VHM0 variable. We also fix the point (NE Atlantic, west of Brittany) used later for the time-series.

In [ ]:

DATASET_ID = 'cmems_mod_glo_wav_my_0.2deg_PT3H-i'
POINT_LAT, POINT_LON = 48.0, -16.0  # NE Atlantic, west of Brittany

ds_meta = Catalog().get_dataset(DATASET_ID)
print(ds_meta)
print(f'domain: {ds_meta.domain}')
print(ds_meta.variables['VHM0'])

DATASET_ID = 'cmems_mod_glo_wav_my_0.2deg_PT3H-i' POINT_LAT, POINT_LON = 48.0, -16.0 # NE Atlantic, west of Brittany ds_meta = Catalog().get_dataset(DATASET_ID) print(ds_meta) print(f'domain: {ds_meta.domain}') print(ds_meta.variables['VHM0'])

Download the storm window over the North Atlantic¶

3-hourly over three days = 25 time steps of a 2-D field on a ~25°×15° box. The window is centred on the storm rather than opening at its peak: at the probe point Hs climbs from ~4.5 m on 02-07, tops out on 02-08, and decays through 02-09, so the series covers the whole passage. We build the request first — source, date window, cadence, dataset, variable, bounding box, output path, and the Copernicus Marine credentials.

In [ ]:

el = EarthLens(
    data_source='cmems',
    start='2014-02-07',
    end='2014-02-10',
    cadence='hourly',
    dataset=DATASET_ID,
    variables=['VHM0'],
    aoi=[-30.0, 40.0, -5.0, 55.0],
    path=OUT_DIR,
    service_username=os.environ.get('COPERNICUSMARINE_SERVICE_USERNAME'),
    service_password=os.environ.get('COPERNICUSMARINE_SERVICE_PASSWORD'),
)

el = EarthLens( data_source='cmems', start='2014-02-07', end='2014-02-10', cadence='hourly', dataset=DATASET_ID, variables=['VHM0'], aoi=[-30.0, 40.0, -5.0, 55.0], path=OUT_DIR, service_username=os.environ.get('COPERNICUSMARINE_SERVICE_USERNAME'), service_password=os.environ.get('COPERNICUSMARINE_SERVICE_PASSWORD'), )

With the request built, download() fetches the subset to the output directory and returns the list of written file paths.

In [ ]:

paths = el.download()
print(paths)

paths = el.download() print(paths)

Open the field and pull a point series¶

Read the downloaded NetCDF with pyramids' NetCDF, which decodes the CF metadata itself, and report the variables and dimensions.

In [ ]:

nc = NetCDF.read_file(paths[0], read_only=True)
print('variables :', nc.variable_names)
print('dimensions:', nc.dimension_sizes)

nc = NetCDF.read_file(paths[0], read_only=True) print('variables :', nc.variable_names) print('dimensions:', nc.dimension_sizes)

Significant wave height at the chosen point¶

Select the VHM0 field and pull the nearest grid point to (48N, 16W), then report the peak significant wave height as the storm passes.

In [ ]:

vhm0 = nc.get_variable('VHM0')

# get_time_variable defaults to a date-only format, which collapses all eight
# steps of a day onto one label and throws away the 3-hourly cadence this
# notebook is about, so ask for the hour and build a real datetime axis.
stamps = nc.get_time_variable(time_format='%Y-%m-%d %H:%M')
times = np.array(stamps, dtype='datetime64[m]')

# VHM0 is stored packed (int16 with a scale_factor); NetCDF.read_array() unpacks
# CF-packed data to physical units by default as of pyramids 0.64.
cube = vhm0.read_array(masked=True).astype('float64')

row, col = vhm0.rowcol(
    POINT_LON, POINT_LAT
)  # nearest cell, from the variable's geotransform
point = cube[:, row, col]
peak_step = int(np.ma.argmax(point))
step_hours = int((times[1] - times[0]) / np.timedelta64(1, 'h'))
print(f'{len(times)} steps, {step_hours} h apart')
print(f'peak Hs at ({POINT_LAT}N, {abs(POINT_LON)}W): {float(point.max()):.1f} m')
print(f'peak step : {peak_step} ({stamps[peak_step]} UTC)')

vhm0 = nc.get_variable('VHM0') # get_time_variable defaults to a date-only format, which collapses all eight # steps of a day onto one label and throws away the 3-hourly cadence this # notebook is about, so ask for the hour and build a real datetime axis. stamps = nc.get_time_variable(time_format='%Y-%m-%d %H:%M') times = np.array(stamps, dtype='datetime64[m]') # VHM0 is stored packed (int16 with a scale_factor); NetCDF.read_array() unpacks # CF-packed data to physical units by default as of pyramids 0.64. cube = vhm0.read_array(masked=True).astype('float64') row, col = vhm0.rowcol( POINT_LON, POINT_LAT ) # nearest cell, from the variable's geotransform point = cube[:, row, col] peak_step = int(np.ma.argmax(point)) step_hours = int((times[1] - times[0]) / np.timedelta64(1, 'h')) print(f'{len(times)} steps, {step_hours} h apart') print(f'peak Hs at ({POINT_LAT}N, {abs(POINT_LON)}W): {float(point.max()):.1f} m') print(f'peak step : {peak_step} ({stamps[peak_step]} UTC)')

Snapshot map + point time-series¶

Left: the wave field at the stormiest step, with Natural Earth land and coastlines under it so the domain reads as the North Atlantic, and the values banded at 1 m so the storm's structure is legible. Right: the 3-hourly Hs series at the probe point through the whole passage. First slice the snapshot out of the cube at the peak step.

In [ ]:

snapshot = Dataset.from_array(
    cube[peak_step].filled(np.nan),
    no_data_value=np.nan,
    geo_ref=GeoReference(geo=vhm0.geotransform, epsg=vhm0.epsg),
)
print(f'snapshot {snapshot.rows} x {snapshot.columns}, epsg {snapshot.epsg}')

snapshot = Dataset.from_array( cube[peak_step].filled(np.nan), no_data_value=np.nan, geo_ref=GeoReference(geo=vhm0.geotransform, epsg=vhm0.epsg), ) print(f'snapshot {snapshot.rows} x {snapshot.columns}, epsg {snapshot.epsg}')

Draw the two panels¶

matplotlib supplies the 1×2 grid and nothing else. pyramids draws the map into the left axes and cleopatra's LineGlyph draws the series into the right one; the probe point is a red star on the map and a dashed rule at the snapshot time on the series.

In [ ]:

fig, axes = plt.subplots(1, 2, figsize=(15, 5.0), width_ratios=[1.35, 1])
fig.subplots_adjust(left=0.05, right=0.98, wspace=0.32, bottom=0.16, top=0.90)

map_glyph = snapshot.plot(
    fig=fig,
    ax=axes[0],
    cmap='viridis',
    contour=Contour(levels=list(np.arange(0.0, 15.1, 1.0))),
    basemap=Basemap(
        relief=False,
        features=[
            Feature('land', facecolor='#e4ded2', edgecolor='none'),
            Feature('coastline', edgecolor='#3a3a3a', linewidth=0.6),
        ],
    ),
    colorbar=ColorBar(label='significant wave height (m)'),
    title=f'MFWAM VHM0 — {stamps[peak_step]} UTC',
)
map_glyph.ax.plot(
    [POINT_LON],
    [POINT_LAT],
    'r*',
    markersize=15,
    markeredgecolor='white',
    markeredgewidth=0.9,
)
map_glyph.ax.xaxis.set_ticks_position('bottom')
map_glyph.ax.xaxis.set_label_position('bottom')
map_glyph.ax.set_xlabel('longitude')
map_glyph.ax.set_ylabel('latitude')

series = LineGlyph(
    times,
    point.filled(np.nan),
    fig=fig,
    ax=axes[1],
    marker='o',
    line_width=1.8,
    color_1='#1f4e79',
)
_, series_ax, _ = series.line(title=f'Hs at {POINT_LAT:.0f}N, {abs(POINT_LON):.0f}W')
# line() applies only the title; the axis labels and grid it declares as options
# are ignored (serapeum-org/cleopatra#347), so they are set on the returned axes.
series_ax.set_xlabel('time (UTC)')
series_ax.set_ylabel('significant wave height (m)')
series_ax.grid(alpha=0.3)
series_ax.axvline(times[peak_step], color='firebrick', lw=1.0, ls='--')
series_ax.tick_params(axis='x', rotation=30)
for label in series_ax.get_xticklabels():
    label.set_horizontalalignment('right')

nc.close()

fig, axes = plt.subplots(1, 2, figsize=(15, 5.0), width_ratios=[1.35, 1]) fig.subplots_adjust(left=0.05, right=0.98, wspace=0.32, bottom=0.16, top=0.90) map_glyph = snapshot.plot( fig=fig, ax=axes[0], cmap='viridis', contour=Contour(levels=list(np.arange(0.0, 15.1, 1.0))), basemap=Basemap( relief=False, features=[ Feature('land', facecolor='#e4ded2', edgecolor='none'), Feature('coastline', edgecolor='#3a3a3a', linewidth=0.6), ], ), colorbar=ColorBar(label='significant wave height (m)'), title=f'MFWAM VHM0 — {stamps[peak_step]} UTC', ) map_glyph.ax.plot( [POINT_LON], [POINT_LAT], 'r*', markersize=15, markeredgecolor='white', markeredgewidth=0.9, ) map_glyph.ax.xaxis.set_ticks_position('bottom') map_glyph.ax.xaxis.set_label_position('bottom') map_glyph.ax.set_xlabel('longitude') map_glyph.ax.set_ylabel('latitude') series = LineGlyph( times, point.filled(np.nan), fig=fig, ax=axes[1], marker='o', line_width=1.8, color_1='#1f4e79', ) _, series_ax, _ = series.line(title=f'Hs at {POINT_LAT:.0f}N, {abs(POINT_LON):.0f}W') # line() applies only the title; the axis labels and grid it declares as options # are ignored (serapeum-org/cleopatra#347), so they are set on the returned axes. series_ax.set_xlabel('time (UTC)') series_ax.set_ylabel('significant wave height (m)') series_ax.grid(alpha=0.3) series_ax.axvline(times[peak_step], color='firebrick', lw=1.0, ls='--') series_ax.tick_params(axis='x', rotation=30) for label in series_ax.get_xticklabels(): label.set_horizontalalignment('right') nc.close()