This post describe the generation of high-resolution DEM from LiDAR data for the county Kall in North Rhine-Westphalia in Germany. The dataset is collected and provided by Geobasis NRW and can be downloaded here. For data processing I use the LidaR package in R.
Introduction
The topography affects various key processes in the earth’s system such as precipitation, erosion, sedimentation, flooding, wind currents, etc. Therefore, topographic information is required for many environmental modelling techniques and Digital Elevation Models (DEM) are usually the choice to go with. Several DEMs with almost global coverage are freely available (e.g. ASTER, AW3D 30m, SRTM). With a spatial resolution of up to 30m these datasets are great for large-scale objectives but much information is lost on the local scale. Using a high-resolution DEM for local-scale objectives (e.g. flood simulations) can significantly increase the quality of the assessment and it is surprisingly simple to generate a very detailed and accurate DEM from LiDAR data.
LiDAR sensors measure the time between emitting and receiving an electromagnetic pulse (the runtime). The position of the reflecting object is then calculated from the runtime, the pulse’s speed and the position of the sensor. Pulses can be reflected at the surface of vegetation but also penetrate it and reflect at lower layers so that this technology is very good for measuring forests. The received points are stored in 3D point clouds.
A bunch of software for the processing of point clouds is available; one user-friendly option is LAStools which can be accessed via QGIS and enables the fast processing of large datasets. However, it’s a commercial tool and needs licensing for commercial usages. Furthermore, also GRASSGIS offers LiDAR processing capabilities (also accessible via QGIS). In this post, I will use the LidR package in R to generate the following datasets:
- Digital Terrain Model (DTM): A topographic model of the bare earth
- Digital Surface Model (DSM): A topographic model of the terrain and all objects (e.g. houses, trees)
- Canopy Height Model (CHM): A model of all objects excluding the terrain.
Exploratory data analysis
Let’s start with getting the hand on the data and examine the dataset. First, we load the package, read the data and print a summary of the dataset.
library(lidR)
# load data
las <- readLAS('C:/DopeData/Data/LiDAR/Kall/Test/Raw/3dm_32_322_5593_1_nw.laz')
# plot points
plot(las)
The point density is an important statistic that describes the level of detail in the dataset and should be considered when choosing the resolution of the DEM. At a given resolution the uncertainty in the DEM increases with decreasing point density. In this dataset, the average number of 13.4 points per sq. m. To investigate the point density in more detail we divide the study area into 1x1m grid cells and calculate the number of points per cell.
# create raster of the point count for 1m resolution
d <- grid_density(las, res = 1)
# Print some statistics
summary(d)
From the summary function we can see that around 8k pixels have no point observation. When increasing the resolution to 0.5 m the number rises to 85k pixels. To get a better picture about the spatial distribution of the point density we can simply plot the density grid.
library(tmap)
library(tmaptools)
# Plot grid density
tm_shape(d) +
tm_raster(style= "cont", palette= viridisLite::viridis(20)) +
tm_layout(legend.outside = TRUE)
# For an overview of the colour palettes run:
library(shinyjs)
palette_explorer()
From the point density map we can see that the density follows the spatial patterns of the land cover. In fact this raster could be used to classify forest, grasslands and water areas.
DTM generation
To generate a Digital Terrain Model we can call the function grid_terrain() and specify a grid resolution and an interpolation method that estimates the elevation of unobserved cells from the neighbouring cells. In this case, the k-nearest inverse distance weighting method is used with the default parameters (k = 10, power = 2). To speed up processing and prevent very large file sizes a resolution of 1m is chosen.
# Create the DTM
dtm <- grid_terrain(las, res = 1, knnidw(k = 10, p = 2), keep_lowest = TRUE)
#Plot it
tm_shape(dtm) +
tm_raster(style= "cont", palette=get_brewer_pal("Greys", plot=FALSE))
The DTM stores the average terrain elevation per 1m grid cell. However, while it is a good format for further processing, it is not great for visualization. Therefore, we will generate a hillshade raster.
# Create Hillshade
slope <- terrain(dtm, opt='slope')
aspect <- terrain(dtm, opt='aspect')
hillshade <- hillShade(slope, aspect, angle=45, direction=315)
# Plot it
tm_shape(hillshade) +
tm_raster(style= "cont", palette=get_brewer_pal("Greys", plot=FALSE))+
tm_layout(legend.show = FALSE)
DSM generation
The DTM contains information about the terrain elevations. Since LiDAR data contains information about the elevation of surface objects too, we can use it to generate a DSM which store the elevation of the terrain + the surface objects. For this, we can use the function grid_canopy() and select an algorithm for the DSM generation (here dsmtin).
# Create DSM
dsm <- grid_canopy(las, res = 1, dsmtin(max_edge = 0))
# Create Hillshade
slope <- terrain(dsm, opt='slope')
aspect <- terrain(dsm, opt='aspect')
hillshade <- hillShade(slope, aspect, angle=45, direction=315)
# Plot it
tm_shape(hillshade)+
tm_raster(style= "cont", palette=get_brewer_pal("Greys", plot=FALSE))+
tm_layout(legend.show = FALSE)
CHM generation
Since we now know the elevation of the terrain and the terrain + surface objects we can easily calculate the height of the surface objects by subtracting the DTM from the DSM. To prevent negative values we set the min elevation to 0.
# Create Canopy Height Model
ndsm <- dsm – dtm
ndsm[ndsm<0] <- 0
# Plot it
tm_shape(ndsm)+
tm_raster(style= "quantile", n=7, palette=get_brewer_pal("Greens", n=7, plot=FALSE))+
tm_layout(legend.outside = TRUE)
Automatization
LiDAR files can become quite big for larger areas which is why the data is usually split into smaller tiles, in this case, 1x1 km tiles. For the county Kall the data is split into 91 tiles of which we only processed one tile so far. To generate a DTM (or DSM) for the whole county we can automate the processing and so that we don’t have to repeat the process 91 times. For this, we list all LiDAR datasets in a directory and iterate through the files. After processing a tile we merge it with the previously processed tiles.
# load libraries
library(rgdal)
library(raster)
library(lidR)
# Define input and output directory
in_dir <- 'D:/DopeData/Data/Lidar/Kall/Raw/'
out_dir <- 'D:/DopeData/Data/Lidar/Kall/DEM_R/'
# list las files in input directory
file_list <- list.files(path = in_dir, pattern = "*.laz")
# Iterate through the files
for (i in 1:length(file_list)) {
print(paste("Processing file: ", i, " of: ",length(file_list), sep = "
"))
# load LiDAR tile
las <- readLAS(paste(in_dir, file_list[i], sep = ""))
# create DTM (1m resolution)
dtm <- grid_terrain(las, res = 1, knnidw(k = 10, p = 2))
# Merge tile with previously processed tiles
if (i == 1) {
dtm_merge <- dtm
} else {
dtm_merge <- merge(dtm_merge, dtm)
}
# Remove files from memory
rm(dtm, las)
}
# write raster file
filename <- paste(out_dir, "DTM_Merge.tif", sep = "")
writeRaster(dtm_merge, filename, overwrite = TRUE)
Commenti