Create Digital Elevation Models using regional data

Create Digital Elevation Models using regional data#

This script shows how to create Digital Elevation Models (DEMs) with publicly available data. Using a couple of vegetation survey plots as an example, we show how to add additional information to the elevation models.

If you have questions, suggestions, spot errors, or want to contribute, get in touch with us through planthub@idiv.de.

Author: David Schellenberger Costa

Requirements#

To run the script, the following is needed:

a couple of publicly available geodata files provided by the German Federal State of Saxony
1. A digital surface model including vegetation and buildings (DOM, Digitales Oberflächenmodell in German)
  To download the tiles, do the following:
  - open this Saxony state geodata page
  - scroll down on the page to the map view and enter “Winkelmühle, Doberschütz (Ortsteil)” into the upper right search box
  - left of the map, click the dropdown titled “Produkt” and select DOM1
  - click on “Mit Suchergebnis auswählen”
  - click on “Auswahl herunterladen”
  - unpack the downloaded file
2. A digital elevation model excluding vegetation and buildings (DGM, Digitales Geländemodell in German)
  To download the tiles, follow the instruction from 1. and select DGM1 instead of DOM1
3. A digital landscape model (DLM, Digitales Landschaftsmodell in German) including information on streets, agriculture, forests, cities, etc.
  - open this Saxony state geodata page
  - select and unpack the file named “Download DLM50 (Shape) (*.zip, 221 MB)”
field-recorded GPX data showing plot locations available here
some R libraries that may need to be installed

Code#

# load in libraries
library(plotKML)
library(rgdal)
library(raster)

# clear workspace
rm(list = ls())

Let’s read in the data.

# set working directory (adapt this!)
setwd(paste0(.brd, "Vegetationsoekologisches Praktikum"))

# get the field-recorded GPX data
plotGPS <- readGPX("Vegetationsoekologisches Praktikum 2021.gpx")
plotGPS <- plotGPS$waypoints # just keep waypoints, we do not need tracks etc.

We add a geographic reference system (GRS) to the coordinates and transform them into the GRS of our digital elevation models. As the GPX data comes in WGS84 format, we select the EPSG code 4326. EPSG codes are shorthand ways to select specific GRS. As we can see on the Saxony state geodata page (provided we now some German), their data uses the GRS with the ESPG code 25833. So we convert the coordinates into this format.

plotCoords <- SpatialPoints(coords = plotGPS[, 1:2], proj4string = CRS("+init=epsg:4326")) # convert to spatial points
plotCoords <- spTransform(plotCoords, CRS("+init=epsg:25833"))

The DLM comes as a collection of files representing landscape features. Let’s load some interesting ones.

# vegetation
fields <- readOGR("Shape_DLM50_20210406_UTM33/veg01_f.shp")
forests <- readOGR("Shape_DLM50_20210406_UTM33/veg02_f.shp")
shrubs <- readOGR("Shape_DLM50_20210406_UTM33/veg03_f.shp")

# houses
houses <- readOGR("Shape_DLM50_20210406_UTM33/sie02_f.shp")

# infrastructure
roads <- readOGR("Shape_DLM50_20210406_UTM33/ver01_l.shp") # loading this takes a while
paths <- readOGR("Shape_DLM50_20210406_UTM33/ver02_l.shp")

# water
water <- readOGR("Shape_DLM50_20210406_UTM33/gew01_f.shp")

We create vectors storing the file paths to the DGM/DOM files, so we can load them in a loop. Depending on where you unpacked the data, this may need to be adapted. We then store the data in lists and convert them into matrices. We print the range extent of the files.

# names of elevation data without consideration of vegetation and buildings
dgmNames <- c("dgm1_3425712/333425712_dgm1.xyz", "dgm1_3425714/333425714_dgm1.xyz")

# names of elevation data with consideration of vegetation and buildings
domNames <- c("dom1_3425712/333425712_dom1.xyz", "dom1_3425714/333425714_dom1.xyz")

# load DGM files and convert tables into matrices for plotting and get ranges
dgms <- list()
dgmMatrices <- list()
for (i in seq_along(dgmNames)) {
	dgms[[i]] <- read.table(dgmNames[i])
	dgmMatrices[[i]] <- matrix(dgms[[i]][, 3], length(unique(dgms[[i]][, 1])),
		length(unique(dgms[[i]][, 2])),
		byrow = TRUE
	)
	print(range(dgms[[i]][, 1]))
	print(range(dgms[[i]][, 2]))
}

