Bathymetry — DEM quickstart¶

Fetch a small GEBCO 2020 topography / bathymetry subset over a deep-ocean AOI, map the sea-floor depth, then compare the two ETOPO1 variants (ice surface vs bedrock) over Antarctica to read off ice thickness. Every DEM is subset on the NOAA ERDDAP griddap server and written to GeoTIFF by pyramids — earthlens never touches a competing array stack.

Setup¶

earthlens provides the unified EarthLens entry point; pyramids reads the written GeoTIFF; downloads go to a temporary directory.

In [ ]:

import tempfile

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

from earthlens import EarthLens
from earthlens.bathymetry import Catalog

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

The DEM catalog¶

Three curated DEMs ship today — GEBCO 2020 (15″) and the two ETOPO1 (1′) variants. Reading the catalog is offline.

In [ ]:

catalog = Catalog()
for dataset_id in sorted(catalog.datasets):
    row = catalog.get(dataset_id)
    print(f'{dataset_id:16} {row.native_resolution:14} band={row.variable!r}')

catalog = Catalog() for dataset_id in sorted(catalog.datasets): row = catalog.get(dataset_id) print(f'{dataset_id:16} {row.native_resolution:14} band={row.variable!r}')

Download a GEBCO subset¶

Build the request — the bathymetry source, the gebco_2020 dataset, and a small AOI off the NW African coast (aoi is [min_lon, min_lat, max_lon, max_lat]). download() returns the written GeoTIFF path(s).

In [ ]:

out = tempfile.mkdtemp()
paths = EarthLens(
    data_source='bathymetry',
    dataset='gebco_2020',
    aoi=[-18.0, 25.0, -17.0, 26.0],
    path=out,
).download(progress_bar=False)
paths

out = tempfile.mkdtemp() paths = EarthLens( data_source='bathymetry', dataset='gebco_2020', aoi=[-18.0, 25.0, -17.0, 26.0], path=out, ).download(progress_bar=False) paths

Read and map the sea floor¶

Read the GeoTIFF with pyramids, mask the no-data fill (32767), and map the elevation — every value here is below sea level, so the depths are negative.

In [ ]:

ds = Dataset.read_file(str(paths[0]))
arr = np.asarray(ds.read_array(), dtype='float32')
arr = arr[0] if arr.ndim == 3 else arr
arr = np.where(arr == 32767, np.nan, arr)
print('epsg', ds.epsg, '| shape', arr.shape,
      '| depth range', round(float(np.nanmin(arr))), '..', round(float(np.nanmax(arr))), 'm')

ds = Dataset.read_file(str(paths[0])) arr = np.asarray(ds.read_array(), dtype='float32') arr = arr[0] if arr.ndim == 3 else arr arr = np.where(arr == 32767, np.nan, arr) print('epsg', ds.epsg, '| shape', arr.shape, '| depth range', round(float(np.nanmin(arr))), '..', round(float(np.nanmax(arr))), 'm')

In [ ]:

plt.figure(figsize=(6, 5))
plt.imshow(arr, cmap='Blues_r')
plt.colorbar(label='elevation (m, negative = below sea level)')
plt.title('GEBCO 2020 (15″) — NW African shelf')
plt.axis('off')
plt.show()

plt.figure(figsize=(6, 5)) plt.imshow(arr, cmap='Blues_r') plt.colorbar(label='elevation (m, negative = below sea level)') plt.title('GEBCO 2020 (15″) — NW African shelf') plt.axis('off') plt.show()

ETOPO1: ice surface vs bedrock¶

ETOPO1 ships two variants. Over an Antarctic AOI, the ice surface is the top of the ice sheet and the bedrock is its base; their difference is the ice thickness. We download both and subtract.

In [ ]:

def fetch_etopo(dataset_id):
    out_dir = tempfile.mkdtemp()
    written = EarthLens(
        data_source='etopo',
        dataset=dataset_id,
        aoi=[150.0, -78.0, 170.0, -74.0],
        path=out_dir,
    ).download(progress_bar=False)
    grid = np.asarray(Dataset.read_file(str(written[0])).read_array(), dtype='float32')
    return grid[0] if grid.ndim == 3 else grid


ice = fetch_etopo('etopo1_ice')
bedrock = fetch_etopo('etopo1_bedrock')
thickness = ice - bedrock
print('max ice thickness', round(float(np.nanmax(thickness))), 'm')

def fetch_etopo(dataset_id): out_dir = tempfile.mkdtemp() written = EarthLens( data_source='etopo', dataset=dataset_id, aoi=[150.0, -78.0, 170.0, -74.0], path=out_dir, ).download(progress_bar=False) grid = np.asarray(Dataset.read_file(str(written[0])).read_array(), dtype='float32') return grid[0] if grid.ndim == 3 else grid ice = fetch_etopo('etopo1_ice') bedrock = fetch_etopo('etopo1_bedrock') thickness = ice - bedrock print('max ice thickness', round(float(np.nanmax(thickness))), 'm')

In [ ]:

plt.figure(figsize=(6, 5))
plt.imshow(thickness, cmap='viridis')
plt.colorbar(label='ice thickness (m)')
plt.title('ETOPO1 ice surface − bedrock — East Antarctica')
plt.axis('off')
plt.show()

plt.figure(figsize=(6, 5)) plt.imshow(thickness, cmap='viridis') plt.colorbar(label='ice thickness (m)') plt.title('ETOPO1 ice surface − bedrock — East Antarctica') plt.axis('off') plt.show()

Attribution¶

• GEBCO — GEBCO Compilation Group (2020) GEBCO 2020 Grid.
• ETOPO1 — Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model.

Both are served by the NOAA ERDDAP (U.S. Government public domain). A DEM is a static grid with no time axis, so download(aggregate=...) is rejected.