Arctic sea-ice thickness — a polar L4 map¶

Demonstrates the surface-only polar-map pattern: pull one day of Arctic sea-ice thickness from the merged CryoSat-2 / SMOS L4 product and render it as a map over the Arctic basin.

Dataset esa_obs-si_arc_phy-sit_nrt_l4-multi_P1D-m is near-real-time and surface-only (no depth axis). Sea-ice thickness retrieval from altimetry is winter-weighted — coverage is densest Oct–Apr — and the NRT window rolls forward, so we probe a recent date and note that a deep-summer re-run may need an earlier (winter) date.

Reads credentials from COPERNICUSMARINE_SERVICE_USERNAME / COPERNICUSMARINE_SERVICE_PASSWORD.

Setup¶

The imports: numpy and matplotlib for the array math and the map, xarray plus pyramids' NetCDF for reading the downloaded L4 file, and earthlens for the unified EarthLens entry point and the CMEMS Catalog. The output directory is created up front.

In [ ]:

import datetime as dt
import os
from pathlib import Path

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

from earthlens import EarthLens
from earthlens.cmems import Catalog

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

import datetime as dt import os from pathlib import Path import matplotlib.pyplot as plt import numpy as np import xarray as xr from pyramids.netcdf import NetCDF from earthlens import EarthLens from earthlens.cmems import Catalog OUT_DIR = Path('data/cmems-seaice') OUT_DIR.mkdir(parents=True, exist_ok=True)

Request parameters¶

The dataset id is fixed; the probe date is computed ~45 days back to stay inside the rolling NRT window. In deep summer (when thickness coverage thins out) bump this to a winter date such as '2024-02-15'.

In [ ]:

DATASET_ID = 'esa_obs-si_arc_phy-sit_nrt_l4-multi_P1D-m'
# ~45 days back keeps us inside the rolling NRT window; in deep summer
# bump this to a winter date (e.g. '2024-02-15') where SIT coverage is dense.
PROBE_DATE = (dt.datetime.now(dt.timezone.utc) - dt.timedelta(days=45)).strftime(
    '%Y-%m-%d'
)

DATASET_ID = 'esa_obs-si_arc_phy-sit_nrt_l4-multi_P1D-m' # ~45 days back keeps us inside the rolling NRT window; in deep summer # bump this to a winter date (e.g. '2024-02-15') where SIT coverage is dense. PROBE_DATE = (dt.datetime.now(dt.timezone.utc) - dt.timedelta(days=45)).strftime( '%Y-%m-%d' )

Inspect the catalog entry¶

Before downloading, look up the dataset in the CMEMS Catalog to confirm its domain and cadence and to see which variables it exposes.

In [ ]:

ds_meta = Catalog().get_dataset(DATASET_ID)
print(f'{DATASET_ID}: domain={ds_meta.domain}, cadence={ds_meta.cadence}')
print('variables:', sorted(ds_meta.variables))
print('probe date:', PROBE_DATE)

ds_meta = Catalog().get_dataset(DATASET_ID) print(f'{DATASET_ID}: domain={ds_meta.domain}, cadence={ds_meta.cadence}') print('variables:', sorted(ds_meta.variables)) print('probe date:', PROBE_DATE)

Download one day over the Arctic basin¶

Build the EarthLens request first — source, date window, cadence, dataset, the single sea_ice_thickness variable, an Arctic-basin bounding box, the output path, and the CMEMS credentials from the environment.

In [ ]:

el = EarthLens(
    data_source='cmems',
    start=PROBE_DATE,
    end=PROBE_DATE,
    cadence='daily',
    dataset=DATASET_ID,
    variables=['sea_ice_thickness'],
    aoi=[-180.0, 65.0, 180.0, 88.0],
    path=str(OUT_DIR),
    service_username=os.environ.get('COPERNICUSMARINE_SERVICE_USERNAME'),
    service_password=os.environ.get('COPERNICUSMARINE_SERVICE_PASSWORD'),
)

el = EarthLens( data_source='cmems', start=PROBE_DATE, end=PROBE_DATE, cadence='daily', dataset=DATASET_ID, variables=['sea_ice_thickness'], aoi=[-180.0, 65.0, 180.0, 88.0], path=str(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 and returns the list of written NetCDF paths.

In [ ]:

paths = el.download()
print(paths)

paths = el.download() print(paths)

Open and inspect the field¶

The L4 grid is a polar stereographic projection, so the array carries xc / yc projected axes (plus 2-D lat/lon). We read the file through pyramids' NetCDF, decode it into an xarray dataset, and close the handle.

In [ ]:

nc = NetCDF.read_file(str(paths[0]), read_only=True)
ds = xr.decode_cf(nc.to_xarray())
nc.close()
print('dims:', dict(ds.sizes))

nc = NetCDF.read_file(str(paths[0]), read_only=True) ds = xr.decode_cf(nc.to_xarray()) nc.close() print('dims:', dict(ds.sizes))

Collapse the single time step to get the 2-D thickness field, then summarise the valid (non-NaN) cells.

In [ ]:

sit = ds['sea_ice_thickness'].isel(time=0)
valid = sit.values[~np.isnan(sit.values)]
print(f'valid cells: {valid.size}')
if valid.size:
    print(
        f'thickness min/mean/max: {valid.min():.2f} / {valid.mean():.2f} / {valid.max():.2f} m'
    )

sit = ds['sea_ice_thickness'].isel(time=0) valid = sit.values[~np.isnan(sit.values)] print(f'valid cells: {valid.size}') if valid.size: print( f'thickness min/mean/max: {valid.min():.2f} / {valid.mean():.2f} / {valid.max():.2f} m' )

Map the thickness field¶

A plain pcolormesh of the projected grid — thicker multi-year ice north of the Canadian Arctic Archipelago, thinner first-year ice toward the marginal seas.

In [ ]:

fig, ax = plt.subplots(figsize=(6, 5))
mesh = ax.pcolormesh(sit.values, cmap='Blues', vmin=0, vmax=5)
fig.colorbar(mesh, ax=ax, label='Sea-ice thickness (m)')
ax.set_title(f'Arctic sea-ice thickness (L4 NRT, {PROBE_DATE})')
ax.set_xlabel('grid x')
ax.set_ylabel('grid y')
ax.set_aspect('equal')
fig.tight_layout()

fig, ax = plt.subplots(figsize=(6, 5)) mesh = ax.pcolormesh(sit.values, cmap='Blues', vmin=0, vmax=5) fig.colorbar(mesh, ax=ax, label='Sea-ice thickness (m)') ax.set_title(f'Arctic sea-ice thickness (L4 NRT, {PROBE_DATE})') ax.set_xlabel('grid x') ax.set_ylabel('grid y') ax.set_aspect('equal') fig.tight_layout()