Export Velodyne PCAP to LAS files in VeloView

VeloView is an open source application for working with Velodyne's popular LiDAR sensors (HDL-64E, HDL-32E, VLP-16 etc.). It is available from this website https://www.paraview.org/veloview/.

VeloView has an export to LAS file format command but you would not know it as it is not available from the graphical user interface. To generate LAS files, it is necessary to use the Python command line interface. The steps below show how it can be done.

1. Run VeloView and load a Velodyne PCAP file.
   
   
2. Choose Tools | Python Console.
   
   The Python Console appears.
   
3. To export out the current frame, type in the following command at the Python prompt.
   
   >>> vv.saveLASCurrentFrame('c:/path/to/output.las', 0)
   
   Note 1: the last argument indicates the coordinate transformation mode. 0 = raw, 1 = absolute, 2 = relative.
   Note 2: Replace c:/path/to/ with your actual file system path
   
4. To export a range of frames by frame numbers, type in the following command at the Python prompt.
   
   >>> vv.saveLASFrames('c:/path/to/output.las', 10, 20, 0)
   
   Note: in this example, 10 is the starting frame, 20 is the ending frame, and 0 is the raw transformation mode.
   
5. To export a list of frames by time steps, type in the following command at the Python prompt.
   
   >>> vv.saveLAS('c:/path/to/output/las', timesteps, 0)
   
   Note: timesteps is an array of integers specifying the time steps e.g. [0,1,2,...10] and 0 is the raw transformation mode.
   
   To get the available time steps in the PCAP file, type in the following command at the Python prompt.
   
   >>> timesteps = vv.getCurrentTimesteps()
   
6. Optional. Display the resultant LAS file in a viewer.
   
   

Use PDAL to export colored LAS file to Point Cloud XYZ format text file

PDAL can be used to export out a colored laser LAS file like the one in the screenshot below into a text file in Point Cloud XYZ format, which is just an ASCII text file with each point in a row of X, Y, Z, R, G, B values.

To use PDAL to perform the conversion, follow the following steps.

1. In Windows, open up an OSGeo4W Shell.
   
   
   
2. In the command prompt, change directory to the folder containing the LAS file e.g. D:\Temp\data\.
   
   
3. Use a text editor to create a PDAL pipeline JSON file, e.g. las2xyz.json, with the following contents.
   
   

   {
     "pipeline":[
       {
         "type":"readers.las",
         "filename":"color.las"
       },
       {
         "type":"writers.text",
         "format":"csv",
         "order":"X,Y,Z,RED:0,GREEN:0,BLUE:0",
         "keep_unspecified":"false",
      "quote_header": "false",
      "delimiter": " ",
         "filename":"outputfile.xyz"
       }
     ]
   }
   

   
   where color.las is the name of the input LAS file
   and outputfile.xyz is obviously the name of the output XYZ file.
   
4. In the command prompt, type in the following command to perform the conversion.
   
   D:\> pdal pipeline las2xyz.json
   
5. Press Enter.
   
   The X, Y, Z, R, G, B values of the LAS file are exported out into the text Point Cloud XYZ file.
   

Perform a 2D geo-registration of a LiDAR LAS file to a GeoTiff image with Cloud Compare

Recently I received LiDAR LAS point clouds generated without using accurate GPS positioning. As a result, the data is not placed at the correct geographic location in the world. To rectify this, I used Cloud Compare to so-called geo-register the point clouds to the world using GeoTiff images as reference.

The following illustrates how it was done:

Load the GeoTiff image as reference

1. Start Cloud Compare.
   
   
   
2. Select File | Open.
   
   The Open file(s) dialog box appears.
   
   
3. Browse and select the reference GeoTiff raster file, e.g. autzen.tif. Click Open.
   
   The Result type prompt appears.
   
   
4. Click No to load as a texture map instead of converting to a point cloud.
   
   The Global shift/scale dialog box appears.
   
   
5. Accept the defaults. Click Yes.
   
   The GeoTiff image is loaded.
   

Load the LAS file

1. Select File | Open.
   
   
   
2. Browse and choose the LAS file, e.g. autzen.las. Click Open.
   
   The Open LAS File dialog box appears.
   
3. Click Apply.
   
   The Global shift/scale dialog box appears.
   
   
4. Click Yes to accept the defaults.
   
   The LAS file is loaded.
   

 Setup the display

