3D terrain was introduced in Mapbox GL JS in v2 (announcement), providing a seamless digital elevation basemap.

I’ve seen some questions about how to: a.) display custom terrain in Mapbox, and b.) how to augment the existing Mapbox terrain with your own elevation data. I outline methods for approximating both below.

Choose an area of interest

This process works best if you choose a single XYZ tile as an area of interest. Eventually, we will swap out this tile, and all of this tile’s children tiles, for our augmented elevation tiles. In this case, I’ve chosen XYZ tile [70, 57, 7], which corresponds to the quadkey 1222112. You can explore the XYZ tile/quadkey coordinate space and find a tile of your choice in Mapbox’s What the Tile service:

Download the tile from AWS Terrain Tiles

We’ll augment the elevation in an existing elevation tile from AWS Terrain Tiles, which, conveniently, are stored according to XYZ coordinates:

aws s3 cp s3://elevation-tiles-prod/v2/geotiff/7/70/57.tif <LOCAL PATH>

Create an augmentation image

We will arithmetically add the values of an image to the existing elevation. In short, we will create an 8-bit image, where dark values are closer to 0 and light values are closer to the maximum 8-bit value, 255. So, dark values in our augmentation image will contribute less new elevation to the new layer than light values. I used a simple clip art figure, but you could just as easily use a more realistic custom elevation model, for example, captured with your own LiDAR sensor.

Below are the original augmentation image, and a blurred version of that image (saved as body_blur.png), which will serve to feather the newly added elevation into the existing elevation. In this case, I used GIMP software to blur the image. If you plan ahead to make the dimensions of the augmentation image a square to scale with the original elevation image, it will nicely overlay when the time comes.

Georeference the augmentation image

I used the Freehand Raster Georeferencer plugin for QGIS to geographically locate the augmentation image over the original elevation tile. My process consisted of loading the georeferenced elevation tile, then moving and scaling the augmentation image approximately overtop using the georeferencer plugin. I saved the georeferenced file as body_blur.tif.

Augment the original elevation

I used the QGIS raster calculator to combine the existing elevation raster (the file was called “57.tif” which corresponds to 57@1 in the raster calculator) with the augmentation image (called “body_blur@1” in the raster calculator). The formula I entered in the raster calculator was:

"57@1" + ("body_blur@1" * 4)

This formula translates to: for each pixel in the first band of the raster 57, add the corresponding pixel from the first band of the raster body_blur multiplied by 4. In this case, the values of body_blur range from 0-255, and I wanted to stretch the output to about a maximum of 1000m, so I multiplied by 4. A pure white pixel (255) will add 1020m to the output, while a pure black pixel (0), will add no elevation. The combined image, saved as dem_plus_body.tif, is shown below.

Encode elevation values to Terrarium format

The ultimate goal here is to display our augmented elevation layer as a Mapbox GL JS terrain layer, which accepts tiles in either mapbox or terrarium format. Because AWS Terrain Tiles are already provided in terrarium format, we will convert our combined augmented layer to terrarium format, using the Python code below:

with rasterio.open('dem_plus_body.tif') as src:
    dem = src.read(1)
    r = np.zeros(dem.shape)
    g = np.zeros(dem.shape)
    b = np.zeros(dem.shape)

    v = dem + 32768
    r += np.floor(v / 256.0)
    g += np.floor(v % 256.0)
    b += np.floor((v - np.floor(v)) * 256.0)

    meta = src.profile
    meta["dtype"] = rasterio.uint8
    meta["nodata"] = None
    meta["count"] = 3

    with rasterio.open('encoded_elevation.tif', 'w', **meta) as dst:
        dst.write_band(1, r.astype(rasterio.uint8))
        dst.write_band(2, g.astype(rasterio.uint8))
        dst.write_band(3, b.astype(rasterio.uint8))

Run Generate XYZ Tiles tools

I used the Generate XYZ Tiles (Directory) tool in QGIS to create an organized tile directory, using encoded_elevation.tif created above. If you intend to match with AWS Terrain Tiles, create tiles for zoom levels 0 – 15.

Make the tiles available somewhere

I uploaded my tiles to AWS S3, but you can host them wherever you like. I will leave it as an exercise for the reader to sort out a secure permission configuration, but of course, the web map will need permissions to read the data.

Display the augmented tiles

Finally, we can load our tiles in a web map! It’s fairly easy in Mapbox GL JS v2, simply add an entry in the tiles parameter of a raster-dem source, and call setTerrain referencing that source.

map.addSource("augmented-dem", {
     type: "raster-dem",
     tiles: [
     "https://<BUCKET>.s3.amazonaws.com/<PREFIX>/{z}/{x}/{y}.png",
     ],
     tileSize: 256,
     maxzoom: 14,
     encoding: "terrarium",
 });
 map.setTerrain({ source: "augmented-dem", exaggeration: 10 });

(Optional) Create Proxy API

The method above will load only the tiles from your new augmented layer, which is great. However, it gets a little tricky if you want to combine terrain sources. For example, it is likely that you will want to display your augmented elevation tiles AND original elevation tiles for the rest of the world.

One option is to add multiple tile sources to the tiles parameter in the raster-dem source (i.e. append a tile source that references AWS Terrain Tiles terrarium product), but this seems to result in a race condition between the tile sources, so you can’t be sure which tile will be loaded.

A more reliable, albeit potentially more costly, solution is to create an API that handles the tile-source decision for you. In this scenario, you create an API that proxies requests to a microservice (e.g. a Lambda function) that fetches and returns augmented tiles if they fall within a given quadkey, and returns an original elevation tile from AWS Terrain Tiles if not.

I use a handler like below in a Lambda function to swap augmented tiles into AWS Terrain Tiles.

s3 = boto3.client("s3")
quadkey = "<YOUR_QUADKEY>"
tile_bucket = "YOUR_BUCKET"
tile_prefix = "<TILES_PREFIX>/{z}/{x}/{y}"

def lambda_handler(event, context):
    path = event["path"]
    splitpath = path.split("/")
    x = splitpath[-2]
    y = splitpath[-1]
    z = splitpath[-3]

    t = mercantile.Tile(int(x), int(y.split(".")[0]), int(z))

    if str(mercantile.quadkey(t)).startswith(quadkey):
        bucket = tile_bucket
        key = tile_prefix
    else:
        bucket = "elevation-tiles-prod"
        key = f"v2/terrarium/{z}/{x}/{y}"

    response = s3.get_object(
        Bucket=bucket,
        Key=key
    )

    image = response["Body"].read()

    return {
        "headers": {
            "Content-Type": "image/png", 
            "Access-Control-Allow-Origin" : "*",
            "Access-Control-Allow-Methods": "OPTIONS,GET",
            "Access-Control-Allow-Headers": "Content-Type"
        },
        "statusCode": 200,
        "body": base64.b64encode(image).decode("utf-8"),
        "isBase64Encoded": True
    }