EFHM quickstart — a 1-in-100-year flood-depth map¶

This notebook downloads the JRC European Flood Hazard Map for a small area of the Rhine delta (Netherlands) at the 100-year return period, reading only the AOI's pixel window over /vsicurl — so it transfers kilobytes, not gigabytes. The EFHM is CC-BY-4.0 (no licence warning).

In [ ]:

import math
import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
from cleopatra.styling.colors import DATA_STYLES
from pyramids.dataset import Dataset

from earthlens.core import EarthLens

# `oslo` from cleopatra's Crameri palettes, reversed so the ramp runs light ->
# dark and deeper water reads as more ink. Picked by measurement, not taste:
# `cleopatra.styling.perceptual.perceptual_uniformity` scores it 0.09 against
# matplotlib's `Blues` at 0.27 (0 == perfectly even steps), and it is the only
# candidate that stays a single blue hue end to end.
DEPTH_CMAP = DATA_STYLES["oslo"]["oslo"]["cmap"].reversed()

out = Path(tempfile.mkdtemp(prefix="efhm-"))
paths = EarthLens(
    data_source="jrc",
    lat_lim=[51.7, 52.0],
    lon_lim=[4.6, 5.1],
    return_periods=[100],
    path=out,
).download()
paths

import math import tempfile from pathlib import Path import matplotlib.pyplot as plt from cleopatra.styling.colors import DATA_STYLES from pyramids.dataset import Dataset from earthlens.core import EarthLens # `oslo` from cleopatra's Crameri palettes, reversed so the ramp runs light -> # dark and deeper water reads as more ink. Picked by measurement, not taste: # `cleopatra.styling.perceptual.perceptual_uniformity` scores it 0.09 against # matplotlib's `Blues` at 0.27 (0 == perfectly even steps), and it is the only # candidate that stays a single blue hue end to end. DEPTH_CMAP = DATA_STYLES["oslo"]["oslo"]["cmap"].reversed() out = Path(tempfile.mkdtemp(prefix="efhm-")) paths = EarthLens( data_source="jrc", lat_lim=[51.7, 52.0], lon_lim=[4.6, 5.1], return_periods=[100], path=out, ).download() paths

Map the water depth¶

Cells carry river-flood water depth in metres; -9999 marks no data (dry / outside the modelled river network).

In [ ]:

depth = Dataset.read_file(paths[0])

# `stats()` excludes the nodata; `approx_ok=False` forces a full read rather
# than an overview estimate, which matters when the value sets a colour range.
vmax = float(depth.stats(approx_ok=False)["max"].iloc[0])

# The figure is built here so it can carry `constrained_layout`, which sizes
# margins from the real text extents rather than matplotlib's fixed fractions
# and keeps the colour bar clear of the longitude labels. matplotlib owns only
# the layout; pyramids still renders the map, on its own lon/lat axes.
fig, ax = plt.subplots(figsize=(12, 7.2), constrained_layout=True)
glyph = depth.plot(
    fig=fig,
    ax=ax,
    cmap=DEPTH_CMAP,
    title="EFHM RP100 river-flood water depth",
    title_size=15,
    colorbar=False,
)

# Tidy the axes through the returned cleopatra glyph. Anchoring the range at 0
# and asking for whole-metre ticks replaces the default labels, which land on
# the data's own min/max and read as 0.100 / 1.601 / 5.102.
glyph.im.set_clim(0, vmax)
glyph.ax.xaxis.set_ticks_position("bottom")
glyph.ax.set_xlabel("longitude", fontsize=12)
glyph.ax.set_ylabel("latitude", fontsize=12)
glyph.ax.tick_params(labelsize=11)

# The bar runs horizontally under the map, so it does not compete with it for
# width -- the dimension that sizes an aspect-locked panel.
bar = fig.colorbar(
    glyph.im,
    ax=ax,
    orientation="horizontal",
    location="bottom",
    shrink=0.5,
    pad=0.02,
    aspect=40,
)
bar.set_label("depth (m)", size=12)
bar.set_ticks(range(0, math.floor(vmax) + 1))
bar.ax.tick_params(labelsize=11)

depth.stats(approx_ok=False)

