Encoding elevation into pixels: DEM terrain tiles in Rust

Modern web maps are expected to do 3D terrain. Hillshade, contour lines, the tilted camera flying over mountains — all of it is powered by something delightfully dumb: PNG images where the color of each pixel is the elevation.

I recently shipped native DEM tile encoding in tileserver-rs (#1220, closing #1008), and the whole thing fits in one Rust file. But that file is full of decisions that will silently produce garbage terrain if you get them wrong, and none of them are documented in one place. So here's that place.

Elevation as RGB

The trick behind raster-dem sources is packing a floating-point elevation into 24 bits of RGB. There are two conventions, and MapLibre GL JS decodes both:

Terrarium: elevation = R*256 + G + B/256 - 32768, precision 1/256 m
Mapbox-RGB: elevation = -10000 + (R*65536 + G*256 + B) * 0.1, precision 0.1 m

Encoding is the reverse: scale the elevation into integer units, split into three bytes. Here's the Terrarium encoder from sources/dem.rs:

/// Terrarium base offset in metres: decoded = packed - 32768.
const TERRARIUM_BASE: f64 = 32768.0;
/// Largest packed value representable in 24 bits (2^24 - 1).
const MAX_U24: f64 = 16_777_215.0;

#[must_use]
pub fn encode_terrarium(elevation: f64) -> [u8; 3] {
    // Work in 1/256-m units so the blue byte is an exact integer; round to
    // nearest before splitting (the round-before-truncate rule), then clamp
    // to the 24-bit packed range [0, 256^3 - 1].
    let units = ((elevation + TERRARIUM_BASE) * 256.0).round();
    let packed = units.clamp(0.0, MAX_U24) as u32;
    let r = (packed >> 16) & 0xFF;
    let g = (packed >> 8) & 0xFF;
    let b = packed & 0xFF;
    [r as u8, g as u8, b as u8]
}

Twelve lines. The interesting part is the two functions hiding inside them: .round() and .clamp().

The round-before-truncate rule

The classic bug in this space comes from rio-rgbify and its descendants: scale the elevation, then let the integer cast truncate it. Floating-point arithmetic being what it is, (1000.5 + 32768.0) * 256.0 can land at 8644735.999999..., and truncation eats a full least-significant bit. Your terrain is now systematically one elevation-interval low, and nobody notices until someone compares contour lines against a survey.

MapLibre's own packDEMData rounds before splitting, for exactly this reason. So does this encoder, and there's a regression test pinning the exact case:

#[test]
fn mapbox_round_before_truncate() {
    // The classic rgbify bug: 1000.5 m packs to exactly 110005.0; a
    // truncating encoder that computed 110004.999.. would lose the LSB.
    // Round-then-split must land on the exact triplet.
    let rgb = encode_mapbox(1000.5);
    let decoded = decode_mapbox(rgb);
    assert!(
        (decoded - 1000.5).abs() < 0.001,
        "mapbox round-before-truncate 1000.5 → {rgb:?} → {decoded}"
    );
}

The clamp matters too. An absurd input — say a corrupt raster cell claiming 40,000 m — must saturate at the 24-bit ceiling, not wrap around a byte and come back as a plausible-looking valley:

#[test]
fn terrarium_clamps_out_of_range() {
    assert_eq!(encode_terrarium(-40000.0), [0, 0, 0], "terrarium underflow");
    assert_eq!(encode_terrarium(40000.0), [255, 255, 255], "terrarium overflow");
}

The alpha channel will betray you

Here's the gotcha I want you to remember if you remember nothing else. Elevation rasters have nodata — ocean, void fills, the edge of coverage. The obvious move is to write those pixels as transparent and call it a day.

MapLibre does not care about your transparency. Its DEMData constructor reads R, G, B and ignores the alpha channel entirely. A transparent pixel still decodes its RGB to an elevation, and if you wrote (0, 0, 0, 0) in Mapbox-RGB encoding, that decodes to −10,000 m. Congratulations: your coastline now has a trench deeper than the Mariana next to every beach.

The fix is a sentinel RGB — a color that decodes to a deliberately out-of-band elevation:

/// Default Terrarium nodata sentinel: (0, 0, 0) decodes to −32768 m.
const TERRARIUM_NODATA: [u8; 3] = [0, 0, 0];
/// Default Mapbox-RGB nodata sentinel: (1, 134, 160) decodes to 0 m
/// (sea level), the Mapbox convention for "no data / ocean".
const MAPBOX_NODATA: [u8; 3] = [1, 134, 160];

That (1, 134, 160) is worth decoding by hand once: 1*65536 + 134*256 + 160 = 100000, and -10000 + 100000 * 0.1 = 0. Sea level. Missing data renders as flat ocean, which is almost always what you want. The encoder still writes alpha 0 on nodata pixels — but as a courtesy to clients that do honor alpha, not as the mechanism.

The pixel loop treats three things as nodata: the band's declared nodata value, NaN, and infinity (GDAL masked reads produce non-finite values):

#[must_use]
pub fn encode_pixels(values: &[f64], params: &EncodeParams) -> Vec<u8> {
    let mut out = Vec::with_capacity(values.len() * 4);
    for &raw in values {
        let is_nodata = !raw.is_finite()
            || params
                .nodata_value
                .is_some_and(|nd| (raw - nd).abs() < f64::EPSILON);
        if is_nodata {
            out.extend_from_slice(&params.nodata_rgba);
        } else {
            let elevation = raw * params.scale + params.offset;
            let [r, g, b] = encode_elevation(elevation, params.encoding);
            out.extend_from_slice(&[r, g, b, 255]);
        }
    }
    out
}

Note scale and offset in there: if your DEM is stored in feet (looking at you, some USGS products), scale = 0.3048 fixes it at encode time instead of requiring a raster rewrite.

This function is deliberately pure — no GDAL types, just &[f64] in and RGBA bytes out. That's what makes every branch unit-testable without a native library in the loop, which is how the module carries 17 unit tests covering sea level, Everest (8848.86 m), the Mariana Trench floor (−10994 m), clamps, NaN, and the feet-to-meters path.

From GeoTIFF to tile

The serving path is GDAL doing what GDAL does: warp the source raster to the tile's Web Mercator bounding box, read floats, encode, PNG:

let px = (maxx - minx) / f64::from(tile_size);
let py = (maxy - miny) / f64::from(tile_size);
warped.set_geo_transform(&[minx, px, 0.0, maxy, 0.0, -py])?;
warped.set_spatial_ref(&web_mercator)?;
reproject(dataset, &warped)?;

let buffer: Buffer<f64> = out_band.read_as::<f64>(
    (0, 0),
    (tile_size as usize, tile_size as usize),
    (tile_size as usize, tile_size as usize),
    Some(resampling),
)?;

let rgba = encode_pixels(buffer.data(), &params);

One decision here saved real precision: the pipeline reads Buffer<f64> straight from GDAL and encodes float-to-RGB directly. No intermediate 8-bit PNG, no quantize-then-requantize round trip. The float leaves the raster exactly once.

The config side wraps any existing raster source, so an operator who already serves a Copernicus COG adds terrain with four lines of TOML:

[[sources]]
id = "terrain"
type = "dem"
input_source = "copernicus-dem" # reuse the already-open COG
dem_encoding = "terrarium"

input_source shares the referenced source's cached GDAL dataset handle — no second file open, no second block cache. The alternative was every operator pre-baking terrain tiles with rio terrain or gdaldem and doubling their storage, or running a second tile server just for elevation. We already had GDAL in the building; the encoder is a couple hundred lines.

"The tests pass" is not the finish line

The part of this feature I'm most confident in isn't the unit tests — it's the QA that happened after them, because encoding correctness and served-tile correctness are different claims. From the PR:

Built with --features dem, served both encodings, fetched real tiles over HTTP: Terrarium 173 KB, Mapbox-RGB 130 KB, both valid PNG.
Decoded the served tiles back to elevations and checked them against the fixture's known range: got 0–818.6 m against a truth of −41 to 1041 m. The decode-what-you-served step is the one that catches a bug the encoder's own round-trip tests structurally cannot — anything that happens between encode_pixels and the HTTP response.
Drove MapLibre GL JS in a real browser with the server as a raster-dem hillshade source: 20+ tiles at HTTP 200, actual relief on the actual screen, zero console errors.

A byte-exact encoder proves your math. Decoding your own served tiles proves your pipeline. A browser rendering hillshade proves the whole claim. They're three different tests, and terrain has enough silent-garbage failure modes that I now consider all three mandatory for anything raster-shaped.

For throughput nerds: the encoder benches at roughly a nanosecond per pixel — a full 256×256 tile encodes at 110–195 Melem/s — so the cost of a DEM tile is dominated by the GDAL warp, not the encoding.

If you want to see the result instead of reading about it, there's a live demo of Everest in 3D terrain — that's a Copernicus DEM going through this exact encoder, rendered by MapLibre with the camera pitched at the Himalayas.

If you're serving elevation from COGs and want terrain without another service in the stack, the feature ships in tileserver-rs v2.34.0+. And if you're writing your own encoder: round before you truncate, clamp before you cast, and never, ever trust the alpha channel.