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 DatasetCollection
from pyramids.plot import Basemap

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

# One spelling of the axis label, shared by the table, the twin axis and the bar.
SST_LABEL = "SST [°C]"

from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from pyramids.dataset import DatasetCollection from pyramids.plot import Basemap from earthlens.core import AggregationConfig, EarthLens from earthlens.ecmwf import Catalog # One spelling of the axis label, shared by the table, the twin axis and the bar. SST_LABEL = "SST [°C]"

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

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

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=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=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¶

DatasetCollection.from_files opens each variable's monthly GeoTIFFs as one cube whose band order follows the sorted file names (the collection carries no time axis of its own), so the rasters keep their georeferencing instead of being flattened into a bare array.

In [ ]:

agg_dir = OUT / "aggregated"

sst = DatasetCollection.from_files(agg_dir, glob="sea_surface_temperature_*_1MS_*.tif")
ice = DatasetCollection.from_files(agg_dir, glob="sea_ice_cover_*_1MS_*.tif")
print("sst cube", sst.shape, "| ice cube", ice.shape)

agg_dir = OUT / "aggregated" sst = DatasetCollection.from_files(agg_dir, glob="sea_surface_temperature_*_1MS_*.tif") ice = DatasetCollection.from_files(agg_dir, glob="sea_ice_cover_*_1MS_*.tif") print("sst cube", sst.shape, "| ice cube", ice.shape)

Convert SST to °C and ice cover to a percentage. The domain mean is a spatial reduction — one number per month — so it comes from each member raster's stats(), which already excludes the declared no-data. The collection's own mean() / min() / max() reduce across time instead, returning one map, so they answer a different question.

In [ ]:

months = pd.date_range("2022-01-01", periods=12, freq="MS")


def domain_mean(collection):
    """Spatial mean of each time step, straight from the rasters' own stats."""
    means = []
    for index in range(len(collection)):
        raster = collection.iloc(index)
        stats = raster.stats(approx_ok=False)
        means.append(float(stats["mean"].iloc[0]))
    return np.array(means)


sst_C = domain_mean(sst) - 273.15
ice_pct = 100.0 * domain_mean(ice)

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

months = pd.date_range("2022-01-01", periods=12, freq="MS") def domain_mean(collection): """Spatial mean of each time step, straight from the rasters' own stats.""" means = [] for index in range(len(collection)): raster = collection.iloc(index) stats = raster.stats(approx_ok=False) means.append(float(stats["mean"].iloc[0])) return np.array(means) sst_C = domain_mean(sst) - 273.15 ice_pct = 100.0 * domain_mean(ice) pd.DataFrame( {SST_LABEL: 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_LABEL, 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_LABEL, 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) patterns side by side. Both panels share one colour scale, so a single bar below serves them both, and basemap=Basemap() fills the land — which the SST rasters leave as no-data — with relief and a coastline, so Iceland reads as Iceland rather than a hole in the data.

In [ ]:

panels = {"February": sst.iloc(1), "August": sst.iloc(7)}
celsius = {
    label: raster.apply(lambda a: a - 273.15) for label, raster in panels.items()
}

bounds = []
for raster in celsius.values():
    stats = raster.stats(approx_ok=False)
    bounds.append((float(stats["min"].iloc[0]), float(stats["max"].iloc[0])))
vmin = min(low for low, _ in bounds)
vmax = max(high for _, high in bounds)
print(f"shared scale: {vmin:.2f} .. {vmax:.2f} C")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.4))
glyphs = []
for ax, (label, raster) in zip(axes, celsius.items()):
    glyphs.append(
        raster.plot(
            fig=fig,
            ax=ax,
            cmap="RdBu_r",
            vmin=vmin,
            vmax=vmax,
            colorbar=False,
            basemap=Basemap(),
            title=f"{label} 2022",
        )
    )

# Both panels share one scale, so one bar serves both -- drawn on the figure
# rather than on either axes, so neither map loses width to it. Either glyph
# would do; take the first rather than whichever the loop happened to leave
# bound.
bar = fig.colorbar(
    glyphs[0].im,
    ax=list(axes),
    orientation="horizontal",
    fraction=0.06,
    pad=0.10,
    shrink=0.7,
)
bar.set_label(SST_LABEL)

panels = {"February": sst.iloc(1), "August": sst.iloc(7)} celsius = { label: raster.apply(lambda a: a - 273.15) for label, raster in panels.items() } bounds = [] for raster in celsius.values(): stats = raster.stats(approx_ok=False) bounds.append((float(stats["min"].iloc[0]), float(stats["max"].iloc[0]))) vmin = min(low for low, _ in bounds) vmax = max(high for _, high in bounds) print(f"shared scale: {vmin:.2f} .. {vmax:.2f} C") fig, axes = plt.subplots(1, 2, figsize=(11, 4.4)) glyphs = [] for ax, (label, raster) in zip(axes, celsius.items()): glyphs.append( raster.plot( fig=fig, ax=ax, cmap="RdBu_r", vmin=vmin, vmax=vmax, colorbar=False, basemap=Basemap(), title=f"{label} 2022", ) ) # Both panels share one scale, so one bar serves both -- drawn on the figure # rather than on either axes, so neither map loses width to it. Either glyph # would do; take the first rather than whichever the loop happened to leave # bound. bar = fig.colorbar( glyphs[0].im, ax=list(axes), orientation="horizontal", fraction=0.06, pad=0.10, shrink=0.7, ) bar.set_label(SST_LABEL)

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.