depth = Dataset.read_file(paths[0]) # `stats()` excludes the nodata; `approx_ok=False` forces a full read rather # than an overview estimate, which matters when the value sets a colour range. vmax = float(depth.stats(approx_ok=False)["max"].iloc[0]) # The figure is built here so it can carry `constrained_layout`, which sizes # margins from the real text extents rather than matplotlib's fixed fractions # and keeps the colour bar clear of the longitude labels. matplotlib owns only # the layout; pyramids still renders the map, on its own lon/lat axes. fig, ax = plt.subplots(figsize=(12, 7.2), constrained_layout=True) glyph = depth.plot( fig=fig, ax=ax, cmap=DEPTH_CMAP, title="EFHM RP100 river-flood water depth", title_size=15, colorbar=False, ) # Tidy the axes through the returned cleopatra glyph. Anchoring the range at 0 # and asking for whole-metre ticks replaces the default labels, which land on # the data's own min/max and read as 0.100 / 1.601 / 5.102. glyph.im.set_clim(0, vmax) glyph.ax.xaxis.set_ticks_position("bottom") glyph.ax.set_xlabel("longitude", fontsize=12) glyph.ax.set_ylabel("latitude", fontsize=12) glyph.ax.tick_params(labelsize=11) # The bar runs horizontally under the map, so it does not compete with it for # width -- the dimension that sizes an aspect-locked panel. bar = fig.colorbar( glyph.im, ax=ax, orientation="horizontal", location="bottom", shrink=0.5, pad=0.02, aspect=40, ) bar.set_label("depth (m)", size=12) bar.set_ticks(range(0, math.floor(vmax) + 1)) bar.ax.tick_params(labelsize=11) depth.stats(approx_ok=False)

Global vs European coverage¶

This is the higher-fidelity, Europe-focused EFHM. For the global JRC flood hazard (which also covers Europe, at ~90 m), use the gee backend with asset="JRC/CEMS_GLOFAS/FloodHazard/v2_1".

Check the plan before you fetch¶

Each return period is a whole-Europe GeoTIFF of roughly 23 GB, so it is worth knowing what a request will do before it does it. count() reports how many products match and preview() shows the first few as plain dicts — neither downloads anything.

In [ ]:

plan = EarthLens(
    data_source="jrc",
    lat_lim=[51.7, 52.0],
    lon_lim=[4.6, 5.1],
    return_periods=[10, 100, 500],
    path=out,
)
print("products that match:", plan.count())
for row in plan.preview(3):
    print("  ", {key: row[key] for key in ("id", "rp") if key in row})

plan = EarthLens( data_source="jrc", lat_lim=[51.7, 52.0], lon_lim=[4.6, 5.1], return_periods=[10, 100, 500], path=out, ) print("products that match:", plan.count()) for row in plan.preview(3): print(" ", {key: row[key] for key in ("id", "rp") if key in row})

download() writes files; load() hands you the object¶

download() returns the paths it wrote. When you only want the data in memory — a quick look, a notebook, a step in a pipeline — load() runs the same request into a temporary directory and returns the object instead: a list of pyramids Datasets here, and a pandas.DataFrame for the coastal summary.

In [ ]:

grids = EarthLens(
    data_source="jrc",
    lat_lim=[51.8, 51.9],
    lon_lim=[4.8, 4.9],
    return_periods=[100],
).load()

print(type(grids).__name__, "of", type(grids[0]).__name__)
print(f"{grids[0].rows} x {grids[0].columns} cells, EPSG:{grids[0].epsg}")

grids = EarthLens( data_source="jrc", lat_lim=[51.8, 51.9], lon_lim=[4.8, 4.9], return_periods=[100], ).load() print(type(grids).__name__, "of", type(grids[0]).__name__) print(f"{grids[0].rows} x {grids[0].columns} cells, EPSG:{grids[0].epsg}")

Other ways to say where¶

lat_lim / lon_lim is the explicit form. aoi= also accepts a (lon, lat) point, which needs a buffer in degrees to become an area.

In [ ]:

point = EarthLens(
    data_source="jrc",
    aoi=(4.9, 51.8),
    buffer=0.1,
    return_periods=[100],
    path=out,
)
print("a point plus a 0.1 deg buffer matches", point.count(), "product(s)")

point = EarthLens( data_source="jrc", aoi=(4.9, 51.8), buffer=0.1, return_periods=[100], path=out, ) print("a point plus a 0.1 deg buffer matches", point.count(), "product(s)")

Where files land when you do not say¶

Omitting path= does not scatter files into the working directory: earthlens writes to a per-source subdirectory of its configured output directory, and set_output_dir() moves that. The cache directory is separate and holds intermediate downloads.

In [ ]:

from earthlens.core import cache_dir, output_dir

print("output dir:", output_dir())
print("cache dir :", cache_dir())

from earthlens.core import cache_dir, output_dir print("output dir:", output_dir()) print("cache dir :", cache_dir())

A key you will probably try first¶

sea-level-forecast looks like the obvious name, and it is deliberately not a valid key. Generic subjects are reserved: several providers could serve "sea level", so a bare topic is ambiguous and earthlens refuses it rather than picking one for you. The error names the qualified keys that do work.

In [ ]:

raises-exception

# NBVAL_RAISES_EXCEPTION
# Deliberately the wrong form -- this raises, and the message is the useful part.
EarthLens(data_source="sea-level-forecast")

# NBVAL_RAISES_EXCEPTION # Deliberately the wrong form -- this raises, and the message is the useful part. EarthLens(data_source="sea-level-forecast")