diff --git a/R/ALS_data_preprocessing.Rmd b/R/ALS_data_preprocessing.Rmd
index 36a770470830394b249dff9a1769a65b730c1c2a..ade9591a2acb22a7336e414ec2dbab84d02819bc 100755
--- a/R/ALS_data_preprocessing.Rmd
+++ b/R/ALS_data_preprocessing.Rmd
@@ -84,9 +84,11 @@ A single LAS file can be loaded in R with the `lidR::readLAS` function, which re
* `data`: point cloud attributes.
* `crs`: projection information.
-The bounding box of the data cona be accessed with the function `sf::st_bbox`.
+The bounding box of the data can be accessed with the function `sf::st_bbox`.
-It is advisable to fill the projection information if missing, as spatial objects computed from the point cloud will inherit the information. A warning is issued when loading this file because unexpected values are present in the scan angle field.
+It is advisable to fill the projection information if missing, as spatial objects computed from the point cloud will inherit the information.
+
+A warning is issued when loading the following example file because unexpected values are present in the scan angle field.
```{r ALS.load, include = TRUE}
# read file
@@ -144,16 +146,16 @@ If one considers that the power transmitted by the emitter $P_t$ is constant and
In case the aircraft trajectory is available, it is possible to normalize amplitude values from the range ($R$) effect. One must however be aware that the resulting distribution of amplitude values might not be homogeneous because in areas with longer range the signal to noise ratio is lower.
-```{r ALS.intensity, eval=TRUE, fig.width=6}
+```{r ALS.intensity, eval=TRUE, fig.width=4}
hist(point_cloud$Intensity, xlab = "Raw value", main = "Histogram of Intensity")
```
**Echo position within returns from the same pulse**: `ReturnNumber` is the order of arrival of the echo among the `NumberOfReturns` echoes associated to a given pulse. Echoes are sometimes referred to as:
-* single if `NumberOfReturns = ReturnNumber = 1`
-* first if `ReturnNumber = 1`
-* last if `NumberOfReturns = ReturnNumber`
-* intermediate if `1 < ReturnNumber < NumberOfReturns`
+* *single* if `NumberOfReturns = ReturnNumber = 1`
+* *first* if `ReturnNumber = 1`
+* *last* if `NumberOfReturns = ReturnNumber`
+* *intermediate* if `1 < ReturnNumber < NumberOfReturns`
The contingency table of `ReturnNumber` and `NumberOfReturns` should have no values below the diagonal, and approximately the same values in each column.
@@ -190,7 +192,7 @@ table(point_cloud$Classification)
**`ScanAngleRank`** is the scan angle associated to the pulse. Origin of values might correspond to the horizontal or vertical. In most cases it is the scan angle relative to the laser scanner, but sometimes values relative to nadir are indicated. In the case of this file values are from 60 to 120: the scan range is ± 30 degrees on both sides of the scanner, with 90 the value when pulses are emitted downwards.
-```{r ALS.scanangle, eval=TRUE, fig.width=6}
+```{r ALS.scanangle, eval=TRUE, fig.width=4}
hist(point_cloud$ScanAngleRank, main = "Histogram of scan angles", xlab = "Angle (degrees)")
```
**`UserData`** is a field to store additional information, but it is limited in capacity. In the latest version of LAS files, additional attributes can be added.
@@ -249,9 +251,9 @@ 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::grid_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. 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.
-```{r catalog.grid_metrics, include = TRUE, fig.width=9, fig.height = 5.5, warning = FALSE}
+```{r catalog.grid_metrics, include = TRUE, fig.width=9, fig.height = 5, warning = FALSE}
# resolution
res <- 5
# build function to compute desired statistics
@@ -277,9 +279,12 @@ p_pulse_5 <- round(sum(terra::values(metrics$npulses)<10) / length(terra::values
From those maps summary statistics can be derived, e.g. the percentage of cells with no ground points (`r p_no_ground`). Checking that the spatial distribution of pulses is homogeneous can be informative.
-```{r catalog.maps, fig.width = 9, fig.height = 4}
-par(mfrow = c(1,2))
+```{r catalog.maps, fig.width = 9, fig.height = 3}
+par(mfrow = c(1,1))
+layout(matrix(c(1,2,2), nrow = 1, ncol = 3, byrow = TRUE))
+# par(fig=c(0,3,0,10)/10)
hist(terra::values(metrics$npulses), main = "Histogram of pulse density", xlab = "Pulse density /m2", ylab = "Number of 5m cells")
+# par(fig=c(3,10,0,10)/10)
terra::plot(metrics$npulses < 5, main = "Areas with pulse density < 5 /m2")
```
@@ -310,7 +315,7 @@ Digital elevation models (DEMs) are raster images where each cell value correspo
### Digital terrain model (DTM)
-It represents the altitude of the bare earth surface. It is computed using points classified as ground. For better representation of certain topographical features, some additional points or lines might be added. The function `lidR::grid_terrain` proposes different algorithms for computation, and works either with catalogs or LAS objects.
+It represents the altitude of the bare earth surface. It is computed using points classified as ground. For better representation of certain topographical features, some additional points or lines might be added. The function `lidR::rasterize_terrain` proposes different algorithms for computation, and works either with catalogs or LAS objects.
```{r dem.terrain, include = TRUE, message = FALSE}
# create dtm at 0.5 m resolution with TIN algorithm
@@ -320,7 +325,7 @@ dtm
Functions from the `terra` package are available to derive topographical rasters (slope, aspect, hillshade, contour lines) from the DTM.
-```{r dem.terrain.derived, include = TRUE, message = FALSE, fig.dim = c(6, 4)}
+```{r dem.terrain.derived, include = TRUE, message = FALSE, fig.width = 6.3, fig.height = 3.2}
dtm_slope <- terra::terrain(dtm, unit = "radians")
dtm_aspect <- terra::terrain(dtm, v = "aspect", unit = "radians")
dtm_hillshade <- terra::shade(dtm_slope, dtm_aspect)
@@ -356,7 +361,7 @@ chm_norm
If high points are present in the point cloud, a threshold can be applied to the CHM values to remove outliers.
-```{r dem.chm.threshold, include = TRUE, message = FALSE, fig.dim = c(6, 4)}
+```{r dem.chm.threshold, include = TRUE, message = FALSE, fig.width = 6.3, fig.height = 3.2}
threshold <- 40
chm[chm > threshold] <- threshold
terra::plot(chm, main = "Canopy Height Model (m)")
diff --git a/R/area-based.1.data.preparation.Rmd b/R/area-based.1.data.preparation.Rmd
index f8369857bcd2b1c0c79bc31bd3c2a4647b3c05c2..2dc165b88e4c35691d7702c9bf613359fa400d89 100755
--- a/R/area-based.1.data.preparation.Rmd
+++ b/R/area-based.1.data.preparation.Rmd
@@ -26,7 +26,7 @@ The code below presents a workflow to prepare inventory data for the calibration
Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.1.data.preparation.Rmd)
-Required `R` packages : `ggplot2`, `sf`, `ggmap`
+Required `R` packages : `ggplot2`, `sf`, `ggmap`, `lidaRtRee` (tested with version `r packageVersion("lidaRtRee")`) and `lidR` (tested with version `r packageVersion("lidR")`)
Many thanks to Pascal Obstétar for checking code and improvement suggestions.
diff --git a/R/area-based.2.model.calibration.Rmd b/R/area-based.2.model.calibration.Rmd
index 816eb5789b6a905a5f44bdd045324dcd82622a63..b6c484e4426049f2dd49e2b8bcaccdfd78c27e0c 100755
--- a/R/area-based.2.model.calibration.Rmd
+++ b/R/area-based.2.model.calibration.Rmd
@@ -88,7 +88,7 @@ Two types of vegetation metrics can be computed.
## 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/Jean-Romain/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.
+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.
```{r computeMetrics, include=TRUE}
# define function for later use
@@ -159,7 +159,7 @@ model_aba$stats
The function computes values predicted in leave-one-out cross-validation, by using the same combination of dependent variables and fitting the regression coefficients with all observations except one. Predicted values can be plotted against field values with the function `lidaRtRee::aba_plot`. It is also informative to check the correlation of prediction errors with other forest or environmental variables.
-The model seems to fail to predict large values, and the prediction errors are positively correlated with basal area.
+The prediction errors seem positively correlated with basal area.
```{r modelCorrelation, include=TRUE, fig.height = 4.5, fig.width = 8}
# check correlation between errors and other variables
@@ -180,7 +180,8 @@ lidaRtRee::aba_plot(model_aba,
legend("topleft", c("public", "private"), col = c("green", "blue"), pch = 1)
plot(model_aba$values$predicted,# plots[subsample, c("G_m2_ha")],
model_aba$values$residual,
- ylab = "Prediction errors", xlab = "Predicted in LOOCV",
+ ylab = "Prediction error", xlab = "Predicted in LOOCV",
+ main = variable,
col = ifelse(plots$stratum == "public", "green", "blue")
)
abline(h = 0, lty = 2)
diff --git a/R/area-based.3.mapping.and.inference.Rmd b/R/area-based.3.mapping.and.inference.Rmd
index a67a867623d1187d1656140f04d7b1c8343629dd..26e5e5dcec3fa233dae2d6cdce90b5ef7243bac1 100755
--- a/R/area-based.3.mapping.and.inference.Rmd
+++ b/R/area-based.3.mapping.and.inference.Rmd
@@ -116,18 +116,18 @@ point_cloud_ground <- lidR::clip_rectangle(cata_altitude, 899600, 6447600, 89990
### Point cloud metrics
-For computation of point cloud metrics, the same function used in `lidaRtRee::clouds_metrics` is now supplied to `lidR::grid_metrics`. Some metrics are displayed hereafter.
+For computation of point cloud metrics, the same function used in `lidaRtRee::clouds_metrics` is now supplied to `lidR::pixel_metrics`. Some metrics are displayed hereafter.
```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
# compute point metrics
-metrics_points <- lidR::grid_metrics(point_cloud_height, aba_point_metrics_fun, res = resolution)
+metrics_points <- lidR::pixel_metrics(point_cloud_height, aba_point_metrics_fun, res = resolution, pkg = "terra")
```
-```{r plotpointMetrics, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 9, fig.height = 3}
+```{r plotpointMetrics, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 3.3}
par(mfrow = c(1, 3))
for (i in c("zmean", "zsd", "imean"))
{
- raster::plot(metrics_points[[i]], main = i)
+ terra::plot(metrics_points[[i]], main = i)
}
```
@@ -143,7 +143,7 @@ It is very important that tree segmentation parameters are the same as in the ca
```{r treeMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 6, fig.height = 6}
# compute chm
-chm <- lidaRtRee::points2DSM(point_cloud_height, res = aba_res_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
# tree top detection (same parameters as used by clouds_tree_metrics in calibration step -> default parameters)
@@ -177,44 +177,44 @@ dummy <- lidaRtRee::raster_metrics(r_treechm,
output = "raster"
)
names(dummy) <- c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
-metrics_trees <- raster::addLayer(metrics_trees, dummy)
+metrics_trees <- c(metrics_trees, dummy)
```
The detected trees are plotted below, with the CHM as background.
-```{r plotTrees, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 6, fig.height = 6}
+```{r plotTrees, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 6, fig.height = 5}
par(mfrow = c(1, 1))
-raster::plot(chm, main = "CHM and detected trees")
-sp::plot(trees, cex = trees$h / 60, add = TRUE)
+terra::plot(chm, main = "CHM and detected trees")
+plot(sf::st_geometry(trees), cex = trees$h / 60, pch = 3, add = TRUE)
```
Some metrics maps are displayed below.
-```{r plottreeMetrics, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 8, fig.height = 3}
+```{r plottreeMetrics, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 3.3}
par(mfrow = c(1, 3))
for (i in c("Tree_meanH", "Tree_sdH", "Tree_density"))
{
- raster::plot(metrics_trees[[i]], main = i)
+ terra::plot(metrics_trees[[i]], main = i)
}
```
### Terrain metrics
-Terrain metrics can be computed from the altitude values of ground points. A wrapping function is required to pass `lidaRtRee::terrain_points_metrics` to `lidR::grid_metrics`.
+Terrain metrics can be computed from the altitude values of ground points. A wrapping function is required to pass `lidaRtRee::terrain_points_metrics` to `lidR::pixel_metrics`.
```{r terrainMetrics, include=TRUE, message=FALSE, warning=FALSE}
# compute terrain metrics
f <- function(x, y, z) {
as.list(lidaRtRee::terrain_points_metrics(data.frame(x, y, z)))
}
-metrics_terrain <- lidR::grid_metrics(point_cloud_ground, ~ f(X, Y, Z), res = resolution)
+metrics_terrain <- lidR::pixel_metrics(point_cloud_ground, ~ f(X, Y, Z), res = resolution)
```
-```{r plotterrainMetrics, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 8, fig.height = 3}
+```{r plotterrainMetrics, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 3.3}
par(mfrow = c(1, 3))
-raster::plot(metrics_terrain[["altitude"]], main = "altitude")
-raster::plot(metrics_terrain$azimut_gr, main = "azimut (gr)", breaks = seq(from = 0, to = 400, length.out = 81), col = rainbow(81))
-raster::plot(metrics_terrain$slope_gr, main = "Slope (gr)")
+terra::plot(metrics_terrain[["altitude"]], main = "altitude")
+terra::plot(metrics_terrain$azimut_gr, main = "azimut (gr)", col = rainbow(128), zlim = c(0,400))
+terra::plot(metrics_terrain$slope_gr, main = "Slope (gr)")
```
### Batch processing of tiles
@@ -249,8 +249,7 @@ 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")])
+ 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,
@@ -273,14 +272,14 @@ metrics <- future.apply::future_lapply(
lidR::filter_poi(a, is.element(Classification, class_points) &
Z <= h_points)
# check number of remaining points
- if (nrow(a@data) == 0)
+ if (nrow(a) == 0)
return(NULL)
# set negative heights to 0
- a@data$Z[a@data$Z < 0] <- 0
+ a$Z[a$Z < 0] <- 0
#
# compute chm
chm <-
- lidaRtRee::points2DSM(a, res = aba_res_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
@@ -303,7 +302,7 @@ metrics <- future.apply::future_lapply(
return(NULL)
# compute raster metrics
metrics_trees <- lidaRtRee::raster_metrics(
- trees[,-1],
+ trees,
res = resolution,
fun = function(x) {
lidaRtRee::std_tree_metrics(x, resolution ^ 2 / 10000)
@@ -317,7 +316,7 @@ metrics <- future.apply::future_lapply(
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])),
+ terra::crop(r_treechm, terra::ext(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,
@@ -332,44 +331,49 @@ metrics <- future.apply::future_lapply(
# remove buffer points
a <- lidR::filter_poi(a, buffer == FALSE)
#
- if (nrow(a@data) == 0)
+ if (nrow(a) == 0)
return(NULL)
# all points metrics
metrics_points <-
- lidR::grid_metrics(a, aba_point_metrics_fun, res = resolution)
+ lidR::pixel_metrics(a, aba_point_metrics_fun, res = resolution)
#
# extend / crop to match metrics_points
metrics_trees <-
- raster::extend(metrics_trees, metrics_points, values = NA)
+ terra::extend(metrics_trees, metrics_points)
metrics_trees2 <-
- raster::extend(metrics_trees2, metrics_points, values = NA)
+ terra::extend(metrics_trees2, metrics_points)
metrics_trees <-
- raster::crop(metrics_trees, metrics_points)
+ terra::crop(metrics_trees, metrics_points)
metrics_trees2 <-
- raster::crop(metrics_trees2, metrics_points)
+ terra::crop(metrics_trees2, metrics_points)
# merge rasterstacks
metrics <- metrics_points
metrics <-
- raster::addLayer(metrics, metrics_trees)
+ c(metrics, metrics_trees)
metrics <-
- raster::addLayer(metrics, metrics_trees2)
- return(metrics)
+ c(metrics, metrics_trees2)
+ if (is.null(metrics))
+ {
+ return(NULL)
+ } else {
+ return(terra::wrap(metrics))
+ }
}
)
# remove NULL rasters
metrics <- metrics[which(!sapply(metrics, is.null))]
+# unwrap rasters
+metrics <- lapply(metrics, terra::rast)
# merge rasters
-names_metrics <- names(metrics[[1]])
-# error might occur if one raster has only one cell
-metrics_map <- do.call(raster::merge, metrics)
-names(metrics_map) <- names_metrics
+# error might occur if one raster has only one cell (check if still relevant in terra package)
+metrics_map <- do.call(terra::merge, metrics)
```
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 grid_metrics
+# 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)))
}
@@ -378,8 +382,8 @@ 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
+ 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]
@@ -388,25 +392,25 @@ metrics_terrain <- future.apply::future_lapply(
if (class(a) == "try-error") {
return(NULL)
}
- lidR::grid_metrics(a, ~ f(X, Y, Z), res = resolution)
+ terra::wrap(lidR::pixel_metrics(a, ~ f(X, Y, Z), res = resolution))
}
)
# remove NULL rasters
-metrics <- metrics[which(!sapply(metrics, is.null))]
+metrics_terrain <- metrics_terrain[which(!sapply(metrics_terrain, is.null))]
+# unwrap rasters
+metrics_terrain <- lapply(metrics_terrain, terra::rast)
# mosaic rasters
-names_metrics <- names(metrics_terrain[[1]])
-metrics_terrain <- do.call(raster::merge, metrics_terrain)
-names(metrics_terrain) <- names_metrics
+metrics_terrain <- do.call(terra::merge, metrics_terrain)
```
Finally all metrics are combined in a single object.
```{r allMetrics, include=TRUE, message=FALSE, warning=FALSE}
-metrics_terrain <- raster::extend(metrics_terrain, metrics_map, values = NA)
-metrics_terrain <- raster::crop(metrics_terrain, metrics_map)
-metrics_terrain <- raster::dropLayer(metrics_terrain, "adjR2.plane")
-# merge rasterstacks
-metrics_map <- raster::addLayer(metrics_map, metrics_terrain)
+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)
```
If the metrics maps are to be displayed in external software, the "tif" format can be used to produce one single with all bands. Meanwhile band names are not retained, so it is useful to save them in a separated `R` archive.
@@ -414,7 +418,7 @@ If the metrics maps are to be displayed in external software, the "tif" format c
```{r saveMetrics, include=TRUE, message=FALSE, warning=FALSE}
metrics_names <- names(metrics_map)
save(metrics_names, file = "../data/aba.model/output/metrics_names.rda")
-raster::writeRaster(metrics_map, file = "../data/aba.model/output/metrics_map.tif", overwrite = TRUE)
+terra::writeRaster(metrics_map, file = "../data/aba.model/output/metrics_map.tif", overwrite = TRUE)
save(metrics_map, file = "../data/aba.model/output/metrics_map.rda")
```
@@ -424,9 +428,9 @@ save(metrics_map, file = "../data/aba.model/output/metrics_map.rda")
Forest parameters are then mapped using the mapped metrics and the previously calibrated model with the function `lidaRtRee::aba_predict`.
-```{r mapForest, include=TRUE, message=FALSE, warning=FALSE, fig.width = 6, fig.height = 4}
+```{r mapForest, include=TRUE, message=FALSE, warning=FALSE, fig.width = 6, fig.height = 3}
prediction_map <- lidaRtRee::aba_predict(model_aba, metrics_map)
-raster::plot(prediction_map, main = "prediction map")
+terra::plot(prediction_map, main = "prediction map")
```
### Stratified models
@@ -435,16 +439,16 @@ In case strata-specific models have been calibrated, it is necessary to add a la
```{r rasterizeStrat, include=TRUE, message=FALSE, warning=FALSE}
# rasterize layer of public forest using value = 1 in polygons
-metrics_map$stratum <- raster::rasterize(public, metrics_map, field = 1)
+metrics_map$stratum <- terra::rasterize(terra::vect(public), metrics_map, field = 1)
# set to 0 in non-public forests areas
metrics_map$stratum[is.na(metrics_map$stratum)] <- 0
# label levels in the raster metadata
-levels(metrics_map$stratum) <- data.frame(ID = c(0, 1), stratum = c("private", "public"))
+levels(metrics_map$stratum) <- data.frame(id = c(0, 1), stratum = c("private", "public"))
# crop public vector
-public_cropped <- sf::st_crop(public, raster::extent(metrics_map))
+public_cropped <- sf::st_crop(public, terra::ext(metrics_map))
```
-Then the function `lidaRtRee::aba_predict` is used to produce the map. If the input is an object containing the stratified model and if the metrics map contains the strata layer, then the model directly maps the output by retaining the corresponding model in each strata.
+Then the function `lidaRtRee::aba_predict` is used to produce the map. If the input is an object containing the stratified model and if the metrics map contains the strata layer, then the model directly maps the output by retaining the corresponding model in each stratum.
```{r mapForestStratified, include=TRUE, message=FALSE, warning=FALSE}
# map with stratified model
@@ -460,18 +464,19 @@ model_public <- list(model = model_aba_stratified_mixed$model$public, stats = mo
prediction_map_public <- lidaRtRee::aba_predict(model_public, metrics_map)
```
-```{r plotmapForestStratified, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 9, fig.height = 6}
+```{r plotmapForestStratified, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 6}
# plot all maps
par(mfrow = c(2, 2))
+# legend range
+limits <- range(c(terra::values(prediction_map_private), terra::values(prediction_map_public)))
# plot raster and vector
-raster::plot(metrics_map$stratum, main = "Ownership", legend = FALSE, col = c("green", "blue"), xaxt = "n", yaxt = "n")
-legend("bottom", legend = c("private", "public"), fill = c("green", "blue"))
-limits <- range(c(raster::values(prediction_map_private), raster::values(prediction_map_public)))
+terra::plot(prediction_map_private, main = "Private model", range = limits, xaxt = "n", yaxt = "n")
+terra::plot(prediction_map_public, main = "Public model", range = limits, xaxt = "n", yaxt = "n")
+terra::plot(metrics_map$stratum, main = "Ownership", col = c("green", "blue"), xaxt = "n", yaxt = "n")
+# legend("bottom", legend = c("private", "public"), fill = c("green", "blue"))
plot(public_cropped, col = NA, add = TRUE)
-raster::plot(prediction_map_mixed, main = "Stratified model", zlim = limits, xaxt = "n", yaxt = "n")
+terra::plot(prediction_map_mixed, main = "Stratified model", range = limits, xaxt = "n", yaxt = "n")
plot(public_cropped, col = NA, add = TRUE)
-raster::plot(prediction_map_private, main = "Private model", zlim = limits, xaxt = "n", yaxt = "n", legend = FALSE)
-raster::plot(prediction_map_public, main = "Public model", zlim = limits, xaxt = "n", yaxt = "n", legend = FALSE)
```
### Forest mask and thresholds
@@ -486,22 +491,25 @@ The first step is to rasterize the forest mask. In the current example, all the
```{r cleanMap, include=TRUE, message=FALSE, warning=FALSE}
# rasterize the forest mask
-forest_mask <- raster::rasterize(forest, metrics_map, field = "FID")
+forest_mask <- terra::rasterize(terra::vect(forest), metrics_map, field = "FID")
# set values to 1 in polygon with ID 27, 0 outside
-raster::values(forest_mask) <- ifelse(raster::values(forest_mask) == 27, 1, NA)
+terra::values(forest_mask) <- ifelse(terra::values(forest_mask) == 27, 1, NA)
+# label levels in the raster metadata
+# forest_mask <- terra::as.factor(forest_mask)
+levels(forest_mask) <- data.frame(id = c(1), forest = c("forest"))
# crop forest vector
-forest_cropped <- sf::st_crop(forest, raster::extent(metrics_map))
+forest_cropped <- sf::st_crop(forest, terra::ext(metrics_map))
# clean raster
# threshold values and set to NA outside of forest
prediction_map_mixed_clean <- lidaRtRee::clean_raster(prediction_map_mixed, c(0, 80), forest_mask)
```
-```{r plotcleanMap, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 2.8}
-par(mfrow = c(1, 3))
-raster::plot(forest_mask, main = "Forest mask")
+```{r plotcleanMap, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 3.02}
+par(mfrow = c(1, 2))
+terra::plot(forest_mask, main = "Forest mask")
plot(sf::st_geometry(forest_cropped), add = TRUE)
-raster::plot(prediction_map_mixed, main = "Stratified predictions map")
-raster::plot(prediction_map_mixed_clean, main = "Masked map")
+# terra::plot(prediction_map_mixed, main = "Stratified predictions map")
+terra::plot(prediction_map_mixed_clean, main = "Masked map")
plot(sf::st_geometry(forest_cropped), add = TRUE)
```
diff --git a/R/coregistration.Rmd b/R/coregistration.Rmd
index 9cc2896b9c05a23469360f7d0809927702e0dcdc..635aef1a9b824bf22f9d9cbad3b4488b2f5e6b72 100755
--- a/R/coregistration.Rmd
+++ b/R/coregistration.Rmd
@@ -87,6 +87,7 @@ The code to extract LAS objects from a directory containing the LAS files is (co
```{r extractlas, include = TRUE, eval=FALSE}
# create catalog of LAS files
cata <- lidR::readALSLAScatalog("/directory_with_classified_laz_files/")
+# "/media/reseau/lessem/ProjetsCommuns/Lidar/data/Franche-Comte/norm.laz/"
# set coordinate system
lidR::projection(cata) <- 2154
# option to read only xyzc attributes (coordinates, intensity, echo order and classification) from files
@@ -124,11 +125,11 @@ s_step <- 0.5
The first step is to compute the canopy height model from the ALS data, and remove artefacts by thresholding extreme values and applying a median filter.
-```{r chm, include = TRUE, out.width='70%', fig.width = 9, fig.height = 4.5}
+```{r chm, include = TRUE, out.width='70%', fig.width = 9, fig.height = 3.5}
# Choose a plot as example
i <- 10
# compute canopy height model
-chm <- lidaRtRee::points2DSM(las[[i]], res = r_res)
+chm <- lidR::rasterize_canopy(las[[i]], res = r_res, algorithm = lidR::p2r())
# apply threshold to canopy height model (CHM)
chm[chm > hmax] <- hmax
# fill NA values with 0
@@ -137,15 +138,15 @@ chm[is.na(chm)] <- 0
chmfilt <- lidaRtRee::dem_filtering(chm, "Median", 3, sigmap = 0)$non_linear_image
# plot canopy height model
par(mfrow = c(1, 2))
-raster::plot(chm, main = "Raw canopy height model")
-raster::plot(chmfilt, main = "Filtered canopy height model")
+terra::plot(chm, main = "Raw canopy height model")
+terra::plot(chmfilt, main = "Filtered canopy height model")
```
### Plot mask from tree inventory
The trees corresponding to the plot are extracted, and a plot mask is computed from the plot center and radius.
-```{r mask, include = TRUE, out.width = '50%', fig.dim=c(5.5, 4.5), warning=FALSE}
+```{r mask, include = TRUE, out.width = '40%', fig.dim=c(4.5, 3), warning=FALSE}
# plot centre
centre <- p[i, c("XGPS", "YGPS")]
# extract plot trees
@@ -165,9 +166,9 @@ binary = T
# replace 0 by NA
r_mask[r_mask == 0] <- NA
# specify projection
-raster::crs(r_mask) <- 2154
+terra::crs(r_mask) <- "epsg:2154"
# display plot mask
-raster::plot(r_mask, main = "Plot mask and tree positions", legend = FALSE)
+terra::plot(r_mask, main = "Plot mask and tree positions", legend = FALSE)
# add tree positions
points(trees[, c("x", "y")], cex = trees[, "dia"] / 40)
# add plot center
@@ -190,13 +191,13 @@ round(coreg$local_max, 2)
```
The maximum of the correlation image is located at (`r coreg$local_max$dx1`, `r coreg$local_max$dy1`), given by attributes `dx1` and `dy1`.
-```{r plotcorrelation, include = TRUE, out.width='100%', fig.width = 12, fig.height = 4.5}
+```{r plotcorrelation, include = TRUE, out.width='100%', fig.width = 12, fig.height = 4}
par(mfrow = c(1, 3))
-raster::plot(chmfilt, main = "Initial tree positions and CHM")
+terra::plot(chmfilt, main = "Initial tree positions and CHM")
# display initial tree positions
graphics::points(trees$x, trees$y, cex = trees$dia / 40)
# display correlation image
-raster::plot(coreg$correlation_raster,
+terra::plot(coreg$correlation_raster,
main = "Correlation image",
col = cm.colors(16)
)
@@ -206,7 +207,7 @@ abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
legend("topleft", "Maximum", pch = 4)
#
-raster::plot(chmfilt, main = "Coregistered tree positions and CHM")
+terra::plot(chmfilt, main = "Coregistered tree positions and CHM")
# display coregistered tree positions
graphics::points(trees$x + coreg$local_max$dx1, trees$y + coreg$local_max$dy1,
cex = trees$dia / 40, col = "red"
@@ -244,7 +245,7 @@ l_chm <- future.apply::future_lapply(
as.list(1:length(las)),
FUN = function(i) {
# compute CHM
- chm <- lidaRtRee::points2DSM(las[[i]], res = r_res)
+ chm <- lidR::rasterize_canopy(las[[i]], res = r_res, algorithm = lidR::p2r(), pkg = "terra")
# apply threshold to canopy height model (CHM)
chm[chm > hmax] <- hmax
# fill NA values with 0
@@ -252,7 +253,7 @@ l_chm <- future.apply::future_lapply(
# apply median filter (3x3 window) to CHM
chmfilt <-
lidaRtRee::dem_filtering(chm, "Median", 3, sigmap = 0)[[1]]
- return(chmfilt)
+ return(terra::wrap(chmfilt))
}
)
```
@@ -280,8 +281,8 @@ l_field <- future.apply::future_lapply(
# replace 0 by NA
r_mask[r_mask == 0] <- NA
# specify projection
- raster::crs(r_mask) <- 2154
- return(list(trees = trees[, c("x", "y", "dia")], r_mask = r_mask))
+ terra::crs(r_mask) <- "epsg:2154"
+ return(list(trees = trees[, c("x", "y", "dia")], r_mask = terra::wrap(r_mask)))
}
)
```
@@ -293,15 +294,19 @@ Then the correlation image and corresponding statistics are computed for each pl
l_coreg <- future.apply::future_lapply(
as.list(1:length(las)),
FUN = function(i) {
- lidaRtRee::coregistration(
- l_chm[[i]],
+ dummy <- lidaRtRee::coregistration(
+ terra::rast(l_chm[[i]]),
l_field[[i]]$trees,
- mask = l_field[[i]]$r_mask,
+ mask = terra::rast(l_field[[i]]$r_mask),
buffer = b_size,
step = s_step,
dm = 2,
plot = FALSE
)
+ # wrap raster
+ dummy$correlation_raster <- terra::wrap(dummy$correlation_raster)
+ # output
+ return(dummy)
}
)
```
@@ -330,11 +335,16 @@ round(head(translations, n = 3L), 2)
### Export graphics as pdf files
```{r batch.plots, include = TRUE, eval=FALSE}
+# unwrap raster files
+l_chm <- lapply(l_chm, terra::rast)
+# l_coreg <- lapply(l_coreg, function(x) list(correlation_raster = terra::rast(x$correlation_raster), local_max = x$local_max))
+# l_field <- lapply(l_field, function(x) list(trees = x$trees, r_mask = terra::rast(x$r_mask)))
+#
for (i in 1:length(las))
{
# CHM
pdf(file = paste("Coregistration_", p$placette[i], ".pdf", sep = ""))
- raster::plot(l_chm[[i]], asp = 1)
+ terra::plot(l_chm[[i]], asp = 1)
# display initial tree positions
graphics::points(l_field[[i]]$trees$x, l_field[[i]]$trees$y,
cex = l_field[[i]]$trees$dia / 40
diff --git a/R/forest.structure.metrics.Rmd b/R/forest.structure.metrics.Rmd
index 8c1640906e1d9fb2a97013ae339935f5e6cba611..8cd0b472362c5ef80edcfa887e50d6939ed6c4fe 100755
--- a/R/forest.structure.metrics.Rmd
+++ b/R/forest.structure.metrics.Rmd
@@ -67,6 +67,25 @@ class_ground <- 2
sigma <- c(0, 1, 2, 4, 8, 16)
# convert vector of sigma values in list
sigma_l <- as.list(sigma)
+# fonction ti computed raster statistics from multi-scale smoothing
+smoothed_raster_stats <- function(x) {
+ data.frame(
+ CHM0_sd = sd(x$smoothed_image_0),
+ CHM1_sd = sd(x$smoothed_image_1),
+ CHM2_sd = sd(x$smoothed_image_2),
+ CHM4_sd = sd(x$smoothed_image_4),
+ CHM8_sd = sd(x$smoothed_image_8),
+ CHM16_sd = sd(x$smoothed_image_16),
+ CHM_mean = mean(x$smoothed_image_0),
+ CHM_PercInf0_5 = sum(x$smoothed_image_0 < 0.5) * mf,
+ CHM_PercInf1 = sum(x$smoothed_image_0 < 1) * mf,
+ CHM_PercSup5 = sum(x$smoothed_image_0 > 5) * mf,
+ CHM_PercSup10 = sum(x$smoothed_image_0 > 10) * mf,
+ CHM_PercSup20 = sum(x$smoothed_image_0 > 20) * mf,
+ CHM_Perc1_5 = (sum(x$smoothed_image_0 < 5) - sum(x$smoothed_image_0 < 1)) * mf,
+ CHM_Perc2_5 = (sum(x$smoothed_image_0 < 5) - sum(x$smoothed_image_0 < 2)) * mf
+ )
+}
# height breaks for penetration ratio and density
breaks_h <- c(-Inf, 0, 0.5, 1, 2, 5, 10, 20, 30, 60, Inf)
# percentiles of height distribution
@@ -95,7 +114,7 @@ lidR::opt_progress(cata) <- FALSE
# option to read only xyzc attributes (coordinates, intensity, echo order and classification) from height files (saves memory)
lidR::opt_select(cata) <- "xyzirnc"
# display
-sp::plot(cata, main = "Bounding box of las files")
+lidR::plot(cata, main = "Bounding box of las files")
```
## Workflow exemplified on one tile
@@ -107,7 +126,7 @@ The tile is loaded, plus a buffer on adjacent tiles. Buffer points have a column
# tile to process
i <- 1
# tile extent
-b_box <- as.numeric(cata@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
+b_box <- sf::st_bbox(cata[i,])
# load tile extent plus buffer
a <- lidR::clip_rectangle(
cata, b_box[1] - buffer_size, b_box[2] - buffer_size, b_box[3] + buffer_size,
@@ -126,7 +145,7 @@ Only points from desired classes are retained. In case some mis-classified, high
# remove unwanted point classes, and points higher than height threshold
a <- lidR::filter_poi(a, is.element(Classification, class_points) & Z <= points_max_h)
# set negative heights to 0
-a@data$Z[a@data$Z < 0] <- 0
+a$Z[a$Z < 0] <- 0
#
summary(a@data)
```
@@ -143,63 +162,46 @@ The CHM is computed and NA values are replaced by 0. A check is performed to mak
```{r computeCHM, include = TRUE, out.width = '50%', fig.dim=c(5, 4.5), message=FALSE}
# compute chm
-chm <- lidaRtRee::points2DSM(a, res = res_chm)
+chm <- lidR::rasterize_canopy(a, res = res_chm, algorithm = lidR::p2r(), pkg = "terra")
# replace NA, low and high values
chm[is.na(chm) | chm < 0 | chm > points_max_h] <- 0
# display CHM
-raster::plot(chm, asp = 1, main = "Canopy height model")
+terra::plot(chm, asp = 1, main = "Canopy height model")
```
## Metrics computation
### 2D CHM metrics
-The CHM is smoothed with a Gaussian filter with different `sigma` values. Smoothed results are stored in a list and then converted to a raster stack.
+The CHM is smoothed with a Gaussian filter with different `sigma` values. Smoothed results are stored in a list and then integrated ito a single raster
-```{r 2dchmMetrics, include = TRUE, fig.dim = c(9, 4.5), out.width='100%', message=FALSE}
+```{r 2dchmMetrics, include = TRUE, fig.dim = c(9, 4), out.width='100%', message=FALSE}
# for each value in list of sigma, apply filtering to chm and store result in list
st <- lapply(sigma_l, FUN = function(x) {
lidaRtRee::dem_filtering(chm, nl_filter = "Closing", nl_size = 5, sigmap = x)$smoothed_image
})
# convert to raster
-st <- raster::stack(st)
+st <- terra::rast(st)
+# modify layer names
+names(st) <- paste0(names(st), "_", unlist(sigma_l))
# display
-raster::plot(st)
+terra::plot(st)
```
-The raster stack is then converted to points. Points outside the tile extent are removed and summary statistics (mean and standard deviation) of smoothed CHM heights are computed for each pixel at the output metrics resolution. Canopy cover in different height layers are also computed for the initial CHM after applying a non-linear filter designed to fill holes.
+The raster is then converted to points. Points outside the tile extent are removed and summary statistics (mean and standard deviation) of smoothed CHM heights are computed for each pixel at the output metrics resolution. Canopy cover in different height layers are also computed for the initial CHM after applying a non-linear filter designed to fill holes.
```{r 2dmetrics2points, include = TRUE, message=FALSE, warning=FALSE}
# crop to bbox
-st <- raster::crop(st, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4]))
+st <- terra::crop(st, terra::ext(b_box[1], b_box[3], b_box[2], b_box[4]))
# multiplying factor to compute proportion
mf <- 100 / (resolution / res_chm)^2
-# modify layer names
-names(st) <- paste0("layer.", 1:6)
+
# compute raster metrics
metrics_2dchm <- lidaRtRee::raster_metrics(st,
res = resolution,
- fun = function(x) {
- c(
- sd(x$layer.1), sd(x$layer.2), sd(x$layer.3),
- sd(x$layer.4), sd(x$layer.5), sd(x$layer.6),
- mean(x$layer.1), sum(x$layer.1 < 0.5) * mf,
- sum(x$layer.1 < 1) * mf,
- sum(x$layer.1 > 5) * mf,
- sum(x$layer.1 > 10) * mf,
- sum(x$layer.1 > 20) * mf,
- (sum(x$layer.1 < 5) - sum(x$layer.1 < 1)) * mf,
- (sum(x$layer.1 < 5) - sum(x$layer.1 < 2)) * mf
- )
- },
+ fun = smoothed_raster_stats,
output = "raster"
)
-names(metrics_2dchm) <- c(
- "CHM0_sd", "CHM1_sd", "CHM2_sd", "CHM4_sd", "CHM8_sd",
- "CHM16_sd", "CHM_mean", "CHM_PercInf0_5", "CHM_PercInf1",
- "CHM_PercSup5", "CHM_PercSup10", "CHM_PercSup20",
- "CHM_Perc1_5", "CHM_Perc2_5"
-)
#
metrics_2dchm
```
@@ -207,7 +209,7 @@ metrics_2dchm
Some outputs are displayed hereafter. On the first line are the standard deviation of CHM smoothed with different `sigma` values. On the second line are the CHM mean and percentages of CHM values below 1 m and above 10 m.
```{r 2dmetricsOutputs, include = TRUE, fig.dim = c(9, 5), out.width='100%', warning=FALSE}
-raster::plot(metrics_2dchm[[c(
+terra::plot(metrics_2dchm[[c(
"CHM0_sd", "CHM2_sd", "CHM8_sd", "CHM_mean",
"CHM_PercInf1", "CHM_PercSup10"
)]])
@@ -216,7 +218,7 @@ raster::plot(metrics_2dchm[[c(
### Gaps and edges metrics
Gaps are computed with the function `gap_detection`.
-```{r gapMetrics, include = TRUE, fig.dim = c(9, 3), out.width='100%', warning=FALSE, message=FALSE}
+```{r gapMetrics, include = TRUE, fig.dim = c(9, 3), warning=FALSE, message=FALSE}
# compute gaps
gaps <- lidaRtRee::gap_detection(chm,
ratio = 2, gap_max_height = 1,
@@ -225,34 +227,34 @@ gaps <- lidaRtRee::gap_detection(chm,
)
# display maps
par(mfrow = c(1, 3))
-raster::plot(chm, main = "CHM")
-raster::plot(log(gaps$gap_surface), main = "Log of gap surface")
+terra::plot(chm, main = "CHM")
+terra::plot(log(gaps$gap_surface), main = "Log of gap surface")
# display gaps
dummy <- gaps$gap_id
dummy[dummy == 0] <- NA
-raster::plot(dummy %% 8, asp = 1, main = "Gap id (random colors)", col = rainbow(8), legend = FALSE)
+terra::plot(dummy %% 8, asp = 1, main = "Gap id (random colors)", col = rainbow(8), legend = FALSE)
```
Summary statistics are then computed: pixel surface occupied by gaps in different size classes. Results are first cropped to the tile extent.
```{r gapStats, include = TRUE, warning=FALSE, message=FALSE}
# crop results to tile size
-gaps_surface <- raster::crop(gaps$gap_surface, raster::extent(
+gaps_surface <- terra::crop(gaps$gap_surface, terra::ext(
b_box[1], b_box[3],
b_box[2], b_box[4]
))
# compute gap statistics at final metrics resolution, in case gaps exist
-if (!all(is.na(raster::values(gaps_surface)))) {
+if (!all(is.na(terra::values(gaps_surface)))) {
metrics_gaps <- lidaRtRee::raster_metrics(gaps_surface,
res = resolution,
fun = function(x) {
- hist(x$layer,
+ hist(x$lyr.1,
breaks = breaks_gap_surface,
plot = F
)$counts * (res_chm / resolution)^2
}, output = "raster"
)
# compute total gap proportion
- metrics_gaps$sum <- raster::stackApply(metrics_gaps, rep(1, length(names(metrics_gaps))),
+ metrics_gaps$sum <- terra::tapp(metrics_gaps, rep(1, length(names(metrics_gaps))),
fun = sum
)
# if gaps are present
@@ -274,24 +276,24 @@ Edges are detected with the function `edge_detection` as the outer envelope of
edges <- lidaRtRee::edge_detection(gaps$gap_id > 0)
# display maps
par(mfrow = c(1, 3))
-raster::plot(chm, main = "CHM")
-raster::plot(gaps$gap_id > 0, main = "Gaps", legend = FALSE)
-raster::plot(edges, main = "Edges", legend = FALSE)
+terra::plot(chm, main = "CHM")
+terra::plot(gaps$gap_id > 0, main = "Gaps", legend = FALSE)
+terra::plot(edges, main = "Edges", legend = FALSE)
```
The percentage of surface occupied by edges is then computed as summary statistic. Results are first cropped to the tile extent.
```{r edgeStats, include = TRUE, warning=FALSE, message=FALSE}
# crop results to tile size
-edges <- raster::crop(
+edges <- terra::crop(
edges,
- raster::extent(b_box[1], b_box[3], b_box[2], b_box[4])
+ terra::ext(b_box[1], b_box[3], b_box[2], b_box[4])
)
# compute gap statistics at final metrics resolution, in case gaps exist
-if (!all(is.na(raster::values(edges)))) {
+if (!all(is.na(terra::values(edges)))) {
metrics_edges <- lidaRtRee::raster_metrics(edges,
res = resolution,
fun = function(x) {
- sum(x$layer) * (res_chm / resolution)^2
+ sum(x$lyr.1) * (res_chm / resolution)^2
}, output = "raster"
)
names(metrics_edges) <- "edges.proportion"
@@ -309,28 +311,28 @@ Some outputs are displayed hereafter. The proportion of surface occupied by gaps
```{r gapDisplay, include = TRUE, fig.dim = c(9, 3.2), out.width='100%', warning=FALSE, message=FALSE}
par(mfrow = c(1, 3))
-# raster::plot(metrics_gaps[["G_s4to16"]], main="Prop of gaps 4-16m2")
-raster::plot(metrics_gaps[["G_s64to256"]], main = "Prop of gaps 64-256m2")
-raster::plot(metrics_gaps[["G_s4toInf"]], main = "Prop of gaps >4m2")
-raster::plot(metrics_edges, main = "Proportion of edges")
+# terra::plot(metrics_gaps[["G_s4to16"]], main="Prop of gaps 4-16m2")
+terra::plot(metrics_gaps[["G_s64to256"]], main = "Prop of gaps 64-256m2")
+terra::plot(metrics_gaps[["G_s4toInf"]], main = "Prop of gaps >4m2")
+terra::plot(metrics_edges, main = "Proportion of edges")
```
### Tree metrics
Tree tops are detected with the function `treeSegmentation` and then extracted with `treeExtraction`.
-```{r treeMetrics, include = TRUE, fig.dim = c(9, 4.5), out.width='100%', , warning=FALSE, message=FALSE}
+```{r treeMetrics, include = TRUE, fig.width = 8, fig.height = 3.3, warning=FALSE, message=FALSE}
# tree top detection (default parameters)
segms <- lidaRtRee::tree_segmentation(chm, hmin = 5)
# extraction to data.frame
trees <- lidaRtRee::tree_extraction(segms$filled_dem, segms$local_maxima, segms$segments_id)
# display maps
par(mfrow = c(1, 2))
-raster::plot(chm, main = "CHM and tree tops")
+terra::plot(chm, main = "CHM and tree tops")
points(trees$x, trees$y, pch = 4, cex = 0.3)
# display segments, except ground segment
dummy <- segms$segments_id
dummy[dummy == 0] <- NA
-raster::plot(dummy %% 8, asp = 1, main = "Segments (random colors)", col = rainbow(8), legend = FALSE)
+terra::plot(dummy %% 8, asp = 1, main = "Segments (random colors)", col = rainbow(8))
points(trees$x, trees$y, pch = 4, cex = 0.3)
```
@@ -355,9 +357,9 @@ r_tree_chm <- segms$filled_dem
# set chm to NA in non segment area
r_tree_chm[segms$segments_id == 0] <- NA
# compute raster metrics
-dummy <- lidaRtRee::raster_metrics(raster::crop(
+dummy <- lidaRtRee::raster_metrics(terra::crop(
r_tree_chm,
- raster::extent(
+ terra::ext(
b_box[1], b_box[3], b_box[2],
b_box[4]
)
@@ -373,21 +375,21 @@ output = "raster"
)
names(dummy) <- c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
#
-dummy <- raster::extend(dummy, metrics_trees)
-metrics_trees <- raster::addLayer(metrics_trees, dummy)
+dummy <- terra::extend(dummy, metrics_trees)
+metrics_trees <- c(metrics_trees, dummy)
#
metrics_trees
```
Some outputs are displayed hereafter.
-```{r displayTreeMetrics, include = TRUE, fig.width = 9, fig.height = 3.2, warning=FALSE, message=FALSE}
+```{r displayTreeMetrics, include = TRUE, fig.width = 12, fig.height = 3.3, warning=FALSE, message=FALSE}
par(mfrow = c(1, 3))
-raster::plot(metrics_trees$Tree_meanH, main = "Mean (detected) tree height (m)")
+terra::plot(metrics_trees$Tree_meanH, main = "Mean (detected) tree height (m)")
points(trees$x, trees$y, cex = trees$h / 40)
-raster::plot(metrics_trees$TreeSup20_density,
+terra::plot(metrics_trees$TreeSup20_density,
main = "Density of (detected) trees > 20m (/ha)"
)
-raster::plot(metrics_trees$Tree_meanCrownVolume,
+terra::plot(metrics_trees$Tree_meanCrownVolume,
main = "Mean crown volume of detected trees (m3)"
)
```
@@ -401,7 +403,7 @@ Buffer points are first removed from the point cloud, then metrics are computed
a <- lidR::filter_poi(a, buffer == 0)
#
# all points metrics
-metrics_1d <- lidR::grid_metrics(a, as.list(c(
+metrics_1d <- lidR::pixel_metrics(a, as.list(c(
as.vector(quantile(Z, probs = percent)),
sum(Classification == 2),
mean(Intensity[ReturnNumber == 1])
@@ -410,7 +412,7 @@ names(metrics_1d) <- c(paste("Hp", percent * 100, sep = ""), "nb_ground", "Imean
#
# vegetation-only metrics
a <- lidR::filter_poi(a, Classification != class_ground)
-dummy <- lidR::grid_metrics(a, as.list(c(
+dummy <- lidR::pixel_metrics(a, as.list(c(
mean(Z),
max(Z),
sd(Z),
@@ -419,18 +421,18 @@ dummy <- lidR::grid_metrics(a, as.list(c(
)), res = resolution)
names(dummy) <- c("Hmean", "Hmax", "Hsd", "nb_veg", n_breaks_h)
# merge rasterstacks
-dummy <- raster::extend(dummy, metrics_1d)
-metrics_1d <- raster::addLayer(metrics_1d, dummy)
+dummy <- terra::extend(dummy, metrics_1d)
+metrics_1d <- c(metrics_1d, dummy)
rm(dummy)
# crop to tile extent
-metrics_1d <- raster::crop(metrics_1d, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4]))
+metrics_1d <- terra::crop(metrics_1d, terra::ext(b_box[1], b_box[3], b_box[2], b_box[4]))
```
Based on those metrics, additional metrics are computed (Simpson index, relative density in height bins and penetration ratio in height bins)
```{r 1dmetricsAdditional, include = TRUE, warnings=FALSE}
# Simpson index
-metrics_1d$Hsimpson <- raster::stackApply(
+metrics_1d$Hsimpson <- terra::tapp(
metrics_1d[[n_breaks_h[c(-1, -length(n_breaks_h))]]],
1,
function(x, ...) {
@@ -448,7 +450,7 @@ for (i in n_breaks_h[c(-1, -length(n_breaks_h))])
cumu_sum <- metrics_1d[["nb_ground"]]
for (i in n_breaks_h)
{
- cumu_sum <- raster::addLayer(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics_1d[[i]])
+ cumu_sum <- c(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics_1d[[i]])
}
names(cumu_sum) <- c("nb_ground", n_breaks_h)
# compute interception ratio of each layer
@@ -459,7 +461,7 @@ for (i in 2:dim(cumu_sum)[3])
}
names(intercep_ratio) <- gsub("nb", "intercepRatio", names(cumu_sum)[-1])
# merge rasterstacks
-metrics_1d <- raster::addLayer(metrics_1d, intercep_ratio)
+metrics_1d <- c(metrics_1d, intercep_ratio)
# rm(cumu_sum, intercep_ratio)
#
metrics_1d
@@ -467,17 +469,17 @@ metrics_1d
Some outputs are displayed hereafter.
-```{r display1dmetrics, include = TRUE, fig.width = 9, fig.height = 3.2, warning=FALSE, message=FALSE}
+```{r display1dmetrics, include = TRUE, fig.width = 12, fig.height = 3.3, warning=FALSE, message=FALSE}
# display results
par(mfrow = c(1, 3))
-raster::plot(metrics_1d$Imean_1stpulse, main = "Mean intensity (1st return)")
-raster::plot(metrics_1d$Hp50, main = "Percentile 50 of heights")
-raster::plot(metrics_1d$intercepRatio_H2to5, main = "Interception ratio 2-5 m")
+terra::plot(metrics_1d$Imean_1stpulse, main = "Mean intensity (1st return)")
+terra::plot(metrics_1d$Hp50, main = "Percentile 50 of heights")
+terra::plot(metrics_1d$intercepRatio_H2to5, main = "Interception ratio 2-5 m")
```
## Batch processing
-The following code uses parallel processing to handle multiple files of a catalog, and arranges all metrics in a raster stack. Code in the "parameters" section has to be run beforehand.
+The following code uses parallel processing to handle multiple files of a catalog, and arranges all metrics in a raster. Code in the "parameters" section has to be run beforehand.
```{r batchProcessing, include = TRUE, eval=TRUE, warning=FALSE, message=FALSE, fig.width = 9, fig.height = 2.6}
# processing by tile
@@ -485,8 +487,7 @@ metrics <- future.apply::future_lapply(
as.list(1:nrow(cata)),
FUN = function(i) {
# tile extent
- b_box <-
- as.numeric(cata@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
+ b_box <- sf::st_bbox(cata[i,])
# load tile extent plus buffer
a <- try(lidR::clip_rectangle(
cata,
@@ -510,13 +511,13 @@ metrics <- future.apply::future_lapply(
lidR::filter_poi(a,
is.element(Classification, class_points) & Z <= points_max_h)
# check number of remaining points
- if (nrow(a@data) == 0)
+ if (nrow(a) == 0)
return(NULL)
# set negative heights to 0
- a@data$Z[a@data$Z < 0] <- 0
+ a$Z[a$Z < 0] <- 0
#
# compute chm
- chm <- lidaRtRee::points2DSM(a, res = res_chm)
+ chm <- lidR::rasterize_canopy(a, res = res_chm, algorithm = lidR::p2r(), pkg = "terra")
# replace NA, low and high values
chm[is.na(chm) | chm < 0 | chm > points_max_h] <- 0
#-----------------------
@@ -528,14 +529,16 @@ metrics <- future.apply::future_lapply(
lidaRtRee::dem_filtering(chm,
nl_filter = "Closing",
nl_size = 5,
- sigmap = x)[[2]]
+ sigmap = x)$smoothed_image
}
)
# convert to raster
- st <- raster::stack(st)
+ st <- terra::rast(st)
+ # modify layer names
+ names(st) <- paste0(names(st), "_", unlist(sigma_l))
# crop to bbox
st <-
- raster::crop(st, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4]))
+ terra::crop(st, terra::ext(b_box[1], b_box[3], b_box[2], b_box[4]))
# multiplying factor to compute proportion
mf <- 100 / (resolution / res_chm) ^ 2
# compute raster metrics
@@ -543,24 +546,7 @@ metrics <- future.apply::future_lapply(
lidaRtRee::raster_metrics(
st,
res = resolution,
- fun = function(x) {
- data.frame(
- CHM0.sd = sd(x$non_linear_image),
- CHM1.sd = sd(x$smoothed_image.1),
- CHM2.sd = sd(x$smoothed_image.2),
- CHM4.sd = sd(x$smoothed_image.4),
- CHM8.sd = sd(x$smoothed_image.8),
- CHM16.sd = sd(x$smoothed_image.16),
- CHM.mean = mean(x$non_linear_image),
- CHM.PercInf0.5 = sum(x$non_linear_image < 0.5) * mf,
- CHM.PercInf1 = sum(x$non_linear_image < 1) * mf,
- CHM.PercSup5 = sum(x$non_linear_image > 5) * mf,
- CHM.PercSup10 = sum(x$non_linear_image > 10) * mf,
- CHM.PercSup20 = sum(x$non_linear_image > 20) * mf,
- CHM.Perc1_5 = (sum(x$non_linear_image < 5) - sum(x$non_linear_image < 1)) * mf,
- CHM.Perc2_5 = (sum(x$non_linear_image < 5) - sum(x$non_linear_image < 2)) * mf
- )
- },
+ fun = smoothed_raster_stats,
output = "raster"
)
#
@@ -577,16 +563,16 @@ metrics <- future.apply::future_lapply(
gap_reconstruct = TRUE
)
# crop results to tile size
- gaps_surface <- raster::crop(gaps$gap_surface,
- raster::extent(b_box[1], b_box[3], b_box[2],
+ gaps_surface <- terra::crop(gaps$gap_surface,
+ terra::ext(b_box[1], b_box[3], b_box[2],
b_box[4]))
# compute gap statistics at final metrics resolution, in case gaps exist
- if (!all(is.na(raster::values(gaps_surface)))) {
+ if (!all(is.na(terra::values(gaps_surface)))) {
metrics_gaps <- lidaRtRee::raster_metrics(
gaps_surface,
res = resolution,
fun = function(x) {
- hist(x$layer,
+ hist(x$lyr.1,
breaks = breaks_gap_surface,
plot = F)$counts * (res_chm / resolution) ^
2
@@ -594,9 +580,10 @@ metrics <- future.apply::future_lapply(
output = "raster"
)
# compute total gap proportion
- metrics_gaps$sum <-
- raster::stackApply(metrics_gaps, rep(1, length(names(metrics_gaps))),
- fun = sum)
+ metrics_gaps$sum <- terra::tapp(metrics_gaps,
+ rep(1, length(names(metrics_gaps))),
+ fun = sum
+ )
# if gaps are present
names(metrics_gaps) <-
c(n_breaks_gap, paste("G_s", min(breaks_gap_surface), "toInf", sep = ""))
@@ -612,15 +599,15 @@ metrics <- future.apply::future_lapply(
# Perform edge detection by erosion
edges <- lidaRtRee::edge_detection(gaps$gap_id > 0)
# crop results to tile size
- edges <- raster::crop(edges, raster::extent(b_box[1], b_box[3],
+ edges <- terra::crop(edges, terra::ext(b_box[1], b_box[3],
b_box[2], b_box[4]))
# compute edges statistics at final metrics resolution, in case edges exist
- if (!all(is.na(raster::values(edges)))) {
+ if (!all(is.na(terra::values(edges)))) {
metrics_edges <- lidaRtRee::raster_metrics(
edges,
res = resolution,
fun = function(x) {
- sum(x$layer) * (res_chm / resolution) ^ 2
+ sum(x$lyr.1) * (res_chm / resolution) ^ 2
},
output = "raster"
)
@@ -657,7 +644,9 @@ metrics <- future.apply::future_lapply(
#
# extend layer in case there are areas without trees
metrics_trees <-
- raster::extend(metrics_trees, metrics_2dchm, values = 0)
+ terra::extend(metrics_trees, metrics_2dchm)
+ # set NA values to 0
+ metrics_trees[is.na(metrics_trees)] <- 0
# compute canopy cover in trees and canopy mean height in trees
# in region of interest, because it is not in previous step.
r_tree_chm <- segms$filled_dem
@@ -665,7 +654,7 @@ metrics <- future.apply::future_lapply(
r_tree_chm[segms$segments_id == 0] <- NA
# compute raster metrics
dummy <- lidaRtRee::raster_metrics(
- raster::crop(r_tree_chm, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4])),
+ terra::crop(r_tree_chm, terra::ext(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 / res_chm) ^ 2,
@@ -675,20 +664,22 @@ metrics <- future.apply::future_lapply(
)
names(dummy) <-
c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
- dummy <- raster::extend(dummy, metrics_2dchm, values = 0)
+ dummy <- terra::extend(dummy, metrics_2dchm)
+ # set NA values to 0
+ dummy[is.na(dummy)] <- 0
#
- metrics_trees <- raster::addLayer(metrics_trees, dummy)
+ metrics_trees <-c(metrics_trees, dummy)
#
# -------------------------
# compute 1D height metrics
# remove buffer points
a <- lidR::filter_poi(a, buffer == 0)
#
- if (nrow(a@data) == 0) {
+ if (nrow(a) == 0) {
metrics_1d <- NULL
} else {
# all points metrics
- metrics_1d <- lidR::grid_metrics(a, as.list(c(
+ metrics_1d <- lidR::pixel_metrics(a, as.list(c(
as.vector(quantile(Z, probs = percent)),
sum(Classification == 2),
mean(Intensity[ReturnNumber == 1])
@@ -700,8 +691,8 @@ metrics <- future.apply::future_lapply(
#
# vegetation-only metrics
a <- lidR::filter_poi(a, Classification != class_ground)
- if (nrow(a@data) != 0) {
- dummy <- lidR::grid_metrics(a, as.list(c(
+ if (nrow(a) != 0) {
+ dummy <- lidR::pixel_metrics(a, as.list(c(
mean(Z),
max(Z),
sd(Z),
@@ -715,24 +706,27 @@ metrics <- future.apply::future_lapply(
)), res = resolution)
names(dummy) <-
c("Hmean", "Hmax", "Hsd", "nb_veg", n_breaks_h)
- # merge rasterstacks
- dummy <- raster::extend(dummy, metrics_1d)
- metrics_1d <- raster::addLayer(metrics_1d, dummy)
+ # merge raster
+ dummy <- terra::extend(dummy, metrics_1d)
+ metrics_1d <- c(metrics_1d, dummy)
rm(dummy)
# -------------
# merge metrics
metrics_1d <-
- raster::extend(metrics_1d, metrics_2dchm, values = 0)
- metrics_1d <- raster::crop(metrics_1d, metrics_2dchm)
+ terra::extend(metrics_1d, metrics_2dchm)
+ # metrics_1d[is.na(metrics_1d)] <- 0
+ metrics_1d <- terra::crop(metrics_1d, metrics_2dchm)
metrics_gaps <-
- raster::extend(metrics_gaps, metrics_2dchm, values = 0)
- metrics_gaps <- raster::crop(metrics_gaps, metrics_2dchm)
+ terra::extend(metrics_gaps, metrics_2dchm)
+ # metrics_gaps[is.na(metrics_gaps)] <- 0
+ metrics_gaps <- terra::crop(metrics_gaps, metrics_2dchm)
metrics_edges <-
- raster::extend(metrics_edges, metrics_2dchm, values = 0)
- metrics_edges <- raster::crop(metrics_edges, metrics_2dchm)
- metrics_trees <- raster::crop(metrics_trees, metrics_2dchm)
+ terra::extend(metrics_edges, metrics_2dchm)
+ # metrics_edges[is.na(metrics_edges)] <- 0
+ metrics_edges <- terra::crop(metrics_edges, metrics_2dchm)
+ metrics_trees <- terra::crop(metrics_trees, metrics_2dchm)
#
- temp <- raster::addLayer(metrics_1d,
+ temp <- c(metrics_1d,
metrics_2dchm,
metrics_gaps,
metrics_edges,
@@ -740,18 +734,20 @@ metrics <- future.apply::future_lapply(
temp[is.na(temp)] <- 0
} # end of vegetation-only points presence check
} # end of buffer-less points presence check
- return(temp)
+ return(terra::wrap(temp))
}
)
# --------------
# mosaic rasters
-names_metrics <- names(metrics[[1]])
-metrics <- do.call(raster::merge, metrics)
-names(metrics) <- names_metrics
+# unpack rasters
+metrics <- lapply(metrics, terra::rast)
+# names_metrics <- names(metrics[[1]])
+metrics <- do.call(terra::merge, metrics)
+# names(metrics) <- names_metrics
# --------------------------
# compute additional metrics
# Simpson index
-metrics$Hsimpson <- raster::stackApply(metrics[[n_breaks_h[c(-1,-length(n_breaks_h))]]], 1,
+metrics$Hsimpson <- terra::tapp(metrics[[n_breaks_h[c(-1,-length(n_breaks_h))]]], 1,
function(x, ...) {
vegan::diversity(x, index = "simpson")
})
@@ -768,7 +764,7 @@ cumu_sum <- metrics[["nb_ground"]]
for (i in n_breaks_h)
{
cumu_sum <-
- raster::addLayer(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics[[i]])
+ c(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics[[i]])
}
names(cumu_sum) <- c("nb_ground", n_breaks_h)
# interception ratio
@@ -779,7 +775,7 @@ for (i in 2:dim(cumu_sum)[3])
}
names(intercep_ratio) <- gsub("nb", "intercepRatio", names(cumu_sum)[-1])
# merge rasterstacks
-metrics <- raster::addLayer(metrics, intercep_ratio)
+metrics <- c(metrics, intercep_ratio)
#
# ----------------------
# export as raster files
@@ -790,7 +786,7 @@ metrics <- raster::addLayer(metrics, intercep_ratio)
# }
# -------
# display
-raster::plot(metrics[[c("Imean_1stpulse", "Hsimpson")]],
+terra::plot(metrics[[c("Imean_1stpulse", "Hsimpson")]],
main = c("Mean intensity of 1st pulse",
"Simpson indice of point heights")
)
diff --git a/R/gaps.edges.detection.Rmd b/R/gaps.edges.detection.Rmd
index ded196db2e38348da2c3976a8773984f5a1b33fd..81a19fd56496d72ef218e606bb5894be87aa6993 100755
--- a/R/gaps.edges.detection.Rmd
+++ b/R/gaps.edges.detection.Rmd
@@ -33,6 +33,7 @@ The study area is a 200m x 200m forest located in the Jura mountain (France). T
```{r extractals, eval=FALSE}
# build catalog from folder with normalized LAS/LAZ files
cata <- lidR::readALSLAScatalog("/las.folder/")
+# "/media/reseau/lessem/ProjetsCommuns/Lidar/data/Franche-Comte/norm.laz/"
# set coordinate system
lidR::projection(cata) <- 2154
# option to read only xyzc attributes (coordinates, intensity, echo order and classification) from height files
@@ -46,8 +47,8 @@ las <- lidR::clip_rectangle(
centre[2] + size
)
# compute canopy height model from normalized point cloud
-chm <- lidaRtRee::points2DSM(
- las, 1, centre[1] - size, centre[1] + size, centre[2] - size,
+chm <- lidR::rasterize_canopy(
+ las, 1, terra::ext(centre[1] - size, centre[1] + size, centre[2] - size,
centre[2] + size
)
# save chm in Rdata file
@@ -58,22 +59,22 @@ If the CHM has been previously calculated.
```{r loadals, include = TRUE}
# load dataset
load(file = "../data/gap.detection/chm.rda")
+chm <- terra::rast(chm)
```
Missing, low and high values are replaced for better visualization and processing.
-```{r checkchm, include = TRUE, out.width='70%', fig.width = 9, fig.height = 4.5}
-# replace NA by 0 # raster::values to avoid error in knitr, otherwise use:
-# chm[is.na(chm[])] <- 0
-raster::values(chm)[is.na(raster::values(chm))] <- 0
+```{r checkchm, include = TRUE, out.width='70%', fig.width = 9, fig.height = 3.5}
+# replace NA by 0
+chm[is.na(chm)] <- 0
# replace negative values by 0
-chm[chm[] < 0] <- 0
+chm[chm < 0] <- 0
# set maximum threshold at 50m
chm[chm > 50] <- 50
# display CHM and areas where vegetation height is lower than 1m
par(mfrow = c(1, 2))
-raster::plot(chm, main = "CHM")
-raster::plot(chm < 1, asp = 1, main = "CHM < 1m", legend = FALSE)
+terra::plot(chm, main = "CHM (m)")
+terra::plot(chm < 1, asp = 1, main = "CHM < 1m")
```
## Gap detection
@@ -93,18 +94,18 @@ gaps2 <- lidaRtRee::gap_detection(chm, ratio = NULL, gap_max_height = 1, min_gap
```
The following figure displays the obtained gaps, colored by size.
-```{r gapsplot, include = TRUE, out.width='70%', fig.width = 9, fig.height = 4.5}
+```{r gapsplot, include = TRUE, out.width='70%', fig.width = 9, fig.height = 3.5}
# display CHM and areas where vegetation height is lower than 1m
par(mfrow = c(1, 2))
# max value of surface (log)
-z_lim <- log(max(raster::values(gaps1$gap_surface)))
+z_lim <- log(max(terra::values(gaps1$gap_surface), na.rm = TRUE))
# plot gap surface
-raster::plot(log(gaps1$gap_surface),
- main = "log(Gap surface) (1)", zlim = c(0, z_lim),
+terra::plot(log(gaps1$gap_surface),
+ main = "log(Gap surface) (1)", range = c(0, ceiling(z_lim)),
col = rainbow(255, end = 5 / 6)
)
-raster::plot(log(gaps2$gap_surface),
- main = "log(Gap surface) (2)", zlim = c(0, z_lim),
+terra::plot(log(gaps2$gap_surface),
+ main = "log(Gap surface) (2)", range = c(0, ceiling(z_lim)),
col = rainbow(255, end = 5 / 6)
)
```
@@ -113,20 +114,20 @@ Gaps can be vectorized for display over the original CHM.
```{r gapsvectorize, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5)}
# convert to vector
-gaps1_v <- raster::rasterToPolygons(gaps1$gap_id, dissolve = TRUE)
+gaps1_v <- terra::as.polygons(gaps1$gap_id)
# display gaps over original CHM
-raster::plot(chm, main = "Gaps over CHM (1)")
-sp::plot(gaps1_v, add = T)
+terra::plot(chm, main = "Gaps over CHM (1)")
+terra::plot(gaps1_v, add = T)
```
Statistics on gaps (number, size) can be derived from the raster objects. Example for gaps (2)
```{r gapsstats, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5) }
# compute gaps surface
-gaps_df <- data.frame(t(table(raster::values(gaps2$gap_id))))[, -1]
+gaps_df <- data.frame(t(table(terra::values(gaps2$gap_id))))[, -1]
# change column names
names(gaps_df) <- c("id", "surface")
# convert surface in pixels to square meters
-gaps_df$surface <- gaps_df$surface * raster::res(chm)[1]^2
+gaps_df$surface <- gaps_df$surface * terra::xres(chm)^2
# data.frame
head(gaps_df, n = 3)
# cumulative distribution of gap surface by surface
@@ -137,8 +138,8 @@ Histogram can be used to compare the distribution of surface in different classe
```{r gapsstats2, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5) }
# extract vector of total gap surface of pixel
-surface1 <- raster::values(gaps1$gap_surface)
-surface2 <- raster::values(gaps2$gap_surface)
+surface1 <- terra::values(gaps1$gap_surface)
+surface2 <- terra::values(gaps2$gap_surface)
# remove large gaps
surface1[surface1 > 500] <- NA
surface2[surface2 > 500] <- NA
@@ -152,9 +153,9 @@ legend("topright", c("1", "2"), fill = c("red", "blue"))
In the previous part, a pixel that fulfills the maximum height criterion but does not fulfill the distance ratio criterion (i.e. it is too close to surrounding vegetation based on the distance / vegetation height ratio) is removed from gap surface. For edge detection, it seems more appropriate to integrate those pixels in the gap surface, otherwise some edges would be detected inside flat areas. The difference is exemplified in the following plots. The second image exhibits more gaps, because some gaps which are not reconstructed do not reach the minimum surface. The gaps in the second image are also larger, because gaps extend to all neighboring pixels complying with the height threshold, even if they do not comply with the distance criterion.
-```{r gapsReconstruct, include = TRUE, out.width='100%', fig.width = 13, fig.height = 4.5, warning = FALSE}
+```{r gapsReconstruct, include = TRUE, out.width='100%', fig.width = 12, fig.height = 3.5, warning = FALSE}
# load canopy height model
-chm <- lidaRtRee::chm_chablais3
+chm <- terra::rast(lidaRtRee::chm_chablais3)
chm[is.na(chm)] <- 0
# Perform gap detection with distance ratio criterion
gaps1 <- lidaRtRee::gap_detection(chm,
@@ -169,9 +170,9 @@ gaps1r <- lidaRtRee::gap_detection(chm,
# display detected gaps
par(mfrow = c(1, 3))
# plot gaps
-raster::plot(chm, main = "Canopy heightmodel")
-raster::plot(gaps1$gap_id > 0, main = "Gaps", legend = FALSE)
-raster::plot(gaps1r$gap_id > 0, main = "Gaps extended by reconstruction", legend = FALSE)
+terra::plot(chm, main = "Canopy heightmodel")
+terra::plot(gaps1$gap_id > 0, main = "Gaps")
+terra::plot(gaps1r$gap_id > 0, main = "Gaps extended by reconstruction")
```
Edge detection is performed by extraction of the difference between a binary image of gaps and the result of a morphological erosion or dilation applied to the same image.
@@ -183,9 +184,9 @@ edge_inside <- lidaRtRee::edge_detection(gaps1r$gap_id > 0)
edge_outside <- lidaRtRee::edge_detection(gaps1r$gap_id > 0, inside = FALSE)
# display zoom on edges
par(mfrow = c(1, 2))
-raster::plot(edge_inside, main = "Edges (detection by erosion)", legend = FALSE)
-raster::plot(edge_outside, main = "Edges (detection by dilation)", legend = FALSE)
+terra::plot(edge_inside, main = "Edges (detection by erosion)", legend = FALSE)
+terra::plot(edge_outside, main = "Edges (detection by dilation)", legend = FALSE)
```
-Percentage of edges can be calculated as the ratio between edge pixels surface and total surface, multiplied by 100. The result is dependent on the raster resolution. Obtained percentage for the "erosion" method is `r round(sum(raster::values(edge_inside))/(nrow(edge_inside)*ncol(edge_inside))*100, 1)`, whereas it is `r round(sum(raster::values(edge_outside))/(nrow(edge_outside)*ncol(edge_outside))*100, 1)` for method "dilation".
+Percentage of edges can be calculated as the ratio between edge pixels surface and total surface, multiplied by 100. The result is dependent on the raster resolution. Obtained percentage for the "erosion" method is `r round(sum(terra::values(edge_inside))/(nrow(edge_inside)*ncol(edge_inside))*100, 1)`, whereas it is `r round(sum(terra::values(edge_outside))/(nrow(edge_outside)*ncol(edge_outside))*100, 1)` for method "dilation".
## References
\ No newline at end of file
diff --git a/R/tree.detection.Rmd b/R/tree.detection.Rmd
index 9c21b499ea5a6c90253a581e899b9d1d6d1bb5aa..486913624d8a51efcfcbd1463a8ee86991b2424e 100755
--- a/R/tree.detection.Rmd
+++ b/R/tree.detection.Rmd
@@ -259,7 +259,7 @@ head(apices, n = 3L)
segments_v <- terra::as.polygons(segms$segments_id)
#
# display initial image
-raster::plot(chm, col = gray(seq(0, 1, 1 / 255)), main = "CHM and detected positions")
+terra::plot(chm, col = gray(seq(0, 1, 1 / 255)), main = "CHM and detected positions")
# display segments border
terra::plot(segments_v, border = "white", add = T)
# display plot mask
@@ -391,11 +391,11 @@ Plotting the reference trees with the mean intensity of lidar points in the segm
# create raster of segment' mean intensity
r_mean_intensity_segm <- segms$segments_id
# match segment id with id in metrics data.frame
-dummy <- match(raster::values(r_mean_intensity_segm), apices$id)
+dummy <- match(terra::values(r_mean_intensity_segm), apices$id)
# replace segment id with corresponding mean intensity
-raster::values(r_mean_intensity_segm) <- apices$meanI[dummy]
+terra::values(r_mean_intensity_segm) <- apices$meanI[dummy]
# display tree inventory with mean intensity in segment background
-raster::plot(r_mean_intensity_segm, main = "Mean intensity in segment")
+terra::plot(r_mean_intensity_segm, main = "Mean intensity in segment")
# display tree inventory
lidaRtRee::plot_tree_inventory(tree_inventory_chablais3[, c("x", "y")],
tree_inventory_chablais3$h,
@@ -532,7 +532,7 @@ To display the classification results, an image of the classified segments is cr
species <- segms$segments_id
# replace segment id by predicted species in image
terra::values(species) <-
- metrics$predicted_s[match(raster::values(segms$segments_id), metrics$seg_id)]
+ metrics$predicted_s[match(terra::values(segms$segments_id), metrics$seg_id)]
# remove ground segment
species[segms$segments_id == 0] <- NA
# convert to factor raster
@@ -550,7 +550,7 @@ levels(species)[[1]] <- rat
# display results
terra::plot(species, col = rat$col, legend = FALSE)
legend("bottomright", legend=rat$Species, fill=rat$col)
-sp::plot(segments_v, add = TRUE, border = "white")
+terra::plot(segments_v, add = TRUE, border = "white")
lidaRtRee::plot_tree_inventory(tree_metrics_h
[, c("x", "y")],
tree_metrics_h$h,
@@ -588,11 +588,11 @@ The point cloud can be displayed colored by segment, with poles at the location
```{r displayPointCloud, include=TRUE, eval=html, webgl=TRUE, fig.width=6, fig.height=6, warning=FALSE}
rgl::par3d(mouseMode = "trackball") # parameters for interaction with mouse
# select segment points and offset them to avoid truncated coordinates in 3d plot
-points.seg <- lasn@data[which(lasn@data$seg_id != 0), c("X", "Y", "Z", "seg_id")]
+points.seg <- lasn[which(lasn$seg_id != 0), c("X", "Y", "Z", "seg_id")]
points.seg$X <- points.seg$X - 974300
points.seg$Y <- points.seg$Y - 6581600
# plot point cloud
-rgl::plot3d(points.seg[, c("X", "Y", "Z")], col = points.seg$seg_id %% 10 + 1, aspect = FALSE)
+rgl::plot3d(points.seg@data[, c("X", "Y", "Z")], col = points.seg$seg_id %% 10 + 1, aspect = FALSE)
#
# add inventoried trees
tree_inventory_chablais3$z <- 0
@@ -855,15 +855,11 @@ resultats <- future.apply::future_mapply(
apices_v <- terra::as.polygons(segms$segments_id)
# remove polygons which treetop is in buffer
apices_v <- apices_v[is.element(apices_v$segments_id, apices$id), ]
- # names(apices_v) <- "id"
- # add attributes
- # errors when using sp::merge so using sp::over even if it is probably slower
# convert to sf
apices_v <- sf::st_as_sf(apices_v)
names(apices_v)[1] <- "id"
- apices_v <- merge(apices_v, sf::st_drop_geometry(apices))
- # merge(apices_v@data, apices, all.x = TRUE)
- #apices_v@data <- cbind(apices_v@data, sp::over(apices_v, apices))
+ # add attributes
+ apices_v <- merge(apices_v, sf::st_drop_geometry(apices), all.x = TRUE)
# save in list
output$apices_v <- apices_v
}
@@ -888,6 +884,7 @@ apices <- do.call(rbind, apices)
#
# chm
if (out_chm) {
+ # extract and unwrap chms
chm <- lapply(resultats, function(x) terra::rast(x[["chm"]]))
# merge chm
# no names in list otherwise do.call returns an error
@@ -927,9 +924,9 @@ The following lines save outputs to files.
# merged chm
terra::writeRaster(chm_all, file = "chm.tif", overwrite = TRUE)
# apices
-sf::st_write(apices, "apices_points.gpkg", delete_layer = TRUE)
+sf::st_write(apices, "apices_points.gpkg")# , delete_layer = TRUE)
# vectorized apices
-if (out_vectorize_apices) sf::st_write(apices_v, "v_apices_points.gpkg", delete_layer = TRUE)
+if (out_vectorize_apices) sf::st_write(apices_v, "v_apices_points.gpkg")#, delete_layer = TRUE)
```
## References
\ No newline at end of file
diff --git a/README.md b/README.md
index 7a52b1ba871e58269e02b621fd66b4379451242f..d9dfed08ebbfc19748ffad323f4ce49eab6261a6 100755
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
# `lidaRtRee` tutorials
-`lidaRtRee_tutorials` is a repository that provides tutorials for forest analysis with airborne laser scanning (ALS or lidar remote sensing) data, using functions from `R` packages [lidR](https://github.com/Jean-Romain/lidR/) (also available on [CRAN](https://cran.r-project.org/package=lidR)) and [lidaRtRee](https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee) (also available on [CRAN](https://cran.r-project.org/package=lidaRtRee)). Tutorials are available as `Rmarkdown` and `html` files. Datasets useful for code running are also provided.
+`lidaRtRee_tutorials` is a repository that provides tutorials for forest analysis with airborne laser scanning (ALS or lidar remote sensing) data, using functions from `R` packages [lidR](https://github.com/r-lidar/lidR/) (also available on [CRAN](https://cran.r-project.org/package=lidR)) and [lidaRtRee](https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee) (also available on [CRAN](https://cran.r-project.org/package=lidaRtRee)). Tutorials are available as `Rmarkdown` and `html` files. Datasets useful for code running are also provided.
# Available tutorials
@@ -8,4 +8,4 @@ They are listed on the [wiki home page](https://gitlab.irstea.fr/jean-matthieu.m
# Acknowledgements
-`lidaRtRee_tutorials` development was funded by the [ADEME](https://www.ademe.fr/en) (french Agency for Ecological Transition) through the [PROTEST](https://protest.inrae.fr/) project (grant 1703C0069 of the GRAINE program).
+`lidaRtRee_tutorials` development (2018-2021) was funded by the [ADEME](https://www.ademe.fr/en) (french Agency for Ecological Transition) through the [PROTEST](https://protest.inrae.fr/) project (grant 1703C0069 of the GRAINE program).