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.

Setup¶

Consolidate the imports up front. earthlens provides the unified EarthLens entry point plus AggregationConfig and the ECMWF Catalog; pyramids provides Dataset for reading the downloaded GeoTIFFs.

In [ ]:

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pyramids.dataset import Dataset

from earthlens import AggregationConfig, EarthLens
from earthlens.ecmwf import Catalog

from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from pyramids.dataset import Dataset from earthlens import AggregationConfig, EarthLens from earthlens.ecmwf import Catalog

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 [ ]:

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}"
    )

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}" )

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.

Output directory¶

Pick a local folder for the GeoTIFFs and create it if needed.

In [ ]:

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

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

Build the request¶

Configure the ECMWF backend with the monthly-means dataset, both variables, the Iceland Sea bounding box, and the output path.

In [ ]:

earthlens = EarthLens(
    data_source="ecmwf",
    cadence="monthly",
    start="2022-01-01",
    end="2022-12-01",
    dataset="reanalysis-era5-single-levels-monthly-means",
    variables=[
        "sea-surface-temperature",
        "sea-ice-cover",
    ],
    aoi=[-30.0, 60.0, -10.0, 70.0],
    path=str(OUT),
)

earthlens = EarthLens( data_source="ecmwf", cadence="monthly", start="2022-01-01", end="2022-12-01", dataset="reanalysis-era5-single-levels-monthly-means", variables=[ "sea-surface-temperature", "sea-ice-cover", ], aoi=[-30.0, 60.0, -10.0, 70.0], path=str(OUT), )

Download and aggregate¶

download() retrieves each variable from the CDS and aggregates it to a monthly mean (op="auto" resolves to mean for these state variables). This may take a few minutes while the CDS processes the request.

In [ ]:

earthlens.download(aggregate=AggregationConfig(freq="1MS", op="auto"))

earthlens.download(aggregate=AggregationConfig(freq="1MS", op="auto"))

Step 3 — extract domain-mean time series¶

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

In [ ]:

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])

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])

Stack both variables, convert SST to °C and ice cover to a percentage, then collapse each cube to a domain-mean monthly series with np.nanmean (which skips the land NaNs).

In [ ]:

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
)

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 )

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 [ ]:

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()

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 [ ]:

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.