From bccb7cd411b0396fbb15781de5425a96810f7eb9 Mon Sep 17 00:00:00 2001
From: jmmonnet <jean-matthieu.monnet@inrae.fr>
Date: Fri, 15 Jul 2022 13:58:37 +0200
Subject: [PATCH] Updated area-based mapping to include catalog_engine
---
R/ALS_data_preprocessing.Rmd | 6 +-
R/area-based.2.model.calibration.Rmd | 23 ++-
R/area-based.3.mapping.and.inference.Rmd | 228 +++++++++++------------
3 files changed, 130 insertions(+), 127 deletions(-)
diff --git a/R/ALS_data_preprocessing.Rmd b/R/ALS_data_preprocessing.Rmd
index ade9591..8c1cc65 100755
--- a/R/ALS_data_preprocessing.Rmd
+++ b/R/ALS_data_preprocessing.Rmd
@@ -108,7 +108,7 @@ point_cloud@header
```
### Basic attributes
-The `data` slot contains point attributes. Some attributes may be present or not or have different names depending on the format version of the LAS file. The main attributes (coordinates, intensity, classification, return number) are always present.
+The `data` slot contains point attributes. Some attributes may be present or not, or have different names depending on the format version of the LAS file. The main attributes (coordinates, intensity, classification, return number) are always present.
**`X`, `Y`, `Z`**: coordinates. The point cloud can be displayed with `lidR::plot`
```{r ALS.coordinates, eval=FALSE}
@@ -251,7 +251,7 @@ For checking purposes one might be interested to map the following features, in
The spatial resolution of the map should be relevant with respect to the announced specifications, and reach a trade-off between the amount of spatial detail and the possibility to evaluate the whole dataset at a glance.
-The `lidR::pixel_metrics` function is very efficient to derive statistics maps from the point cloud. Point statistics are computed for all pixels at the given resolution, based on points contained in the cells. In this case the density variation observed between the tiles is due to the fact that points from a given flight strip where written in the eastern tile but not in the western one.
+The `lidR::pixel_metrics` function is very efficient to derive statistics maps from the point cloud. Point statistics are computed for all pixels at the given resolution, based on points contained in the cells (the function `lidR::rasterize_density` is specifically designed to compute point density). In this example the density variation observed between the tiles is due to the fact that points from a given flight strip where written in the eastern tile but not in the western one.
```{r catalog.grid_metrics, include = TRUE, fig.width=9, fig.height = 5, warning = FALSE}
# resolution
@@ -296,7 +296,7 @@ terra::plot(metrics$npulses < 5, main = "Areas with pulse density < 5 /m2")
Normalization of a point cloud refers to the operation that consists in replacing the altitude coordinate by the height to the ground in the `Z` attribute. This is an important pre-requisite for the analysis of the vegetation 3D structure from the point cloud.
A pre-existing digital terrain model (DTM) can be used for normalizing, or it can be computed on-the-fly for this purpose. Obtaining relevant height values requires a DTM with sufficient resolution or the presence of enough ground points inside the data. To make sure that correct values are obtained at the border of the region of interest, it is advisable to add a buffer area before normalization.
-The next line normalizes the point cloud using a TIN constructed from the points classified as ground (class 2). The height values are stored in the `Z` attribute. The altitude value is copied into a new attribute called `Zref`. When the LAS object is written to a LAS file, this field is lost, unless the file version allows the writing of additional fields.
+The next line normalizes the point cloud using a TIN constructed from the points classified as ground (class 2). The height values are stored in the `Z` attribute. The altitude value is copied into a new attribute called `Zref`. When the LAS object is written to a LAS file, this field is lost, unless the file version allows the writing of additional fields.
```{r normalization.LAS, include = TRUE, fig.dim=c(8, 4)}
# normalize with tin algorithm for computed of surface from ground points
diff --git a/R/area-based.2.model.calibration.Rmd b/R/area-based.2.model.calibration.Rmd
index b6c484e..be517ae 100755
--- a/R/area-based.2.model.calibration.Rmd
+++ b/R/area-based.2.model.calibration.Rmd
@@ -86,6 +86,26 @@ Two types of vegetation metrics can be computed.
* Point cloud metrics are directly computed from the point cloud or from the derived surface model on the whole plot extent. These are the metrics generally used in the area-based approach.
+ Tree metrics are computed from the characteristics of trees detected in the point cloud (or in the derived surface model). They are more CPU-intensive to compute and require ALS data with higher density, but in some cases they allow a slight improvement in models prediction accuracy.
+## Removal of unwanted points
+
+The point cloud contains some points with classification code 7 (Low point) and 12 (Reserved). In this case the provider indicated that the 12 code also contains vegetation points: they are left in the point cloud whereas code 7 are removed. Points higher than a height threshold are also removed. Points with negative heights are set to zero.
+
+```{r removeLowpoints, include=TRUE}
+# define classes to retain
+class_points <- c(0, 1, 2, 3, 4, 5, 12)
+# height threshold (60 m)
+h_points <- 60
+# filter points
+llas_height <- lapply(llas_height, function(x) lidR::filter_poi(x, is.element(Classification, class_points) &
+ Z <= h_points))
+llas_height <- lapply(llas_height, function(x)
+{
+ x$Z[x$Z<0] <- 0
+ x
+}
+)
+```
+
## Point cloud metrics
Point cloud metrics are computed with the function `lidaRtRee::clouds_metrics`, which applies the function `lidR::cloud_metrics` to all point clouds in a list. Default computed metrics are those proposed by the function [`lidR::stdmetrics`](https://github.com/r-lidar/lidR/wiki/stdmetrics). Additional metrics are available with the function `lidaRtRee::aba_metrics`. The buffer points, which are located outside of the plot extent inventoried on the field, should be removed before computing those metrics.
@@ -324,7 +344,8 @@ The following lines save the data required for the [area-based mapping step](htt
```{r saveModels, eval=FALSE}
save(model_aba_stratified_mixed, model_aba, aba_point_metrics_fun, aba_res_chm,
- file = "../data/aba.model/output/models.rda"
+ class_points, h_points,
+ file = "../data/aba.model/output/models.rda"
)
```
diff --git a/R/area-based.3.mapping.and.inference.Rmd b/R/area-based.3.mapping.and.inference.Rmd
index 26e5e5d..e574148 100755
--- a/R/area-based.3.mapping.and.inference.Rmd
+++ b/R/area-based.3.mapping.and.inference.Rmd
@@ -65,19 +65,25 @@ out <- unzip(temp, exdir = "../data/aba.model/ALS/")
unlink(temp)
```
+Filters are directly specified in the catalog options, in order to load only relevant points. In particular, ground and vegetation points are loaded for height catalog, whereas only ground and water points are loaded for the altitude catalog. Points with heights larger than a threshold are alos filtered out.
+
```{r loadALS, include = TRUE, warning = FALSE}
# build catalogs
-cata_height <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.norm.laz/")
-cata_altitude <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.laz/")
-lidR::opt_progress(cata_altitude) <- lidR::opt_progress(cata_height) <- FALSE
+# option to read only xyzirnc attributes (coordinates, intensity, echo order and classification) from altitude files
+# build filter to keep certain classes and discard high points
+filter_cata_height <- paste(c("-keep_class", class_points, "-drop_z_above", h_points), collapse = " ")
+cata_height <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.norm.laz/",
+ progress = FALSE,
+ select = "xyzirnc",
+ filter = filter_cata_height)
+# option to read only xyzc attributes (coordinates and classification) from height files
+# option to read only ground and water points from altitude files
+cata_altitude <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.laz/",
+ progress = FALSE,
+ select = "xyzc",
+ filter = "-keep_class 2 9")
# set projection info
lidR::projection(cata_height) <- lidR::projection(cata_altitude) <- 2154
-# option to read only xyzc attributes (coordinates, intensity, echo order and classification) from height files
-lidR::opt_select(cata_height) <- "xyzirnc"
-# option to read only ground and water points from altitude files
-lidR::opt_filter(cata_altitude) <- "-keep_class 2 9"
-# option to read only xyzc attributes (coordinates and classification) from altitude files
-lidR::opt_select(cata_altitude) <- "xyzc"
cata_height
```
@@ -145,13 +151,13 @@ It is very important that tree segmentation parameters are the same as in the ca
# compute chm
chm <- lidR::rasterize_canopy(point_cloud_height, res = aba_res_chm, algorithm = lidR::p2r())
# replace NA, low and high values
-chm[is.na(chm) | chm < 0 | chm > 60] <- 0
+chm[is.na(chm) | chm < 0] <- 0
# tree top detection (same parameters as used by clouds_tree_metrics in calibration step -> default parameters)
segms <- lidaRtRee::tree_segmentation(chm)
# extraction to data.frame
-trees <- lidaRtRee::tree_extraction(segms$filled_dem, segms$local_maxima, segms$segments_id)
+trees <- lidaRtRee::tree_extraction(segms)
# compute raster metrics
-metrics_trees <- lidaRtRee::raster_metrics(trees[, -1],
+metrics_trees <- lidaRtRee::raster_metrics(trees,
res = resolution,
fun = function(x) {
lidaRtRee::std_tree_metrics(x, resolution^2 / 10000)
@@ -219,7 +225,7 @@ terra::plot(metrics_terrain$slope_gr, main = "Slope (gr)")
### Batch processing of tiles
-For the batch processing of tiles, parallel processing with packages `future` and `future.apply` is used. The processing capabilities of `lidR` are not yet used.
+For the batch processing of tiles, the catalog engine of `lidR` is used. It allows parallel processing with package `future`.
```{r setupCluster, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
# specify to use two parallel sessions
@@ -228,84 +234,87 @@ future::plan("multisession", workers = 2L)
options(future.rng.onMisuse = "ignore")
```
-A buffering procedure is designed to handle the border effects when detecting trees at tile edges. A 10 m buffer is enough for tree metrics. One can also specify points classes to use, or apply a threshold to high points. Some parameters specified in the previous paragraphs are also required for the batch processing (output map resolution, chm resolution for tree segmentation, functions for metrics computation).
+The catalog engine encompasses a buffering to handle the border effects when detecting trees at tile edges. A 10 m buffer is enough for tree metrics. One can also specify points classes to use, or apply a threshold to high points. Some parameters specified in the previous paragraphs are also required for the batch processing (output map `resolution`, chm resolution for tree segmentation `aba_res_chm`, functions for metrics computation `aba_point_metrics_fun`). Points classes to retain (`class_points`) and height threshold (`h_points`) have already been set as options for catalog processing.
```{r paramCluster, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
+rm(list = setdiff(ls(), c("resolution", "aba_res_chm", "aba_point_metrics_fun",
+ "cata_altitude", "cata_height", "class_points", "h_points",
+ "model_aba", "model_aba_stratified_mixed", "public", "forest")))
# processing parameters
b_size <- 10
-# height threshold (m) for high points removal (points above this threshold are considered
-# as outliers)
-h_points <- 60
-# points classes to retain for analysis (vegetation + ground)
+# point filter already in catalog options
+lidR::opt_filter(cata_height)
+# "-drop_z_above": height threshold (m) for high points removal
+# "-keep_class 0 1 2 3 4 5 12" points classes to retain for analysis (vegetation + ground)
# ground should be 2
-class_points <- c(0, 1, 2, 3, 4, 5)
+# here 12 code is also vegetation
```
This first chunk of code computes tree and point metrics from the normalized ALS tiles.
```{r batchProcessing, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
-# processing by tile: apply to all tiles a function that loads a tile plus buffer and then processes the data)
-metrics <- future.apply::future_lapply(
- as.list(1:nrow(cata_height)),
- FUN = function(i) {
- # tile extent
- b_box <- sf::st_bbox(cata_height[i,])
- # load tile extent plus buffer
- a <- try(lidR::clip_rectangle(cata_height,
- b_box[1] - b_size,
- b_box[2] - b_size,
- b_box[3] + b_size,
- b_box[4] + b_size))
- #
- # check if points are successfully loaded
- if (class(a) == "try-error") {
+# SET CATALOG ENGINE OPTIONS
+# process by file
+lidR::opt_chunk_size(cata_height) <- 0
+# buffer size
+lidR::opt_chunk_buffer(cata_height) <- b_size
+#
+# FUNCTION to apply to chunks
+routine <-
+ function(chunk,
+ res_map = resolution,
+ points_class = class_points,
+ points_h = h_points,
+ res_chm = aba_res_chm,
+ aba_fun = aba_point_metrics_fun) {
+ # extract bounding box (without buffer)
+ b_box <- sf::st_bbox(chunk)
+ # load chunk (including buffer)
+ a <- lidR::readLAS(chunk)
+ # return NULL if empty
+ if (lidR::is.empty(a))
return(NULL)
- }
# add 'buffer' flag to points in buffer with TRUE value in this new attribute
+ # the buffer flag is already present in chunk but it seems points lying on
+ # the northern and eastern borders are not considered "buffer"
a <- lidR::add_attribute(a,
a$X < b_box[1] |
a$Y < b_box[2] |
a$X >= b_box[3] | a$Y >= b_box[4],
"buffer")
- # remove unwanted point classes, and points higher than height threshold
- a <-
- lidR::filter_poi(a, is.element(Classification, class_points) &
- Z <= h_points)
+ # remove unwanted points: not necessary if specified in catalog options
+ # a <-
+ # lidR::filter_poi(a, is.element(Classification, points_class) &
+ # Z <= points_h)
# check number of remaining points
- if (nrow(a) == 0)
- return(NULL)
+ if (lidR::is.empty(a)) return(NULL)
# set negative heights to 0
a$Z[a$Z < 0] <- 0
- #
# compute chm
chm <-
lidR::rasterize_canopy(a, res = aba_res_chm, algorithm = lidR::p2r())
- # replace NA, low and high values by 0
- chm[is.na(chm) |
- chm < 0 | chm > h_points] <- 0
+ # replace NA by 0
+ # there should be no values below 0 and above h_points
+ chm[is.na(chm)] <- 0
#
# compute tree metrics
# tree top detection (default parameters)
segms <- lidaRtRee::tree_segmentation(chm)
# extraction to data.frame
- trees <-
- lidaRtRee::tree_extraction(segms$filled_dem, segms$local_maxima, segms$segments_id)
+ trees <- lidaRtRee::tree_extraction(segms)
+ # if no trees return null
+ if (is.null(trees)) return(NULL)
# remove trees outside of tile
- trees <-
- trees[trees$x >= b_box[1] &
- trees$x < b_box[3] &
- trees$y >= b_box[2] & trees$y < b_box[4],]
+ # no trees should be on a tile border so st_crop can be used
+ trees <- sf::st_crop(trees, sf::st_bbox(chunk))
# if no trees return null
- if (is.null(trees))
- return(NULL)
- if (nrow(trees) == 0)
- return(NULL)
- # compute raster metrics
+ if (nrow(trees) == 0) return(NULL)
+ # compute raster tree metrics
metrics_trees <- lidaRtRee::raster_metrics(
trees,
- res = resolution,
+ res = res_map,
fun = function(x) {
- lidaRtRee::std_tree_metrics(x, resolution ^ 2 / 10000)
+ lidaRtRee::std_tree_metrics(x, res_map ^ 2 / 10000)
},
output = "raster"
)
@@ -314,12 +323,14 @@ metrics <- future.apply::future_lapply(
r_treechm <- segms$filled_dem
# set chm to NA in non segment area
r_treechm[segms$segments_id == 0] <- NA
- # compute raster metrics
+ # crop to bbox
+ r_treechm <- terra::crop(r_treechm, terra::ext(chunk))
+ # compute raster metrics in bbox
metrics_trees2 <- lidaRtRee::raster_metrics(
- terra::crop(r_treechm, terra::ext(b_box[1], b_box[3], b_box[2], b_box[4])),
- res = resolution,
+ r_treechm,
+ res = res_map,
fun = function(x) {
- c(sum(!is.na(x$filled_dem)) / (resolution / aba_res_chm) ^ 2,
+ c(sum(!is.na(x$filled_dem)) / (res_map / aba_res_chm) ^ 2,
mean(x$filled_dem, na.rm = T))
},
output = "raster"
@@ -329,78 +340,50 @@ metrics <- future.apply::future_lapply(
#
# compute 1D height metrics
# remove buffer points
- a <- lidR::filter_poi(a, buffer == FALSE)
+ a <- lidR::filter_poi(a, buffer == 0)
#
- if (nrow(a) == 0)
- return(NULL)
+ if (lidR::is.empty(a)) return(NULL)
# all points metrics
- metrics_points <-
- lidR::pixel_metrics(a, aba_point_metrics_fun, res = resolution)
+ metrics_points <- lidR::pixel_metrics(a, aba_fun, res = res_map)
#
- # extend / crop to match metrics_points
- metrics_trees <-
- terra::extend(metrics_trees, metrics_points)
- metrics_trees2 <-
- terra::extend(metrics_trees2, metrics_points)
- metrics_trees <-
- terra::crop(metrics_trees, metrics_points)
- metrics_trees2 <-
- terra::crop(metrics_trees2, metrics_points)
- # merge rasterstacks
- metrics <- metrics_points
- metrics <-
- c(metrics, metrics_trees)
- metrics <-
- c(metrics, metrics_trees2)
- if (is.null(metrics))
- {
- return(NULL)
- } else {
- return(terra::wrap(metrics))
- }
+ if (is.null(metrics_points)) return(NULL)
+ # extend / crop to match bbox
+ b_box <- terra::ext(chunk)
+ metrics_trees <- terra::extend(metrics_trees, b_box)
+ metrics_trees2 <- terra::extend(metrics_trees2, b_box)
+ metrics_points <- terra::extend(metrics_points, b_box)
+ metrics_trees <- terra::crop(metrics_trees, b_box)
+ metrics_trees2 <- terra::crop(metrics_trees2, b_box)
+ metrics_points <- terra::crop(metrics_points, b_box)
+ # merge rasters
+ metrics <- c(metrics_points, metrics_trees, metrics_trees2)
+ #
+ return(metrics)
}
-)
-# remove NULL rasters
-metrics <- metrics[which(!sapply(metrics, is.null))]
-# unwrap rasters
-metrics <- lapply(metrics, terra::rast)
-# merge rasters
-# error might occur if one raster has only one cell (check if still relevant in terra package)
-metrics_map <- do.call(terra::merge, metrics)
+#
+# APPLY TO CATALOG, specify to merge results
+metrics_map <- lidR::catalog_apply(cata_height, routine, .options = list(automerge = TRUE))
```
-
The second chunk of code computes terrain metrics from the ALS tiles with altitude values.
```{r batchProcessingTerrain, include=TRUE, message=FALSE, warning=FALSE}
-# function to compute terrain metrics to input to pixel_metrics
+# SET CATALOG ENGINE OPTIONS
+# process by file
+lidR::opt_chunk_size(cata_altitude) <- 0
+# no buffer additional buffer
+lidR::opt_chunk_buffer(cata_altitude) <- 0
+# enable auto-merge of results
+lidR::opt_merge(cata_altitude) <- TRUE
+#
+# FUNCTION to compute terrain metrics to input to pixel_metrics
f <- function(x, y, z) {
as.list(lidaRtRee::terrain_points_metrics(data.frame(x, y, z)))
}
-# apply function to all tiles specified by indices
-metrics_terrain <- future.apply::future_lapply(
- as.list(1:nrow(cata_height)),
- FUN = function(i) {
- # tile extent
- b_box <- sf::st_bbox(cata_altitude[i,])
- # load tile extent
- a <- try(lidR::clip_rectangle(
- cata_altitude, b_box[1], b_box[2], b_box[3],
- b_box[4]
- ))
- # check if points are successfully loaded
- if (class(a) == "try-error") {
- return(NULL)
- }
- terra::wrap(lidR::pixel_metrics(a, ~ f(X, Y, Z), res = resolution))
- }
-)
-# remove NULL rasters
-metrics_terrain <- metrics_terrain[which(!sapply(metrics_terrain, is.null))]
-# unwrap rasters
-metrics_terrain <- lapply(metrics_terrain, terra::rast)
-# mosaic rasters
-metrics_terrain <- do.call(terra::merge, metrics_terrain)
+# APPLY TO CATALOG
+metrics_terrain <- lidR::pixel_metrics(cata_altitude, ~ f(X, Y, Z), res = resolution)
+# remove unused layer
+metrics_terrain <- terra::subset(metrics_terrain, which(names(metrics_terrain) != "adjR2_plane"))
```
Finally all metrics are combined in a single object.
@@ -408,7 +391,6 @@ Finally all metrics are combined in a single object.
```{r allMetrics, include=TRUE, message=FALSE, warning=FALSE}
metrics_terrain <- terra::extend(metrics_terrain, metrics_map)
metrics_terrain <- terra::crop(metrics_terrain, metrics_map)
-metrics_terrain <- terra::subset(metrics_terrain, which(names(metrics_terrain) != "adjR2_plane"))
# merge rasters
metrics_map <- c(metrics_map, metrics_terrain)
```
--
GitLab