CZO Lidar Workshop 2012
Lidar (Light Detection and Ranging) is a remote sensing technology that measures distance to an object (and other properties) using laser pulses; the time it takes a laser pulse to reflect from an object and return to a sensor is used to calculate the distance to that object. A single pulse can have multiple returns (up to four in most cases) as it penetrates different levels of an object (such as a tree). The first level to return is called the first return and can be considered the top most layer of an object. The last level to be returned is called the last-return and used to help identify points considered to be the ground. Several pulses can be collected to form a point cloud and used to create useful products that explain the surface being recorded.
Lidar point cloud from forest
Lidar data is usually recorded to a binary format (for efficiency) with the extension “.las” known as an LAS file. The first part of an LAS file is called the header and contains information about how to read the data (e.g. how many points the file contains, point offset and scaling, etc.).
LAS file header
The second part of the LAS file is where individual point entries are stored and their attributes. Some attributes include (besides x, y and z) classification (e.g. ground, building, water, etc), intensity, number of return and other attributes.
LAS file (individual point entry) - Version 2
All programs and Lidar data for these exercises are contained in the “C:\workshop2012” directory on the lab computers. The raw Lidar data being used is located in “C:\workshop2012\data\SS_LAS”. All output should go into the “C:\workshop2012\data\output” directory. The programs are located in “C:\workshop2012\programs” unless otherwise stated. Alternatively, you can download all this from the “Resources” section.
For the CZO, Lidar was delivered in several tiles that make up each area. This is because it is much easier to deal with the smaller file sizes of the tiles rather than one large file for an entire area. If a particular study area is only a small extent of the entire area, it may be okay to deal with one tile; the study area may be extracted from the one tile. However, it is most likely that a study area will spread across several tiles and would require extracting the study area from each individual tile, than merging them.
Several tools exist to extract subsets from each LAS tile that overlap with the desired study area. For this exercise, we will use tools from LAStools and the scripting language Python. LAStools are free (for non-commercial) tools written in C++ for efficiently processing lidar. Python allows scripting of these tools for batch programming.
The steps required for this process include:
A script has already been written to accomplish this called “ClipLAS.py”. Open “ClipLAS.py” (located in the “programs\clip_lastiles” directory) in a text editor if you disire to see how it utilizes LAStools. “clip_lastiles_GUI.py” is the graphical user interface version of “ClipLAS.py”.
There are 2 ways to run ClipLAS: through the command line and through a graphical user interface (GUI). For this exercise, we will be working with the GUI.
Note: All Exercises build on each other!
Open the clip lastiles GUI by entering the “programs” directory, and double clicking “clip_lastiles_GUI”. Point the input directory to where the LAS files are saved. Give an output (including filename) to the output field (call it “studyarea.las”). Fill out the extent parameters using the following values.
Your study area will have the following extents:
Click OK. Once it completes, close the program. Check the output directory for your new LAS file.
Saga GIS is an open source program with many GIS tools including importing and viewing Lidar files. Saga can easily import, manipulate and visualize LAS files.
Open Saga from the programs directory (or start menu if you installed manually). Many tool boxes filled with modules are available in the Workspace window (usually on the left side of the program) under the modules tab. Scroll down to the toolbox named “Import/Export - LAS” and expand it. Double click the “Import LAS Files” module.
This opens a dialog where you can specify the file as well as the different attributes to import (like classification, return number, etc). Click on the “Input File” field, then click on the ellipse that appears on the right to navigate to the LAS file you created from exercise 1. Under “Attributes to import...”, check intensity and classification.
Click OK. The Message window toward the bottom will indicate that you are running the import module and the status bar at the bottom will display the current progress for the module. Allow it finish.
In the workspace window in Saga, click the “Data” tab. Notice that the data you imported is listed under the PointCloud category. Double click your data, you will be presented with a 2D view of your data.
In the workspace window in Saga, click the “Modules” tab. Scroll to the “Shapes - Point Clouds Viewer” Toolbox and expand it. Double click “Point Cloud Viewer”. In the dialog, choose your data for “>> Points”.
The Point Cloud Viewer will start and you will see a three dimensional version of your data. Using the mouse, hold the left mouse button down on your data and move the it around - you can view all angles of the data.
Use the Extent box in the lower left corner to view smaller parts of the data.
Change the Color Attribute drop menu to intensity - Notice how the view changed. Now change the Color Attribute to classification. Experiment with different angles with the different color attributes.
Boundary polygons are useful for representing lidar data. For this we will use another program from LAStools called “lasboundary”. “lasboundary” creates a shapefile representing the extents of an las file. It can be used from the command prompt by issuing the following command:
lasboundary -i input.las -o output.shp
Enter the programs directory. Double click on “Boundary Command Prompt”. This will open the command prompt and add an environment variable to allow you to use “lasboundary.exe” from any directory. First change to the directory that contains your data by typing the following command (your path may vary depending on where you output your file):
To create the boundary, type the following:
lasboundary -i studyarea.las -o studyarea.shp
Close the command prompt. Your output directory should contain “studyarea.shp” which is the boundary of the study area. To view this in Saga, from within Saga select File -> Shapes -> Load.
Navigate to the “studyarea.shp” file and click open. From the workspace window, select the “Data” tab and double click studyarea under the Shapes category.
MCC-lidar is an open source software that classifies lidar data as ground (indicated by a 2) and non-ground using the Multiscale Curvature Classification algorithm (see Acknowledgements section for more information) useful for forest applications.
This program will be ran from the command line and requires 2 parameters in addition to the input file and output file. The parameters are scale (‘s’) and the curvature threshold (‘t’). The scale is a function of the size of objects (e.g., trees, rocks) on the ground and the resolution of the Lidar data. The threshold is a function of the expected slope of an area.
Since areas can be heterogeneous (e.g., trees of differing sizes, downed trees, rocks, etc.), selecting values for the parameters may become trial and error. So, it is recommended that you start with values of 0.9 for scale and increment it up or down by 0.5 and 0.3 for curvature and increment it up or down by 0.1.
The following is an example of how to use MCC-LIDAR:
mcc-lidar -s 0.9 -t 0.3 input_filename.las output_filename.las
See the Acknowledgements section and visit the MCC-LIDAR website for more information or type “mcc-lidar --help” for more options.
Open MCC-LIDAR Command Prompt from the programs directory (or start menu if you installed manually) - This opens the command prompt and temporarily adds mcc-lidar as a system environment variable allowing you to execute it from any directory. Change directory to where your output from exercise 1 is saved - type the following in the command prompt (your path may differ):
Now, execute mcc-lidar using the data from exercise one as the input file and the recommended starting parameters:
Type in the following into the command prompt (with the appropriate file names):
mcc-lidar -s 0.9 -t 0.3 studyarea.las studyarea_reclass.las
Import the new classified lidar file into Saga GIS (be sure to include the classification attribute) and view it in the point cloud viewer. Change the “Color Attribute” to classification. Notice the distinction between the trees and the ground.
Note: This process may take a while. Since some points have already been classified in the original dataset, please continue with the rest of the exercises while this process continues in the background. This exercise was to demonstrate how to classify ground points in a data set where they have not yet been classified or if you need to identify more ground points in your dataset or study area.
There may be a need to extract points from Lidar that satisfy some criteria. For example, when looking at the intensity of Lidar points, you may find that paved road points give a much lower intensity than other points in a lidar data set (this is because the pavement absorbs more of the laser pulse). To extract points for the paved road, you could determine a threshold on the intensity that distinguishes the road points from the other points using intensity. Than simply extract all points whose intensity fall below the threshold.
Another example might be in creating a digital elevation model (DEM) where you would need to extract only ground points from a dataset.
In Saga, in the workspace window make sure you are on the modules tab, scroll to the “Shapes - Point Clouds” toolbox and expand it. Double click the module “Point Cloud Reclassifier / Subset Extractor”.
You will be presented with the subset extractor tool (which is also the reclassifier tool and can be used to reclass points) where you can extract a subset of points based on some attribute and criteria. In this case we want to extract ground points (which are classified as 2 under classification).
For the point cloud, select “studyarea.las” which was imported earlier. Select classification as the attribute. Leave result as “[create]” to create a new point cloud with the extracted subset. For “Mode of operation” select “Extract Subset” and leave “Method” as “single”. Set “old value” to 2 (points matching the value we wish to extract) and the operator to “=”.
Press OK and allow the module to finish. View the result in the point cloud viewer - Notice how it appears to be only the ground points (no vegetation).
Sometimes it may be useful to know the point density when producing the Lidar products. This gives us an idea of the quality of the data which could aid in determining things like grid resolution for Lidar products.
To determine the quality of the ground points, from within the workspace window under the modules tab, expand the “Shapes - Point Cloud” toolbox and double click the “Point cloud to Grid” module.
In the dialog box that appears, select “studyarea_subset_classification” for “>>points”. Leave “only z” for “Output” and change the “Aggregation” to “mean value”. Set “Cellsize” to 1 (meter).
Press OK and allow the tool to finish. When it is complete you will have 2 grids : one representing the mean Z value of points that fall in each grid cell and the other representing how many points fall in each grid cell. In the workspace window, change to the “data” tab, you will see the 2 grids under the “Grids” category.
At this point, we still need to fill in the “No Data” cells with a value of 0 for correct statistics on the grid, otherwise the statistics will only be performed on cells that contain data. In the workspace window under the modules tab, scroll up to the “Grid - Tools” toolbox and expand it. Double click the “Create Constant Grid” module.
In the window that appears, for “Grid System”, select the grid system that contains the 2 grids and for the “Base grid” select the “Study area Subset Classification [Points Per Cell]” grid. For “value”, use 0 for the constant value for all grid cells.
Click OK. This produces a grid that is the same size and resolution as the point density grid whose value for all cells is 0. Now, double click the “Merging” module in the same toolbox, which allows you to merge grids together. For the “Grids to Merge” field, click on the input field then click on the ellipse that appears all the way to the right.
In this window, add the points per cell grid and the constant grid you created (order does matter, so add them in that order).
Click OK. For Grid system select the same grid system your point density grid is contained in. Select “[Create]” for the target grid. For “Interpolation” and “Overlapping cells” select “Nearest Neighbor” and “first value in order of grid list” respectively.
Click Ok. This will create a grid called “Merged Grid” in which the no data cells from the point density grid will be filled with a 0.
From the data tab in the workspace window, select the “Merged Grid” grid by clicking on it once. On the right side of Saga is the properties window. At the bottom, select the “Description” tab. Scroll all the way down and notice the Arithmetic Mean (should be about 1.5 points per cell. Also notice the Standard Deviation and other grid statistics).
Generally it is good to have at least 1 point per grid cell when interpolating a grid (more is better). For this exercise, we can assume from this that a grid product with a cell resolution of 1 is acceptable.
Here, the digital elevation model (DEM) is a bare earth surface where vegetation and other objects (buildings, rocks, etc) have been removed such that you have a continuous surface representing a bare landscape. In this sense, it is also called the digital terrain model (DTM). When using Lidar to produce the DEM, vegetation points are removed which leaves holes or no data under the vegetation. These areas often have to be interpolated to create a continuous surface.
First, let us create our own target grid where we specify the size and resolution. This is useful if you have existing grids and you want to match their exact size and resolution. In the workspace window, expand the “Grid -Tools” module. Double click the “Create Grid System” Module. See the image below for parameters.
(xMax and yMax will be ignored because of the option we chose for “Set extent by”). Click OK. This will produce your grid system and dummy grid that you will output your target grid to.
For this exercise, we will use the triangulation interpolation method to create the continuous DEM surface. In the workspace window, expand the “Grid - Gridding” toolbox and double click the “Triangulation” module. For “>>Points” select “studyarea_subset_classification” and for “Attribute” select “Z”. Select “grid” for “Target grid”.
Click OK. After a few seconds, a “Choose grid” dialog with appear. Select the grid system you produced earlier as the “Grid system” and for “<< Grid” select “[create]”.
Click OK. Wait for the module to finish, then double click the new grid it to view it (Should be called “studyarea_subset_classification (Triangulation)”).
Optional: For further visual analysis on the new DEM, create a hillshade by using the “Analytical Hillshade” tool found in the “Terrain analysis - Lighting, Visibility” toolbox. Load up the triangulation DEM and accept default parameters. View the new hillshade grid. You may find a different interpolation method yields better DEMs.
The digital surface model (DSM) is a continuous grid of the top most surface including all objects across a landscape. This is most commonly produced by accepting the points with the highest Z value in each grid cell. With Lidar, the top most points are usually the most abundant (because they are easier to attain from above considering that penetrating the vegetation is unnecessary), and thus can usually be produced with no interpolation; the point density is usually high enough such that interpolation is unnecessary.
In the “Grid - Gridding” toolbox, double click the “Shapes to Grid” Module. For “>> Shapes”, select “studyarea”, which contains all points (including ground and vegetation). Select “Z” for the “Attribute” parameter. For “Method for Multiple Values” and “Target Grid” select “maximum” and “grid” respectively.
Click OK. In the “Choose Grid” dialog, select the grid system that contains your DEM and “[create]” for the “<< Grid” parameter.
Click OK. Double click the new grid to look at it.
The next exercise demonstrates how a Lidar product can be used.
The canopy height model (CHM) is the normalized height model where each pixel represents the real height for each cell (eliminates elevation). This is simply the DEM subtracted from the DSM.
Expand the “Grid - Calculus” toolbox and select the “Grid Difference” module. Select the grid system that contains your DEM and DSM. For “>> A” select your DSM (studyarea [Z]) and for “>> B” select your DEM (studyarea_subset_classification (Triangulation)).
Click OK. Double click the new grid to view it (Notice everywhere near the ground is consistently colored due to the normalization).
Next we will filter the CHM using a Gaussian filter to smooth out noise to easier identify trees.
Expand the “Grid - Filter” toolbox and double click the “Gaussian filter”. Select the grid system that contains the CHM then select the CHM for the “>> Grid” parameter. For “< Filtered Grid” select “[create]”. Set “Search Radius” to 5 (this should roughly be the width of the trees).
Expand the “Imagery - Segmentation” toolbox and double click the “Watershed Segmentation” module. Choose the grid system that contains the filtered grid. For “>> Grid” choose the filtered grid you just produced. Select “[create]” for “< Segmentation” and “Seed Value” for “Output”.
For this exercise, we only want to identify larger trees above a certain value such that the follow equation is satisfied:
where a is the segmentation grid you produced and g is the output grid.
Open the grid calculator by expanding the “Grid Calculus” toolbox and double clicking on “Grid Calculator”. Select the grid system that contains the segmented grid you just created. Click the “Grids >>” input field then click the ellipse to the right. In the dialog that appears, add the segmented grid.
Click Ok. Select “[create]” for “<< Result”. In the “Formula” field, type the following equation:
where “ifelse” is a boolean statement (where 1 is true and 0 is false) and “gt” is a greater than function (that returns 1 if a condition is true and 0 otherwise)
Now the the trees are separated, the final step is vectorize the tree segments to points.
Expand the “Shapes - Grid” toolbox and double click the “Vectorizing Grid Classes” module. Choose the grid system from the threshold calculation you just performed and select that calculation grid for the “>> Grid” parameter.
Click OK. You will now have a vector layer with polygons that represent the segmented trees. Now, expand the “Shapes - Polygons” toolbox and select the “Polygon Centroids” module. Select the polygon layer you just created for the “>> Polygon” parameter. For “<< Centroids” select “[create]”.
Click OK. Compare the CHM with the segmented tree positions by double clicking the “Different (A-B)” layer and selecting a new map. Then double click the tree location layer “Calculation [ifelse(gt(a,5),a,-99999)]” under points, and add to the map you just created (the one just above new).
Click OK to see how they compare.
For tree number estimate, in the workspace window under the data tab, click once on the tree location layer “Calculation [ifelse(gt(a,5),a,-99999)]” under points. In the properties window on the right side of Saga, click the description tab and notice the value for “Number Of Shapes”. This is the estimated number of trees.
Saga can export data to common formats to be used in other software such as ArcGIS.
Save the DEM layer by expanding the “Import/Export - GDAL/ORG” toolbox and double clicking the “GDAL: Export Raster” Module. Select the grid system that contains the DEM. Click on the “>> Grid(s)” input field and click the ellipse that appears to the right. In the dialog that appears, add the DEM grid.
Click Ok. Now click the “File” input field and click the ellipse that appears. Name and Save the file somewhere with the extension “.img”, e.g. DEM.img. For format, choose “Erdas Imagine Images (.img)”.
Click Ok. The grid was saved in the image format (*.img). You can repeat the same process to save the DSM and CHM.