dtm.qmd

```{r,echo=FALSE,message=FALSE,warning=FALSE}
library(lidR)
library(ggplot2)
library(raster)
library(stars)
#library(forcats)

source("function_plot_crossection.R")

r3dDefaults = rgl::r3dDefaults
m = structure(c(0.921, -0.146, 0.362, 0, 0.386, 0.482, -0.787, 0, 
-0.06, 0.864, 0.5, 0, 0, 0, 0, 1), .Dim = c(4L, 4L))
r3dDefaults$FOV = 50
r3dDefaults$userMatrix = m
r3dDefaults$zoom = 0.75

knitr::opts_chunk$set(
  comment =  "#>", 
  collapse = TRUE,
  fig.align = "center")

rgl::setupKnitr(autoprint = TRUE)
```

# Digital terrain model {#sec-dtm}

Generating a Digital Terrain Model (DTM) is usually the second step in processing that follows classification of ground points (@sec-gnd). Put simply, a DTM can be described as an "image" of the ground. Methods to generate DTMs have been intensively studied and several algorithms have been proposed for various terrain situations. DTMs are used for a variety of purposes in practice, such as determination of the catchment basins of water retention and stream flow, or the identification of drivable roads to access resources. It also enables users to normalize point clouds i.e. subtract the local terrain from the elevation of points to allow a manipulation of point clouds as if they were acquired on a flat surface (@sec-norm).

The construction of a DTM starts with known or sampled ground points and uses various spatial interpolation techniques to infer ground points at unsampled locations. Accuracy of the DTM is very important because errors will propagate to future processing stages like tree height estimation. A wide range of methods exist for spatial interpolation of points. In the following section we will use the classified `Topography.laz` data set, which is included internally within `lidR` to create reproducible examples.

```{r plot-las-topo-dtm, warning=FALSE, rgl=TRUE}
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzc")
plot(las, size = 3, bg = "white")
```

## Triangular irregular network {#sec-tin}

This method is based on triangular irregular network (TIN) of ground point data to derive a bivariate function for each triangle, which is then used to estimate the values at unsampled locations (between known ground points).

```{r triangulation, echo = FALSE, snapshot = TRUE}
las2 <- filter_ground(las)
dxyz <- with(las2@data, deldir::deldir(X, Y, z = Z, suppressMsge = TRUE))
col <- with(las2@data, height.colors(20)[1 + round(19*(Z - min(Z))/diff(range(Z)))])
rgl::persp3d(dxyz, col = col, smooth = FALSE, add = T, front="lines", axes = F)
rgl::axes3d(c("x-", "y+", "z-"))
rgl::grid3d(side=c('x-', 'y+', 'z-'), col="gray")
```

Planar facets of each generated triangle are used to interpolate. Used with a Delaunay triangulation, this is the most simple solution because it involves no parameters. The Delaunay triangulation is unique and the linear interpolation is parameter-free. The drawbacks of the method are that it creates a non-smooth DTM and that it cannot extrapolate the terrain outside the convex hull delimited by the ground points since there are no triangle facets outside the convex hull. Moreover, the interpolation is weak at the edges because large irrelevant triangles are often created. It's therefore important to compute the triangulation with a buffer to be able to crop the DTM and clear the edge artifacts (see @sec-engine).

To generate a DTM model with the TIN algorithm we use `rasterize_terrain()` where `algorithm = tin()`.

```{r tin-dtm, rgl = TRUE}
dtm_tin <- rasterize_terrain(las, res = 1, algorithm = tin())
plot_dtm3d(dtm_tin, bg = "white") 
```

Notice the ugly edge interpolations. This occurs because we didn't process with a buffer.

## Invert distance weighting {#sec-idw}

Invert distance weighting (IDW) is one of the simplest and most readily available methods that can be applied to create DTMs. It is based on an assumption that the value at an unsampled point can be approximated as a weighted average of values at points within a certain cut-off distance *d*, or from a given number *k* of closest neighbours. Weights are usually inversely proportional to a power *p* of the distance between the location and the neighbour, which leads to the computing of an estimator.

Compared to `tin()` this method is more robust to edge artifacts because it uses a more relevant neighbourhood but generates terrains that are "bumpy" and probably not as realistic as those generated using TINs. There are always trade-offs to different methods!

To generate a DTM model with the IDW algorithm we use `rasterize_terrain()` where `algorithm = knnidw()`.

```{r idw-dtm, rgl = TRUE}
dtm_idw <- rasterize_terrain(las, algorithm = knnidw(k = 10L, p = 2))
plot_dtm3d(dtm_idw, bg = "white") 
```

Notice the bumpy nature of the DTM compared to the previous one generated with `tin()`. In 1D and IDW interpolation looks like:

