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.

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

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

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

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

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

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

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

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.

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.
This blog post is better than practically any journal article on the topic! Thanks for posting all the details!
Can’t agree more. And most of all it can be easily replicated too! Thanks.
Very nice indeed. I wonder how high was the flight (and how good is good enough) to get a resolution of these canopy for 100m x 100m sub-window is? What could be possible to use quantify or visualize pulse density with flight height to get good CHM – in LASTools. Would the pits also be explained say if the leaves were partially wet from water for example or too moist from mist or haze and are not really artifacts? Maybe morphological operation of fill holes can also work – I’m guessing.
very helpful, Thanks.
Any suggestions for computing CHMs in densely vegetated areas with very little points reaching the ground? Too few ground points apparently result in very large triangles.
Leo, this year at SilviLaser 2015 we will introduce a new technique that may help you with computing vegetation raster products when only poor ground information is available …
Pingback: Two ASPRS awards for “pit-free” CHM algorithm | rapidlasso GmbH
Pingback: Generating Spike-Free Digital Surface Models from LiDAR | rapidlasso GmbH
Where can I find the „subcircle“ attribute??
The option ‚-subcircle 0.1‘ for 20 cm disks can be used with either the lasthin or the lasgrid tool. Also see the corresonding README files.
Hi. I’m having a problem here. I already have the .las file with only the highest returns within an area. I have to interpolate these highest returns into a 0.5m grid (to create the CHM) and keep the correspondind Z value for each point of my 0.5m grid raster. I’ve been trying all las2dem options but it just doesn’t keep the Z value for each pixel, can somebody help me on how to interpolate those returns and still keep the Z value information for each pixel of the grid? (Z value = elevation)
Thanks in advance.
Pingback: FANR5640/7640: LAStools2 – GIS Note
Pingback: LiDAR Canopy Height Model Resolution – Does it Matter? | Interpine Innovation
Hello,
I have a question to the pit-free method of Khosravipour et al.
The code shows:
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
My question is how can I merge .bil-files? Another question is how should .bil-files fit into the tool „lasgrid“? It says: The input LiDAR points can be in LAS, LAZ, BIN, SHP, or simple XYZ ASCII format.
I’m quiet a beginner, am I just dumb or is there any error?
Greetings
Our lasgrid tool can read BIL ASC, and DEM rasters by on-the-fly conversion from grid to points that is implemented in form of a LASreader that is part of the LASlib library. But you can also use a proper raster processing software such as GDAL to merge the rasters.
Can you re-link the data? The links lead to 404 Not Found
Done. Here is the new link: drawno.laz
Any suggestions on how to speed this process up? I like the quality of the CHMs that it produces, but can get decent CHMs using the mean- or median-smoothing with a rolling window using FUSION’s CanopyModel much more quickly, particularly if you run this commands in sequence. It looks like you could execute the various height-layer las2dem calls in parallel, but am wondering if there’s anything else you’d recommend?
I suggest you have a look at the improvement on the pit-free algorithm (called the spike-free algorithm) that does not need the initial normalization step and computes only a single spike-free TIN from which a single spike-free DSM is rasterized.
https://rapidlasso.com/2016/02/03/generating-spike-free-digital-surface-models-from-lidar/