GLORYS12 — thermocline temperature time-series¶

Demonstrates the depth-axis pattern in the CMEMS backend: pull one year of daily GLORYS12 (thetao, potential temperature) at three fixed depths and plot the time-series at one ocean point.

The CMEMS toolbox returns thetao as a 4-D (time, depth, lat, lon) NetCDF; the minimum_depth / maximum_depth kwargs let you clip server-side so you only pay for the levels you want. This notebook downloads the surface-to-500 m slab once and slices three target depths client-side.

Reads credentials from COPERNICUSMARINE_SERVICE_USERNAME / COPERNICUSMARINE_SERVICE_PASSWORD (see Authentication).

Setup¶

Consolidate the imports up front. pyramids provides NetCDF (reading the downloaded slab); earthlens provides the unified EarthLens entry point and the CMEMS Catalog. xarray and matplotlib are used later for the depth-axis slicing and the plot.

In [ ]:

import os
from pathlib import Path

import matplotlib.pyplot as plt
import xarray as xr
from pyramids.netcdf import NetCDF

from earthlens import EarthLens
from earthlens.cmems import Catalog

import os from pathlib import Path import matplotlib.pyplot as plt import xarray as xr from pyramids.netcdf import NetCDF from earthlens import EarthLens from earthlens.cmems import Catalog

Request parameters¶

The output directory, the GLORYS12 dataset id, the three target depths, and the ocean point (mid-Atlantic, off the Azores) with a 1° × 1° box around it. These constants drive every cell below.

In [ ]:

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

DATASET_ID = 'cmems_mod_glo_phy_my_0.083deg_P1D-m'
TARGET_DEPTHS_M = (20.0, 100.0, 500.0)
POINT_LAT, POINT_LON = 35.0, -25.0  # mid-Atlantic, off the Azores
BBOX = dict(lat_lim=[34.5, 35.5], lon_lim=[-25.5, -24.5])

OUT_DIR = Path('data/cmems-glorys') OUT_DIR.mkdir(parents=True, exist_ok=True) DATASET_ID = 'cmems_mod_glo_phy_my_0.083deg_P1D-m' TARGET_DEPTHS_M = (20.0, 100.0, 500.0) POINT_LAT, POINT_LON = 35.0, -25.0 # mid-Atlantic, off the Azores BBOX = dict(lat_lim=[34.5, 35.5], lon_lim=[-25.5, -24.5])

Inspect the catalog entry¶

Before downloading, look the dataset up in the CMEMS Catalog to confirm its cadence, domain, and the units of thetao.

In [ ]:

ds = Catalog().get_dataset(DATASET_ID)
print(f'{DATASET_ID}: cadence={ds.cadence}, domain={ds.domain}')
print('thetao units:', ds.variables['thetao'].units)

ds = Catalog().get_dataset(DATASET_ID) print(f'{DATASET_ID}: cadence={ds.cadence}, domain={ds.domain}') print('thetao units:', ds.variables['thetao'].units)

Download — one year of daily thetao, surface to 500 m¶

1° × 1° box, 366 days, depths 0-500 m. The toolbox returns a NetCDF a few MB in size.

Build the EarthLens request first — source, date window, cadence, dataset, variable, the bounding box, the output path, the depth clip, and the CMEMS credentials from the environment.

In [ ]:

earthlens = EarthLens(
    data_source='cmems',
    start='2020-01-01',
    end='2020-12-31',
    cadence='daily',
    dataset=DATASET_ID,
    variables=['thetao'],
    **BBOX,
    path=str(OUT_DIR),
    minimum_depth=0.0,
    maximum_depth=500.0,
    service_username=os.environ.get('COPERNICUSMARINE_SERVICE_USERNAME'),
    service_password=os.environ.get('COPERNICUSMARINE_SERVICE_PASSWORD'),
)

