Making a Topographical Map of Alaska

Objective

To create a topographical (shaded colour relief map) image of Alaska (not publication quality).

Assumptions, Constraints and Considerations

  1. Open source tools.
    In this article I primarily use the GDAL (Geospatial Data Abstraction Library) tools from the Open Source Geospatial Foundation (OSGeo) project.
  2. Open data.
    Data freely available to the general public, without onerous terms and conditions.
  3. Command line.
    I generally tend to use command line tools rather than graphical tools, partly because of a smaller footprint, faster processing and general convenience, but also because of the ability to script and automate the workflow. In this article I use a Bourne-shell compatible shell on a POSIX (Unix, Linux) system and the reader is assumed to be familiar with shell operation and programming.
  4. File formats.
    For raster geospatial data, I'm using the GeoTIFF format, since this is what most of the tools use by default. GeoTIFF is an extension of the TIFF format that adds attributes to define the projection, extent of the file, interval, coordinate system etc. For vector geospatial data I tend to use the ESRI Shapefile format for processing. This stores geometries of a single type (linestring, point, polygon) actually as a set of files with the same base name but different extensions, and the set of files is typically specified by the filename with the .shp extension.

Issues

Immediate problems:

  1. Alaska straddles the anti-meridian (180°) which some tools can struggle to handle.
  2. DEM coverages. One source of DEM data, SRTM [1], is high resolution but only covers land masses up to 60°N. The GEBCO database [2] covers the whole earth (including bathymetry) but is lower resolution. I wanted the highest resolution at any part of the map, so these had to be combined.

Steps

The data processing steps are roughly as follows:

  1. Data preparation.
  2. Create DEM files.
  3. Reproject.
  4. Create hillshade files.
  5. Create hyposometric tints (relief shade).
  6. Combine the relief shade with the hillshade to produce an image.

Filename Conventions

With geospatial processing and any other processing that involves multiple steps with a number of intermediate files, it's important to impose a consistent naming convention.

Initially, I started with a roughly rectangular bounding area containing all of Alaska. Because of issues to do with different DEM data sources and the anti-meridian, I divided the area into four parts, and used as a naming convention the prefixes alaska_ne, alaska_nw, alaska_se and alaska_sw, shown in the table below. Note that with this nomenclature, the "west" files are actually in the eastern hemisphere!

Filename convention: prefixes and corresponding extents
Prefix DEM Extent
alaska_ne GEBCO 60°N 180°W to 71°30’N 129°30’W
alaska_nw GEBCO 60°N 172°E to 71°N 180°E
alaska_se SRTM 51°N 180°W to 60°N 129°30’W
alaska_sw SRTM 51°N 172°E to 60°N 180°E

(Notice that the alaska_nw part does not actually contain Alaska but the landmass across the Bering Strait, which is in the Russian Federation. This will be dealt with later.)

Extending the naming convention to denote the purposes/data of each file, I designated the suffixes shown in the following table. "Projected" coordinates means latitude/longitude coordinates that have been projected onto a Cartesian plane by some transformation, and the units are metres.

Filename suffixes and file types
Suffix Type Spatial Coordinates Values
None DEM Latitude/longitude Height (ft)
proj DEM Projected (m) Height (ft)
hillshade Hillshade Projected (m) Hillshade [0, 255]
relief Relief tints Projected (m) Colour (rgb or rgba)

1. Data preparation

The data for generating the shaded relief map are Digital Elevation Model (DEM) files obtained from two sources: SRTM (Shuttle Radar Topography Mission) [1] and GEBCO [2] databases. SRTM data are sampled at 30m intervals and cover land masses but not the polar regions (about 60N). GEBCO (General Bathymetric Chart of the Ocean) provides bathymetric data as well as land surface elevations for the whole globe at 15 arc-second intervals.

SRTM can be downloaded as files in HGT format covering 1x1 degrees, so you will need to download those sufficient to cover the planned area and put them in a directory. I converted the SRTM hgt files into GeoTIFF format for easier general handling, with the following command (the dollar sign is the shell command line prompt and not part of the command).

$ gdalwarp -srcnodata 32767 -dstalpha <hgtfile> <tiff_file>

This replaces "no data" values (32767) that cover ocean areas and voids in the data set with an "alpha mask". In the discussion below, I put the tiff files in a directory labelled <srtm> in the text.

The complete GEBCO data can be downloaded in GeoTIFF (4 Gbytes compressed!) and NetCDF formats. Otherwise, a web interface is available to download user-specified areas. In the discussion below, I merged all the GEBCO data from the 2019 dataset into a single file GEBCO_2019.tiff.

For cutting to country borders, one must provide a polygon (or polygons) that define the country's land outlines. I used the Natural Earth dataset [3] at the highest available resolution.

2. Create DEM files of the required area

Alaska has three bits to consider: the panhandle, the main body of Alaska on the continental landmass, and the Aleutian Islands.

Because the SRTM DEM only goes up to 60°N, but I want higher resolutions where available, we're going to have to use the SRTM for south of 60°N and GEBCO for north of 60°N. In addition, to avoid some of the issues with the anti-meridian, we are going to need to extract the western and eastern hemispheres separately.

I therefore initially create four DEM files that will be merged later. Using the command line tool gdal_merge.py:

$ gdal_merge.py -v -o alaska_se.tif -ul_lr -180 60 -129.5 51 <srtm>/N5?W1??.tiff
$ gdal_merge.py -v -o alaska_sw.tif -ul_lr 172 60 180 51 <srtm>/N5?E1??.tiff
$ gdal_merge.py -v -o alaska_ne.tif ul_lr -180 71.5 -129.5 60 GEBCO_2019.tiff
$ gdal_merge.py -v -o alaska_nw.tif ul_lr 172 71.5 180 60 GEBCO_2019.tiff

where:

-v
Species the 'verbose' option. (Leave it out if you don't want diagnostic information printed, or if you're making a script.)
-o <filename>
Specifies the output filename.
-ul_lr <lx> <uy> <rx> <by>
Specifies the bounds to extract (short for 'upper-left lower-right').

The last argument(s) are a list of input files from which to extract the required data. There's a "gotcha" here when using wildcard expansion like <srtm>/*.tiff — if it includes too many files, gdal_merge might take a long time to scan a lot of files it doesn't use, and worse, the expanded wildcard might overflow the command line length limit and the processing may not include all the files! The examples above attempt to use a more restricted set of files.

(If you're scripting, it's advisable to define the limits as variables, so as if you need to fine-tune the limits you need only make a change in one place in the file.)

3. Reproject

The DEM files we have use latitude and longitude for the horizontal plane and feet for the vertical plane, so some scaling needs to be done when we compute the hillshade to match the horizontal and vertical scales. Unfortunately, the scaling of latitude and longitude are not uniform, particularly this far from the equator. It would be therefore advisable to reproject the DEM files onto a "flat" surface with uniform horizontal scaling.

Cartographers have come up with a number of defined "standard" projections of different types with characteristics to meet the needs of map users. Some of these are specific to certain parts of the world to meet requirements for accuracy within a given area. However, all projections of a spheriod onto a flat plane necessarily involve a distortion, and the choice of projection is not necessarily straightforward [4]. A quick search of projections for Alaska threw up NAD83/Alaska Albers (EPSG:3338) as a candidate [5]. This is an Albers equal-area conic projection that covers our area and projects to metres.

Three projections are illustrated below. "WGS84" just plots latitude and longitude values on a rectilinear grid. "Webmercator" is widely used by web map apps.

Projections
WGS84 Webmercator NAD 83/Alaska Albers
proj_wgs84 proj_webmercator proj_epsg3338

The basic command to reproject is:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3338 -r cubicspline -dstalpha <infile> <outfile>

where:

-s_srs EPSG:4326
Specifies the spatial reference set of the input (source) file. In this case, the DEM files use the WGS84 system, which is referenced by its EPSG (European Petroleum Survey Group) reference code of EPSG:4326.
-t_srs EPSG:3338
Specifies the target spatial reference set of the output file. In this case, NAD 83/Alaska Albers (EPSG:3338).
-r cubicspline
Specifes the resampling method. The default (if not specified) is quick and dirty (worst interpolation quality). There are several to choose from but I picked Cubic Spline interpolation.
-dstalpha
Tells gdalwarp to create an alpha band to identify no-data pixels. This is important for later processing when we want to merge several files, since the rectangular output raster will contain areas that are not present in the input.

You can either type in the command repeatedly for each file but I used the following fragment of shell script:

for FILE in alaska_??.tif; do
  OUTFILE=`echo ${FILE} | sed -e 's/\.tif/_proj\.tif/'`
  gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3338 -r cubicspline -dstalpha ${FILE} ${OUTFILE}
done

For example, this takes our alaska_ne.tif DEM file in WGS84 coordinates (latitude, longitude) as input and outputs the reprojected file in Cartesian coordinates (metres) as alaska_ne_proj.tif.

4. Create hillshade files

Using the gdaldem command. The basic command is:

gdaldem hillshade -z 4 -s 3.28084 <input_file> <output_file>

where:

hillshade
Specifies that gdaldem computes a hillshade file
-z 4
Specifies the vertical exaggeration. I used 4 in this case.
-s 3.28084
Specifies the ratio of horizontal to vertical scales; in this case, the ratio between metres (horizontal) and feet (vertical)

Again, you can run the commands individually or use a shell command

for FILE in alaska_??_proj.tif; do
  OUTFILE=`echo $FILE|sed -e 's/proj/hillshade/'`
  gdaldem hillshade -z 4 -s 3.28084 $FILE $OUTFILE
done

5. Create Hypsometric tints

We now want to create colour relief tints. Again, we can use the gdaldem command but we need an additional file that maps elevation values to colours: a "colour ramp". This is specified in a simple plain-text file with four or five columns per line in the order:

<elevation> <red> <green> <blue> [<alpha>]

where:

<elevation>
Is the elevation value (integer or floating point value)
<red> <green> <blue>
are three integer values in the range [0, 255] that specify the colour corresponding to the elevation value
<alpha>
is an optional alpha value ranging from 0 (transparent) to 255 (full opacity)

I created a plain text file called colour_ramp.txt for this purpose, with the following contents:

0    255 255 255 0
1      6 100  51
50    15 122  47
500  205 202 112
1200 169  80  12
1700 135  26  24
2800 111 105 105
4000 246 246 246
5000 255 255 255

(Excess spaces are not significant.) Note that this will map elevation values of 0 or less (mostly sea) to transparent. However, this is a pretty crude way of extracting coastlines and won't of course produce good results if there are any land elevations below sea level.

Then to create a relief map from a DEM file with the colour ramp:

gdaldem color-relief -alpha <input_dem_file> <colour_ramp_file> <output_file>

or to batch process:

for FILE in alaska_??_proj.tif; do OUTFILE=`echo $FILE|sed -e 's/proj/relief/'`; gdaldem color-relief -alpha $FILE colour_ramp.txt $OUTFILE; done

6. Combine

The shaded relief files may now be combined with the hillshade file to produce output images. Hillshade files are greyscale files with pixel values in the range [0, 255]. The pixel values of the coloured relief file are "blended" (combined mathematically) pixel by pixel with the corresponding pixel values in the hillshade file. Parameters include the blend mode (mathematical function that combines the relief colour and hillshade value) and functions that can modify the hillshade layer, e.g. adjusting its brightness and contrast, to adjust the blended effect. The figures below show renderings with hypsometric tints only and with a hillshade overlay, at x4 vertical exaggeration.

Hypsometric tints only With Hillshade overlay
Hypsometric tints Hypsometric tints and hillshade overlay

Intermediate Result

I loaded the files using the QGIS tool, with the hillshade layer blend mode set to Multiply and brightness set to 80, to get the following result:

Alaska relief map with hillshade

At fairly low resolution it's broadly acceptable. There are some artefacts in the land and ocean at the seams of the various files. Artefacts in the ocean can be edited out with a pixel editor. (Save as a JPEG with a sidecar file specifying the projection, and edit the image.)

Also, although the ocean floor is not particularly steep in the Arctic Ocean, there exists some hillshading from the bathymetry data in the sea parts north of 60N.

Closer inspection reveals difference between the GEBCO (15 arc second resolution) and SRTM models, and the discontinuous border between them.

Image extract showing different resolutions of GEBCO and SRTM data

Cutting to the State Border

In some cases, one may wish to restrict the relief map to a certain geographical area; for example, only the state of Alaska. This can be achieved by defining a polygon and using it to mask the desired areas, for example when creating the projected DEM file with gdalwarp.

I downloaded land polygon data at the highest resolution (1:10m) from the Natural Earth website. The data is available in ESRI Shapefile, geoDB or GeoTIFF formats. I assume below a set of ESRI Shapefiles ne_10m_admin_countries.shp.

One must then extract the polygons for the USA from the datafile to use as a "cutline" in the gdalwarp operation. This could be carried out directly in gdalwarp, for example:

$ gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3338 -r cubicspline -dstalpha \
  -cutline ne_10m_admin_countries.shp -crop_to_cutline \
  -cwhere "ADMIN='United States of America'" \
  <infile> <outfile>

However, in this case the polygons straddle the anti-meridian and cause an extremely large output file because the extreme western coordinates are expressed with negative longitude values less than -180 degrees. I found it better to extract the individual Alaska polygon without the portion of the Aleutian islands in the eastern hemisphere and to use a simple bounding box for the Aleutian islands in the eastern hemisphere. This can be done with the ogr2ogr tool:

$ ogr2ogr \
  -sql "select * from ne_10m_admin_0_countries WHERE name='United States of America" \
  -f "ESRI Shapefile" \
  -clipsrc -180 51 -129.5 71 \
  alaska_east_land.shp ne_10m_admin_countries.shp

where:

-sql <statement>
is an SQL select statement that selects the desired data from the input file
-clipsrc <xmin> <ymin> <xmax> <ymax>
Specifies the "bounding box" extents to clip the file to; in this case, the extents of Alaska in the western hemisphere

This produces a polygon in set of Shapefiles alaska_east_land.shp, as shown below (projected to EPSG:3338).

Alaska land polygon, except Aleutian islands in the eastern hemisphere

Then using them at the stage of creating the projected DEM files for the western hemisphere (parts of Alaska east of 180W):

for FILE in alaska_?e.tif; do
  OUTFILE=`echo $FILE | sed -e 's/\.tif/_proj\.tif/'`
  gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3338 \
           -r cubicspline -dstalpha \
           -cutline alaska_east_land.shp -crop_to_cutline $FILE $OUTFILE
done

For the part of the Aleutian islands in the eastern hemisphere, a simple bounding box from 51N 172E to 55.1N 180E will suffice, and the alaska_nw portion is redundant (it is in Russia rather than Alaska). We can specify the bounding area directly to gdalwarp using the -te and -te_srs options:

$ gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3338 \
  -te 172 51 180 55.1 -te_srs EPSG:4326 \
  -r cubicspline -dstalpha \
  alaska_sw.tif alaska_sw_proj.tif

Hillshade and colour relief files can then be created in the same way as above. The final result with just the state of Alaska is like this:

State of Alaska relief map with hillshade

References

[1](1, 2) Shuttle Radar Topography Mission: https://www2.jpl.nasa.gov/srtm
[2](1, 2) General Bathymetric Chart of the Oceans: https://www.gebco.net
[3]Natural Earth: https://www.naturalearthdata.com
[4]Corey, M.: Choosing the Right Map Projection: https://source.opennews.org/articles/choosing-right-map-projection/
[5]EPSG:3338 (NAD83/Alaska Albers): https://epsg.io/3338/

social