Making a Beautiful Terrain Map

Introduction

Because my research focuses on the dynamics of civil war, my papers often include detailed maps showing the distribution of violence over space and time. There are easy ways to add base-layers from Google and other services to your maps, but often these have too many, or the wrong kind, of details that distract from the purpose of the map. This tutorial explains how to make your own terrain map that will look great in publications and presentations.

Download the data

The first step is to get access to digital elevation model (DEM) data. These are available at a variety of levels of granularity, so which is best depends on your purpose.

  • NASA Shuttle Radar Topographic Mission (SRTM) (1km / 15 arc seconds) here.
  • SRTM through CGIAR-CSI GeoPortal (90m / 3 arc seconds) here or here.
  • ASTER Global DEM (30m / 1 arc secomd) here or or here.

You should probably also be able to access these data using an API or FTP. This stackexchange post suggests the following code but I usually use local files kept on my external harddrive because they can take a bit to download.

curl -f -L 'http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_[01-72]_[01-24].zip' -o 'SRTM_V41_#1_#2.zip'

Aggregate the data

The geographic area you’re interested in will usually be split across several DEM tiles so we need to combine several seperate tile pieces into a single file. Here’s some code to search through a directory of DEM files and create a single mosaic out of them.

#load libraries
library(dplyr)
library(raster)
library(rgdal)
library(maptools)
library(beepr)

#locate the directory of the files and create a list
setwd("/Users/vincentbauer/Downloads/DEM")

#some projection settings
unproj <- CRS("+proj=longlat +datum=WGS84") #default WGS84 projection
proj <- CRS("+init=epsg:32638")  #projected to UTM 38N (Iraq)

#four different options for selecting files
files <- list.files(recursive=TRUE) #select all files in the directory
files <- list.files(recursive=TRUE, pattern="^I")  #select files in a particular row, such as those that start (^) with the letter I
files <- list.files(recursive=TRUE, pattern=".hgt$") #files that end ($) with the .hgt extension

#turn each file into a raster, saved as a list
files <- files[c(1,2)] #to speed up the tutorial, only use two DEM tiles
rasters.list <- sapply(files, raster)

#run the mosaic function on the raster list, which is a little complicated
#solution here, https://stackoverflow.com/questions/22109774/r-raster-mosaic-from-list-of-rasters
names(rasters.list) <- NULL
rasters.list$fun <- mean
mosaic <- do.call(mosaic, rasters.list)

#create a plot
plot(mosaic, col = gray.colors(20, start = 0, end = 1))

#notify you when the code completes
beep()

You should probably also project your mosaic file to the specific area you’re mapping

mosaic <- projectRaster(mosaic, crs=proj) #project the combined raster, can take a while

Create a hillshade

Once you have the DEM raster, you can compute the average elevation in different administrative regions and other related quantities, but it is just a visual representation of elevation, not a map with perspective. To make a map with perspective and show, we need to estimate what the shadows would look like if the sun was at a 45 degree angle to the plane.

slope <- terrain(mosaic, opt="slope", unit='radians')
aspect <- terrain(mosaic, opt="aspect", unit='radians')
hillshade <- hillShade(slope, aspect, angle=45, direction=315)

plot(hillshade, col = gray.colors(20, start = 0, end = 1))

For better presentation

That’s all you need to create a DEM raster that covers your area of interest. You could now compute the average elevation in different administrative regions or produce a simple map. But to make a map that really looks good, there are a couple of key steps.

Exaggerate the elevation

The world is actually pretty flat – mountains are not nearly as tall as the world is round – but this can create some really underwhelming maps with very little contrast. Almost all visual displays exaggerate the change in elevation versus distance on the ground. In the technical literature, this is called setting the z-factor (see articles here and here) but all you’re actually doing is multiplying the elevations by a set factor. Ten is usually pretty good for country-level maps but you might want to increase it more if you’re making a map of the entire globe.

#set the z-factor, which increases contrast
mosaic <- mosaic * 10 

#then create the hillshade
slope <- terrain(mosaic, opt="slope", unit='radians')
aspect <- terrain(mosaic, opt="aspect", unit='radians')
hillshade <- hillShade(slope, aspect, angle=45, direction=315)

#plot again
plot(hillshade, col = gray.colors(20, start = 0, end = 1))

Add some smoothing

Often the hillshade will be too detailed, where we only want to pick out major landfeatures. Use the following code to slightly smooth the data and make it better for presentations. Whether and how much smoothing to add depends on how large of an area you want to cover. In this example, the smoothing is not actually necessary because I’m only loading two of the many DEM tiles. If I was actually creating a mosaic of many tiles I would want to aggregate more.

