How to draw elevation profile for your path

The MapTiler Terrain-RGB tileset contains elevations encoded in a simple png image format as RGB values. You might have wondered if it is possible to access the elevation values and how to get the elevation profile in an HTML/javascript application. It is possible indeed! However, it is not straightforward and one can get sidetracked by many GIS-related challenges. This article attempts to offer guidance and outline steps that need to be taken in order to get an elevation profile from the Terrain-RGB tileset.

The Interactive Elevation Profile Tool

The application we are going to build is a web application that allows us to draw a path and retrieve elevation profiles along that path

The process

Before we dive into the code, it is beneficial to understand the process of retrieving elevation from Terrain-RGB tiles. Let’s go through it now.

Draw a path on the map

The first thing we need to do is to create a simple web app with a map and add the ability to draw paths. We will use maplibre-gl-js library along with maplibre-gl-draw plugin which allows us to do it. By the path, we mean a linestring with several definition points.

Path sampling

In the next step, we need to generate a set of points along the line. Why is this necessary? Imagine that the user draws a line with only two definition points (red dots). Without sampling, the elevation profile would be just a line connecting two elevations.

Later, we will retrieve the elevation for each point. Distance between points should be regular so that we get a nice and smooth elevation profile but it should not be smaller than Terrain-RGB pixel size.

Retrieving Tiles

In order to access the elevation data, we need to fetch Terrain-RGB tiles from the server. Which tiles? Map tiles are organized in a pyramid indexed by (x, y, zoom) indices where the top contains a single tile that covers the entire world, the next level contains 2 x 2 tiles, n-th level contains 2^zoom x 2^zoom tiles etc. We also need a projection that will convert geographic coordinates (latitude, longitude) to meters. We will use Web Mercator projection.

For the given geographic coordinate (latitude, longitude), we need to find 3 values: (x, y, zoom).

zoom: we will use MapTiler Terrain-RGB tileset whose most detailed zoom level is 12. We expect users will draw short profiles so we will the most detailed data. (For long profiles spanning many kilometers, it would be better to use lower resolution data though).

x,y: we will calculate x and y coordinates using the globalmaptiles library. This library allows us to translate geographic coordinates to the Web Mercator coordinates and find the pixel coordinates for a given zoom level. The conversion process is as follows

(latitude, longitude) [degrees]
  -> (x, y) [meters in web mercator]
  -> (px, py) [pixels at a given zoom level]
  -> (tx, ty) [tile index]

Once we have the tile index, we can download the tile from MapTiler cloud.

Finding relevant pixels

We have downloaded the tile successfully. The tile is a png image that we can decode to a byte array where pixels are organized into rows and columns, each represented as an RGBA value. The byte array has origin in top-left corner and there are 512 rows and 512 columns.

The next task is to find the correct pixel in the byte array which covers the requested geographic coordinates. We are working on zoom level 12 and our tile size is 512 pixels. With these parameters, we calculate:

• The pixel coordinates (px, py) for the requested geographic coordinate

• The pixel coordinate (tx, ty) for the tile lower-left corner

• The offsets (x-offset, y-offset) to retrieve the pixel and decode the elevation

Decoding the elevation

Once we have offsets within the RGBA byte array, we can retrieve the 3 values (we don't need A) and decode the elevation:

elevation = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)

Draw the chart

When we determined elevation for each sample along the path, the last step is to draw a chart so that the user can see the profile. We will use an open-source Chartist library.

The implementation

Let’s dive into the code. We are going to explain the most important parts. The full source code is available on GitHub.

The app skeleton

The application user interface is very simple - just a container with two placeholders - one for the map and one for the elevation profile chart. The challenge is to make sure that the map and the chart fill browser window height, do not overflow, and resizes correctly when the user resizes the window. We use CSS Flex Box Layout to build this layout.

Here is the HTML fragment which implements the layout:

<div id="container">
    <div id="map"></div>
    <div id="chart">
        <p>Please draw a line on the map to get the elevation profile.</p>
    </div>
</div>

And the corresponding CSS stylesheet