# repeat for DOM files
doms <- list()
domMatrices <- list()
for (i in seq_along(domNames)) {
	doms[[i]] <- read.table(domNames[i])
	domMatrices[[i]] <- matrix(doms[[i]][, 3], length(unique(doms[[i]][, 1])),
		length(unique(doms[[i]][, 2])),
		byrow = TRUE
	)
	print(range(doms[[i]][, 1]))
	print(range(doms[[i]][, 2]))
}

As we downloaded adjacent tiles, it makes sense to join them.

dgm <- rbind(dgmMatrices[[1]], dgmMatrices[[2]])
dim(dgm)
colnames(dgm) <- unique(dgms[[1]][, 1])
rownames(dgm) <- c(unique(dgms[[1]][, 2]), unique(dgms[[2]][, 2]))

dom <- rbind(domMatrices[[1]], domMatrices[[2]])
dim(dom)
colnames(dom) <- unique(doms[[1]][, 1])
rownames(dom) <- c(unique(doms[[1]][, 2]), unique(doms[[2]][, 2]))

We can now remove the data used to create the DGM/DOM files

rm(dgms)
rm(dgmMatrices)
rm(doms)
rm(domMatrices)

The DLM data covers a much larger area than the DOM and DGM data. In fact, we downloaded and loaded data for whole Saxony. We will use the crop() function to reduce the area to the one covered by the DOM and DGM. Before this, we will extract the spatial extent of the DGM by creating a rectangle of the outermost points.

areaExtent <- SpatialPoints(coords = expand.grid(
	range(as.numeric(colnames(dgm))),
	range(as.numeric(rownames(dgm)))
), proj4string = CRS("+init=epsg:25833"))

forests <- crop(forests, extent(areaExtent))
fields <- crop(fields, extent(areaExtent))
shrubs <- crop(shrubs, extent(areaExtent))
houses <- crop(houses, extent(areaExtent))
roads <- crop(roads, extent(areaExtent))
paths <- crop(paths, extent(areaExtent))
water <- crop(water, extent(areaExtent))

Our data is now ready for plotting. Let’s save it into an RData file, so we can use it later.

save.image("Digital terrain model start data.RData")

The DGM/DOM have a very high resolution (1m x 1m). This means that any attempt to plot them as they are will take a lot of time. We will therefore reduce the resolution of the data to speed up computation. This is done by selecting only every xth row or column from the dataset. It is always a good idea to start with a very low resolution and increase it stepwise until we find a good trade-off between computation time and accuracy. As the extent of the DGM/DOM in the x dimension is half of the y dimension (because we selected two tiles that are located one above the other), we should always select half as many columns as rows to avoid distance distortions. We recommend choosing the following values for the res (=resolution) variable: 1000 for printing as a PDF, and 100 for viewing on the screen.

res <- 100
dgmSmall <- dgm[seq(1, nrow(dgm), l = res), seq(1, ncol(dgm), l = res / 2)]
domSmall <- dom[seq(1, nrow(dom), l = res), seq(1, ncol(dom), l = res / 2)]

We only need a part of the area covered by the DGM/DOM to show the locations of our research plots. We select the needed area using the rownames and colnames of the DGM/DOM which are essentially the coordinates.

dgmSmall <- dgmSmall[, as.numeric(colnames(dgmSmall)) >= 342400 & as.numeric(colnames(dgmSmall)) < 343600]
dgmSmall <- dgmSmall[as.numeric(rownames(dgmSmall)) >= 5713800 & as.numeric(rownames(dgmSmall)) < 5715000, ]
domSmall <- domSmall[, as.numeric(colnames(domSmall)) >= 342400 & as.numeric(colnames(domSmall)) < 343600]
domSmall <- domSmall[as.numeric(rownames(domSmall)) >= 5713800 & as.numeric(rownames(domSmall)) < 5715000, ]

We repeat the reduction process for the DLM data

areaExtent <- SpatialPoints(coords = expand.grid(
	range(as.numeric(colnames(dgmSmall))),
	range(as.numeric(rownames(dgmSmall)))
), proj4string = CRS("+init=epsg:25833"))

forests <- crop(forests, extent(areaExtent))
fields <- crop(fields, extent(areaExtent))
shrubs <- crop(shrubs, extent(areaExtent))
houses <- crop(houses, extent(areaExtent))
roads <- crop(roads, extent(areaExtent))
paths <- crop(paths, extent(areaExtent))
water <- crop(water, extent(areaExtent))

