Thinning a LAS file with FME 2011

FME Desktop 2011 comes with a bunch of revamped transformers designed to work with point clouds - specifically LiDAR datasets. In previous versions, FME can only read in LiDAR LAS files but I'm thrilled with this new version now since it comes with a LAS file writer.

I tried out the new point cloud transformers by performing one of the standard LiDAR processing production tasks - thinning the data by resampling at every nth point. It is simple to build a transformation workspace to do the job. Here are the steps:

1. Start up FME Workbench 2011.
   
   
2. Click Blank Workspace.
   
   
3. From the Windows Explorer, drag a LAS file and drop it onto the FME Workbench Main pane.
   The Add Reader dialog box appears.
   
4.  Click OK.
5. Activate the Transformer Gallery tab on the left.
6. Find and select the PointCloudThinner transformer in the Transformer Gallery tree list. Drag and drop it onto the Main pane.
7. Select Writers | Add Writer.
   
   The Add Writer dialog box appears.
8. Choose the ASPRS Lidar Data Exchange Format (LAS). Browse and select the output dataset.
9. Click OK. Click Yes and OK to the next few prompts.
10. Connect the datasets and transformer as shown.
11. Open up the PointCloudThinner Parameters.
    
    The PointCloudThinner Parameters dialog box appears.
12. Choose the Keep Every Nth Point Thinning Type. In the Amount field, type in a value e.g. 10.
    
    
13. Click OK.
    
    
14. Run the translation.
    
    The LAS file is thinned
    

Split an area polygon with gvSIG

gvSIG's data editing functions are among the best I have ever encountered in a free and open source GIS software application. One of the useful command is the Split command. Here is how it works.

1. Start up gvSIG OADE 2010. Add in a polygon shape file layer, e.g. polygons.shp.
   
   The layer is displayed in the View.
   
2. In the legend, click the polygon layer.
   
   The layer is selected.
3. Select Layer | Start Editing.
   
   The legend name is highlighted in red and the command pane appears.
   
4. Select Geometry | Select. Click on the polygon.
   
   The polygon fill color changes to red and the vertex handles appear.
   
5. Select Geometry | Modify | Split.
   
   The prompt Insert first point appears in the Command pane.
6. Click a few points to digitize the splitting line.
   
   
7. To complete, type in e in the Command pane prompt. Press RETURN.
   
   The polygon is split.
   
8. Select Layer | Stop editing.
   
   You are prompted to save the edits.
9. Click Yes.
   
   The changes are saved.

Removing holes from a polygon with gvSIG

I wanted to simplify a polygon of a LiDAR flight coverage strip by removing all the holes (inner polygons in gvSIG terminology) within the polygon. I could not find a ready function other than the Delete Vertex command. Using this would be very tedious for hundreds of holes. So I thought of using gvSIG's Geoprocessing Tools to do the job. The idea is (1) to decompose the polygon and its holes into lines, (2) to remove all the short lines except the longest line and optionally to simplify the line work, and (3)  to make a polygon from the remaining line.

Decompose the Polygon into Lines

1. Start up gvSIG OADE 2010. Add a polygon with holes layer e.g. polygonholes.shp to the view.
   
   
2. Select View | Geoprocessing Tools.
   
   The Geoprocessing tools dialog box appears.
   
   
3. Expand the Geoprocessing tools node and select Topology | Reduce to lines.
4. Click Open Tool.
   
   The Analysis tools dialog box appears.
   
5. In the Input layer field, choose the polygon with holes layer e.g. polygonholes.shp.
6. Click Choose and type in the name of the Output layer e.g. lines.shp.
7. Click OK.
   
   The polygon and holes are reduced to lines in the output layer lines.shp.
   

Find and remove short lines

1. In the legend (or table of contents in gvSIG terminology) pane, click the decomposed lines layer, e.g. lines.shp.
   
   The selected layer becomes the active layer.
2. Select Layer | Start Editing.
   
   
   
   The layer name changes to bold red in the legend, and the console pane appears.
   
3. Select Layer | Add Geometry Info.
   
   
   The Add geometry info dialog box appears.
   
4. In the Select geometry information group box, select Length. Click the > button. Click OK.
   
   The length attribute field is created for the lines.shp layer.
   
5. Select Layer | Show attribute table.
   
   The Table: Attribute table for lines.shp dialog box appears.
   
6. Click the Length column header.
7. Select Table | Sort descending.
   
   The table is sorted with the longest length in the first row.
   
8. Select the second row. Scroll to the last row. Press SHIFT and click the last row.
   
   All the rows between the second row to the last row are selected.
   
9. Press Delete on the keyboard.
   
   All the selected records are deleted.
   
10. Close the table. Select Layer | Stop editing.
    
    The Save prompt message appears.
    
11. Click Yes.
    
    The changes are saved.

Simplify the line work

1. Select View | Geoprocessing tools.
   
   The Geoprocessing tools dialog box appears.
   
2. Expand the Geoprocessing tools node. Select Computational Data conversion | Generalize.
3. Click Open tool.
   
   The Analysis tools dialog box appears.
   
4. In the Input layer field, select lines.shp.
5. Toggle on Douglas-Peucker method.
6. In the Distance tolerance field, type in a value e.g. 20.
7. In the Output layer field, type in the path and name of the output file, e.g. c:\temp\generalize.shp.
8. Click OK.
   
   The generalized line layer is created.
   

Convert the line work to a polygon

1. Select View | Geoprocessing tools.
   
   The Geoprocessing tools dialog box appears.
2. Expand the Geoprocessing tools node. Select Topology | Build polygons.
   
   
3. Click Open tool.
   
   The Analysis tools dialog box appears.
   
4. In the Input layer field, choose generalize.shp.
5. In the Output layer field, choose or type in the path and name of the output layer e.g. C:\Temp\polygons.shp.
6. Click OK.
   
   The polygon layer without holes is created.
   

Tip for developing a liblas C++ program with Microsoft Visual Studio 2008

Recently I tried to write a C++ program for reading a LiDAR LAS file using Microsoft Visual Studio and the open source library liblas. I followed the sample tutorial from the liblas.org web site. The simple code below simply opens up a LAS file and creates a LASReader class. Compiling the simple code is okay but when I ran it, the program crashed when trying to create the LASReader class. What was the problem?

 // liblas  
 #include <liblas/liblas.hpp>  
 #include <liblas/laspoint.hpp>  
 #include <liblas/lasreader.hpp>  
 #include <liblas/cstdint.hpp>  
 #include <fstream>  
 #include <iostream>  
 #include <string>  
 int main(int argc, char* argv[])  
 {  
   std::ifstream ifs;  
   ifs.open( "c:\\temp\\input.las", std::ios::in | std::ios::binary);   
   liblas::LASReader reader(ifs);  
 }  

The source code is simple enough and the syntax is correct. After some research I found out that the liblas.lib file I was linking was a released (non-debug) version while I was compiling a debug version of my C++ executable. Due to licensing issues, libLAS can only be distributed in released version. If you want to compile a debug executable, then you have to download the liblas source code and compile your own debug liblas.lib. Otherwise, you can only build release copies of your executable.