html { height: 100%; }
body { height: 100%; margin: 0; padding: 0;}
#container {display: flex; flex-flow: column; height: 100%;}
#map { flex: 1 1 65%;}
#chart { flex: 1 1 35%;}
#chart p { text-align: center; vertical-align: center; }

Drawing the path

After we prepared the layout, we can create a new instance of the MapLibre map control and tell it to use #map placeholder.

this.map = new maplibre.Map({
    container: 'map',
    style: `https://api.maptiler.com/maps/outdoor/style.json?key=${this.apiKey}`,
    center: [10.988677069124009, 46.88158715828973],
    zoom: 12
});

In order to draw paths, we are going to use maplibre-gl-draw plugin. The plugin can do many things but we will use it just for drawing paths (line strings). We create the plugin instance and add it to the collection of MapLibre map controls.

this.draw = new MaplibreDraw({
    displayControlsDefault: false,
    controls: {
        line_string: true,
        trash: true
    },
    defaultMode: 'draw_line_string'
});
this.map.addControl(this.draw);

There are a few events the plugin eventually triggers. For us the most important one is draw.create event which is triggered when the user finished the path. The callback which handles the event receives one argument and it contains the GeoJSON representation of the path the user created. We only need coordinates which we can access as follows.

this.map.on('draw.create', async (e) => {
    const feature = e.features[0];
    const coordinates = feature.geometry.coordinates;
    await this.showChart(coordinates);
});

Generating sampling points

As mentioned in the intro section, we need to generate a set of points along the path. This is easy thanks to the Turf.js library. But before we calculate these points we have to determine the step size. To display the elevation profile in the browser window 200 samples should be sufficient. So we will simply calculate the length of the path and divide it by the number of samples. We will also calculate the pixel size for the given latitude and will make sure the step size is not smaller than the pixel size, otherwise, our profile would have couple of samples with equal elevation.

const options = {units: 'meters'};
const line = turf.lineString(coordinates);
const lineLength = turf.length(line, options);
let numSamples = 200;
const metersPerPx = this.getZoomLevelResolution(coordinates[0][1], 12);

const stepSize = Math.max(metersPerPx, lineLength / numSamples);
numSamples = lineLength / stepSize;

To generate points, we use turf.along method in a simple loop.

const samples = [];
for (let i = 0; i <= numSamples; i++) {
    const along = turf.along(line, i * stepSize, options);
    const coords = along.geometry.coordinates;
    samples.push(coords);
}

Downloading tiles

In the GitHub sample, the code for downloading tiles is written in the ElevationProvider class. The class uses GlobalMercator class for coordinate system and projection calculations. Terrain-RGB tiles have 512x512 pixels so we pass this size into the GlobalMercator constructor. We will use a lazy loading pattern to download the tile from the server once we encounter a first point that intersects that tile and we cache its data to tiles dictionary so that we don't have to download again for the following points intersecting that tile.

class ElevationProvider {

    constructor(apiKey) {
        this.apiKey = apiKey;
        this.gm = new GlobalMercator(512);
        this.tiles = {}
    }
    ...
}

Tile index is calculated simply by converting latitude, longitude to meters using Web Mercator projection, then converting meters to pixels at the given zoom level. GlobalMercator class has a handy LatLonToTile method for it.

getTileIndex(lat, lon, zoom) {
    const tms = this.gm.LatLonToTile(lat, lon, zoom);
    const google = this.gm.GoogleTile(tms.tx, tms.ty, zoom);
    return {
        x: google.tx,
        y: google.ty,
        zoom: zoom
    };
}

With the tile index, we can download the tile from MapTiler cloud and read the tile data.

async fetchTile(tileIndex) {
    const url = `https://api.maptiler.com/tiles/terrain-rgb-v2/${tileIndex.zoom}/${tileIndex.x}/${tileIndex.y}.png?key=${this.apiKey}`
    const image = await this.loadImage(url);
    return this.getImageData(image);
}

We will load the tile using Image class. The reason is that the tile is a PNG image and we want to use the browser API to decode it and access raw bytes. So we will package the code for downloading into the Promise so that we can call it from our asynchronous function.