1. In the vertical tool bar on the left, click the Set Current View Mode icon
   and choose Orthographic projection.
   
   Note: we do not want the perspective projection display to mislead us as to the location of the features
   
2. Optional. Select Edit | Colors | Height Ramp. Click OK in the Gradient color dialog box to display the point cloud with a height color gradient for better visualization.
   
   
   

Perform the rectification

1.  In the DB Tree pane, select the point cloud to rectify, e.g. autzen.las.
   
2. In the top tool bar, click the Translate/Rotate icon .
   
   The Translate/Rotate widget appears in the top right corner.
   
3. In the Rotation field, select Z. Toggle off Tz.
   
   Note: we only want to move horizontally in the 2D plane.
   
4. Look for matching features in the point cloud and the GeoTiff image.
   
   
   
   
5. To translate the point cloud, press down the mouse right button and drag the point cloud to the new location. To rotate the point cloud, press down the mouse left button and rotate the point cloud.
   
   Note 1: To toggle the display of the point cloud, toggle on/off the Visible property in the Properties pane on the left.
   Note 2: To zoom in/out/pane around, pause the Translate/Rotate widget by clicking the widget's Pause icon.
   
6. In the Translate/Rotate widget, click the Tick icon to save the changes.
   
   The point cloud is rectified.

Use PDAL to apply a vertical datum geoid correction to a LAS file

PDAL can be used to apply a vertical datum correction to LAS files, e.g. to convert LiDAR data in ellipsoidal heights to mean sea level using a geoid such as the EGM2008. PDAL uses geoids in gtx format and can be downloaded from http://download.osgeo.org/proj/vdatum/

The following steps show how to perform a vertical datum correction with the reprojection filter of PDAL:

1. Optional. Download an appropriate geoid file e.g. egm08_25.gtx for EGM2008 from http://download.osgeo.org/proj/vdatum/ and place it in a folder e.g. D:\Temp\PDAL\
   
2. Open up the OSGeo4W Shell.
   
3. In the Command Prompt, type in the following command:
   
   D:\> pdal translate -i input.las -o output.laz reprojection --filters.reprojection.in_srs="EPSG:32610+4326" --filters.reprojection.out_srs="+init=EPSG:32610 +geoidgrids=D:/Temp/PDAL/egm08_25.gtx" --writers.las.compression="true" --writers.las.a_srs="EPSG:32610+3855" -v 4
   
   where
   --filters.reprojection.in_srs specifies the source coordinate reference system EPSG:32610 or UTM 10 North  and the source vertical datum of EPSG:4326 which is the WGS84 Ellipsoid
   
   --filters.reprojection.out_srs specifies the destination coordinate reference system EPSG:32610 (the same as the input) and the EGM2008 grid file to apply
   
   --writers.las.a_srs tells PDAL to write the destination coordinate reference system of EPSG:32610 and vertical datum EPSG:3855 (which is the EGM2008 datum) to the output file.
   
   
4. Press RETURN.
   
   Processing messages appear. The input las file is reprojected to the EGM2008 vertical datum.
   
   
5. Optional. Overlay the input and output LAS files in a viewer and observe the vertical offset in a profile as shown below.
   
   

Use PDAL to register a point cloud to control points

Sometimes, some 3D LiDAR point cloud data are not at the correct location or they are not aligned with other point clouds. When this happens, at least 3 control points from the point cloud and the reference need to be identified and a transformation matrix (usually a 4x4 homogeneous matrix) be calculated.

Once the matrix has been calculated, PDAL's transformation filter can then be used to perform the registration or transformation of the point cloud to the reference.

While PDAL can apply the transformation, it does not have a tool for calculating the matrix. This has to be done externally. This post uses a simple WebApp calculator at this site https://dominoc925-pages.appspot.com/webapp/calc_transf3d/default.html to calculate a rigid 4x4 transformation matrix given 2 corresponding sets of 3 or more control points.

Identifying 2 sets of corresponding control points

1. Using a point cloud viewer such as FugroViewer, identify 3 or more source control points from the unregistered point cloud., e.g. p1, p2, p3 and p4.
   
   
   
   
   
2. Using a GPS receiver or from a database, identify the corresponding points from the destination reference system.
   
   

Calculating the 4x4 transformation matrix

1. Open up url https://dominoc925-pages.appspot.com/webapp/calc_transf3d/default.html in an Internet browser.
   