We now apply a trick to increase the 3D effect of our models. As they are now, the edges of the models would be transparent. This does not look very nice. To give the model a more solid appearance, we add points with less than the lowest elevation found in the DGM/DOM at the borders of our DGM/DOM. As we only add one row/column of points at each edge, they will not be visible, but the edges will appear to be closed. In this example, the elevation of the outermost points is the minimum of the DGM/DOM minus the elevation range of the DGM/DOM divided by 100. You may experiment with other values.

dgmSmall <- rbind(
	min(dgmSmall) - diff(range(dgmSmall)) / 100,
	cbind(min(dgmSmall) - diff(range(dgmSmall)) / 100, dgmSmall, min(dgmSmall) -
		diff(range(dgmSmall)) / 100), min(dgmSmall) - diff(range(dgmSmall)) / 100
)
domSmall <- rbind(
	min(domSmall) - diff(range(domSmall)) / 100,
	cbind(min(domSmall) - diff(range(domSmall)) / 100, domSmall, min(domSmall) -
		diff(range(domSmall)) / 100), min(domSmall) - diff(range(domSmall)) / 100
)

It’s now time to choose colors for elevation levels of the DGM (it makes no sense for the DOM, as that one includes buildings and vegetation). First, we choose a color gradient. Then, we transfer the colors to a matrix of the same dimension as the DGM.

cols <- terrain.colors(15) # typical elevation colors predefined in R
cols <- cols[1:10] # reduce to the lower ten because the upper one's refer to snowy peaks
colSeq <- seq(min(dgmSmall), max(dgmSmall), l = length(cols) + 1) # choose elevation intervals for color plotting

# create a matrix of the same dimension as the DGM with grey base color
fill <- matrix("gray", nrow(dgmSmall), ncol(dgmSmall))

# fill the matrix depending on elevation
for (i in seq_along(cols)) {
	fill[dgmSmall >= colSeq[i]] <- cols[i]
}
fill <- fill[-(nrow(fill) - 1), -(ncol(fill) - 1)] # fills are one less than grid points in each direction

Let’s create a simple 2D overview plot of the plot coordinates first.

latDgmSmall <- as.numeric(rownames(dgmSmall)) # y range of later plots
lonDgmSmall <- as.numeric(colnames(dgmSmall)) # x range of later plots

plot(plotCoords@coords[, 1], plotCoords@coords[, 2],
	xlim = range(lonDgmSmall, na.rm = TRUE), ylim = range(latDgmSmall, na.rm = TRUE)
)
points(plotCoords@coords[19, 1], plotCoords@coords[19, 2],
	col = "red", pch = 16
) # use other color for base camp than actual plots

We need to create relative coordinates for the DGM/DOM models that basically give the distance between the elevation points. As we have our fake points at the borders, we want to keep the distance of those and the next real points small. We therefore first create equal-spaced sequences of the number of points in x and y directions, subtract two because of the fake points at the borders and at two very tiny distances at the beginning and end. The size of the the fake intervals needs to be adapted depending on the distances between the real points. In this example, 1e-5 seems to be a good choice.

xlab <- seq(0, 1, l = ncol(dgmSmall) - 2)
xlab <- c(-1e-5, xlab, max(xlab) + 1e-5)
ylab <- seq(0, 1, l = nrow(dgmSmall) - 2)
ylab <- c(-1e-5, ylab, max(ylab) + 1e-5)

We can now choose whether we want to plot the DGM or DOM. If we choose the DOM, the fill we be uniform grey, but this can of course be changed.

plotChoice <- 2 # 1=DGM, 2=DOM

if (plotChoice == 1) {
	demSmall <- dgmSmall
	plotCol <- t(fill)
} else {
	demSmall <- domSmall
	plotCol <- "grey"
}

Let’s do the real plotting now! Uncomment the pdf() and dev.off() commands if you want to print your results to a PDF file. We first plot the DGM/DOM. If plotting on the screen, you may wish to check the result already. We then add the landscape features from the DLM. We do this with a loop for each landscape class separately, as each class contains many objects, e.g. houses, roads, etc. Feel free to experiment with the colors or select only some of the features for plotting by commenting out the others. Finally, we add the coordinates of the research plots.

# pdf("Digital terrain model.pdf")  # uncomment this for PDF creation
par(mfrow = c(2, 2))
par(mar = c(1, 1, 1, 1))
par(oma = c(0, 0, 0, 0))

# plotting of the DGM/DOM
q1 <- persp(xlab, ylab, t(demSmall),
	theta = -20, phi = 25, col = plotCol,
	expand = 0.1, border = NA, scale = TRUE, box = FALSE, ltheta = -120, lphi = 45,
	shade = 0.45, axes = FALSE, main = "Winkelm\ufchle"
)

