02-io-display.Rmd

```{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)
old_hooks <- fansi::set_knit_hooks(knitr::knit_hooks, which = c("output", "message", "error"))
```

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

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

Discrete return ALS sensors record a number of pieces of data. First and foremost, 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, or the beam incidence angle of each point. Reading, writing, and efficient storage of these ALS data is a critical step prior to any subsequent analysis. 

ALS data is 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 [American Society for Photogrammetry & Remote Sensing (ASPRS)](http://www.asprs.org/wp-content/uploads/2019/07/LAS_1_4_r15.pdf). LAS files do however 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, which has been designed to process LAS and LAZ files both as input and output, taking advantage of 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 depth in a [dedicated vignette](https://cran.r-project.org/web/packages/lidR/vignettes/lidR-LAS-class.html). To briefly summarize, a LAS file is made of two parts:

1. The header that 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` {#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. We call these data *attributes*. In pratice many attributes are not actually useful but they are loaded anyway by default. This can take up a lot of processing memory because R is a language that 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` {#filter}

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

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

It is important to understand that the option `filter` in `readLAS()` keeps or discards point **at read time** i.e. while reading at the C++ level without implying any R code. For example the R function to filter points of interest (POI) is `filter_poi()` may return the exact 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 above example we are (1) reading only the first returns or (2) Reading all the points then filtering the first returns in R. Both outputs are strictly identical but the first one is faster and more memory efficient because it doesn't load the whole file in R and does not use extra processing memory. It should always be preferred when possible. Multiple filter commands can be used at once to e.g. read only first return between 5 and 50 m.

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

The full list of available commands is given by `readLAS(filter = "-help")`. Users of `LAStools` may recognize these commands because both `LAStools` and `lidR` use the same library (`LASlib` and `LASzip`) to read and write LAS and LAZ files.

## Validating lidar data {#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 deep inspection of `LAS` objects. This function checks if a `LAS` object meets the ASPRS LAS specifications and whether it is valid for processing, giving warnings if otherwise.

A simple example that happens fairly often is that a `LAS` file contains duplicate points. This may lead to problems like trees being detected twice, to invalid metrics, or to errors in DTM generation. We can also encounter invalid return numbers, incoherent return numbers and number of returns attributes, invalid coordinate reference system etc. Always make sure to run the `las_check()` function before digging deep 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 check-las}
las_check(las)
```

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 {#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 coloured by Z coordinates on a black background as default.

### Basic 3D rendering {#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}
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}
# 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}
plot(las, color = "Intensity", breaks = "quantile", bg = "white")
```

### Overlays {#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 \@ref(dtm)) and `add_treetops3d()` to visualize the output of an individual tree detection (section \@ref(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}
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}
x <- plot(las, bg = "white", size = 3)
add_treetops3d(x, ttops)
```

It is also possible to combine two point clouds with different colour 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}
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 {#plot-advanced}

With `lidR` being based on `rgl` it is easy to add objects in the main rendering using `rgl` functions such as `rgl::point3d()`, `rgl::text()`, `rgl::surface3d()` and so on to produce publication ready rendering.  However `lidR` introduced an additional challenge because it does not display the points with their actual coordinates. The points are shifted to be rendered close to (0, 0) (a matter of accuracy because `rgl` uses `float` (decimal numbers on 32 bits) instead of `double` (decimal numbers on 64 bits)). When `plot()` is used it invisibly returns the shift values that can be used later 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 section \@ref(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}
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 {#plot-crossection}

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

For example, to create a 100 m long cross section we may define the beginning and the end and then use `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(273357, 5274357)
p2 <- c(273542, 5274542)
las_tr <- clip_transect(las, p1, p2, width = 4, xz = TRUE)
```

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

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

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

The two steps required to create a cross section (clipping the point cloud and plotting) can be combined. Below we create a simple function that will become handy at multiple occasions throughout this book. To make this function even easier to use we will specify the default values for `p1` and `p2` so that the cross section is located in the centre of the point cloud, along the X-axis. The default width will be 4 m.

```{r, code=readLines("function_plot_crossection.R")}

```

Then we can used the function:

```{r ggplot-transect-2, fig.height=1.5, fig.width=8}
plot_crossection(las, colour_by = factor(Classification))
```