09_solutions.qmd

---
title: "Excercise Solutions"
---

```{r, echo = FALSE, warnings = FALSE}
library(rgl)

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) 
options(lidR.progress = FALSE)
```

## Resources

[Code](https://github.com/tgoodbody/lidRtutorial/tree/main/code/solutions)

## 1-LAS

```{r, warning=FALSE, message=FALSE}
# Load packages
library(lidR)
library(sf)
library(terra)
```

#### E1.

What are withheld points? Where are they in our pointcloud?

```{r}
#| code-fold: true
# According to ASPRS LAS specification http://www.asprs.org/wp-content/uploads/2019/07/LAS_1_4_r15.pdf page 18 "a point # that should not be included in processing (synonymous with Deleted)"

# They are on the edges. It looks like they correspond to a buffer. LAStools makes use of the withheld bit to flag some # points. Without more information on former processing step it is hard to say.
```

#### E2.

Read the file dropping the withheld points.

```{r}
#| code-fold: true
las <- readLAS("data/MixedEucaNat_normalized.laz", filter = "-drop_withheld")
plot(las)
```

#### E3.

The withheld points seem to be legitimate points that we want to keep.

Try to load the file including the withheld points but get rid of the warning (without using `suppressWarnings()`). Hint: Check available `-set_withheld` filters in `readLAS(filter = "-h")`

```{r}
#| code-fold: true
las <- readLAS("data/MixedEucaNat_normalized.laz", filter = "-set_withheld_flag 0")
plot(las, color = "Withheld_flag")
```

#### E4.

Load only the ground points and plot the point-cloud coloured by the returnnumber of the point. Do it loading the strict minimal amount of memory (4.7 Mb). Hint: use `?lidR::readLAS` and see what `select` options might help.

```{r}
#| code-fold: true
las <- readLAS("data/MixedEucaNat_normalized.laz", filter = "-keep_class 2 -set_withheld_flag 0", select = "r")
plot(las, color = "ReturnNumber", legend = T)
format(object.size(las), "Mb")
```

## 2-ROI

```{r}
plots <- st_read("data/shapefiles/MixedEucaNatPlot.shp")
plot(las@header, map = FALSE)
plot(plots, add = TRUE)
```

#### E1.

Clip the 5 plots with a radius of 11.3 m,

```{r}
#| code-fold: true
inventory <- clip_roi(las, plots, radius = 11.3)
plot(inventory[[2]])
```

#### E2.

Clip a transect from A c(203850, 7358950) to B c(203950, 7959000).

```{r}
#| code-fold: true
tr <- clip_transect(las, c(203850, 7358950), c(203950, 7359000), width = 5)
plot(tr, axis = T)
```

#### E3.

Clip a transect from A c(203850, 7358950) to B c(203950, 7959000) but reorient it so it is no longer on the XY diagonal. Hint = ?clip_transect

```{r, eval = FALSE}
#| code-fold: true
ptr <- clip_transect(las, c(203850, 7358950), c(203950, 7359000), width = 5, xz = TRUE)
plot(tr, axis = T)
plot(ptr, axis = T)
plot(ptr$X, ptr$Z, cex = 0.25, pch = 19, asp = 1)
```

## 3-ABA

```{r}
las <- readLAS("data/MixedEucaNat_normalized.laz", select = "*",  filter = "-set_withheld_flag 0")
```

#### E1.

Assuming that biomass is estimated using the equation `B = 0.5 * mean Z + 0.9 * 90th percentile of Z` applied on first returns only, map the biomass.

```{r}
#| code-fold: true
B <- pixel_metrics(las, ~0.5*mean(Z) + 0.9*quantile(Z, probs = 0.9), 10, filter = ~ReturnNumber == 1L)
plot(B, col = height.colors(50))

B <- pixel_metrics(las, .stdmetrics_z, 10)
B <- 0.5*B[["zmean"]] + 0.9*B[["zq90"]]
plot(B, col = height.colors(50))

pixel_metrics(las, ~as.list(quantile(Z), 10))
```

#### E2.

Map the density of ground returns at a 5 m resolution with `pixel_metrics(filter = ~Classification == LASGROUND)`.

```{r}
#| code-fold: true
GND <- pixel_metrics(las, ~length(Z)/25, res = 5, filter = ~Classification == LASGROUND)
plot(GND, col = heat.colors(50))
```

#### E3.

Map pixels that are flat (planar) using `stdshapemetrics`. These could indicate potential roads.

```{r}
#| code-fold: true
m <- pixel_metrics(las, .stdshapemetrics, res = 3)
plot(m[["planarity"]], col = heat.colors(50))
flat <- m[["planarity"]] > 0.85
plot(flat)
```

## 5-DTM

#### E1.

Plot and compare these two normalized point-clouds. Why do they look different? Fix that. Hint: filter.

Some non ground points are below 0. It can be slightly low noise point not classified as ground by the data provider. This low points not being numerous and dark blue we hardly see them

```{r}
#| code-fold: true
las1 <- readLAS("data/MixedEucaNat.laz", filter = "-set_withheld_flag 0")
nlas1 <- normalize_height(las1, tin())
nlas2 <- readLAS("data/MixedEucaNat_normalized.laz", filter = "-set_withheld_flag 0")
plot(nlas1)
plot(nlas2)

nlas1 <- filter_poi(nlas1, Z > -0.1)
plot(nlas1)
```

#### E2.

Clip a plot somewhere in `MixedEucaNat.laz` (the non-normalized file).

```{r}
#| code-fold: true
circ <- clip_circle(las, 203930, 7359000, 25)
plot(circ)
```

#### E3.

Compute a DTM for this plot. Which method are you choosing and why?

```{r}
#| code-fold: true
dtm <- rasterize_terrain(circ, 0.5, kriging())
plot_dtm3d(dtm)
```

#### E4.

Compute a DSM (digital surface model). Hint: Look back to how you made a CHM.

```{r}
#| code-fold: true
dsm <- rasterize_canopy(circ, 1, p2r(0.1))
plot(dsm, col = height.colors(50))
```

#### E5.

Normalize the plot.

```{r}
#| code-fold: true
ncirc <- circ - dtm
plot(ncirc)
```

#### E6.

Compute a CHM.

```{r}
#| code-fold: true
chm <- rasterize_canopy(ncirc, 1, p2r(0.1))
plot(chm, col = height.colors(50))
```

#### E7.

Estimate some metrics of interest in this plot with cloud_metric()

```{r}
#| code-fold: true
metrics <- cloud_metrics(ncirc, .stdmetrics_z)
metrics
```

## 6-ITS

Using:

```{r}
las <- readLAS("data/example_corrupted.laz", select = "xyz")
col1 <- height.colors(50)
```

#### E1.

Run `las_check()` and fix the errors.

```{r}
#| code-fold: true
las_check(las)

las <- filter_duplicates(las = las)

las_check(las)
```

#### E2.

Find the trees and count the trees.

```{r}
#| code-fold: true
ttops <- locate_trees(las = las, algorithm = lmf(ws = 3, hmin = 5))
x <- plot(las)
add_treetops3d(x = x, ttops = ttops)
```

#### E3.

Compute and map the density of trees with a 10 m resolution.

```{r}
#| code-fold: true
r <- terra::rast(x = ttops)
terra::res(r) <- 10
r <- terra::rasterize(x = ttops, y = r, "treeID", fun = 'count')
plot(r, col = viridis::viridis(20))
```

#### E4.

Segment the trees.

```{r}
#| code-fold: true
chm <- rasterize_canopy(las = las, res = 0.5, algorithm = p2r(subcircle = 0.15))
plot(chm, col = col1)
ttops <- locate_trees(las = chm, algorithm = lmf(ws = 2.5))
las <- segment_trees(las = las, dalponte2016(chm = chm, treetops = ttops))

plot(las, color = "treeID")
```

#### E5.

Assuming that a value of interest of a tree can be estimated using the crown area and the mean Z of the points with the formula `2.5 * area + 3 * mean Z`. Estimate the value of interest of each tree.

```{r}
#| code-fold: true
value_of_interest <- function(x,y,z)
{
  m <- stdtreemetrics(x,y,z)
  avgz <- mean(z)
  v <- 2.5*m$convhull_area + 3 * avgz
  return(list(V = v))
}

V <- crown_metrics(las = las, func = ~value_of_interest(X,Y,Z))
plot(x = V["V"])

# 6. Map the total biomass at a resolution of 10 m. The output is a mixed of ABA and ITS

Vtot <- rasterize(V, r, "V", fun = "sum")
plot(Vtot, col = viridis::viridis(20))
```

## 7-LASCTALOG

This exercise is complex because it involves options not yet described. Be sure to use the lidRbook and package documentation.

https://cran.r-project.org/web/packages/lidR/lidR.pdf https://r-lidar.github.io/lidRbook/index.html

Using:

```{r}
ctg <- readLAScatalog(folder = "data/Farm_A/")
```

#### E1.

Generate a raster of point density for the provided catalog. Hint: Look through the documentation for a function that will do this!

```{r}
#| code-fold: true
ctg <- readLAScatalog("data/Farm_A/", filter = "-drop_withheld -drop_z_below 0 -drop_z_above 40")
D1 <- rasterize_density(las = ctg, res = 4)
plot(D1, col = heat.colors(50))
```

#### E2.

Modify the catalog to have a point density of 10 pts/m2 using the `decimate_points()` function. If you get an error make sure to read the documentation for `decimate_points()` and try: using `opt_output_file()` to write files to a temporary directory.

https://r-lidar.github.io/lidRbook/engine.html#engine-dtm-ondisk

```{r, eval = FALSE}
#| code-fold: true
newctg <- decimate_points(las = ctg, algorithm = homogenize(density = 10, res = 5))
#>  Error: This function requires that the LAScatalog provides an output file template.
```

```{r}
#| code-fold: true
opt_filter(ctg) <- "-drop_withheld"
opt_output_files(ctg) <- paste0(tempdir(), "/{ORIGINALFILENAME}")
newctg <- decimate_points(las = ctg, algorithm = homogenize(density = 10, res = 5))
```

#### E3.

Generate a raster of point density for this new decimated dataset.

```{r}
#| code-fold: true
opt_output_files(newctg) <- ""
D2 <- rasterize_density(las = newctg, res = 4)
plot(D2, col = heat.colors(50))
```

#### E4.

Read the whole decimated catalog as a single las file. The catalog isn't very big - not recommended for larger data sets!

```{r}
#| code-fold: true
las <- readLAS(newctg)
plot(las)
```

#### E5.

Read documentation for the catalog_retile() function and merge the dataset into larger tiles. Use `ctg` metadata to align new chunks to the lower left corner of the old ones. Hint: Visualize the chunks and use `opt_chunk_*` options.

```{r}
#| code-fold: true
opt_chunk_size(ctg) <- 280
opt_chunk_buffer(ctg) <- 0
opt_chunk_alignment(ctg) <- c(min(ctg$Min.X), min(ctg$Min.Y))
plot(ctg, chunk = T)

opt_output_files(ctg) <- "{tempdir()}/PRJ_A_{XLEFT}_{YBOTTOM}"
newctg <- catalog_retile(ctg = ctg)
plot(newctg)
```

## 8-ENGINE

#### E1.

In example 2 (section B) what does last line `m <- m[m$treeID %in% p$treeID,]` do? Adjust the function to not include that line to see what happens (use `catalog_select()` to select 4 tiles to test on).

```{r, eval = FALSE}
#| code-fold: true
# Subset catalog
subctg <- catalog_select(ctg)

# without line
routine_trees_test <- function(chunk) {
  # Read in check, check NULL status, get bbox
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)
  bbox <- st_bbox(chunk)
  
  # Filter surface points and generate chm
  las <- filter_surfacepoints(las, res = 0.5)
  chm <- rasterize_canopy(las = las, res = 0.5, algorithm = p2r())
  
  # Tree detection, segmentation, metrics
  ttops <- locate_trees(las = las, algorithm = lmf(ws = 3, hmin = 5))
  las_trees <- segment_trees(las = las, algorithm = dalponte2016(chm = chm, treetops = ttops))
  p <- crown_metrics(las = las_trees, func = .stdtreemetrics)
  
  # Remove buffer
  p <- sf::st_crop(x = p, y = bbox)
  
  # Delineate crowns
  output <- delineate_crowns(las_trees)
  
  #output <- m[m$treeID %in% p$treeID,]
  
  return(output)
}

options <-  list(automerge = TRUE)
m <-  catalog_apply(subctg, routine_trees_test, .options = options)
plot(m, col = rgb(0,0,1,0.3))
```

```{r, eval = FALSE}
#| code-fold: true
ctg <-  readLAScatalog("data/Farm_A/")
opt_select(ctg) <- "xyz"
opt_filter(ctg) <- "-drop_withheld -drop_z_below 0 -drop_z_above 40"
opt_chunk_buffer(ctg) <- 15
opt_chunk_size(ctg) <- 0
subctg <-  catalog_select(ctg)
options <-  list(automerge = TRUE)
m <- catalog_apply(subctg, routine_trees_test, .options = options)

plot(m, col = rgb(0,0,1,0.3))
```

#### E2.

The following is a simple (and a bit naive) function to remove high noise points. - Explain what this function does - Create a user-defined function to apply using `catalog_apply()` - Hint: Dont forget about buffered points... remember `lidR::filter_*` functions.

```{r}
#| code-fold: true
filter_noise <- function(las, sensitivity)
{
  p95 <- pixel_metrics(las, ~quantile(Z, probs = 0.95), 10)
  las <- merge_spatial(las, p95, "p95")
  las <- filter_poi(las, Z < 1+p95*sensitivity, Z > -0.5)
  las$p95 <- NULL
  return(las)
}

filter_noise_collection = function(cl, sensitivity)
{
  las <- readLAS(cl)
  if (is.empty(las)) return(NULL)
  las <- filter_noise(las, sensitivity)
  las <- filter_poi(las, buffer == 0L)
  return(las)
}

ctg = readLAScatalog("data/Farm_A/")
opt_select(ctg) <- "*"
opt_filter(ctg) <- "-drop_withheld -drop_"
opt_output_files(ctg) <- "{tempdir()}/*"
opt_chunk_buffer(ctg) <- 20
opt_chunk_size(ctg) <- 0

options <- list(automerge = TRUE)
output <- catalog_apply(ctg, filter_noise_collection, sensitivity = 1.2, .options = options)

las <- readLAS(output)
plot(las)
```

#### E3.

Design an application that retrieves the convex hull of each flightline (hard). Use the `concaveman::concaveman()` function, adn functions from `sf`. Start by designing a test function that works on a LAS object and later apply on the collection. The output should look like:

```{r}
flightlines <- st_read("data/flightlines.shp")
plot(flightlines, col = sf.colors(6, alpha = 0.5))
plot(flightlines[3,])
```

```{r}
#| code-fold: true
# Read the catalog
ctg <- readLAScatalog("data/Farm_A/")

# Read a single file to perform tests
las <- readLAS(ctg$filename[16], select = "xyzp", filter = "-drop_withheld -drop_z_below 0 -drop_z_above 40")

# Define a function capable of building the hull from the XY of a given PointSourceID
enveloppes <- function(x,y, psi)
{
  hull <- concaveman::concaveman(cbind(x,y), length_threshold = 10)
  hull <- sf::st_polygon(list(hull))
  hull <- sf::st_sfc(hull)
  hull <- sf::st_simplify(hull, dTolerance = 1)
  hull <- sf::st_sf(hull)
  hull$ID <- psi[1]
  list(hull = list(hull = hull))
}

# Define a function that apply the previous function to each PointSourceID from a LAS object
flighline_polygons <- function(las)
{
  u <- las@data[ , enveloppes(X,Y, PointSourceID), by = PointSourceID]
  hulls <- Reduce(rbind, u$hull)
  return(hulls)
}

# Test this function on a LAS
hulls <- flighline_polygons(las)
plot(hulls, col = sf.colors(3, alpha = 0.5))


# It works so let make a function that works with a LAScatalog
flighline_polygons <- function(las)
{
  if (is(las, "LAS"))  {
    u <- las@data[ , enveloppes(X,Y, PointSourceID), by = PointSourceID]
    hulls <- Reduce(rbind, u$hull)
    return(hulls)
  }
  
  if (is(las, "LAScluster")) {
    las <- readLAS(las)
    if (is.empty(las)) return(NULL)
    hulls <- flighline_polygons(las)
    return(hulls)
  }
  
  if (is(las, "LAScatalog")) {
    opt_select(las) <-  "xyzp"
    options <- list(
      need_output_file = FALSE,
      need_buffer = TRUE,
      automerge = TRUE)
    output <- catalog_apply(las, flighline_polygons, .options = options)
    hulls <- dplyr::summarise(dplyr::group_by(output, ID), ID = ID[1])
    return(hulls)
  }
  
  stop("Invalid input")
}

library(future)
future::plan(multisession)
opt_chunk_buffer(ctg) <- 5
opt_filter(ctg) <- "-drop_withheld -drop_z_below 0 -drop_z_above 40"
flightlines <- flighline_polygons(ctg)
plot(flightlines, col = sf.colors(6, alpha = 0.5))
```