---
title: "Pipelines"
author: "Jean-Romain Roussel"
output:
  html_document:
    toc: true
    toc_float:
      collapsed: false
      smooth_scroll: false
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{3. Pipelines}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup, echo=F}
suppressPackageStartupMessages(library(lasR))
```

After reading the [tutorial](tutorial.html), one may have noticed that writing certain pipelines can be tedious, especially when normalizing a dataset, which always implies using the same code.

```{r}
tri <- triangulate(filter = keep_ground())
trans <- transform_with(tri)
norm <- tri + trans
```

This vignette shows a list of pipeline factories, i.e., functions that generate useful pieces of pipelines that can be reused. It also defines some functions that execute a pipeline, which may be regularly useful instead of writing it out. The lists are not comprehensive and will be updated if good tools come to mind. Some are already installed natively in the package.

## Pipeline factories

### Normalize

<blockquote style="background-color: #d6e9f9; border-left: 5px solid #428bca; padding: 10px;font-size: 14px; border-radius: 5px;">
This pipeline is already installed in the package.
</blockquote>

```{r}
normalize = function(extrabytes = FALSE)
{
  tri <- triangulate(filter = keep_ground())
  pipeline <- tri
  
  if (extrabytes)
  {
    extra = add_extrabytes("int", "HAG", "Height Above Ground")
    trans = transform_with(tri, store_in_attribute = "HAG")
    pipeline = pipeline + extra + trans
  }
  else
  {
    trans = transform_with(tri)
    pipeline = pipeline + trans
  }
  
  return(pipeline)
}
```

It can be used this way

```{r, eval=FALSE}
pipeline = reader(f) + normalize() + write_las(o)
```

### Classify ground

#### With CSF

```r
ground_csf = function(smooth = FALSE, threshold = 0.5, resolution = 0.5, rigidness = 1L, iterations = 500L, step = 0.65)
{
  csf = function(data, smooth, threshold, resolution, rigidness, iterations, step)
  {
    id = RCSF::CSF(data, smooth, threshold, resolution, rigidness, iterations, step)
    class = integer(nrow(data))
    class[id] = 2L
    data$Classification <- class
    return(data)
  }
  
  classify = callback(csf, expose = "xyz", smooth = smooth, threshold = threshold, resolution = resolution, rigidness = rigidness, iterations = iterations, step = step)
  return(classify)
}
```

```{r, eval=FALSE}
pipeline = reader(f) + ground_csf() + write_las(o)
```

#### With MCC

```r
ground_mcc = function(s = 1.5, t = 0.3)
{
  csf = function(data, s, t)
  {
    id = RMCC::MCC(data, s, t)
    class = integer(nrow(data))
    class[id] = 2L
    data$Classification <- class
    return(data)
  }
  
  classify = callback(csf, expose = "xyz", s = s, t = t)
  return(classify)
}
```

```{r, eval=FALSE}
pipeline = reader(f) + ground_mcc() + write_las(o)
```

<blockquote style="background-color: #fff3cd; border-left: 5px solid #ffc107; padding: 10px; font-size: 14px; border-radius: 5px;">
These pipelines use `callback()` that exposes the point cloud as a `data.frame`. One of the reasons why `lasR` is more memory-efficient and faster than `lidR` is that it **does not** expose the point cloud as a `data.frame`. Thus, these pipelines are not very different from the `classify_ground()` function in `lidR`. The advantage of using `lasR` here is the ability to pipe different stages.
</blockquote>

### Canopy Heigh Model

<blockquote style="background-color: #d6e9f9; border-left: 5px solid #428bca; padding: 10px;font-size: 14px; border-radius: 5px;">
These two pipelines are natively installed in the package under the name `chm()`.
</blockquote>

```{r}
chm_p2r = function(res, filter = "", ofile = tempfile(fileext = ".tif"))
{
  return(rasterize(res, "max", filter = filter, ofile = ofile))
}
```

```{r}
chm_tin = function(res, ofile = tempfile(fileext = ".tif"))
{
  tin = triangulate(filter = keep_first())
  chm = rasterize(res, tin, ofile = ofile)
  return(tin+chm)
}
```

### Digital Terrain Model 

<blockquote style="background-color: #d6e9f9; border-left: 5px solid #428bca; padding: 10px;font-size: 14px; border-radius: 5px;">
This one is also natively installed in the package.
</blockquote>

`add_class` can be used to add a class used as ground such as 9 for water. 

```{r}
dtm = function(res, ofile = tempfile(fileext = ".tif"), add_class = NULL)
{
  filter = keep_ground()
  if (!is.null(add_class)) filter = filter + keep_class(add_class)
  tin = triangulate(filter = filter)
  chm = rasterize(res, tin, ofile = ofile)
  return(tin+chm)
}
```

## Useful functions

### Read LAS

```{r}
read_las = function(files, select = "xyzi", filter = "")
{
  load = function(data) { return(data) }
  read = reader_las(filter = filter)
  call = callback(load, expose = select, no_las_update = TRUE)
  return(exec(read+call, on = f))
}
```

### Buffer tiles

```{r}
buffer_tiles = function(files, buffer, ofiles = paste0(tempdir(), "/*_buffered.las"))
{
  read = reader_las()
  write = write_las(ofiles, keep_buffer = TRUE)
  return(exec(read+write, on = files, buffer = buffer))
}
```

### Clip circle

Writes LAS files or returns `data.frame`s. Supports `sf` objects as input.

```{r}
clip_circle = function(files, geometry, radius, ofiles = paste0(tempdir(), "/*_clipped.las"))
{
  if (sf::st_geometry_type(geometry, FALSE) != "POINT") 
    stop("Expected POINT geometry type")

  coordinates <- sf::st_coordinates(geometry)
  xcenter <- coordinates[,1]
  ycenter <- coordinates[,2]
  
  read = reader_las(xc = xcenter, yc = ycenter, r = radius)
  
  if (length(ofiles) == 1L && ofiles == "")
    stage = callback(function(data) { return(data) }, expose = "*", no_las_update = T)
  else
    stage = write_las(ofiles)
  
  ans = exec(read+stage, on = files)
  return(ans)
}
```

### CRS

A CRS as `sf` object. The cost of applying `hulls()` is virtually null.

```{r crs}
crs = function(files)
{
  pipeline = reader_las() + hulls()
  ans = exec(pipeline, on = files)
  return(sf::st_crs(ans))
}
```

```{r, echo=FALSE, results='hide'}
f = paste0(system.file(package="lasR"), "/extdata/bcts/")
f = list.files(f, pattern = "(?i)\\.la(s|z)$", full.names = TRUE)
crs(f)
```

### Inventory metrics

Using an `sf` object to provide plot centers and offering the option to normalize on-the-fly.  It returns the same `sf` object with extra attributes.

```{r}
inventory_metrics = function(files, geometry, radius, fun, normalize = FALSE)
{
  if (sf::st_geometry_type(geometry, FALSE) != "POINT") 
    stop("Expected POINT geometry type")

  coordinates <- sf::st_coordinates(geometry)
  xcenter <- coordinates[,1]
  ycenter <- coordinates[,2]
  
  pipeline <- reader_las(xc = xcenter, yc = ycenter, r = radius)
  
  if (normalize)
  {
    tri <- triangulate(filter = keep_ground())
    trans <- transform_with(tri)
    pipeline <- pipeline + tri + trans
  }

  pipeline <- pipeline + callback(fun, expose = "*")
  ans <- exec(pipeline, on = files)
  ans <- lapply(ans, as.data.frame)
  ans <- do.call(rbind, ans)
  return(cbind(geometry, ans))
}
```

```{r, echo=FALSE, results='hide'}
f = paste0(system.file(package="lasR"), "/extdata/bcts/")
f = list.files(f, pattern = "(?i)\\.la(s|z)$", full.names = TRUE)
f = f[1:2]

x = c(885100, 885100)
y = c(629200, 629600)
df = data.frame(x,y)
geometry = sf::st_as_sf(df, coords = c("x", "y"))
radius = 20
fun = function(data) { list(metric1 = median(data$Z), metric2 = mean(data$Intensity)) }

inventory_metrics(f, geometry, radius, fun)
```

### Virtual point cloud

```{r}
build_vpc = function(files, ofile)
{
  read = reader_las()
  write = write_vpc(ofile)
  exec(read+write, on = files)
}
```

```{r, echo=FALSE, results='hide'}
f = paste0(system.file(package="lasR"), "/extdata/bcts/")
build_vpc(f, tempfile(fileext = ".vpc"))
```