Copernicus DEM — quickstart (anonymous, no account)¶

Fetch a Copernicus DEM GLO-30 tile over the Nile Delta from the public AWS Open Data bucket, then read + plot it with pyramids. No credentials, no SDK login, no API key — the bucket is anonymous. earthlens.dem just hands you the raw COG; the cropping / mosaicking / hillshade is pyramids' job.

Then a second bbox spanning two land tiles shows the multi-tile behaviour and how to mosaic them into one continuous raster.

Setup¶

earthlens provides the unified EarthLens entry point; pyramids reads the downloaded COG; downloads go to a per-notebook temp directory.

In [1]:

import tempfile
from pathlib import Path

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

from earthlens import EarthLens
from earthlens.dem import Catalog

download_root = Path(tempfile.mkdtemp(prefix='dem-quickstart-'))
print(f'downloads will land under: {download_root}')

import tempfile from pathlib import Path import matplotlib.pyplot as plt import numpy as np from pyramids.dataset import Dataset from earthlens import EarthLens from earthlens.dem import Catalog download_root = Path(tempfile.mkdtemp(prefix='dem-quickstart-')) print(f'downloads will land under: {download_root}')

2026-07-05 14:07:38 | INFO | pyramids.base.config | Logging is configured.

downloads will land under: C:\Users\main\AppData\Local\Temp\dem-quickstart-7u4xer7t

The DEM catalog¶

Two datasets ship today — the ~30 m GLO-30 and the ~90 m GLO-90 grids. Reading the catalog is offline.

In [2]:

catalog = Catalog()
for dataset_id in sorted(catalog.datasets):
    row = catalog.get_dataset(dataset_id)
    print(
        f'{dataset_id:16} bucket={row.bucket:22} '
        f'token={row.resolution_token} ≈{row.native_resolution_m} m'
    )

catalog = Catalog() for dataset_id in sorted(catalog.datasets): row = catalog.get_dataset(dataset_id) print( f'{dataset_id:16} bucket={row.bucket:22} ' f'token={row.resolution_token} ≈{row.native_resolution_m} m' )

cop-dem-glo-30   bucket=copernicus-dem-30m     token=10 ≈30 m
cop-dem-glo-90   bucket=copernicus-dem-90m     token=30 ≈90 m

Download one GLO-30 tile over the Nile Delta¶

The Copernicus DEM grid is 1° x 1° tiles keyed by their SW corner. The bbox here (lat=[30.2, 30.8], lon=[31.2, 31.8]) lies entirely inside the tile at (30 N, 31 E), so one COG is returned.

In [3]:

single_out = download_root / 'nile_delta'
paths = EarthLens(
    data_source='dem',
    dataset='cop-dem-glo-30',
    lat_lim=[30.2, 30.8],
    lon_lim=[31.2, 31.8],
    path=single_out,
).download()

for path in paths:
    print(f'  {path.name}  ({path.stat().st_size / 1024**2:.1f} MB)')

single_out = download_root / 'nile_delta' paths = EarthLens( data_source='dem', dataset='cop-dem-glo-30', lat_lim=[30.2, 30.8], lon_lim=[31.2, 31.8], path=single_out, ).download() for path in paths: print(f' {path.name} ({path.stat().st_size / 1024**2:.1f} MB)')

dem/cop-dem-glo-30:   0%|          | 0/1 [00:00<?, ?tile/s]

dem/cop-dem-glo-30: 100%|██████████| 1/1 [00:00<00:00,  1.23tile/s]

dem/cop-dem-glo-30: 100%|██████████| 1/1 [00:00<00:00,  1.23tile/s]

  Copernicus_DSM_COG_10_N30_00_E031_00_DEM.tif  (47.2 MB)

Read the tile back with pyramids¶

The returned COG carries a WGS84 CRS and one elevation band (metres above the EGM2008 geoid). Read it into a pyramids.Dataset and plot the raw elevation.

In [4]:

dem = Dataset.read_file(str(paths[0]))
print(f'EPSG        : {dem.epsg}')
print(f'shape (y, x): {dem.rows} x {dem.columns}')
print(f'pixel size  : {dem.cell_size}')
print(f'no-data     : {dem.no_data_value}')

elevation = np.asarray(dem.read_array(), dtype=float)
# Mask on the raster's declared no-data sentinel, not a heuristic threshold.
nodata = dem.no_data_value
if nodata is not None and nodata[0] is not None:
    elevation = np.where(elevation == nodata[0], np.nan, elevation)

