08_engine2.qmd

---
title: "LAScatalog processing engine"
---

```{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)
```

## Relevant resources:

[Code](https://github.com/tgoodbody/lidRtutorial/blob/main/code/tutorials/08_engine2.R)

-   [lidRbook section: Engine](https://r-lidar.github.io/lidRbook/engine2.html)

-   [lidRbook section: Thinking outside the box](https://r-lidar.github.io/lidRbook/outbox.html)

## Overview

This code showcases the LASCATALOG PROCESSING ENGINE, which efficiently and in parallel applies various functions to LiDAR catalogs. It introduces the `catalog_apply()` function for processing LiDAR data in a catalog. The code includes routines to detect trees and calculate metrics on the LiDAR catalog.

## Environment

```{r clear warnings, warnings = FALSE, message=FALSE}
# Clear environment
rm(list = ls(globalenv()))

# Load packages
library(lidR)
library(terra)
library(future)
```

## Basic Usage

In this section, we will cover the basic usage of the `lidR` package, including reading LiDAR data, visualization, and inspecting metadata.

### Basic Usage of `lidR` Package

This section introduces the basic usage of the `lidR` package for reading and visualizing LiDAR data, as well as inspecting metadata.

### Reading and Visualizing LiDAR Data

We start by reading a LAS catalog and inspecting one of its LAS files.

```{r read_las_catalog}
# Read a LAS catalog
ctg <- readLAScatalog(folder = "data/Farm_A/")

# Inspect the first LAS file in the catalog
las_file <- ctg$filename[1]
las <- readLAS(las_file)
las
```

### Visualizing LiDAR Data

We visualize the LiDAR data from the selected LAS file using a 3D plot.

```{r visualize_las}
# Visualize the LiDAR data in 3D
plot(las, bg = "white")
```

### `catalog_apply()` Function

This section demonstrates the use of the `catalog_apply()` function for efficient processing of LiDAR data within a LAS catalog.

### Problem Statement

We start by addressing a common problem - how can we apply operations to `LAS` data in a catalog?

```{r problem}
# Read a LAS file from the catalog and filter surface points
las_file <- ctg$filename[16]
las <- readLAS(files = las_file, filter = "-drop_withheld -drop_z_below 0 -drop_z_above 40")
surflas <- filter_surfacepoints(las = las, res = 1)
```

### Visualizing LiDAR Data

We visualize the selected LiDAR data, including both the original data and the surface points.

``` r
# Visualize the LiDAR data with a default color palette
plot(las)
```

```{r visualize_las_ctg2, echo = FALSE, rgl = TRUE, fig.width = 8, fig.height = 6}
# Visualize the LiDAR data with a default color palette
plot(las, bg = "white")
```

``` r
# Visualize the surface points using a default color palette
plot(surflas)
```

```{r visualize_surflas_ctg2, echo = FALSE, rgl = TRUE, fig.width = 8, fig.height = 6}
# Visualize the surface points using a default color palette
plot(surflas, bg = "white")
```

### Calculating Rumple Index

We calculate the rumple index using the `pixel_metrics()` function.

```{r calculate_rumple}
# Generate Area-based metrics
ri <- pixel_metrics(las = las, ~rumple_index(X,Y,Z), res = 10)
plot(ri)
```

### Solution: `LAScatalog` Processing Engine

This section introduces the `LAScatalog` processing engine, a powerful tool for efficient processing of LAS data within a catalog.

### Basic Usage of the `catalog_apply()` Function

We demonstrate the basic usage of the `catalog_apply()` function with a simple user-defined function.

```{r basic_usage}
# User-defined function for processing chunks
routine <- function(chunk){
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)
  
  # Perform computation
  m <- pixel_metrics(las = las, func = ~max(Z), res = 20)
  output <- terra::crop(x = m, terra::ext(chunk))
  
  return(output)
}

# Initialize parallel processing
plan(multisession)

# Specify catalog options
opt_filter(ctg) <- "-drop_withheld"

# Apply routine to catalog
out <- catalog_apply(ctg = ctg, FUN = routine)

# Inspect the output list
out[1:5]

# Use the engine-supported method for merging
options <- list(automerge = TRUE)
out <- catalog_apply(ctg = ctg, FUN = routine, .options = options)
print(out)
```

### User-Defined Functions for Processing

We demonstrate the use of user-defined functions to process LiDAR data within a catalog.

```{r user_defined_functions}
# User-defined function for rumple index calculation
routine_rumple <- function(chunk, res1 = 10, res2 = 1){
  las <-  readLAS(chunk)
  if (is.empty(las)) return(NULL)
  bbox <- terra::ext(chunk)
  
  las <- filter_surfacepoints(las = las, res = res2)
  ri  <- pixel_metrics(las = las, ~rumple_index(X,Y,Z), res1)
  
  output <- terra::crop(x = ri, y = bbox)
  return(output)
}

# Set catalog options
opt_select(ctg) <- "xyz"
opt_filter(ctg) <- "-drop_withheld -drop_z_below 0 -drop_z_above 40"
opt_chunk_buffer(ctg) <- 0
opt_chunk_size(ctg) <- 0

# Specify options for merging
options <- list(automerge = TRUE, alignment = 10)

# Apply the user-defined function to the catalog
ri <- catalog_apply(ctg = ctg, FUN = routine_rumple, res1 = 10, res2 = 0.5, .options = options)

# Plot the output
plot(ri, col = height.colors(50))
```

## `catalog_apply()` - Example 2

In this section, we provide another example of using the `catalog_apply()` function to detect trees and calculate metrics on a catalog.

### Defining Routines for Tree Detection and Metrics

We define a routine that detects trees, calculates metrics, and returns relevant data.

```{r tree_detection_metrics}
# User-defined routine for tree detection and metrics
routine_trees <- function(chunk) {
  # Read in the chunk and check for emptiness
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)
  
  # Get the chunk bounding box
  bbox <- sf::st_bbox(obj = chunk)

  # Filter surface points and create canopy height model (CHM)
  las <- filter_surfacepoints(las, res = 0.5)
  chm <- rasterize_canopy(las = las, res = 0.5, algorithm = p2r())

  # Detect and segment trees
  ttops <- locate_trees(las = las, algorithm = lmf(ws = 3, hmin = 5))
  las_trees <- segment_trees(las = las, algorithm = dalponte2016(chm = chm, treetops = ttops))
  
  # Generate metrics for each tree
  p <- crown_metrics(las = las_trees, func = .stdtreemetrics)
  p <- sf::st_crop(x = p, y = bbox)

  # Delineate convex hulls
  m <- delineate_crowns(las_trees)
  output <- m[m$treeID %in% p$treeID,]

  return(output)
}

# Set options for the catalog
opt_chunk_buffer(ctg) <- 15
options <- list(automerge = TRUE) # Merge all outputs

# Apply the function to the catalog
m <- catalog_apply(ctg = ctg, FUN = routine_trees, .options = options)

# View and visualize the output
m
plot(m)

# End parallel processing
future::plan(sequential)
```

::: callout-note
## Thinking outside the box

The LAScatalog engine is versatile! The functions that can be applied to LiDAR data are infinite - leverage the flexibility of `lidR` and create software that pushes the boundaries of research in forest inventory and management!
:::

## Exercises

#### E1.

Understanding Chunk Filtering

On line 111 (routine_trees function) `m <- m[m$treeID %in% p$treeID,]` is used. Explain the purpose of this line. To understand its impact, modify the function to exclude this line and observe the results. You can use the `catalog_select()` function to choose a subset of tiles for testing.

``` r
subctg <- catalog_select(ctg)
```

#### E2.

Implement Noise Filtering

-   Explain the purpose of the `filter_noise()` function.
-   Create a user-defined function to apply noise filtering using the `catalog_apply()` function.
-   Make sure to consider buffered points when using lidR's `filter_*` functions.

#### E3.

Flightline Convex Hull Application

Design an application to retrieve the convex hull of each flightline using the `concaveman::concaveman()` function and functions from the `sf` package.

-   Begin by designing a test function that works on a single LAS object.
-   Apply the function to a collection of LAS files.
-   Visualize the results using the flightlines' shapefile.

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

## Conclusion

This concludes the tutorial on using the `catalog_apply` function in the `lidR` package to efficiently process LAS data within a catalog.

```{r, echo = FALSE, results = FALSE}
# Instructions for cleaning up any existing .lax files
# (Note: Please replace 'path' with the appropriate path)
path <- "data/Farm_A/"
file_list <- list.files(path)
delete_lax <- file_list[grep("\\.lax$", file_list)]
file.remove(file.path(path, delete_lax))
```