earthlens = EarthLens( data_source='cmems', start='2020-01-01', end='2020-12-31', cadence='daily', dataset=DATASET_ID, variables=['thetao'], **BBOX, path=str(OUT_DIR), minimum_depth=0.0, maximum_depth=500.0, service_username=os.environ.get('COPERNICUSMARINE_SERVICE_USERNAME'), service_password=os.environ.get('COPERNICUSMARINE_SERVICE_PASSWORD'), )

Run the download as its own step; it returns the list of written file paths.

In [ ]:

paths = earthlens.download()
print(paths)

paths = earthlens.download() print(paths)

Slice three depths at one ocean point¶

Open the returned NetCDF with pyramids.netcdf.NetCDF, then index the depth axis to the three target levels and the lat/lon axes to the central pixel.

Open the slab as labelled xarray¶

Read the NetCDF and convert it to an xarray dataset. decode_cf turns the CF "hours since 1950" time axis into datetime64, giving labelled (time, depth, latitude, longitude) axes.

In [ ]:

nc = NetCDF.read_file(str(paths[0]), read_only=True)
# decode_cf turns the CF "hours since 1950" time axis into datetime64.
ds = xr.decode_cf(nc.to_xarray())  # labelled (time, depth, latitude, longitude)
nc.close()
print('dims:', dict(ds.sizes))
print('depth levels in slab:', [round(float(d), 1) for d in ds['depth'].values])

nc = NetCDF.read_file(str(paths[0]), read_only=True) # decode_cf turns the CF "hours since 1950" time axis into datetime64. ds = xr.decode_cf(nc.to_xarray()) # labelled (time, depth, latitude, longitude) nc.close() print('dims:', dict(ds.sizes)) print('depth levels in slab:', [round(float(d), 1) for d in ds['depth'].values])

Select the point and the three nearest levels¶

Take the nearest grid point to the target location, then the nearest model level to each target depth — xarray label-based selection handles the snapping.

In [ ]:

# Nearest grid point to the target location, nearest model level to each
# target depth — xarray label-based selection handles the snapping.
point = ds['thetao'].sel(latitude=POINT_LAT, longitude=POINT_LON, method='nearest')
series = {d_m: point.sel(depth=d_m, method='nearest') for d_m in TARGET_DEPTHS_M}
for d_m, da in series.items():
    print(
        f'{int(d_m):>4d} m -> nearest level {float(da["depth"]):.1f} m, '
        f'{da.size} daily values'
    )

# Nearest grid point to the target location, nearest model level to each # target depth — xarray label-based selection handles the snapping. point = ds['thetao'].sel(latitude=POINT_LAT, longitude=POINT_LON, method='nearest') series = {d_m: point.sel(depth=d_m, method='nearest') for d_m in TARGET_DEPTHS_M} for d_m, da in series.items(): print( f'{int(d_m):>4d} m -> nearest level {float(da["depth"]):.1f} m, ' f'{da.size} daily values' )

Plot the seasonal cycle at three depths¶

Surface follows the seasonal cycle clearly; below the thermocline (here at 100 m and 500 m) the temperature is much smoother and the seasonal signal is damped.

In [ ]:

fig, ax = plt.subplots(figsize=(8, 4))
for d_m, da in series.items():
    ax.plot(da['time'].values, da.values, label=f'{int(d_m)} m')
ax.set_xlabel('Date (2020)')
ax.set_ylabel('Potential temperature (degrees_C)')
ax.set_title('GLORYS12 thetao at 35N, 25W (2020)')
ax.legend()
ax.grid(alpha=0.3)
fig.tight_layout()

fig, ax = plt.subplots(figsize=(8, 4)) for d_m, da in series.items(): ax.plot(da['time'].values, da.values, label=f'{int(d_m)} m') ax.set_xlabel('Date (2020)') ax.set_ylabel('Potential temperature (degrees_C)') ax.set_title('GLORYS12 thetao at 35N, 25W (2020)') ax.legend() ax.grid(alpha=0.3) fig.tight_layout()