Tips and tutorials for GIS and Remote Sensing…

Tag Archives: DEM

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:

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:

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:

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.

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:

Numba is used to compile the Python code – which is a lot faster than pure Python

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:

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.