Tips and tutorials for GIS and Remote Sensing…

Category Archives: LiDAR

The first post of this series was written quite a while ago now. Apologies it has taken so long for a follow up. Since the first post has been written there have been two exciting developments:

The methods described for generating full waveform metrics have been used to perform the LiDAR analysis for a paper led by Chloe Brown, University of Nottingham.Brown, C.; Boyd, D.S.; Sjögersten, S.; Clewley, D.; Evers, S.L.; Aplin, P. Tropical Peatland Vegetation Structure and Biomass: Optimal Exploitation of Airborne Laser Scanning. Remote Sens. 2018, 10, 671. https://doi.org/10.3390/rs10050671

SPDLib is now available on Windows, macOS and Linux through conda-forge (as is RSGISLib). See part one for updated install instructions.

At the end of part one we had imported the LAS 1.3 file into SPDLib and decomposed the waveforms. This next section will cover ground classification and metrics generation.

Spatially index data

In part one we had been working with an SPD file without a spatial index (UPD file). However, for subsequent processing steps a spatial index is needed so a spatially indexed file is generated using the spdtranslate command.

As the gridding can use a lot of RAM we are going to process in tiles and then stitch them together. We create a temporary directory to store the tiles using:

The temporary directory can be removed after processing has completed:

rm -fr spd_tmp

Classify ground returns and populate height

Many metrics use the height about ground rather than absolute elevation so this must be defined. To derive heights from LiDAR data it is first necessary to determine the ground elevation so heights can be calculated above this. Within SPDLib the ground classification results are achieved using a combination of two classification algorithms: a Progressive Morphology Filter (PMF; [1]) followed by the Multi-Scale Curvature algorithm (MCC; [2]). Both these algorithms use only the discrete points rather than the waveform information.

The final step of the SPD processing is to attribute each pulse with heights above ground level. An interpolation is used for ground points, similar to generating a Digital Terrain Model (DTM), but rather than using a regular grid the ground height is calculated for the position of each point.

After all the pre-processing steps to convert the LAS 1.3 file into a gridded SPD format file with a defined height it is possible to generate a number of metrics from the waveform data. The command to calculate metrics within SPDLib (`spdmetrics`) takes an XML file in which the metrics are defined. There are a large number of metrics available and operators (addition, subtraction etc.,) allowing existing metrics to be combined to implement new metrics. The full list of metrics is available in the The full list of metrics is available in the SPDMetrics.xml file, distributed with the source of SPDLib. Most metrics have an option to specify the minimum number of returns (`minNumReturns`), setting this to 0 will use the waveform information to calculate the metric, setting to 1 (default) or above will use the discrete data. In this way full waveform and discrete metrics can be created at the same time.For this exercise we will be calculating Height of Medium Energy (HOME) and waveform distance (WD), a detailed description of these metrics is given in [3].

First, create a file containing these metrics. Create a text file called ‘spd_metrics.xml’ and paste the text below into it:

More metrics can be added to the ‘spd_metrics.xml’ file as needed, it is also possible to define new metrics using the operator tags.

This post was derived from the LiDAR practical given as part of the NERC-ARF workshop held at BAS, Cambridge in March 2018. If you have any questions about working with NERC-ARF data contact the NERC-ARF Data Analysis Node (NERC-ARF-DAN) see https://nerc-arf-dan.pml.ac.uk/ or follow us on twitter: @NERC_ARF_DAN.

The SPD file format was designed around storing LiDAR pulses with digitised waveforms and associated points. The most recent version (3.3) has the ability to import waveform data from LAS files using LASlib, which is part of LAStools. Binaries are available for Linux, macOS and Windows, they can be installed through conda.

You can also follow through with any of the other NERC-ARF datasets or other waveform LAS files you have.

Convert LAS to SPD format
First you need a WKT file to define the projection. This step is optional but is more reliable than reading from a LAS file. For the example the projection is UTM50N, you can download a WKT file using:

One of the limitations of discrete systems is there is are only a given number of ‘points’ recorded (normally 2 – 4) and the rest of the information is lost. As full waveform data records the entire waveform it is possible to extract more returns after data are acquired. A common approach to this ‘Gaussian Decomposition’ which involves fitting Gaussian distributions to the peaks, within SPDLib this is available as the ‘spddecomp’ command.

A good way of visualising LiDAR point clouds is to attribute them with colour information.
LAS files of point type 2 and above have the ability to store RGB colour information. These can be viewed with programs such as plas.io and CloudCompare.

The script uses a combination of laspy (to read and write a LAS file) and GDAL (to read the raster) and is capable of working with images with more than 3 bands such as multispectral or hyperspectral files without creating a three band copy of the image first.

For this post I’ll demonstrate running the script on the JASMIN/CMEMS system. In a nutshell it is a Linux computing facility with lots of fast storage and direct access to many datasets, including the ARSF archive of airborne LiDAR and Hyperspectral data (apply for access to ARSF data here). Researchers in the UK who are part of the NERC or National Centre for Earth Observation (NCEO) scientific communities can apply for access. If you don’t have access to JASMIN you can run on your own machine just ignore the JASMIN specific steps.

