Terrain

Generate contours, hillshade, Terrain RGB, and slope angle shading map tiles from Digital Elevation Models (DEMs).

Overview

I use OpenMapTiles to create self-hosted vector map tiles. However, I'm interested in building a topographic outdoors-oriented map. Contours, hillshading, and slope-angle shading are instrumental in intuitively understanding the terrain.

Source data

Since my map is focused on the continental United States, I use data from the US Geological Survey (USGS), which has high accuracy but is limited to the US. If you're interested in creating a map with international scope, check out 30-meter SRTM data, which is generally the most accurate worldwide source available. Most of the code in this repository is applicable to other DEM sources with few modifications.

Regarding the USGS data, they have a few sources available:

• 1 arc-second seamless DEM. This has ~30m horizontal accuracy, which is accurate enough for many purposes, and gives it the smallest file sizes, making it easy to work with.
• 1/3 arc-second seamless DEM. This dataset has the best precision available (~10m horizontal accuracy) for a seamless dataset. Note that the file sizes are about 9x bigger than the 1 arc-second data, making each 1x1 degree cell about 450MB unzipped.
• 1/9 arc-second project-based DEM and 1-meter project-based DEM. These have very high horizontal accuracy, but aren't available for the entire US yet. If you want to use these datasets, go to the National Map download page, check "Elevation Products (3DEP)", and then click "Show Availability" under the layer you're interested in, so that you can see if they exist for the area you're interested in.

For my purposes, I use the 1 arc-second data initially for testing, but then the 1/3 arc-second data for production use. In a couple years, once more of the United States is mapped at 1 meter resolution, I'd suggest considering using 1 meter data where available.

Terrain RGB vs Raster Hillshade

Historically, the way to make a hillshade is to generate raster images where each cell stores the level of grayscale to display.

With Mapbox GL, you have a new option: Terrain RGB tiles. Instead of encoding the grayscale in the raster, it encodes the raw elevation value in 0.1m increments. This enables a whole host of cool things to do client-side, like retrieving elevation for a point, or generating the viewshed from a point.

I currently exclusively use Terrain RGB tiles in my projects because it allows for the possibility of doing cool things in the future if I have time.

Terrarium dataset

If you want to use Terrain RGB tiles, but don't want to create them yourself, it's also currently possible to use the publicly-hosted Terrarium dataset on AWS Public Datasets for free. It isn't even in a requester-pays bucket. I don't know how long this will be available for free, so I figured I'd just generate my own.

If you want to go that route, set this as your source in your style.json (note "encoding": "terrarium"):

"terrarium": {
  "type": "raster-dem",
  "tiles": [
  	"https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png"
  ],
  "minzoom": 0,
  "maxzoom": 15,
  "tileSize": 256,
  "encoding": "terrarium"
}

Note that the Terrarium dataset uses a different encoding than Mapbox's RGB tiles:

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

Contours

