Those are my notes on taking some LiDAR data and testing different commands in PDAL. If you like, you can follow along using your own data.
Most probably I’ll de adding or deleting stuff, but you can take a look to see some of the main operations that can be made on point clouds using PDAL.
This is an adaptation of the original workshop from the PDAL website, Point Cloud Processing and Analysis with PDAL, the single difference being that I use the same data for every operation and show the commands output, too, while in the workshop there are used different data sources. Also, because the files were large, I had to make tiles and use batch processing for some operations.
You can use Fugro Viewer or LASTools to view LiDAR files on desktop, but only in Windows! The Linux users, although forgotten, are lucky to have Potree and Plasio that can consume Entwine Point Tile (EPT) datasets in the browser. But let’s just approach LiDAR with baby steps before digging too deep.
LiDAR refers to a remote sensing technology that emits intense, focused beams of light and measures the time it takes for the reflections to be detected by the sensor. This information is used to compute ranges, or distances, to objects. In this manner, lidar is analogous to radar (radio detecting and ranging), except that it is based on discrete pulses of laser light. The three-dimensional coordinates (e.g., x,y,z or latitude, longitude, and elevation) of the target objects are computed from 1) the time difference between the laser pulse being emitted and returned, 2) the angle at which the pulse was “fired,” and 3) the absolute location of the sensor on or above the surface of the Earth.
You can use your own files, but I’ll be using free data for Yosemite Valley.
There are three point clouds of LAS type for Yosemite Valley on OpenTopography. In order to conduct this study in the area, the files were downloaded, merged and transformed into raster files. One file was too large, it was splitted in two for download, resulting four files for processing. We will convert the files to LAZ format, too, in order to save space.
The original files and the job ids used to download and split the data are:
Yosemite National Park, CA: Rockfall Studies. Distributed by OpenTopography. https://doi.org/10.5069/G9D798B8. Accessed: 2019-12-10
We create a working folder, yose and a subfolder named raw_data, then move all the downloaded point cloud data into it.
# go home, create subdirectory and parent directory cd&&mkdir -pv yose/raw_data
Install PDAL
You can read this other post to find out how to install PDAL. Every time you want to access the environment where PDAL was installed, useconda activate followed by the environment name, in our case new_clouds.
I will refer to Linux, though the commands for Windows will be mentioned, too. NOTE: PDAL commands, except batch processing, are the same in Windows and Linux, but the line separator differs: ^ in Windows and \ in Linux.
# List the LAZ filescd ~/yose
ls -lah raw_data
# Query the las data using pdal info
pdal info raw_data/points_little_yose.las -p 0
pdal info raw_data/points_extra_yose.las -p 0
pdal info raw_data/points_rockfall_1.las -p 0
Create a new folder named temp_data, decompress the file from LAZ to LAS or vice-versa using pdal translate, then compare their size. Note that it would take some time to complete the transformation, if the files were large, but after the compression the storage decreases significantly.
mkdir temp_data && pdal translate raw_data/points_little_yose.las temp_data/points_little_yose.laz
pdal translate raw_data/points_extra_yose.las temp_data/points_extra_yose.laz
pdal translate raw_data/points_rockfall_1.las temp_data/points_rockfall_1.laz
pdal translate raw_data/points_rockfall_2.las temp_data/points_rockfall_2.laz
ls -lah raw_data/points_little_yose.*
ls -lah raw_data/points_extra_yose.*
ls -lah raw_data/points_rockfall_1.*
ls -lah temp_data/points_rockfall_2.*
(new_clouds) myuser@comp:~/yose$ pdal translate raw_data/points_little_yose.las raw_data/points_little_yose.laz
(new_clouds) myuser@comp:~/yose$ pdal translate raw_data/points_extra_yose.las raw_data/points_extra_yose.laz
(new_clouds) myuser@comp:~/yose$ pdal translate raw_data/points_rockfall_1.las raw_data/points_rockfall_1.laz
(new_clouds) myuser@comp:~/yose$ pdal translate raw_data/points_rockfall_2.las raw_data/points_rockfall_2.laz
(new_clouds) myuser@comp:~/yose$ ls -lah raw_data/points_little_yose.*
-rwxrwxrwx 1 myuser myuser 800M dec 2022:33 raw_data/points_little_yose.las
-rw-r--r-- 1 myuser myuser 78M feb 1622:41 raw_data/points_little_yose.laz
(new_clouds) myuser@comp:~/yose$ ls -lah raw_data/points_extra_yose.*
-rwxrwxrwx 1 myuser myuser 4,6G dec 2011:25 raw_data/points_extra_yose.las
-rw-r--r-- 1 myuser myuser 643M feb 1622:43 raw_data/points_extra_yose.laz
(new_clouds) myuser@comp:~/yose$ ls -lah raw_data/points_rockfall_1.*
-rwxrwxrwx 1 myuser myuser 3,6G nov 2917:54 raw_data/points_rockfall_1.las
-rw-r--r-- 1 myuser myuser 613M feb 1622:44 raw_data/points_rockfall_1.laz
(new_clouds) myuser@comp:~/yose$ ls -lah raw_data/points_rockfall_2.*
-rwxrwxrwx 1 myuser myuser 4,4G nov 2919:28 raw_data/points_rockfall_2.las
-rw-r--r-- 1 myuser myuser 708M feb 1622:46 raw_data/points_rockfall_2.laz
(new_clouds) myuser@comp:~/yose$
From now on we will use the compressed files in our operations.
We will use Entwine to create a new folder named ept_yose, then rearrange points into a lossless octree structure known as EPT.
You can read this other post to find out how to install Entwine.
The build command is used to generate an Entwine Point Tile (EPT) dataset from point cloud data. Because we want to build the visualization from all of our data, for input we provide the path to raw_data folder the command. Note that we select only the LAZ files.
Now we can statically serve the entwine directory with an HTTP server and visualize it with the WebGL-based Potree and Plasio projects.
Potree and Plasio
We can visualize our data in Potree and Plasio, directly in the browser.
# Install NodeJS and http server from nodejs
conda install nodejs -y
npminstall http-server -g
# Make sure you are in the working folder, and the EPT folder# is a subfolder# Start the server using --cors to allow serving data to the# remote Potree/Plasio pages
http-server -p 8066 --cors
(new_clouds) myuser@comp:~/yose$ http-server -p 8066 --cors
Starting up http-server, serving ./
Available on:
http://127.0.0.1:8066
http://192.168.0.102:8066
http://172.18.0.1:8066
http://192.168.61.1:8066
http://192.168.63.1:8066
Hit CTRL-C to stop the server
# In another tab (Shift+Ctrl+T), visualise out local point cloud data in Potreexdg-open"http://potree.entwine.io/data/view.html?r=http://localhost:8066/ept_yose"
Assign a colour scale to the points based on elevation, classification and other criteria, then play with some terrain profiles.
first view
elevation
el capitan profile
inner profile
the dome profile
inner profile
elevation
ground only
ground & unclassified
about
Visualize the points by elevation:
top view
top zoom
front view
left view
back view
right view
the dome
the dome
the dome
the dome
the dome
Glacier Point apron
camp 4
the sentinel
el capitan top
el capitan & three bro.
el capitan & cathedral sp
cathedral rocks & sp.
el capitan
el capitan & bridalveil
The commonly used versions of the LAS format (1.2 and 1.3) have 8 classification categories pre-defined and can handle up to 32. For 1.4 see this document.
Visualize the points by point classification:
ground
low vegetation
overlap
ground & unclassified
LiDAR data can contain significantly more information than x, y, and z values and may include, among others, the intensity of the returns, the point classification (if done), the number of returns, the time, and the source (flight line) of each point.
Intensity can be used as a substitute for aerial imagery when none is available.
Visualize the points by other attributes (intensity data, return number):
level of detail
elevation
index
intensity
intensity gradient
return number
source
Visualize the points in Plasio, too.
# Visualise our local point cloud data in Plasioxdg-open"http://dev.speck.ly/?s=0&r=ept://localhost:8066/ept_yose&c0s=local://color"
Or use Hexbin to output the boundary of actual points in point buffer, not just rectangular extents. Create a file named hexbin-pipeline.json and populate it with the code below, then run pdal pipeline.
Before running the next commands, we will have to split the files into tiles. Then run the command in batches, because we can run out of memory if the file was too large, with the error PDAL: std::bad_allocfor example.
We will start from the beginning, and use the LAS files from the raw_data folder. First, we reproject points_little_yose.las from “EPSG:32611” to “EPSG:26911” (from “WGS 84 / UTM zone 11N” to “NAD83 / UTM zone 11N”).
WGS 84 / UTM zone 11N - EPSG:32611 - EPSG.io
epsg.io › 32611
Traducerea acestei pagini
EPSG:32611 Projected coordinate system for Between 120°W and 114°W, northern hemisphere between equator and 84°N, onshore and offshore.
NAD83 / UTM zone 11N - EPSG:26911 - EPSG.io
epsg.io › 26911
Traducerea acestei pagini
EPSG:26911 Projected coordinate system for North America - between 120°W and 114°W - onshore and offshore.
Then we move the files we need for making tiles in a new folder, las_files.
# List the LAS files in folderls -la raw_data/*.las
# Create the folder where we will gather our datamkdir -pv las_files
# Move the datamv raw_data/points_little_yose_nad83.las raw_data/points_extra_yose.las raw_data/points_rockfall_1.las raw_data/points_rockfall_2.las las_files
# List the files in the newly created folderls -la las_files
(new_clouds) myuser@comp:~/yose$ ls -la raw_data/*.las
-rwxrwxrwx 1 myuser myuser 4858578967 dec 2011:25 raw_data/points_extra_yose.las
-rwxrwxrwx 1 myuser myuser 838430789 dec 2022:33 raw_data/points_little_yose.las
-rw-r--r-- 1 myuser myuser 838430703 feb 1920:58 raw_data/points_little_yose_nad83.las
-rwxrwxrwx 1 myuser myuser 3853695669 nov 2917:54 raw_data/points_rockfall_1.las
-rwxrwxrwx 1 myuser myuser 4627696169 nov 2919:28 raw_data/points_rockfall_2.las
(new_clouds) myuser@comp:~/yose$ mkdir -pv las_files
mkdir: created directory 'las_files'(new_clouds) myuser@comp:~/yose$ mv raw_data/points_little_yose_nad83.las raw_data/points_extra_yose.las raw_data/points_rockfall_1.las raw_data/points_rockfall_2.las las_files
(new_clouds) myuser@comp:~/yose$ ls -la las_files
total 13846136
drwxr-xr-x 2 myuser myuser 4096 feb 1921:14 .
drwxr-xr-x 7 myuser myuser 4096 feb 1921:13 ..
-rwxrwxrwx 1 myuser myuser 4858578967 dec 2011:25 points_extra_yose.las
-rw-r--r-- 1 myuser myuser 838430703 feb 1920:58 points_little_yose_nad83.las
-rwxrwxrwx 1 myuser myuser 3853695669 nov 2917:54 points_rockfall_1.las
-rwxrwxrwx 1 myuser myuser 4627696169 nov 2919:28 points_rockfall_2.las
(new_clouds) myuser@comp:~/yose$
And now we create LAZ tiles using the pdal tilecommand. If we don’t specify --length, each output file contains points in the 1000×1000 square units represented by the tile.
# Create the folder for tilesmkdir -pv laz_tiles
# Create the tiles from all files in the specified folder
pdal tile -i "las_files/*.las"\
-o "laz_tiles/yose_#.laz"\
--buffer="20"\
--out_srs="EPSG:26911"# List the files in the newly created folderls -la laz_tiles |head -10
ls -la laz_tiles |tail -10
(new_clouds) myuser@comp:~/yose$ mkdir -pv laz_tiles
mkdir: created directory 'laz_tiles'(new_clouds) myuser@comp:~/yose$ pdal tile -i "las_files/*.las"\> -o "laz_tiles/yose_#.laz"\> --out_srs="EPSG:26911"(new_clouds) myuser@comp:~/yose$ ls -la laz_tiles |head -10
total 2095836
drwxr-xr-x 2 myuser myuser 4096 feb 1921:23 .
drwxr-xr-x 7 myuser myuser 4096 feb 1921:16 ..
-rw-r--r-- 1 myuser myuser 426303 feb 1921:24 yose_0_0.laz
-rw-r--r-- 1 myuser myuser 6070318 feb 1921:24 yose_0_-1.laz
-rw-r--r-- 1 myuser myuser 9484799 feb 1921:24 yose_0_-2.laz
-rw-r--r-- 1 myuser myuser 8178030 feb 1921:24 yose_0_-3.laz
-rw-r--r-- 1 myuser myuser 7145497 feb 1921:24 yose_0_-4.laz
-rw-r--r-- 1 myuser myuser 6998980 feb 1921:24 yose_0_-5.laz
-rw-r--r-- 1 myuser myuser 1387969 feb 1921:24 yose_0_-6.laz
(new_clouds) myuser@comp:~/yose$ ls -la laz_tiles |tail -10
-rw-r--r-- 1 myuser myuser 46019612 feb 1921:24 yose_9_0.laz
-rw-r--r-- 1 myuser myuser 27848492 feb 1921:24 yose_9_-1.laz
-rw-r--r-- 1 myuser myuser 39658296 feb 1921:24 yose_9_1.laz
-rw-r--r-- 1 myuser myuser 33373584 feb 1921:24 yose_9_-2.laz
-rw-r--r-- 1 myuser myuser 29441795 feb 1921:24 yose_9_2.laz
-rw-r--r-- 1 myuser myuser 13122580 feb 1921:24 yose_9_-3.laz
-rw-r--r-- 1 myuser myuser 9873921 feb 1921:24 yose_9_3.laz
-rw-r--r-- 1 myuser myuser 1817222 feb 1921:24 yose_9_-4.laz
-rw-r--r-- 1 myuser myuser 27836 feb 1921:24 yose_9_4.laz
-rw-r--r-- 1 myuser myuser 167052 feb 1921:24 yose_9_-5.laz
(new_clouds) myuser@comp:~/yose$
To aquire a more accurate representation of the ground returns, you can remove the non-ground data to just view the ground data, then use the translate command to stack the filters.outlier and filters.smrf.
Now do the same for all the files in laz_tiles, providing the input file name and the output file name on the fly in batch processing. Parallel, as its name implies, allows paralell operations. After that, merge all the resulting tiles in a new LAZ file named denoised_grnd_f.laz.
(new_clouds) myuser@comp:~/yose$ ls laz_tiles/*.laz |\> parallel -I{} pdal translate {}\> -o denoised_grnd/denoised_grnd{/.}.laz \> outlier smrf range \> --filters.outlier.method="statistical"\> --filters.outlier.mean_k=8 --filters.outlier.multiplier=3.0\> --filters.smrf.ignore="Classification[7:7]"\> --filters.range.limits="Classification[2:2]"\> --writers.las.compression=true \> --verbose 4
Academic tradition requires you to cite works you base your article on.
When using programs that use GNU Parallel to process data for publication
please cite:
O. Tange (2011): GNU Parallel - The Command-Line Power Tool,
;login: The USENIX Magazine, February 2011:42-47.
This helps funding further development; AND IT WON'T COST YOU A CENT.
If you pay 10000 EUR you should feel free to use GNU Parallel without citing.
To silence this citation notice: run 'parallel --citation'.
(PDAL Debug) Debugging...
(pdal translate Debug) Executing pipeline in standard mode.
(pdal translate filters.smrf Debug) progressiveFilter: radius =140177 ground 995 non-ground (2.42%)(pdal translate filters.smrf Debug) progressiveFilter: radius =133026 ground 8146 non-ground (19.79%)(pdal translate filters.smrf Debug) progressiveFilter: radius =231082 ground 10090 non-ground(24.51%)(pdal translate filters.smrf Debug) progressiveFilter: radius =330493 ground 10679 non-ground(25.94%)(pdal translate filters.smrf Debug) progressiveFilter: radius =430315 ground 10857 non-ground(26.37%)(pdal translate filters.smrf Debug) progressiveFilter: radius =530276 ground 10896 non-ground(26.46%)(pdal translate filters.smrf Debug) progressiveFilter: radius =630265 ground 10907 non-ground(26.49%)[...](pdal translate filters.smrf Debug) progressiveFilter: radius =16843723 ground 237877 non-ground(21.99%)(pdal translate filters.smrf Debug) progressiveFilter: radius =17843470 ground 238130 non-ground(22.02%)(pdal translate filters.smrf Debug) progressiveFilter: radius =18843232 ground 238368 non-ground(22.04%)(pdal translate writers.las Debug) Wrote 6703951 points to the LAS file(new_clouds) myuser@comp:~/yose$
Generate an elevation model surface using the output from above. Create a file named dtm.json and populate it with the code below, then run pdal pipeline. The GDAL writer creates a raster from a point cloud using an interpolation algorithm. (src)
Generate a raster surface using a fully classified point cloud. Create a file named classification.json and populate it with the below code, then run pdal pipeline.
(new_clouds) myuser@comp:~/yose$ gdal_translate temp_data/intensity.tif temp_data/intensity.png -of PNG
Input file size is 22747, 9618
Warning 6: PNG driver doesn't support data type Float32. Only eight bit (Byte) and sixteen bit (UInt16) bands supported. Defaulting to Byte
0...10...20...30...40...50...60...70...80...90...100 - done.
(new_clouds) myuser@comp:~/yose$
# import numpyimport numpy as np
# import matplotlib stuff and make sure to use the# AGG renderer.import matplotlib
matplotlib.use('Agg')import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
# This only works for Python 3. Use# StringIO for Python 2.from io import BytesIO
# The make_plot function will do all of our work. The# filters.programmable filter expects a function name in the# module that has at least two arguments -- "ins" which# are numpy arrays for each dimension, and the "outs" which# the script can alter/set/adjust to have them updated for# further processing.defmake_plot(ins, outs):# figure position and row will increment
figure_position =1
row =1
fig = plt.figure(figure_position, figsize=(6,8.5), dpi=300)for key in ins:
dimension = ins[key]
ax = fig.add_subplot(len(ins.keys()),1, row)# histogram the current dimension with 30 bins
n, bins, patches = ax.hist( dimension,30,
facecolor='grey',
alpha=0.75,
align='mid',
histtype='stepfilled',
linewidth=None)# Set plot particulars
ax.set_ylabel(key, size=10, rotation='horizontal')
ax.get_xaxis().set_visible(False)
ax.set_yticklabels('')
ax.set_yticks((),)
ax.set_xlim(min(dimension),max(dimension))
ax.set_ylim(min(n),max(n))# increment plot position
row = row +1
figure_position = figure_position +1# We will save the PNG bytes to a BytesIO instance# and the nwrite that to a file.
output = BytesIO()
plt.savefig(output,format="PNG")# a module global variable, called 'pdalargs' is available# to filters.programmable and filters.predicate modules that contains# a dictionary of arguments that can be explicitly passed into# the module by the user. We passed in a filename arg in our `pdal pipeline` callif'filename'in pdalargs:
filename = pdalargs['filename']else:
filename ='histogram.png'# open up the filename and write out the# bytes of the PNG stored in the BytesIO instance
o =open(filename,'wb')
o.write(output.getvalue())
o.close()# filters.programmable scripts need to# return True to tell the filter it was successful.returnTrue
(new_clouds) myuser@comp:~/yose$ pdal pipeline histogram.json
anything:46: UserWarning: Attempting to set identical left == right ==0 results in singular transformations; automatically expanding.
(new_clouds) myuser@comp:~/yose$
Example syntax. There is another example in Identifying ground section.
# Create directory unless existsmkdir -pv gtiles
# List files in Windows dos and iterate them# use verbose to view logsfor %f in(*.*)do pdal pipeline denoise_ground.json ^
--readers.las.filename=%f ^
--writers.las.filename=../gtiles/%f
--verbose 4# Create directory unless existsmkdir -pv denoised_grnd
# List files in Linux shell and iterate themls laz_tiles/*.laz |\
parallel -I{} pdal pipeline denoise_ground.json \
--readers.las.filename={}\
--writers.las.filename=denoised_grnd/denoised_grnd{/.}.laz \
--verbose 4
Repeated the steps for Yosemite Valley
What’s next? Repeat the steps for an urban area Then try LasTools, too And, of course, the 3D in Blender.
PS: You can admire each dataset directly in Potree using Generate 3D point cloud browser visualization on OpenTopography corresponding page. Reminder: The colored point cloud data that we have used throughout this post can be viewed at https://potree.entwine.io/data/view.html?r=”https://maptheclouds.com/ept_yose_colored/”. Below, I leave you a Potree mix of all the Yosemite Valley, from my personal computer.
Yosemite Valley
The Dome
El Capitan
Published by MapTheClouds
Hi all! Here’s Cristina, founder @MapTheClouds.
I’m an independent cartographer from România, based in Cluj-Napoca.
I all 3D and maps.
Feel free to check out my projects at maptheclouds.com, and reach out if you need a custom map.
View all posts by MapTheClouds
Repeated the steps for Yosemite Valley 
Repeat the steps for an urban area 
all 3D and maps.
Feel free to check out my projects at maptheclouds.com, and reach out if you need a custom map.