# add terrain information in loops for each component separately
for (i in seq_along(forests@polygons)) {
	hp <- vapply(seq_len(nrow(forests@polygons[[i]]@Polygons[[1]]@coords)), function(x) {
		demSmall[
			order(abs(forests@polygons[[i]]@Polygons[[1]]@coords[x, 2] - latDgmSmall))[1],
			order(abs(forests@polygons[[i]]@Polygons[[1]]@coords[x, 1] - lonDgmSmall))[1]
		]
	}, 1)
	xp <- (forests@polygons[[i]]@Polygons[[1]]@coords[, 1] - min(lonDgmSmall, na.rm = TRUE)) /
		(max(lonDgmSmall, na.rm = TRUE) - min(lonDgmSmall, na.rm = TRUE))
	yp <- (forests@polygons[[i]]@Polygons[[1]]@coords[, 2] - min(latDgmSmall, na.rm = TRUE)) /
		(max(latDgmSmall, na.rm = TRUE) - min(latDgmSmall, na.rm = TRUE))
	temp <- trans3d(x = xp, y = yp, z = hp, pmat = q1)
	polygon(temp, col = "#00880088", border = NA)
}
for (i in seq_along(fields@polygons)) {
	hp <- vapply(seq_len(nrow(fields@polygons[[i]]@Polygons[[1]]@coords)), function(x) {
		demSmall[
			order(abs(fields@polygons[[i]]@Polygons[[1]]@coords[x, 2] - latDgmSmall))[1],
			order(abs(fields@polygons[[i]]@Polygons[[1]]@coords[x, 1] - lonDgmSmall))[1]
		]
	}, 1)
	xp <- (fields@polygons[[i]]@Polygons[[1]]@coords[, 1] - min(lonDgmSmall, na.rm = TRUE)) /
		(max(lonDgmSmall, na.rm = TRUE) - min(lonDgmSmall, na.rm = TRUE))
	yp <- (fields@polygons[[i]]@Polygons[[1]]@coords[, 2] - min(latDgmSmall, na.rm = TRUE)) /
		(max(latDgmSmall, na.rm = TRUE) - min(latDgmSmall, na.rm = TRUE))
	temp <- trans3d(x = xp, y = yp, z = hp, pmat = q1)
	polygon(temp, col = "#88880088", border = NA)
}
for (i in seq_along(shrubs@polygons)) {
	hp <- vapply(seq_len(nrow(shrubs@polygons[[i]]@Polygons[[1]]@coords)), function(x) {
		demSmall[
			order(abs(shrubs@polygons[[i]]@Polygons[[1]]@coords[x, 2] - latDgmSmall))[1],
			order(abs(shrubs@polygons[[i]]@Polygons[[1]]@coords[x, 1] - lonDgmSmall))[1]
		]
	}, 1)
	xp <- (shrubs@polygons[[i]]@Polygons[[1]]@coords[, 1] - min(lonDgmSmall, na.rm = TRUE)) /
		(max(lonDgmSmall, na.rm = TRUE) - min(lonDgmSmall, na.rm = TRUE))
	yp <- (shrubs@polygons[[i]]@Polygons[[1]]@coords[, 2] - min(latDgmSmall, na.rm = TRUE)) /
		(max(latDgmSmall, na.rm = TRUE) - min(latDgmSmall, na.rm = TRUE))
	temp <- trans3d(x = xp, y = yp, z = hp, pmat = q1)
	polygon(temp, col = "#00888888", border = NA)
}
for (i in seq_along(houses@polygons)) {
	hp <- vapply(seq_len(nrow(houses@polygons[[i]]@Polygons[[1]]@coords)), function(x) {
		demSmall[
			order(abs(houses@polygons[[i]]@Polygons[[1]]@coords[x, 2] - latDgmSmall))[1],
			order(abs(houses@polygons[[i]]@Polygons[[1]]@coords[x, 1] - lonDgmSmall))[1]
		]
	}, 1)
	xp <- (houses@polygons[[i]]@Polygons[[1]]@coords[, 1] - min(lonDgmSmall, na.rm = TRUE)) /
		(max(lonDgmSmall, na.rm = TRUE) - min(lonDgmSmall, na.rm = TRUE))
	yp <- (houses@polygons[[i]]@Polygons[[1]]@coords[, 2] - min(latDgmSmall, na.rm = TRUE)) /
		(max(latDgmSmall, na.rm = TRUE) - min(latDgmSmall, na.rm = TRUE))
	temp <- trans3d(x = xp, y = yp, z = hp, pmat = q1)
	polygon(temp, col = "#66666688", border = NA)
}
for (i in seq_along(water@polygons)) {
	hp <- vapply(seq_len(nrow(water@polygons[[i]]@Polygons[[1]]@coords)), function(x) {
		demSmall[
			order(abs(water@polygons[[i]]@Polygons[[1]]@coords[x, 2] - latDgmSmall))[1],
			order(abs(water@polygons[[i]]@Polygons[[1]]@coords[x, 1] - lonDgmSmall))[1]
		]
	}, 1)
	xp <- (water@polygons[[i]]@Polygons[[1]]@coords[, 1] - min(lonDgmSmall, na.rm = TRUE)) /
		(max(lonDgmSmall, na.rm = TRUE) - min(lonDgmSmall, na.rm = TRUE))
	yp <- (water@polygons[[i]]@Polygons[[1]]@coords[, 2] - min(latDgmSmall, na.rm = TRUE)) /
		(max(latDgmSmall, na.rm = TRUE) - min(latDgmSmall, na.rm = TRUE))
	temp <- trans3d(x = xp, y = yp, z = hp, pmat = q1)
	polygon(temp, col = "#00008888", border = NA)
}
for (i in seq_along(roads@lines)) {
	hp <- vapply(seq_len(nrow(roads@lines[[i]]@Lines[[1]]@coords)), function(x) {
		demSmall[
			order(abs(roads@lines[[i]]@Lines[[1]]@coords[x, 2] - latDgmSmall))[1],
			order(abs(roads@lines[[i]]@Lines[[1]]@coords[x, 1] - lonDgmSmall))[1]
		]
	}, 1)
	xp <- (roads@lines[[i]]@Lines[[1]]@coords[, 1] - min(lonDgmSmall, na.rm = TRUE)) /
		(max(lonDgmSmall, na.rm = TRUE) - min(lonDgmSmall, na.rm = TRUE))
	yp <- (roads@lines[[i]]@Lines[[1]]@coords[, 2] - min(latDgmSmall, na.rm = TRUE)) /
		(max(latDgmSmall, na.rm = TRUE) - min(latDgmSmall, na.rm = TRUE))
	temp <- trans3d(x = xp, y = yp, z = hp, pmat = q1)
	lines(temp, col = "#000000AA", border = NA)
}
for (i in seq_along(paths@lines)) {
	hp <- vapply(seq_len(nrow(paths@lines[[i]]@Lines[[1]]@coords)), function(x) {
		demSmall[
			order(abs(paths@lines[[i]]@Lines[[1]]@coords[x, 2] - latDgmSmall))[1],
			order(abs(paths@lines[[i]]@Lines[[1]]@coords[x, 1] - lonDgmSmall))[1]
		]
	}, 1)
	xp <- (paths@lines[[i]]@Lines[[1]]@coords[, 1] - min(lonDgmSmall, na.rm = TRUE)) /
		(max(lonDgmSmall, na.rm = TRUE) - min(lonDgmSmall, na.rm = TRUE))
	yp <- (paths@lines[[i]]@Lines[[1]]@coords[, 2] - min(latDgmSmall, na.rm = TRUE)) /
		(max(latDgmSmall, na.rm = TRUE) - min(latDgmSmall, na.rm = TRUE))
	temp <- trans3d(x = xp, y = yp, z = hp, pmat = q1)
	lines(temp, col = "#444444AA", border = NA)
}