```{r plot-1d-idw, echo = FALSE, fig.height=2, fig.width=6, message=FALSE}
set.seed(43)
x <- runif(20, 0, 20)
z <- 2*sin(x)+2
y <- rep(0, 20)

llas = LAS(data.frame(X=x, Y = y, Z = z))

X <- seq(0,20,0.1)
Y <- rep(201, 0)

Z <- lidR:::C_knnidw(llas, X, Y, 7, 2.5, 100, 1)

opar <- par(mar = c(0,0,0,0))
plot(x,z, asp = 1)
lines(X,Z)
par(opar)
```

## Kriging {#sec-kriging}

Kriging is the most advanced approach and utilizes advanced geostatistical interpolation methods that take into account the relationships between the returns and their respective distances from each other. `lidR` uses the package `gstat` to perform the kriging. This method is very advanced, difficult to manipulate, and extremely slow to compute, but probably provides the best results with minimal edge artifacts.

To generate a DTM model with the kriging algorithm we use `rasterize_terrain()` where `algorithm = kriging()`.

```{r kri-dtm, message=F, warning=F, rgl = TRUE}
dtm_kriging <- rasterize_terrain(las, algorithm = kriging(k = 40))
plot_dtm3d(dtm_kriging, bg = "white") 
```

Notice that the algorithm has issues interpolating regions with missing point such as lakes.

## Pros and cons {#sec-dtm-pros-cons}

-   **Triangulation** is a very fast and efficient method that generates very good DTMs and is robust to empty regions inside the point cloud. It is however weak at edges. Although `lidR` uses the nearest neighbour to complete the missing pixel out of the convex hull of the ground points the interpolation remains poor. This algorithm must therefore always be used with a buffer of extra points to ensure that the region of interest is not on an edge. The TIN method is recommended for broad DTM computation but should be avoided for small regions of interest loaded without buffers.

-   **Invert distance weighting** is fast, but approximately twice as slower than TIN. The terrain is not very realistic, but edges are likely to be free of strong edge artifacts. IDW is a compromise between TIN and KRIGING. It is recommended if you want a simple method, if you cannot load a buffer, and if edge regions are important.

-   **Kriging** is very slow because it is computationally demanding. It is not recommended for use on medium to large areas. It can be used for small plots without buffers to get a nice DTM without strong edges artifact.

Whatever the method used, edges are critical. Results will always be weak if the method needs to guess the local topography with only partial information on the neighborhood. Though different methods provide better and worse estimates in these regions, best practice is to **always** use a buffer to obtain some information about the neighborhood and remove the buffer once the terrain is computed.

## Other methods {#sec-dtm-other}

Spatial interpolation is not limited to the 3 methods described above. Many more have been presented and described in the literature. In @sec-plugins we will learn how to create a plugin algorithm compatible with `rasterize_terrain()` based on a multilevel B-spline approximation (MBA) using the [`MBA`](https://cran.r-project.org/web/packages/MBA/index.html) package.

```{r, echo = FALSE}
mba <- function(n = 1, m = 1, h = 8, extend = TRUE)
{
  f <- function(las, where, scales = c(0,0), offsets = c(0,0)) {
    res <- MBA::mba.points(las@data, where, n, m , h, extend)
    return(res$xyz.est[,3])
  }
  class(f) <- lidR:::LIDRALGORITHMSPI
  return(f)
}
```

```{r plot-dtm-3d, rgl = TRUE}
dtm_mba <- rasterize_terrain(las, algorithm = mba())
plot_dtm3d(dtm_mba, bg = "white")
```

## Render shaded DTM {#sec-hillshade}

```{r, echo=FALSE, print=FALSE}
suppressMessages(library(terra))
```

Generating a hillshade layer in R is relatively straight forward and is done using functions from the `terra` package. The `terrain()` and `hillShade()` functions can be combined to take the DTM raster layers as input and return a hillshade raster:

```{r plot-dtm-shaded, fig.height=5.85, fig.width=6, warning=FALSE}
library(terra)
dtm <- rasterize_terrain(las, algorithm = tin())
dtm_prod <- terrain(dtm, v = c("slope", "aspect"), unit = "radians")
dtm_hillshade <- shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
plot(dtm_hillshade, col =gray(0:30/30), legend = FALSE)
```

The [`rayshader`](https://github.com/tylermorganwall/rayshader) package also provides interesting tools to generate shaded DTM. The dtm must be a `RasterLayer`

```{r, message=FALSE}
library(rayshader)
dtm <- raster::raster(dtm)
elmat <- raster_to_matrix(dtm)
map <- elmat %>%
  sphere_shade(texture = "imhof1", progbar = FALSE) %>%
  add_water(detect_water(elmat), color = "imhof1") %>%
  add_shadow(ray_shade(elmat, progbar = FALSE), 0.5) %>%
  add_shadow(ambient_shade(elmat, progbar = FALSE), 0)
```

2D plot

```{r plot-2d-dtm-rayshaded}
plot_map(map)
```

3D plot

```{r plot-3d-dtm-rayshaded, rgl = TRUE, warning  = FALSE, fig.height=15}
plot_3d(map, elmat, zscale = 1, windowsize = c(800, 800))
```