Oceanography — North Atlantic SST and sea-ice cover¶

Pull a year of monthly sea-surface temperature and sea-ice cover for a North Atlantic box. SST is a state variable (instantaneous K samples); sea-ice cover is a fractional area (0–1, also state). Both use op="auto" → mean for monthly aggregation.

Domain context. North Atlantic SST and sea-ice extent are the classical inputs to studies of the AMOC, NAO teleconnections, and Arctic–subarctic climate. ERA5 single-levels carries both, so one retrieve gives you a coherent monthly time-series ready for visual inspection.

Step 1 — catalog inspection¶

Both variables live on reanalysis-era5-single-levels. Note is_flux=False for both — the right reduction is the time-window mean.

In [1]:

from earthlens.ecmwf import Catalog

cat = Catalog()
for code in ("sea-surface-temperature", "sea-ice-cover"):
    spec = cat.get_variable("reanalysis-era5-single-levels", code)
    print(
        f"{code:25s}  nc={spec.nc_variable:6s}  units={spec.units:25s}  is_flux={spec.is_flux}"
    )

from earthlens.ecmwf import Catalog cat = Catalog() for code in ("sea-surface-temperature", "sea-ice-cover"): spec = cat.get_variable("reanalysis-era5-single-levels", code) print( f"{code:25s} nc={spec.nc_variable:6s} units={spec.units:25s} is_flux={spec.is_flux}" )

sea-surface-temperature    nc=sst     units=K                          is_flux=False
sea-ice-cover              nc=siconc  units=(0 - 1)                    is_flux=False

Step 2 — retrieve a year of monthly means¶

Domain: Iceland Sea (60°–70°N, 30°W–10°W). Just on the border of the seasonal sea-ice edge so we get visible variability in both fields. Pulling 12 months keeps the retrieve small.

In [2]:

from pathlib import Path
from earthlens import EarthLens, AggregationConfig

OUT = Path("data/era5-iceland-sea")
OUT.mkdir(parents=True, exist_ok=True)

earthlens = EarthLens(
    data_source="ecmwf",
    temporal_resolution="monthly",
    start="2022-01-01",
    end="2022-12-01",
    variables={
        "reanalysis-era5-single-levels-monthly-means": [
            "sea-surface-temperature",
            "sea-ice-cover",
        ],
    },
    lat_lim=[60.0, 70.0],
    lon_lim=[-30.0, -10.0],
    path=str(OUT),
)
earthlens.download(aggregate=AggregationConfig(freq="1MS", op="auto"))

from pathlib import Path from earthlens import EarthLens, AggregationConfig OUT = Path("data/era5-iceland-sea") OUT.mkdir(parents=True, exist_ok=True) earthlens = EarthLens( data_source="ecmwf", temporal_resolution="monthly", start="2022-01-01", end="2022-12-01", variables={ "reanalysis-era5-single-levels-monthly-means": [ "sea-surface-temperature", "sea-ice-cover", ], }, lat_lim=[60.0, 70.0], lon_lim=[-30.0, -10.0], path=str(OUT), ) earthlens.download(aggregate=AggregationConfig(freq="1MS", op="auto"))

2026-05-10 01:41:55.633 | INFO     | earthlens.ecmwf.backend:download:536 - Download ECMWF reanalysis-era5-single-levels-monthly-means/sea-surface-temperature data for period 2022-01-01 00:00:00 till 2022-12-01 00:00:00

2026-05-10 01:41:56.418 | INFO     | earthlens.ecmwf.backend:_api:724 - Requesting reanalysis-era5-single-levels-monthly-means from CDS; this may take several minutes

2026-05-10 01:41:56,759 INFO Request ID is 9a6c26d6-b484-4025-964d-2dd1d75accec

2026-05-10 01:41:56,923 INFO status has been updated to accepted

2026-05-10 01:42:29,652 INFO status has been updated to successful

