Tag Archives: SPDLib

Working with full waveform LiDAR data in SPDLib (part 1)

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. Currently binaries are only available for Linux (although binaries for OS X are planned), they can be installed through conda.

conda install -c rios -c conda-forge spdlib

Windows users are advised to install using a Linux Virtual Machine. This post also has more details on installing conda.

For this example LAS 1.3 files acquired by the NERC Airborne Research Facility (NERC-ARF, previously ARSF) over Borneo using a Leica ALS50-II instrument with full waveform digitiser will be used.

Once you have registered with NEODC and applied for access to the ARSF archive these data can be downloaded from: http://browse.ceda.ac.uk/browse/neodc/arsf/2014/RG13_06/RG13_06-2014_303_Maliau_Basin/LiDAR/flightlines/fw_laser/las1.3

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

  1. 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:

    wget https://bitbucket.org/petebunting/rsgis_scripts/raw/c8cf94528cdb58b753029df3bc631a2509740ad1/WKT/UTM_WGS84/UTM_WGS84_Z50_N.wkt
    

    Then convert to an unsorted SPD file (UPD).

    spdtranslate --if LAS --of UPD \
                 -x LAST_RETURN \
                 --input_proj UTM_WGS84_Z50_N.wkt \
                 -i LDR-FW-RG13_06-2014-303-05.LAS \
                 -o LDR-FW-RG13_06-2014-303-05.spd
    
  2. Subset SPD file (optional)

    As full waveform processing is quite intensive it is recommended to subset the data for the purpose of running though this tutorial, you can do this using the spdsubset command.

    spdsubset --xmin 494400 --ymin 524800 \
              --xmax 494800 --ymax 525000 \
              -i LDR-FW-RG13_06-2014-303-05.spd \
              -o LDR-FW-RG13_06-2014-303-05_subset.spd
    
  3. Decompose waveform

    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.

    spddecomp --all --noise --threshold 25 \
              -i  LDR-FW-RG13_06-2014-303-05_subset.spd \
              -o  LDR-FW-RG13_06-2014-303-05_subset_decomp.spd
    

    This will still take around 5 minutes to run. If you decide to decompose the full dataset after, expect it to take an hour or so.

  4. Export returns to LAS file
    The final step for this part of the tutorial is to export the returns to a LAS file, using the spdtranslate command.

    spdtranslate --if SPD --of LAS \
                 -i LDR-FW-RG13_06-2014-303-05_subset_decomp.spd \
                 -o LDR-FW-RG13_06-2014-303-05_subset_decomp.las
    

The next part of the tutorial will cover more tools in SPDLib for classifying ground returns and calculating LiDAR metrics.

