JRC sea-level (Total Water Level) forecasts¶
The earthlens.jrc backend serves the JRC / Copernicus-EMS probabilistic,
data-driven sea-level forecasts — storm surge + tide + wave-derived coastal
Total Water Level (TWL) — alongside the European Flood Hazard Map (EFHM).
There are two products, both global 0.25° NetCDF-4 cubes:
- medium-term — issued twice daily, 15-day horizon;
- subseasonal — issued weekly, ~46-day horizon, with a global per-country coastal-summary CSV alongside the gridded cube.
This notebook fetches the latest medium-term forecast for the North Sea, maps
it, and reads the global coastal summary. It downloads live from the JRC
open-data server (no credentials — CC-BY-4.0); only the small area of interest
is read over /vsicurl, not the whole multi-gigabyte cube.
import tempfile
from pathlib import Path
from cleopatra.styling.colors import DATA_STYLES
from pyramids.dataset import Dataset
from pyramids.plot import ColorBar
from earthlens.core import EarthLens
# The facade key for the gridded TWL forecast. Named once here because every
# fetch below asks for the same product; `product=` is what varies.
SEA_LEVEL = "jrc:sea-level-forecast"
# The per-country coastal summary that rides with the subseasonal run.
COASTAL = "jrc:coastal-forecast"
# Same `oslo` (Crameri) ramp reversed as the flood notebooks: TWL is the same
# job -- a strictly positive water height in metres -- so it gets the same
# sequential single-hue treatment. It is deliberately NOT the `sea_level_anomaly`
# preset: that one diverges about zero, and TWL75 is an absolute level (this
# North Sea window runs 0.27 to 3.54 m), never negative, so a diverging map
# would invent a polarity the data does not have.
TWL_CMAP = DATA_STYLES["oslo"]["oslo"]["cmap"].reversed()
OUT = Path(tempfile.mkdtemp(prefix="jrc_twl_"))
OUT
1. Gridded TWL forecast for a coastal area¶
Select product="medium_term" and a bounding box. reference_time defaults to
"latest", which resolves the newest complete forecast cycle (the backend
walks the dated directory tree and honours the endFls sentinel). The default
field is TWL75 — the 75th-percentile total water level — and every forecast
time step becomes a band of the written GeoTIFF.
paths = EarthLens(
data_source=SEA_LEVEL,
product="medium_term",
lat_lim=[50.0, 58.0], # North Sea
lon_lim=[-2.0, 10.0],
path=OUT,
).download()
paths
2. Map the first forecast step¶
Read the cropped GeoTIFF back with pyramids. It is georeferenced (EPSG:4326) and
carries one band per forecast day; land cells are NaN. We plot the first
forecast step.
grid = Dataset.read_file(paths[0])
print("bands (forecast steps):", grid.band_count, "| epsg:", grid.epsg)
print("geotransform:", grid.geotransform)
vmax = float(grid.stats(band=0, approx_ok=False)["max"].iloc[0])
# The cube is georeferenced, so `plot()` places it on its own lon/lat axes --
# no extent arithmetic from the geotransform is needed.
glyph = grid.plot(
band=0,
cmap=TWL_CMAP,
title="JRC medium-term TWL75 forecast — first step",
title_size=12,
figsize=(8, 5),
colorbar=ColorBar(label="total water level (m)", label_size=10, length=0.85),
)
# Tidy the axes through the returned cleopatra glyph. Half-metre ticks replace
# the defaults, which land on the data's own min/max and read to 3 decimals
# (0.270 / 1.576 / 3.536).
glyph.im.set_clim(0, vmax)
glyph.cbar.set_ticks([t / 2 for t in range(0, int(vmax * 2) + 1)])
glyph.ax.xaxis.set_ticks_position("bottom")
glyph.ax.set_xlabel("longitude", fontsize=10)
glyph.ax.set_ylabel("latitude", fontsize=10)
glyph.ax.tick_params(labelsize=9)
3. Global coastal summary¶
The subseasonal product also publishes a global, per-country coastal summary as
a table. The jrc:coastal-forecast key returns it directly as a pandas.DataFrame
(exceedance probabilities against return-period thresholds, plus a 1–10 severity
summary per country).
summary = EarthLens(data_source=COASTAL).download()
print("countries:", len(summary))
summary[["GID_0", "NAME_0", "summary_TWL_1_10"]].head(10)
Notes¶
- Distinct from the EFHM. The
efhm/jrc-floodkeys serve the static European river-flood depth map (per return period); the sea-level keys here serve the coastal TWL forecast. Both live in the oneearthlens.jrcbackend. - Windowed reads. Only the AOI window is transferred over
/vsicurl, so a small area costs little of the 13–38 GB cube. A forecast cycle is chosen byreference_time, soaggregate=is rejected. - Licence. CC-BY-4.0 (Copernicus EMS / EC JRC).
Asking for a specific cycle¶
Every example so far used the default reference_time="latest", which resolves
the newest complete cycle. A forecast archive is more useful when you can ask for
a named one — to reproduce a figure, or to compare what was predicted on two
different days.
The format is "YYYY-MM-DDTHH", and the hour must be a real issue time (00 or 12
for the twice-daily medium-term product).
named = EarthLens(
data_source=SEA_LEVEL,
product="medium_term",
reference_time="2026-08-26T12",
lat_lim=[51.0, 53.0],
lon_lim=[3.0, 5.0],
path=OUT,
)
print("cycle 2026-08-26T12 matches", named.count(), "product(s)")
When a cycle has aged out¶
The server keeps a rolling window of recent cycles, not the whole archive, so a date inside the product's stated period is not automatically retrievable. Asking for one that is gone fails immediately with a message naming the cycle, rather than returning an empty or silently substituted result.
# NBVAL_RAISES_EXCEPTION
# Well inside the medium-term period (2022 -> present), but long since rotated out.
EarthLens(
data_source=SEA_LEVEL,
product="medium_term",
reference_time="2020-01-01T00",
lat_lim=[51.0, 53.0],
lon_lim=[3.0, 5.0],
path=OUT,
).count()
The coastal summary in memory¶
The coastal product is a table, not a raster, so load() returns a
pandas.DataFrame directly — no file is written and nothing needs reading back.
summary = EarthLens(data_source=COASTAL).load()
print(type(summary).__name__, summary.shape)
summary.head()
Choosing a field, and what it costs¶
A cycle file holds 68 variables, not one. field= picks which becomes the
written raster, and the choice has a large cost implication, because the fields
differ in shape, not just meaning:
| field | bands | what it is |
|---|---|---|
TWL75 (default) |
16 | 75th-percentile total water level, one band per forecast step |
episWL75 |
16 | the same percentile for the epistemic-uncertainty band |
TWL, episWL |
800 | the full ensemble — 16 steps x 50 members |
TWLcoast, episWLcoast |
17247 | per-station coastal series on a 16 x 50 index grid, not a map |
The percentile summaries are what you usually want. Asking for TWL reads fifty
times the data for the same window, and the *coast variables are not
georeferenced grids at all — the backend refuses them rather than writing a
raster with meaningless axes.
import time
started = time.time()
episodic = EarthLens(
data_source=SEA_LEVEL,
product="medium_term",
field="episWL75",
lat_lim=[51.0, 53.0],
lon_lim=[3.0, 5.0],
path=OUT,
).download()
print(f"episWL75 fetched in {time.time() - started:.0f}s -> {episodic[0].name}")
The bands carry their forecast valid times¶
A 16-band raster is unreadable if the bands are anonymous, so each one is named
for the time step it holds, read from the cube's CF time coordinate. That is
what makes the written GeoTIFF self-describing once it leaves this notebook.
band_grid = Dataset.read_file(episodic[0])
print("bands:", band_grid.band_count)
for name in band_grid.band_names[:5]:
print(" ", name)
print(" ...")
print(" ", band_grid.band_names[-1])