Apply for access and logon to JASMIN
If you are using JASMIN the first stage is to get an account and log on by following the steps in this guide: http://jasmin.ac.uk/workflow/

Set up scripts
Once you have logged onto one of the shared VMs (e.g., cems-sci2.cems.rl.ac.uk) you need to load laspy, which isn’t installed by default. You can do this using:

module load contrib/arsf/laspy/1.3.0

Then checkout the repository containing the ‘colour_las_file.py’ script using:

git clone https://github.com/pmlrsg/arsf_tools.git

Check the script runs OK by typing:

python2.7 arsf_tools/colour_las_file.py -h

Note, you need to specify Python 2.7 as the default system Python is 2.6 but most of the Python libraries (e.g., GDAL) are build against 2.7. If you are not running the scripts on JASMIN just using ‘python’ should be fine.

Colour using ARSF hyperspectral dataFor this example data from Mont Blanc flown as part of a EUFAR campaign in 2015 on Julian day 175a (24/06/2015) will be used. This data can be downloaded from NEODC (direct link).

This will attribute using bands 38, 25 and 12 (true colour) and will scale the pixel values from 16 bit to 8 bit (0 – 255) using a standard deviation stretch.To read the hyperspectral data without unzipping a GDAL virtual file system is used, as described in a previous post.

LASzip files and Download
To reduce the size of the files before downloading you can use LASzip. On JASMIN LAStools can be loaded using:

module load contrib/arsf/lastools/20150925

To compress the LAS file and drop points flagged as noise the following command can be used:

Open http://plas.io/ in a modern browser and click ‘Browse’ to load in the coloured LAS file.
The nice thing about using plas.io and JASMIN is no specialist software needs to be installed on your local machine, just an ssh client and a web browser.

View using Cloud Compare Viewer

If you want to view files locally you can use Cloud Compare, which is available for Windows, OS X and Linux. To open the coloured LAS/LAZ file in the cloud compare viewer just drag it into the viewer window.

Additional – Colour using an existing three band image
If you already have a three band image, with pixel values 0 – 255 a simpler command can be used:

The Environment Agency have just made their high resolution LiDAR-derived Digital Surface Model (DSM) and Digital Terrain Model (DTM) data available under the Open Government Licence through http://environment.data.gov.uk/ds/survey

The files are downloaded as zipped archives containing ASCII format files. To create a mosaic from them I have created a script using:

The archive reader functionality from TuiView, which can get a list of files within a zip archive and convert them to paths which can be read with GDAL using a virtual filesytem as described in an earlier post

The ‘createImageMosaic’ function in RSGISLib to create the mosaic, changing the no data value from -9999 to 0

The ‘assignProj’ function in RSGISLib to assign the correct projection to the mosaic

The ‘popImageStats’ function in RSGISLib to calculate stats and overviews for fast display in TuiView and other programs

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.

It is common to get LiDAR data (from airborne or terrestrial systems) in ASCII format. However, the format of data within the file (which data is in which column, delimiter etc.,) often varies. To account for these differences SPDLib uses an XML scheme to define which data is in which column and the delimiter. There are some example schemas provided with in the ‘schemas’ folder with the SPDLib source

A common product of LiDAR data is high resolution Digital Elevation Models (DEM) which are a raster (gridded) product. There are two types of DEM: a Digital Terrain Model (DTM) is a model of the bare earth and doesn’t contain trees or buildings, a Digital Surface Model (DSM) is a model of the surface which includes the top of buildings and trees etc. The difference between these can provide the height of trees and buildings.

The Sorted Pulse Data (SPD) library is a format for storing, and a set of tools for processing, full waveform and discrete return LiDAR data. More details on SPDLib are available from spdlib.org and in the following publications:

The spdtranslate command can also be used with data in ASCII format. For full waveform data in LAS 1.3 format a separate script is needed – there will be more details on working with waveform data in SPDLib in a future post.

If the LAS file doesn’t have the projection defined properly you can specify it by setting –input_proj and –output_proj to point to WKT file with the correct projection.

Using the -b option creates a spatially indexed point cloud, which is needed for display and many of the other algorithms in SPDLib. In this example a bin size of 10 m is used.

Classify ground returns

There are multiple algorithms in SPDLib to classify ground returns, one of the recommenced algorithms is a progressive morphology filter [1]. This is available through the command spdpmfgrd.

Within SPDLib there is also an implementation of the multi-scale curvature algorithm [2], created at the US Forest Service. This does a good job at classifying ground returns under a forest canopy while retaining the terrain but it does not differentiate the buildings.

spdmccgrd -r 50 --overlap 10 -i LiDAR_10m.spd -o LiDAR_10m_mccgrd.spd

More details on both these algorithms and the required parameters are available in the SPDLib Documentation and their respective references

Interpolate to DTM and DSM

Once the points have been classified they need to be interpolated to create a raster DEM. A key parameter is the resolution of the raster which is generated (-b) which needs to be a whole number multiple of the SPD index. For example, if the SPD file has a bin size of 10 m then the the output raster file resolution can be 1, 2 or 5 m but not 3 m. There are different interpolation algorithms available in SPDLib, the Natural Neighbour algorithm is recommended.