If you have further questions about using SPDLib please contact the mailing list (details available from https://lists.sourceforge.net/lists/listinfo/spdlib-develop). 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).

Import ASCII format LiDAR data to SPDLib

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

As an example to demonstrate reading an ASCII file into SPDLib, LiDAR data from NERC’s Airborne Research & Survey Facility (ARSF) is available to download from http://neodc.nerc.ac.uk. You just need to sign up for a NEODC account and then apply for access to the ARSF data.

The ASCII lidar files supplied by ARSF are space separated and contain the following columns:

Time, Easting, Northing, elevation, intensity, classification, return number, number of returns for given pulse, scan angle rank

Note: the ASCII files are only used as an example here, ARSF also supply data in LAS format which is recommended for importing into SPDLib

To import the points with associated intensity and classification the following schema is used:

<?xml version="1.0" encoding="UTF-8" ?>
<line delimiter=" " comment="#" ignorelines="0" >
    <field name="X" type="spd_double" index="1" />
    <field name="Y" type="spd_double" index="2" />
    <field name="Z" type="spd_float" index="3" />
    <field name="AMPLITUDE_RETURN" type="spd_uint" index="4" />
    <field name="CLASSIFICATION" type="spd_uint" index="5" />
</line>

This is saved as ‘arsf_ascii_lidar.xml’ and is passed into SPDLib using the ‘–schema’ flag:

spdtranslate -i LDR-BAS13_01-2014-063-20.txt \
             --if ASCII --schema arsf_ascii_lidar.xml \
             -o LDR-BAS13_01-2014-063-20.spd \
             --of SPD \
             -b 2

Create a DTM and DSM from LAS Format LiDAR data using SPDLib

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:

  • Bunting, P.J., Armston, J., Lucas, R. M., Clewley, D. 2013. Sorted pulse data (SPD) library. Part I: A generic file format for LiDAR data from pulsed laser systems in terrestrial environments. Computers & Geosciences, 56, pp.197–206. (link)
  • Bunting, P.J., Armston, J., Clewley, D., Lucas, R. M. 2013. Sorted pulse data (SPD) – Part II: A processing framework for LiDAR data from pulsed laser systems in terrestrial environments. Computers & geosciences, 56, pp.207–215. (link)

For Linux the software can be installed through conda. Once anaconda / miniconda has been installed (Python 3.5) SPDLib and the prerequisites can be installed using:

conda create -n spdlib_env \
   -c conda-forge -c rios -c osgeo \
   spdlib spd3dpointsviewer tuiview
source activate spdlib_env

To create the DTM and DSM the following steps are required.

  1. Convert to SPD Format
  2. Assuming discrete LiDAR data in LAS 1.2 format:

    spdtranslate --if LAS --of SPD -b 10 -x LAST_RETURN \
    -i LiDAR.las -o LiDAR_10m.spd
    

    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.

  3. Classify ground returns
  4. 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.

    spdpmfgrd -r 50 --overlap 10 --initelev 0.1 --maxfilter 14 -b 0.5 \
    -i LiDAR_10m.spd -o LiDAR_10m_pmfgrd.spd
    

    To view the classified points open the file in SPDPoints viewer, and select ‘Classification’ as the point colour.

    SPDPointsViewer_ground_classification

    Note – you might have problems running SPDPointsViewer if you’re running Linux through VirtualBox.

    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

  5. Interpolate to DTM and DSM
  6. 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.

    # DTM
    spdinterp --dtm --topo -r 50 --overlap 10 --in NATURAL_NEIGHBOR \
    -f GTiff -b 1 -i LiDAR_10m_pmfgrd.spd -o LiDAR_1m_dtm.tif
    
    # DSM
    spdinterp --dsm --topo -r 50 --overlap 10 --in NATURAL_NEIGHBOR \
    -f GTiff -b 1 -i LiDAR_10m_pmfgrd.spd -o LiDAR_1m_dsm.tif
    

    The DTM and DSM can be opened in TuiView. To see the differences you can create a colour composite and look at a profile (the DSM and DTM will be shown in a different colour)

    # Create virtual raster stack
    gdalbuildvrt -separate LiDAR_1m_dtm_dms_stack.vrt \
        LiDAR_1m_dtm.tif LiDAR_1m_dsm.tif
    
    # Open composite in TuiView
    tuiview --rgb -b 2,1,1 --stddev LiDAR_1m_dtm_dms_stack.vrt
    

    DTM_DSM_screengrab

    You can also create a hillshade image using GDAL:

    gdaldem hillshade -of GTiff LiDAR_1m_dtm.tif \
    LiDAR_1m_dtm_hillshade.tif
    

    And view in TuiView:

    Hillshade_screengrab

    Using the spdinterp command it is also possible to interpolate a canopy height model (CHM) to provide the height of a forest canopy using –chm.

    This post was based on the SPDLib documentation written by Pete Bunting. A PDF of the documentation and example datasets can be downloaded from bitbucket.org/petebunting/spdlib-documentation/

    [1] Zhang, K., Chen, S., Whitman, D., Shyu, M., Yan, J., Zhang, C., 2003. A progressive morphological filter for removing nonground measurements from airborne LIDAR data. IEEE Transactions on Geoscience and Remote Sensing 41 (4), pp. 872–882.
    [2] Evans, J. S., Hudak, A. T., 2007. A multiscale curvature algorithm for classifying discrete return lidar in forested environments. IEEE Transactions on Geoscience and Remote Sensing 45 (4), pp. 1029–1038.