Skip to content
GitLab
Commit eef15d73 authored by Monnet Jean-Matthieu's avatar Monnet Jean-Matthieu
Browse files

future_lapply used in aba3

parent 601fe438
master dev
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 222 additions and 205 deletions
+222 -205
R/area-based.3.mapping.and.inference.Rmd
+ 132
− 109
View file @ eef15d73
......@@ -21,7 +21,7 @@ knitr::opts_chunk$set(fig.align = "center")
```
---
The code below presents a workflow to map forest parameters using calibrated prediction models (see [previous tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.2.model.calibration.Rmd)) and wall-to-wall coverage of the region of interest with airborne laser scanning. This tutorial is the last step of the so-called area-based approach. It is based on functions from `R` packages `lidaRtRee` and `lidR`.
The code below presents a workflow to map forest parameters using calibrated prediction models (see [previous tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.2.model.calibration.Rmd)) and wall-to-wall coverage of the region of interest with airborne laser scanning. This tutorial is the last step of the so-called area-based approach. It is based on functions from `R` packages [lidaRtRee](https://cran.r-project.org/package=lidaRtRee) (`r packageVersion("lidaRtRee")`) and [lidR](https://cran.r-project.org/package=lidR) (`r packageVersion("lidR")`).
Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.3.mapping.and.inference.Rmd)
......@@ -212,13 +212,13 @@ raster::plot(metrics_terrain$slope_gr, main = "Slope (gr)")
### Batch processing of tiles
For the batch processing of tiles, parallel processing with packages `foreach` and `doFuture` is used. The processing capabilities of `lidR` are not yet used.
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.
```{r setupCluster, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
library(foreach)
# create parallel frontend, specify to use two parallel sessions
doFuture::registerDoFuture()
# specify to use two parallel sessions
future::plan("multisession", workers = 2L)
# remove warning when using random numbers in parallel sessions
options(future.rng.onMisue = "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).
......@@ -237,95 +237,115 @@ class_points <- c(0, 1, 2, 3, 4, 5)
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
metrics <- foreach::foreach(i = 1:nrow(cata_height), .errorhandling = "remove") %dopar% {
# tile extent
b_box <- as.numeric(cata_height@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# 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") {
return(NULL)
}
# add 'buffer' flag to points in buffer with TRUE value in this new attribute
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)
# check number of remaining points
if (length(a) == 0) {
return(NULL)
}
# set negative heights to 0
a@data$Z[a@data$Z < 0] <- 0
#
# compute chm
chm <- lidaRtRee::points2DSM(a, res = aba_res_chm)
# replace NA, low and high values by 0
chm[is.na(chm) | chm < 0 | chm > h_points] <- 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)
# 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], ]
# compute raster metrics
metrics_trees <- lidaRtRee::raster_metrics(trees[, -1],
res = resolution,
fun = function(x) {
lidaRtRee::std_tree_metrics(x, resolution^2 / 10000)
},
output = "raster"
)
# compute canopy cover in trees and canopy mean height in trees
# in region of interest, because it is not in previous step.
r_treechm <- segms$filled_dem
# set chm to NA in non segment area
r_treechm[segms$segments_id == 0] <- NA
# compute raster metrics
metrics_trees2 <- lidaRtRee::raster_metrics(
raster::crop(r_treechm, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4])),
res = resolution,
fun = function(x) {
c(
sum(!is.na(x$filled_dem)) / (resolution / aba_res_chm)^2,
mean(x$filled_dem, na.rm = T)
)
},
output = "raster"
)
names(metrics_trees2) <- c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
#
# compute 1D height metrics
# remove buffer points
a <- lidR::filter_poi(a, buffer == 0)
#
if (length(a) == 0) {
return(NULL)
# 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 <-
as.numeric(cata_height@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# 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") {
return(NULL)
}
# add 'buffer' flag to points in buffer with TRUE value in this new attribute
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)
# check number of remaining points
if (length(a) == 0) {
return(NULL)
}
# set negative heights to 0
a@data$Z[a@data$Z < 0] <- 0
#
# compute chm
chm <-
lidaRtRee::points2DSM(a, res = aba_res_chm)
# replace NA, low and high values by 0
chm[is.na(chm) |
chm < 0 | chm > h_points] <- 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)
# 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],]
# compute raster metrics
metrics_trees <- lidaRtRee::raster_metrics(
trees[,-1],
res = resolution,
fun = function(x) {
lidaRtRee::std_tree_metrics(x, resolution ^ 2 / 10000)
},
output = "raster"
)
# compute canopy cover in trees and canopy mean height in trees
# in region of interest, because it is not in previous step.
r_treechm <- segms$filled_dem
# set chm to NA in non segment area
r_treechm[segms$segments_id == 0] <- NA
# compute raster metrics
metrics_trees2 <- lidaRtRee::raster_metrics(
raster::crop(r_treechm, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4])),
res = resolution,
fun = function(x) {
c(sum(!is.na(x$filled_dem)) / (resolution / aba_res_chm) ^ 2,
mean(x$filled_dem, na.rm = T))
},
output = "raster"
)
names(metrics_trees2) <-
c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
#
# compute 1D height metrics
# remove buffer points
a <- lidR::filter_poi(a, buffer == 0)
#
if (length(a) == 0) {
return(NULL)
}
# all points metrics
metrics_points <-
lidR::grid_metrics(a, aba_point_metrics_fun, res = resolution)
#
# extend / crop to match metrics_points
metrics_trees <-
raster::extend(metrics_trees, metrics_points, values = NA)
metrics_trees2 <-
raster::extend(metrics_trees2, metrics_points, values = NA)
metrics_trees <-
raster::crop(metrics_trees, metrics_points)
metrics_trees2 <-
raster::crop(metrics_trees2, metrics_points)
# merge rasterstacks
metrics <- metrics_points
metrics <-
raster::addLayer(metrics, metrics_trees)
metrics <-
raster::addLayer(metrics, metrics_trees2)
return(metrics)
}
# all points metrics
metrics_points <- lidR::grid_metrics(a, aba_point_metrics_fun, res = resolution)
#
# extend / crop to match metrics_points
metrics_trees <- raster::extend(metrics_trees, metrics_points, values = NA)
metrics_trees2 <- raster::extend(metrics_trees2, metrics_points, values = NA)
metrics_trees <- raster::crop(metrics_trees, metrics_points)
metrics_trees2 <- raster::crop(metrics_trees2, metrics_points)
# merge rasterstacks
metrics <- metrics_points
metrics <- raster::addLayer(metrics, metrics_trees)
metrics <- raster::addLayer(metrics, metrics_trees2)
return(metrics)
}
)
# mosaic rasters
names_metrics <- names(metrics[[1]])
metrics_map <- do.call(raster::merge, metrics)
......@@ -340,23 +360,26 @@ The second chunk of code computes terrain metrics from the ALS tiles with altitu
f <- function(x, y, z) {
as.list(lidaRtRee::terrain_points_metrics(data.frame(x, y, z)))
}
#
metrics_terrain <- foreach::foreach(i = 1:nrow(cata_altitude), .errorhandling = "remove") %dopar% {
# tile extent
b_box <- as.numeric(cata_altitude@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# load tile extent plus buffer
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)
# 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 <- as.numeric(cata_altitude@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# load tile extent plus buffer
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)
}
# keep only ground points
a <- lidR::filter_ground(a)
lidR::grid_metrics(a, ~ f(X, Y, Z), res = resolution)
}
# keep only ground points
a <- lidR::filter_ground(a)
lidR::grid_metrics(a, ~ f(X, Y, Z), res = resolution)
}
)
# mosaic rasters
names_metrics <- names(metrics_terrain[[1]])
metrics_terrain <- do.call(raster::merge, metrics_terrain)
......
export/area-based.3.mapping.and.inference.html
+ 90
− 96
View file @ eef15d73
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment