Segment individual trees from TLS point clouds with 3DForest

Identifying individual trees from point clouds scanned using Terrestrial LiDAR System (TLS can be easily segmented using this freely downloable software 3DForest. For more information about 3DForest, visit http://www.3dforest.eu/.

Note: although 3DForest has its own ground classification filter, I found it better to perform the ground / non-ground classification using alternative software like PDAL.

To use 3DForest to segment individual trees, the following steps can be done:

1. Start up 3DForest. Select Project | New Project to create a new project.
   
   The 3DForest application appears.
   
2. In the menu, choose Project | Import | Import Terrain Cloud to load in a LAS file of only ground points.
   
   
3. In the dialog box that pops up, choose a ground class only LAS file, e.g. sample_ground.las.
   
   The LAS file is converted into a PCD format file and displayed.
   
4. Next, choose Project | Import Vegetation Cloud. Choose a non-ground class only LAS file, e.g. sample_vegetation.las.
   
   The non-ground class LAS file is converted into PCD format and displayed.
   
5. In the menu, select Vegetation | Automatic Segmentation to run the individual tree segmentation algorithm.
   
   The Automatic Segmentation dialog box appears.
   
   
6. In the Input Vegetation Cloud combo box, choose the newly loaded sample_vegetation LAS file.
   
7. In the Input Terrain Cloud combo box, choose the loaded sample_ground LAS file.
   
8. If necessary, change the parameters, e.g. Prefix of Clouds value from ID to TreeID. Click OK.
   
   The individual trees are segmented into individual point clouds.
   
   
    
9. If required, select Project | Export | Export Clouds to export the segmented trees to text or PLY formatted files.

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.
   

Use CMake to build Visual Studio C++ projects with PDAL on Windows

The tutorial on the PDAL website https://pdal.io/development/writing.html provides the details on developing and compiling a simple C++ program using the PDAL API on *nux. I followed the steps but encountered errors trying to compile the example on Windows 10. The following screenshot shows the compilation errors - the Microsoft Visual Studio linker could not find some external symbols such as __imp_htonl_ among others.


 
To fix the compilation errors, the Windows library ws2_32.lib needs to be linked to the final executable. This can be done in the following steps:

1. Add the library ws2_32 to the CMakeLists.txt file, as shown below.
   
   
2. Run CMake to build the Visual Studio project files.
   
   C:\> cmake -G"Visual Studio 15 2017 Win64" ..
   
   
   The project files are generated.
   
   
3. Open up the solution project in Visual Studio. Then build the solution.
   
   The project is successfully built and the executable is generated.
   

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.
   
   

Using PDAL to perform simple translation transformations on LAS files

Sometimes for whatever reason you may want to apply a simple linear translation (dx, dy, dz) transformation shift on a LiDAR las file. PDAL has a transformation filter that you can use to perform the shifting, but you have to pass in the translation amounts as part of a 4x4 homogeneous transformation matrix.

A 4x4 transformation matrix that performs translations has the following form:
[1 0 0 dx]
[0 1 0 dy]
[0 0 1 dz]
[0 0 0 1]

So to use PDAL to perform a simple translation, just replace the dx, dy, dz with the translation amounts and pass the matrix to the PDAL transformation filter.

The following example command performs a vertical shift of +5.0 to the input las file.
C:> pdal translate -i input.las -o output.laz -f transformation   --filters.transformation.matrix="1 0 0 0 0 1 0 0 0 0 1 5.0 0 0 0 1"  --writers.las.compression="true" -v 4


The following screen shot shows a profile of the input and output las files with the output las file visibly shifted 5.0 above the input las file.

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.
   

Generating LiDAR flight line boundaries with PDAL's tindex

A common task after an airborne LiDAR data collection is to determine the boundaries of the LiDAR LAS file strips, to visualize where they cover approximately. The screenshot below shows an example of how these LiDAR strips look like, colored by elevation. 

Example of LiDAR data strips

The boundaries can be generated in various formats using PDAL's tindex command. Using this command is simple, as shown in the example steps below:

1. Open up the OSGeo4W Command Prompt. Change the directory to the folder containing the *.laz files.
   
   

   

   A folder containing LiDAR *.las file strips

   
   
2. At the prompt, type in the command:
   
   C:\> pdal tindex boundary.sqlite *.laz -f SQLite --lyr_name "bnd"
   
   where boundary.sqlite is the output file
   *.laz specifies the input files using a wildcard
   -f SQLite is the output format
   --lyr_name "bnd" specifies the layer name to append to the output name
3. Press RETURN.
   
   The boundaries are generated.
   
   
4. Optional. After the process is completed, drag and drop the output file e.g. boundary.sqlite onto QGIS.
   
   The LiDAR strips boundary_bnd_polygon are displayed.
   
   
   The output boundaries have the following attributes as shown below.
   

PDAL: Colorize a LAS file with multiple GeoTiff images

A LiDAR LAS file may not necessarily share the same bounds as a single raster image file. More often the case, to cover the LAS file more than one raster image file is needed. One way to resolve this is to merge all the raster images into a single mosaic file. Another way may be to keep the raster images as they are and load them in one by one (just-in-time fashion) to colorize the LAS file. This workflow is a little more complicated but it can be done using PDAL's pipeline processing mechanism.

In this post, PDAL will be used to colorize a LAS file using multiple tiles of GeoTIFF image files - four to be exact. You just need to create a pipeline that loads in the LAS file, apply the RGB values from each image file individually and write the colored points into an output LAS file. The following sections illustrate this workflow.

Here's how the source LAS and GeoTIFF files look like.

The source LiDAR LAS file displayed by elevation

The source GeoTIFF raster image files

Create a PDAL pipeline JSON file

1. Using a text editor, type in the JSON syntax to colorize a LAS file using PDAL's colorization filter, as shown below.
   
   
   where
   uncompahgre.laz is the source LAS file,
   tile0.tif, tile1.tif, tile2.tif and tile3.tif are the source image files,
   and color.laz is the output LAS file name
   
2. Save the text into a file e.g. colorLas.json.
   

A sample pipeline JSON is show in the listing below.

{
    "pipeline": [
        "uncompahgre.laz",
        {
            "type": "filters.colorization",
            "raster": "tile0.tif"
        },
        {
            "type": "filters.colorization",
            "raster": "tile1.tif"
        },
        {
            "type": "filters.colorization",
            "raster": "tile2.tif"
        },
        {
            "type": "filters.colorization",
            "raster": "tile3.tif"
        },
        {
            "type": "writers.las",
            "compression": "true",
            "minor_version": "2",
            "dataformat_id": "3",
            "filename":"color.laz"
        }
    ]
}

Run the colorization process

1. Open up a OSGeo4W Command Prompt.
2. At the prompt, type in the pdal pipeline command:
   
   c:\> pdal pipeline colorLas.json
   where colorLas.json is the pipeline JSON file created in the previous section.
3. Run the command.
   
   Processing messages appear.
   
   
   The output colorized LAS file color.laz is generated.
4. Optional. Display the resultant colored LAS file in a LAS Viewer.
   
   

Use PDAL to generate a DEM from a LiDAR LAS file

PDAL can be used to generate a GeoTIFF digital elevation model (DEM) from a LiDAR point cloud LAS file via the GDAL writer driver. The resultant elevation GeoTIFF file may not be displayed as what is expected in some software such as Global Mapper because the PDAL generated file can have one or more bands of data including an alpha channel, depending on the PDAL options used; refer to https://www.pdal.io/stages/writers.gdal.html for the types of output bands that can be written to the DEM GeoTIFF file. For normal display, it might be necessary to separate the bands into individual files - in this post, the OTB (Orfeo Tool Box) otbcli_splitImage executable is used.

The following steps show how to generate a GeoTIFF DEM and split the resultant file into separate bands.

Generate the DEM

1. Open up a text editor. Type in the following JSON text. Save the text into a file e.g. pdal_dtm.json
   
   
   
   Note: the JSON pipeline text will do the following:
   (a) Load in a LAS file e.g. autzen.laz
   (b) Filter away non-ground classified point
   (c) Write only the mean Z values to a GeoTIFF file dem.tif with a resolution of 10 meters
   
2. Open up the OSGeo4W Command Prompt.
   
   The Command Prompt appears.
   
3. Type in the PDAL command to process the pipeline. Run the command.
   
   C:\> pdal pipeline pdal_dtm.json
   
   The GeoTIFF file dem.tif is generated.
   

Splitting the bands

Displaying the resultant dem.tif containing only the mean Z band in QGIS will show something resembling a grayscale terrain with masked areas, as shown in the screenshot below.

However, displaying the resultant GeoTIFF in Global Mapper will show a black rectangle, as shown below, due to the presence of an alpha channel. This need to be removed for proper interpretation of the elevation data.

To split the bands, do the following:

1. In the OSGeo4W Command Window, type in the command:
   
   C:\> otbcli_splitimage -in dem.tif -out split.tif
   
   
   
   The input file dem.tif is split into separate files with a numbered suffix starting from 0 e.g. split_0.tif.
   
2. Now, the file can be displayed properly in Global Mapper.
   

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.