# add plots
hp <- vapply(seq_len(nrow(plotCoords@coords)), function(x) {
	demSmall[
		order(abs(plotCoords@coords[x, 2] - latDgmSmall))[1],
		order(abs(plotCoords@coords[x, 1] - lonDgmSmall))[1]
	]
}, 1)
xp <- (plotCoords@coords[, 1] - min(lonDgmSmall, na.rm = TRUE)) /
	(max(lonDgmSmall, na.rm = TRUE) - min(lonDgmSmall, na.rm = TRUE))
yp <- (plotCoords@coords[, 2] - min(latDgmSmall, na.rm = TRUE)) /
	(max(latDgmSmall, na.rm = TRUE) - min(latDgmSmall, na.rm = TRUE))
temp <- trans3d(x = xp, y = yp, z = hp, pmat = q1)
for (i in 1:6) {
	points(temp$x[1:3 + (i - 1) * 3], temp$y[1:3 + (i - 1) * 3],
		col = rainbow(6)[i], pch = 16
	)
} # use rainbow colors
points(temp$x[19], temp$y[19], pch = 16)
points(temp$x, temp$y)

# dev.off()  # uncomment this for PDF creation

Create Digital Elevation Models using regional data

Contents

Create Digital Elevation Models using regional data#

Requirements#

Code#