Snow & cryosphere — Alpine winter snowpack¶

Track the build-up and ablation of a winter snowpack at a mountain site using ERA5-Land's snow variables. Two complementary fields:

• snow-depth — geometric depth of the snowpack in metres. Convenient for visual / engineering work.
• snow-depth-water-equivalent — water content of the snowpack in metres of water equivalent. The hydrologically meaningful quantity (snowmelt feeds runoff during spring).

Both are state variables (instantaneous), so monthly aggregation uses op="auto" → mean.

Domain context. Snow water equivalent (SWE) is the operational metric for spring melt forecasting, reservoir management, and avalanche risk. Mountain SWE typically peaks in March–April depending on elevation, then drops to zero by mid-summer.

Setup¶

All the imports up front. earthlens provides the unified EarthLens entry point plus AggregationConfig for the monthly aggregation; pyramids provides Dataset for reading the written GeoTIFFs; numpy / pandas / matplotlib handle the time-series maths and plotting.

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

Step 1 — pull a winter at monthly resolution¶

Box: 1° around the central Alps (46°–47°N, 10°–11°E) — Ortler / Engadine area. Range Oct 2021–Jun 2022 covers a full snow season.

Choose an output directory¶

Each ERA5-Land variable is written as a GeoTIFF under OUT, with the aggregated monthly fields landing in an aggregated/ subdirectory.

In [ ]:

OUT = Path("data/era5-land-alps-snow")
OUT.mkdir(parents=True, exist_ok=True)

OUT = Path("data/era5-land-alps-snow") OUT.mkdir(parents=True, exist_ok=True)

Build the ECMWF request¶

Configure the backend with the dataset, the three snow / temperature variables, the bounding box, and the date window. Nothing is fetched until download() runs.

In [ ]:

earthlens = EarthLens(
    data_source="ecmwf",
    cadence="monthly",
    start="2021-10-01",
    end="2022-06-01",
    dataset="reanalysis-era5-land-monthly-means",
    variables=[
        "snow-depth",
        "snow-depth-water-equivalent",
        "2m-temperature",
    ],
    aoi=[10.0, 46.0, 11.0, 47.0],
    path=str(OUT),
)

earthlens = EarthLens( data_source="ecmwf", cadence="monthly", start="2021-10-01", end="2022-06-01", dataset="reanalysis-era5-land-monthly-means", variables=[ "snow-depth", "snow-depth-water-equivalent", "2m-temperature", ], aoi=[10.0, 46.0, 11.0, 47.0], path=str(OUT), )

Download and aggregate to monthly means¶

download() fetches each variable from the CDS and aggregates it to a 1MS monthly mean on a 0.1° grid (op="auto" resolves to mean for these state variables).

In [ ]:

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

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

Step 2 — assemble the seasonal time series¶

Extract the spatial mean over the box for each month.

Stack the per-month rasters¶

stack() reads every aggregated GeoTIFF for one variable (sorted by month) and stacks them into a single array, ready for a spatial average.

In [ ]:

agg = OUT / "aggregated"


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

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

Reduce each field to a box mean¶

Average over the spatial dimensions to get one value per month for snow depth, SWE, and 2 m temperature (converted from kelvin to °C).

In [ ]:

sde = np.nanmean(stack("snow_depth"), axis=(1, 2))  # m (geometric)
swe = np.nanmean(stack("snow_depth_water_equivalent"), axis=(1, 2))  # m of water eq
T2m = np.nanmean(stack("2m_temperature"), axis=(1, 2)) - 273.15  # °C

sde = np.nanmean(stack("snow_depth"), axis=(1, 2)) # m (geometric) swe = np.nanmean(stack("snow_depth_water_equivalent"), axis=(1, 2)) # m of water eq T2m = np.nanmean(stack("2m_temperature"), axis=(1, 2)) - 273.15 # °C

Build the monthly table¶

Assemble the box means into a DataFrame indexed by month, adding the bulk snow density (swe / sde) as a derived column.

In [ ]:

months = pd.date_range("2021-10-01", periods=len(sde), freq="MS")
df = pd.DataFrame(
    {
        "Snow depth [m]": sde.round(3),
        "SWE [m]": swe.round(3),
        "Snow density": (swe / np.where(sde > 0, sde, np.nan)).round(2),
        "T_2m [°C]": T2m.round(1),
    },
    index=months,
)
df

months = pd.date_range("2021-10-01", periods=len(sde), freq="MS") df = pd.DataFrame( { "Snow depth [m]": sde.round(3), "SWE [m]": swe.round(3), "Snow density": (swe / np.where(sde > 0, sde, np.nan)).round(2), "T_2m [°C]": T2m.round(1), }, index=months, ) df

Step 3 — plot snowpack evolution¶

SWE rises through fall as cold-season precipitation accumulates, peaks in late winter / early spring, and drops to zero through the melt season. The 0 °C 2 m temperature line is a useful reference for when the snowpack is energetically melting vs accumulating.

In [ ]:

fig, ax1 = plt.subplots(figsize=(9, 5))
ax1.fill_between(months, 0, sde, alpha=0.4, color="tab:blue", label="Snow depth [m]")
ax1.plot(months, swe, marker="o", color="tab:cyan", label="SWE [m]")
ax1.set_ylabel("Snow [m]", color="tab:blue")

ax2 = ax1.twinx()
ax2.plot(months, T2m, marker="^", color="tab:red", label="T_2m [°C]")
ax2.axhline(0, color="gray", lw=0.5, ls="--")
ax2.set_ylabel("T_2m [°C]", color="tab:red")

ax1.set_title("Central Alps snowpack evolution — Oct 2021 to Jun 2022")
fig.legend(loc="upper right", bbox_to_anchor=(0.98, 0.95))
ax1.grid(alpha=0.3)
plt.tight_layout()
plt.show()

fig, ax1 = plt.subplots(figsize=(9, 5)) ax1.fill_between(months, 0, sde, alpha=0.4, color="tab:blue", label="Snow depth [m]") ax1.plot(months, swe, marker="o", color="tab:cyan", label="SWE [m]") ax1.set_ylabel("Snow [m]", color="tab:blue") ax2 = ax1.twinx() ax2.plot(months, T2m, marker="^", color="tab:red", label="T_2m [°C]") ax2.axhline(0, color="gray", lw=0.5, ls="--") ax2.set_ylabel("T_2m [°C]", color="tab:red") ax1.set_title("Central Alps snowpack evolution — Oct 2021 to Jun 2022") fig.legend(loc="upper right", bbox_to_anchor=(0.98, 0.95)) ax1.grid(alpha=0.3) plt.tight_layout() plt.show()

Notes¶

• Snow depth vs SWE. Geometric depth (sde) and water equivalent (sd) differ by snow density (~0.1–0.5). New snow is fluffy (~0.1); compacted late-winter snowpack is denser (~0.3–0.5). The ratio swe / sde recovers the bulk density.
• ERA5-Land vs single-levels. ERA5-Land at 0.1° captures Alpine topography much better than the 0.25° single-levels product. The bbox-mean still smooths over high peaks vs valleys; for site-level work, slice to a specific pixel rather than averaging.
• Snowfall vs snowpack. ERA5 also reports snowfall (a flux — monthly accumulation). For the seasonal accumulation story, sum monthly snowfall; for the snowpack at a moment, use SWE.
• Glacier-relevant variables. Combine SWE with ERA5 2m-temperature and surface-net-solar-radiation for a temperature-index melt model (degree-day approach).