Derived Terrain Protocol

demdiff:// (lib/demdiff-protocol.ts) is the engine behind the Difference of two sources terrain type — the feature side is documented in nDSM and Comparison. It takes two elevation sources and emits a third that is their difference, as an ordinary Terrain-RGB tile, so every downstream consumer (3D terrain, hillshade, hypsometric tint, slope, contours, the elevation picker, the GeoTIFF export) treats the result exactly like any other DEM.

The URL

demdiff://<encA>/<encB>/<tileSize>/<offsetM>/<template A>/<template B>/{z}/{x}/{y}

Both operand templates are percent-encoded inside the URL, so MapLibre's literal {z}/{x}/{y} substitution only touches the trailing coordinates and never the embedded ones. buildDemDiffUrl() assembles it.

Because the whole configuration lives in the URL, the tile result cache keys on it for free: change the offset, or either operand, and it is a different key — a stale difference can never be served.

Why both operands are fetched at the same z/x/y

The two sources are fetched through fetchDecodedTile — the same decoded-tile cache the viz-mode protocols use, so cog://, titiler, plain TMS and float32dem:// all work as operands — and at the same tile coordinate. That is what makes the subtraction meaningful: two tiles at the same z/x/y cover the same ground and the same pixel grid, whatever the sources' native resolutions. No reprojection, no resampling to a common grid, no alignment step.

Ancestor-tile upsampling

The interesting case is a 0.5 m DSM minus a 30 m reference. Above the coarse source's maxzoom there is no tile to fetch, and the naive result is that the difference simply stops existing exactly where it is most useful.

fetchOperand walks up instead: it tries z, then z-1, z-2, … up to 6 levels, and returns the first ancestor that exists together with the sub-window to read:

type Operand = { tile: DecodedTile; scale: number; ox: number; oy: number }

sampleOperand then reads that window bilinearly, so the coarse side is smoothly upsampled on the fly rather than blocking into 30 m steps. Six levels is 64× — enough to put a GLO-30 reference under a half-metre drone survey.

Holes, and why they are alpha 254

A pixel is a hole when either operand is missing it — nodata, masked, outside coverage, or nothing found within six zoom levels. The obvious encoding is "write transparent and let the consumer work it out", and it is wrong in a way that is worth stating plainly, because the same trap is laid in three different places in this app.

MapLibre's DEM decoder ignores alpha, and the browser premultiplies it. A transparent pixel is premultiplied to RGB 0,0,0 on decode, which in Terrain-RGB is the floor, −10 000 m. A single transparent pixel is a 10 km pit; a nodata edge is a 10 km cliff along the whole boundary.

So a hole is written as 0 m in the colour channels with alpha 254:

const hole = !(Number.isFinite(va) && Number.isFinite(vb))
const [r, g, bl] = elevationToTerrainrgb(hole ? 0 : va - vb + offset)
out[i] = r; out[i + 1] = g; out[i + 2] = bl; out[i + 3] = hole ? 254 : 255

MapLibre reads the colour channels, sees flat ground and draws nothing dramatic. This app's own decoders read the alpha and know the cell is invalid. 254 is chosen because it is fully opaque for every practical purpose while still being distinguishable from the 255 that every real pixel carries.

The same convention, three places

If you add a fourth path that produces elevation tiles, it has to honour this or it will either dig pits or invent ground.

Titiler does it server-side

Titiler-backed sources do not need the convention at all, because titiler's terrainrgb algorithm can encode masked pixels at a chosen height:

&algorithm_params={"nodata_height":0}&return_mask=false

This is the same flattening, done once on the server, for free — it replaced an earlier client-side demfix:// wrapper that decoded and re-encoded every tile in the browser. Dropping the now-pointless alpha channel also makes the tiles about 9% smaller. See TITILER_FLAT_NODATA in lib/source-builder.ts.

The client-side viz modes and demdiff:// itself still request the mask (forClientDecode: true), because they read alpha as a validity flag and are decoding the tile anyway.

Zoom range

A derived source's zoom range is not its own — it is derived from the operands, and the ceiling follows the finer operand + 2. Overzooming past the finer source by a couple of levels is what lets 3D terrain stay smooth when you push in; going further would only invent bilinear detail on both sides at once.

A titiler-served operand contributes no zoom range of its own, which once silently capped a difference at the coarse side's maxzoom — worth remembering when a difference mysteriously stops refining.

Output

The per-pixel loop writes into a Uint8ClampedArray, which the shared toTileImage tail (lib/tile-image.ts) turns into an ImageBitmap. That is what MapLibre receives — a bog-standard Terrain-RGB tile with no trace of how it was made. It is handed over decoded rather than PNG-encoded, which is worth roughly 99 ms a tile; see Tile Caches.

The offset, and how it is measured

diffOffsetM is added to every difference before encoding. It exists because a DSM minus a DTM from two producers, two dates, or an ellipsoidal source minus an orthometric one rarely centres on zero: there is a local constant between them — a datum, a co-registration bias, a different notion of "ground".

The Difference offset row in the Elevation Color section (only shown when the source on screen is a difference) both holds that number and can measure it. It lives there rather than beside the source's name because the hypsometric ramp is where the offset is visible: a diverging ramp centred on zero is what it exists to make honest.

It is not a least-squares fit and not an SSD minimiser. It is an interquartile mean of the raw per-pixel difference:

1. Take the tiles under the viewport at min(zoom, 14), halving the zoom until there are at most 16 of them. The z14 cap keeps a WMS operand from being asked for a screenful of z18 GetMaps just for a statistic.
2. For each tile, fetch both operands through the protocol's own fetchOperand — so an operand with no tile at that zoom is read from its nearest ancestor, up to 6 levels. This matters more than it sounds: Mapterhorn declares maxzoom: 18 but falls back to Copernicus GLO-30 over most of the world and 404s from z13 there, so sampling it directly returned nothing at all against a fine COG.
3. Sample a 24 × 24 grid per tile at fractional positions — fractional so operands with different tile sizes still line up — skipping anything the validity mask or the DEM sentinel guards reject.
4. Sort, drop the top and bottom quartiles, take the mean of the middle half. The offset is the negation of that, since the protocol adds it.

The interquartile trim is the whole point. Buildings, canopy, a quarry — the very things a difference exists to show — sit in the tails, so a plain mean would chase them and a median alone would throw away the agreeing half. The toast reports the sample count, the tile count, the zoom and the median alongside the result, so a suspicious answer is visible as one.

The size of the offset tells you what it is

A useful diagnostic, borrowed from uw-cryo/groundcontrol's vertical-datum guide: a measured offset of geoid magnitude means the datum assumption is wrong, not that the sources are misregistered. Co-registration bias between two surveys of the same ground is metres. Ellipsoidal against orthometric is tens of metres, and up to ±100 m — −105 m south of India, +85 m around Iceland. If the button reports 40 m over the Himalaya, the answer is not to accept it; it is that one source is on the ellipsoid and the other is not.

The library ships both geoid undulation grids (EGM96 15′ and EGM2008 2.5′) as ordinary terrain sources, so that conversion is itself a difference: subtract the geoid from an ellipsoidal DEM and the result is orthometric.

Worth knowing for the sources that are explicitly ellipsoidal: PGC's ArcticDEM and REMA, Re:Earth Terrain, and — when access to it ever arrives — Vantor Precision3D, which is stated on WGS84 G1674 specifically rather than "WGS84" generally. The WGS84 ensemble is about 2 m of deliberate ambiguity across its realizations, which is why a source's realization is worth recording where a producer states one.