Rasterizing Perfect Canopy Height Models from LiDAR

In literature you sometimes read “we generated a Canopy Height Models (CHM) and then did this and that” without the process that was used to create the CHM being described in detail. One approach computes the CHM as a difference between DSM and DTM: create a DTM from the ground returns and a DSM from the first returns and subtract the two rasters. Also here is still a question left to be answered: how exactly are the DTM and the DSM generated. A different approach computes the CHM directly from height-normalized LiDAR points. And again there are many ways of doing so and we want to look at the possibilities in more detail.

the 100 by 100 meter sample plot 'drawno.laz'

the 100 by 100 meter sample plot ‘drawno.laz’

In the following we demonstrate different alternatives for CHM generation on a 100 by 100 meter sample LiDAR tile ‘drawno.laz‘ from a forest near Drawno in Poland (that you can download here), slowly converging towards the CHM generation method that we recommend using. We start with ground classifying the LiDAR using lasground:

lasground -i drawno.laz ^
          -wilderness ^
          -o ground.laz

Then we height-normalize the LiDAR using lasheight. As we know that there are no trees higher than 28 meters in this plot we drop all LiDAR points that are higher than 30 meters that may be bird hits or other noise.

lasheight -i ground.laz ^
          -drop_above 30 ^
          -replace_z ^
          -o normalized.laz

As the sample plot has an average pulse spacing of around 0.3 meters we decide to use a step size of 0.33333 meters to create a 300 by 300 pixel raster. We produce a false color visualization instead of a height raster for our CHMs so we can include the results here. The simplest method uses lasgrid with option ‘-highest’ that uses for each pixel the highest z coordinate among all LiDAR returns falling into the corresponding 0.33333 by 0.33333 meter area.

lasgrid -i normalized.laz ^
        -step 0.33333 ^
        -highest ^
        -false -set_min_max 0 25 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o chm_grd.png
Gridding the highest point that falls into each 0.33333 meter by 0.33333 meter cell.

Gridding the highest point that falls into each 0.33333 meter by 0.33333 meter cell.

The resulting CHM (shown at a 200 % zoom) is full of empty pixels and so called “pits” that will hamper subsequent analysis for single tree detection, height and crown diameter computation, and the like. A simple improvement can be obtained by replacing each LiDAR return with a small disk. After all, the laser beam has – depending on the flying height – a diameter of 10 to 50 centimeter and approximating this with a single point of zero area seems overly conservative. In lasgrid there is the option ‘-subcircle 0.1’, which replaces each return by a circle with a radius of 10 centimeter or a diameter of 20 centimeter.

lasgrid -i normalized.laz ^
        -subcircle 0.1 ^
        -step 0.33333 ^
        -highest ^
        -false -set_min_max 0 25 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o chm_grd_d20.png
Gridding the highest z value after turning each point into a circle with 20 cm diameter.

Gridding the highest z value after turning each point into a circle with 20 cm diameter.

The resulting CHM (shown again at a 200 % zoom) is much improved but there are still empty pixels and “pits”. We could simply widen the circles further with ‘-subcircle 0.15’, ‘-subcircle 0.2’, or ‘-subcircle 0.25’. As you can see below this produces increasingly smooth CHMs with widening tree crowns. But this “splats” the LiDAR returns into circles that are growing larger and larger than the laser beam diameter and thus have less and less in common with reality.

Gridding the highest returns will often leave empty pixels in the data even when “splatting” the points. Another popular approach avoids this by interpolating all first returns with a triangulated irregular network (TIN) and then rasterizing it onto a grid to create the CHM. This can be implemented with las2dem as shown below:

las2dem -i normalized.laz ^
        -first_only ^
        -step 0.33333 ^
        -false -set_min_max 0 25 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o chm_tin.png
Rasterizing the TIN that interpolates all first returns onto a  0.33333 meter grid.

Rasterizing the TIN that interpolates all first returns onto a 0.33333 meter grid.

The result has no more empty pixels but is full of pits because many laser pulses manage to deeply penetrate the canopy before producing the first return. When combining multiple flight lines some laser pulses may have an unobstructed view of the ground under the canopy without hitting any branches. These “pits” and how to avoid them is discussed in great length in the September 2014 edition of the ASPRS PE&RS journal by a paper of Khosravipour et al. We build upon these ideas in the following. But first we combine the ‘-highest’ gridding with TIN interpolation with a two step approach: (1) keep only one highest return per grid cell with lasthin (2) interpolate all these highest returns with las2dem.