hillshade2 <- aggregate(hillshade , fact = 5 , method = "bilinear" )
hillshade2 <- focal(hillshade2, w=matrix(1/9, nc=3, nr=3), mean)

plot(hillshade2, col = gray.colors(20, start = 0, end = 1))

Add some slopeshade

The hillshade by itself is a bit sharp looking. You can create something that looks a little bit more natural by also including the slope shade in the final image. We already made the slopeshade in the process of making the hillshade. You can tradeoff between the very sharp hillshade and the very smooth slope shade using the transparencies. I think that the best results come from 80% hillshade and 20% slope shade.

You also want to reverse the grey-scale for the slope shade because you want areas with more slope to be darker, not lighter. Notice that the greyscale is flipped below.

plot(slope, col = gray.colors(20, start = 0, end = 1))

plot(slope, col = gray.colors(20, start = 1, end = 0))

If you smoothed the hillshade, also smooth the slope otherwise they will not match up with each other.

slope2 <- aggregate(slope , fact = 5 , method = "bilinear" )
slope2 <- focal(slope2, w=matrix(1/9, nc=3, nr=3), mean)

Add some color

The DEM is just a representation of the elevation, so it can only express a single color scale. There are color-scales which are designed to show how colors often changed based on elevation (low levels are lush green, intermediate levels are brown or grey, and the highest levels are snow white) but these may not be accurate for the region you’re interested in (they’ll do a much better job with Europe and North America). You can add more realistic color by downloading data on land coverage and overlaying this with the DEM.

Here, I have downloaded a large landcover raster from Natural Earth that covers the entire globe.

landcover <- "/Users/vincentbauer/Downloads/NE1_LR_LC.tif"
landcover <- brick(landcover)  #brick for multi-band rasters
landcover <- crop(landcover, extent(mosaic))  #crop it to the area of interest

#view the color image
plotRGB(landcover,r=1,g=2,b=3) 

You can also download the landcover raster through the Natural Earth API but I do not know the steps to address the multiple bands.

library(rnaturalearth)

#low resolution version, replace LR with HR for high resolution
landcover <- ne_download(scale = 10, type = 'NE1_LR_LC', category="raster")  

Combine these elements

For the final map, we want to combine all of these elements: the landcover, the hillshade, and the slope shade.

plotRGB(landcover,r=1,g=2,b=3, alpha=255)
plot(hillshade, col=grey.colors(100, start=0, end=1), legend=F, alpha=.8, add=TRUE)
plot(slope, col=grey.colors(100, start=1, end=0), legend=F, alpha=.2, add=TRUE)

Here’s the ggplot version

library(ggplot)
library(ggsn) #scale bars

#ggplot makes me turn the raster into points
hills.df <- rasterToPoints(hillshade)
hills.df <-  data.frame(hills.df)
colnames(hills.df) <- c("lon", "lat", "hills")

#the key to a pretty map is merging the slope shade with the hillshade
slope.df <- rasterToPoints(slope)
slope.df <-  data.frame(slope.df)
colnames(slope.df) <- c("lon", "lat", "slope")
slope.df$slope <- 1- slope.df$slope #invert the scale so that more slope is darker

plot <- ggplot(aes(x=long, y=lat, group=group)) +
    coord_quickmap(xlim = c(extent@xmin, extent@xmax),ylim = c(extent@ymin, extent@ymax), expand=FALSE) +
    geom_raster(data = hills.df, aes(lon, lat, fill = hills, group=1), alpha = .6, interpolate=TRUE) +
    geom_raster(data = slope.df, aes(lon, lat, fill = slope, group=2), alpha = .2, interpolate=TRUE) +
    scale_fill_gradientn(colours = grey.colors(100, start = 0, end=.95), guide='none') +
    theme_map() +
    north(symbol = 12, scale = .05,  location = 'topright',anchor=c(x=70.48, y = 33.23),
          x.min =extent@xmin, x.max=extent@xmax, y.min=extent@ymin, y.max=extent@ymax) +
    ggsn::scalebar(dist=10, height=.015, st.size=2, dd2km=TRUE, model="WGS84", location = 'bottomleft', anchor=c(x=69.55, y = 32.77),
                   x.min =extent@xmin, x.max=extent@xmax, y.min=extent@ymin, y.max=extent@ymax)

ggsave(file="../Charts/example_map.png", plot = plot, width = 5, height=3, units="in")