Commit 6a271dc4 authored by Monnet Jean-Matthieu's avatar Monnet Jean-Matthieu
Browse files

work in progress ABA3

No related merge requests found
Showing with 266 additions and 92 deletions
+266 -92
...@@ -242,7 +242,7 @@ Terrain statistics can be extracted from the cloud of ground points with altitud ...@@ -242,7 +242,7 @@ Terrain statistics can be extracted from the cloud of ground points with altitud
```{r terrainStats, include=TRUE, warning=FALSE, message=FALSE, fig.width = 4, fig.height = 6} ```{r terrainStats, include=TRUE, warning=FALSE, message=FALSE, fig.width = 4, fig.height = 6}
# use mapply to apply points2terrainStats function to each point cloud # use mapply to apply points2terrainStats function to each point cloud
# while providing the coordinates of each center # while providing the coordinates of each center
terrain.stats <- mapply(function(x, y) metrics.terrain <- mapply(function(x, y)
{ {
lidaRtRee::points2terrainStats(x, plots[y, c("X", "Y")], p.radius) lidaRtRee::points2terrainStats(x, plots[y, c("X", "Y")], p.radius)
}, },
...@@ -250,13 +250,13 @@ x = llas.ground, ...@@ -250,13 +250,13 @@ x = llas.ground,
y = as.list(1:length(llas.ground)), y = as.list(1:length(llas.ground)),
SIMPLIFY = FALSE) SIMPLIFY = FALSE)
# bind results in data.frame # bind results in data.frame
terrain.stats <- do.call(rbind, terrain.stats) metrics.terrain <- do.call(rbind, metrics.terrain)
# #
# compute terrain stats without specifying centre # compute terrain stats without specifying centre
terrain.stats2 <- lapply(llas.ground, lidaRtRee::points2terrainStats) metrics.terrain2 <- lapply(llas.ground, lidaRtRee::points2terrainStats)
terrain.stats2 <- do.call(rbind, terrain.stats2) metrics.terrain2 <- do.call(rbind, metrics.terrain2)
# display results for comparison # display results for comparison
head(cbind(terrain.stats, terrain.stats2), n=3) head(cbind(metrics.terrain, metrics.terrain2), n=3)
``` ```
## Check field inventory data ## Check field inventory data
...@@ -529,7 +529,7 @@ The following lines save the data required for the [area-based model calibration ...@@ -529,7 +529,7 @@ The following lines save the data required for the [area-based model calibration
```{r save.plots, include=TRUE} ```{r save.plots, include=TRUE}
# save(plots, file="../data/aba.model/output/plots.rda") # save(plots, file="../data/aba.model/output/plots.rda")
# save(llas.height, file="../data/aba.model/output/llas.height.rda", compress="bzip2") # save(llas.height, file="../data/aba.model/output/llas.height.rda", compress="bzip2")
# save(terrain.stats, file="../data/aba.model/output/terrain.stats.rda") # save(metrics.terrain, file="../data/aba.model/output/metrics.terrain.rda")
``` ```
...@@ -61,15 +61,15 @@ The first lines of the terrain statistics are displayed hereafter. ...@@ -61,15 +61,15 @@ The first lines of the terrain statistics are displayed hereafter.
```{r loadTerrainStats, echo = FALSE} ```{r loadTerrainStats, echo = FALSE}
# terrain statistics previously computed with (non-normalized) ground points inside each plot extent # terrain statistics previously computed with (non-normalized) ground points inside each plot extent
load("../data/aba.model/output/terrain.stats.rda") load("../data/aba.model/output/metrics.terrain.rda")
head(terrain.stats[, 1:3], n=3) 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. The following lines ensure that the plots are ordered in the same way in the three data objects.
```{r orderData, include=TRUE} ```{r orderData, include=TRUE}
llas.height <- llas.height[plots$plotId] llas.height <- llas.height[plots$plotId]
terrain.stats <- terrain.stats[plots$plotId,] metrics.terrain <- metrics.terrain[plots$plotId,]
``` ```
# ALS metrics computation # ALS metrics computation
...@@ -87,8 +87,8 @@ Point cloud metrics are computed with the function `lidaRtRee::cloudMetrics`, wh ...@@ -87,8 +87,8 @@ Point cloud metrics are computed with the function `lidaRtRee::cloudMetrics`, wh
# define function for later use # define function for later use
aba.pointMetricsFUN <- ~lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2) aba.pointMetricsFUN <- ~lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2)
# apply function on each point cloud in list # apply function on each point cloud in list
point.metrics <- lidaRtRee::cloudMetrics(llas.height, aba.pointMetricsFUN) metrics.points <- lidaRtRee::cloudMetrics(llas.height, aba.pointMetricsFUN)
round(head(point.metrics[, 1:8], n = 3),2) round(head(metrics.points[, 1:8], n = 3),2)
``` ```
## Tree metrics ## Tree metrics
...@@ -101,14 +101,14 @@ aba.resCHM <- 0.5 ...@@ -101,14 +101,14 @@ aba.resCHM <- 0.5
# specify plot radius to exclude trees located outside plots # specify plot radius to exclude trees located outside plots
plot.radius <- 15 plot.radius <- 15
# compute tree metrics # compute tree metrics
tree.metrics <- lidaRtRee::cloudTreeMetrics(llas.height, plots[, c("X", "Y")], metrics.tree <- lidaRtRee::cloudTreeMetrics(llas.height, plots[, c("X", "Y")],
plot.radius, res = aba.resCHM, plot.radius, res = aba.resCHM,
func=function(x) func=function(x)
{ {
lidaRtRee::stdTreeMetrics(x, lidaRtRee::stdTreeMetrics(x,
area.ha=pi*plot.radius^2/10000) area.ha=pi*plot.radius^2/10000)
}) })
round(head(tree.metrics[, 1:5], n = 3), 2) round(head(metrics.tree[, 1:5], n = 3), 2)
``` ```
## Other metrics ## Other metrics
...@@ -116,9 +116,9 @@ round(head(tree.metrics[, 1:5], n = 3), 2) ...@@ -116,9 +116,9 @@ round(head(tree.metrics[, 1:5], n = 3), 2)
In case terrain metrics have been computed from the cloud of ground points only, they can also be added as variables, and so do other environmental variables which might be relevant in modeling. In case terrain metrics have been computed from the cloud of ground points only, they can also be added as variables, and so do other environmental variables which might be relevant in modeling.
```{r bindMetrics, include=TRUE} ```{r bindMetrics, include=TRUE}
metrics <- cbind(point.metrics[plots$plotId, ], metrics <- cbind(metrics.points[plots$plotId, ],
tree.metrics[plots$plotId, ], metrics.tree[plots$plotId, ],
terrain.stats[plots$plotId, 1:3]) metrics.terrain[plots$plotId, 1:3])
``` ```
# Model calibration # Model calibration
...@@ -147,7 +147,7 @@ In this example, only tree metrics are selected in the basal area prediction mod ...@@ -147,7 +147,7 @@ In this example, only tree metrics are selected in the basal area prediction mod
```{r modelPlot, include=TRUE, fig.height = 4.5, fig.width = 8} ```{r modelPlot, include=TRUE, fig.height = 4.5, fig.width = 8}
# check correlation between errors and other variables # 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")], terrain.stats[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 # 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 # plot predicted VS field values
...@@ -158,20 +158,20 @@ abline(h = 0, lty = 2) ...@@ -158,20 +158,20 @@ 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. 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 point.metricsOnly, include=TRUE, fig.height = 4.5, fig.width = 8} ```{r metrics.pointsOnly, include=TRUE, fig.height = 4.5, fig.width = 8}
model.ABA.point.metrics <- lidaRtRee::ABAmodel(plots[subsample,variable], point.metrics[subsample,], transform="boxcox", nmax=4, xy = plots[subsample, c("X", "Y")]) model.ABA.metrics.points <- lidaRtRee::ABAmodel(plots[subsample,variable], metrics.points[subsample,], transform="boxcox", nmax=4, xy = plots[subsample, c("X", "Y")])
# renames outputs # renames outputs
row.names(model.ABA.point.metrics$stats) <- names(model.ABA.point.metrics$model) <- variable row.names(model.ABA.metrics.points$stats) <- names(model.ABA.metrics.points$model) <- variable
# model.ABA.point.metrics$model[[variable]] # model.ABA.metrics.points$model[[variable]]
model.ABA.point.metrics$stats model.ABA.metrics.points$stats
# cor.test(model.ABA.point.metrics$values$residual, plots[subsample, variable]) # cor.test(model.ABA.metrics.points$values$residual, plots[subsample, variable])
par(mfrow=c(1,2)) par(mfrow=c(1,2))
# plot predicted VS field values # plot predicted VS field values
lidaRtRee::ABAmodelPlot(model.ABA.point.metrics, main = variable, lidaRtRee::ABAmodelPlot(model.ABA.metrics.points, main = variable,
col = ifelse(plots$stratum == "public", "green", "blue")) col = ifelse(plots$stratum == "public", "green", "blue"))
legend("topleft", c("public", "private"), col = c("green", "blue"), pch = 1) legend("topleft", c("public", "private"), col = c("green", "blue"), pch = 1)
plot(plots[subsample, c("G.m2.ha")], plot(plots[subsample, c("G.m2.ha")],
model.ABA.point.metrics$values$residual, model.ABA.metrics.points$values$residual,
ylab = "Prediction errors", xlab = "Field values", ylab = "Prediction errors", xlab = "Field values",
col = ifelse(plots$stratum == "public", "green", "blue")) col = ifelse(plots$stratum == "public", "green", "blue"))
abline(h = 0, lty = 2) abline(h = 0, lty = 2)
...@@ -279,7 +279,7 @@ for (i in levels(plots[,strat])) ...@@ -279,7 +279,7 @@ for (i in levels(plots[,strat]))
subsample <- which(plots[,strat]==i) subsample <- which(plots[,strat]==i)
if (length(subsample)>0) if (length(subsample)>0)
{ {
model.ABA.stratified.none[[i]] <- lidaRtRee::ABAmodel(plots[subsample, variable], point.metrics[subsample,], transform="none", xy = plots[subsample,c("X", "Y")]) model.ABA.stratified.none[[i]] <- lidaRtRee::ABAmodel(plots[subsample, variable], metrics.points[subsample,], transform="none", xy = plots[subsample,c("X", "Y")])
} }
} }
# combine list of models into single object # combine list of models into single object
......
...@@ -68,25 +68,35 @@ sf::st_crs(public) <- sf::st_crs(forest) <- 2154 ...@@ -68,25 +68,35 @@ sf::st_crs(public) <- sf::st_crs(forest) <- 2154
## ALS metrics mapping ## ALS metrics mapping
First the ALS metrics used as independent variables in the ABA prediction models have to be mapped on the whole area. For this the SAME functions used for computing the metrics from the plot-extent cloud metrics should be used. They are used inside the function `lidR::grid_metrics` which performs the computation the area of interest divided in square pixels. Pixel surface should be of the same size as calibration plots. Using round values is advisable in case other remote sensing data are considered for use. First the ALS metrics used as independent variables in the ABA prediction models have to be mapped on the whole area. For this the SAME functions used for computing the metrics from the plot-extent cloud metrics should be used. Metrics are computed in each cell of the area of interest divided in square pixels. Pixel surface should be of the same size as calibration plots. Using round values is advisable in case other remote sensing data are considered for use.
```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6} ```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
# resolution for metrics map # resolution for metrics map
resolution <- 25 resolution <- 25
``` ```
### Example on one tile The following paragraph show how to compute metrics for a 100 x 100 m^2^ subsample of the first ALS tile.
The following code computes point metrics for a 100 x 100 m^2^ subsample of the first ALS tile. ```{r loadTile, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 6}
```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
# load ALS tile # load ALS tile
point.cloud.height <- lidR::clip_rectangle(cata.height, 900100, 6447100, 900400, 6447400) point.cloud.height <- lidR::clip_rectangle(cata.height, 899600, 6447100, 899900, 6447400)
point.cloud.ground <- lidR::filter_ground(lidR::clip_rectangle(cata.altitude, 900100, 6447100, 900400, 6447400)) point.cloud.ground <- lidR::filter_ground(lidR::clip_rectangle(cata.altitude, 899600, 6447100, 899900, 6447400))
```
### 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.
```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
# compute point metrics # compute point metrics
point.metrics <- lidR::grid_metrics(point.cloud.height, ~lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2), res=resolution) metrics.points <- lidR::grid_metrics(point.cloud.height, aba.pointMetricsFUN, res=resolution)
raster::plot(point.metrics[c("")]) par(mfrow=c(1,3))
for (i in c("zmean", "zsd", "imean"))
{
raster::plot(metrics.points[[i]], main = i)
}
``` ```
### Tree metrics
For tree metrics, three steps are required: For tree metrics, three steps are required:
...@@ -94,7 +104,9 @@ For tree metrics, three steps are required: ...@@ -94,7 +104,9 @@ For tree metrics, three steps are required:
* summary statistics of tree features by target pixels (values calculated based on trees which apices are inside the pixel); * summary statistics of tree features by target pixels (values calculated based on trees which apices are inside the pixel);
* additional statistics about the percentage of cover and the mean canopy height (values calculated based on segmented tree crowns which are inside the target pixel). * additional statistics about the percentage of cover and the mean canopy height (values calculated based on segmented tree crowns which are inside the target pixel).
```{r treeMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6} It is very important that tree segmentation parameters are the same as in the calibration steps, and that the function for aggregating tree attributes into raster metrics is the also the same as in the calibration step (except for the area which might differ slightly between the field plots and target cells).
```{r treeMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 6, fig.height = 6}
# compute chm # compute chm
chm <- lidaRtRee::points2DSM(point.cloud.height, res = aba.resCHM) chm <- lidaRtRee::points2DSM(point.cloud.height, res = aba.resCHM)
# replace NA, low and high values # replace NA, low and high values
...@@ -129,77 +141,239 @@ dummy <- lidaRtRee::rasterMetrics(r.treechm, ...@@ -129,77 +141,239 @@ dummy <- lidaRtRee::rasterMetrics(r.treechm,
names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot") names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot")
metrics.trees <- raster::addLayer(metrics.trees, dummy) metrics.trees <- raster::addLayer(metrics.trees, dummy)
# #
metrics.trees par(mfrow=c(1,1))
raster::plot(chm, main = "CHM and detected trees")
sp::plot(trees, cex = trees$h/60, add = TRUE)
#
```
```{r plottreeMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
par(mfrow=c(1,3))
for (i in c("Tree.meanH", "Tree.sdH", "Tree.density"))
{
raster::plot(metrics.trees[[i]], main = i)
}
``` ```
Terrain statistics can be computed from the altitude values of ground points. ### 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`.
```{r terrainMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6} ```{r terrainMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
# compute terrain metrics # compute terrain metrics
f <- function(x, y, z) f <- function(x, y, z)
{ {
as.list(lidaRtRee::points2terrainStats(data.frame(x, y, z))) as.list(lidaRtRee::points2terrainStats(data.frame(x, y, z)))
} }
terrain.metrics <- 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)
par(mfrow=c(1,3))
#for (i in c("altitude", "azimut.gr", "slope.gr"))
#{
# raster::plot(metrics.terrain[[i]], main = i)
#}
raster::plot(metrics.terrain[["altitude"]], main = "altitude")
raster::plot(metrics.terrain$azimut.gr, main = "azimut (gr)", breaks = seq(from = 0, to = 400, length.out = 40), col=rainbow(40))
raster::plot(metrics.terrain$slope.gr, main = "Slope (gr)")
``` ```
```{r old, include = FALSE} ### Batch processing of tiles
# workflow associated with package lidaRtRee
# Copyright Irstea For the batch processing of tiles, parallel processing with packages `foreach` and `doFuture` is used. The processing capabilities of `lidR` are not yet used.
# Author(s): Jean-Matthieu Monnet
# Licence: CeCill 2.1 ```{r setupCluster, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
#######################" library(foreach)
# workflow for computation of ALS metrics on LAZ dataset # create parallel frontend, specify to use two parallel sessions
# application of prediction model to obtain forest parameter map doFuture::registerDoFuture()
# inference function to compute estimates from whole map future::plan("multisession", workers = 2L)
# ```
# resolution for metrics map
resolution <- 25 A buffering procedure is designed to handle the border effects when detecting trees at tile edges. A 10 m buffer is enough for tree metrics. One can also specify points classes to use or apply a threshold to high points. Some parameters specified in the previous paragraphs are also required for the batch processing (output map resolution, chm resolution for tree segmentation, functions for metrics computation).
taille.dalle <- 1000
################ ```{r paramCluster, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
# metrics computation # processing parameters
# b.size <- 10
# compute metrics if not already done # height threshold (m) for high points removal (points above this threshold are considered
if (1)#!file.exists("./data/map.metrics.rda")) # as outliers)
{ h.points <- 60
# create catalog of normalized LAZ files : should be rectangular tiles, with edges aligned on multiple values of the resolution # points classes to retain for analysis (vegetation + ground)
cata <- lidR::catalog("/media/reseau/services/LESSEM/ProjetsCommuns/Lidar/data/73-Bauges/norm.laz.pnr") # ground should be 2
# cata <- lidR::catalog("/media/data/las/4montagnes/laznorm/") class.points <- c(0, 1, 2, 3, 4, 5)
library(foreach) ```
doFuture::registerDoFuture()
future::plan("multisession", workers = 2L) This first chunk of code computes tree and point metrics from the normalized ALS tiles.
#
map.metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar% ```{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%
{ {
centroid.x <- mean(as.numeric(cata@data[i,c("Min.X", "Max.X")])) # tile extent
centroid.y <- mean(as.numeric(cata@data[i,c("Min.Y", "Max.Y")])) b.box <- as.numeric(cata.height@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
origine.x <- floor(centroid.x/taille.dalle)*taille.dalle # load tile extent plus buffer
origine.y <- floor(centroid.y/taille.dalle)*taille.dalle a <- try(lidR::clip_rectangle(cata.height, b.box[1]-b.size, b.box[2]-b.size, b.box[3]+b.size,
a <- lidR::readLAS(cata$filename[i]) b.box[4]+b.size))
a <- lidR::lasfilter(a, X>=origine.x & X<origine.x+taille.dalle & Y>=origine.y & Y<origine.y+taille.dalle) #
# check if points are successfully loaded
if (class(a)=="try-error") return(NULL)
# add 'buffer' flag to points in buffer with TRUE value in this new attribute
a <- lidR::add_attribute(a,
a$X < b.box[1] | a$Y < b.box[2] | a$X >= b.box[3] | a$Y >= b.box[4],
"buffer")
# remove unwanted point classes, and points higher than height threshold
a <- lidR::filter_poi(a, is.element(Classification, class.points) & Z<=h.points)
# check number of remaining points
if (length(a)==0) return(NULL)
# set negative heights to 0 # set negative heights to 0
a@data$Z[a@data$Z<0] <- 0 a@data$Z[a@data$Z<0] <- 0
dummy <- lidR::grid_metrics(a, lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2), res=resolution) #
# compute chm
return(dummy) chm <- lidaRtRee::points2DSM(a, res=aba.resCHM)
# replace NA, low and high values by 0
chm[is.na(chm) | chm<0 | chm>h.points] <- 0
#
# compute tree metrics
# tree top detection (default parameters)
segms <- lidaRtRee::treeSegmentation(chm)
# 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")
# compute canopy cover in trees and canopy mean height in trees
# in region of interest, because it is not in previous step.
r.treechm <- segms$filled.dem
# set chm to NA in non segment area
r.treechm[segms$segments.id==0] <- NA
# compute raster metrics
metrics.trees2 <- lidaRtRee::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/aba.resCHM)^2,
mean(x$filled.dem,na.rm=T))
}
, output="raster")
names(metrics.trees2) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot")
#
# compute 1D height metrics
# remove buffer points
a <- lidR::filter_poi(a, buffer==0)
#
if (length(a)==0) return(NULL)
# all points metrics
metrics.points <- lidR::grid_metrics(a, aba.pointMetricsFUN, res = resolution)
#
# extend / crop to match metrics.points
metrics.trees <- raster::extend(metrics.trees, metrics.points, values = NA)
metrics.trees2 <- raster::extend(metrics.trees2, metrics.points, values = NA)
metrics.trees <- raster::crop(metrics.trees, metrics.points)
metrics.trees2 <- raster::crop(metrics.trees2, metrics.points)
# merge rasterstacks
metrics <- metrics.points
metrics <- raster::addLayer(metrics, metrics.trees)
metrics <- raster::addLayer(metrics, metrics.trees2)
return(metrics)
} }
# rbind data.frame with option "fill" to handle data.frames with 0 rows and less columns # mosaic rasters
names.metrics <- names(map.metrics[[1]]) names.metrics <- names(metrics[[1]])
test <- Reduce(raster::merge, map.metrics) metrics.map <- do.call(raster::merge, metrics)
map.metrics <- do.call(raster::merge, map.metrics) names(metrics.map) <- names.metrics
names(map.metrics) <- names.metrics ```
} else {load("./data/map.metrics.rda")}
The second chunk of code computes terrain metrics from the ALS tiles with altitude values.
```{r batchProcessingTerrain, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
# function to compute terrain metrics to input to grid_metrics
f <- function(x, y, z)
{
as.list(lidaRtRee::points2terrainStats(data.frame(x, y, z)))
}
# #
# ADD TERRAIN METRICS <- REQUIRES NON NORMALIZED POINT CLOUD metrics.terrain <- foreach::foreach (i=1:nrow(cata.altitude), .errorhandling="remove") %dopar%
# ADD TREE METRICS <- REQUIRES ADDITIONNAL CODE TO COMPUTE THEM ON GRID {
################# # tile extent
# rasterize strata b.box <- as.numeric(cata.altitude@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
public <- rgdal::readOGR("./data/", "Public4Montagnes") # load tile extent plus buffer
map.metrics@crs <- public@proj4string a <- try(lidR::clip_rectangle(cata.altitude, b.box[1], b.box[2], b.box[3],
map.metrics$strata <- raster::rasterize(public, map.metrics) b.box[4]))
map.metrics$strata[is.na(map.metrics$strata)] <- 0 # check if points are successfully loaded
levels(map.metrics$strata) <- data.frame(ID=c(0,1), propriete=c("private", "public")) if (class(a)=="try-error") return(NULL)
# keep only ground points
a <- lidR::filter_ground(a)
lidR::grid_metrics(a, ~f(X, Y, Z), res=resolution)
}
# mosaic rasters
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, fig.width = 12, fig.height = 4}
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)
```
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, fig.width = 12, fig.height = 4}
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`.
```{r mapForest, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
prediction.map <- lidaRtRee::ABApredict(model.ABA, metrics.map)
```
### Stratified models
In case strata-specific models have been calibrated, it is necessary to add a layer containing the strata information in the metrics raster. For this, the existing vector layer of public forests has to be rasterized and added to the metrics raster. The levels in the layer metadata have to be filled, both in the ID (code, numeric) and stratum (label, character) fields.
```{r rasterizeStrat, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
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
# label levels in the raster metadata
levels(metrics.map$stratum) <- data.frame(ID=c(0,1), stratum=c("private", "public"))
# plot raster and vector
raster::plot(metrics.map$stratum, legend=FALSE, col=c("green", "blue"),xaxt='n', yaxt='n')
legend("bottom", legend=c("private","public"), fill=c("green","blue"))
# crop public vector
public.cropped <- sf::st_crop(public, raster::extent(metrics.map))
plot(public.cropped, col = NA, add = TRUE)
```
```{r rasterizeStrat, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4}
prediction.map.mixed <- lidaRtRee::ABApredict(model.ABA.stratified.mixed, metrics.map, stratum="stratum")
# verification par comparaison avec model prive / public
model.private <- list(model = model.ABA.stratified.mixed$model$private, stats = model.ABA.stratified.mixed$stats["private",])
prediction.map.private <- lidaRtRee::ABApredict(model.private, metrics.map)
model.public <- list(model = model.ABA.stratified.mixed$model$public, stats = model.ABA.stratified.mixed$stats["public",])
prediction.map.public <- lidaRtRee::ABApredict(model.public, metrics.map)
```
```{r old, include = FALSE}
} else {load("./data/map.metrics.rda")}
# #
################ ################
# requires previous calibration of models (script workflow.prediction.model.calibration.R) # requires previous calibration of models (script workflow.prediction.model.calibration.R)
......
...@@ -467,7 +467,7 @@ metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar% ...@@ -467,7 +467,7 @@ metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar%
a@data$Z[a@data$Z<0] <- 0 a@data$Z[a@data$Z<0] <- 0
# #
# compute chm # compute chm
chm <- lidaRtRee::points2DSM(a,res=resCHM) chm <- lidaRtRee::points2DSM(a, res=resCHM)
# replace NA, low and high values # replace NA, low and high values
chm[is.na(chm) | chm<0 | chm>h.points] <- 0 chm[is.na(chm) | chm<0 | chm>h.points] <- 0
# ########################## # ##########################
......
File deleted
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment