FABDEM quickstart — bare-earth DEM, and how it differs from the surface DEM¶

This notebook downloads a small FABDEM V1-2 subset, visualises the bare-earth terrain, then downloads the Copernicus GLO-30 surface DEM over the same area and shows the difference: FABDEM has forest canopy and building heights removed, which is exactly what you want for flood routing.

Licence. FABDEM is CC-BY-NC-SA 4.0 (non-commercial) — download() emits a LicenseWarning. For commercial use, obtain a licence from Fathom. FABDEM ships as 0.8–2.4 GB 10° bundles, so this live download takes a few minutes; keep the bounding box small.

In [ ]:

import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from pyramids.dataset import Dataset, GeoReference
from pyramids.plot import ColorBar, ColorScaling

from earthlens.core import EarthLens

out = Path(tempfile.mkdtemp(prefix="fabdem-"))
# A small coastal AOI (SE England) with a mix of open ground and settlement.
lat_lim = [50.85, 50.95]
lon_lim = [0.05, 0.15]
print("downloads land under:", out)

# Most of this AOI lies below 30 m with a thin tail to ~150 m, so the class
# breaks crowd the low ground, where the terrain actually varies.
ELEVATION_BREAKS = [0, 5, 10, 20, 30, 50, 75, 100, 150]

import tempfile from pathlib import Path import matplotlib.pyplot as plt import numpy as np from pyramids.dataset import Dataset, GeoReference from pyramids.plot import ColorBar, ColorScaling from earthlens.core import EarthLens out = Path(tempfile.mkdtemp(prefix="fabdem-")) # A small coastal AOI (SE England) with a mix of open ground and settlement. lat_lim = [50.85, 50.95] lon_lim = [0.05, 0.15] print("downloads land under:", out) # Most of this AOI lies below 30 m with a thin tail to ~150 m, so the class # breaks crowd the low ground, where the terrain actually varies. ELEVATION_BREAKS = [0, 5, 10, 20, 30, 50, 75, 100, 150]

Download the bare-earth DEM¶

download() emits a LicenseWarning — FABDEM is CC-BY-NC-SA 4.0, so it may not be used commercially without a licence from Fathom. It is left to surface as a real warning rather than being captured and reprinted.

In [ ]:

fabdem_paths = EarthLens(
    data_source="fabdem",
    lat_lim=lat_lim,
    lon_lim=lon_lim,
    path=out / "fabdem",
).download()
fabdem_paths

fabdem_paths = EarthLens( data_source="fabdem", lat_lim=lat_lim, lon_lim=lon_lim, path=out / "fabdem", ).download() fabdem_paths

Visualise the terrain¶

In [ ]:

fabdem = Dataset.read_file(fabdem_paths[0])
bare_masked = fabdem.read_array(masked=True)
bare = bare_masked.filled(np.nan)

stats = fabdem.stats(approx_ok=False)
print(
    f"elevation range: {float(stats['min'].iloc[0]):.1f} .. "
    f"{float(stats['max'].iloc[0]):.1f} m"
)

fabdem.plot(
    cmap="terrain",
    color=ColorScaling.boundary(bounds=ELEVATION_BREAKS),
    colorbar=ColorBar(label="elevation (m)"),
    title="FABDEM V1-2 bare-earth elevation",
)

fabdem = Dataset.read_file(fabdem_paths[0]) bare_masked = fabdem.read_array(masked=True) bare = bare_masked.filled(np.nan) stats = fabdem.stats(approx_ok=False) print( f"elevation range: {float(stats['min'].iloc[0]):.1f} .. " f"{float(stats['max'].iloc[0]):.1f} m" ) fabdem.plot( cmap="terrain", color=ColorScaling.boundary(bounds=ELEVATION_BREAKS), colorbar=ColorBar(label="elevation (m)"), title="FABDEM V1-2 bare-earth elevation", )

Compare against the Copernicus surface DEM¶

The Copernicus GLO-30 DEM is a surface model — it includes tree canopy and building tops. Subtracting FABDEM (bare-earth) from it leaves the removed above-ground height (canopy + buildings), which should be roughly zero over open ground and positive over woodland and built-up areas.

In [ ]:

copdem_paths = EarthLens(
    data_source="dem",
    dataset="cop-dem-glo-30",
    lat_lim=lat_lim,
    lon_lim=lon_lim,
    path=out / "copdem",
).download()
copdem_paths

copdem_paths = EarthLens( data_source="dem", dataset="cop-dem-glo-30", lat_lim=lat_lim, lon_lim=lon_lim, path=out / "copdem", ).download() copdem_paths

In [ ]:

from pyramids.dataset.merge import merge_rasters

# Align the Copernicus tile(s) onto the FABDEM grid so the arrays subtract cleanly.
aligned = out / "copdem_aligned.tif"
# merge_rasters defaults no_data_value to 0. This AOI is coastal, so sea level
# sits at exactly 0 m over most of the tile; the default would mark every one of
# those cells no-data. The Copernicus tiles declare none, so none is inherited.
merge_rasters(
    src=list(copdem_paths),
    dst=aligned,
    dst_crs=None,
    resampling="bilinear",
    no_data_value="none",
)
cop_raw = Dataset.read_file(aligned)
cop = cop_raw.align(fabdem)
surface_masked = cop.read_array(masked=True)
surface = surface_masked.filled(np.nan)
# `no_data_value="none"` above means the aligned raster declares nothing, so
# masked=True masks nothing. Copernicus carries its voids as large negatives,
# and an unmasked void would subtract into the difference map as hundreds of
# metres of "canopy" — drop anything below a physical floor for this coast.
surface[surface < -1000] = np.nan

diff = surface - bare  # removed above-ground height (canopy + buildings)

geo_ref = GeoReference(geo=fabdem.geotransform, epsg=fabdem.epsg)
bare_grid = Dataset.from_array(bare, no_data_value=np.nan, geo_ref=geo_ref)
surface_grid = Dataset.from_array(surface, no_data_value=np.nan, geo_ref=geo_ref)
diff_grid = Dataset.from_array(diff, no_data_value=np.nan, geo_ref=geo_ref)

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
elevation_scale = ColorScaling.boundary(bounds=ELEVATION_BREAKS)
for ax, grid, title, cmap, scale in [
    (axes[0], bare_grid, "FABDEM (bare-earth)", "terrain", elevation_scale),
    (axes[1], surface_grid, "Copernicus GLO-30 (surface)", "terrain", elevation_scale),
    (axes[2], diff_grid, "surface - bare  (canopy + buildings)", "viridis", None),
]:
    grid.plot(fig=fig, ax=ax, cmap=cmap, color=scale, title=title)

from pyramids.dataset.merge import merge_rasters # Align the Copernicus tile(s) onto the FABDEM grid so the arrays subtract cleanly. aligned = out / "copdem_aligned.tif" # merge_rasters defaults no_data_value to 0. This AOI is coastal, so sea level # sits at exactly 0 m over most of the tile; the default would mark every one of # those cells no-data. The Copernicus tiles declare none, so none is inherited. merge_rasters( src=list(copdem_paths), dst=aligned, dst_crs=None, resampling="bilinear", no_data_value="none", ) cop_raw = Dataset.read_file(aligned) cop = cop_raw.align(fabdem) surface_masked = cop.read_array(masked=True) surface = surface_masked.filled(np.nan) # `no_data_value="none"` above means the aligned raster declares nothing, so # masked=True masks nothing. Copernicus carries its voids as large negatives, # and an unmasked void would subtract into the difference map as hundreds of # metres of "canopy" — drop anything below a physical floor for this coast. surface[surface < -1000] = np.nan diff = surface - bare # removed above-ground height (canopy + buildings) geo_ref = GeoReference(geo=fabdem.geotransform, epsg=fabdem.epsg) bare_grid = Dataset.from_array(bare, no_data_value=np.nan, geo_ref=geo_ref) surface_grid = Dataset.from_array(surface, no_data_value=np.nan, geo_ref=geo_ref) diff_grid = Dataset.from_array(diff, no_data_value=np.nan, geo_ref=geo_ref) fig, axes = plt.subplots(1, 3, figsize=(15, 4.5)) elevation_scale = ColorScaling.boundary(bounds=ELEVATION_BREAKS) for ax, grid, title, cmap, scale in [ (axes[0], bare_grid, "FABDEM (bare-earth)", "terrain", elevation_scale), (axes[1], surface_grid, "Copernicus GLO-30 (surface)", "terrain", elevation_scale), (axes[2], diff_grid, "surface - bare (canopy + buildings)", "viridis", None), ]: grid.plot(fig=fig, ax=ax, cmap=cmap, color=scale, title=title)

Takeaway¶

FABDEM gives you the terrain surface with vegetation and buildings stripped out — the bare earth that water actually flows over. The difference map highlights woodland and settlement, exactly the features that a surface DEM would wrongly treat as ground in a flood model.