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. pyramids 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 numpy as np
import pandas as pd
from pyramids.netcdf import NetCDF

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

import os from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from pyramids.netcdf import NetCDF from earthlens.cmems import Catalog from earthlens.core import EarthLens

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

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

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(ds)
print(f'domain: {ds.domain}')
print(ds.variables['thetao'])

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

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'],
    aoi=[-25.5, 34.5, -24.5, 35.5],  # 1-degree box on the Azores point (W, S, E, N)
    path=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'], aoi=[-25.5, 34.5, -24.5, 35.5], # 1-degree box on the Azores point (W, S, E, N) path=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 and read its dimensions¶

Read the NetCDF and inspect its variable and dimension names. The depth levels come from NetCDF.get_dimension_values('depth'); the (time, depth, latitude, longitude) cube itself is reshaped by hand from the band-major layout further down.

In [ ]:

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

# GLORYS12's 50-level grid, clipped here to the levels above 500 m.
# NetCDF.get_dimension_values reads any dimension's coordinate values directly
# (added by pyramids serapeum-org/pyramids#1133, present since the 0.64.0 floor).
DEPTH_LEVELS_M = np.asarray(nc.get_dimension_values('depth'))
print('depth levels in slab:', [round(float(x), 1) for x in DEPTH_LEVELS_M])

nc = NetCDF.read_file(paths[0], read_only=True) print('variables :', nc.variable_names) print('dimensions:', nc.dimension_sizes) # GLORYS12's 50-level grid, clipped here to the levels above 500 m. # NetCDF.get_dimension_values reads any dimension's coordinate values directly # (added by pyramids serapeum-org/pyramids#1133, present since the 0.64.0 floor). DEPTH_LEVELS_M = np.asarray(nc.get_dimension_values('depth')) print('depth levels in slab:', [round(float(x), 1) for x in DEPTH_LEVELS_M])

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 — nearest-level selection handles the snapping.

In [ ]:

# thetao is stored packed (int16 with a scale_factor and add_offset); NetCDF.read_array()
# unpacks CF-packed data to physical units by default as of pyramids 0.64.
thetao = nc.get_variable('thetao')

n_time = nc.dimension_sizes['time']
n_depth = nc.dimension_sizes['depth']
cube = thetao.read_array(masked=True).astype('float64')
cube = cube.reshape(n_time, n_depth, thetao.rows, thetao.columns)

row, col = thetao.rowcol(POINT_LON, POINT_LAT)  # nearest cell from the geotransform
dates = pd.to_datetime(nc.get_time_variable())
nc.close()

series = {}
for target in TARGET_DEPTHS_M:
    level = int(np.abs(DEPTH_LEVELS_M - target).argmin())
    series[target] = cube[:, level, row, col]
    print(
        f'{int(target):>4d} m -> nearest level {DEPTH_LEVELS_M[level]:.1f} m, '
        f'{series[target].count()} daily values'
    )

# thetao is stored packed (int16 with a scale_factor and add_offset); NetCDF.read_array() # unpacks CF-packed data to physical units by default as of pyramids 0.64. thetao = nc.get_variable('thetao') n_time = nc.dimension_sizes['time'] n_depth = nc.dimension_sizes['depth'] cube = thetao.read_array(masked=True).astype('float64') cube = cube.reshape(n_time, n_depth, thetao.rows, thetao.columns) row, col = thetao.rowcol(POINT_LON, POINT_LAT) # nearest cell from the geotransform dates = pd.to_datetime(nc.get_time_variable()) nc.close() series = {} for target in TARGET_DEPTHS_M: level = int(np.abs(DEPTH_LEVELS_M - target).argmin()) series[target] = cube[:, level, row, col] print( f'{int(target):>4d} m -> nearest level {DEPTH_LEVELS_M[level]:.1f} m, ' f'{series[target].count()} 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, values in series.items():
    ax.plot(dates, 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, values in series.items(): ax.plot(dates, 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()