Creating a DEM from LiDAR data using the ARSF DEM Scripts

The Airborne Survey and Research Facility Data Analysis Node (ARSF-DAN) based at Plymouth Marine Laboratory (PML) process airborne hyperspectral and lidar data acquired by the ARSF. As part of the LiDAR pre-processing a Digital Surface Model (DSM) is produced from discrete lidar returns, patched with a lower resolution DSM (normally ASTER) suitable for use in APL for hyperspectral data processing. The scripts used to produce these DSM use GRASS. Updated versions, which use the GRASS Python bindings, have recently been made available on GitHub under a GPLv3 license:

https://github.com/pmlrsg/arsf_dem_scripts

Installation

There are two main pre-requisites for the ARSF DEM Scripts: GRASS and the open source LAStools (namely las2txt), both are available under Windows, Linux and OS X. Once these are installed, download the scripts from GitHub (direct link) and install using:

python setup.py install

As the scripts use the GRASS Python bindings they need to be run using the same version of Python used by GRASS, this will likely be Python 2.7. If you have Python 3 installed as your default Python you should be able to specify they use Python 2.7 by installing using:

python2 setup.py install

For more detailed installation instructions see here.

Create a DSM Mosaic

To create a mosaic of all LAS files the following command is used:

create_dem_from_lidar.py --in_projection UTM33N \
                      --outdem EUFAR11_02-2011-187_dsm.dem \
                      las1.2

This will:

1. Set up a GRASS database.
2. For each line, convert the LAS file to a temporary text file using las2txt, keeping only the first returns and dropping points flagged as noise (class 7).
3. Import the text file into GRASS using r.in.xyz.
4. Patch all LAS files together using r.patch.
5. Export the mosaic using r.out.gdal.
6. Remove the GRASS database (unless –keepgrassdb is specified) and any other temp files created.

The library supports horizontal and vertical transforms by setting the ‘–out_projection’ flag. For transforms from British National Grid some additional files are required to ensure and accurate transform.

• The OSTN02 transform file, which can be downloaded from Ordnance Survey for horizontal transforms.
• A vertical difference file between the WGS-84 ellipsoid and the Ordnance Survey datum – ARSF can provide a copy of this ready for use in the DEM scripts.

The location of these needs to be set in the arsf_dem.cfg file. This is installed with the library but can be overridden by placing a copy in the current working directory (using the current name) or the home folder (stored with the name ‘.arsf_dem’ or ‘.arsf_dem.cfg’ so it is hidden).

DSM / DTM Utility Scripts

In addition to the scripts used as part of the standard ARSF processing, based on GRASS, there are two utility scripts to create a DSM and Digital Terrain Model (DTM). These two scripts (las_to_dsm and las_to_dtm) provide a common interface to a number of open source (e.g., SPDLib, FUSION) and paid packages (e.g., LAStools). If SPDLib is installed a DTM/DSM can be created using:

las_to_dsm -o LDR-EUFAR11_02-2011-187-01_spdlib_dsm.tif \
           --hillshade LDR-EUFAR11_02-2011-187-01_spdlib_dsm_hillshade.tif \
           --projection UTM33N  \
           --method SPDLib \
           LDR-EUFAR11_02-2011-187-01.LAS

las_to_dtm -o LDR-EUFAR11_02-2011-187-01_spdlib_dtm.tif \
           --hillshade LDR-EUFAR11_02-2011-187-01_spdlib_dtm_hillshade.tif \
           --projection UTM33N  \
           --method SPDLib \
           LDR-EUFAR11_02-2011-187-01.LAS

Note, these utility scripts create a DSM/DTM using the default settings. For more control over output accessing the programs directly is advised. See previous post for an example of creating a DSM/DTM using SPDLib.

For more details on the use of the ARSF DEM scripts, including creating patched DSM for use in APL see the tutorial.

Archived LiDAR data from ARSF flights is available to download from NEODC. Registered users can apply for access to the ARSF archive here.

Create a slope layer from a DEM

To calculate slope from a DEM the gdaldem utility can be used:

gdaldem slope -of KEA in_dem.kea out_slope.kea

There is an option to pass in a scale to convert the horizontal spacing from degrees to metres. However, if the DEM covers a range of latitudes it is desirable to convert the horizontal spacing on a per-pixel basis. To accomplish this I wrote a Python script, with RIOS used to read and write out data in blocks, making it efficient to process large datasets.

Implementation
As the slope calculation requires looping through a block of data, a task which is slow in pure Python, two strategies are used to improve the speed:

1. Numba is used to compile the Python code – which is a lot faster than pure Python
2. Fortran code is used to perform the actual slope calculation, compiled as a Python module using f2py

To provide a comparison of the time required for each implementation slope was calculated from a 1800 x 1800 pixel DEM using each method:

Implementation	Time (s)
Pure Python	136.7
Numba	2.6
Fortran	2.1

Tests using Python 3.4 with Numba 0.13.4 running under OS X with 1.8 GHz Intel i5 processor. The gfortran Fortran compiler was used. Average runtime over 3 runs.

So the Fortran version is slightly faster than the Numba version but only slightly and required a lot more effort to implement. Both are significantly faster than the pure Python version. The code provides a nice example of applying more complicated algorithms to images with RIOS used to handle the I/O.

To account for noise in the DEM there is also a version of the slope calculation which will use least squares fitting to fit a plane over a window of pixels and calculate the slope from this. As the Python version requires calls to the NumPy linear fitting code there is no improvement using Numba. The Fortran version uses Accelerate under OS X or ATLAS under Linux for the plane fitting.

The code is available to download from https://github.com/MiXIL/calcSlopeDegrees

Usage
To calculate slope from a DEM where the horizontal sloping is in degrees:

calcSlopeDegrees.py --spacing_degrees inDEM.kea outSlope.kea

To calculate slope by fitting a plane over 9 x 9 window:

calcSlopeDegrees.py --plane_ls --window_size 9 inDEM.kea outSlope.kea