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, so coverage is densest Oct–Apr. Its
coverage starts 2023-10-18; we pin a mid-winter date inside that window
rather than computing one from today's date.
Reads credentials from COPERNICUSMARINE_SERVICE_USERNAME /
COPERNICUSMARINE_SERVICE_PASSWORD.
Setup¶
The imports: numpy for the array math, pyramids' NetCDF for reading the
downloaded L4 file and Dataset / GeoReference / ColorBar for the
reprojected plot (plotting goes entirely through polar.plot(), not
matplotlib directly), and earthlens for the unified EarthLens entry
point and the CMEMS Catalog. The output directory is created up front.
import os
from pathlib import Path
import numpy as np
from pyramids.dataset import Dataset, GeoReference
from pyramids.netcdf import NetCDF
from pyramids.plot import ColorBar
from earthlens.cmems import Catalog
from earthlens.core import EarthLens
OUT_DIR = Path('data/cmems-seaice')
OUT_DIR.mkdir(parents=True, exist_ok=True)
Request parameters¶
The dataset id and the probe date are both fixed. A docs example should give the
same answer whichever day it runs — and "near-real-time" is not a promise that
the window keeps rolling. This product has not advanced past 2026-04-12, so
a date derived from datetime.now() walks off the end of the coverage and the
subset is rejected with CoordinatesOutOfDatasetBounds. Mid-February also sits
in the dense part of the altimetry retrieval season.
DATASET_ID = 'esa_obs-si_arc_phy-sit_nrt_l4-multi_P1D-m'
# Pinned, not computed from today: the NRT window has not advanced past
# 2026-04-12, so a `now - N days` probe falls outside the coverage
# (2023-10-18 onward) and the subset fails. Assumes this NRT product's
# retention only ever grows forward from 2023-10-18 -- if it were instead
# a rolling window, a pinned date could eventually fall out of the front
# of the window too, not just risk falling past the end.
PROBE_DATE = '2026-02-15'
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.
ds_meta = Catalog().get_dataset(DATASET_ID)
print(ds_meta)
print(f'domain: {ds_meta.domain}')
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.
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=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.
paths = el.download()
print(paths)
Open and inspect the field¶
The subset comes back on a regular geographic grid — latitude / longitude
in EPSG:4326, not the polar-stereographic xc / yc axes the native L4 product
is distributed on, because the CMEMS subset service reprojects. We read the file
through pyramids' NetCDF and close the handle when done.
nc = NetCDF.read_file(paths[0], read_only=True)
print('variables :', nc.variable_names)
print('dimensions:', nc.dimension_sizes)
read_array already collapses the single time step, so the thickness comes back
as a 2-D (rows, columns) masked array — there is no leading time axis to index
away. The variable is stored unpacked (scale_factor 1, add_offset 0), so the
values are metres as they are read.
sit_var = nc.get_variable('sea_ice_thickness')
sit = sit_var.read_array(masked=True).astype('float64')
print(f'grid {sit.shape[0]} x {sit.shape[1]}, epsg {sit_var.epsg}')
print(f'valid cells: {sit.count()} of {sit.size}')
print(
f'thickness min/mean/max: {float(sit.min()):.2f} / '
f'{float(sit.mean()):.2f} / {float(sit.max()):.2f} m'
)
Map the thickness field¶
The subset arrives on a lat/lon grid, which draws the Arctic basin as a 15:1 strip. Reprojecting to EPSG:3413 (NSIDC north polar stereographic) at the 25 km cell size conventional for sea ice turns it back into a map: thicker multi-year ice banked against the Canadian Arctic Archipelago and northern Greenland, thinner first-year ice out toward the marginal seas. The hole at the pole is the 88°N edge of the bounding box requested above, not missing data.
field = Dataset.from_array(
sit.filled(np.nan),
no_data_value=np.nan,
geo_ref=GeoReference(geo=sit_var.geotransform, epsg=sit_var.epsg),
)
# 25 km is the conventional sea-ice grid spacing, and it keeps the reprojected
# raster small enough to plot -- a native-resolution warp is ~60 million cells.
polar = field.to_crs(3413, cell_size=25000)
print(f'reprojected: {polar.rows} x {polar.columns}, epsg {polar.epsg}')
glyph = polar.plot(
cmap='Blues',
vmin=0,
vmax=5,
colorbar=ColorBar(label='sea-ice thickness (m)'),
title=f'Arctic sea-ice thickness (L4 NRT, {PROBE_DATE})',
)
glyph.ax.xaxis.set_ticks_position('bottom')
glyph.ax.xaxis.set_label_position('bottom')
glyph.ax.set_xlabel('x (m, EPSG:3413)')
glyph.ax.set_ylabel('y (m)')
glyph.ax.set_facecolor('#eef2f5')
nc.close()