I originally generated contours vector tiles straight from USGS vector data (https://github.com/nst-guide/contours). These downloads have contour lines pregenerated in vector format, and so are quite easy to work with; just download, run ogr2ogr to convert to GeoJSON, and then run tippecanoe to convert to vector tiles.

There are a couple drawbacks of using the USGS vector contour data:

1. Inconsistent vertical spacing. In some areas, I found that the data included 10' contours, while others had a precision of only 80'. Since I desired a 40' contour, this meant that there were occasionally visual discontinuities in the map between tiles with source data of at least 40' precision and source data of less than 40' precision.

2. Metric contours. Although I'm currently making maps of the United States, where imperial measurements are the standard, it's desirable to provide metric measurements as well. With Mapbox GL, it's quite easy to write a style expression that converts feet into meters, but the spacing of the contour lines will still be in feet. If neighboring imperial contour lines are 2000 feet and 2040 feet, displaying those same lines on a metric map would represent unintuitive values of 609.6 meters and 621.8 meters.

   In order to display metric contour lines, it is necessary to generate a separate set of contours with lines that represent 100 meters, 110 meters, etc. This must be generated from the original DEMs.

Slope-angle shading

Slope angle shading displays polygons representing the slope of the terrain in degrees. This is just a couple of GDAL functions joined with a color scheme that roughly matches the one from Caltopo.

Integration with style.json

The style JSON spec tells Mapbox GL how to style your map. Add the hillshade tiles as a source to overlay them with the other sources.

Within sources, each object key defines the name by which the later parts of style.json should refer to the layer. Note the difference between a normal raster layer and the terrain RGB layer.

"sources": {
  "contours": {
    "type": "vector",
    "url": "https://example.com/contours/tile.json"
  },
  "slope_angle": {
    "type": "raster",
    "url": "https://example.com/slope_angle/tile.json",
  	"tileSize": 512
  },
  "terrain-rgb": {
    "type": "raster-dem",
    "url": "https://example.com/terrain_rgb/tile.json",
    "tileSize": 512,
    "minzoom": 0,
    "maxzoom": 12,
    "encoding": "mapbox"
  }
}

Where the tile.json for a raster layer should be something like:

{
    "attribution": "<a href=\"https://www.usgs.gov/\" target=\"_blank\">© USGS</a>",
    "description": "Hillshade generated from 1 arc-second USGS DEM",
    "format": "png",
    "id": "hillshade",
    "maxzoom": 16,
    "minzoom": 0,
    "name": "hillshade",
    "scheme": "tms",
    "tiles": ["https://example.com/url/to/tiles/{z}/{x}/{y}.png"],
    "version": "2.2.0"
}

Later in the style JSON, refer to the hillshade to style it. Example for terrain RGB:

{
  "id": "terrain-rgb",
  "source": "terrain-rgb",
  "type": "hillshade",
  "paint": {
    "hillshade-shadow-color": "hsl(39, 21%, 33%)",
    "hillshade-illumination-direction": 315,
    "hillshade-exaggeration": 0.8
  }
}

Installation

Clone the repository:

git clone https://github.com/nst-guide/terrain
cd terrain

Then install dependencies

conda env create -f environment.yml
source activate terrain

Code Overview

download.py

Downloads USGS elevation data for a given bounding box.

> python download.py --help
Usage: download.py [OPTIONS]

Options:
  --bbox TEXT  Bounding box to download data for. Should be west, south, east,
               north.  [required]
  --overwrite  Re-download and overwrite existing files.
  --high_res   Download high-res 1/3 arc-second DEM.
  --help       Show this message and exit.

This script calls the National Map API and finds all the 1x1 degree elevation products that intersect the given bounding box. Right now, this uses 1 arc-second data, which has about a 30 meter resolution. It would also be possible to use the 1/3 arc-second seamless data, which is the best seamless resolution available for the continental US, but those file sizes are 9x bigger, so for now I'm just going to generate from the 1 arc-second.

The script then downloads each of these files to data/raw/. By default, it doesn't re-download and overwrite a file that already exists. If you wish to overwrite an existing file, use --overwrite.

unzip.sh

Takes downloaded DEM data from data/raw/ or data/raw_hr/, unzips it, and places it in data/unzipped/ or data/unzipped_hr/.

Usage

Preparation

Download desired DEM tiles, then unzip them, build a VRT (Virtual Dataset), and optionally download my fork of gdal2tiles which allows for creating 512x512 pngs.

If you're not using data from the USGS, you'll need to figure out which DEMs to download yourself. The files should be geospatially adjacent so that there are no holes in the generated map data.

# Download for Washington state
python download.py --bbox="-126.7423, 45.54326, -116.9145, 49.00708"
# Or, download high-resolution 1/3 arc-second tiles
python download.py --bbox="-126.7423, 45.54326, -116.9145, 49.00708" --high_res
bash unzip.sh
# Create seamless DEM:
gdalbuildvrt data/dem.vrt data/unzipped/*.img
gdalbuildvrt data/dem_hr.vrt data/unzipped_hr/*.img
# Download my fork of gdal2tiles.py
# I use my own gdal2tiles.py fork for retina 2x 512x512 tiles
git clone https://github.com/nst-guide/gdal2tiles
cp gdal2tiles/gdal2tiles.py ./

Terrain RGB

# Create a new VRT specifically for the terrain RGB tiles, manually setting the
# nodata value to be -9999
gdalbuildvrt \
    `# Set the VRT's nodata value to -9999` \
    -vrtnodata -9999 \
    data/dem_hr_9999.vrt data/unzipped_hr/*.img
# Reproject DEM to web mercator
gdalwarp \
    `# Use cubicspline averaging` \
    -r cubicspline \
    `# Source projection is EPSG 4269