Finding max elevation of a DEM tile

There are several sources of Digital Elevation Model (DEM) tiles, sometimes called “terrain tiles” online. Instead of containing cartographic features such as roads, rivers and cities, each pixel of a DEM tile encodes the elevation found in the corresponding location. Recently, I became interested in finding the maximum elevation for each one of these tiles in order to accelerate shadow calculation on shademap.app

DEM tile for zoom level 0

At zoom level 0, the entire world map is encoded on just one tile. I would expect the max elevation of any pixel of this tile to equal the height of Mount Everest (8848m), since this is the highest point on earth. In reality, I found that the maximum value was actually (5887m). My theory is that there is some averaging or interpolation that is done when reducing the data to a single 256x256 pixel tile. Each pixel of a tile at zoom level 0 technically represents ~156 square kilometers, so a large amount of other elevation values would be used to calculate the pixel where Mount Everest sits.

Getting the tiles #

I used the Terrain Tiles data set hosted by AWS Open Data because it covers the entire planet. I used the following command to download tiles for a given [zoom] level into a tiles/ directory on my machine:

aws s3 sync --no-sign-request s3://elevation-tiles-prod/terrarium/[zoom] tiles/[zoom] --exclude "*" --include "*.png"

I believe I installed the aws command line utility using brew install awscli. You will need a ~/.aws/credentials and ~/.aws/config file for the CLI to work. There are other formats beside .png available if you adjust the --exclude/include flag. Below the amount of disk space required to download each corresponding zoom level:

Calculating max elevation per tile #

Each pixel of a .png tile uses the red, green and blue channels to encode the elevation. One byte for red, one byte for green and one byte for blue, each containing a value from 0-255. The following formula converts the color values into meters above sea level:

meters_above_sea_level  = (red * 256 + green + blue / 256) - 32768

The GDAL library (installed with brew install gdal) contains a utility for displaying metadata for DEM tiles, including the maximum red, green and blue encountered in the image. A single pixel might not contain the maximum red, blue and green value, however, adding the maximum red, blue and green values for a tile yields an upper bound for the maximum elevation.

$> gdalinfo 0/0/0.png -stats

Band 1 Block=256x1 Type=Byte, ColorInterp=Red
  Min=119.000 Max=135.000 
  Minimum=119.000, Maximum=135.000, Mean=127.557, StdDev=2.453

Band 2 Block=256x1 Type=Byte, ColorInterp=Green
  Min=0.000 Max=255.000 
  Minimum=0.000, Maximum=255.000, Mean=126.723, StdDev=73.106

Band 3 Block=256x1 Type=Byte, ColorInterp=Blue
  Min=0.000 Max=0.000 
  Minimum=0.000, Maximum=0.000, Mean=0.000, StdDev=0.000

Max elevation of each tile for given zoom level #

I wrote a bash script. The only argument is the zoom level to do max elevation calculations for. Here is the output produced when running on the four tiles included in zoom level 1:

$> ./max-height.sh 1
4095 1/0/0.png
5375 1/0/1.png
6143 1/1/0.png
4095 1/1/1.png

Here is the script:

#!/bin/bash

let ZOOM=$1
let FILE_NUM=2**$ZOOM-1

for FOLDER in $(seq 0 $FILE_NUM)
do
  for FILE in $(seq 0 $FILE_NUM)
  do
    # print Maximum=1234.5 value for R,G,B channels
    gdalinfo $ZOOM/$FOLDER/$FILE.png -stats | grep Maximum | \
    # put all values on a single line seperated by whitespace 12.3 45.6 78.9
    cut -d "," -f 2 | cut -d "=" -f 2 | tr '\n' ' ' | \
    # calculate and print elevation in meters along with filename
    awk -v zoom="$ZOOM" -v folder="$FOLDER" -v file="$FILE" '{print ($1 * 256.0 + $2 + $3 / 256.0) - 32768 " " zoom "/" folder "/" file ".png"}'
  done
done

The script runs in under a minute for zoom levels < 5. After that it is significantly slower.

Conclusion #

It appears that there is some errors in the DEM data since values over 10,000 meters exist on some tiles