Using PDAL to clean legacy LiDAR dataPete Gadomski January 31, 2017 [how-to] #riegl #lidar #snow #pdal
The CRREL/UCSV Energy Site (CUES) is a snow research site on Mammoth Mountain, CA:
In February 2011 a Riegl LMS-Z390 LiDAR scanner came online at the site, mounted in a glass enclosure on the central structure.
This scanner regularly captures snow depths over a portion of the study area, and stores these snow depths on the CUES web server as
While these data have been used with some research, in general they have been under-utilized, partially due to their file format.
.3dd is an old Riegl file format, so old that the company no longer produces binaries to read from the format — the only binaries are available in
riscanlib for Windows 32 systems.
And because Riegl doesn't release the layout of their file formats to the public, we can't1 write our own multi-platform reader.
riscanlib, which requires a Riegl website account to download, comes with an example executable named
conv2asc.exe that reads a
.3dd file and prints the X, Y, Z, and Intensity coordinates to standard output.
For now, this is what we're using to extract information from
However, there's a couple of hangups:
- Our trusty ol' Z390 scanner produces a lot of bogus data at
(0., 0., 0.)that should be filtered out.
.3ddintensities are floats in [0, 1], but other point cloud formats (e.g. las) work best when intensity is scaled from [0, 65536].
- Our scanner shifts and moves through time, both intentionally when someone makes and adjustment and naturally due to shifts in the ground. We need to apply a rigid transform to our data to ensure all of our scans are in the same reference frame.
We'll tackle parts 1 and 2 in this post, and save part 3 for a later date.
If you want to follow along at home with the same data, download it from the CUES website.
We'll assume you've already used
conv2asc.exe on the
3dd file and redirected the output to
PDAL (with the Python plugin enabled) is the only required software for this exercise (assuming you have a unix-y system that includes
Filtering and scaling
PDAL's text reader expects reasonably well-formed text files.
conv2asc.exe output looks like this, which isn't good enough:
$ head txt/original/20161229-1315-30.Z390.frame.txt 0.0000, 0.0000, 0.0000, 0.0000 0.0000, 0.0000, 0.0000, 0.0000 0.0000, 0.0000, 0.0000, 0.0000 0.0000, 0.0000, 0.0000, 0.0000 0.0000, 0.0000, 0.0000, 0.0000 0.0000, 0.0000, 0.0000, 0.0000 0.0000, 0.0000, 0.0000, 0.0000 0.0000, 0.0000, 0.0000, 0.0000 0.0000, 0.0000, 0.0000, 0.0000 0.0000, 0.0000, 0.0000, 0.0000
We need to remove those spaces between the comma and the next number, and add a header. This is a simple two-step process:
mkdir -p txt/cleaned echo "X,Y,Z,Intensity" > txt/cleaned/20161229-1315-30.Z390.frame.txt sed 's/, /,/g' txt/original/20161229-1315-30.Z390.frame.txt >> txt/cleaned/20161229-1315-30.Z390.frame.txt
We'll use PDAL to do the
(0., 0., 0.) filtering and intensity scaling.
The filtering can be done out of the box with range filter, but we need to write some code to do the intensity scaling.
We'll use the programmable filter to do the scaling work.
Create a file called
scripts/scale_intensity.py with the following code:
= = * 65536 = return True
This takes our
[0, 1] intensity values and scales them to
[0, 65536], which happens to be the maximum range of intensity values supported by the las format.2
We then can created a "cleaned" laz file with no origin points3 and scaled intensity values with this one-liner:
By combining simple external tools (e.g.
sed) with PDAL's powerful filters, we can transform uncomfortable data formats into more usable layouts.
I've gisted a Makefile and the scale intensity script that are useful for batch processing:
In a later post I'll discuss how to use transformation matrices to bring all of these data into the same coordinate system.
Inquiring minds might wonder if we're corrupting useful data by scaling these values all willy-nilly. By convention, when LiDAR systems report "intensity" values, they are sensor-specific values that have no absolute relationship to a real-world radiometric property of the reflective object. Riegl does attempt to provide an absolute measure of an object's reflectivity in the laser's wavelength with the "Reflectivity" attribute, but even this is fraught with issues.
As of this writing, the negation operator ("!") is an undocumented feature of the range filter.
Back to top