Reading 9 national lidar terrain models over HTTP, without downloading a country

Most countries in Europe now publish a lidar-derived terrain model (DTM) at 1–5 m as open data. On paper that's amazing. In practice, each one sits behind a different web service with its own idea of what a "bounding box" is, what "no data" looks like, and how long you're allowed to wait.

I build TopoLines, a web tool that traces contour lines for a zone you draw on a map. To trace a zone in France on 5 m lidar instead of a 30 m global model, I need that zone's elevations, in a few seconds, without downloading a département. So I wired nine national models plus two global ones, all read over HTTP, all in plain Python (rasterio + numpy, no GDAL CLI, no Docker).

This post is the field notes. The full list — endpoints, native resolution, licence, native CRS — is open on GitHub as Markdown and CSV.

The four shapes of "give me elevations for this box"

Every provider I've met falls into one of four access patterns:

1. A global COG: read just a window

GEDTM30 (OpenGeoHub, CC BY 4.0) is a 30 m bare-earth model of the whole planet, published as a single Cloud-Optimized GeoTIFF. You never download it: rasterio opens it through /vsicurl/ and only fetches the byte ranges for your window.

The GDAL settings matter more than you'd think. Without them GDAL tries to list the bucket directory before opening the file, which is slow and pointless:

import rasterio
from rasterio.windows import from_bounds

URL = ("https://opentopography.s3.sdsc.edu/raster/GEDTM30/GEDTM30_be/"
       "gedtm_rf_m_30m_s_20060101_20151231_go_epsg.4326.3855_v1.2.tif")

env = {
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",   # don't list the bucket
    "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif",    # don't probe sidecar files
    "GDAL_HTTP_TIMEOUT": "20",
    "VSI_CACHE": "TRUE",
}

def gedtm_window(west, south, east, north):
    with rasterio.Env(**env), rasterio.open(f"/vsicurl/{URL}") as ds:
        win = from_bounds(west, south, east, north, ds.transform)
        win = win.round_offsets().round_lengths()
        arr = ds.read(1, window=win).astype("float32")
        arr[arr >= 3.4e38] = float("nan")          # float32-max = sea / void
        return arr, ds.window_transform(win)

The same file is mirrored on OpenGeoHub's own S3, so I try one host and fall through to the other.

2. WMS and ImageServer that return numbers, not pictures

People think of WMS as "the thing that returns map images". But IGN's Géoplateforme will hand you raw float32 elevations as a GeoTIFF if you ask for the right layer:

params = {
    "SERVICE": "WMS", "VERSION": "1.3.0", "REQUEST": "GetMap",
    "LAYERS": "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES",
    "STYLES": "",
    "CRS": "EPSG:4326",
    "BBOX": f"{south},{west},{north},{east}",   # WMS 1.3.0 + 4326 = lat,lon order!
    "WIDTH": width, "HEIGHT": height,
    "FORMAT": "image/geotiff",
}
r = requests.get("https://data.geopf.fr/wms-r/wms", params=params, timeout=30)

Two things bit me here:

• Axis order. WMS 1.3.0 with EPSG:4326 wants lat,lon. Get it wrong and you get a valid GeoTIFF of somewhere else.
• "200 OK" isn't success. WMS can answer 200 with an XML ServiceException in the body. Check the Content-Type before you write anything to disk.

I also measured how big one GetMap can be: 0.8 Mpx answered in 4.7 s, 2.4 Mpx in 7.8 s, and 4.2 Mpx returned a 504 no matter how long I waited. Splitting into tiles works (the seams are invisible), but it multiplies the failure surface. So I cap a request at 3 Mpx — roughly 75 km² at 5 m — and above that I use a coarser source rather than promise 5 m and time out.

USGS 3DEP works the same way through an ArcGIS ImageServer exportImage. The nice surprise: the mosaic picks the best source per pixel, so a 5 m request returns lidar-derived heights wherever lidar has been flown and the 1/3″ (~10 m) product elsewhere. Force your own noData so you know what a void looks like:

params = {
    "bbox": f"{west},{south},{east},{north}", "bboxSR": 4326, "imageSR": 4326,
    "size": f"{width},{height}", "format": "tiff", "pixelType": "F32",
    "noData": -9999, "noDataInterpretation": "esriNoDataMatchAny",
    "interpolation": "RSP_BilinearInterpolation", "f": "image",
}

3. WCS 2.0.1: same standard, four dialects

WCS is the "proper" way to get coverage data, and four of my providers use it. They do not agree on much.

Spain (IGN España) is the polite one: you subset in lat/long, SCALESIZE pins the output grid, the server reprojects for you.

params = [
    ("SERVICE", "WCS"), ("VERSION", "2.0.1"), ("REQUEST", "GetCoverage"),
    ("COVERAGEID", "Elevacion4258_5"),
    ("SUBSET", f"long({west},{east})"),
    ("SUBSET", f"lat({south},{north})"),
    ("SCALESIZE", f"long({width}),lat({height})"),
    ("FORMAT", "image/tiff"),
]

The catch: there is no void value. The sea and the land outside Spain (Portugal, France) both come back as a flat 0. If you trace contours on that, you get a 0 m contour drawn along the Portuguese border. I remap 0 to a void sentinel before anything else.

England (Environment Agency), Italy (INGV TINITALY) and North Rhine-Westphalia only behave in their native projected CRS. On the EA's GeoServer, a lon/lat subset errors or hangs for minutes. Asking in British National Grid eastings/northings with SCALEFACTOR returns a 5 m grid for a 6 × 6 km box in about a second:

from rasterio.warp import transform_bounds

e0, n0, e1, n1 = transform_bounds("EPSG:4326", "EPSG:27700", west, south, east, north)
pad = 50  # metres, so the warp back to lon/lat has data at the edges
params = [
    ("SERVICE", "WCS"), ("VERSION", "2.0.1"), ("REQUEST", "GetCoverage"),
    ("COVERAGEID", COVERAGE_ID),
    ("SUBSET", f"E({e0 - pad:.0f},{e1 + pad:.0f})"),
    ("SUBSET", f"N({n0 - pad:.0f},{n1 + pad:.0f})"),
    ("SCALEFACTOR", f"{1.0 / 5.0:g}"),   # native 1 m → 5 m
    ("FORMAT", "image/tiff"),
]

Then I warp the result onto my lon/lat grid myself. The same function serves England (EPSG:27700), Italy (EPSG:32632 — for the whole country, eastern regions included) and NRW (EPSG:25832).

One more English trap: a subset entirely off the data (Wales, Scotland) also comes back as a flat 0 with no nodata tag. Real lidar ground in England is never exactly 0.0 — even the Fens read −11.42 m — so there, 0 means "no data".

4. Tiles: find them, read them, mosaic them

Switzerland (swisstopo swissALTI3D) publishes 1 km Cloud-Optimized GeoTIFFs behind a STAC API. Ask for items in your bbox, and you'll often get several survey years for the same tile, so keep the newest:

url = ("https://data.geo.admin.ch/api/stac/v0.9/collections/ch.swisstopo.swissalti3d/items"
       f"?bbox={west},{south},{east},{north}&limit=100")
newest = {}
while url:
    data = requests.get(url, timeout=30).json()
    for feat in data["features"]:
        _, year, tile = feat["id"].rsplit("_", 2)          # swissalti3d_<year>_<E>-<N>
        href = next(a["href"] for a in feat["assets"].values()
                    if a["href"].endswith("_2_2056_5728.tif"))   # the 2 m asset
        if int(year) > newest.get(tile, (0, ""))[0]:
            newest[tile] = (int(year), href)
    url = next((l["href"] for l in data["links"] if l["rel"] == "next"), None)

The tiles are kilometre-aligned in LV95, so they mosaic without resampling at the seams. Then I warp once to lon/lat.

Austria (BEV) publishes 50 km COG tiles in EPSG:3035 with internal overviews. Asking rasterio for a window at a coarser out_shape makes GDAL pick the matching overview, so a 5–20 m read only pulls a few MB. There are exactly 55 published tiles; anything else 404s, so I keep the list in code instead of guessing.

Bavaria is the hard one: no open WCS (the INSPIRE one needs a login), only 1 km GeoTIFF tiles that are not cloud-optimized, about 2.5 MB each. Every tile is a full download, so I read them in parallel and cap the count.

And Germany as a whole has no free nationwide DGM1 — the federal one is licensed — so each Land is its own provider with its own service. Baden-Württemberg's INSPIRE WCS answers, but over part of the Land it returns whole-metre integers with no CRS, which terraces every contour. It's not wired yet.

The part nobody warns you about: is that void the sea, or France?

National models stop at the border and at the coast, and both show up as "no data". For contours, these mean opposite things:

• Void = sea → keep it, it's your coastline.
• Void = a neighbouring country → this zone isn't fully covered; use another source.

My heuristic: look at the ring of valid pixels around each void. If the median elevation of that ring is below 60 m, it's a shoreline. If it's higher, the void is foreign land the national model doesn't cover, and the zone falls back to a worldwide source. Above 0.5 % of void pixels inside the box, I don't trust the national model for that zone at all.

This also leads to a rule I'd recommend to anyone mixing sources: one dataset per map. A national model is only used when the zone lies entirely inside its coverage. Stitching a 1 m lidar model to a 30 m global one at a border gives you a visible step in the contours exactly where you least want it.

What I'd tell myself a year ago

1. Check the Content-Type, always. Every OGC server here returns errors as XML or JSON with a 200 status.
2. Ask in the native CRS. Reprojection on the server is where WCS implementations break.
3. Know what "no data" looks like for each source. It's NaN, -9999, -99999, float32-max, or 0, depending on who you ask.
4. Measure the server, then cap below it. Falling back fast is better than a 30 s timeout.
5. Licences differ, attribution doesn't. CC BY 4.0, OGL v3.0, Etalab 2.0 and swisstopo's terms all require a credit; I write it into the metadata of every export.

The full list, with endpoints, native resolution, licence and CRS for each model, is on GitHub: open national elevation models (CSV). If you know an open national DTM with a usable web service that isn't there, open an issue and I'll test it.

And if you just want the contours: TopoLines traces them for any zone on Earth, on whichever of these models covers it, and exports GeoJSON or DXF in your CRS.