Digital Elevation Models (DEMs) are maps of the bare topographic surface of the Earth, excluding features like trees and buildings etc, and usually collected by satellites. DEMs can be used to map topographic elevation, relief, slope, and aspect, and are useful sources for a range of purposes such as flood risk modelling, contour mapping, correcting aerial and satellite imagery, and land cover classification.
A well known source of this data is the USGS/NASA Shuttle Radar Topography Mission (SRTM), which provides tiled DEM data for the entire globe within latitudes 60° N and 56° S, at a resolution of 1 arc-second (~30m).
In this post, I describe some methods to source and process SRTM imagery. Before you can download the data though, you’ll first need to register for a free Earthdata account.
PEC Sample DEM
As an example, we can first create an elevation map for a smaller region, in this case Prince Edward County, Ontario (the peninsular in north eastern lake Ontario). For downloading a handful of tiles, this useful tool created by Derek Watkins can be used to select tiles from an interactive map: http://dwtkns.com/srtm30m/.
Download and unzip the four tiles that cover Prince Edward County, N44W078 (pictured below) and the three others east and south of it: N44W077; N43W077; and N43W078, entering your Earthdata username and password when prompted. Also, download Canadian census divisions data from Statistics Canada (Cartographic Boundary File), which we can use to clip the SRTM data.
With the four tiles downloaded and unzipped to .hgt files, we can merge them using GDAL, building a VRT file, a virtual dataset useful for intermediate processing steps. The below example uses Docker and the osgeo/gdal image, but can be modified for other configurations.
docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdalbuildvrt /mnt/data/pec-srtm.vrt /mnt/data/N44W077.hgt /mnt/data/N44W078.hgt /mnt/data/N43W077.hgt /mnt/data/N43W078.hgt
We can then clip the VRT file to the census divisions layer with gdalwarp as below, selecting the ‘Prince Edward’ feature as a cutline. This creates a clipped DEM GeoTiff like that displayed at the head of this post.
docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdalwarp -cutline /mnt/data/lcd_000b16a_e/lcd_000b16a_e.shp -cwhere "CDNAME = 'Prince Edward'" -crop_to_cutline /mnt/data/pec-srtm.vrt /mnt/data/pec-dem.tif
Some alternative calculations and visualizations can be created with the gdaldem tool as below, for example hillshade (shaded relief to represent 3d terrain), slope, and aspect (the direction of the slope).
docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdaldem hillshade -s 111120 /mnt/data/pec-dem.tif /mnt/data/pec-hillshade.tif docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdaldem slope -s 111120 /mnt/data/pec-dem.tif /mnt/data/pec-slope.tif docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdaldem aspect /mnt/data/pec-dem.tif /mnt/data/pec-aspect.tif
Entire Ontario DEM
Needless to say, it would be tedious and time consuming to manually download individual tiles to cover the entirety of Ontario (or beyond) in this way. So to scale up to a larger area we can write a script to automate the process using e.g. Python, like shown in the notebook below (available for download here).
This script uses a layer of the SRTM grid extent and a Provinces/territories boundaries layer, again from Statistics Canada, overlaying the two to determine which tiles to download. Note that the entire selection consists of 170 files, totalling ~750mb, which make take some time to download.
With the tiles downloaded we can again use GDAL to merge and crop the tiles in a similar way to the PEC sample e.g.
docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest sh -c 'gdalbuildvrt /mnt/data/on-srtm.vrt /mnt/data/*.hgt' docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdalwarp -cutline /mnt/data/lpr_000b16a_e/lpr_000b16a_e.shp -cwhere "PRNAME = 'Ontario'" -crop_to_cutline -co COMPRESS=LZW -co BIGTIFF=YES -co TILED=YES /mnt/data/on-srtm.vrt /mnt/data/on-dem.tif
Addendum: Web Map
I you want to use the results for a web map, web tiles can be created using gdal2tiles.py, after first converting the data type via gdal_translate. The below commands generate web map tiles at zoom levels 4-10, along with the HTML content for rendering a web map, such as leaflet.
docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdal_translate -of VRT -ot Byte -scale /mnt/data/on-dem.tif /mnt/data/temp.vrt docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdal2tiles.py -e -x -t "Ontario DEM" --zoom 4-10 --processes 2 /mnt/data/temp.vrt /mnt/data/web-tiles