LiDAR Digital Elevation Model (DEM) using QGIS and GRASS GIS

Overview

Modeling terrain is a powerful function of a GIS. With the emergence of Airborne LiDAR, being able to turn the millions of highly accurate laser point readings into a surface that can be used to perform analysis on a large scale is invaluable.

As LiDAR data becomes more accessible, many industries and consultants will be able to utilize the data for analysis. It seems, however, that many people shy away from using LiDAR or do not use LiDAR to its full potential due to lack of experience in working with the data or the assumption that the cost of the data is too expensive or the cost of the software required to process it is too complicated or too expensive as well.

As a first step in breaking down the accessibility barriers to LiDAR, this article will walk through the common task of creating a Digital Elevation Model (DEM) or more specifically a Digital Surface Model (DSM) which is a raster/image/gridded representation of the ground topography with the vegetation and buildings removed. The DEM/DSM is the base data that can be used in many types of environmental and terrain modeling such as calculating watersheds, flood areas, viewsheds, slopes and terrain feature identification. In order to accomplish this we will use QGIS 3.4 which has access to the GRASS GIS 7.4 tools.

There are several ways to process LiDAR with free and open source solutions. I like using this method because all of the processing, analysis, visualization and mapping can be done with one environment. This is adventitious as analysis becomes more complicated at which point you can start creating scripts and tools to automate a lot of the time consuming and repetitive processes (a topic for a future tutorial). QGIS is ideal because comes with GRASS tools integrated into it (possibly not for macOS)

The data set we will use for this tutorial is from Open Topography and is available at this link.

Open Topography is an awesome source for elevation information numerous reasons. Not only does it provide LiDAR data at no cost but it will provide you with many of the computationally intensive elevation products, such as Hillshades, DEMs, Point Cloud Viewer and much much more!

However, for the purpose of this tutorial we will choose to download only a small portion of the LiDAR data, unselecting all of the other options that Open Topography provides.

The image below shows the parameters used to download that LiDAR data be sure to select .las format!

Step 1

Inspect the LiDAR metadata to find out what Spatial Reference System (SRS) the LiDAR data is in, to get an idea of how any points per square meter there are on average and see how the LiDAR data has been classified (i.e. ground returns and non-ground returns). To get a better understanding of how LiDAR data is structured take a look at the GISGeography website has an excellent description of how LiDAR works and ASPRS - The Imaging and Geospatial Information Society has set the standard for how LiDAR should be classified.

Below is the meta data that can be downloaded along with the point cloud data from Open Topography. The most important details to note at this point is the horizontal and vertical coordinate systems (specifically the EPSG codes).

Step 2

In this next step we will use some statistics go gain a little more insight into the point cloud data. Our goal is to see what what we can use as a maximum resolution for our DEM. To do this we will calculate how many LiDAR ground points there are per square meter using GRASS GIS' r.in.lidar tool and the built in QGIS Raster Information... tool.

Open a new QGIS project (download and install a free copy now if you have not done so already! Once the program is open make sure you set the default spatial reference system/projection/EPSG to match that of the LiDAR dataset. Next, click on the gear shaped "Toolbox" button in the tool bar to open the "Processing Toolbox". In the search box enter r.in.lidar

The key parameters to enter into the tool are:
LAS input file=points.las
Statistics to use for raster values=n (i.e. number of lidar points per raster cell)
Storage type for resultant raster map=CELL
Filter range for z data=-10000, 10000 (use large values to ensure all points fall within the range)
Filter range for intensity values=-10000, 10000
Output raster resolution=1
Only import points of selected class(es)=2
Lidar Raster=/output/file/path/lidar_density.tif
Advanced parameters
Use the extent of the input for the raster extent=Checked
Set the computation region to match the new raster map=Checked
Override projection check=Checked

The output raster should look similar to the image below and at first look will not give us much information. In step 3 we will perform some simple statistics on this raster to help us choose the optimal resolution for our elevation products.

Step 3

Using the built in QGIS Raster Information... tool. The tool is located in the drop-down menu at the top of the QGIS window: Raster -> Miscellaneous -> Raster Information...
Using the raster generated in Step 2 as the input layer, with Force computation of the actual min/max values for each band checked and Read and display image statistics (force computation if necessary) checked we will compute the raster statistics.
The tool will output a .html file that should be saved along with the other outputs and data from this tutorial for future reference.
A snippet of the output can be seen below

The image below is the LiDAR density raster that has been reclassified to only the cells that do not have any point measurements for that area. This is important for two reasons. First, we see what areas of the terrain are not returning any LiDAR points so we can get a sense of where the DEM might be less accurate. Second, we can use this visualization to adjust our DEM resolution to ensure there are an adequate number of LiDAR points per raster cell throughout the study area.

Step 4

Convert .las to GeoPackage

In order for a DEM to be created from LiDAR points we have to convert the data from .las format to GeoPackage format. This step is the most computationally intensive process and will require some additional hard drive space as the GeoPackage format does not store the point cloud as compactly as LAS. For example, the point cloud used for this tutorial is 650 MB in size in .las format and is 808 MB. This could partially be due to the fact that a topology is built during the point cloud import and is also likely why the process is computationally intensive. This does make sense in the long run as future calculations will be a lot more efficient once the GeoPackage is created.

Import point cloud with v.in.lidar

In your open QGIS project search for v.in.lidar in the Processing Toolbox panel or browse for it under GRASS -> Vector (v.*) -> v.in.lidar. There are input parameters that are key for the import. First is the Filter range for z data [optional]. This is not optional! If you do not set these values above and below the elevation range of the point data not all of the point data will be imported! Second, Only import points of selected class(es) (comma separated integers) [optional] should be set to 2 to only import the ground classified point. The rest of the parameters can be left as default (see image below).

Step 5

In this step we will use the GRASS GIS v.surf.idw tool to create a DEM. We know from the statistics performed in the previous steps that the majority of the cells have enough point readings per cell (~4) to use a 1 meter resolution for the DEM. Also, because the points are irregularly dispersed across the study area and there are significantly sized areas that do not have elevation readings, we will use a mathematical method called Inverse Distance Weighting (IDW) to calculate an appropriate elevation value for each raster cell. GISGeography is an excellent resource for understanding how IDW works.

There are several parameters that can be tweaked to affect the resulting DEM depending on what the data will be used for. The most significant variable is the Power parameter. A higher power will increase the weight of elevation points closer to a given cell meaning better accuracy for cells with lots of LiDAR points close by but results in a much rougher surface. A lower power will result in a smoother surface which might be contagious when producing hillshade rasters and surface flow models.

The next key parameter is cell size. The resulting DEM will likely be the base for lots of future analysis and models, so the cell size must be chosen carefully to try and find the balance of highest possible resolution, processing demands, project scale and accuracy demands. A 1m resolution is quite common and is flexible in that it can be down sampled to a lower resolution later and is appropriate for most industries.

The parameters to generate a good all round DEM can be seen in the image below. The key parameters that were changed from the default values are:
Number of interpolation points:Delete Default value of 12 leaving value of "Not Set"
Attribute Table Column with values to interpolate:z_coord
GRASS GIS 7 region cellsize:1 *This is the resulting raster cell size