Rasterizing perfect canopy height models
Document updated on February 13th 2020 and up-to-date with
lidR v3.0.0
The following page shows how to compute a canopy height model in the R environment using the lidR package. The file used in this work can be downloaded here. This work is the "lidR version" of the blog article published on rapidlasso.com by Martin Isenburg and is therefore deeply inspired by it. We will use the same file to do (almost) the same things.
As Martin said in his article, in the literature, authors often compute a canopy height model (CHM) without describing the process that was used to create the CHM. There are many ways of achieving such a task and we will describe here what lidR proposes.
library(lidR)
las = readLAS("drawno.laz")
las = normalize_height(las, tin())
plot(las)
The first and simplest method consists of returning the highest z coordinates among all LiDAR returns falling into the corresponding n by n meter area (e.g. 0.5 by 0.5 m). It is a point-to-raster method.
chm1 = grid_canopy(las, 0.5, p2r())
plot(chm1, col = height.colors(50))
At this resolution, the resulting CHM is full of empty pixels and so-called "pits". Theses empty pixels and pits are harmful for consistency of subsequent analyses. A simple improvement proposed by Martin Isenburg in lastools can be obtained by replacing each LiDAR return with a small disk. Since a laser has a width and a footprint, this tweak may simulate this fact. The lidR package implements such a feature too.
chm2 = grid_canopy(las, 0.5, p2r(0.2))
plot(chm2, col = height.colors(50))
The resulting CHM no longer contains empty pixels. However, the resolution is still coarse. Let's try with a higher resolution (e.g. 0.25 by 0.25 m)
chm3 = grid_canopy(las, 0.25, p2r(0.2))
plot(chm3, col = height.colors(50))
Here there is a trade-off between the subcircle parameter and the resolution to avoid the few occurring empty pixels. An interpolation of empty pixels can also be done (look at parameter na.fill). However, gridding the highest returns will often leave empty pixels in the CHM at high resolution, even after "subcircling" the points. Moreover the CHM will still have pits.
Another popular approach to avoid the empty pixel problem consists of the interpolation of first returns with a triangulated irregular network (TIN) and then rasterizing it onto a grid to create the CHM. This can be implemented in lidR with grid_tincanopy. Then a linear interpolation is done within each triangle.
chm4 = grid_canopy(las, 0.25, dsmtin())
plot(chm4, col = height.colors(50))
As explained by Martin Isenburg, 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. In this context, Khosravipour et al. (2014) [1] developed a "pit-free" algorithm for lastools, and this algorithm is also implemented in lidR. Basically, it consists of several layers of triangulation at different elevations. In the following example at z = 0 (classical triangulation and at 2, 5, 10 and 15 meters.
chm5 = grid_canopy(las, 0.25, pitfree(thresholds = c(0,2,5,10,15), max_edge = c(0,1)))
plot(chm5, col = height.colors(50))
This CHM looks nice. No empty pixels and all the pits are gone. However it is still not perfect. Some shallow pits still remain. This maybe can be solved by playing with the max_edge parameter, but let's try another approach.
The pit-free algorithm combined with the subcircling tweak works and is available in lidR and we get a really nice CHM without empty pixels or pits!
chm6 = grid_canopy(las, 0.25, pitfree(c(0,2,5,10,15), c(0,1), subcircle = 0.2))
plot(chm6, col = height.colors(50))
However, one may guess that this last one is computationally demanding.
[1] Khosravipour, Anahita & Skidmore, Andrew & Isenburg, Martin & Wang, Tiejun & Hussin, Yousif. (2014). Generating Pit-free Canopy Height Models from Airborne Lidar. Photogrammetric Engineering & Remote Sensing. 80. 863-872. 10.14358/PERS.80.9.863.