2. In the Source 3D points field, type or paste in the source control point coordinates from the previous section.
   
   
   
3. In the Destination 3D points field, type or paste in the corresponding control points from the destination reference system from the previous section.
   
4. Click Calculate.
   
   The 4x4 transformation matrix and root mean square error is calculated.
   
   
   Note: if the RMSE is large, then there may be some error in the input coordinates.
   
5. Using a text editor, create a PDAL pipeline JSON file which specifies the input, filter, and output. Copy and paste the 4x4 transformation matrix into the matrix field in the JSON file as shown below. Replace the commas with blanks. 

{
  "pipeline":[
    "input.laz",
    {
      "type":"filters.transformation",
      "matrix":"

0.21095277662308629 0.9774781497556848 -0.00594918357528593 23238.318780198064 

-0.9774332258657229 0.21100438075203337 0.01007175640974158 30100.055196603218 

0.011100225616467532 0.0036902647131527543 0.9999315811282344 32.53070599039722 

0 0 0 1"
    },
    {
      "type":"writers.las",
   "compression": "laszip",
      "filename":"output.laz"
    }
  ]
}

Running the transformation

1. In Windows, open up the OSGeo4W command prompt.
   
   The Administrator: OSGeo4W Shell appears.
   
   
   
2. In the Command Prompt, type in the following and press RETURN:
   
   C:> pdal pipeline transform_pipeline.json
   
   where transform_pipeline.json is the JSON file created in the previous section.
   
   The point cloud is transformed.
   

displaz - A small but fast LiDAR las file viewer

I'm happy with this LiDAR las file viewer displaz I stumbled upon while doing some searching on Google. It's free, open-source, no-frills, fast, and relatively easy to use - just perfect to drag and drop some las files for a quick look. It can be downloaded from this web site http://c42f.github.io/displaz/.

The following are some sample screenshots of the viewer.

1. Upon start up, the viewer shows a blank view with a shader parameters and Log or Datasets panes on the right.
   
   
   
   Note: The points are colored according to intensity values. This can be changed by choosing an alternative coloring method in the Shader parameters pane.
   
2. Las files and be opened by simply dragging it(them) and dropping onto the application window from the operating system's file managers.
   
   
   Note: once opened, the file(s) are listed in the Datasets pane. To view the point cloud from a different orientation, drag and move the mouse left button. To zoom, roll the mouse wheel.
   
3. This is a little confusing at first, but to pan or move the view, you have to press the middle button of the mouse to center the view around the chosen point. At the same time, it will show the attributes of the chosen point in the Log pane, as shown below.
   
   
   
4. Finally, to unload a las file, select one or more files in the Datasets pane and press the DELETE button on the keyboard.
   
   
   

Mask a LAS file using PDAL and QGIS

The Point Data Abstraction Library (PDAL) comes with a useful function to mask or crop out a LiDAR LAS file with one or more polygons. The example screenshot below shows a point cloud overlaid with a red polygon in the upper right corner, which outlines the desired area of the point cloud to be cropped.

The cropping can be done using PDAL's crop filter but it requires the cropping polygon to be specified in the Well Known Text (WKT) string format. This is a bit of pain but can be overcome using a few methods, one of which is shown below using QGIS and the Plain Geometry Editor plugin.

Define the cropping polygon

1. In QGIS, draw a new polygon e.g. mask, as shown below.
   
   
   
   Note: The mask should be created in the same coordinate system as the LAS file. 
2. Click the Plain Geometry Editor icon in the toolbar (red circle above). Click on the clipping polygon.
   
   The Plain Geometry Editor dialog box appears.
   
   
   Note: Install the Plain Geometry Editor plugin if the icon is not available.
3. In the text field, select and copy all the polygon WKT text string into the Windows clipboard.  

Create a PDAL processing pipeline JSON file

1. In a text editor, type in something similar to the example below.
   
2. From the Windows Clipboard, paste the WKT string from the previous section to the "polygon" attribute and surround it with double quote '"' characters.
   
3. Save the JSON file e.g. process.json.

