io.qmd

```{r,echo=FALSE,message=FALSE,warning=FALSE}
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")

library(lidR)
library(sf)
library(stars)
library(ggplot2)

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile)

options(crayon.enabled = TRUE)
rgl::setupKnitr(autoprint = TRUE)
```

# Reading, Plotting, Querying & Validating {#sec-io}

## Reading LiDAR data using `readLAS` {#sec-read}

Discrete return ALS sensors record various types of data. Primarily, they capture positional data in three dimensions (X, Y, Z), followed by additional information like the intensity for each point, the position of each point in the return sequence, and the beam incidence angle of each point. Reading, writing, and efficient storage of these ALS data are critical steps prior to any subsequent analysis.

ALS data are most commonly distributed in LAS format, which is specifically designed to store ALS data in a standardized way. These data are officially documented and maintained by the [ASPRS](http://www.asprs.org/wp-content/uploads/2019/07/LAS_1_4_r15.pdf). However, LAS files require a large amount of memory because they are not compressed. The LAZ format has become the standard compression scheme because it is free and open-source.

The widespread use, standardization, and open-source nature of the LAS and LAZ formats promoted the development of the `lidR` package. This package is designed to process LAS and LAZ files both as input and output, leveraging the LASlib and LASzip C++ libraries via the [`rlas`](https://cran.r-project.org/package=rlas) package.

The function `readLAS()` reads a LAS or LAZ file and returns an object of class `LAS`. The `LAS` formal class is documented in detail in a [dedicated vignette](https://cran.r-project.org/web/packages/lidR/vignettes/lidR-LAS-class.html). To briefly summarize, a LAS file consists of two parts:

1. The header, which stores summary information about its content, including the bounding box of the file, coordinate reference system, and point format.
2. The payload, i.e., the point cloud itself.

The function `readLAS()` reads and creates an object that contains both the header and the payload.

``` r
las <- readLAS("files.las")
```

When printed it displays a summary of its content.

```{r print-las}
print(las)
```

For a more in-depth print out of the data use the function `summary()` instead of `print()`.

### Parameter `select` {#sec-select}

A LAS file stores the `X Y Z` coordinates of each point as well as many other data such as intensity, incidence angle, and return sequence position. These data are called *attributes*. In practice, many attributes are not actually useful but are loaded by default. This can consume a lot of processing memory because R does not allow for choosing data storage modes (see [this vignette](https://cran.r-project.org/web/packages/lidR/vignettes/lidR-LAS-class.html) for more details).

To save memory, `readLAS()` can take an optional parameter `select`, which enables the user to selectively load the attributes of interest. For example, one can choose to load only the `X Y Z` attributes.

```r
las <- readLAS("file.las", select = "xyz")  # load XYZ only
las <- readLAS("file.las", select = "xyzi") # load XYZ and intensity only
```

Examples of other attribute abbreviations are: `t` - gpstime, `a` - scan angle, `n` - number of returns, `r` - return number, `c` - classification, `s` - synthetic flag, `k` - keypoint flag, `w` - withheld flag, `o` - overlap flag (format 6+), `u` - user data, `p` - point source ID, `e` - edge of flight line flag, `d` - direction of scan flag

### Parameter `filter` {#sec-filter}

While `select` enables the user to choose "columns" (or attributes) while reading files, `filter` allows selection of "rows" (or points) during the reading process. Removing superfluous data at read time saves memory and increases computation speed. For example, it's common practice in forestry to process only the first returns.

```r
las <- readLAS("file.las", filter = "-keep_first") # Read only first returns
```

It is important to understand that the `filter` option in `readLAS()` keeps or discards points **at read time**, i.e., while reading at the C++ level, without involving any R code. For example, the R function `filter_poi()` may return the same output as the `filter` option in `readLAS()`:


```r
las1 <- readLAS("file.las", filter = "-keep_first")

las2 <- readLAS("file.las")
las2 <- filter_poi(las2, ReturnNumber == 1L)
```

In the example above, we are (1) reading only the first returns or (2) reading all the points and then filtering the first returns in R. Both outputs are strictly identical, but the first method is faster and more memory-efficient because it doesn't load the entire file into R and avoids using extra processing memory. It should always be preferred when possible. Multiple filter commands can be used simultaneously to, for example, read only the first returns between 5 and 50 meters.


```r
las <-  readLAS("file.las", filter = "-keep_first -drop_z_below 5 -drop_z_above 50")
```

The full list of available commands can be obtained by using `readLAS(filter = "-help")`. Users of `LAStools` may recognize these commands, as both `LAStools` and `lidR` use the same libraries (`LASlib` and `LASzip`) to read and write LAS and LAZ files.

## Validating LiDAR Data {#sec-asprs-compliance}

An important first step in ALS data processing is ensuring that your data is complete and valid according to the [ASPRS LAS specifications](http://www.asprs.org/wp-content/uploads/2019/07/LAS_1_4_r15.pdf). Users commonly report bugs arising from invalid data. This is why we introduced the `las_check()` function to perform a thorough inspection of `LAS` objects. This function checks whether a `LAS` object meets the ASPRS LAS specifications and whether it is valid for processing, providing warnings if it does not.

A common issue is that a `LAS` file contains duplicate points. This can lead to problems such as trees being detected twice, invalid metrics, or errors in DTM generation. We may also encounter invalid return numbers, incoherent return numbers and number of returns attributes, and invalid coordinate reference systems, among other issues. Always make sure to run the `las_check()` function before delving deeply into your data.

```{r, echo = FALSE}
las$X[2] <- las$X[1]
las$Y[2] <- las$Y[1]
las$Z[2] <- las$Z[1]
las$Classification[1:2] <- 2L
las$ReturnNumber[3] <- 0L
```

```r
las_check(las)
#>  Checking the data
#>   - Checking coordinates... ✓
#>   - Checking coordinates type... ✓
#>   - Checking coordinates range... ✓
#>   - Checking coordinates quantization... ✓
#>   - Checking attributes type... ✓
#>   - Checking ReturnNumber validity...
#>     ⚠ Invalid data: 1 points with a return number equal to 0 found.
#>  [...]
```

A check is performed at read time regardless, but the read time check is not as thorough as `las_check()` for computation time reasons. For example duplicated points are not checked at read time.

```{r read-corrupted, warning = TRUE}
las <- readLAS("data/chap1/corrupted.laz")
```

## Plotting {#sec-plot}

```{r, echo = FALSE, warning = FALSE}
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile)
```

The `lidR` package takes advantage of the [`rgl`](https://cran.r-project.org/package=rgl) package to provide a versatile and interactive 3D viewer with points colored by Z coordinates on a black background as default.

### Basic 3D rendering {#sec-plot-3d}

The very basic way to render a point cloud is the function `plot()`.

``` r
plot(las)
```

```{r plot-las, echo = FALSE, rgl = TRUE, fig.width = 4, fig.height = 3}
#| label: plot-las
#| fig-cap: "Simple 3D interactive plot of a point cloud"
plot(las, size = 3)
```

Users can change the attributes used for coloring by providing the name of the attribute used to colorize the points. The background color of the viewer can also be changed by assigning a color using the `bg` argument. Axes can also be added and point sizes can be changed.

```{r plot-las-custom, rgl = TRUE, fig.width = 4, fig.height = 3}
#| label: plot-las-custom
#| fig-cap: "Simple 3D interactive plot of a point cloud colored by scan angle, with a white background and a color scale"

# Plot las object by scan angle, 
# make the background white, 
# display XYZ axis and  scale colors
plot(las, color = "ScanAngleRank", bg = "white", axis = TRUE, legend = TRUE)
```

Note that if your file contains RGB data the string `"RGB"` is supported:

``` r
plot(las, color = "RGB")
```

The argument `breaks` enables to defined more adequate breaks in the color palette for example when intensity contains large outliers. Otherwise the palette range would be too large and most of the values would be considered as "very low", so everything would appear in the same color.

```{r plot-las-fail, rgl = TRUE, fig.width = 4, fig.height = 3}
#| label: plot-las-fail
#| fig-cap: "Using quantiles to color the intensity provides a clear display despite outliers."
plot(las, color = "Intensity", breaks = "quantile", bg = "white")
```

### Overlays {#sec-plot-overlay}

The package also provides some easy to use functions for common overlay. For example `add_dtm3d()` to add a digital terrain model (section @sec-dtm)) and `add_treetops3d()` to visualize the output of an individual tree detection (section @sec-itd))

```{r, echo = FALSE, warning=FALSE}
dtm <- rasterize_terrain(las, 2, tin())
las <- clip_circle(las, 273516, 5274496, 100)
```

```{r plot-las-dtm, rgl = TRUE, fig.width = 4, fig.height = 3}
#| label: plot-las-dtm
#| fig-cap: "3D interactive rendering of a point cloud with a digital terrain model overlaid."
x <- plot(las, bg = "white", size = 3)
add_dtm3d(x, dtm)
```

```{r, echo = FALSE}
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzc")

ttops <- locate_trees(las, lmf(ws = 5))
```

```{r plot-las-trees, rgl = TRUE, fig.width = 4, fig.height = 3}
#| label: plot-las-trees
#| fig-cap: "3D interactive rendering of a point cloud with segmented tree tops overlaid."
x <- plot(las, bg = "white", size = 3)
add_treetops3d(x, ttops)
```

It is also possible to combine two point clouds with different color palettes. In the following example we are using a previously classified point cloud. We first separate the vegetation and non vegetation points using `filter_poi()` and then plot both on top of each other with different colour schemes using `add` options in `plot()`

```{r, echo = FALSE}
r3dDefaults$zoom = 0.3
las = readLAS("data/chap11/building_WilliamsAZ_Urban_normalized.laz", filter = "-thin_random_fraction 0.4")
```

```{r plot-las-add, rgl = TRUE, fig.width=8, fig.height=3}
#| label: plot-las-add
#| fig-cap: "3D interactive rendering of two point clouds overlaid with different color palettes."
nonveg <- filter_poi(las, Classification != LASHIGHVEGETATION)
veg <- filter_poi(las, Classification == LASHIGHVEGETATION)

x <- plot(nonveg, color = "Classification", bg = "white", size = 3)
plot(veg, add = x)
```

### Advanced 3D Rendering {#sec-plot-advanced}

Since `lidR` is based on `rgl`, it is easy to add objects to the main rendering using `rgl` functions such as `rgl::point3d()`, `rgl::text()`, `rgl::surface3d()`, and so on to produce publication-ready renderings. However, `lidR` introduces an additional challenge: it does not display the points with their actual coordinates. Instead, the points are shifted to be rendered close to (0, 0) due to accuracy issues, as `rgl` uses `float` (32-bit decimal numbers) rather than `double` (64-bit decimal numbers). When `plot()` is used, it invisibly returns the shift values, which can later be used to realign other objects.


```{r print-offset}
offsets <- plot(las)
print(offsets)
```

The coordinates of the objects must be corrected to align with the point cloud. In the following we will add lines to render the trunks. We read a file, we locate the trees (see @sec-itd)), we extract the coordinates and sizes of the trees and plot lines with `rgl::segment3d()`.

```{r, echo=FALSE}
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
```

```{r plot-las-truncks, rgl = TRUE, webgl = FALSE, fig.width = 4, fig.height = 3, snapshot = TRUE}
#| label: plot-las-trunks
#| code-fold: true
#| code-summary: "Show the code"
#| fig-cap: "3D interactive rendering of a point cloud with segmented tree tops and trunks overlaid."
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzc")

# get the location of the trees
ttops <- locate_trees(las, lmf(ws = 5)) 

# plot the point cloud
offsets <- plot(las, bg = "white", size = 3)
add_treetops3d(offsets, ttops)

# extract the coordinates of the trees and
# apply the shift to display the lines
# in the rendering coordinate system
x <- sf::st_coordinates(ttops)[,1] - offsets[1] 
y <- sf::st_coordinates(ttops)[,2] - offsets[2] 
z <- ttops$Z

# Build a GL_LINES matrix for fast rendering
x <- rep(x, each = 2)
y <- rep(y, each = 2)
tmp <- numeric(2*length(z)) 
tmp[2*1:length(z)] <- z
z <- tmp
M <- cbind(x,y,z)

# Display lines
rgl::segments3d(M, col = "black", lwd = 2)
```

### Voxel rendering

```{r, echo=F}
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile)
```

It is possible to render voxels. This is useful to render the output of the function `voxelise_points()` or `voxel_metrics()` for examples.

```{r plot-voxels, rgl = TRUE, fig.width = 4, fig.height = 3}
vox <- voxelize_points(las, 6)
plot(vox, voxel = TRUE, bg = "white")
```

### Cross Sections 2D Rendering {#sec-plot-crossection}

To better visualize the vertical structure of a point cloud, investigate classification results, or compare the results of different interpolation routines, a cross section can be plotted. To do this, we first need to decide where the cross section will be located (i.e., define the beginning and end) and specify its width. The point cloud can then be clipped, and the `X` and `Z` coordinates used to create the plot.

For example, to create a 200 m long cross section, we might define the beginning and end, and then use the `clip_transect()` function to subset the point cloud.


```{r, echo = FALSE}
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las <- readLAS(LASfile)
```

```{r}
p1 <- c(273457, 5274357)
p2 <- c(273542, 5274542)
las_tr <- clip_transect(las, p1, p2, width = 5, xz = TRUE)
```

Rendering can be achieved with base plot or `ggplot2`. Notice the use of `payload()` to extract the `data.frame` from the `LAS` object.

```{r ggplot-transect, fig.height=1.5, fig.width=8}
library(ggplot2)

ggplot(payload(las_tr), aes(X,Z, color = Z)) + 
  geom_point(size = 0.5) + 
  coord_equal() + 
  theme_minimal() +
  scale_color_gradientn(colours = height.colors(50))
```