2026-05-10 01:42:32 | INFO | pyramids.base.config | Logging is configured.

C:\gdrive\algorithms\remote-sensing\earthlens\src\earthlens\aggregate.py:385: RuntimeWarning: Mean of empty slice
  result = reducer(arr, axis=0)
2026-05-10 01:42:32.756 | INFO     | earthlens.ecmwf.backend:download:536 - Download ECMWF reanalysis-era5-single-levels-monthly-means/sea-ice-cover data for period 2022-01-01 00:00:00 till 2022-12-01 00:00:00

2026-05-10 01:42:32.757 | INFO     | earthlens.ecmwf.backend:_api:724 - Requesting reanalysis-era5-single-levels-monthly-means from CDS; this may take several minutes

2026-05-10 01:42:33,144 INFO Request ID is 770b3cfc-3eea-4f6a-8bce-8e5672d093d3

2026-05-10 01:42:33 | INFO | ecmwf.datastores.legacy_client | Request ID is 770b3cfc-3eea-4f6a-8bce-8e5672d093d3

2026-05-10 01:42:33,227 INFO status has been updated to accepted

2026-05-10 01:42:33 | INFO | ecmwf.datastores.legacy_client | status has been updated to accepted

2026-05-10 01:42:54,662 INFO status has been updated to successful

2026-05-10 01:42:54 | INFO | ecmwf.datastores.legacy_client | status has been updated to successful

2026-05-10 01:42:54 | INFO | multiurl.base | Downloading https://object-store.os-api.cci2.ecmwf.int:443/cci2-prod-cache-2/2026-05-09/5c9afca3b315767458544ea452a42ee2.nc

C:\gdrive\algorithms\remote-sensing\earthlens\src\earthlens\aggregate.py:385: RuntimeWarning: Mean of empty slice
  result = reducer(arr, axis=0)
2026-05-10 01:42:55.522 | INFO     | earthlens.ecmwf.backend:download:575 - ECMWF download summary: all 2 variables succeeded ([('reanalysis-era5-single-levels-monthly-means', 'sea-surface-temperature'), ('reanalysis-era5-single-levels-monthly-means', 'sea-ice-cover')])

Step 3 — extract domain-mean time series¶

Stack each variable's monthly GeoTIFFs and average over space (mask land NaNs out).

In [3]:

import numpy as np
import pandas as pd
from pyramids.dataset import Dataset

agg_dir = OUT / "aggregated"

def stack(cds_variable: str) -> np.ndarray:
    paths = sorted(agg_dir.glob(f"{cds_variable}_1MS_*.tif"))
    return np.stack([Dataset.read_file(str(p)).read_array() for p in paths])

sst = stack("sea_surface_temperature")  # K, NaN over land
ice = stack("sea_ice_cover")             # 0..1 fraction

months = pd.date_range("2022-01-01", periods=12, freq="MS")
sst_C = np.nanmean(sst, axis=(1, 2)) - 273.15
ice_pct = 100.0 * np.nanmean(ice, axis=(1, 2))

pd.DataFrame({"SST [°C]": sst_C.round(2), "Ice cover [%]": ice_pct.round(1)}, index=months)

import numpy as np import pandas as pd from pyramids.dataset import Dataset agg_dir = OUT / "aggregated" def stack(cds_variable: str) -> np.ndarray: paths = sorted(agg_dir.glob(f"{cds_variable}_1MS_*.tif")) return np.stack([Dataset.read_file(str(p)).read_array() for p in paths]) sst = stack("sea_surface_temperature") # K, NaN over land ice = stack("sea_ice_cover") # 0..1 fraction months = pd.date_range("2022-01-01", periods=12, freq="MS") sst_C = np.nanmean(sst, axis=(1, 2)) - 273.15 ice_pct = 100.0 * np.nanmean(ice, axis=(1, 2)) pd.DataFrame({"SST [°C]": sst_C.round(2), "Ice cover [%]": ice_pct.round(1)}, index=months)