lasthin -i normalized.laz ^
        -step 0.33333 ^
        -highest ^
        -o temp.laz
las2dem -i temp.laz ^
        -step 0.33333 ^
        -false -set_min_max 0 25 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o chm_tin_his.png
Rasterizing the TIN that interpolates only the highest points falling into each  0.33333 meter by 0.33333 meter grid cell.

Rasterizing the TIN that interpolates only the highest points falling into each 0.33333 meter by 0.33333 meter grid cell.

Next we integrate the idea of “splatting” the points into circles with a diameter of 20 centimeter to account for the laser beam diameter by adding option ‘-subcircle 0.1’ to lasthin.

lasthin -i normalized.laz ^
        -subcircle 0.1 ^
        -step 0.33333 ^
        -highest ^
        -o temp.laz
las2dem -i temp.laz ^
        -step 0.33333 ^
        -false -set_min_max 0 25 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o chm_tin_his_d20.png
Rasterizing the TIN that interpolates only the highest points of a 0.33333 meter grid  after first splatting points into circles with 20 cm in diameter.

Rasterizing the TIN that interpolates only the highest points of a 0.33333 meter grid after first splatting points into circles with 20 cm in diameter.

The results are much nicer but there are still pits. However, one may argue at this points that thinning with a step size of 0.33333 is too agressive and removes too many points for subsequent interpolation and rasterization with the exact same step size. So we double the resolution of the temporary point cloud by thinning with half the step size of 0.16667.

lasthin -i normalized.laz ^
        -subcircle 0.1 ^
        -step 0.16667 ^
        -highest ^
        -o temp.laz
las2dem -i temp.laz ^
        -step 0.33333 ^
        -false -set_min_max 0 25 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o chm_tin_hhs_d20.png
Rasterizing the TIN that interpolates only the highest points of a 0.16667 meter grid  after first splatting points into circles with 20 cm in diameter onto a raster with step size 0.33333.

Rasterizing the TIN that interpolates only the highest points of a 0.16667 meter grid after first splatting points into circles with 20 cm in diameter onto a raster with step size 0.33333.

We now have more detail but also many more pits. Furthermore the interpolation of highest returns makes a big error across areas were we do not have any LiDAR returns that are flanked by canopy returns on both sides. This happens in the area circled where wrong canopy is created because the TIN is interpolating canopy returns across a small wet area without any returns, We show now how the pit-free method of Khosravipour et al. can be used to generate the perfect CHM:

rmdir tmp_chm /s /q
mkdir tmp_chm
las2dem -i normalized.laz ^
        -drop_z_above 0.1 ^
        -step 0.33333 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o tmp_chm/chm_ground.bil
lasthin -i normalized.laz ^
        -subcircle 0.1 ^
        -step 0.16667 ^
        -highest ^
        -o temp.laz
las2dem -i temp.laz ^
        -step 0.33333 -kill 1.0 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o tmp_chm/chm_00.bil
las2dem -i temp.laz ^
        -drop_z_below 2 ^
        -step 0.33333 -kill 1.0 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o tmp_chm/chm_02.bil
las2dem -i temp.laz ^
        -drop_z_below 5 ^
        -step 0.33333 -kill 1.0 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o tmp_chm/chm_05.bil
las2dem -i temp.laz ^
        -drop_z_below 10 ^
        -step 0.33333 -kill 1.0 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o tmp_chm/chm_10.bil
las2dem -i temp.laz ^
        -drop_z_below 15 ^
        -step 0.33333 -kill 1.0 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o tmp_chm/chm_15.bil
las2dem -i temp.laz ^
        -drop_z_below 20 ^
        -step 0.33333 -kill 1.0 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o tmp_chm/chm_20.bil
las2dem -i temp.laz ^
        -drop_z_below 25 ^
        -step 0.33333 -kill 1.0 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o tmp_chm/chm_25.bil
lasgrid -i tmp_chm/chm*.bil -merged ^
        -step 0.33333 ^
        -highest ^
        -false -set_min_max 0 25 ^
        -ll 278200 602200 -ncols 300 -nrows 300 ^
        -o chm_pit_free_d20.png
rmdir tmp_chm /s /q
Running the pit-free algorithm on the highest LiDAR returns in a 0.16667 meter grid (after splatting them to circles 20 cm in diameter) and producing a 0.33333 meter raster CHM.