loadImage(url) {
    return new Promise((resolve, reject) => {
        const img = new Image();
        img.crossOrigin = "Anonymous"
        img.addEventListener('load', () => resolve(img));
        img.addEventListener('error', reject);
        img.src = url;
    });
}

After the Image has been loaded, we can use the browser’s HTMLCanvasElement and CanvasRenderingContext2D to turn the image into an array of raw bytes. The most important method is context.getImageData which will return the 2D array of RGBA values.

getImageData(image) {
    const canvas = document.createElement('canvas')
    canvas.width = image.width
    canvas.height = image.height
    const context = canvas.getContext('2d')
    context.drawImage(image, 0, 0)
    return context.getImageData(0, 0, image.width, image.height);
}

Extracting pixels

Now we need to calculate offsets as described in the intro section. We will again use GlobalMercator class to calculate tile extent in meters (TileBoundsmethod) and then covert it to pixels (MetersToPixels method)

getTileExtentPixels(x, y, zoom) {
    const tms = this.gm.TMSTile(x, y, zoom);
    const tileBounds = this.gm.TileBounds(tms.tx, tms.ty, zoom);
    return {
        lowerLeft: this.gm.MetersToPixels(tileBounds.minx, tileBounds.miny, zoom),
        upperRight: this.gm.MetersToPixels(tileBounds.maxx, tileBounds.maxy, zoom)
    }
}

Similarly, we calculate the pixel coordinates of the sample point. With pixel coordinates of both the sample point and the tile extent, we can evaluate indices that we will use to locate the correct but chunk in the tile image data. The “noise“ in the code (min and max methods) just ensures that we always get the correct index in interval <0,511) even for corner cases caused by rounding.

const meters = this.gm.LatLonToMeters(lat, lon);
const pixels = this.gm.MetersToPixels(meters.mx, meters.my, tileIndex.zoom);

const tilePixelExtent = this.getTileExtentPixels(tileIndex.x, tileIndex.y, tileIndex.zoom);

let xOffset = Math.floor(pixels.px - tilePixelExtent.lowerLeft.px);
xOffset = Math.max(0, Math.min(tileData.width- 1, xOffset));
let yOffset = tileData.height - Math.floor(pixels.py - tilePixelExtent.lowerLeft.py);
yOffset = Math.max(0, Math.min(tileData.height - 1, yOffset));

Decoding the elevation

The next step is super-easy. Just use calculated indices and retrieve RGB values. Then decode it into the elevation.

const imageDataIndex = yOffset * (tileData.width * 4) + xOffset * 4;
const red = tileData.data[imageDataIndex];
const green = tileData.data[imageDataIndex + 1];
const blue = tileData.data[imageDataIndex + 2];

const elevation = -10000 + ((red * 256 * 256 + green * 256 + blue) * 0.1);

Calculating the profile

With all the code written above, we can just loop over samples and calculate elevation for each sample.

const samples = this.sampleProfileLine(coordinates);
const elevationProfile = [];
for (const c of samples) {
    const elevation = await this.elevationProvider.getElevation(c[1], c[0]);
    elevationProfile.push(elevation);
}

Rendering the chart

And finally, show the chart. In order to get a nice chart, we determine the minimum elevation in the profile to configure the chart component and along with some padding we set it as a low limit in the chart component.

const minElevation = Math.min(...elevationProfile);

this.chart = new Chartist.Line('#chart', {
    series: [elevationProfile]
}, {
    low: minElevation - 100,
    showArea: true,
    showPoint: false,
    fullWidth: true,
    lineSmooth: Chartist.Interpolation.cardinal({
        tension: 4,
        fillHoles: false
    })
});

Conclusion

With a little effort, it is quite easy to generate terrain profiles using MapTiler Terrain-RGB tileset. Access to the elevation data opens a wide range of applications such as microwave link planning, line-of-sight analysis, pipelines, power lines, etc.

Useful links

RGB-encoded elevation Data in MapTiler Cloud