Tag Archives: ARSF

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).

Useful GDAL commands for working with hyperspectral data

The Geospatial Data Abstraction Library (GDAL) contains a number of utilities for working with remote sensing data. The following are some commands which are particularly useful when working with airborne hyperspectral data.

Create subset and change format

To create a band and spatial subset of a larger hyperspectral dataset the following command can be used:

gdal_translate -of GTiff -co COMPRESS=LZW \
 -projwin 481500 766000 482500 765000 \
 -b 190 -b 40 -b 15  \
 outputs/f281013b_utm_wgs84N50_mapped_subsection.bil \
 outputs/f281013b_utm_wgs84N50_mapped_subsection_3band.tif
  • -of: Output Format – we are using GeoTiff file. For all available options see http://www.gdal.org/formats_list.html
  • -co: Creation Option – specifies the creation options, these are specific to the format. For GeoTiff we use compression to reduce the file size. For ENVI we can use the creation options to specify the Interleave (BSQ or BIL)
  • -projwin: A rectangle delineating an area to be extracted, in the format min X min Y max X max Y
  • -b: A given band to extract
  • final two arguments: input and output filenames

Create mosaic of multiple flightlines

Often a mosaic of all individual flightlines is required. Due to the large sizes of the mapped files this can often be time consuming. One option is to first make a ‘Virtial Raster’ using the ‘gdalbuildvrt’ command:

gdalbuildvrt -srcnodata 0 outputs/f2810_mapped_mosaic.vrt \
 outputs/f281013b_utm_wgs84N50_mapped_3band.bil \
 outputs/f281023b_utm_wgs84N50_mapped_3band.bil
  • -srcnodata: No data value of input lines
  • first argument: output VRT file
  • final arguments: list of input files

Rather than creating a new file a small .vrt text file is created which contains links to the individual lines.

Many programs which use GDAL (e.g., TuiView) can open this file directly. To make opening the file easier we can do some optimisation. First pre-calculate statistics using:

gdalinfo -stats outputs/f2810_mapped_mosaic.vrt
  • -stats: Calculate statistics and save to .aux.xml file
  • final argument: file name

The gdalinfo command can also be used to print information about a raster. Note for the .vrt file it shows the individual files which are included.

After creating statistics generate overviews (lower resolution copies at different resolutions) using gdaladdo:

gdaladdo -r average outputs/f2810_mapped_mosaic.vrt 2 4 8 16
  • -r average: Use average for resampling overviews
  • first argument: file name
  • numbers: overview levels to build (2x, 4x, 8x and 16x)

We can also make a real raster using gdal_translate (as shown previously).

Extract pixel values for a given location

gdallocationinfo -valonly -geoloc \
 outputs/f281013b_utm_wgs84N50_mapped_subsection.bil \
 481801 765289
  • -valonly: Only print values to the screen
  • -geoloc: Interprerate x and y values as coordinates rather than pixel and line
  • first argument: file name
  • x y: Easting and Northing of pixel to extract statistics from

This will print the values to the screen. To save as a text file we can redirect the output using ‘>’

gdallocationinfo -valonly -geoloc \
 outputs/f281013b_utm_wgs84N50_mapped_subsection.bil \
 481801 765289 > outputs/extracted_values.txt

This text file can be opened in other programs for further analysis.

This post is based on the hyperspectral practical session from the NERC Airborne Survey and Research Facility workshop held at Plymouth Marine Laboratory, March 2016. For more information about the ARSF Data Analysis Node (ARSF-DAN) and details about future training opportunities follow @ARSF_DAN on twitter. Hyperspectral data acquired by ARSF are available to download from NEODC.

Copy information between ENVI header files

The header files (.hdr) of ENVI files are text files which contain the size of the file (pixels and lines) geospatial information and other meta data such as wavelengths, data units etc.,

Not all programs are aware of all items in then header file and some may become lost during processing. ARSF-DAN have made a script available which can copy attributes from one header to another. The script can be downloaded from https://github.com/pmlrsg/arsf_tools and is run using:

copy_header_info.py f221013_mapped_osng.bil.hdr \
                      f221013_nadir_mapped_osng.hdr \
                      -k "map info" "projection info" "wavelength"

To show the available fields the following script (also from the arsf_tools repo) can be used:

get_info_from_header.py f221013_mapped_osng.bil.hdr

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.

Viewing hyperspectral data in TuiView

As TuiView is based on GDAL, it has always been able to open the ENVI files commonly used for hyperspectral data. However, there have been some recent enhancements to TuiView which have improved handling of hyperspectral data.

TuiView can be installed for Linux, OS X and Windows through conda using:

conda install -c rios tuiview

The source is available to download from https://bitbucket.org/chchrsc/tuiview/.

As hyperspectral data contain many bands, specifying the ones to use for display when you load the file in makes things easier. You can do this from the command line using the ‘–bands’ argument. For example:

tuiview --rgb --stddev --bands 280,460,160 f159223b_mapped_osng.bil

For the Specim Fenix dataset used as an example the command will display a colour composite of NIR, SWIR and Red at approximately the centre wavelengths of the Landsat 8 bands with a standard deviation stretch applied to the pixel values. Opening the image may take a while as the statistics need to be calculated for the the three bands selected for display in order to apply a standard deviation stretch.

Once the file has opened you can use the ‘Query tool’ to display spectral profiles for each pixel.

tuiview_screengrab

If you get an error that TuiView is unable to open the file it could be that ‘data ignore value = 0’ is set in the header and the GDAL driver is unable to handle it. You can export the environmental variable:

export GDAL_PAM_ENABLED=ON

on Linux/OS X, or:

set GDAL_PAM_ENABLED=ON

on Windows. This will create a ‘*.aux.xml’ metadata file with additional features the driver doesn’t support when TuiView opens the file. More information is available in the GDALPamDataset documentation.

If you have hyperspectral data which isn’t mapped (e.g., level1b / level2 data in the NASA Data processing levels) TuiView will not open these by default and will give an error about only allowing north-up images. If you know the data don’t have geospatial information you can force TuiView to open them by setting:

export TUIVIEW_ALLOW_NOGEO=YES

on Linux/OS X, or:

set TUIVIEW_ALLOW_NOGEO=YES

on Windows.

More details about available environmental variables used by TuiView are available on the TuiView Wiki.

Hyperspectral data from the NERC ARSF used for this post were provided courtesy of NERC via the NERC Earth Observation Data Centre (NEODC).

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