Running the pit-free algorithm on the highest LiDAR returns in a 0.16667 meter grid (after splatting them to circles 20 cm in diameter) and producing a 0.33333 meter raster CHM.

With such a perfectly pit-free output we can now be even more conservative and lower the diameter of the circles that we replace each LiDAR return with from 20 cm to 10 cm by replacing ‘-subcircle 0.1’ with ‘-subcircle 0.05’ in the above script.

Running the pit-free algorithm on the highest LiDAR returns in a 0.16667 meter grid (after splatting them to circles only 10 cm in diameter) and producing a 0.33333 meter raster CHM.

Running the pit-free algorithm on the highest LiDAR returns in a 0.16667 meter grid (after splatting them to circles only 10 cm in diameter) and producing a 0.33333 meter raster CHM.

Compared to the original pit-free algorithm published by Khosravipour et al. there are two minor differences: (1) instead of using all first returns as input we use one highest returns per grid cell after splatting all returns as small circles (instead of area-less points) onto a grid that has twice the output resolution. (2) we also have a ‘-kill 1.0’ threshold to also generate a partial CHM from all points for ‘chm_00.bil’ and add a new ‘chm_ground.bil’ to fill the potential holes. This prevents that higher up canopy returns are wrongly connected across water bodies where there are no LiDAR returns at all.

Reference:
Khosravipour, A., Skidmore, A.K., Isenburg, M., Wang, T.J., Hussin, Y.A., 2014. Generating pit-free Canopy Height Models from Airborne LiDAR. PE&RS = Photogrammetric Engineering and Remote Sensing 80, 863-872.

Tutorial: editing LAS or LAZ files “by hand” with lasview