fig, ax = plt.subplots(figsize=(6, 6))
im = ax.imshow(
    elevation,
    cmap='terrain',
    extent=[31.0, 32.0, 30.0, 31.0],
    origin='upper',
)
ax.set_title('Copernicus DEM GLO-30 — Nile Delta tile')
ax.set_xlabel('longitude (°E)')
ax.set_ylabel('latitude (°N)')
cbar = fig.colorbar(im, ax=ax, shrink=0.8)
cbar.set_label('elevation (m, EGM2008)')
plt.tight_layout()
plt.show()

dem = Dataset.read_file(str(paths[0])) print(f'EPSG : {dem.epsg}') print(f'shape (y, x): {dem.rows} x {dem.columns}') print(f'pixel size : {dem.cell_size}') print(f'no-data : {dem.no_data_value}') elevation = np.asarray(dem.read_array(), dtype=float) # Mask on the raster's declared no-data sentinel, not a heuristic threshold. nodata = dem.no_data_value if nodata is not None and nodata[0] is not None: elevation = np.where(elevation == nodata[0], np.nan, elevation) fig, ax = plt.subplots(figsize=(6, 6)) im = ax.imshow( elevation, cmap='terrain', extent=[31.0, 32.0, 30.0, 31.0], origin='upper', ) ax.set_title('Copernicus DEM GLO-30 — Nile Delta tile') ax.set_xlabel('longitude (°E)') ax.set_ylabel('latitude (°N)') cbar = fig.colorbar(im, ax=ax, shrink=0.8) cbar.set_label('elevation (m, EGM2008)') plt.tight_layout() plt.show()

EPSG        : 4326
shape (y, x): 3600 x 3600
pixel size  : 0.0002777777777777778
no-data     : (None,)

Crop the tile to the exact bbox¶

earthlens returns the whole 1° tile. To keep only the bbox pixels, hand the raster to pyramids.Dataset.crop with the bbox in [west, south, east, north] order.

In [5]:

cropped = dem.crop(bbox=(31.2, 30.2, 31.8, 30.8), epsg=4326)
print(f'cropped shape (y, x): {cropped.rows} x {cropped.columns}')
cropped_array = np.asarray(cropped.read_array(), dtype=float)
nodata = cropped.no_data_value
if nodata is not None and nodata[0] is not None:
    cropped_array = np.where(cropped_array == nodata[0], np.nan, cropped_array)

fig, ax = plt.subplots(figsize=(6, 6))
im = ax.imshow(
    cropped_array,
    cmap='terrain',
    extent=[31.2, 31.8, 30.2, 30.8],
    origin='upper',
)
ax.set_title('Copernicus DEM GLO-30 — cropped bbox')
ax.set_xlabel('longitude (°E)')
ax.set_ylabel('latitude (°N)')
cbar = fig.colorbar(im, ax=ax, shrink=0.8)
cbar.set_label('elevation (m)')
plt.tight_layout()
plt.show()

cropped = dem.crop(bbox=(31.2, 30.2, 31.8, 30.8), epsg=4326) print(f'cropped shape (y, x): {cropped.rows} x {cropped.columns}') cropped_array = np.asarray(cropped.read_array(), dtype=float) nodata = cropped.no_data_value if nodata is not None and nodata[0] is not None: cropped_array = np.where(cropped_array == nodata[0], np.nan, cropped_array) fig, ax = plt.subplots(figsize=(6, 6)) im = ax.imshow( cropped_array, cmap='terrain', extent=[31.2, 31.8, 30.2, 30.8], origin='upper', ) ax.set_title('Copernicus DEM GLO-30 — cropped bbox') ax.set_xlabel('longitude (°E)') ax.set_ylabel('latitude (°N)') cbar = fig.colorbar(im, ax=ax, shrink=0.8) cbar.set_label('elevation (m)') plt.tight_layout() plt.show()

cropped shape (y, x): 3600 x 3600

Multi-tile bbox and mosaic¶

A wider bbox spans several tiles — earthlens returns one COG per intersected tile in row-major order. pyramids.dataset.merge.merge_rasters composites them into one continuous raster.

The bbox below (lat=[45.4, 45.6], lon=[7.4, 8.6]) intersects two Alpine tiles at (45 N, 7 E) and (45 N, 8 E).

In [6]:

alps_out = download_root / 'alps'
alp_paths = EarthLens(
    data_source='dem',
    dataset='cop-dem-glo-90',
    lat_lim=[45.4, 45.6],
    lon_lim=[7.4, 8.6],
    path=alps_out,
).download()

for path in alp_paths:
    print(f'  {path.name}  ({path.stat().st_size / 1024**2:.2f} MB)')

print(f'{len(alp_paths)} tiles fetched')

alps_out = download_root / 'alps' alp_paths = EarthLens( data_source='dem', dataset='cop-dem-glo-90', lat_lim=[45.4, 45.6], lon_lim=[7.4, 8.6], path=alps_out, ).download() for path in alp_paths: print(f' {path.name} ({path.stat().st_size / 1024**2:.2f} MB)') print(f'{len(alp_paths)} tiles fetched')

dem/cop-dem-glo-90:   0%|          | 0/2 [00:00<?, ?tile/s]

dem/cop-dem-glo-90:  50%|█████     | 1/2 [00:00<00:00,  2.54tile/s]

dem/cop-dem-glo-90: 100%|██████████| 2/2 [00:00<00:00,  3.44tile/s]

dem/cop-dem-glo-90: 100%|██████████| 2/2 [00:00<00:00,  3.27tile/s]

  Copernicus_DSM_COG_30_N45_00_E007_00_DEM.tif  (5.00 MB)
  Copernicus_DSM_COG_30_N45_00_E008_00_DEM.tif  (4.90 MB)
2 tiles fetched

In [7]:

from pyramids.dataset.merge import merge_rasters

mosaic_path = alps_out / 'alps_mosaic.tif'
merge_rasters([str(p) for p in alp_paths], str(mosaic_path))

mosaic = Dataset.read_file(str(mosaic_path))
print(f'mosaic shape (y, x): {mosaic.rows} x {mosaic.columns}')
print(f'mosaic bbox        : {mosaic.bbox}')

mb = mosaic.bbox
arr = np.asarray(mosaic.read_array(), dtype=float)
nodata = mosaic.no_data_value
if nodata is not None and nodata[0] is not None:
    arr = np.where(arr == nodata[0], np.nan, arr)

fig, ax = plt.subplots(figsize=(8, 4))
im = ax.imshow(
    arr,
    cmap='terrain',
    extent=[mb[0], mb[2], mb[1], mb[3]],
    origin='upper',
)
ax.set_title('Copernicus DEM GLO-90 — two-tile Alpine mosaic')
ax.set_xlabel('longitude (°E)')
ax.set_ylabel('latitude (°N)')
cbar = fig.colorbar(im, ax=ax, shrink=0.8)
cbar.set_label('elevation (m, EGM2008)')
plt.tight_layout()
plt.show()

from pyramids.dataset.merge import merge_rasters mosaic_path = alps_out / 'alps_mosaic.tif' merge_rasters([str(p) for p in alp_paths], str(mosaic_path)) mosaic = Dataset.read_file(str(mosaic_path)) print(f'mosaic shape (y, x): {mosaic.rows} x {mosaic.columns}') print(f'mosaic bbox : {mosaic.bbox}') mb = mosaic.bbox arr = np.asarray(mosaic.read_array(), dtype=float) nodata = mosaic.no_data_value if nodata is not None and nodata[0] is not None: arr = np.where(arr == nodata[0], np.nan, arr) fig, ax = plt.subplots(figsize=(8, 4)) im = ax.imshow( arr, cmap='terrain', extent=[mb[0], mb[2], mb[1], mb[3]], origin='upper', ) ax.set_title('Copernicus DEM GLO-90 — two-tile Alpine mosaic') ax.set_xlabel('longitude (°E)') ax.set_ylabel('latitude (°N)') cbar = fig.colorbar(im, ax=ax, shrink=0.8) cbar.set_label('elevation (m, EGM2008)') plt.tight_layout() plt.show()

mosaic shape (y, x): 1200 x 2400
mosaic bbox        : [6.999583333333334, 45.000416666666666, 8.999583333333334, 46.000416666666666]

What just happened — the account-free path¶

Every step above spoke to s3://copernicus-dem-30m and s3://copernicus-dem-90m anonymously: no ~/.aws/credentials, no AWS profile, no environment variable, no signer. Copernicus DEM is the only globally-consistent DEM reachable through earthlens without an account, which is the entire justification for the dem backend.