Commit 195d658a authored by Monnet Jean-Matthieu's avatar Monnet Jean-Matthieu
Browse files

WIP 4.0.0

No related merge requests found
Showing with 315 additions and 297 deletions
+315 -297
......@@ -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)")
......
......@@ -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.
......
......@@ -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)
......
......@@ -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)
```
......
......@@ -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
......
This diff is collapsed.
......@@ -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
......@@ -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
# `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).
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