Out[3]:

	SST [°C]	Ice cover [%]
2022-01-01	4.52	8.0
2022-02-01	4.22	7.0
2022-03-01	4.29	4.7
2022-04-01	4.55	6.2
2022-05-01	5.38	6.6
2022-06-01	6.94	4.0
2022-07-01	8.16	0.5
2022-08-01	8.71	0.0
2022-09-01	8.42	0.0
2022-10-01	6.71	0.1
2022-11-01	6.12	1.9
2022-12-01	5.52	6.4

Step 4 — plot the seasonal cycle on twin axes¶

SST and sea-ice cover trade off seasonally — winter ice maximum coincides with the SST minimum, summer melt with the SST peak.

In [4]:

import matplotlib.pyplot as plt

fig, ax1 = plt.subplots(figsize=(9, 5))
color1, color2 = "tab:red", "tab:blue"

ax1.plot(months, sst_C, marker="o", color=color1, label="SST")
ax1.set_ylabel("SST [°C]", color=color1)
ax1.tick_params(axis="y", labelcolor=color1)

ax2 = ax1.twinx()
ax2.plot(months, ice_pct, marker="s", color=color2, label="Sea-ice cover")
ax2.set_ylabel("Sea-ice cover [%]", color=color2)
ax2.tick_params(axis="y", labelcolor=color2)

ax1.set_title("Iceland Sea — monthly SST and sea-ice cover, 2022")
ax1.grid(alpha=0.3)
plt.tight_layout()
plt.show()

import matplotlib.pyplot as plt fig, ax1 = plt.subplots(figsize=(9, 5)) color1, color2 = "tab:red", "tab:blue" ax1.plot(months, sst_C, marker="o", color=color1, label="SST") ax1.set_ylabel("SST [°C]", color=color1) ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() ax2.plot(months, ice_pct, marker="s", color=color2, label="Sea-ice cover") ax2.set_ylabel("Sea-ice cover [%]", color=color2) ax2.tick_params(axis="y", labelcolor=color2) ax1.set_title("Iceland Sea — monthly SST and sea-ice cover, 2022") ax1.grid(alpha=0.3) plt.tight_layout() plt.show()

Step 5 — winter vs summer SST maps¶

Compare the February (typical ice maximum) and August (ice minimum) spatial patterns side-by-side.

In [5]:

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
for ax, idx, label in zip(axes, (1, 7), ("February", "August")):
    img = sst[idx] - 273.15
    im = ax.imshow(img, cmap="RdBu_r", origin="upper")
    ax.set_title(f"{label} 2022 SST [°C]")
    ax.set_xlabel("lon (pixels)")
    ax.set_ylabel("lat (pixels)")
    fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 2, figsize=(11, 4)) for ax, idx, label in zip(axes, (1, 7), ("February", "August")): img = sst[idx] - 273.15 im = ax.imshow(img, cmap="RdBu_r", origin="upper") ax.set_title(f"{label} 2022 SST [°C]") ax.set_xlabel("lon (pixels)") ax.set_ylabel("lat (pixels)") fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04) plt.tight_layout() plt.show()

Notes¶

• NaN over land. SST is undefined over land; ERA5 fills with NaN. np.nanmean correctly excludes those pixels from the domain mean. Same for sea-ice cover.
• For a proper ocean reanalysis, look at ORAS5. ERA5's SST is the atmospheric model's surface boundary condition (interpolated from HadISST/OSTIA), not a full ocean state. ORAS5 is on the catalog as "reanalysis-oras5" and exposes proper 3-D fields like potential-temperature and salinity with vertical_resolution: all_levels. ORAS5 carries request_kind: oceanic_monthly so the request shape strips day/time/area automatically.
• Daily SST is also available. Pass temporal_resolution="daily" and the catalog's daily dataset name (reanalysis-era5-single-levels) for finer-resolution time series.