diff --git a/R/area-based.1.data.preparation.Rmd b/R/area-based.1.data.preparation.Rmd
index 5e99a9c7cfae57c07a389580a1111a4440cdde60..ee2a9b096868919df5a863568a2ec236a0867de2 100644
--- a/R/area-based.1.data.preparation.Rmd
+++ b/R/area-based.1.data.preparation.Rmd
@@ -3,8 +3,8 @@ title: "R workflow for ABA data preparation"
author: "Jean-Matthieu Monnet"
date: "`r Sys.Date()`"
output:
- pdf_document: default
html_document: default
+ pdf_document: default
papersize: a4
bibliography: "../bib/bibliography.bib"
---
@@ -44,61 +44,60 @@ The first step is to import tree and plot information in two separate data.frame
```{r treeInventory, include=TRUE}
# set inventory parameters
# plot radius (m)
-p.radius <- 15
+p_radius <- 15
# DBH threshold (cm) for inventory, trees smaller than this value are not inventoried
-dbh.min <- 7.5
+dbh_min <- 7.5
# list tree files by pattern matching
-files.t <- as.list(dir(
+files_t <- as.list(dir(
path = "../data/aba.model/field/",
pattern = "Verc-[[:alnum:]]{2}-[[:digit:]]{1}_ArbresTerrain.csv",
full.names = TRUE
))
# load content of all files with lapply
-trees <- lapply(files.t, function(x) {
+trees <- lapply(files_t, function(x) {
# read table
dummy <- read.table(x, sep = ";", header = T, stringsAsFactors = F)
- # add one column with plotId from file name
- cbind(dummy, data.frame(plotId = rep(substr(x, 31, 34), nrow(dummy))))
+ # add one column with plot_id from file name
+ cbind(dummy, data.frame(plot_id = rep(substr(x, 31, 34), nrow(dummy))))
})
# bind elements of list into a single data.frame
trees <- do.call(rbind, trees)
-# add study area info in plotId
-trees$plotId <- paste0("Verc-", trees$plotId)
+# add study area info in plot_id
+trees$plot_id <- paste0("Verc-", trees$plot_id)
```
Fields are:
-* `treeId`: tree id in the plot
-+ `poleId`: plot id inside the cluster (1 to 4)
-+ `Species`: tree species abreviated as GESP (GEnus SPecies)
-+ `Azimuth.gr`: azimuth in grades from the plot center to tree center
-+ `Slope.gr`: slope in grades from the plot center to tree center
-+ `Diameter.cm`: tree diameter at breast height (1.3 m, upslope), in cm
-+ `Ground.Distance.m`: ground distance between the plot center and tree edge at breast height, in m
-+ `Height.m`: tree height, in m
-+ `Appearance`:
+* `tree_id`: tree id in the plot
++ `pole_id`: plot id inside the cluster (1 to 4)
++ `species`: tree species abreviated as GESP (GEnus SPecies)
++ `azimuth_gr`: azimuth in grades from the plot center to tree center
++ `slope_gr`: slope in grades from the plot center to tree center
++ `diameter_cm`: tree diameter at breast height (1.3 m, upslope), in cm
++ `ground_distance_m`: ground distance between the plot center and tree edge at breast height, in m
++ `height_m`: tree height, in m
++ `appearance`:
+ 0: lying or missing,
+ 1: live,
+ 2: live with broken treetop,
+ 3: dead with branches,
+ 4: dead without branches (snag)
-+ `Tilted`: 0: no, 1:yes
-+ `Remark`
-+ `plotId`: plot id
-
-The following lines set column names to English and convert the `Appearance` column to factor with adequate values for the levels.
++ `tilted`: 0: no, 1:yes
++ `remark`
++ `plot_id`: plot id
+The following lines set column names to English and convert the `appearance` column to factor with adequate values for the levels.
```{r treeInventory2, include=TRUE}
# change column names to English
names(trees) <- c(
- "treeId", "poleId", "Species", "Azimuth.gr", "Slope.gr", "Diameter.cm",
- "Ground.Distance.m", "Height.m", "Appearance", "Tilted", "Remark", "plotId"
+ "tree_id", "pole_id", "species", "azimuth_gr", "slope_gr", "diameter_cm",
+ "ground_distance_m", "height_m", "appearance", "tilted", "remark", "plot_id"
)
-# add factor levels to column Appearance
-trees$Appearance <- factor(trees$Appearance, levels = 0:4)
-levels(trees$Appearance) <- c("missing or lying", "live", "live with broken treetop", "dead with branches", "dead without branches (snag)")
+# add factor levels to column appearance
+trees$appearance <- factor(trees$appearance, levels = 0:4)
+levels(trees$appearance) <- c("missing or lying", "live", "live with broken treetop", "dead with branches", "dead without branches (snag)")
# convert species column to factor
-trees$Species <- factor(trees$Species)
+trees$species <- factor(trees$species)
head(trees, n = 3)
```
@@ -108,7 +107,7 @@ The next step is to import plot coordinates.
```{r plotLocation, include=TRUE, fig.width = 4, fig.height = 6}
# list location files by pattern matching
-files.p <- dir(
+files_p <- dir(
path = "../data/aba.model/field/",
pattern = "Verc-[[:alnum:]]{2}_PiquetsTerrain.csv",
full.names = TRUE
@@ -116,26 +115,26 @@ files.p <- dir(
# initialize data.frame
plots <- NULL
# load all plot files with a loop
-for (i in files.p)
+for (i in files_p)
{
# read file
dummy <- read.table(i, sep = ";", header = T, stringsAsFactors = F)
- # append to data.frame and add plotId info
- plots <- rbind(plots, cbind(dummy, data.frame(plotId = rep(substr(i, 31, 32), nrow(dummy)))))
+ # append to data.frame and add plot_id info
+ plots <- rbind(plots, cbind(dummy, data.frame(plot_id = rep(substr(i, 31, 32), nrow(dummy)))))
}
# keep only necessary data in plots data.frame (remove duplicate position measurements)
plots <- plots[is.element(plots$Id, c("p1", "p2", "p3", "p4")), ]
-# add plotId to clusterId
-plots$clusterId <- paste0("Verc-", substr(plots$plotId, 1, 2))
-plots$plotId <- paste(plots$clusterId, substr(plots$Id, 2, 2), sep = "-")
+# add plot_id to cluster_id
+plots$cluster_id <- paste0("Verc-", substr(plots$plot_id, 1, 2))
+plots$plot_id <- paste(plots$cluster_id, substr(plots$Id, 2, 2), sep = "-")
# keep only coordinates and Id in data.frame
-plots <- plots[, c("X", "Y", "plotId", "clusterId")]
+plots <- plots[, c("X", "Y", "plot_id", "cluster_id")]
# convert to spatial sf object
-plots.sf <- sf::st_as_sf(plots, coords = c("X", "Y"))
+plots_sf <- sf::st_as_sf(plots, coords = c("X", "Y"))
# set coordinate system
-sf::st_crs(plots.sf) <- 2154
+sf::st_crs(plots_sf) <- 2154
# add coordinates in data.frame of sf object
-plots.sf <- cbind(plots.sf, data.frame(sf::st_coordinates(plots.sf)))
+plots_sf <- cbind(plots_sf, data.frame(sf::st_coordinates(plots_sf)))
```
The following lines maps the plot locations with a background from OpenStreetMap :
@@ -146,9 +145,9 @@ The following lines maps the plot locations with a background from OpenStreetMap
```{r plotLocation.osm, include=TRUE, fig.width = 4, fig.height = 6, warning=FALSE, message=FALSE}
# transform coordinates to lat / lon
-plots.sf.transform <- sf::st_transform(plots.sf, 4326)
+plots_sf_transform <- sf::st_transform(plots_sf, 4326)
# extract bounding box
-ext <- sf::st_bbox(plots.sf.transform)
+ext <- sf::st_bbox(plots_sf_transform)
# set buffer
buffer <- 0.02
# define location extent
@@ -160,7 +159,7 @@ map <- ggmap::get_map(location = loc, source = "osm", zoom = 11)
# draw map
ggmap::ggmap(map) +
ggplot2::labs(title = "Field plot location", x = "Longitude", y = "Latitude") +
- ggplot2::geom_sf(data = plots.sf.transform, inherit.aes = FALSE)
+ ggplot2::geom_sf(data = plots_sf_transform, inherit.aes = FALSE)
```
In case tree positions have to be precisely calculated in a given projected coordinates system, it is necessary to correct magnetic azimuth with magnetic declination at the time of inventory and meridian convergence at the location. The values of magnetic declination and meridian convergence at each plot location are loaded from an additional table.
@@ -171,28 +170,28 @@ meta <- read.table(file = "../data/aba.model/field/plot.metadata.csv", sep = ";"
# keep only required attributes
meta <- meta[, c("Id", "Convergence_gr", "Declinaison_gr")]
# rename attributes
-names(meta) <- c("clusterId", "Convergence.gr", "Declination.gr")
+names(meta) <- c("cluster_id", "convergence_gr", "declination_gr")
# merge to add convergence and declination in plots data.frame
-plots.sf <- merge(plots.sf, meta)
-head(plots.sf, n = 3)
+plots_sf <- merge(plots_sf, meta)
+head(plots_sf, n = 3)
```
### Compute tree coordinates
Tree coordinates can be computed from the plot center coordinates and from the azimuth, slope and ground distance measurements. In case ground distance is measured to the tree edge, the tree diameter has to be taken into account to compute the position of the tree center.
```{r treeCoordinates, include=TRUE}
# merge tree and plots data.frames to import plot center
-trees <- merge(trees, plots.sf[, c("X", "Y", "Convergence.gr", "Declination.gr", "plotId")],
- by = "plotId"
+trees <- merge(trees, plots_sf[, c("X", "Y", "convergence_gr", "declination_gr", "plot_id")],
+ by = "plot_id"
)
# compute projected coordinates
dummy <- lidaRtRee::polar2Projected(
- trees$X, trees$Y, 0, trees$Azimuth.gr / 200 * pi,
- trees$Ground.Distance.m, trees$Slope.gr / 200 * pi,
- trees$Declination.gr / 200 * pi,
- trees$Convergence.gr / 200 * pi, trees$Diameter.cm / 100
+ trees$X, trees$Y, 0, trees$azimuth_gr / 200 * pi,
+ trees$ground_distance_m, trees$slope_gr / 200 * pi,
+ trees$declination_gr / 200 * pi,
+ trees$convergence_gr / 200 * pi, trees$diameter_cm / 100
)
# add coordinates to trees data.frame
-trees[, c("X", "Y", "Horiz.Distance.m")] <- dummy[, c("x", "y", "d")]
+trees[, c("X", "Y", "horiz_dist_m")] <- dummy[, c("x", "y", "d")]
head(trees, n = 3)
```
@@ -215,26 +214,26 @@ lidR::projection(cata) <- 2154
# disable display of catalog processing
lidR::opt_progress(cata) <- FALSE
# extract from LAZ files in catalog
-llas.height <- lidR::clip_circle(cata, plots.sf$X, plots.sf$Y, p.radius + 5)
+llas_height <- lidR::clip_circle(cata, plots_sf$X, plots_sf$Y, p_radius + 5)
# add "buffer" attribute equal to TRUE to points outside plot
-for (i in 1:length(llas.height))
+for (i in 1:length(llas_height))
{
- llas.height[[i]] <- lidR::add_attribute(
- llas.height[[i]],
- (llas.height[[i]]$X - plots.sf$X[i])^2 +
- (llas.height[[i]]$Y - plots.sf$Y[i])^2 > p.radius^2,
+ llas_height[[i]] <- lidR::add_attribute(
+ llas_height[[i]],
+ (llas_height[[i]]$X - plots_sf$X[i])^2 +
+ (llas_height[[i]]$Y - plots_sf$Y[i])^2 > p_radius^2,
"buffer"
)
}
# set names of point clouds in list
-names(llas.height) <- plots.sf$plotId
+names(llas_height) <- plots_sf$plot_id
# set negative height values to 0
-for (i in 1:length(llas.height)) {
- llas.height[[i]]@data$Z[llas.height[[i]]@data$Z < 0] <- 0
+for (i in 1:length(llas_height)) {
+ llas_height[[i]]@data$Z[llas_height[[i]]@data$Z < 0] <- 0
}
# write plot clouds
-# for (i in names(llas.height))
-# {lidR::writeLAS(llas.height[[i]], file=paste0(lazdir, i, ".laz"))}
+# for (i in names(llas_height))
+# {lidR::writeLAS(llas_height[[i]], file=paste0(lazdir, i, ".laz"))}
```
Point clouds with altitude values can be used to compute terrain statistics. The folder `plots.laz` contains the point clouds with altitude value extracted on the extent of each field plot. No buffer is added when loading those point clouds. The points not classified as ground are removed.
@@ -249,34 +248,34 @@ lidR::projection(cata) <- 2154
# disable plot of tile processing
cata@processing_options$progress <- FALSE
# extract las point cloud without buffer
-llas.ground <- lidR::clip_circle(cata, plots$X, plots$Y, p.radius)
-names(llas.ground) <- plots$plotId
+llas_ground <- lidR::clip_circle(cata, plots$X, plots$Y, p_radius)
+names(llas_ground) <- plots$plot_id
# keep only ground points in point clouds
-llas.ground <- lapply(llas.ground, lidR::filter_ground)
+llas_ground <- lapply(llas_ground, lidR::filter_ground)
```
### Computation of terrain statistics
-Terrain statistics can be extracted from the cloud of ground points with altitude values. The function `points2terrainStats` computes the aspect, altitude and slope of a point cloud by fitting a plane to points. The altitude is extracted by interpolating values at user-specified coordinates, or as the mean of the range of altitude values if no coordinates are specified.
+Terrain statistics can be extracted from the cloud of ground points with altitude values. The function `terrain_points_metrics` computes the aspect, altitude and slope of a point cloud by fitting a plane to points. The altitude is extracted by interpolating values at user-specified coordinates, or as the mean of the range of altitude values if no coordinates are specified.
```{r terrainStats, include=TRUE, warning=FALSE, message=FALSE, fig.width = 4, fig.height = 6}
# use mapply to apply points2terrainStats function to each point cloud
# while providing the coordinates of each center
-metrics.terrain <- mapply(function(x, y) {
- lidaRtRee::points2terrainStats(x, plots[y, c("X", "Y")], p.radius)
+metrics_terrain <- mapply(function(x, y) {
+ lidaRtRee::terrain_points_metrics(x, plots[y, c("X", "Y")], p_radius)
},
-x = llas.ground,
-y = as.list(1:length(llas.ground)),
+x = llas_ground,
+y = as.list(1:length(llas_ground)),
SIMPLIFY = FALSE
)
# bind results in data.frame
-metrics.terrain <- do.call(rbind, metrics.terrain)
+metrics_terrain <- do.call(rbind, metrics_terrain)
#
# compute terrain stats without specifying centre
-metrics.terrain2 <- lapply(llas.ground, lidaRtRee::points2terrainStats)
-metrics.terrain2 <- do.call(rbind, metrics.terrain2)
+metrics_terrain2 <- lapply(llas_ground, lidaRtRee::terrain_points_metrics)
+metrics_terrain2 <- do.call(rbind, metrics_terrain2)
# display results for comparison
-head(cbind(metrics.terrain, metrics.terrain2), n = 3)
+head(cbind(metrics_terrain, metrics_terrain2), n = 3)
```
## Check field inventory data
@@ -291,19 +290,19 @@ To ensure that all trees were correctly inventoried according to the protocol :
- Diameters can also be checked to avoid that trees below the DBH limit remain in the inventory.
```{r checkData, include=TRUE}
-summary(trees[, c("Diameter.cm", "Horiz.Distance.m")])
+summary(trees[, c("diameter_cm", "horiz_dist_m")])
```
Once values have been checked for potential writing errors, remaining errors should be removed.
```{r checkProtocol, include=TRUE, }
-nb.trees <- nrow(trees)
+nb_trees <- nrow(trees)
# keep only trees inside plot
-trees <- trees[trees$Horiz.Distance.m <= p.radius, ]
+trees <- trees[trees$horiz_dist_m <= p_radius, ]
# keep trees above the DBH limit
-trees <- trees[trees$Diameter.cm >= dbh.min, ]
+trees <- trees[trees$diameter_cm >= dbh_min, ]
# number of removed trees
-# nb.trees - nrow(trees)
+# nb_trees - nrow(trees)
```
-In this case `r nb.trees - nrow(trees)` trees were removed.
+In this case `r nb_trees - nrow(trees)` trees were removed.
### Consistency of trees azimuth and slope to center (by plot)
@@ -311,12 +310,12 @@ If tree 3D spherical coordinates have been recorded, azimuth can be checked agai
```{r checkSlopeAzimuth, include=TRUE, fig.width = 5, fig.height = 4}
# example plot to test
-plot.test <- "Verc-01-1"
+plot_test <- "Verc-01-1"
# plot slope as a function of azimuth,
# symbol size proportional to horizontal distance from center
-plot(Slope.gr ~ Azimuth.gr,
- data = trees, subset = which(trees$plotId == plot.test),
- cex = Horiz.Distance.m / p.radius, main = "Trees slope and azimuth to plot center"
+plot(slope_gr ~ azimuth_gr,
+ data = trees, subset = which(trees$plot_id == plot_test),
+ cex = horiz_dist_m / p_radius, main = "Trees slope and azimuth to plot center"
)
```
@@ -324,9 +323,9 @@ The following lines create a spatial polygon corresponding to one plot extent.
```{r plotExtent, include = TRUE}
# extract example plot
-plots.sf.t <- plots.sf[plots.sf$plotId == plot.test, ]
+plots_sf.t <- plots_sf[plots_sf$plot_id == plot_test, ]
# create polygon with buffer of plot radius
-plots.extent.t <- sf::st_buffer(plots.sf.t, p.radius, nQuadSegs = 10)
+plots_extent_t <- sf::st_buffer(plots_sf.t, p_radius, nQuadSegs = 10)
```
### Consistency of field operator path (by plot)
@@ -335,20 +334,20 @@ In case trees have been numbered in the order of inventory on the field, plottin
```{r checkPath, include=TRUE, fig.width = 5, fig.height = 4}
# extraction of trees of plot to test
-trees.t <- trees[trees$plotId == plot.test, ]
+trees_t <- trees[trees$plot_id == plot_test, ]
# display plot limits
-plot(sf::st_geometry(plots.extent.t), border = "red", axes = TRUE)
+plot(sf::st_geometry(plots_extent_t), border = "red", axes = TRUE)
# display plot center
-plot(sf::st_geometry(plots.sf.t), col = "red", pch = 3, add = TRUE)
+plot(sf::st_geometry(plots_sf.t), col = "red", pch = 3, add = TRUE)
# add inventoried trees
-lidaRtRee::plotTreeInventory(trees.t[, c("X", "Y")],
- trees.t$Height.m,
- species = as.character(trees.t$Species),
- pch = as.numeric(trees.t$Appearance) - 1,
+lidaRtRee::plot_tree_inventory(trees_t[, c("X", "Y")],
+ trees_t$height_m,
+ species = as.character(trees_t$species),
+ pch = as.numeric(trees_t$appearance) - 1,
add = TRUE
)
# draw lines between trees
-lines(trees.t[order(trees.t$treeId), c("X", "Y")])
+lines(trees_t[order(trees_t$tree_id), c("X", "Y")])
```
### Height / diameter allometry (all trees)
@@ -357,15 +356,15 @@ The allometry of trees can also be checked in case heights were also measured. P
```{r checkAllometry, include=TRUE, fig.width = 6, fig.height = 4.5}
# symbol size proportional to horizontal distance from center
-plot(Height.m ~ Diameter.cm,
+plot(height_m ~ diameter_cm,
data = trees,
- col = lidaRtRee::speciesColor()[as.character(trees$Species), "col"],
- pch = as.numeric(trees$Appearance) - 1,
+ col = lidaRtRee::species_color()[as.character(trees$species), "col"],
+ pch = as.numeric(trees$appearance) - 1,
main = "Height VS diameter, all trees"
)
legend("bottomright",
- legend = levels(trees$Appearance),
- pch = (1:length(levels(trees$Appearance)) - 1)
+ legend = levels(trees$appearance),
+ pch = (1:length(levels(trees$appearance)) - 1)
)
```
@@ -374,14 +373,14 @@ The `ggplot2` package provides useful syntax and functions for producing informa
```{r checkAllometry.ggplot2, include=TRUE, fig.width = 12, fig.height = 5, warning=FALSE, message = FALSE}
g <- ggplot(
- trees[is.element(trees$Species, c("ABAL", "PIAB", "FASY")), ],
- aes(x = Diameter.cm, y = Height.m, shape = Appearance)
+ trees[is.element(trees$species, c("ABAL", "PIAB", "FASY")), ],
+ aes(x = diameter_cm, y = height_m, shape = appearance)
)
-g + geom_point(aes(color = Species)) +
+g + geom_point(aes(color = species)) +
scale_x_sqrt() +
- geom_smooth(aes(group = Species), color = "black", linetype = "dashed") +
- facet_grid(. ~ Species) +
- labs(x = "Diameter.cm (sqrt scale)")
+ geom_smooth(aes(group = species), color = "black", linetype = "dashed") +
+ facet_grid(. ~ species) +
+ labs(x = "diameter_cm (sqrt scale)")
```
### Diameter distribution (by plot)
@@ -391,35 +390,35 @@ For each plot, displaying the histogram of trees, colored by species is also inf
```{r plotDiaDistribution, include=TRUE, fig.width = 5, fig.height = 5, warning=FALSE, message = FALSE}
# plot diameter distribution
# remove unused levels in species factor column
-trees.t$Species <- factor(trees.t$Species)
+trees_t$species <- factor(trees_t$species)
# stacked histogram : data.frame split by species
-dummy <- split(trees.t$Diameter.cm, trees.t$Species)
-lidaRtRee::histStack(dummy,
- col = lidaRtRee::speciesColor()[names(dummy), "col"],
+dummy <- split(trees_t$diameter_cm, trees_t$species)
+lidaRtRee::hist_stack(dummy,
+ col = lidaRtRee::species_color()[names(dummy), "col"],
breaks = seq(
- from = dbh.min, by = 5,
- to = 5 * ceiling(max(trees.t$Diameter.cm - 2.5) / 5) + 2.5
+ from = dbh_min, by = 5,
+ to = 5 * ceiling(max(trees_t$diameter_cm - 2.5) / 5) + 2.5
),
- main = paste0("Diameter distribution of plot ", plot.test),
+ main = paste0("Diameter distribution of plot ", plot_test),
ylab = "Number of trees",
xlab = "Diameter (cm)",
)
-legend("topright", names(dummy), fill = lidaRtRee::speciesColor()[names(dummy), "col"])
+legend("topright", names(dummy), fill = lidaRtRee::species_color()[names(dummy), "col"])
```
```{r plotDiaDistribution.ggplot2, include=TRUE, fig.width = 4, fig.height = 4, warning=FALSE, message = FALSE}
# plot diameter distribution
# change order of levels to display more abundant species at the bottom
-trees.t$Species <- factor(trees.t$Species, levels = names(sort(table(trees.t$Species))))
+trees_t$species <- factor(trees_t$species, levels = names(sort(table(trees_t$species))))
# load custom species table with associated colors
-table.species <- lidaRtRee::speciesColor()
+table.species <- lidaRtRee::species_color()
# extract species present on the plot
-table.species <- table.species[levels(trees.t$Species), "col"]
+table.species <- table.species[levels(trees_t$species), "col"]
# stacked histogram : data.frame split by species
-ggplot(trees.t, aes(x = Diameter.cm, fill = Species)) +
- geom_histogram(breaks = seq(from = 7.5, to = max(trees.t$Diameter.cm) + 5, by = 5), color = "black") +
+ggplot(trees_t, aes(x = diameter_cm, fill = species)) +
+ geom_histogram(breaks = seq(from = 7.5, to = max(trees_t$diameter_cm) + 5, by = 5), color = "black") +
scale_fill_manual(values = table.species) +
labs(
- title = paste0("Diameter distribution of plot ", plot.test),
+ title = paste0("Diameter distribution of plot ", plot_test),
y = "Number of trees", x = "Diameter (cm)"
)
```
@@ -434,21 +433,21 @@ The tree positions can be plotted with the canopy height model computed with the
```{r plotRS, include = TRUE, out.width = '70%', fig.dim=c(6.5, 4.5), warning=FALSE}
# compute background
# compute canopy height model
-chm <- lidaRtRee::points2DSM(llas.height[[plot.test]], res = 0.5)
+chm <- lidaRtRee::points2DSM(llas_height[[plot_test]], res = 0.5)
# display CHM
raster::plot(chm,
col = gray(seq(0, 1, 1 / 255)),
main = "Canopy Height Model and tree positions"
)
# add inventoried trees
-lidaRtRee::plotTreeInventory(trees.t[, c("X", "Y")],
- trees.t$Height.m,
- species = as.character(trees.t$Species),
- pch = as.numeric(trees.t$Appearance) - 1,
+lidaRtRee::plot_tree_inventory(trees_t[, c("X", "Y")],
+ trees_t$height_m,
+ species = as.character(trees_t$species),
+ pch = as.numeric(trees_t$appearance) - 1,
add = TRUE
)
# display plot limits
-plot(plots.extent.t, border = "red", color = NA, add = TRUE)
+plot(plots_extent_t, border = "red", color = NA, add = TRUE)
```
Similar output with `ggplot2`.
@@ -465,19 +464,19 @@ g1 <- ggplot() +
# background scale in black and white
scale_fill_gradientn(colours = c("black", "white"), na.value = "white") +
# add trees
- geom_point(data = trees.t, aes(x = X, y = Y, shape = Appearance, size = Diameter.cm, colour = Species)) +#
+ geom_point(data = trees_t, aes(x = X, y = Y, shape = appearance, size = diameter_cm, colour = species)) +#
# shape radius (not area) proportional to diameter
- scale_radius(name = "Diameter.cm") +
+ scale_radius(name = "diameter_cm") +
# custom colors for species
scale_colour_manual(values = table.species) +
# shapes are empty
scale_shape(solid = FALSE) +
# add plot extent
- geom_sf(data = plots.extent.t, fill = NA, color = "red") +
+ geom_sf(data = plots_extent_t, fill = NA, color = "red") +
# specify coordinate system
coord_sf(datum = 2154) +
# add title
- labs(title = paste0("Map of plot ", plot.test))
+ labs(title = paste0("Map of plot ", plot_test))
print(g1)
```
@@ -487,48 +486,48 @@ print(g1)
Three stand level parameters are computed by aggregating the tree level information at the plot level:
-* basal area (`G.m2.ha`, m^2^/ha)
-+ mean diameter (`D.mean.cm`, cm)
-+ stem density (`N.ha`, /ha).
+* basal area (`G_m2_ha`, m^2^/ha)
++ mean diameter (`D_mean_cm`, cm)
++ stem density (`N_ha`, /ha).
Only live trees are considered in this analysis.
```{r standParameters, include=TRUE, fig.width = 12, fig.height = 4}
# keep only live trees
-trees <- trees[is.element(trees$Appearance, c("live", "live with broken treetop")), ]
+trees <- trees[is.element(trees$appearance, c("live", "live with broken treetop")), ]
# Basal area in m2/ha
-dummy <- aggregate(Diameter.cm ~ plotId, trees,
+dummy <- aggregate(diameter_cm ~ plot_id, trees,
FUN = function(x) {
- sum(pi * (x / 200)^2, na.rm = T) * (10000 / (pi * p.radius^2))
+ sum(pi * (x / 200)^2, na.rm = T) * (10000 / (pi * p_radius^2))
}
)
# rename colum of result
-names(dummy)[2] <- "G.m2.ha"
+names(dummy)[2] <- "G_m2_ha"
# add information to plots data.frame
plots <- merge(plots, dummy)
# Stem density in m2/ha
-dummy <- aggregate(Diameter.cm ~ plotId, trees,
+dummy <- aggregate(diameter_cm ~ plot_id, trees,
FUN = function(x) {
- length(x) * (10000 / (pi * p.radius^2))
+ length(x) * (10000 / (pi * p_radius^2))
}
)
-names(dummy)[2] <- "N.ha"
+names(dummy)[2] <- "N_ha"
plots <- merge(plots, dummy)
# Mean diameter in cm
-dummy <- aggregate(Diameter.cm ~ plotId, trees,
+dummy <- aggregate(diameter_cm ~ plot_id, trees,
FUN = function(x) {
mean(x)
}
)
-names(dummy)[2] <- "D.mean.cm"
+names(dummy)[2] <- "D_mean_cm"
plots <- merge(plots, dummy)
# Summary statistics
-summary(plots[, c("G.m2.ha", "N.ha", "D.mean.cm")])
+summary(plots[, c("G_m2_ha", "N_ha", "D_mean_cm")])
# display histograms
par(mfrow = c(1, 3))
-hist(plots$G.m2.ha, main = "Basal area", xlab = "(m2/ha)", ylab = "Plot number")
-hist(plots$N.ha, main = "Stem density", xlab = "(/ha)", ylab = "Plot number")
-hist(plots$D.mean.cm, main = "Mean diameter", xlab = "(cm)", ylab = "Plot number")
+hist(plots$G_m2_ha, main = "Basal area", xlab = "(m2/ha)", ylab = "Plot number")
+hist(plots$N_ha, main = "Stem density", xlab = "(/ha)", ylab = "Plot number")
+hist(plots$D_mean_cm, main = "Mean diameter", xlab = "(cm)", ylab = "Plot number")
```
### Add stratum information from external data
@@ -545,25 +544,25 @@ public <- sf::st_read("../data/aba.model/GIS/Public4Montagnes.shp", stringsAsFac
# set coordinate system
sf::st_crs(public) <- 2154
# spatial query
-plots.sf <- sf::st_join(plots.sf, public)
+plots_sf <- sf::st_join(plots_sf, public)
# rename column and levels
-plots.sf$stratum <- factor(ifelse(is.na(plots.sf$FID), "private", "public"))
+plots_sf$stratum <- factor(ifelse(is.na(plots_sf$FID), "private", "public"))
# also add to data.frame
-plots <- merge(plots, sf::st_drop_geometry(plots.sf[, c("plotId", "stratum")]))
+plots <- merge(plots, sf::st_drop_geometry(plots_sf[, c("plot_id", "stratum")]))
# remove column
-plots.sf$FID <- NULL
+plots_sf$FID <- NULL
```
The following map displays the plots colored based on ownership, after spatial extraction of this information from the polygons of public forests (black borders).
```{r stratumPlot, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
# project points into map system
-plots.sf.transform <- sf::st_transform(plots.sf, 4326)
-public.transform <- sf::st_transform(public, 4326)
+plots_sf_transform <- sf::st_transform(plots_sf, 4326)
+public_transform <- sf::st_transform(public, 4326)
# display plots
ggmap::ggmap(map) +
- ggplot2::geom_sf(data = sf::st_geometry(public.transform), inherit.aes = FALSE, fill = NA, colour = "black") +
- ggplot2::geom_sf(data = plots.sf.transform, inherit.aes = FALSE, aes(color = stratum))
+ ggplot2::geom_sf(data = sf::st_geometry(public_transform), inherit.aes = FALSE, fill = NA, colour = "black") +
+ ggplot2::geom_sf(data = plots_sf_transform, inherit.aes = FALSE, aes(color = stratum))
```
## Save data before next tutorial
@@ -572,9 +571,9 @@ The following lines save the data required for the [area-based model calibration
```{r save.plots, include=TRUE}
save(plots, file = "../data/aba.model/output/plots.rda")
-save(llas.height,
- file = "../data/aba.model/output/llas.height.rda",
+save(llas_height,
+ file = "../data/aba.model/output/llas_height.rda",
compress = "bzip2"
)
-save(metrics.terrain, file = "../data/aba.model/output/metrics.terrain.rda")
-```
+save(metrics_terrain, file = "../data/aba.model/output/metrics_terrain.rda")
+```
\ No newline at end of file
diff --git a/R/area-based.2.model.calibration.Rmd b/R/area-based.2.model.calibration.Rmd
index 518d56f97084e436fecda1ebfd0d75dbb2e9188f..b63f1d7b347cfbc005c5ef1e620c35420dce1b9f 100644
--- a/R/area-based.2.model.calibration.Rmd
+++ b/R/area-based.2.model.calibration.Rmd
@@ -3,8 +3,8 @@ title: "R workflow for ABA prediction model calibration"
author: "Jean-Matthieu Monnet"
date: "`r Sys.Date()`"
output:
- pdf_document: default
html_document: default
+ pdf_document: default
papersize: a4
bibliography: "../bib/bibliography.bib"
---
@@ -33,13 +33,13 @@ The "Quatre Montagnes" dataset from France, prepared as described in the [data p
## Field data
-The file "plots.rda" contains the field data, organized as a data.frame named `plots`. For subsequent use in the workflow, the data.frame should contain at least two fields: `plotId` (unique plot identifier) and a forest stand parameter. Each line in the data.frame corresponds to a field plot. A factor variable is required to calibrate stratified models. Plot coordinates are required for subsequent inference computations.
+The file "plots.rda" contains the field data, organized as a data.frame named `plots`. For subsequent use in the workflow, the data.frame should contain at least two fields: `plot_id` (unique plot identifier) and a forest stand parameter. Each line in the data.frame corresponds to a field plot. A factor variable is required to calibrate stratified models. Plot coordinates are required for subsequent inference computations.
The provided data set includes one categorical variable: `stratum`, which corresponds to forest ownership, XY coordinates and three forest stand parameters :
-* basal area in m^2^/ha (`G.m2.ha`),
-+ stem density in /ha (`N.ha`),
-+ mean diameter at breast height in cm (`D.mean.cm`).
+* basal area in m^2^/ha (`G_m2_ha`),
++ stem density in /ha (`N_ha`),
++ mean diameter at breast height in cm (`D_mean_cm`).
Scatterplots of stand parameters are presented below, colored by ownership (green for public forest, blue otherwise).
@@ -48,7 +48,7 @@ Scatterplots of stand parameters are presented below, colored by ownership (gree
load(file = "../data/aba.model/output/plots.rda")
summary(plots)
# display forest variables
-plot(plots[, c("G.m2.ha", "N.ha", "D.mean.cm")],
+plot(plots[, c("G_m2_ha", "N_ha", "D_mean_cm")],
col = ifelse(plots$stratum == "public", "green", "blue")
)
```
@@ -59,24 +59,24 @@ Normalized ALS point clouds extracted over each plot, as well as terrain statist
```{r loadALSdData, include=TRUE, message = FALSE, warning = FALSE}
# list of LAS objects: normalized point clouds inside plot extent
-load("../data/aba.model/output/llas.height.rda")
+load("../data/aba.model/output/llas_height.rda")
# display one point cloud # lidR::plot(llasn[[1]])
-llas.height[[1]]
+llas_height[[1]]
```
The first lines of the terrain statistics are displayed hereafter.
```{r loadTerrainStats, echo = FALSE}
# terrain statistics previously computed with (non-normalized) ground points inside each plot extent
-load("../data/aba.model/output/metrics.terrain.rda")
-head(metrics.terrain[, 1:3], n = 3)
+load("../data/aba.model/output/metrics_terrain.rda")
+head(metrics_terrain[, 1:3], n = 3)
```
The following lines ensure that the plots are ordered in the same way in the three data objects.
```{r orderData, include=TRUE}
-llas.height <- llas.height[plots$plotId]
-metrics.terrain <- metrics.terrain[plots$plotId, ]
+llas_height <- llas_height[plots$plot_id]
+metrics_terrain <- metrics_terrain[plots$plot_id, ]
```
# ALS metrics computation
@@ -88,36 +88,36 @@ Two types of metrics can be computed.
## Point cloud metrics
-Point cloud metrics are computed with the function `lidaRtRee::cloudMetrics`, which applies the `lidR::cloud_metrics` to all point clouds in the 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::ABAmodelMetrics`.
+Point cloud metrics are computed with the function `lidaRtRee::clouds_metrics`, which applies the `lidR::cloud_metrics` to all point clouds in the 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::ABAmodelMetrics`.
```{r computeMetrics, include=TRUE}
# define function for later use
-aba.pointMetricsFUN <- ~ lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2)
+aba_point_metrics_fun <- ~ lidaRtRee::aba_metrics(Z, Intensity, ReturnNumber, Classification, 2)
# apply function on each point cloud in list
-metrics.points <- lidaRtRee::cloudMetrics(llas.height, aba.pointMetricsFUN)
-round(head(metrics.points[, 1:8], n = 3), 2)
+metrics_points <- lidaRtRee::clouds_metrics(llas_height, aba_point_metrics_fun)
+round(head(metrics_points[, 1:8], n = 3), 2)
```
## Tree metrics
-Tree metrics rely on a preliminary detection of trees, which is performed with the `lidaRtRee::treeSegmentation` function. For more details, please refer to the [tree detection tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/tree.detection.Rmd). Tree segmentation requires point clouds or canopy height models with an additional buffer in order to avoid border effects when computing tree characteristics. Once trees are detected, metrics are derived with the function `lidaRtRee::stdTreeMetrics`. A user-specific function can be specified to compute other metrics from the features of detected trees. Plot radius has to be specified as it is required to exclude trees detected outside of the plot, and to compute the plot surface. Tree segmentation is not relevant when the point cloud density is too low, typically below five points per m^2^. The function first computes a canopy height model which default resolution is 0.5 m, but this should be set to 1 m with low point densities.
+Tree metrics rely on a preliminary detection of trees, which is performed with the `lidaRtRee::tree_segmentation` function. For more details, please refer to the [tree detection tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/tree.detection.Rmd). Tree segmentation requires point clouds or canopy height models with an additional buffer in order to avoid border effects when computing tree characteristics. Once trees are detected, metrics are derived with the function `lidaRtRee::std_tree_metrics`. A user-specific function can be specified to compute other metrics from the features of detected trees. Plot radius has to be specified as it is required to exclude trees detected outside of the plot, and to compute the plot surface. Tree segmentation is not relevant when the point cloud density is too low, typically below five points per m^2^. The function first computes a canopy height model which default resolution is 0.5 m, but this should be set to 1 m with low point densities.
-```{r computeTreeMetrics, include=TRUE}
+```{r computeTreeMetrics, include=TRUE, warning = FALSE}
# resolution of canopy height model (m)
-aba.resCHM <- 0.5
+aba_res_chm <- 0.5
# specify plot radius to exclude trees located outside plots
-plot.radius <- 15
+plot_radius <- 15
# compute tree metrics
-metrics.tree <- lidaRtRee::cloudTreeMetrics(llas.height, plots[, c("X", "Y")],
- plot.radius,
- res = aba.resCHM,
+metrics_trees <- lidaRtRee::clouds_tree_metrics(llas_height, plots[, c("X", "Y")],
+ plot_radius,
+ res = aba_res_chm,
func = function(x) {
- lidaRtRee::stdTreeMetrics(x,
- area.ha = pi * plot.radius^2 / 10000
+ lidaRtRee::std_tree_metrics(x,
+ area_ha = pi * plot_radius^2 / 10000
)
}
)
-round(head(metrics.tree[, 1:5], n = 3), 2)
+round(head(metrics_trees[, 1:5], n = 3), 2)
```
## Other metrics
@@ -126,9 +126,9 @@ In case terrain metrics have been computed from the cloud of ground points only,
```{r bindMetrics, include=TRUE}
metrics <- cbind(
- metrics.points[plots$plotId, ],
- metrics.tree[plots$plotId, ],
- metrics.terrain[plots$plotId, 1:3]
+ metrics_points[plots$plot_id, ],
+ metrics_trees[plots$plot_id, ],
+ metrics_terrain[plots$plot_id, 1:3]
)
```
@@ -139,52 +139,52 @@ metrics <- cbind(
Once a dependent variable (forest parameter of interest) has been chosen, the function `lidaRtRee::ABAmodel` is used to select the linear regression model that yields the highest adjusted-R^2^ with a defined number of independent variables, while checking linear model assumptions. A Box-Cox transformation of the dependent variable can be applied to normalize its distribution, or a log transformation of all variables (parameter `transform`). Model details and cross-validation statistics are available from the returned object.
```{r modelCalibration, include=TRUE, message = FALSE, warning = FALSE}
-variable <- "G.m2.ha"
+variable <- "G_m2_ha"
# no subsample in this case
subsample <- 1:nrow(plots)
# model calibration
-model.ABA <- lidaRtRee::ABAmodel(plots[subsample, variable], metrics[subsample, ], transform = "boxcox", nmax = 4, xy = plots[subsample, c("X", "Y")])
+model_aba <- lidaRtRee::aba_build_model(plots[subsample, variable], metrics[subsample, ], transform = "boxcox", nmax = 4, xy = plots[subsample, c("X", "Y")])
# renames outputs with variable name
-row.names(model.ABA$stats) <- variable
+row.names(model_aba$stats) <- variable
# display selected linear regression model
-model.ABA$model
+model_aba$model
# display calibration and validation statistics
-model.ABA$stats
+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::ABAmodelPlot`. It is also informative to check the correlation of prediction errors with other forest or environmental variables.
+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.
In this example, only tree metrics are selected in the basal area prediction model. The model seems to fail to predict large values. The prediction errors are positively correlated with basal area because large values are under-estimated.
```{r modelPlot, include=TRUE, fig.height = 4.5, fig.width = 8}
# check correlation between errors and other variables
-round(cor(cbind(model.ABA$values$residual, plots[subsample, c("G.m2.ha", "N.ha", "D.mean.cm")], metrics.terrain[subsample, 1:3])), 2)[1, ]
+round(cor(cbind(model_aba$values$residual, plots[subsample, c("G_m2_ha", "N_ha", "D_mean_cm")], metrics_terrain[subsample, 1:3])), 2)[1, ]
# significance of correlation value
-cor.test(model.ABA$values$residual, plots[subsample, variable])
+cor.test(model_aba$values$residual, plots[subsample, variable])
# plot predicted VS field values
par(mfrow = c(1, 2))
-lidaRtRee::ABAmodelPlot(model.ABA, main = variable)
-plot(plots[subsample, c("G.m2.ha")], model.ABA$values$residual, ylab = "Prediction errors", xlab = "Field values")
+lidaRtRee::aba_plot(model_aba, main = variable)
+plot(plots[subsample, c("G_m2_ha")], model_aba$values$residual, ylab = "Prediction errors", xlab = "Field values")
abline(h = 0, lty = 2)
```
In case only point cloud metrics are used as potential inputs, the errors are hardly better distributed. Coloring points by ownership shows that plots located in private forests have the largest basal area values which tend to be under-estimated.
-```{r metrics.pointsOnly, include=TRUE, fig.height = 4.5, fig.width = 8}
-model.ABA.metrics.points <- lidaRtRee::ABAmodel(plots[subsample, variable], metrics.points[subsample, ], transform = "boxcox", nmax = 4, xy = plots[subsample, c("X", "Y")])
+```{r metrics_pointsOnly, include=TRUE, fig.height = 4.5, fig.width = 8}
+model_aba_metrics_points <- lidaRtRee::aba_build_model(plots[subsample, variable], metrics_points[subsample, ], transform = "boxcox", nmax = 4, xy = plots[subsample, c("X", "Y")])
# renames outputs
-row.names(model.ABA.metrics.points$stats) <- names(model.ABA.metrics.points$model) <- variable
-# model.ABA.metrics.points$model[[variable]]
-model.ABA.metrics.points$stats
-# cor.test(model.ABA.metrics.points$values$residual, plots[subsample, variable])
+row.names(model_aba_metrics_points$stats) <- names(model_aba_metrics_points$model) <- variable
+# model_aba_metrics_points$model[[variable]]
+model_aba_metrics_points$stats
+# cor.test(model_aba_metrics_points$values$residual, plots[subsample, variable])
par(mfrow = c(1, 2))
# plot predicted VS field values
-lidaRtRee::ABAmodelPlot(model.ABA.metrics.points,
+lidaRtRee::aba_plot(model_aba_metrics_points,
main = variable,
col = ifelse(plots$stratum == "public", "green", "blue")
)
legend("topleft", c("public", "private"), col = c("green", "blue"), pch = 1)
-plot(plots[subsample, c("G.m2.ha")],
- model.ABA.metrics.points$values$residual,
+plot(plots[subsample, c("G_m2_ha")],
+ model_aba_metrics_points$values$residual,
ylab = "Prediction errors", xlab = "Field values",
col = ifelse(plots$stratum == "public", "green", "blue")
)
@@ -196,13 +196,13 @@ abline(h = 0, lty = 2)
The following code calibrates models for several forest parameters. In case different transformations have to be performed on the parameters, models have to be calibrated one by one.
```{r multipleModels, include=TRUE, warning = FALSE, message = FALSE}
-models.ABA <- list()
-for (i in c("G.m2.ha", "D.mean.cm", "N.ha"))
+models_aba <- list()
+for (i in c("G_m2_ha", "D_mean_cm", "N_ha"))
{
- models.ABA[[i]] <- lidaRtRee::ABAmodel(plots[, i], metrics, transform = "boxcox", nmax = 4, xy = plots[, c("X", "Y")])
+ models_aba[[i]] <- lidaRtRee::aba_build_model(plots[, i], metrics, transform = "boxcox", nmax = 4, xy = plots[, c("X", "Y")])
}
# bind model stats in a data.frame
-model.stats <- do.call(rbind, lapply(models.ABA, function(x) {
+model_stats <- do.call(rbind, lapply(models_aba, function(x) {
x[["stats"]]
}))
```
@@ -217,20 +217,20 @@ The obtained models are presented below. The table columns correspond to:
```{r multipleModelsTable, echo = FALSE, fig.width = 12, fig.height = 4.5}
# prepare output for report
-table.output <- cbind(
- model.stats[, c("n", "formula")],
- round(model.stats[, c("adjR2", "looR2", "cvrmse")] * 100, 1),
- data.frame(rmse = round(model.stats[, "rmse"], 1))
+table_output <- cbind(
+ model_stats[, c("n", "formula")],
+ round(model_stats[, c("adjR2", "looR2", "cvrmse")] * 100, 1),
+ data.frame(rmse = round(model_stats[, "rmse"], 1))
)
-names(table.output) <- c("n", "metrics", "adj-R2.%", "CV-R2.%", "CV-RMSE.%", "CV-RMSE")
-knitr::kable(table.output)
+names(table_output) <- c("n", "metrics", "adj-R2.%", "CV-R2.%", "CV-RMSE.%", "CV-RMSE")
+knitr::kable(table_output)
#
par(mfrow = c(1, 3))
-for (i in names(models.ABA))
+for (i in names(models_aba))
{
- lidaRtRee::ABAmodelPlot(models.ABA[[i]], main = i)
+ lidaRtRee::aba_plot(models_aba[[i]], main = i)
}
-rm(models.ABA, model.stats)
+rm(models_aba, model_stats)
```
# Stratified models
@@ -241,7 +241,7 @@ When calibrating a statistical relationship between forest stand parameters, whi
## Calibration of stratum-specific models
-Stratum-specific models are computed and stored in a list during a `for` loop. The function `lidaRtRee::ABAmodelCombineStrata` then combines the list of models corresponding to each stratum to compute aggregated statistics for all plots, making it easier to compare stratified with non-stratified models.
+Stratum-specific models are computed and stored in a list during a `for` loop. The function `lidaRtRee::aba_combine_strata` then combines the list of models corresponding to each stratum to compute aggregated statistics for all plots, making it easier to compare stratified with non-stratified models.
In this example, the model for "private" yields a large error on the plot "Verc-C5-1", which considerably lowers the accuracy of the stratified approach.
@@ -249,37 +249,37 @@ In this example, the model for "private" yields a large error on the plot "Verc-
# stratification variable
strat <- "stratum"
# create list of models
-model.ABA.stratified <- list()
+model_aba_stratified <- list()
# calibrate each stratum model
for (i in levels(plots[, strat]))
{
subsample <- which(plots[, strat] == i)
if (length(subsample) > 0) {
- model.ABA.stratified[[i]] <- lidaRtRee::ABAmodel(plots[subsample, variable], metrics[subsample, ], transform = "boxcox", nmax = 4, xy = plots[subsample, c("X", "Y")])
+ model_aba_stratified[[i]] <- lidaRtRee::aba_build_model(plots[subsample, variable], metrics[subsample, ], transform = "boxcox", nmax = 4, xy = plots[subsample, c("X", "Y")])
}
}
# backup list of models for later use
-model.ABA.stratified.boxcox <- model.ABA.stratified
+model_aba_stratified_boxcox <- model_aba_stratified
# combine list of models into single object
-model.ABA.stratified <- lidaRtRee::ABAmodelCombineStrata(model.ABA.stratified, plots$plotId)
-# model.ABA.stratified$stats
+model_aba_stratified <- lidaRtRee::aba_combine_strata(model_aba_stratified, plots$plot_id)
+# model_aba_stratified$stats
```
```{r stratifiedmodelTable, echo=FALSE, fig.height = 4.5, fig.width = 8}
# bind model stats in a data.frame for comparison
-model.stats <- rbind(model.ABA$stats, model.ABA.stratified$stats)
-row.names(model.stats)[1] <- "NOT.STRATIFIED"
+model_stats <- rbind(model_aba$stats, model_aba_stratified$stats)
+row.names(model_stats)[1] <- "NOT.STRATIFIED"
# prepare output for report
-table.output <- cbind(
- model.stats[, c("n", "formula")],
- round(model.stats[, c("adjR2", "looR2", "cvrmse")] * 100, 1),
- data.frame(rmse = round(model.stats[, "rmse"], 1))
+table_output <- cbind(
+ model_stats[, c("n", "formula")],
+ round(model_stats[, c("adjR2", "looR2", "cvrmse")] * 100, 1),
+ data.frame(rmse = round(model_stats[, "rmse"], 1))
)
-names(table.output) <- c("n", "metrics", "adj-R2.%", "CV-R2.%", "CV-RMSE.%", "CV-RMSE")
-knitr::kable(table.output)
+names(table_output) <- c("n", "metrics", "adj-R2.%", "CV-R2.%", "CV-RMSE.%", "CV-RMSE")
+knitr::kable(table_output)
par(mfrow = c(1, 2))
-lidaRtRee::ABAmodelPlot(model.ABA, main = paste0(variable, ", not stratified"))
-lidaRtRee::ABAmodelPlot(model.ABA.stratified, main = paste0(variable, ", stratified"))
+lidaRtRee::aba_plot(model_aba, main = paste0(variable, ", not stratified"))
+lidaRtRee::aba_plot(model_aba_stratified, main = paste0(variable, ", stratified"))
```
## Stratified models with stratum-specific variable tranformations
@@ -291,35 +291,35 @@ In case one wants to apply different variable transformations, or use different
```{r stratifiedmodelCalibrationTransformation, include=TRUE, warning = FALSE}
# create list of models for no transformation
-model.ABA.stratified.none <- list()
+model_aba_stratified.none <- list()
# calibrate each stratum model
for (i in levels(plots[, strat]))
{
subsample <- which(plots[, strat] == i)
if (length(subsample) > 0) {
- model.ABA.stratified.none[[i]] <- lidaRtRee::ABAmodel(plots[subsample, variable], metrics.points[subsample, ], transform = "none", xy = plots[subsample, c("X", "Y")])
+ model_aba_stratified.none[[i]] <- lidaRtRee::aba_build_model(plots[subsample, variable], metrics_points[subsample, ], transform = "none", xy = plots[subsample, c("X", "Y")])
}
}
# combine list of models into single object
-model.ABA.stratified.mixed <- lidaRtRee::ABAmodelCombineStrata(list(private = model.ABA.stratified.none[["private"]], public = model.ABA.stratified.boxcox[["public"]]), plots$plotId)
+model_aba_stratified_mixed <- lidaRtRee::aba_combine_strata(list(private = model_aba_stratified.none[["private"]], public = model_aba_stratified_boxcox[["public"]]), plots$plot_id)
# bind model stats in a data.frame for comparison
-model.stats <- rbind(model.ABA$stats, model.ABA.stratified.mixed$stats)
-row.names(model.stats)[1] <- "NOT.STRATIFIED"
+model_stats <- rbind(model_aba$stats, model_aba_stratified_mixed$stats)
+row.names(model_stats)[1] <- "NOT.STRATIFIED"
```
```{r stratifiedmodelCalibrationTransformationTable, echo = FALSE, fig.height = 4.5, fig.width = 8}
# prepare output for report
-table.output <- cbind(
- model.stats[, c("n", "formula", "transform")],
- round(model.stats[, c("adjR2", "looR2", "cvrmse")] * 100, 1),
- data.frame(rmse = round(model.stats[, "rmse"], 1))
+table_output <- cbind(
+ model_stats[, c("n", "formula", "transform")],
+ round(model_stats[, c("adjR2", "looR2", "cvrmse")] * 100, 1),
+ data.frame(rmse = round(model_stats[, "rmse"], 1))
)
-names(table.output) <- c("n", "metrics", "transform", "adj-R2.%", "CV-R2.%", "CV-RMSE.%", "CV-RMSE")
-knitr::kable(table.output)
+names(table_output) <- c("n", "metrics", "transform", "adj-R2.%", "CV-R2.%", "CV-RMSE.%", "CV-RMSE")
+knitr::kable(table_output)
# graphics
par(mfrow = c(1, 2))
-lidaRtRee::ABAmodelPlot(model.ABA, main = paste0(variable, ", not stratified"))
-lidaRtRee::ABAmodelPlot(model.ABA.stratified.mixed, main = paste0(variable, ", stratified"))
+lidaRtRee::aba_plot(model_aba, main = paste0(variable, ", not stratified"))
+lidaRtRee::aba_plot(model_aba_stratified_mixed, main = paste0(variable, ", stratified"))
```
# Save data before next tutorial
@@ -327,7 +327,10 @@ lidaRtRee::ABAmodelPlot(model.ABA.stratified.mixed, main = paste0(variable, ", s
The following lines save the data required for the [area-based mapping step](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.3.mapping.Rmd).
```{r saveModels, eval=FALSE}
-save(model.ABA.stratified.mixed, model.ABA, aba.pointMetricsFUN, aba.resCHM,
+save(model_aba_stratified_mixed, model_aba, aba_point_metrics_fun, aba_res_chm,
file = "../data/aba.model/output/models.rda"
)
-```
+# save data for lidaRtRee package
+# quatre_montagnes <- cbind(plots, metrics)
+# save(quatre_montagnes, file = "quatre_montagnes.rda")
+```
\ No newline at end of file
diff --git a/R/area-based.3.mapping.and.inference.Rmd b/R/area-based.3.mapping.and.inference.Rmd
index bdeb0d0f222db87065878e186581604dd1815297..2c586c7d122767c8bc8c39cb6a1c5903626fdf07 100644
--- a/R/area-based.3.mapping.and.inference.Rmd
+++ b/R/area-based.3.mapping.and.inference.Rmd
@@ -3,8 +3,8 @@ title: "R workflow for forest mapping and inference from ABA prediction models"
author: "Jean-Matthieu Monnet"
date: "`r Sys.Date()`"
output:
- pdf_document: default
html_document: default
+ pdf_document: default
papersize: a4
bibliography: "../bib/bibliography.bib"
---
@@ -67,12 +67,12 @@ unlink(temp)
```{r loadALS, include = TRUE, warning = FALSE}
# build catalogs
-cata.height <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.norm.laz/")
-cata.altitude <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.laz/")
-lidR::opt_progress(cata.altitude) <- lidR::opt_progress(cata.height) <- FALSE
+cata_height <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.norm.laz/")
+cata_altitude <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.laz/")
+lidR::opt_progress(cata_altitude) <- lidR::opt_progress(cata_height) <- FALSE
# set projection info
-lidR::projection(cata.height) <- lidR::projection(cata.altitude) <- 2154
-cata.height
+lidR::projection(cata_height) <- lidR::projection(cata_altitude) <- 2154
+cata_height
```
## GIS data
@@ -103,24 +103,24 @@ The following paragraph show how to compute metrics for a 300 x 300 m^2^ subsamp
```{r loadTile, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 6}
# load ALS tile
-point.cloud.height <- lidR::clip_rectangle(cata.height, 899600, 6447600, 899900, 6447900)
-point.cloud.ground <- lidR::filter_ground(lidR::clip_rectangle(cata.altitude, 899600, 6447600, 899900, 6447900))
+point_cloud_height <- lidR::clip_rectangle(cata_height, 899600, 6447600, 899900, 6447900)
+point_cloud_ground <- lidR::filter_ground(lidR::clip_rectangle(cata_altitude, 899600, 6447600, 899900, 6447900))
```
### Point cloud metrics
-For computation of point cloud metrics, the same function used in `lidaRtRee::cloudMetrics` 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::grid_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.pointMetricsFUN, res = resolution)
+metrics_points <- lidR::grid_metrics(point_cloud_height, aba_point_metrics_fun, res = resolution)
```
```{r plotpointMetrics, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 9, fig.height = 3}
par(mfrow = c(1, 3))
for (i in c("zmean", "zsd", "imean"))
{
- raster::plot(metrics.points[[i]], main = i)
+ raster::plot(metrics_points[[i]], main = i)
}
```
@@ -136,41 +136,41 @@ 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.resCHM)
+chm <- lidaRtRee::points2DSM(point_cloud_height, res = aba_res_chm)
# replace NA, low and high values
chm[is.na(chm) | chm < 0 | chm > 60] <- 0
-# tree top detection (same parameters as used by cloudTreeMetrics in calibration step -> default parameters)
-segms <- lidaRtRee::treeSegmentation(chm)
+# tree top detection (same parameters as used by clouds_tree_metrics in calibration step -> default parameters)
+segms <- lidaRtRee::tree_segmentation(chm)
# extraction to data.frame
-trees <- lidaRtRee::treeExtraction(segms$filled.dem, segms$local.maxima, segms$segments.id)
+trees <- lidaRtRee::tree_extraction(segms$filled_dem, segms$local_maxima, segms$segments_id)
# compute raster metrics
-metrics.trees <- lidaRtRee::rasterMetrics(trees[, -1],
+metrics_trees <- lidaRtRee::raster_metrics(trees[, -1],
res = resolution,
fun = function(x) {
- lidaRtRee::stdTreeMetrics(x, resolution^2 / 10000)
+ lidaRtRee::std_tree_metrics(x, resolution^2 / 10000)
}, output = "raster"
)
# remove NAs and NaNs
-metrics.trees[!is.finite(metrics.trees)] <- 0
+metrics_trees[!is.finite(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.treechm <- segms$filled.dem
+r_treechm <- segms$filled_dem
# set chm to 0 in non segment area
-r.treechm[segms$segments.id == 0] <- NA
+r_treechm[segms$segments_id == 0] <- NA
# compute raster metrics
-dummy <- lidaRtRee::rasterMetrics(r.treechm,
+dummy <- lidaRtRee::raster_metrics(r_treechm,
res = resolution,
fun = function(x) {
c(
- sum(!is.na(x$filled.dem) / (resolution / aba.resCHM)^2),
- mean(x$filled.dem, na.rm = T)
+ sum(!is.na(x$filled_dem) / (resolution / aba_res_chm)^2),
+ mean(x$filled_dem, na.rm = T)
)
},
output = "raster"
)
-names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot")
-metrics.trees <- raster::addLayer(metrics.trees, dummy)
+names(dummy) <- c("TreeCanop_cover_in_plot", "TreeCanopy_meanH_in_plot")
+metrics_trees <- raster::addLayer(metrics_trees, dummy)
```
The detected trees are plotted below, with the CHM as background.
@@ -185,29 +185,29 @@ Some metrics maps are displayed below.
```{r plottreeMetrics, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 9, fig.height = 3}
par(mfrow = c(1, 3))
-for (i in c("Tree.meanH", "Tree.sdH", "Tree.density"))
+for (i in c("Tree_meanH", "Tree_sdH", "Tree_density"))
{
- raster::plot(metrics.trees[[i]], main = i)
+ raster::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::points2terrainStats` 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::grid_metrics`.
```{r terrainMetrics, include=TRUE, message=FALSE, warning=FALSE}
# compute terrain metrics
f <- function(x, y, z) {
- as.list(lidaRtRee::points2terrainStats(data.frame(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::grid_metrics(point_cloud_ground, ~ f(X, Y, Z), res = resolution)
```
```{r plotterrainMetrics, include=TRUE, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 9, fig.height = 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)")
+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)")
```
### Batch processing of tiles
@@ -225,26 +225,26 @@ A buffering procedure is designed to handle the border effects when detecting tr
```{r paramCluster, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
# processing parameters
-b.size <- 10
+b_size <- 10
# height threshold (m) for high points removal (points above this threshold are considered
# as outliers)
-h.points <- 60
+h_points <- 60
# points classes to retain for analysis (vegetation + ground)
# ground should be 2
-class.points <- c(0, 1, 2, 3, 4, 5)
+class_points <- c(0, 1, 2, 3, 4, 5)
```
This first chunk of code computes tree and point metrics from the normalized ALS tiles.
```{r batchProcessing, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
# processing by tile
-metrics <- foreach::foreach(i = 1:nrow(cata.height), .errorhandling = "remove") %dopar% {
+metrics <- foreach::foreach(i = 1:nrow(cata_height), .errorhandling = "remove") %dopar% {
# tile extent
- b.box <- as.numeric(cata.height@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
+ b_box <- as.numeric(cata_height@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# load tile extent plus buffer
a <- try(lidR::clip_rectangle(
- cata.height, b.box[1] - b.size, b.box[2] - b.size, b.box[3] + b.size,
- b.box[4] + b.size
+ cata_height, b_box[1] - b_size, b_box[2] - b_size, b_box[3] + b_size,
+ b_box[4] + b_size
))
#
# check if points are successfully loaded
@@ -254,11 +254,11 @@ metrics <- foreach::foreach(i = 1:nrow(cata.height), .errorhandling = "remove")
# add 'buffer' flag to points in buffer with TRUE value in this new attribute
a <- lidR::add_attribute(
a,
- a$X < b.box[1] | a$Y < b.box[2] | a$X >= b.box[3] | a$Y >= b.box[4],
+ a$X < b_box[1] | a$Y < b_box[2] | a$X >= b_box[3] | a$Y >= b_box[4],
"buffer"
)
# remove unwanted point classes, and points higher than height threshold
- a <- lidR::filter_poi(a, is.element(Classification, class.points) & Z <= h.points)
+ a <- lidR::filter_poi(a, is.element(Classification, class_points) & Z <= h_points)
# check number of remaining points
if (length(a) == 0) {
return(NULL)
@@ -267,43 +267,43 @@ metrics <- foreach::foreach(i = 1:nrow(cata.height), .errorhandling = "remove")
a@data$Z[a@data$Z < 0] <- 0
#
# compute chm
- chm <- lidaRtRee::points2DSM(a, res = aba.resCHM)
+ chm <- lidaRtRee::points2DSM(a, res = aba_res_chm)
# replace NA, low and high values by 0
- chm[is.na(chm) | chm < 0 | chm > h.points] <- 0
+ chm[is.na(chm) | chm < 0 | chm > h_points] <- 0
#
# compute tree metrics
# tree top detection (default parameters)
- segms <- lidaRtRee::treeSegmentation(chm)
+ segms <- lidaRtRee::tree_segmentation(chm)
# extraction to data.frame
- trees <- lidaRtRee::treeExtraction(segms$filled.dem, segms$local.maxima, segms$segments.id)
+ trees <- lidaRtRee::tree_extraction(segms$filled_dem, segms$local_maxima, segms$segments_id)
# remove trees outside of tile
- trees <- trees[trees$x >= b.box[1] & trees$x < b.box[3] & trees$y >= b.box[2] & trees$y < b.box[4], ]
+ trees <- trees[trees$x >= b_box[1] & trees$x < b_box[3] & trees$y >= b_box[2] & trees$y < b_box[4], ]
# compute raster metrics
- metrics.trees <- lidaRtRee::rasterMetrics(trees[, -1],
+ metrics_trees <- lidaRtRee::raster_metrics(trees[, -1],
res = resolution,
fun = function(x) {
- lidaRtRee::stdTreeMetrics(x, resolution^2 / 10000)
+ lidaRtRee::std_tree_metrics(x, resolution^2 / 10000)
},
output = "raster"
)
# compute canopy cover in trees and canopy mean height in trees
# in region of interest, because it is not in previous step.
- r.treechm <- segms$filled.dem
+ r_treechm <- segms$filled_dem
# set chm to NA in non segment area
- r.treechm[segms$segments.id == 0] <- NA
+ r_treechm[segms$segments_id == 0] <- NA
# compute raster metrics
- metrics.trees2 <- lidaRtRee::rasterMetrics(
- raster::crop(r.treechm, raster::extent(b.box[1], b.box[3], b.box[2], b.box[4])),
+ metrics_trees2 <- lidaRtRee::raster_metrics(
+ raster::crop(r_treechm, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4])),
res = resolution,
fun = function(x) {
c(
- sum(!is.na(x$filled.dem)) / (resolution / aba.resCHM)^2,
- mean(x$filled.dem, na.rm = T)
+ sum(!is.na(x$filled_dem)) / (resolution / aba_res_chm)^2,
+ mean(x$filled_dem, na.rm = T)
)
},
output = "raster"
)
- names(metrics.trees2) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot")
+ names(metrics_trees2) <- c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
#
# compute 1D height metrics
# remove buffer points
@@ -313,23 +313,23 @@ metrics <- foreach::foreach(i = 1:nrow(cata.height), .errorhandling = "remove")
return(NULL)
}
# all points metrics
- metrics.points <- lidR::grid_metrics(a, aba.pointMetricsFUN, res = resolution)
+ metrics_points <- lidR::grid_metrics(a, aba_point_metrics_fun, res = resolution)
#
- # extend / crop to match metrics.points
- metrics.trees <- raster::extend(metrics.trees, metrics.points, values = NA)
- metrics.trees2 <- raster::extend(metrics.trees2, metrics.points, values = NA)
- metrics.trees <- raster::crop(metrics.trees, metrics.points)
- metrics.trees2 <- raster::crop(metrics.trees2, metrics.points)
+ # extend / crop to match metrics_points
+ metrics_trees <- raster::extend(metrics_trees, metrics_points, values = NA)
+ metrics_trees2 <- raster::extend(metrics_trees2, metrics_points, values = NA)
+ metrics_trees <- raster::crop(metrics_trees, metrics_points)
+ metrics_trees2 <- raster::crop(metrics_trees2, metrics_points)
# merge rasterstacks
- metrics <- metrics.points
- metrics <- raster::addLayer(metrics, metrics.trees)
- metrics <- raster::addLayer(metrics, metrics.trees2)
+ metrics <- metrics_points
+ metrics <- raster::addLayer(metrics, metrics_trees)
+ metrics <- raster::addLayer(metrics, metrics_trees2)
return(metrics)
}
# mosaic rasters
-names.metrics <- names(metrics[[1]])
-metrics.map <- do.call(raster::merge, metrics)
-names(metrics.map) <- names.metrics
+names_metrics <- names(metrics[[1]])
+metrics_map <- do.call(raster::merge, metrics)
+names(metrics_map) <- names_metrics
```
@@ -338,16 +338,16 @@ The second chunk of code computes terrain metrics from the ALS tiles with altitu
```{r batchProcessingTerrain, include=TRUE, message=FALSE, warning=FALSE}
# function to compute terrain metrics to input to grid_metrics
f <- function(x, y, z) {
- as.list(lidaRtRee::points2terrainStats(data.frame(x, y, z)))
+ as.list(lidaRtRee::terrain_points_metrics(data.frame(x, y, z)))
}
#
-metrics.terrain <- foreach::foreach(i = 1:nrow(cata.altitude), .errorhandling = "remove") %dopar% {
+metrics_terrain <- foreach::foreach(i = 1:nrow(cata_altitude), .errorhandling = "remove") %dopar% {
# tile extent
- b.box <- as.numeric(cata.altitude@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
+ b_box <- as.numeric(cata_altitude@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# load tile extent plus buffer
a <- try(lidR::clip_rectangle(
- cata.altitude, b.box[1], b.box[2], b.box[3],
- b.box[4]
+ cata_altitude, b_box[1], b_box[2], b_box[3],
+ b_box[4]
))
# check if points are successfully loaded
if (class(a) == "try-error") {
@@ -358,39 +358,39 @@ metrics.terrain <- foreach::foreach(i = 1:nrow(cata.altitude), .errorhandling =
lidR::grid_metrics(a, ~ f(X, Y, Z), res = resolution)
}
# mosaic rasters
-names.metrics <- names(metrics.terrain[[1]])
-metrics.terrain <- do.call(raster::merge, metrics.terrain)
-names(metrics.terrain) <- names.metrics
+names_metrics <- names(metrics_terrain[[1]])
+metrics_terrain <- do.call(raster::merge, metrics_terrain)
+names(metrics_terrain) <- names_metrics
```
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")
+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_map <- raster::addLayer(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.
```{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)
-save(metrics.map, file = "../data/aba.model/output/metrics.map.rda")
+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)
+save(metrics_map, file = "../data/aba.model/output/metrics_map.rda")
```
## Forest parameters mapping
### Single (not stratified) model
-Forest parameters are then mapped using the mapped metrics and the previously calibrated model with the function `lidaRtRee::ABApredict`.
+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 = 5.1, fig.height = 3}
-prediction.map <- lidaRtRee::ABApredict(model.ABA, metrics.map)
-raster::plot(prediction.map, main = "prediction map")
+prediction_map <- lidaRtRee::aba_predict(model_aba, metrics_map)
+raster::plot(prediction_map, main = "prediction map")
```
### Stratified models
@@ -399,48 +399,48 @@ 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 <- raster::rasterize(public, metrics_map, field = 1)
# set to 0 in non-public forests areas
-metrics.map$stratum[is.na(metrics.map$stratum)] <- 0
+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, raster::extent(metrics_map))
```
-Then the function `lidaRtRee::ABApredict` 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 strata.
```{r mapForestStratified, include=TRUE, message=FALSE, warning=FALSE}
# map with stratified model
-prediction.map.mixed <- lidaRtRee::ABApredict(model.ABA.stratified.mixed, metrics.map, stratum = "stratum")
+prediction_map_mixed <- lidaRtRee::aba_predict(model_aba_stratified_mixed, metrics_map, stratum = "stratum")
# for checking purposes
# extracted "private" stratum model from combined model
-model.private <- list(model = model.ABA.stratified.mixed$model$private, stats = model.ABA.stratified.mixed$stats["private", ])
+model_private <- list(model = model_aba_stratified_mixed$model$private, stats = model_aba_stratified_mixed$stats["private", ])
# produce corresponding map
-prediction.map.private <- lidaRtRee::ABApredict(model.private, metrics.map)
+prediction_map_private <- lidaRtRee::aba_predict(model_private, metrics_map)
# extracted "public" stratum model from combined model
-model.public <- list(model = model.ABA.stratified.mixed$model$public, stats = model.ABA.stratified.mixed$stats["public", ])
+model_public <- list(model = model_aba_stratified_mixed$model$public, stats = model_aba_stratified_mixed$stats["public", ])
# produce corresponding map
-prediction.map.public <- lidaRtRee::ABApredict(model.public, metrics.map)
+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 = 5}
# plot all maps
par(mfrow = c(2, 2))
# plot raster and vector
-raster::plot(metrics.map$stratum, main = "Ownership", legend = FALSE, col = c("green", "blue"), xaxt = "n", yaxt = "n")
+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)))
-plot(public.cropped, col = NA, add = TRUE)
-raster::plot(prediction.map.mixed, main = "Stratified model", zlim = 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)
+limits <- range(c(raster::values(prediction_map_private), raster::values(prediction_map_public)))
+plot(public_cropped, col = NA, add = TRUE)
+raster::plot(prediction_map_mixed, main = "Stratified model", zlim = 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
-To avoid applying models in non-forest areas and to remove the extremes values that may have been predicted due to outliers in the ALS and metrics values, the function `lidaRtRee::cleanRaster` can be applied:
+To avoid applying models in non-forest areas and to remove the extremes values that may have been predicted due to outliers in the ALS and metrics values, the function `lidaRtRee::clean_raster` can be applied:
* it applies a lower and upper threshold to map values,
* it sets to 0 the NA values in the map,
@@ -450,25 +450,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 <- raster::rasterize(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)
+raster::values(forest_mask) <- ifelse(raster::values(forest_mask) == 27, 1, NA)
# crop forest vector
-forest.cropped <- sf::st_crop(forest, raster::extent(metrics.map))
+forest_cropped <- sf::st_crop(forest, raster::extent(metrics_map))
# clean raster
# threshold values and set to NA outside of forest
-prediction.map.mixed.clean <- lidaRtRee::cleanRaster(prediction.map.mixed, c(0, 80), forest.mask)
+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.5}
par(mfrow = c(1, 3))
-raster::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")
-plot(sf::st_geometry(forest.cropped), add = TRUE)
+raster::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")
+plot(sf::st_geometry(forest_cropped), add = TRUE)
```
# Inference
-Work in progress...
+Work in progress...
\ No newline at end of file
diff --git a/R/coregistration.Rmd b/R/coregistration.Rmd
index d82ccda70de8b697c8bc35201803f62cd1ec77d1..305760f79cf42892cebc2b777cc68e52b0e5115a 100644
--- a/R/coregistration.Rmd
+++ b/R/coregistration.Rmd
@@ -3,8 +3,8 @@ title: "R workflow for field plot coregistration with ALS data"
author: "Jean-Matthieu Monnet"
date: "`r Sys.Date()`"
output:
- pdf_document: default
html_document: default
+ pdf_document: default
papersize: a4
bibliography: "../bib/bibliography.bib"
---
@@ -90,7 +90,7 @@ cata <- lidR::readALSLAScatalog("/directory_with_classified_laz_files/")
# set coordinate system
lidR::projection(cata) <- 2154
# extract LAS data on plot extent
-las <- lidR::clip_circle(cata, p$XGPS, p$YGPS, p.radius + b.size + 5)
+las <- lidR::clip_circle(cata, p$XGPS, p$YGPS, p_radius + b_size + 5)
# normalize heights if point cloud are not already normalized
las <- lapply(las, function(x) {
lidR::normalize_height(x, lidR::tin())
@@ -107,13 +107,13 @@ Several parameters have to be specified for optimal results.
# vegetation height threshold for removal of high values
hmax <- 50
# field plot radius for computation of pseudo-CHM on the inventory area
-p.radius <- 14.105
+p_radius <- 14.105
# raster resolution for image processing
-r.res <- 0.5
+r_res <- 0.5
# maximum distance around initial position to search for coregistration
-b.size <- 5
+b_size <- 5
# increment step to search for coregistration (should be a multiple of resolution)
-s.step <- 0.5
+s_step <- 0.5
```
## Coregistration of one plot
@@ -126,13 +126,13 @@ The first step is to compute the canopy height model from the ALS data, and remo
# Choose a plot as example
i <- 10
# compute canopy height model
-chm <- lidaRtRee::points2DSM(las[[i]], res = r.res)
+chm <- lidaRtRee::points2DSM(las[[i]], res = r_res)
# apply threshold to canopy height model (CHM)
chm[chm > hmax] <- hmax
# fill NA values with 0
chm[is.na(chm)] <- 0
# apply median filter (3x3 window) to CHM
-chmfilt <- lidaRtRee::demFiltering(chm, "Median", 3, sigmap = 0)$non.linear.image
+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")
@@ -149,23 +149,23 @@ centre <- p[i, c("XGPS", "YGPS")]
# extract plot trees
trees <- ap[ap$plac == p$placette[i], ]
# create raster with plot extent
-r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p.radius, resolution = r.res)
+r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p_radius, resolution = r_res)
# keep only trees with diameter information
trees <- trees[!is.na(trees[, "dia"]), ]
# create plot mask
-r.mask <- lidaRtRee::rasterXYMask(rbind(
+r_mask <- lidaRtRee::raster_xy_mask(rbind(
c(centre$XGPS, centre$YGPS),
c(centre$XGPS, centre$YGPS)
),
-c(p.radius, p.radius), r,
+c(p_radius, p_radius), r,
binary = T
)
# replace 0 by NA
-r.mask[r.mask == 0] <- NA
+r_mask[r_mask == 0] <- NA
# specify projection
-raster::crs(r.mask) <- 2154
+raster::crs(r_mask) <- 2154
# display plot mask
-raster::plot(r.mask, main = "Plot mask and tree positions", legend = FALSE)
+raster::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
@@ -180,13 +180,13 @@ First the function computes the correlation for different possible translations
```{r correlation, include = TRUE}
# compute correlation image
coreg <- lidaRtRee::coregistration(chmfilt, trees[, c("x", "y", "dia")],
- mask = r.mask,
- buffer = b.size, step = s.step, dm = 2, plot = F
+ mask = r_mask,
+ buffer = b_size, step = s_step, dm = 2, plot = FALSE
)
# correlation image statistics
-round(coreg$local.max, 2)
+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`.
+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}
par(mfrow = c(1, 3))
@@ -194,25 +194,25 @@ raster::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,
+raster::plot(coreg$correlation_raster,
main = "Correlation image",
col = cm.colors(16)
)
# plot local maximum
-graphics::points(coreg$local.max$dx1, coreg$local.max$dy1, pch = 4)
+graphics::points(coreg$local_max$dx1, coreg$local_max$dy1, pch = 4)
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
legend("topleft", "Maximum", pch = 4)
#
raster::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,
+graphics::points(trees$x + coreg$local_max$dx1, trees$y + coreg$local_max$dy1,
cex = trees$dia / 40, col = "red"
)
# display initial plot center
graphics::points(p[i, c("XGPS", "YGPS")], pch = 3)
# display coregistered plot center
-graphics::points(p$XGPS[i] + coreg$local.max$dx1, p$YGPS[i] + coreg$local.max$dy1,
+graphics::points(p$XGPS[i] + coreg$local_max$dx1, p$YGPS[i] + coreg$local_max$dy1,
pch = 3,
col = "red"
)
@@ -237,15 +237,15 @@ library(foreach)
First CHMs are calculated for each point cloud contained in the list.
```{r batch.chm, include = TRUE, eval=TRUE}
-l.chm <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% {
+l_chm <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% {
# compute CHM
- chm <- lidaRtRee::points2DSM(las[[i]], res = r.res)
+ chm <- lidaRtRee::points2DSM(las[[i]], res = r_res)
# apply threshold to canopy height model (CHM)
chm[chm > hmax] <- hmax
# fill NA values with 0
chm[is.na(chm)] <- 0
# apply median filter (3x3 window) to CHM
- chmfilt <- lidaRtRee::demFiltering(chm, "Median", 3, sigmap = 0)[[1]]
+ chmfilt <- lidaRtRee::dem_filtering(chm, "Median", 3, sigmap = 0)[[1]]
return(chmfilt)
}
```
@@ -253,28 +253,28 @@ l.chm <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar%
### Trees lists extraction and plot masks computation
```{r batch.mask, include = TRUE, eval=TRUE}
-l.field <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% {
+l_field <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% {
# plot centre
centre <- p[i, c("XGPS", "YGPS")]
# extract plot trees
trees <- ap[ap$plac == p$placette[i], ]
# create raster with plot extent
- r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p.radius, resolution = r.res)
+ r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p_radius, resolution = r_res)
# keep only trees with diameter information
trees <- trees[!is.na(trees[, "dia"]), ]
# create plot mask
- r.mask <- lidaRtRee::rasterXYMask(rbind(
+ r_mask <- lidaRtRee::raster_xy_mask(rbind(
c(centre$XGPS, centre$YGPS),
c(centre$XGPS, centre$YGPS)
),
- c(p.radius, p.radius), r,
+ c(p_radius, p_radius), r,
binary = T
)
# replace 0 by NA
- r.mask[r.mask == 0] <- 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))
+ raster::crs(r_mask) <- 2154
+ return(list(trees = trees[, c("x", "y", "dia")], r_mask = r_mask))
}
```
@@ -282,10 +282,10 @@ l.field <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar
Then the correlation image and corresponding statistics are computed for each plot.
```{r batch.correlation, include = TRUE, eval=TRUE}
-l.coreg <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% {
- lidaRtRee::coregistration(l.chm[[i]], l.field[[i]]$trees,
- mask = l.field[[i]]$r.mask,
- buffer = b.size, step = s.step, dm = 2, plot = F
+l_coreg <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% {
+ lidaRtRee::coregistration(l_chm[[i]], l_field[[i]]$trees,
+ mask = l_field[[i]]$r_mask,
+ buffer = b_size, step = s_step, dm = 2, plot = FALSE
)
}
```
@@ -297,8 +297,8 @@ translations <- foreach::foreach(
) %dopar% {
# create data.frame with coregistration results and new plot coordinates
data.frame(
- plotid = p$placette[i], X.cor = p$XGPS[i] + l.coreg[[i]]$local.max$dx1,
- Y.cor = p$YGPS[i] + l.coreg[[i]]$local.max$dy1, l.coreg[[i]]$local.max
+ plotid = p$placette[i], X_cor = p$XGPS[i] + l_coreg[[i]]$local_max$dx1,
+ Y_cor = p$YGPS[i] + l_coreg[[i]]$local_max$dy1, l_coreg[[i]]$local_max
)
}
#
@@ -312,21 +312,21 @@ for (i in 1:length(las))
{
# CHM
pdf(file = paste("Coregistration_", p$placette[i], ".pdf", sep = ""))
- raster::plot(l.chm[[i]], asp = 1)
+ raster::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
+ graphics::points(l_field[[i]]$trees$x, l_field[[i]]$trees$y,
+ cex = l_field[[i]]$trees$dia / 40
)
# display coregistered tree positions
- graphics::points(l.field[[i]]$trees$x + l.coreg[[i]]$local.max$dx1,
- l.field[[i]]$trees$y + l.coreg[[i]]$local.max$dy1,
- cex = l.field[[i]]$trees$dia / 40, col = "red"
+ graphics::points(l_field[[i]]$trees$x + l_coreg[[i]]$local_max$dx1,
+ l_field[[i]]$trees$y + l_coreg[[i]]$local_max$dy1,
+ cex = l_field[[i]]$trees$dia / 40, col = "red"
)
# display initial plot center
graphics::points(p[i, c("XGPS", "YGPS")], pch = 3)
# display coregistered plot center
- graphics::points(p$XGPS[i] + l.coreg[[i]]$local.max$dx1,
- p$YGPS[i] + l.coreg[[i]]$local.max$dy1,
+ graphics::points(p$XGPS[i] + l_coreg[[i]]$local_max$dx1,
+ p$YGPS[i] + l_coreg[[i]]$local_max$dy1,
pch = 3, col = "red"
)
graphics::legend("topleft", c("Initial", "Coregistered"),
@@ -342,7 +342,7 @@ for (i in 1:length(las))
Optimized plot center positions are added to the original data.
```{r batch.merge, include = TRUE, eval=TRUE}
-p <- base::merge(p, translations[, c("plotid", "X.cor", "Y.cor")],
+p <- base::merge(p, translations[, c("plotid", "X_cor", "Y_cor")],
by.x = "placette",
by.y = "plotid"
)
@@ -351,15 +351,15 @@ round(head(p, n = 3L), 2)
Tree positions are corrected to account for new center position.
```{r batch.correct, include = TRUE, eval=TRUE}
# add plot center coregistered coordinates to trees data.frame
-ap <- base::merge(ap, p[, c("placette", "X.cor", "Y.cor")], by.x = "plac", by.y = "placette")
+ap <- base::merge(ap, p[, c("placette", "X_cor", "Y_cor")], by.x = "plac", by.y = "placette")
# compute new tree coordinates from coregistered plot center
dummy <- lidaRtRee::polar2Projected(
- ap$X.cor, ap$Y.cor, ap$ZGPS, ap$azimutG / 200 * pi,
+ ap$X_cor, ap$Y_cor, ap$ZGPS, ap$azimutG / 200 * pi,
ap$distR, atan(ap$pente. / 100), 1.55 / 180 * pi,
-2.2 / 180 * pi, 0
)
-ap$X.cor <- dummy$x
-ap$Y.cor <- dummy$y
+ap$X_cor <- dummy$x
+ap$Y_cor <- dummy$y
# save new table
# write.table(round(p.cor,3),file="coregistered_plots.csv", row.names=F,sep=";")
```
@@ -367,26 +367,26 @@ ap$Y.cor <- dummy$y
### Analysis of GPS errors
```{r batch.distance, include = FALSE, eval=TRUE}
-distance <- sqrt((p$X.cor - p$XGPS)^2 + (p$Y.cor - p$YGPS)^2)
+distance <- sqrt((p$X_cor - p$XGPS)^2 + (p$Y_cor - p$YGPS)^2)
```
-Graphics on the difference between initial and corrected positions are displayed. The mean difference is `r round(mean(p$X.cor- p$XGPS), 2)`m in the X axis and `r round(mean(p$Y.cor- p$YGPS),2)`m in the Y axis. `p.values` of a t.test are respectively `r round(t.test(p$X.cor- p$XGPS)$p.value,2)` and `r round(t.test(p$Y.cor- p$YGPS)$p.value, 2)`. Mean absolute distance is `r round(mean(distance), 2)`m with a standard deviation of `r round(sd(distance),2)`m.
+Graphics on the difference between initial and corrected positions are displayed. The mean difference is `r round(mean(p$X_cor- p$XGPS), 2)`m in the X axis and `r round(mean(p$Y_cor- p$YGPS),2)`m in the Y axis. `p.values` of a t.test are respectively `r round(t.test(p$X_cor- p$XGPS)$p.value,2)` and `r round(t.test(p$Y_cor- p$YGPS)$p.value, 2)`. Mean absolute distance is `r round(mean(distance), 2)`m with a standard deviation of `r round(sd(distance),2)`m.
```{r batch.distanceplot, include = TRUE, out.width='80%', fig.width = 8, fig.height = 4}
par(mfrow = c(1, 2))
# plot position difference with additionnal jitter to visualize same points
-plot(jitter(p$X.cor - p$XGPS), jitter(p$Y.cor - p$YGPS),
+plot(jitter(p$X_cor - p$XGPS), jitter(p$Y_cor - p$YGPS),
asp = 1, col = "black",
main = "Corrected-Initial position", xlab = "X difference", ylab = "Y difference"
)
abline(v = 0, lty = 2)
abline(h = 0, lty = 2)
# Distance between initial and corrected
-hist(sqrt((p$X.cor - p$XGPS)^2 + (p$Y.cor - p$YGPS)^2),
+hist(sqrt((p$X_cor - p$XGPS)^2 + (p$Y_cor - p$YGPS)^2),
main = "Histogram of distances",
xlab = "Absolute distance corrected-initial", ylab = "Number of plots"
)
```
-## References
+## References
\ No newline at end of file
diff --git a/R/forest.structure.metrics.Rmd b/R/forest.structure.metrics.Rmd
index 315dfdf93f6db58ede9c0e3fb5dbc9c6b171b5c0..57b043eb5e7a887bae929ce3a58aecdebc4899a3 100644
--- a/R/forest.structure.metrics.Rmd
+++ b/R/forest.structure.metrics.Rmd
@@ -3,8 +3,8 @@ title: "R workflow for forest structure metrics computation from ALS data"
author: "Jean-Matthieu Monnet, with contributions from B. Reineking and A. Glad"
date: "`r Sys.Date()`"
output:
- pdf_document: default
html_document: default
+ pdf_document: default
papersize: a4
bibliography: "../bib/bibliography.bib"
---
@@ -52,37 +52,36 @@ Numerous parameters have to be set for processing.
# output metrics map resolution (m)
resolution <- 10
# canopy height model resolution (m)
-resCHM <- 0.5
+res_chm <- 0.5
# buffer size (m) for tile processing (20 m is better for gaps metrics, 10 m is enough for
# tree metrics)
-b.size <- 20
+buffer_size <- 20
# height threshold (m) for high points removal (points above this threshold are considered
# as outliers)
-h.points <- 60
+points_max_h <- 60
# points classes to retain for analysis (vegetation + ground)
-class.points <- c(0, 1, 2, 3, 4, 5)
+class_points <- c(0, 1, 2, 3, 4, 5)
# ground class
-class.ground <- 2
+class_ground <- 2
# Gaussian filter sigma values for multi-scale smoothing
sigma <- c(0, 1, 2, 4, 8, 16)
# convert vector of sigma values in list
-l.sigma <- as.list(sigma)
+sigma_l <- as.list(sigma)
# height breaks for penetration ratio and density
-breaksH <- c(-Inf, 0, 0.5, 1, 2, 5, 10, 20, 30, 60, Inf)
+breaks_h <- c(-Inf, 0, 0.5, 1, 2, 5, 10, 20, 30, 60, Inf)
# percentiles of height distribution
percent <- c(0.10, 0.25, 0.5, 0.75, 0.9)
# surface breaks for gap size (m2)
-breaksGapSurface <- c(4, 16, 64, 256, 1024, 4096, Inf)
+breaks_gap_surface <- c(4, 16, 64, 256, 1024, 4096, Inf)
#
# gap surface names
-n.breaksG <- gsub("-", "", paste("G.s", breaksGapSurface[c(-length(breaksGapSurface))], "_",
- breaksGapSurface[c(-1)],
- sep = ""
+n_breaks_gap <- gsub("-", "", paste0("G_s", breaks_gap_surface[c(-length(breaks_gap_surface))], "to",
+ breaks_gap_surface[c(-1)]
))
# height bin names
-n.breaksH <- gsub("-", "", paste("H.nb", breaksH[c(-length(breaksH))], "_", breaksH[c(-1)],
- sep = ""
+n_breaks_h <- gsub("-", "", paste0("nb_H", breaks_h[c(-length(breaks_h))], "to", breaks_h[c(-1)]
))
+#
```
The first step is to create a catalog of LAS files (should be normalized, non-overlapping rectangular tiles). Preferably, tiles should be aligned on a multiple of resolution, and points should not lie on the northern or eastern border when such borders are common with adjacent tiles.
@@ -106,15 +105,15 @@ 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 <- as.numeric(cata@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# load tile extent plus buffer
a <- lidR::clip_rectangle(
- cata, b.box[1] - b.size, b.box[2] - b.size, b.box[3] + b.size,
- b.box[4] + b.size
+ cata, b_box[1] - buffer_size, b_box[2] - buffer_size, b_box[3] + buffer_size,
+ b_box[4] + buffer_size
)
# add 'buffer' flag to points in buffer with TRUE value in this new attribute
a <- lidR::add_attribute(
- a, a$X < b.box[1] | a$Y < b.box[2] | a$X >= b.box[3] | a$Y >= b.box[4],
+ a, a$X < b_box[1] | a$Y < b_box[2] | a$X >= b_box[3] | a$Y >= b_box[4],
"buffer"
)
a
@@ -123,7 +122,7 @@ a
Only points from desired classes are retained. In case some mis-classified, high or low points remain in the data set, the point cloud is filtered and negative heights are replaced by 0.
```{r filterPointCloud, include = TRUE, message=FALSE, warnings=FALSE}
# remove unwanted point classes, and points higher than height threshold
-a <- lidR::filter_poi(a, is.element(Classification, class.points) & Z <= h.points)
+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
#
@@ -142,9 +141,9 @@ 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 = resCHM)
+chm <- lidaRtRee::points2DSM(a, res = res_chm)
# replace NA, low and high values
-chm[is.na(chm) | chm < 0 | chm > h.points] <- 0
+chm[is.na(chm) | chm < 0 | chm > points_max_h] <- 0
# display CHM
raster::plot(chm, asp = 1, main = "Canopy height model")
```
@@ -157,8 +156,8 @@ The CHM is smoothed with a Gaussian filter with different `sigma` values. Smooth
```{r 2dchmMetrics, include = TRUE, fig.dim = c(9, 4.5), out.width='100%', message=FALSE}
# for each value in list of sigma, apply filtering to chm and store result in list
-st <- lapply(l.sigma, FUN = function(x) {
- lidaRtRee::demFiltering(chm, nlFilter = "Closing", nlSize = 5, sigmap = x)[[2]]
+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)
@@ -170,64 +169,64 @@ The raster stack is then converted to points. Points outside the tile extent are
```{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 <- raster::crop(st, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4]))
# multiplying factor to compute proportion
-mf <- 100 / (resolution / resCHM)^2
+mf <- 100 / (resolution / res_chm)^2
# modify layer names
names(st) <- paste0("layer.", 1:6)
# compute raster metrics
-metrics.2dchm <- lidaRtRee::rasterMetrics(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
- )
- },
- output = "raster"
+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
+ )
+ },
+ 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"
+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
+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(
- "CHM0.sd", "CHM2.sd", "CHM8.sd", "CHM.mean",
- "CHM.PercInf1", "CHM.PercSup10"
+raster::plot(metrics_2dchm[[c(
+ "CHM0_sd", "CHM2_sd", "CHM8_sd", "CHM_mean",
+ "CHM_PercInf1", "CHM_PercSup10"
)]])
```
### Gaps and edges metrics
-Gaps are computed with the function `gapDetection`.
+Gaps are computed with the function `gap_detection`.
```{r gapMetrics, include = TRUE, fig.dim = c(9, 3), out.width='100%', warning=FALSE, message=FALSE}
# compute gaps
-gaps <- lidaRtRee::gapDetection(chm,
- ratio = 2, gap.max.height = 1,
- min.gap.surface = min(breaksGapSurface), closing.height.bin = 1,
- nlFilter = "Median", nlsize = 3, gapReconstruct = TRUE
+gaps <- lidaRtRee::gap_detection(chm,
+ ratio = 2, gap_max_height = 1,
+ min_gap_surface = min(breaks_gap_surface), closing_height_bin = 1,
+ nl_filter = "Median", nl_size = 3, gap_reconstruct = TRUE
)
# display maps
par(mfrow = c(1, 3))
raster::plot(chm, main = "CHM")
-raster::plot(log(gaps$gap.surface), main = "Log of gap surface")
+raster::plot(log(gaps$gap_surface), main = "Log of gap surface")
# display gaps
-dummy <- gaps$gap.id
+dummy <- gaps$gap_id
dummy[dummy == 0] <- NA
raster::plot(dummy %% 8, asp = 1, main = "Gap id (random colors)", col = rainbow(8), legend = FALSE)
```
@@ -235,46 +234,46 @@ Summary statistics are then computed: pixel surface occupied by gaps in differen
```{r gapStats, include = TRUE, warning=FALSE, message=FALSE}
# crop results to tile size
-gaps.surface <- raster::crop(gaps$gap.surface, raster::extent(
- b.box[1], b.box[3],
- b.box[2], b.box[4]
+gaps_surface <- raster::crop(gaps$gap_surface, raster::extent(
+ 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)))) {
- metrics.gaps <- lidaRtRee::rasterMetrics(gaps.surface,
- res = resolution,
- fun = function(x) {
- hist(x$layer,
- breaks = breaksGapSurface,
- plot = F
- )$counts * (resCHM / resolution)^2
- }, output = "raster"
+if (!all(is.na(raster::values(gaps_surface)))) {
+ metrics_gaps <- lidaRtRee::raster_metrics(gaps_surface,
+ res = resolution,
+ fun = function(x) {
+ hist(x$layer,
+ 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))),
- fun = sum
+ metrics_gaps$sum <- raster::stackApply(metrics_gaps, rep(1, length(names(metrics_gaps))),
+ fun = sum
)
# if gaps are present
- names(metrics.gaps) <- c(n.breaksG, paste("G.s", min(breaksGapSurface), "_Inf", sep = ""))
+ names(metrics_gaps) <- c(n_breaks_gap, paste("G_s", min(breaks_gap_surface), "toInf", sep = ""))
# replace possible NA values by 0 (NA values are present in lines of the target raster
# where no values >0 are present in the input raster)
- metrics.gaps[is.na(metrics.gaps)] <- 0
+ metrics_gaps[is.na(metrics_gaps)] <- 0
} else {
- metrics.gaps <- NULL
+ metrics_gaps <- NULL
}
#
-metrics.gaps
+metrics_gaps
```
-Edges are detected with the function `edgeDetection` as the outer envelope of the previously delineated gaps.
+Edges are detected with the function `edge_detection` as the outer envelope of the previously delineated gaps.
```{r edgeDetection, include = TRUE, fig.dim = c(9, 3), out.width='100%', warning=FALSE, message=FALSE}
# Perform edge detection by erosion
-edges <- lidaRtRee::edgeDetection(gaps$gap.id > 0)
+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(gaps$gap_id > 0, main = "Gaps", legend = FALSE)
raster::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.
@@ -283,35 +282,35 @@ The percentage of surface occupied by edges is then computed as summary statisti
# crop results to tile size
edges <- raster::crop(
edges,
- raster::extent(b.box[1], b.box[3], b.box[2], b.box[4])
+ raster::extent(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)))) {
- metrics.edges <- lidaRtRee::rasterMetrics(edges,
- res = resolution,
- fun = function(x) {
- sum(x$layer) * (resCHM / resolution)^2
- }, output = "raster"
+ metrics_edges <- lidaRtRee::raster_metrics(edges,
+ res = resolution,
+ fun = function(x) {
+ sum(x$layer) * (res_chm / resolution)^2
+ }, output = "raster"
)
- names(metrics.edges) <- "edges.proportion"
+ names(metrics_edges) <- "edges.proportion"
# replace possible NA values by 0 (NA values are present in lines of the target raster
# where no values >0 are present in the input raster)
- metrics.edges[is.na(metrics.edges)] <- 0
+ metrics_edges[is.na(metrics_edges)] <- 0
} else {
- metrics.edges <- NULL
+ metrics_edges <- NULL
}
#
-metrics.edges
+metrics_edges
```
Some outputs are displayed hereafter. The proportion of surface occupied by gaps of various sizes is depicted as well as the proportion of surface occupied by edges.
```{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.s4_16"]], main="Prop of gaps 4-16m2")
-raster::plot(metrics.gaps[["G.s64_256"]], main = "Prop of gaps 64-256m2")
-raster::plot(metrics.gaps[["G.s4_Inf"]], main = "Prop of gaps >4m2")
-raster::plot(metrics.edges, main = "Proportion of edges")
+# 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")
```
### Tree metrics
@@ -319,75 +318,75 @@ raster::plot(metrics.edges, main = "Proportion of edges")
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}
# tree top detection (default parameters)
-segms <- lidaRtRee::treeSegmentation(chm, hmin = 5)
+segms <- lidaRtRee::tree_segmentation(chm, hmin = 5)
# extraction to data.frame
-trees <- lidaRtRee::treeExtraction(segms$filled.dem, segms$local.maxima, segms$segments.id)
+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")
points(trees$x, trees$y, pch = 4, cex = 0.3)
# display segments, except ground segment
-dummy <- segms$segments.id
+dummy <- segms$segments_id
dummy[dummy == 0] <- NA
raster::plot(dummy %% 8, asp = 1, main = "Segments (random colors)", col = rainbow(8), legend = FALSE)
points(trees$x, trees$y, pch = 4, cex = 0.3)
```
-Results are cropped to the tile extent. Summary statistics from `stdTreeMetrics` are then computed for each pixel. Canopy cover in trees and mean canopy height in trees are computed in a second step because they can not be computed from the data of the extracted trees which crowns may span several output pixels. They are calculated from the images obtained during the previous tree segmentation.
+Results are cropped to the tile extent. Summary statistics from `std_tree_metrics` are then computed for each pixel. Canopy cover in trees and mean canopy height in trees are computed in a second step because they can not be computed from the data of the extracted trees which crowns may span several output pixels. They are calculated from the images obtained during the previous tree segmentation.
```{r cropTreeMetrics, include = TRUE, warning=FALSE, message=FALSE}
# remove trees outside of tile
-trees <- trees[trees$x >= b.box[1] & trees$x < b.box[3] & trees$y >= b.box[2] & trees$y < b.box[4], ]
+trees <- trees[trees$x >= b_box[1] & trees$x < b_box[3] & trees$y >= b_box[2] & trees$y < b_box[4], ]
# compute raster metrics
-metrics.trees <- lidaRtRee::rasterMetrics(trees[, -1],
- res = resolution,
- fun = function(x) {
- lidaRtRee::stdTreeMetrics(x, resolution^2 / 10000)
- }, output = "raster"
+metrics_trees <- lidaRtRee::raster_metrics(trees[, -1],
+ res = resolution,
+ fun = function(x) {
+ lidaRtRee::std_tree_metrics(x, resolution^2 / 10000)
+ }, output = "raster"
)
# remove NAs and NaNs
-metrics.trees[!is.finite(metrics.trees)] <- 0
+metrics_trees[!is.finite(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.treechm <- segms$filled.dem
+r_tree_chm <- segms$filled_dem
# set chm to NA in non segment area
-r.treechm[segms$segments.id == 0] <- NA
+r_tree_chm[segms$segments_id == 0] <- NA
# compute raster metrics
-dummy <- lidaRtRee::rasterMetrics(raster::crop(
- r.treechm,
+dummy <- lidaRtRee::raster_metrics(raster::crop(
+ r_tree_chm,
raster::extent(
- b.box[1], b.box[3], b.box[2],
- b.box[4]
+ 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 / resCHM)^2,
- mean(x$filled.dem, na.rm = T)
+ sum(!is.na(x$filled_dem)) / (resolution / res_chm)^2,
+ mean(x$filled_dem, na.rm = T)
)
},
output = "raster"
)
-names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot")
+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 <- raster::extend(dummy, metrics_trees)
+metrics_trees <- raster::addLayer(metrics_trees, dummy)
#
-metrics.trees
+metrics_trees
```
Some outputs are displayed hereafter.
```{r displayTreeMetrics, include = TRUE, fig.width = 9, fig.height = 3.2, warning=FALSE, message=FALSE}
par(mfrow = c(1, 3))
-raster::plot(metrics.trees$Tree.meanH, main = "Mean (detected) tree height (m)")
+raster::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,
- main = "Density of (detected) trees > 20m (/ha)"
+raster::plot(metrics_trees$TreeSup20_density,
+ main = "Density of (detected) trees > 20m (/ha)"
)
-raster::plot(metrics.trees$Tree.meanCrownVolume,
- main = "Mean crown volume of detected trees (m3)"
+raster::plot(metrics_trees$Tree_meanCrownVolume,
+ main = "Mean crown volume of detected trees (m3)"
)
```
@@ -400,68 +399,68 @@ 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::grid_metrics(a, as.list(c(
as.vector(quantile(Z, probs = percent)),
sum(Classification == 2),
mean(Intensity[ReturnNumber == 1])
)), res = resolution)
-names(metrics.1d) <- c(paste("H.p", percent * 100, sep = ""), "nbSol", "I.mean.1stpulse")
+names(metrics_1d) <- c(paste("Hp", percent * 100, sep = ""), "nb_ground", "Imean_1stpulse")
#
# vegetation-only metrics
-a <- lidR::filter_poi(a, Classification != class.ground)
+a <- lidR::filter_poi(a, Classification != class_ground)
dummy <- lidR::grid_metrics(a, as.list(c(
mean(Z),
max(Z),
sd(Z),
length(Z),
- hist(Z, breaks = breaksH, right = F, plot = F)$counts
+ hist(Z, breaks = breaks_h, right = F, plot = F)$counts
)), res = resolution)
-names(dummy) <- c("H.mean", "H.max", "H.sd", "nbVeg", n.breaksH)
+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 <- raster::extend(dummy, metrics_1d)
+metrics_1d <- raster::addLayer(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 <- raster::crop(metrics_1d, raster::extent(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$H.simpson <- raster::stackApply(
- metrics.1d[[n.breaksH[c(-1, -length(n.breaksH))]]],
+metrics_1d$Hsimpson <- raster::stackApply(
+ metrics_1d[[n_breaks_h[c(-1, -length(n_breaks_h))]]],
1,
function(x, ...) {
vegan::diversity(x, index = "simpson")
}
)
-# Relative density
-for (i in n.breaksH[c(-1, -length(n.breaksH))])
+# Relative density
+for (i in n_breaks_h[c(-1, -length(n_breaks_h))])
{
- metrics.1d[[paste0(i, "relative_density")]] <-
- metrics.1d[[i]] / (metrics.1d$nbVeg + metrics.1d$nbSol)
+ metrics_1d[[gsub("nb", "relativeDensity", i)]] <-
+ metrics_1d[[i]] / (metrics_1d$nb_veg + metrics_1d$nb_ground)
}
# Penetration ratio
# cumulative sum
-cumu.sum <- metrics.1d[["nbSol"]]
-for (i in n.breaksH)
+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 <- raster::addLayer(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics_1d[[i]])
}
-names(cumu.sum) <- c("nbSol", n.breaksH)
+names(cumu_sum) <- c("nb_ground", n_breaks_h)
# compute interception ratio of each layer
-intercep.ratio <- cumu.sum[[-1]]
-for (i in 2:dim(cumu.sum)[3])
+intercep_ratio <- cumu_sum[[-1]]
+for (i in 2:dim(cumu_sum)[3])
{
- intercep.ratio[[i - 1]] <- 1 - cumu.sum[[i - 1]] / cumu.sum[[i]]
+ intercep_ratio[[i - 1]] <- 1 - cumu_sum[[i - 1]] / cumu_sum[[i]]
}
-names(intercep.ratio) <- paste0(names(cumu.sum)[-1], "ratio")
+names(intercep_ratio) <- gsub("nb", "intercepRatio", names(cumu_sum)[-1])
# merge rasterstacks
-metrics.1d <- raster::addLayer(metrics.1d, intercep.ratio)
-# rm(cumu.sum, intercep.ratio)
+metrics_1d <- raster::addLayer(metrics_1d, intercep_ratio)
+# rm(cumu_sum, intercep_ratio)
#
-metrics.1d
+metrics_1d
```
Some outputs are displayed hereafter.
@@ -469,9 +468,9 @@ Some outputs are displayed hereafter.
```{r display1dmetrics, include = TRUE, fig.width = 9, fig.height = 3.2, warning=FALSE, message=FALSE}
# display results
par(mfrow = c(1, 3))
-raster::plot(metrics.1d$I.mean.1stpulse, main = "Mean intensity (1st return)")
-raster::plot(metrics.1d$H.p50, main = "Percentile 50 of heights")
-raster::plot(metrics.1d$H.nb2_5ratio, main = "Interception ratio 2-5 m")
+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")
```
## Batch processing
@@ -482,285 +481,283 @@ The following code uses parallel processing to handle multiple files of a catalo
# processing by tile
metrics <- foreach::foreach(i = 1:nrow(cata), .errorhandling = "remove") %dopar% {
# tile extent
- b.box <- as.numeric(cata@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
+ b_box <- as.numeric(cata@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# load tile extent plus buffer
a <- try(lidR::clip_rectangle(
- cata, b.box[1] - b.size, b.box[2] - b.size, b.box[3] + b.size,
- b.box[4] + b.size
+ cata, b_box[1] - buffer_size, b_box[2] - buffer_size, b_box[3] + buffer_size,
+ b_box[4] + buffer_size
))
#
- if (class(a) == "try-error") # check if points are successfully loaded
- {
- # else NULL output
- temp <- NULL
- } else {
- # add 'buffer' flag to points in buffer with TRUE value in this new attribute
- a <- lidR::add_attribute(
- a,
- a$X < b.box[1] | a$Y < b.box[2] | a$X >= b.box[3] | a$Y >= b.box[4],
- "buffer"
- )
- # remove unwanted point classes, and points higher than height threshold
- a <- lidR::filter_poi(a, is.element(Classification, class.points) & Z <= h.points)
- # check number of remaining points
- if (length(a) == 0) {
- # output NULL
- temp <- NULL
- } else {
- # set negative heights to 0
- a@data$Z[a@data$Z < 0] <- 0
- #
- # compute chm
- chm <- lidaRtRee::points2DSM(a, res = resCHM)
- # replace NA, low and high values
- chm[is.na(chm) | chm < 0 | chm > h.points] <- 0
- # ##########################
- # compute 2D CHM metrics
- # for each value in list of sigma, apply filtering to chm and store result in list
- st <- lapply(l.sigma, FUN = function(x) {
- lidaRtRee::demFiltering(chm, nlFilter = "Closing", nlSize = 5, sigmap = x)[[2]]
- })
- # convert to raster
- st <- raster::stack(st)
- # crop to bbox
- st <- raster::crop(st, raster::extent(b.box[1], b.box[3], b.box[2], b.box[4]))
- # multiplying factor to compute proportion
- mf <- 100 / (resolution / resCHM)^2
- # compute raster metrics
- metrics.2dchm <-
- lidaRtRee::rasterMetrics(
- 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
- )
- },
- output = "raster"
- )
- #
- # ##########################
- # compute gap metrics
- gaps <- lidaRtRee::gapDetection(chm,
- ratio = 2, gap.max.height = 1,
- min.gap.surface = min(breaksGapSurface),
- closing.height.bin = 1,
- nlFilter = "Median", nlsize = 3,
- gapReconstruct = TRUE
- )
- # crop results to tile size
- gaps.surface <- raster::crop(gaps$gap.surface, raster::extent(
- 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)))) {
- metrics.gaps <- lidaRtRee::rasterMetrics(gaps.surface,
- res = resolution,
- fun = function(x) {
- hist(x$layer,
- breaks = breaksGapSurface,
- plot = F
- )$counts * (resCHM / resolution)^2
- },
- output = "raster"
- )
- # compute total gap proportion
- metrics.gaps$sum <- raster::stackApply(metrics.gaps, rep(1, length(names(metrics.gaps))),
- fun = sum
- )
- # if gaps are present
- names(metrics.gaps) <- c(n.breaksG, paste("G.s", min(breaksGapSurface), "_Inf", sep = ""))
- # replace possible NA values by 0 (NA values are present in lines of the target raster
- # where no values >0 are present in the input raster)
- metrics.gaps[is.na(metrics.gaps)] <- 0
- } else {
- metrics.gaps <- NULL
- }
- #
- #############################
- # compute edge metrics
- # Perform edge detection by erosion
- edges <- lidaRtRee::edgeDetection(gaps$gap.id > 0)
- # crop results to tile size
- edges <- raster::crop(edges, raster::extent(
- 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)))) {
- metrics.edges <- lidaRtRee::rasterMetrics(edges,
- res = resolution,
- fun = function(x) {
- sum(x$layer) * (resCHM / resolution)^2
- }, output = "raster"
+ # check if points are successfully loaded
+ if (class(a) == "try-error") {
+ return(NULL)
+ }
+ # add 'buffer' flag to points in buffer with TRUE value in this new attribute
+ a <- lidR::add_attribute(
+ a,
+ a$X < b_box[1] | a$Y < b_box[2] | a$X >= b_box[3] | a$Y >= b_box[4],
+ "buffer"
+ )
+ # remove unwanted point classes, and points higher than height threshold
+ a <- lidR::filter_poi(a, is.element(Classification, class_points) & Z <= points_max_h)
+ # check number of remaining points
+ if (length(a) == 0) {
+ return(NULL)
+ }
+ # set negative heights to 0
+ a@data$Z[a@data$Z < 0] <- 0
+ #
+ # compute chm
+ chm <- lidaRtRee::points2DSM(a, res = res_chm)
+ # replace NA, low and high values
+ chm[is.na(chm) | chm < 0 | chm > points_max_h] <- 0
+ #-----------------------
+ # compute 2D CHM metrics
+ # 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)[[2]]
+ })
+ # convert to raster
+ st <- raster::stack(st)
+ # crop to bbox
+ st <- raster::crop(st, raster::extent(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
+ metrics_2dchm <-
+ 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
)
- names(metrics.edges) <- "edges.proportion"
- # replace possible NA values by 0 (NA values are present in lines of the target raster
- # where no values >0 are present in the input raster)
- metrics.edges[is.na(metrics.edges)] <- 0
- } else {
- metrics.edges <- NULL
- }
- #
- # ###########################
- # compute tree metrics
- # tree top detection (default parameters)
- segms <- lidaRtRee::treeSegmentation(chm, hmin = 5)
- # extraction to data.frame
- trees <- lidaRtRee::treeExtraction(segms$filled.dem, segms$local.maxima, segms$segments.id)
- # remove trees outside of tile
- trees <- trees[trees$x >= b.box[1] & trees$x < b.box[3] & trees$y >= b.box[2] & trees$y < b.box[4], ]
- # compute raster metrics
- metrics.trees <- lidaRtRee::rasterMetrics(trees[, -1],
- res = resolution,
- fun = function(x) {
- lidaRtRee::stdTreeMetrics(x, resolution^2 / 10000)
- },
- output = "raster"
+ },
+ output = "raster"
+ )
+ #
+ # -------------------
+ # compute gap metrics
+ gaps <- lidaRtRee::gap_detection(chm,
+ ratio = 2, gap_max_height = 1,
+ min_gap_surface = min(breaks_gap_surface),
+ closing_height_bin = 1,
+ nl_filter = "Median", nl_size = 3,
+ 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],
+ b_box[4]
+ ))
+ # compute gap statistics at final metrics resolution, in case gaps exist
+ if (!all(is.na(raster::values(gaps_surface)))) {
+ metrics_gaps <- lidaRtRee::raster_metrics(gaps_surface,
+ res = resolution,
+ fun = function(x) {
+ hist(x$layer,
+ 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))),
+ fun = sum
+ )
+ # if gaps are present
+ names(metrics_gaps) <- c(n_breaks_gap, paste("G_s", min(breaks_gap_surface), "toInf", sep = ""))
+ # replace possible NA values by 0 (NA values are present in lines of the target raster
+ # where no values >0 are present in the input raster)
+ metrics_gaps[is.na(metrics_gaps)] <- 0
+ } else {
+ metrics_gaps <- NULL
+ }
+ #
+ # --------------------
+ # compute edge metrics
+ # 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],
+ b_box[2], b_box[4]
+ ))
+ # compute edges statistics at final metrics resolution, in case edges exist
+ if (!all(is.na(raster::values(edges)))) {
+ metrics_edges <- lidaRtRee::raster_metrics(edges,
+ res = resolution,
+ fun = function(x) {
+ sum(x$layer) * (res_chm / resolution)^2
+ }, output = "raster"
+ )
+ names(metrics_edges) <- "edges.proportion"
+ # replace possible NA values by 0 (NA values are present in lines of the target raster
+ # where no values >0 are present in the input raster)
+ metrics_edges[is.na(metrics_edges)] <- 0
+ } else {
+ metrics_edges <- NULL
+ }
+ #
+ # --------------------
+ # compute tree metrics
+ # 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)
+ # remove trees outside of tile
+ trees <- trees[trees$x >= b_box[1] & trees$x < b_box[3] & trees$y >= b_box[2] & trees$y < b_box[4], ]
+ # compute raster metrics
+ metrics_trees <- lidaRtRee::raster_metrics(trees[, -1],
+ res = resolution,
+ fun = function(x) {
+ lidaRtRee::std_tree_metrics(x, resolution^2 / 10000)
+ },
+ output = "raster"
+ )
+ # remove NAs and NaNs
+ metrics_trees[!is.finite(metrics_trees)] <- 0
+ #
+ # extend layer in case there are areas without trees
+ metrics_trees <- raster::extend(metrics_trees, metrics_2dchm, values = 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
+ # 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(r_tree_chm, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4])),
+ res = resolution,
+ fun = function(x) {
+ c(
+ sum(!is.na(x$filled_dem)) / (resolution / res_chm)^2,
+ mean(x$filled_dem, na.rm = T)
)
- # remove NAs and NaNs
- metrics.trees[!is.finite(metrics.trees)] <- 0
+ },
+ output = "raster"
+ )
+ names(dummy) <- c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
+ dummy <- raster::extend(dummy, metrics_2dchm, values = 0)
+ #
+ metrics_trees <- raster::addLayer(metrics_trees, dummy)
+ #
+ # -------------------------
+ # compute 1D height metrics
+ # remove buffer points
+ a <- lidR::filter_poi(a, buffer == 0)
+ #
+ if (length(a) == 0) {
+ metrics_1d <- NULL
+ } else {
+ # all points metrics
+ metrics_1d <- lidR::grid_metrics(a, as.list(c(
+ as.vector(quantile(Z, probs = percent)),
+ sum(Classification == 2),
+ mean(Intensity[ReturnNumber == 1])
+ )), res = resolution)
+ names(metrics_1d) <- c(paste("Hp", percent * 100, sep = ""), "nb_ground", "Imean_1stpulse")
+ #
+ # vegetation-only metrics
+ a <- lidR::filter_poi(a, Classification != class_ground)
+ if (length(a) != 0) {
+ dummy <- lidR::grid_metrics(a, as.list(c(
+ mean(Z),
+ max(Z),
+ sd(Z),
+ length(Z),
+ hist(Z, breaks = breaks_h, right = F, plot = F)$counts
+ )), 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)
+ rm(dummy)
+ # -------------
+ # merge metrics
+ metrics_1d <- raster::extend(metrics_1d, metrics_2dchm, values = 0)
+ metrics_1d <- raster::crop(metrics_1d, metrics_2dchm)
+ metrics_gaps <- raster::extend(metrics_gaps, metrics_2dchm, values = 0)
+ metrics_gaps <- raster::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)
#
- # extend layer in case there are areas without trees
- metrics.trees <- raster::extend(metrics.trees, metrics.2dchm, values = 0)
- # compute canopy cover in trees and canopy mean height in trees
- # in region of interest, because it is not in previous step.
- r.treechm <- segms$filled.dem
- # set chm to NA in non segment area
- r.treechm[segms$segments.id == 0] <- NA
- # compute raster metrics
- dummy <- lidaRtRee::rasterMetrics(
- raster::crop(r.treechm, raster::extent(b.box[1], b.box[3], b.box[2], b.box[4])),
- res = resolution,
- fun = function(x) {
- c(
- sum(!is.na(x$filled.dem)) / (resolution / resCHM)^2,
- mean(x$filled.dem, na.rm = T)
- )
- },
- output = "raster"
+ temp <- raster::addLayer(
+ metrics_1d, metrics_2dchm, metrics_gaps,
+ metrics_edges, metrics_trees
)
- names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot")
- dummy <- raster::extend(dummy, metrics.2dchm, values = 0)
- #
- metrics.trees <- raster::addLayer(metrics.trees, dummy)
- #
- ############################
- # compute 1D height metrics
- # remove buffer points
- a <- lidR::filter_poi(a, buffer == 0)
- #
- if (length(a) == 0) {
- metrics.1d <- NULL
- } else {
- # all points metrics
- metrics.1d <- lidR::grid_metrics(a, as.list(c(
- as.vector(quantile(Z, probs = percent)),
- sum(Classification == 2),
- mean(Intensity[ReturnNumber == 1])
- )), res = resolution)
- names(metrics.1d) <- c(paste("H.p", percent * 100, sep = ""), "nbSol", "I.mean.1stpulse")
- #
- # vegetation-only metrics
- a <- lidR::filter_poi(a, Classification != class.ground)
- if (length(a) != 0) {
- dummy <- lidR::grid_metrics(a, as.list(c(
- mean(Z),
- max(Z),
- sd(Z),
- length(Z),
- hist(Z, breaks = breaksH, right = F, plot = F)$counts
- )), res = resolution)
- names(dummy) <- c("H.mean", "H.max", "H.sd", "nbVeg", n.breaksH)
- # merge rasterstacks
- dummy <- raster::extend(dummy, metrics.1d)
- metrics.1d <- raster::addLayer(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)
- metrics.gaps <- raster::extend(metrics.gaps, metrics.2dchm, values = 0)
- metrics.gaps <- raster::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)
- #
- temp <- raster::addLayer(
- metrics.1d, metrics.2dchm, metrics.gaps,
- metrics.edges, metrics.trees
- )
- temp[is.na(temp)] <- 0
- } # end of vegetation-only points presence check
- } # end of buffer-less points presence check
- } # end of buffered area points presence check
- } # success of file reading check
- return(temp)
+ temp[is.na(temp)] <- 0
+ } # end of vegetation-only points presence check
+ } # end of buffer-less points presence check
+return(temp)
}
-# ###########################
+# --------------
# mosaic rasters
-names.metrics <- names(metrics[[1]])
+names_metrics <- names(metrics[[1]])
metrics <- do.call(raster::merge, metrics)
-names(metrics) <- names.metrics
-# ###########################
+names(metrics) <- names_metrics
+# --------------------------
# compute additional metrics
# Simpson index
-metrics$H.simpson <- raster::stackApply(
- metrics[[n.breaksH[c(-1, -length(n.breaksH))]]], 1,
+metrics$Hsimpson <- raster::stackApply(
+ metrics[[n_breaks_h[c(-1, -length(n_breaks_h))]]], 1,
function(x, ...) {
vegan::diversity(x, index = "simpson")
}
)
#
# Relative density
-for (i in n.breaksH[c(-1, -length(n.breaksH))])
+for (i in n_breaks_h[c(-1, -length(n_breaks_h))])
{
- metrics[[paste0(i, "relative_density")]] <- metrics[[i]] / (metrics$nbVeg + metrics$nbSol)
+ metrics[[gsub("nb", "relativeDensity", i)]] <-
+ metrics[[i]] / (metrics$nb_veg + metrics$nb_ground)
}
# Penetration ratio
# compute cumulative sum
-cum.sum <- metrics[["nbSol"]]
-for (i in n.breaksH)
+cumu_sum <- metrics[["nb_ground"]]
+for (i in n_breaks_h)
{
- cum.sum <- raster::addLayer(cum.sum, cum.sum[[dim(cum.sum)[3]]] + metrics[[i]])
+ cumu_sum <- raster::addLayer(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics[[i]])
}
-names(cum.sum) <- c("nbSol", n.breaksH)
+names(cumu_sum) <- c("nb_ground", n_breaks_h)
# interception ratio
-intercep.ratio <- cum.sum[[-1]]
-for (i in 2:dim(cum.sum)[3])
+intercep_ratio <- cumu_sum[[-1]]
+for (i in 2:dim(cumu_sum)[3])
{
- intercep.ratio[[i - 1]] <- 1 - cum.sum[[i - 1]] / cum.sum[[i]]
+ intercep_ratio[[i - 1]] <- 1 - cumu_sum[[i - 1]] / cumu_sum[[i]]
}
-names(intercep.ratio) <- paste0(names(cum.sum)[-1], "ratio")
+names(intercep_ratio) <- gsub("nb", "intercepRatio", names(cumu_sum)[-1])
# merge rasterstacks
-metrics <- raster::addLayer(metrics, intercep.ratio)
+metrics <- raster::addLayer(metrics, intercep_ratio)
#
-# ###########################
+# ----------------------
# export as raster files
# for (i in 1:names(metrics))
# {
# print(i)
# writeRaster(metrics[[i]],file=paste("output/raster_",i,"_",resolution,"m.tif",sep="")
# }
-# ###########################
+# -------
# display
-raster::plot(metrics[[c("I.mean.1stpulse", "H.simpson")]],
- main = c("Mean intensity of 1st pulse", "Simpson indice of point heights")
+raster::plot(metrics[[c("Imean_1stpulse", "Hsimpson")]],
+ main = c("Mean intensity of 1st pulse",
+ "Simpson indice of point heights")
)
```
-## References
+## References
\ No newline at end of file
diff --git a/R/gaps.edges.detection.Rmd b/R/gaps.edges.detection.Rmd
index 85f5f608be68de7d912adcc2e22ee7fe4207c029..a6d8809117435b95ec42f6174a9a5652fd567fd8 100644
--- a/R/gaps.edges.detection.Rmd
+++ b/R/gaps.edges.detection.Rmd
@@ -85,9 +85,9 @@ The procedure of gap detection is exemplified in @Glad20.
```{r gaps, include = TRUE, fig.width = 6.5, fig.height = 6, warning=FALSE}
# Perform gap detection with gap width criterion
-gaps1 <- lidaRtRee::gapDetection(chm, ratio = 2, gap.max.height = 1, min.gap.surface = 100)
+gaps1 <- lidaRtRee::gap_detection(chm, ratio = 2, gap_max_height = 1, min_gap_surface = 100)
# Perform gap detection without gap width nor min gap surface criteria
-gaps2 <- lidaRtRee::gapDetection(chm, ratio = NULL, gap.max.height = 1, min.gap.surface = 0)
+gaps2 <- lidaRtRee::gap_detection(chm, ratio = NULL, gap_max_height = 1, min_gap_surface = 0)
```
The following figure displays the obtained gaps, colored by size.
@@ -95,14 +95,14 @@ The following figure displays the obtained gaps, colored by size.
# 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(raster::values(gaps1$gap_surface)))
# plot gap surface
-raster::plot(log(gaps1$gap.surface),
- main = "log(Gap surface) (1)", zlim = c(0, z.lim),
+raster::plot(log(gaps1$gap_surface),
+ main = "log(Gap surface) (1)", zlim = c(0, z_lim),
col = rainbow(255, end = 5 / 6)
)
-raster::plot(log(gaps2$gap.surface),
- main = "log(Gap surface) (2)", zlim = c(0, z.lim),
+raster::plot(log(gaps2$gap_surface),
+ main = "log(Gap surface) (2)", zlim = c(0, z_lim),
col = rainbow(255, end = 5 / 6)
)
```
@@ -111,32 +111,32 @@ 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 <- raster::rasterToPolygons(gaps1$gap_id, dissolve = TRUE)
# display gaps over original CHM
raster::plot(chm, main = "Gaps over CHM (1)")
-sp::plot(gaps1.v, add = T)
+sp::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(raster::values(gaps2$gap_id))))[, -1]
# change column names
-names(gaps.df) <- c("id", "surface")
+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 * raster::res(chm)[1]^2
# data.frame
-head(gaps.df, n = 3)
+head(gaps_df, n = 3)
# cumulative distribution of gap surface by surface
-Hmisc::Ecdf(log(gaps.df$surface), weights = gaps.df$surface, xlab = "log(surface)")
+Hmisc::Ecdf(log(gaps_df$surface), weights = gaps_df$surface, xlab = "log(surface)")
```
Histogram can be used to compare the distribution of surface in different classes of gap surface. For better visualization of the differences between the set of criteria, gaps bigger than 500 m^2^ are not taken into account.
```{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 <- raster::values(gaps1$gap_surface)
+surface2 <- raster::values(gaps2$gap_surface)
# remove large gaps
surface1[surface1 > 500] <- NA
surface2[surface2 > 500] <- NA
@@ -152,38 +152,38 @@ In the previous part, a pixel that fulfills the maximum height criterion but doe
```{r gapsReconstruct, include = TRUE, out.width='100%', fig.width = 13, fig.height = 4.5, warning = FALSE}
# load canopy height model
-chm <- lidaRtRee::chmchablais3
+chm <- lidaRtRee::chm_chablais3
chm[is.na(chm)] <- 0
# Perform gap detection with distance ratio criterion
-gaps1 <- lidaRtRee::gapDetection(chm,
- ratio = 2, gap.max.height = 1, min.gap.surface = 50,
- nlFilter = "None"
+gaps1 <- lidaRtRee::gap_detection(chm,
+ ratio = 2, gap_max_height = 1, min_gap_surface = 50,
+ nl_filter = "None"
)
# Perform gap detection with distance ratio criterion and gap reconstruction
-gaps1r <- lidaRtRee::gapDetection(chm,
- ratio = 2, gap.max.height = 1, min.gap.surface = 50,
- gapReconstruct = TRUE, nlFilter = "None"
+gaps1r <- lidaRtRee::gap_detection(chm,
+ ratio = 2, gap_max_height = 1, min_gap_surface = 50,
+ gap_reconstruct = TRUE, nl_filter = "None"
)
# 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)
+raster::plot(gaps1$gap_id > 0, main = "Gaps", legend = FALSE)
+raster::plot(gaps1r$gap_id > 0, main = "Gaps extended by reconstruction", legend = FALSE)
```
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.
```{r edgeDetection, include = TRUE, out.width='70%', fig.width = 9, fig.height = 4.5}
# Perform edge detection by erosion
-edgeInside <- lidaRtRee::edgeDetection(gaps1r$gap.id > 0)
+edge_inside <- lidaRtRee::edge_detection(gaps1r$gap_id > 0)
# Perform edge detection by dilation
-edgeOutside <- lidaRtRee::edgeDetection(gaps1r$gap.id > 0, inside = FALSE)
+edge_outside <- lidaRtRee::edge_detection(gaps1r$gap_id > 0, inside = FALSE)
# display zoom on edges
par(mfrow = c(1, 2))
-raster::plot(edgeInside, main = "Edges (detection by erosion)", legend = FALSE)
-raster::plot(edgeOutside, main = "Edges (detection by dilation)", legend = FALSE)
+raster::plot(edge_inside, main = "Edges (detection by erosion)", legend = FALSE)
+raster::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(edgeInside))/(nrow(edgeInside)*ncol(edgeInside))*100, 1)`, whereas it is `r round(sum(raster::values(edgeOutside))/(nrow(edgeOutside)*ncol(edgeOutside))*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(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".
-## References
+## References
\ No newline at end of file
diff --git a/R/tree.detection.Rmd b/R/tree.detection.Rmd
index 174fce288895584e94ee8d3a36b0fc080c12c195..7d82c85b6d259752ecf59d990af2cdd02ec8ce7e 100644
--- a/R/tree.detection.Rmd
+++ b/R/tree.detection.Rmd
@@ -3,8 +3,8 @@ title: "R workflow for tree segmentation from ALS data"
author: "Jean-Matthieu Monnet"
date: "`r Sys.Date()`"
output:
- pdf_document: default
html_document: default
+ pdf_document: default
papersize: a4
bibliography: "../bib/bibliography.bib"
---
@@ -44,16 +44,16 @@ The field inventory corresponds to a 50m x 50m plot located in the Chablais moun
```{r loadTreeInventory, include = TRUE}
# load dataset from package (default)
-data(treeinventorychablais3, package = "lidaRtRee")
+data(tree_inventory_chablais3, package = "lidaRtRee")
```
Otherwise you can load your own data provided positions and heights are measured.
```{r prepareTreeInventory, eval=FALSE}
# import field inventory
fichier <- "chablais3_listeR.csv"
-tree.inventory <- read.csv(file = fichier, sep = ";", header = F)
-names(tree.inventory) <- c("x", "y", "d", "h", "n", "s", "e", "t")
+tree.inventory <- read.csv(file = fichier, sep = ";", header = F, stringsAsFactors = TRUE)
+names(tree_inventory_chablais3) <- c("x", "y", "d", "h", "n", "s", "e", "t")
# save as rda for later access
-# save(tree.inventory,file="tree.inventory.rda")
+# save(tree_inventory_chablais3,file="tree_inventory_chablais3.rda")
```
Attributes are:
@@ -66,26 +66,26 @@ names(tree.inventory) <- c("x", "y", "d", "h", "n", "s", "e", "t")
+ `e` appearance (0: missing or lying, 1: normal, 2: broken treetop, 3: dead with branches, 4: snag)
+ `t` tilted (0: no, 1: yes)
```{r displaytreeinventory, echo=FALSE}
-head(treeinventorychablais3, n = 3L)
+head(tree_inventory_chablais3, n = 3L)
```
-Function `plotTreeInventory` is designed to plot forest inventory data.
+Function `plot_tree_inventory` is designed to plot forest inventory data.
```{r plotTreeInventory, include = TRUE, out.width = '50%', fig.dim=c(6.5, 5.5)}
# display inventoried trees
-lidaRtRee::plotTreeInventory(
- treeinventorychablais3[, c("x", "y")],
- treeinventorychablais3$h,
- species = as.character(treeinventorychablais3$s)
+lidaRtRee::plot_tree_inventory(tree_inventory_chablais3
+ [, c("x", "y")],
+ tree_inventory_chablais3$h,
+ species = as.character(tree_inventory_chablais3$s)
)
```
The `ggplot2` package also provides nice outputs.
```{r ggplot2, include = TRUE, out.width = '50%', fig.dim=c(5.5, 5.5)}
# use table of species of package lidaRtRee to always use the same color for a given species
-plot.species <- lidaRtRee::speciesColor()[levels(treeinventorychablais3$s), "col"]
+plot.species <- lidaRtRee::species_color()[levels(tree_inventory_chablais3$s), "col"]
library(ggplot2)
-ggplot(treeinventorychablais3, aes(x = x, y = y, group = s)) +
+ggplot(tree_inventory_chablais3, aes(x = x, y = y, group = s)) +
geom_point(aes(color = s, size = d)) +
coord_sf(datum = 2154) +
scale_color_manual(values = plot.species) +
@@ -98,9 +98,9 @@ We define the region of interest (ROI) to crop ALS data to corresponding extent
```{r roi, include = TRUE}
# buffer to apply around ROI (meters)
-ROI.buffer <- 10
+ROI_buffer <- 10
# ROI limits
-ROI.range <- data.frame(round(apply(treeinventorychablais3[, c("x", "y")], 2, range)))
+ROI_range <- data.frame(round(apply(tree_inventory_chablais3[, c("x", "y")], 2, range)))
```
### Airborne Laser Scanning data
@@ -109,13 +109,13 @@ In this tutorial, ALS data available in the `lidaRtRee` package is used.
```{r loadALS, include = TRUE, message = FALSE}
# load data in package lidaRtRee (default)
-data(laschablais3, package = "lidaRtRee")
-laschablais3
+data(las_chablais3, package = "lidaRtRee")
+las_chablais3
```
Otherwise, ALS data is loaded with functions of package `lidR`. First a catalog of files containing ALS data is built. Then points located inside our ROI are loaded.
-```{r prepareALS, eval=FALSE, message = FALSE}
+```{r prepareALS, eval=FALSE, message = FALSE, warning = FALSE}
# directory for laz files
lazdir <- "../data/tree.detection"
# build catalog of files
@@ -124,15 +124,15 @@ cata <- lidR::readALSLAScatalog(lazdir)
# set coordinate system
lidR::projection(cata) <- 2154
# extract points in ROI plus additional buffer
-laschablais3 <- lidR::clip_rectangle(
+las_chablais3 <- lidR::clip_rectangle(
cata,
- ROI.range$x[1] - ROI.buffer - 5,
- ROI.range$y[1] - ROI.buffer - 5,
- ROI.range$x[2] + ROI.buffer + 5,
- ROI.range$y[2] + ROI.buffer + 5
+ ROI_range$x[1] - ROI_buffer - 5,
+ ROI_range$y[1] - ROI_buffer - 5,
+ ROI_range$x[2] + ROI_buffer + 5,
+ ROI_range$y[2] + ROI_buffer + 5
)
# save as rda for easier access:
-# save(laschablais3, file="laschablais3.rda", compress = "bzip2")
+# save(las_chablais3, file="las_chablais3.rda", compress = "bzip2")
```
## Data preparation
@@ -140,15 +140,15 @@ laschablais3 <- lidR::clip_rectangle(
From the ALS point cloud, digital elevation models are computed [@Monnet11c, pp. 43-46]. The Digital Terrain Model (DTM) represents the ground surface, it is computed by bilinear interpolation of points classified as ground. The Digital Surface Model (DSM) represents the canopy surface, it is computed by retaining in each raster's pixel the value of the highest point included in that pixel. The Canopy Height Model (CHM) is the normalized height of the canopy. It is computed by subtracting the DTM to the DSM.
-```{r computeDEMs, include=TRUE}
+```{r computeDEMs, include=TRUE, warning = FALSE}
# terrain model computed from points classified as ground
-dtm <- lidaRtRee::points2DTM(lidR::filter_ground(laschablais3),
- res = 0.5, ROI.range$x[1] - ROI.buffer,
- ROI.range$x[2] + ROI.buffer, ROI.range$y[1] - ROI.buffer,
- ROI.range$y[2] + ROI.buffer
+dtm <- lidaRtRee::points2DTM(lidR::filter_ground(las_chablais3),
+ res = 0.5, ROI_range$x[1] - ROI_buffer,
+ ROI_range$x[2] + ROI_buffer, ROI_range$y[1] - ROI_buffer,
+ ROI_range$y[2] + ROI_buffer
)
# surface model
-dsm <- lidaRtRee::points2DSM(laschablais3,
+dsm <- lidaRtRee::points2DSM(las_chablais3,
res = 0.5, dtm@extent@xmin,
dtm@extent@xmax, dtm@extent@ymin, dtm@extent@ymax
)
@@ -172,19 +172,19 @@ A plot mask is computed from the inventoried positions, using a height-dependent
```{r computeMask, include=TRUE, message=FALSE}
# select trees with height measures
-selec <- which(!is.na(treeinventorychablais3$h))
+selec <- which(!is.na(tree_inventory_chablais3$h))
# plot mask computation based on inventoried positions
# convex hull of plot
-ChullMask <- lidaRtRee::rasterChullMask(treeinventorychablais3[selec, c("x", "y")], dsm)
+mask_chull <- lidaRtRee::raster_chull_mask(tree_inventory_chablais3[selec, c("x", "y")], dsm)
# union of buffers around trees
-TreeMask <- lidaRtRee::rasterXYMask(
- treeinventorychablais3[selec, c("x", "y")],
- 2.1 + 0.14 * treeinventorychablais3$h[selec], dsm
+mask_tree <- lidaRtRee::raster_xy_mask(
+ tree_inventory_chablais3[selec, c("x", "y")],
+ 2.1 + 0.14 * tree_inventory_chablais3$h[selec], dsm
)
# union of convexHull and tree buffers
-plotMask <- max(ChullMask, TreeMask)
+mask_plot <- max(mask_chull, mask_tree)
# vectorize plot mask
-v.plotMask <- raster::rasterToPolygons(plotMask, function(x) (x == 1), dissolve = T)
+mask_plot_v <- raster::rasterToPolygons(mask_plot, function(x) (x == 1), dissolve = T)
```
Displaying inventoried trees on the CHM shows a pretty good correspondance of crowns visible in the CHM with trunk locations and sizes.
@@ -195,17 +195,17 @@ raster::plot(chm,
main = "Canopy Height Model and tree positions"
)
# add inventoried trees
-lidaRtRee::plotTreeInventory(treeinventorychablais3[, c("x", "y")],
- treeinventorychablais3$h,
- species = as.character(treeinventorychablais3$s), add = TRUE
+lidaRtRee::plot_tree_inventory(tree_inventory_chablais3[, c("x", "y")],
+ tree_inventory_chablais3$h,
+ species = as.character(tree_inventory_chablais3$s), add = TRUE
)
# display plot mask
-raster::plot(v.plotMask, border = "red", add = TRUE)
+raster::plot(mask_plot_v, border = "red", add = TRUE)
```
## Tree delineation
### Tree segmentation
-Tree segmentation is performed on the Canopy Height Model by using a general function which consists in the following steps:
+Tree segmentation is performed on the Canopy Height Model by using a general function (`tree_segmentation`) which consists in the following steps:
* image pre-processing (non-linear filtering and smoothing for noise removal)
+ local maxima filtering and selection for tree top detection
@@ -214,27 +214,27 @@ Tree segmentation is performed on the Canopy Height Model by using a general fun
The first two steps are documented in @Monnet11c, pp. 47-52.
```{r plotSegm, include=TRUE, fig.width = 12, fig.height = 4.5, out.width='100%'}
# tree detection (default settings), applied on canopy height model
-segms <- lidaRtRee::treeSegmentation(chm)
+segms <- lidaRtRee::tree_segmentation(chm)
#
par(mfrow = c(1, 3))
# display pre-processed chm
-raster::plot(segms$smoothed.dem, main = "Pre-processed CHM")
+raster::plot(segms$smoothed_dem, main = "Pre-processed CHM")
# display selected local maxima
-raster::plot(segms$local.maxima, main = "Selected local maxima")
+raster::plot(segms$local_maxima, main = "Selected local maxima")
# display segments, except ground segment
-dummy <- segms$segments.id
+dummy <- segms$segments_id
dummy[dummy == 0] <- NA
raster::plot(dummy %% 8, main = "Segments (random colors)", col = rainbow(8), legend = FALSE)
```
### Tree extraction
-A data.frame of detected trees located in the plot mask is then extracted with the function `treeExtraction.` Segments can be converted from raster to polygons but the operation is quite slow. Attributes are :
+A data.frame of detected trees located in the plot mask is then extracted with the function `tree_extraction.` Segments can be converted from raster to polygons but the operation is quite slow. Attributes are :
* `id`: tree id
+ `x`: easting coordinate of tree top
+ `y`: northing coordinate of tree top
+ `h`: height of tree top
-+ `dom.radius`: distance of tree top to nearest crown of neighbouring tree with larger height
++ `dom_radius`: distance of tree top to nearest crown of neighbouring tree with larger height
+ `s`: crown surface
+ `v`: crown volume
+ `sp`: crown surface inside plot
@@ -243,21 +243,21 @@ A data.frame of detected trees located in the plot mask is then extracted with t
```{r plotSegmTrees, include=TRUE, out.width = '50%', fig.dim=c(4.5, 4.5)}
# tree extraction only inside plot mask for subsequent comparison
-trees <- lidaRtRee::treeExtraction(
- segms$filled.dem,
- segms$local.maxima,
- segms$segments.id, plotMask
+trees <- lidaRtRee::tree_extraction(
+ segms$filled_dem,
+ segms$local_maxima,
+ segms$segments_id, mask_plot
)
head(trees, n = 3L)
# convert segments from raster to polygons
-v.segments <- raster::rasterToPolygons(segms[[2]], dissolve = T)
+segments_v <- raster::rasterToPolygons(segms$segments_id, dissolve = T)
#
# display initial image
raster::plot(chm, col = gray(seq(0, 1, 1 / 255)), main = "CHM and detected positions")
# display segments border
-sp::plot(v.segments, border = "white", add = T)
+sp::plot(segments_v, border = "white", add = T)
# display plot mask
-sp::plot(v.plotMask, border = "red", add = T)
+sp::plot(mask_plot_v, border = "red", add = T)
# display inventoried trees
graphics::points(trees$x, trees$y, col = "blue", cex = trees$h / 20, pch = 2)
```
@@ -269,14 +269,14 @@ To assess detection accuracy, reference (field) trees should be linked to detect
```{r plotMached, include=TRUE, out.width = '70%', fig.dim=c(6.5, 4.5)}
# match detected treetops with field trees based on relative distance of apices
-matched <- lidaRtRee::treeMatching(
- treeinventorychablais3[selec, c("x", "y", "h")],
+matched <- lidaRtRee::tree_matching(
+ tree_inventory_chablais3[selec, c("x", "y", "h")],
cbind(trees@coords, trees$h)
)
# display matching results
-lidaRtRee::plot2Dmatched(
- treeinventorychablais3[selec, c("x", "y", "h")],
- cbind(trees@coords, trees$h), matched, chm, v.plotMask
+lidaRtRee::plot_matched(
+ tree_inventory_chablais3[selec, c("x", "y", "h")],
+ cbind(trees@coords, trees$h), matched, chm, mask_plot_v
)
```
@@ -284,18 +284,18 @@ lidaRtRee::plot2Dmatched(
```{r detectionStats, include=FALSE}
# height histogram of detections
-detection.stats <- lidaRtRee::histDetection(
- treeinventorychablais3[selec, c("x", "y", "h")],
+detection_stats <- lidaRtRee::hist_detection(
+ tree_inventory_chablais3[selec, c("x", "y", "h")],
cbind(trees@coords, trees$h), matched
)
```
-Detection accuracy is evaluated thanks to the number of correct detections (`r detection.stats$true.detections`), false detections (`r detection.stats$false.detections`) and omissions (`r detection.stats$omissions`). In heterogeneous stands, detection accuracy is height-dependent, it is informative to display those categories on a height histogram.
+Detection accuracy is evaluated thanks to the number of correct detections (`r detection_stats$true_detections`), false detections (`r detection_stats$false_detections`) and omissions (`r detection_stats$omissions`). In heterogeneous stands, detection accuracy is height-dependent, it is informative to display those categories on a height histogram.
```{r plotDetection, include=TRUE, out.width = '40%', fig.dim=c(4.5, 4.5)}
# height histogram of detections
-detection.stats <- lidaRtRee::histDetection(
- treeinventorychablais3[selec, c("x", "y", "h")],
+detection_stats <- lidaRtRee::hist_detection(
+ tree_inventory_chablais3[selec, c("x", "y", "h")],
cbind(trees@coords, trees$h), matched
)
```
@@ -303,44 +303,44 @@ detection.stats <- lidaRtRee::histDetection(
### Height estimation accuracy
```{r heighRegression, include=FALSE}
-heightReg <- lidaRtRee::heightRegression(treeinventorychablais3[selec, c("x", "y", "h")],
+height_reg <- lidaRtRee::height_regression(tree_inventory_chablais3[selec, c("x", "y", "h")],
cbind(trees@coords, trees$h),
matched,
- species = treeinventorychablais3$s
+ species = tree_inventory_chablais3$s
)
```
-For true detections, estimated heights can be compared to field-measured heights. Here, the root mean square error is `r round(heightReg$stats$rmse,1)`m, while the bias is `r round(heightReg$stats$bias,1)`m. The linear regression is displayed hereafter.
+For true detections, estimated heights can be compared to field-measured heights. Here, the root mean square error is `r round(height_reg$stats$rmse,1)`m, while the bias is `r round(height_reg$stats$bias,1)`m. The linear regression is displayed hereafter.
```{r plotRegression, include=TRUE, out.width = '40%', fig.dim=c(4.5, 4.5)}
# linear regression between reference height and estimated height
-heightReg <- lidaRtRee::heightRegression(treeinventorychablais3[selec, c("x", "y", "h")],
+height_reg <- lidaRtRee::height_regression(tree_inventory_chablais3[selec, c("x", "y", "h")],
cbind(trees@coords, trees$h),
matched,
- species = treeinventorychablais3$s
+ species = tree_inventory_chablais3$s
)
```
## Species Classification
+
### Points in segments
Before computation of point cloud metrics in each segment, the whole point cloud is normalized to convert point altitude to height above ground. Points are then labeled with the id of the segment they belong to. A list of LAS objects corresponding to the segments is then created.
```{r pointsSegments, include=TRUE, warning=FALSE, message=FALSE}
# normalize point cloud
-lasn <- lidR::normalize_height(laschablais3, lidR::tin())
+lasn <- lidR::normalize_height(las_chablais3, lidR::tin())
# add segment id in LAS object
-lasn <- lidR::add_attribute(lasn, raster::extract(segms[["segments.id"]], lasn@data[, 1:2]), "seg.id")
+lasn <- lidR::add_attribute(lasn, raster::extract(segms$segments_id, lasn@data[, 1:2]), "seg_id")
# split las object by segment id
-lasl <- split(lasn@data, lasn@data$seg.id)
+las_l <- split(lasn@data, lasn@data$seg_id)
# convert list of data.frames to list of las objects
-lasl <- lapply(lasl, function(x) {
+las_l <- lapply(las_l, function(x) {
lidR::LAS(x, lasn@header)
})
# set coordinate system
-dummy <- sf::st_crs(2154)
-for (i in 1:length(lasl)) {
- lidR::projection(lasl[[i]]) <- dummy
+for (i in 1:length(las_l)) {
+ lidR::projection(las_l[[i]]) <- 2154
}
```
@@ -349,13 +349,13 @@ for (i in 1:length(lasl)) {
Basic point cloud metrics are computed in each segment (maximum, mean, standard deviation of heights, mean and standard deviation of intensity).
```{r metrics, include=TRUE, eval=TRUE}
# compute basic las metrics in each segment
-metrics <- lidaRtRee::cloudMetrics(lasl, func = ~ list(
+metrics <- lidaRtRee::clouds_metrics(las_l, func = ~ list(
maxZ = max(Z), meanZ = mean(Z),
sdZ = sd(Z), meanI = mean(Intensity),
sdI = sd(Intensity)
))
# create segment id attribute
-metrics$seg.id <- row.names(metrics)
+metrics$seg_id <- row.names(metrics)
head(metrics, n = 3L)
```
@@ -364,33 +364,33 @@ head(metrics, n = 3L)
Computed metrics are merged with reference and detected trees, thanks to the segment id.
```{r mergeMetrics, include=TRUE, eval=TRUE}
# associate each reference tree with the segment that contains its trunk.
-treeinventorychablais3$seg.id <- raster::extract(
- segms[["segments.id"]],
- treeinventorychablais3[, c("x", "y")]
+tree_inventory_chablais3$seg_id <- raster::extract(
+ segms$segments_id,
+ tree_inventory_chablais3[, c("x", "y")]
)
# create new data.frame by merging metrics and inventoried trees based on segment id
-tree.metrics <- base::merge(treeinventorychablais3, metrics)
+tree_metrics <- base::merge(tree_inventory_chablais3, metrics)
# remove non-tree segment
-tree.metrics <- tree.metrics[tree.metrics$seg.id != 0, ]
+tree_metrics <- tree_metrics[tree_metrics$seg_id != 0, ]
# add metrics to extracted tree data.frame
-trees <- base::merge(trees, metrics, by.x = "id", by.y = "seg.id", all.x = T)
+trees <- base::merge(trees, metrics, by.x = "id", by.y = "seg_id", all.x = T)
```
Plotting the reference trees with the mean intensity of lidar points in the segments shows that when they are dominant, broadleaf trees tend to have higher mean intensity than coniferous trees.
.
```{r Intensity, include=TRUE, out.width = '70%', fig.dim=c(6.5, 4.5)}
# create raster of segment' mean intensity
-r.mean.intensity.segm <- segms[["segments.id"]]
+r_mean_intensity_segm <- segms$segments_id
# match segment id with id in metrics data.frame
-dummy <- match(raster::values(r.mean.intensity.segm), trees$id)
+dummy <- match(raster::values(r_mean_intensity_segm), trees$id)
# replace segment id with corresponding mean intensity
-raster::values(r.mean.intensity.segm) <- trees$meanI[dummy]
+raster::values(r_mean_intensity_segm) <- trees$meanI[dummy]
# display tree inventory with mean intensity in segment background
-raster::plot(r.mean.intensity.segm, main = "Mean intensity in segment")
+raster::plot(r_mean_intensity_segm, main = "Mean intensity in segment")
# display tree inventory
-lidaRtRee::plotTreeInventory(treeinventorychablais3[, c("x", "y")],
- treeinventorychablais3$h,
- species = as.character(treeinventorychablais3$s), add = T
+lidaRtRee::plot_tree_inventory(tree_inventory_chablais3[, c("x", "y")],
+ tree_inventory_chablais3$h,
+ species = as.character(tree_inventory_chablais3$s), add = T
)
```
@@ -400,28 +400,28 @@ A boxplot of mean intensity in segments per species shows that mean intensity di
```{r metricsByHighestTree, include=FALSE, eval=TRUE}
# create new variable orderer by segment id then decreasing height
-tree.metrics.h <- tree.metrics[order(tree.metrics$seg.id, -tree.metrics$h), ]
+tree_metrics_h <- tree_metrics[order(tree_metrics$seg_id, -tree_metrics$h), ]
i <- 2
# leave only first (highest) tree per segment id
-while (i < nrow(tree.metrics.h)) {
- if (tree.metrics.h$seg.id[i] == tree.metrics.h$seg.id[i + 1]) {
- tree.metrics.h <- tree.metrics.h[-(i + 1), ]
+while (i < nrow(tree_metrics_h)) {
+ if (tree_metrics_h$seg_id[i] == tree_metrics_h$seg_id[i + 1]) {
+ tree_metrics_h <- tree_metrics_h[-(i + 1), ]
} else {
i <- i + 1
}
}
-tree.metrics.h$s <- factor(tree.metrics.h$s)
+tree_metrics_h$s <- factor(tree_metrics_h$s)
```
```{r boxplot2, include=TRUE, out.width = '80%', fig.dim=c(9.5, 5.5)}
par(mfrow = c(1, 2))
boxplot(meanI ~ s,
- data = tree.metrics[, c("s", "maxZ", "meanZ", "sdZ", "meanI", "sdI")],
+ data = tree_metrics[, c("s", "maxZ", "meanZ", "sdZ", "meanI", "sdI")],
ylab = "Mean intensity in segment", xlab = "Specie",
main = "All inventoried trees", las = 2, varwidth = TRUE
)
boxplot(meanI ~ s,
- data = tree.metrics.h, ylab = "Mean intensity in segment",
+ data = tree_metrics_h, ylab = "Mean intensity in segment",
xlab = "Specie", main = "Highest inventoried tree in segment", las = 2,
varwidth = TRUE
)
@@ -430,8 +430,8 @@ boxplot(meanI ~ s,
A linear discriminant analysis shows that main species can be discriminated, thanks to a combination of height and intensity variables.
```{r AFD, include=TRUE, out.width = '60%', fig.dim=c(6, 6)}
-acp2 <- ade4::dudi.pca(tree.metrics.h[, c("maxZ", "meanZ", "sdZ", "meanI", "sdI")], scannf = F)
-afd <- ade4::discrimin(acp2, tree.metrics.h$s, scannf = F, nf = 2)
+acp2 <- ade4::dudi.pca(tree_metrics_h[, c("maxZ", "meanZ", "sdZ", "meanI", "sdI")], scannf = F)
+afd <- ade4::discrimin(acp2, tree_metrics_h$s, scannf = F, nf = 2)
plot(afd)
```
@@ -442,28 +442,28 @@ 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@data[which(lasn@data$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[, c("X", "Y", "Z")], col = points.seg$seg_id %% 10 + 1, aspect = FALSE)
#
# add inventoried trees
-treeinventorychablais3$z <- 0
-for (i in 1:nrow(treeinventorychablais3))
+tree_inventory_chablais3$z <- 0
+for (i in 1:nrow(tree_inventory_chablais3))
{
rgl::lines3d(
rbind(
- treeinventorychablais3$x[i] - 974300,
- treeinventorychablais3$x[i] - 974300
+ tree_inventory_chablais3$x[i] - 974300,
+ tree_inventory_chablais3$x[i] - 974300
),
rbind(
- treeinventorychablais3$y[i] - 6581600,
- treeinventorychablais3$y[i] - 6581600
+ tree_inventory_chablais3$y[i] - 6581600,
+ tree_inventory_chablais3$y[i] - 6581600
),
rbind(
- treeinventorychablais3$z[i],
- treeinventorychablais3$z[i] + treeinventorychablais3$h[i]
+ tree_inventory_chablais3$z[i],
+ tree_inventory_chablais3$z[i] + tree_inventory_chablais3$h[i]
)
)
}
@@ -480,35 +480,35 @@ rgl::rgl.bg(color = "white")
for (j in 1:length(selec))
{
i <- selec[j]
- if (is.na(treeinventorychablais3$h[i]) |
- is.na(treeinventorychablais3$d[i])) {
+ if (is.na(tree_inventory_chablais3$h[i]) |
+ is.na(tree_inventory_chablais3$d[i])) {
next
}
if (!is.element(
- as.character(treeinventorychablais3$s[i]),
+ as.character(tree_inventory_chablais3$s[i]),
c("ABAL", "PIAB", "TABA")
)) {
rLiDAR::LiDARForestStand(
crownshape = "halfellipsoid",
- CL = 0.6 * treeinventorychablais3$h[i],
- CW = treeinventorychablais3$h[i] / 4,
- HCB = 0.4 * treeinventorychablais3$h[i],
- dbh = treeinventorychablais3$d[i] / 50,
+ CL = 0.6 * tree_inventory_chablais3$h[i],
+ CW = tree_inventory_chablais3$h[i] / 4,
+ HCB = 0.4 * tree_inventory_chablais3$h[i],
+ dbh = tree_inventory_chablais3$d[i] / 50,
resolution = "high",
- X = treeinventorychablais3$x[i],
- Y = treeinventorychablais3$y[i],
+ X = tree_inventory_chablais3$x[i],
+ Y = tree_inventory_chablais3$y[i],
mesh = F
)
} else {
rLiDAR::LiDARForestStand(
crownshape = "cone",
- CL = 0.5 * treeinventorychablais3$h[i],
- CW = treeinventorychablais3$h[i] / 4,
- HCB = 0.5 * treeinventorychablais3$h[i],
- dbh = treeinventorychablais3$d[i] / 50,
+ CL = 0.5 * tree_inventory_chablais3$h[i],
+ CW = tree_inventory_chablais3$h[i] / 4,
+ HCB = 0.5 * tree_inventory_chablais3$h[i],
+ dbh = tree_inventory_chablais3$d[i] / 50,
resolution = "high",
- X = treeinventorychablais3$x[i],
- Y = treeinventorychablais3$y[i],
+ X = tree_inventory_chablais3$x[i],
+ Y = tree_inventory_chablais3$y[i],
mesh = F,
stemcolor = "burlywood4",
crowncolor = "darkgreen"
@@ -582,7 +582,7 @@ The steps for processing a batch of las/laz files are :
Be aware that in case tree segments are vectorized, some obtained polygons might overlap. The segmentation algorithm might not be deterministic and borders are sometimes not consistent when adjacent polygons are pasted from different tiles.
-```{r batch, include=TRUE, out.width = '100%', fig.dim=c(8, 8)}
+```{r batch, include=TRUE, out.width = '100%', fig.dim=c(8, 8), warning = FALSE}
rm(list = ls())
# BUILD CATALOG OF FILES
# folder containing the files
@@ -595,11 +595,11 @@ lidR::projection(cata) <- 2154
# BATCH PROCESSING PARAMETERS
# tile size to split area into chunks to process
# trade-off between RAM capacity VS total number of tiles to process
-tile.size <- 70 # here 70 for example purpose with small area
+tile_size <- 70 # here 70 for example purpose with small area
# buffer size: around each tile to avoid border effects on segmentation results
# trade-off between making sure a whole tree crown is processed in case its top
# is on the border VS duplicate processing
-buffer.size <- 10 # 5 m is minimum, 10 is probably better depending on tree size
+buffer_size <- 10 # 5 m is minimum, 10 is probably better depending on tree size
#
# TREE SEGMENTATION PARAMETERS
# set raster resolution
@@ -607,11 +607,11 @@ resolution <- 1
#
# OUTPUT PARAMETERS
# option to vectorize tree crowns (set to FALSE if too slow)
-vectorize.trees <- TRUE
+out_vectorize_trees <- TRUE
# output canopy height models ? (set to FALSE if too much RAM used)
-output.chm <- TRUE
+out_chm <- TRUE
# save chms on disk
-save.chm <- FALSE
+out_save_chm <- FALSE
#
# CLUSTER PARAMETERS
# working directory
@@ -624,20 +624,20 @@ parallel::clusterExport(cl = clust, ls(), envir = .GlobalEnv)
# COORDINATES OF TILES
# low left corner
x <- seq(
- from = floor(cata@bbox["x", "min"] / tile.size) * tile.size,
- by = tile.size,
- to = ceiling(cata@bbox["x", "max"] / tile.size - 1) * tile.size
+ from = floor(cata@bbox["x", "min"] / tile_size) * tile_size,
+ by = tile_size,
+ to = ceiling(cata@bbox["x", "max"] / tile_size - 1) * tile_size
) # [1:2]
-length.x <- length(x)
+length_x <- length(x)
y <- seq(
- from = floor(cata@bbox["y", "min"] / tile.size) * tile.size,
- by = tile.size,
- to = ceiling(cata@bbox["y", "max"] / tile.size - 1) * tile.size
+ from = floor(cata@bbox["y", "min"] / tile_size) * tile_size,
+ by = tile_size,
+ to = ceiling(cata@bbox["y", "max"] / tile_size - 1) * tile_size
) # [1:2]
-length.y <- length(y)
+length_y <- length(y)
# repeat coordinates for tiling while area
-x <- as.list(rep(x, length.y))
-y <- as.list(rep(y, each = length.x))
+x <- as.list(rep(x, length_y))
+y <- as.list(rep(y, each = length_x))
#
# PARALLEL PROCESSING
# function takes coordinates of tile as arguments
@@ -650,36 +650,36 @@ resultats <- parallel::mcmapply(
output <- list()
output$name <- paste0(i, "_", j)
# extract tile plus buffer
- las.tile <- lidR::clip_rectangle(
+ las_tile <- lidR::clip_rectangle(
cata,
- i - buffer.size,
- j - buffer.size,
- i + tile.size + buffer.size,
- j + tile.size + buffer.size
+ i - buffer_size,
+ j - buffer_size,
+ i + tile_size + buffer_size,
+ j + tile_size + buffer_size
)
# check if points are present
- if (nrow(las.tile@data) > 0) {
+ if (nrow(las_tile@data) > 0) {
# normalization if required
- # las.tile <- lidR::normalize_height(las.tile, lidR::tin())
+ # las_tile <- lidR::normalize_height(las_tile, lidR::tin())
# in this example LAS tiles are already normalized
# compute canopy height model
- chm <- lidaRtRee::points2DSM(las.tile)
+ chm <- lidaRtRee::points2DSM(las_tile)
# check all chm is not NA
if (!all(is.na(raster::values(chm)))) {
# output chm if asked for
- if (output.chm) {
+ if (out_chm) {
output$chm <- raster::crop(chm, raster::extent(
- i, i + tile.size,
- j, j + tile.size
+ i, i + tile_size,
+ j, j + tile_size
))
}
# save on disk
- if (save.chm) {
+ if (out_save_chm) {
raster::writeRaster(raster::crop(
chm,
raster::extent(
- i, i + tile.size,
- j, j + tile.size
+ i, i + tile_size,
+ j, j + tile_size
)
),
file = paste0("chm_", i, "_", j, ".tif"),
@@ -688,31 +688,31 @@ resultats <- parallel::mcmapply(
}
#
# tree detection
- segms <- lidaRtRee::treeSegmentation(chm)
+ segms <- lidaRtRee::tree_segmentation(chm)
# tree extraction
- trees <- lidaRtRee::treeExtraction(
- segms$filled.dem,
- segms$local.maxima,
- segms$segments.id
+ trees <- lidaRtRee::tree_extraction(
+ segms$filled_dem,
+ segms$local_maxima,
+ segms$segments_id
)
# remove trees in buffer area
- trees <- trees[trees$x >= i & trees$x < i + tile.size &
- trees$y >= j & trees$y < j + tile.size, ]
+ trees <- trees[trees$x >= i & trees$x < i + tile_size &
+ trees$y >= j & trees$y < j + tile_size, ]
# add tile id to trees to avoid duplicates in final file
trees$tile <- paste0(i, "_", j)
# convert to vectors if option is TRUE
- if (vectorize.trees) {
+ if (out_vectorize_trees) {
# vectorize
- v.trees <- raster::rasterToPolygons(segms$segments.id, dissolve = T)
+ trees_v <- raster::rasterToPolygons(segms$segments_id, dissolve = T)
# remove polygons which treetop is in buffer
- v.trees <- v.trees[is.element(v.trees$segments.id, trees$id), ]
- # names(v.trees) <- "id"
+ trees_v <- trees_v[is.element(trees_v$segments_id, trees$id), ]
+ # names(trees_v) <- "id"
# add attributes
# errors when using sp::merge so using sp::over even if it is probably slower
- # merge(v.trees@data, trees, all.x = TRUE)
- v.trees@data <- cbind(v.trees@data, sp::over(v.trees, trees))
+ # merge(trees_v@data, trees, all.x = TRUE)
+ trees_v@data <- cbind(trees_v@data, sp::over(trees_v, trees))
# save in list
- output$v.trees <- v.trees
+ output$trees_v <- trees_v
}
# save trees in list
output$trees <- trees
@@ -735,20 +735,20 @@ trees[sapply(trees, is.null)] <- NULL
trees <- do.call(rbind, trees)
#
# chm
-if (output.chm) {
+if (out_chm) {
chm <- lapply(resultats, function(x) x[["chm"]])
# merge chm
# no names in list otherwise do.call returns an error
- chm.all <- do.call(raster::merge, chm)
+ chm_all <- do.call(raster::merge, chm)
names(chm) <- id
}
-# v.trees
-if (vectorize.trees) {
- v.trees <- lapply(resultats, function(x) x[["v.trees"]])
+# trees_v
+if (out_vectorize_trees) {
+ trees_v <- lapply(resultats, function(x) x[["trees_v"]])
# remove NULL elements
- v.trees[sapply(v.trees, is.null)] <- NULL
- v.trees <- do.call(rbind, v.trees)
- # 1-pixel overlapping in v.trees might be present because image segmentation
+ trees_v[sapply(trees_v, is.null)] <- NULL
+ trees_v <- do.call(rbind, trees_v)
+ # 1-pixel overlapping in trees_v might be present because image segmentation
# is not fully identical in overlap areas of adjacent tiles.
}
```
@@ -757,14 +757,14 @@ The following image displays the results for the whole area.
```{r batch.plot, include=TRUE, out.width = '90%', fig.dim=c(8.5, 5.5)}
# threshold outsiders in chm
-chm.all[chm.all > 40] <- 40
-chm.all[chm.all < 0] <- 0
+chm_all[chm_all > 40] <- 40
+chm_all[chm_all < 0] <- 0
# display chm
-raster::plot(chm.all,
+raster::plot(chm_all,
main = "Canopy Height Model and segmented trees"
)
# display segments border
-sp::plot(v.trees, border = "white", add = T)
+sp::plot(trees_v, border = "white", add = T)
# add trees
sp::plot(trees, cex = trees$h / 40, add = TRUE, pch = 2)
```
@@ -773,11 +773,11 @@ The following lines save outputs to files.
```{r batch.export, eval = FALSE}
# merged chm
-raster::writeRaster(chm.all, file = "chm.tif", overwrite = TRUE)
+raster::writeRaster(chm_all, file = "chm.tif", overwrite = TRUE)
# trees
sf::st_write(sf::st_as_sf(trees), 'trees_points.gpkg', 'trees_points', delete_dsn = T)
# vectorized trees
-if (vectorize.trees) sf::st_write(sf::st_as_sf(v.trees), 'v_trees_points.gpkg', 'v_trees_points', delete_dsn = T)
+if (out_vectorize_trees) sf::st_write(sf::st_as_sf(trees_v), 'v_trees_points.gpkg', 'v_trees_points', delete_dsn = T)
```
-## References
+## References
\ No newline at end of file