EMODnet Bathymetry over OGC WCS — a guided tour¶

The bathymetry backend serves two kinds of Digital Terrain Model (DTM): the global GEBCO / ETOPO grids over NOAA ERDDAP griddap, and the European high-resolution EMODnet DTM over OGC WCS. This notebook focuses on the WCS side — what dataset="emodnet" gives you, how the request is shaped, how it compares to the global grids, how the release history evolved, and how the European-domain guard protects you from silent out-of-coverage data.

By the end you will be able to pull an EMODnet subset, read its depths, reason about resolution and coverage, and pick between EMODnet and a global DEM for a given area. The read path is pyramids.Dataset.from_wcs — earthlens supplies only the coverage id, AOI, CRS, and protocol version.

API: earthlens.core.EarthLens · earthlens.bathymetry.Catalog.

Setup¶

EarthLens is the unified entry point; pyramids reads the written GeoTIFF; pandas renders the catalog tables and matplotlib the maps. Downloads go to a temporary directory — every dataset here is fetched live, no local fixtures.

In [ ]:

import tempfile

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pyramids.dataset import Dataset

from earthlens.bathymetry import Catalog
from earthlens.core import EarthLens

import tempfile import matplotlib.pyplot as plt import numpy as np import pandas as pd from pyramids.dataset import Dataset from earthlens.bathymetry import Catalog from earthlens.core import EarthLens

We will fetch several small subsets, so a tiny helper keeps the teaching cells about the feature rather than the plumbing. It downloads one DEM subset, reads the elevation band, and masks the GEBCO no-data fill (32767) so land/edges do not skew the colour scale. aoi is [min_lon, min_lat, max_lon, max_lat].

In [ ]:

def fetch_depth(dataset, aoi):
    out = tempfile.mkdtemp()
    paths = EarthLens(
        data_source='bathymetry', dataset=dataset, aoi=aoi, path=out
    ).download(progress_bar=False)
    grid = np.asarray(Dataset.read_file(str(paths[0])).read_array(), dtype='float32')
    grid = grid[0] if grid.ndim == 3 else grid
    return np.where(grid == 32767, np.nan, grid)

def fetch_depth(dataset, aoi): out = tempfile.mkdtemp() paths = EarthLens( data_source='bathymetry', dataset=dataset, aoi=aoi, path=out ).download(progress_bar=False) grid = np.asarray(Dataset.read_file(str(paths[0])).read_array(), dtype='float32') grid = grid[0] if grid.ndim == 3 else grid return np.where(grid == 32767, np.nan, grid)

Two transports, one backend¶

