DEM Tutorial

Matthias Schultheiss

Introduction

This short tutorial will show you how to create a hillshaded topographic map in Rstudio. We will be using Shuttle Radar Topography Mission (SRTM) data, a near global Digital Elevation Model (DEM) collected in February 2000 aboard NASA’s Space Shuttle Endeavour (mission STS-99, link: https://dwtkns.com/srtm/). It was globally released by late 2015. The mission used a X-Band mapping radar to measure the Earth’s topography, built in collaboration with the U.S. Jet Propulsion Laboratory, the U.S. National Imagery and Mapping Agency (now the National Geospatial-Intelligence Agency), and the German and Italian space agencies. For more Information click here: https://www2.jpl.nasa.gov/srtm/.

Setup

To get started, the used packages and paths must be defined and loaded. It is recommended to use as few packages as possible to ensure the reproducibility of the code in the future and under a different R version. For this reason, only packages that are regularly updated and used by the majority of the R community will be used here.

# Loading packages
if (!require(sp)){install.packages("sp"); library(sp)}
if (!require(sf)){install.packages("sf"); library(sf)}
if (!require(raster)){install.packages("raster"); library(raster)}
if (!require(dplyr)){install.packages("dplyr"); library(dplyr)}
if (!require(tidyverse)){install.packages("tidyverse"); library(tidyverse)}

Creating the map

Since the goal is to create a DEM with a Hillshade for a very specific region (Baden-Württemberg, Germany), the boundaries of the space to be considered are needed in addition to the raster data already mentioned in the introduction. While there are various open source providers (e.g. OSM), this project decided to use the publicly available data from the government of the respective federal state for convenience reasons. (Source: https://www.lgl-bw.de/lgl-internet/opencms/de/07_Produkte_und_Dienstleistungen/Open_Data_Initiative/)

# Import vector data
bw.landesgrenzen <- st_transform(
  st_read("data/rohdaten/NUTS/AX_Gebiet_Bundesland.shp"),
  4326)
## Reading layer `AX_Gebiet_Bundesland' from data source `V:\github\nils_website\blogdown\content\showcase\schultheiss\data\rohdaten\NUTS\AX_Gebiet_Bundesland.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 388337.8 ymin: 5265159 xmax: 610069.4 ymax: 5515632
## epsg (SRID):    NA
## proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
# Import raster data
tif1<-raster("data/rohdaten/SRTM/srtm_38_03.tif")
tif2<-raster("data/rohdaten/SRTM/srtm_39_03.tif")

# Merge both raster layer
bw.mosaic<-mosaic(tif1,
                  tif2,
                  fun=mean)
# Check properties (CRS) and set MinMax 
setMinMax(bw.mosaic)
## class      : RasterLayer 
## dimensions : 6001, 12001, 72018001  (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333  (x, y)
## extent     : 4.999584, 15.00042, 44.99958, 50.00042  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : -32, 4783  (min, max)
# Crop raster to boundaries of Baden-Württemberg
bw.crop <- crop(bw.mosaic, bw.landesgrenzen)
bw.SRTM <- mask(bw.crop, bw.landesgrenzen)

# Query lowest and highest point in BaWü
cellStats(bw.SRTM, min) 
## [1] 26
cellStats(bw.SRTM, max)
## [1] 1491
# Look at the intermediate product
bw.geom.landesgrenzen <- st_geometry(bw.landesgrenzen)
plot(bw.SRTM)
plot(bw.geom.landesgrenzen, add=TRUE)

# Add Hillshade
dem.aspect<-terrain(bw.SRTM,
                    opt="aspect")
dem.slope<-terrain(bw.SRTM,
                   opt="slope")
hill<-hillShade(dem.slope,
                dem.aspect,
                40, 270)

# Look at the final product
plot(hill,
     col = grey(0:100/100),
     legend = FALSE,
     main = "DEM Baden-Württemberg")
plot(bw.SRTM,
     col = terrain.colors(100, alpha = 0.50),
     add = TRUE)

With the help of this plot, it can now be very well observed that Baden-Württemberg is a federal state with (by German standards) great topographic contrasts. The Black Forest and the Swabian Alb, which are defined as low mountain ranges, are particularly easy to recognize.The highest point is 1491 m above sea level (southern Black Forest) and the lowest point is only 26 m above sea level (Rhine valley). Try this on your region.The highest point is 1491 m above sea level (southern Black Forest) and the lowest point is only 26 m above sea level (Rhine valley). What are the values of your region?

References:

Horn, B.K., 1981. Hill shading and the reflectance map. Proceedings of the IEEE, 69, 14–47.

Zevenbergen, L.W. & Thorne, C.R., 1987. Quantitative analysis of land surface topography. Earth surface processes and landforms, 12, 47–56.

Previous
Next