{
  "pipeline":[
    "autzen.laz",
    {
      "type":"filters.crop",
      "polygon": "Polygon ((638500.66904077248182148 853359.34703735215589404, 638869.71793351718224585 853365.15883093187585473, 638881.34152069012634456 853208.24040409235749394, 638677.92874516302254051 853199.52271371346432716, 638878.43562390014994889 852818.85023379256017506, 638733.14078423334285617 852842.09740813483949751, 638611.09311891463585198 852975.76866063219495118, 638416.39803376235067844 853211.14630088687408715, 638500.66904077248182148 853359.34703735215589404))"
    },
    {
      "type":"writers.las",
      "filename":"file-cropped.las"
    }
  ]
}

Note: 

• The pipeline JSON file stores the processes to be done in sequence in an array bracketed by the '[' and ']' characters. 
• autzen.laz is the input LAS file for this example
• filters.crop is the process to apply using the "polygon" attribute.
• file-cropped.las is the output LAS file. 

Run the cropping process

1. Open up the OSGeo4W Shell.
   
2. At the prompt, type in the pdal pipeline command:
   
   C:\> pdal pipeline process.json -v 4
   Processing messages appear. The file is cropped.
   
   
3. Optional. Using your preferred LAS Viewer, open up the resultant cropped LAS file.
   
   The cropped file showing only the cropped area is displayed.
   

Using PDAL to classify isolated LiDAR points as noise

LiDAR data often contains noise and it is necessary to identify and/or remove them. An example of a LAS file containing noise in the form of low isolated points beneath the ground is shown in the screen shot below.

Isolated points can be easily identified by using statistical filtering methods, which the PDAL open source software has.

To filter out these points using PDAL, perform the following steps.

1. Open up the OSGeo4W Shell.
   
   The OSGeo4W Shell prompt appears.
2. In the prompt, type in the command:
   
   C:\> pdal translate -i in_noisy.las -o out_filtered.las outlier --filters.outlier.method="statistical" --filters.outlier.mean_k=8 --filters.outlier.multiplier=3.0 -v 4
   
   Note:
   -i in_noisy.las is the input LAS file
   -o out_filtered.las specifies the output LAS file
   outlier tells PDAL to apply the outlier filter
   --filters.outlier.**** options specify the various outlier parameters
   -v 4 indicates the processing messages verbosity level
   
   
3. After running the command, the isolated points are classified as Class 7 - Low Noise points in the output LAS file.
   
   The point cloud colored by classification.
   
   
   The resultant LAS file colored by elevation and with the class 7 - Low noise points turned off.
   

Simple LiDAR ground points classification and segmentation using PDAL

PDAL (Point Data Abstraction Library) comes with a couple of options to segment point clouds by classifying LiDAR ground points (an example unclassified point cloud is shown below) - Simple Morphological Filter (SMRF) or Progressive Morphological Filter (PMF).

I have found the SMRF method to be fast and produce reasonable results while the PMF method seems to take a much longer time to do the job. The steps to run ground classification on a LAS file are describe below.

1. In Windows, open up the OSGeo4W Shell.
   
   The OSGeo4W Shell is displayed.
2. In the OSGeo4W prompt, type in and run the command:
   
   C:\> pdal translate -i unclassified.las -o ground.las smrf -v 4
   
   Notes:
   -i unclassified.las is the input file
   -o ground.las specifies the output file
   smrf is the option to apply the Simple Morphological Filter
   -v 4 is the processing messages verbosity level
   
   Processing messages appear.
   
3. Display the ground classified LAS file in a viewer.
   
   
   
   
4. To use the Progressive Morphological Filter to perform the ground classification, type in the following command:
   
   C:\> pdal translate -i unclassified.las -o ground.las pmf -v 4
   
   Notes:
   -i unclassified.las is the input file
   -o ground.las specifies the output file
   pmf is the option to apply the Progressive Morphological Filter
   -v 4 is the processing messages verbosity level

Assign a coordinate reference system to a LiDAR las file using PDAL

Some of the LiDAR las files I receive do not contain embedded coordinate reference system (CRS) tags. It would be nice to be able to set a CRS tag to the las file so that I don't have to choose a coordinate reference system everytime I load the file; and I found the Point Data Abstraction Library (PDAL) to have the tools to do exactly that.

To use PDAL to assign a CRS to a las file, do the following:

1. Open up a Command prompt. Type in the pdal command:
   
   C:\> pdal translate -i input.las -o output.las --writers.las.a_srs="EPSG:32750"
   
   Note: where a_srs is the option to assign a CRS e.g. EPSG:32750
   
   
2. To double check whether the CRS tag has been assigned to the las file, you can use the lasinfo tool from liblas, as shown below.