Smoothing of water areas of DTM using ALS points

When creating a Digital Terrain Models from LIDAR data, there are discontinuities (absence) of points over areas with water. Algorithms try to approximate missing points, but the result is generally inaccurate. In the case of the method based on triangular irregular network (TIN), such areas are usually represented by relatively large triangles (see Figure 2). Below is a method to smooth out such areas.

My idea was to find a water areas with help of lidRplugins::delineate_lakes(), calculate the mean (or min) values and merge it with original DTM file. The workflow is shown on Figure 1:

Figure 1: Algorithm for smoothing the water surface areas over DTM; r - rasters, v - vectors

Let’s read the DTM already crated. It’s shown on Figure 2.

r <- terra::vrt("data/dtm/rasterize_terrain.vrt")

Figure 2: DTM obtained from LAS with rasterize_terrain() function.

To find and delineate/vectorize the water areas we will use delineate_lakes() function from lidRplugins package (Roussel 2022). In result we will get a vector layer with areas of similar altitude. The tol and tol2 parameters are used to set the tolerances.

l <- lidR::readLAScatalog("data/las/") |>
  lidRplugins::delineate_lakes(tol = 1/1000, tol2 = c(1/10^6, 2/10^4)) |>
  terra::vect(l) |>
  terra::project(terra::crs(r))

As the output of delineate_lakes() can contain the flat terrain areas as well, we will filter out just water (-keep_class 9) points from LAScatalog, rasterize them and convert raster to polygon.

v <- lidR::readLAScatalog("data/las/", filter = "-keep_class 9") |>
  lidR::rasterize_density() |>
  terra::as.polygons() |>
  terra::subset(density > 0, NSE = TRUE) |>
  terra::project(terra::crs(r))

Now, having two polygons: 1st with delineated areas, and 2 only with water areas, lets intersect them to have only delineated water:

w <- l |>
  subset(terra::is.related(l, v, relation = "intersects"))

Now we are going to extract the mean (or you can use any other function) values of altitude for the vector areas and bind it to our vector layer.

w <- terra::extract(r, w, fun=min, bind = TRUE, na.rm=TRUE)
w <- setNames(w, "mean_height")
terra::writeVector(w, "data/w.gpkg", overwrite = TRUE)

terra::plot(r, col = grDevices::gray.colors(50))
terra::plot(w, add = TRUE)

Figure 3: DTM with vectorized water

And finally produce the output raster:

ra <- r
ra[] <- NA
ra <- terra::rasterize(w,ra, field = w$mean_height)

indx <- !is.na(ra)
rb <- r
rb[indx] <- ra[indx]
terra::writeRaster(rb, file = "data/rb.tif", overwrite = TRUE)

Figure 4: Raster layer with rasterized aggregaded values for water areas

Figure 5: DTM with smoothed water areas.

References

Roussel, Jean-Romain. 2022. lidRplugins: Extra Functions and Algorithms for lidR Package.