The bathymetry catalog mixes ERDDAP-griddap rows (GEBCO / ETOPO) with OGC-wcs rows (EMODnet). Reading the catalog is offline. The transport column is what the backend dispatches on; the WCS rows also carry a wcs_version and a native_bbox (the coverage's advertised extent, used by the domain guard).

In [ ]:

catalog = Catalog()
rows = [
    {
        'dataset': d,
        'transport': r.transport,
        'coverage/id': r.dataset_id,
        'native_resolution': r.native_resolution,
    }
    for d in sorted(catalog.datasets)
    for r in [catalog.get(d)]
]
pd.DataFrame(rows).set_index('dataset')

catalog = Catalog() rows = [ { 'dataset': d, 'transport': r.transport, 'coverage/id': r.dataset_id, 'native_resolution': r.native_resolution, } for d in sorted(catalog.datasets) for r in [catalog.get(d)] ] pd.DataFrame(rows).set_index('dataset')

The three global rows are erddap-griddap; the five emodnet* rows are wcs. emodnet is the latest release; the year-stamped ids pin older DTMs. Everything below uses the WCS rows.

Quickstart — an EMODnet subset over the North Sea¶

The shortest end-to-end example: pick dataset="emodnet" and a small North Sea AOI. download() issues the WCS GetCoverage, crops to the AOI, and writes a GeoTIFF; fetch_depth hands back the depth array (metres, negative = below sea level).

In [ ]:

north_sea = [2.0, 53.0, 4.0, 55.0]
emodnet_depth = fetch_depth('emodnet', north_sea)
emodnet_depth.shape

north_sea = [2.0, 53.0, 4.0, 55.0] emodnet_depth = fetch_depth('emodnet', north_sea) emodnet_depth.shape

A quick numeric sanity check before mapping — the range should be plausible shelf-sea depths, all below sea level.

In [ ]:

print(
    'depth range',
    round(float(np.nanmin(emodnet_depth))),
    '..',
    round(float(np.nanmax(emodnet_depth))),
    'm',
)

print( 'depth range', round(float(np.nanmin(emodnet_depth))), '..', round(float(np.nanmax(emodnet_depth))), 'm', )

In [ ]:

plt.figure(figsize=(6, 5))
plt.imshow(emodnet_depth, cmap='Blues_r')
plt.colorbar(label='mean depth (m, negative = below sea level)')
plt.title('EMODnet DTM (~3.75″) — southern North Sea')
plt.axis('off')
plt.show()

plt.figure(figsize=(6, 5)) plt.imshow(emodnet_depth, cmap='Blues_r') plt.colorbar(label='mean depth (m, negative = below sea level)') plt.title('EMODnet DTM (~3.75″) — southern North Sea') plt.axis('off') plt.show()

The southern North Sea is shallow (tens of metres) and the ~3.75-arc-second grid resolves fine seabed structure — sandbanks and channels — that a global DEM would blur.

How the WCS request is shaped¶

earthlens never talks WCS dialects itself — it reads four fields off the catalog row and hands them to pyramids.Dataset.from_wcs. Inspect them for emodnet:

In [ ]:

row = catalog.get('emodnet')
{
    'endpoint': row.endpoint,
    'coverage': row.dataset_id,
    'wcs_version': row.wcs_version,
    'crs': row.crs,
    'native_bbox': row.native_bbox,
}

row = catalog.get('emodnet') { 'endpoint': row.endpoint, 'coverage': row.dataset_id, 'wcs_version': row.wcs_version, 'crs': row.crs, 'native_bbox': row.native_bbox, }

Field	Meaning	Value for emodnet
endpoint	The OGC WCS service URL	ows.emodnet-bathymetry.eu/wcs
coverage	WCS coverage id sent as GetCoverage	emodnet:mean
wcs_version	Protocol version GDAL negotiates	1.0.0
crs	Request / native CRS	EPSG:4326
native_bbox	Advertised extent (W, S, E, N) for the domain guard	(-70.5, 11, 43, 90)

The 1.0.0 version is deliberate: this GeoServer only subsets cleanly at WCS 1.0.0 with the colon coverage id emodnet:mean (the 2.0.1 emodnet__mean id fails GDAL subsetting).

Depth distribution¶

A histogram reads the seabed better than a single min/max — it shows where most of the AOI sits depth-wise.

In [ ]:

plt.figure(figsize=(6, 4))
plt.hist(emodnet_depth[np.isfinite(emodnet_depth)].ravel(), bins=60, color='steelblue')
plt.xlabel('depth (m)')
plt.ylabel('pixel count')
plt.title('EMODnet depth distribution — southern North Sea')
plt.show()

plt.figure(figsize=(6, 4)) plt.hist(emodnet_depth[np.isfinite(emodnet_depth)].ravel(), bins=60, color='steelblue') plt.xlabel('depth (m)') plt.ylabel('pixel count') plt.title('EMODnet depth distribution — southern North Sea') plt.show()

The mass sits in the shallow shelf range with a tail toward deeper channels — the classic North Sea profile.

Regional detail vs the global GEBCO grid¶

EMODnet is the coastal / shelf complement to the global DEMs. Fetch GEBCO over the same AOI and compare the grid shapes: same footprint, very different pixel counts.

In [ ]:

gebco_depth = fetch_depth('gebco_2020', north_sea)
pd.DataFrame(
    {
        'rows': [emodnet_depth.shape[0], gebco_depth.shape[0]],
        'cols': [emodnet_depth.shape[1], gebco_depth.shape[1]],
    },
    index=['emodnet (~3.75″)', 'gebco_2020 (15″)'],
)

gebco_depth = fetch_depth('gebco_2020', north_sea) pd.DataFrame( { 'rows': [emodnet_depth.shape[0], gebco_depth.shape[0]], 'cols': [emodnet_depth.shape[1], gebco_depth.shape[1]], }, index=['emodnet (~3.75″)', 'gebco_2020 (15″)'], )

In [ ]:

fig, axes = plt.subplots(1, 2, figsize=(11, 5))
axes[0].imshow(emodnet_depth, cmap='Blues_r')
axes[0].set_title('EMODnet (~3.75″)')
axes[0].axis('off')
axes[1].imshow(gebco_depth, cmap='Blues_r')
axes[1].set_title('GEBCO 2020 (15″)')
axes[1].axis('off')
plt.show()

fig, axes = plt.subplots(1, 2, figsize=(11, 5)) axes[0].imshow(emodnet_depth, cmap='Blues_r') axes[0].set_title('EMODnet (~3.75″)') axes[0].axis('off') axes[1].imshow(gebco_depth, cmap='Blues_r') axes[1].set_title('GEBCO 2020 (15″)') axes[1].axis('off') plt.show()

Same 2°×2° box, but EMODnet packs ~4× the pixels per axis, so seabed features sharpen. For an open-ocean or global request you would still reach for GEBCO / ETOPO; for European coastal detail, EMODnet.

Release history — the DTM improves over time¶

emodnet tracks the latest published DTM; the year-stamped rows pin older releases. Successive releases refine both resolution and coverage. Compare the 2016 release with the latest over the same North Sea box:

In [ ]:

emodnet2016_depth = fetch_depth('emodnet_2016', north_sea)
pd.DataFrame(
    {
        'grid_rows': [emodnet_depth.shape[0], emodnet2016_depth.shape[0]],
        'min_depth': [
            float(np.nanmin(emodnet_depth)),
            float(np.nanmin(emodnet2016_depth)),
        ],
        'mean_depth': [
            float(np.nanmean(emodnet_depth)),
            float(np.nanmean(emodnet2016_depth)),
        ],
    },
    index=['emodnet (latest)', 'emodnet_2016'],
).round(1)

emodnet2016_depth = fetch_depth('emodnet_2016', north_sea) pd.DataFrame( { 'grid_rows': [emodnet_depth.shape[0], emodnet2016_depth.shape[0]], 'min_depth': [ float(np.nanmin(emodnet_depth)), float(np.nanmin(emodnet2016_depth)), ], 'mean_depth': [ float(np.nanmean(emodnet_depth)), float(np.nanmean(emodnet2016_depth)), ], }, index=['emodnet (latest)', 'emodnet_2016'], ).round(1)

The depths agree closely in this stable shelf area, but the grids differ in resolution — the newer release is finer. Pick a specific release when you need reproducibility against a published DTM version; otherwise emodnet gives you the current best.

The European-domain guard¶

EMODnet is regional, and its WCS server returns a zero-filled grid (not an error) for an area outside coverage. To stop that becoming silent bad data, each row carries its own native_bbox and the backend guards every request against it. Crucially, the older releases cover a smaller domain than the latest:

In [ ]:

releases = ['emodnet', 'emodnet_2022', 'emodnet_2020', 'emodnet_2018', 'emodnet_2016']
pd.DataFrame(
    {'native_bbox (W, S, E, N)': [catalog.get(d).native_bbox for d in releases]},
    index=releases,
)

releases = ['emodnet', 'emodnet_2022', 'emodnet_2020', 'emodnet_2018', 'emodnet_2016'] pd.DataFrame( {'native_bbox (W, S, E, N)': [catalog.get(d).native_bbox for d in releases]}, index=releases, )

The guard has three behaviours, keyed on the request bbox vs the row's native_bbox:

Request vs coverage	Behaviour
Fully inside	proceeds silently
Partially outside	proceeds, but logs a warning that out-of-coverage cells come back as 0.0 fill
Fully outside	raises ValueError, pointing you at the global gebco_2020 / etopo1_ice DEMs

This is why the per-release native_bbox above must be exact: the latest and 2022 run west into the Atlantic (-70.5), while 2016 stops at -36 west and 85°N. Using one extent for all releases would let an older-release request slip into the gap and return silent zeros.

Takeaway¶

• dataset="emodnet" fetches the European high-resolution DTM over OGC WCS through pyramids.Dataset.from_wcs, on the same download() -> GeoTIFF contract as the global griddap DEMs.
• It is ~4× finer than GEBCO over European seas; reach for GEBCO / ETOPO outside the EMODnet domain.
• Year-stamped ids pin older releases, which differ in resolution and extent.
• The domain guard uses each row's native_bbox to reject fully-outside requests and warn on partial overlap, so you never silently map zero-fill as real depth.

Attribution¶

EMODnet Digital Bathymetry (DTM 2024), EMODnet Bathymetry Consortium (doi:10.12770/cf51df64-56f9-4a99-b1aa-36b8d7b743a1); EMODnet use conditions, attribution required. GEBCO Compilation Group (2020) GEBCO 2020 Grid. A DTM is a static grid, so download(aggregate=...) is rejected.