This tutorial describes how to manually edit LiDAR using the new inspection and editing functionality available in ‘lasview.exe’ with the latest release of LAStools (version 140301). We will work with the familiar ‘fusa.laz’ sample LiDAR data set from the LAStools distribution that was recently reported to have shown strange symptoms assumed to be side-effects of the LAZ cloning experiments in the ESRI labs … (-;

Inspecting LiDAR files with cross sections

Copy ‘fusa.laz’ from the folder ‘lastools\data’ to the folder ‘lastools\bin’. Run ‘lasview.exe’ so that it loads ‘fusa.laz’. Either do this via the GUI by double-clicking ‘lasview.exe’, loading ‘fusa.laz’ via the ‘browse …’ menu, and then clicking the ‘VIEW’ button or by entering the command below:

C:\lastools\bin>lasview -i fusa.laz

Press the <x> key to toggle to “select cross” where you can pick a rectangular “cross” section. The default cross section is a profile extending across the bottom of the bounding box.

tutorial4_lasview_01_select_cross

By pressing the <x> key again you toggle back to actually view the cross section. Holding down <ALT> you can rotate the view to look at the cross section from the side. Holding down <CTRL> you can zoom in and out. Holding down <SHIFT> you can translate up and down or left and right. Increase or decrease the size of the points pressing <=> or <->. Hover with the mouse over a point and press <i> to inspect its coordinates and attributes.

tutorial4_lasview_02_cross_inspect_point

Traverse the LiDAR file visually by moving the cross section with the arrow keys <UP> <DOWN> <LEFT> and <RIGHT>. You can move either in the “select cross” view and see the picked rectangle move or in the “cross” view and “walk” through the LiDAR. Hold down the <SHIFT> key simultaneously to take bigger steps or the <ALT> key to take smaller steps. Inspect other points by hovering over them with the cursor and pressing <i>. The point information disappears when pressing <i> with the cursor over the background.

tutorial4_lasview_03_traverse_with_arrow_keysToggle back to the “select cross” view with <x> and pick approximately the same rectangle as shown below:

tutorial4_lasview_04_pick_missclassified_roof

Changing Classifications and Deleting Points

Continuing the steps above, toggle back to the “cross” view by pressing <x>. Note that part of the roof of the house has been miss-classified as vegetation while others are left unclassified. Press <e> to turn on the “EDIT” mode and right-click to select “reclassify points as building (6)” via the pop-up menu.

tutorial4_lasview_05_reclassify_menu

Now use the cross-hair cursor to draw a polygonal fence around all points that should be reclassified. Press <ESC> to remove the last vertex of the polygon if you miss-placed it by a mistake.

tutorial4_lasview_06_reclassify_fenceOnce you are happy with your polygon press <r> to register the edit. A note appears informing you how many points had their classification changed. In the top right corner an “undo” counter appears informing you how many changes you can undo by pressing <CTRL-u>. Try it. Immediataly the changes disappear and a “redo” counter appears instead. Press <CTRL-o> to redo the change you have just undone.

tutorial4_lasview_07_edit_result

Press <CTRL-s> to save this edit as a tiny LAY file using the recently introduced LASlayers concept. In case there was already an existing LAY file (that was not applied with ‘-ilay’ when starting lasview) you will be warned and have to press <CTRL-f> to force overwriting it as shown below.

tutorial4_lasview_07a_first_save

Press the key sequence <SHIFT-b>, <t>, and <a> to get the same visuals above.

tutorial4_lasview_07a_second_save

Press <SHIFT-t> to remove the triangulation again. After saving an edit it can no longer be undone via <CTRL-u>. Instead you will have to strip off this particular layer with the layer management available through “laslayers.exe” as described here. Now press <x> to toggle to the “cross select” view.

tutorial4_lasview_08_move_select_cross

Use the <DOWN> arrow to move the selected cross section to the area shown above that has a few unclassified points in the middle of the roof. Press <x> to go back to the “cross” view and try to understand why these points are not part of the roof. Looks like they are from the top of a chimney, and antenna, or a satellite dish as they do not fit the otherwise planar roof.

Assume we need to remove them for some reason. Pan, translate, and zoom the view such that these points can be easily surrounded by a polygon. Now press <d> to enter the “DELETE” mode, fence in these points, and press <r> to register the deletion.

tutorial4_lasview_09_delete_points

It can be tricky to place a clearly seperating polygon and you may be worried about deleting a few orange building points as well. Press <u> to only display the unclassified points before pressing <r> to register the deletion.

tutorial4_lasview_10_delete_points_unclassified_only

Press <a> to see all points again, then delete the other two points by finding a good view point, pressing <d>, and drawing a polygon.

tutorial4_lasview_11_delete_other_points

After registering this deletion of two points your “undo” counter should be at two. Press <CTRL-u> twice to undo this and the last deletion, then press <CTRL-o> twice to redo them both.

tutorial4_lasview_12_delete_other_points_undo_count_2

Now press <CTRL-s> to save this deletion as another layer. It will be appended to the LAY file that already contains one layer with the roof re-classification edit we did first. Press <t> to triangulate the points in the “cross” view. See how nicly flat the triangulated roofs are now that we deleted these 6 chimney points.

tutorial4_lasview_13_second_save

Look at the size of the tiny LAY file called ‘fusa.lay’ that is in the same folder as the ‘fusa.laz’ file. It contains all the edits we have done so far and mine is only 681 bytes in size. The original LAZ file has not changed. Maybe this is all you want for now. You could send only this tiny LAY file to a colleague elsewhere and he or she could apply those changes locally when needed using the ‘-ilay’ switches. For more on this see the LASlayers page.

C:\lastools\bin>lasview -i fusa.laz
C:\lastools\bin>lasview -i fusa.laz -ilay 1
C:\lastools\bin>lasview -i fusa.laz -ilay 2

However, you may want to eventually apply the changes and produce a new LAZ file. This will be a lot slower as it requires rewriting the entire file. It will also make changes permanent. Press <CTRL-a> and a new file is produced called ‘fusa_1.laz’ that has 6 points less than ‘fusa.laz’ and 69 points with a different classification as “building”. One more thing, press <CTRL-x> if you want to toggle between the “cross” section view and the default view.

tutorial4_lasview_14_applied_laslayers You need to have a license to LAStools to save edits for file that contain 1 million points or more.

LAStools’ BLAST extension can process billions of LiDAR points

Often I get asked about the difference between las2dem.exe and blastdem.exe – the latter being part of the BLAST extension of LAStools. In the academic video from 2006 below, I outline the core technology details behind BLAST, which can seamlessly compute gigantic Delaunay Triangulations from billions of LiDAR points – using very little main memory. The blast2dem.exe module does not store the gigantic TINs that BLAST produces, but immediately rasters the finalized triangles streaming out of BLAST into equally massive DEMs. This avoids to ever having to store such huge triangulations – that are only needed temporarily – to disk. Currently blast2dem.exe can “only” process up to 2 billion points per file – a limitation that shall be lifted soon.