02_roi.qmd

---
title: "Regions of Interest"
---

```{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/02_roi.R)

[lidRbook section](https://r-lidar.github.io/lidRbook/engine.html#engine-clip)

## Overview

This code demonstrates the selection of regions of interest from LiDAR data. Simple geometries like circles and rectangles are selected based on coordinates. Complex geometries are extracted from shapefiles to clip specific areas.

## Environment

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

# Load packages
library(lidR)
library(sf)
```

## Simple Geometries

### Load LiDAR Data and Inspect

We start by loading the LiDAR point cloud data and inspecting its header and the number of point records.

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

# Inspect the header and the number of point records
las@header
las@header$`Number of point records`
```

### Select Circular and Rectangular Areas

We can select circular and rectangular areas from the LiDAR data based on specified coordinates and radii or dimensions.

```{r select_circular_areas}
# Establish coordinates
x <- 203890
y <- 7358935

# Select a circular area
circle <- clip_circle(las = las, xcenter = x, ycenter = y, radius = 30)

# Inspect the circular area and the number of point records
circle
circle@header$`Number of point records`
```

``` r
# Plot the circular area
plot(circle)
```

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

```{r select_rectangular_areas}
# Select a rectangular area
rect <- clip_rectangle(las = las, xleft = x, ybottom = y, xright = x + 40, ytop = y + 30)
```

``` r
# Plot the rectangular area
plot(rect)
```

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

```{r select_multiple}
# Select multiple random circular areas
x <- runif(2, x, x)
y <- runif(2, 7358900, 7359050)

plots <- clip_circle(las = las, xcenter = x, ycenter = y, radius = 10)
```

``` r
# Plot each of the multiple circular areas
plot(plots[[1]])
```

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

``` r
# Plot each of the multiple circular areas
plot(plots[[2]])
```

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

## Extraction of Complex Geometries from Shapefiles

In this section, we demonstrate how to extract complex geometries from shapefiles using the `clip_roi()` function from the `lidR` package.

::: callout-caution
## Legacy packages

`maptools`, `rgdal`, and `rgeos`, underpinning the `sp` package, will retire in October 2023. Please refer to R-spatial evolution reports for details, especially https://r-spatial.org/r/2023/05/15/evolution4.html.
:::

```{r extraction_complex_geometries_sf, warning=FALSE, message=FALSE}
# Load the shapefile using sf
planting <- sf::st_read(dsn = "data/shapefiles/MixedEucaNat.shp", quiet = TRUE)

# Plot the LiDAR header information without the map
plot(las@header, map = FALSE)

# Plot the planting areas on top of the LiDAR header plot
plot(planting, add = TRUE, col = "#08B5FF39")

# Extract points within the planting areas using clip_roi()
eucalyptus <- clip_roi(las = las, geometry = planting)
```

``` r
# Plot the extracted points within the planting areas
plot(eucalyptus)
```

```{r plot_euc, echo = FALSE, rgl = TRUE, fig.width = 8, fig.height = 6}
# Plot the extracted points within the planting areas
plot(eucalyptus, bg = "white")
```

## Exercises and Questions

Now, let's read a shapefile called `MixedEucaNatPlot.shp` using `sf::st_read()` and plot it on top of the LiDAR header plot.

``` r
# Read the shapefile "MixedEucaNatPlot.shp" using st_read()
plots <- sf::st_read(dsn = "data/shapefiles/MixedEucaNatPlot.shp", quiet = TRUE)

# Plot the LiDAR header information without the map
plot(las@header, map = FALSE)

# Plot the extracted points within the planting areas
plot(plots, add = TRUE)
```

#### E1.

Clip the 5 plots with a radius of 11.3 m.

#### E2.

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

#### 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`

## Conclusion

This concludes our tutorial on selecting simple geometries and extracting complex geometries from shapefiles using the `lidR` package in R.