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

Added forest structure metrics

parent ddb5729d
...@@ -70,7 +70,7 @@ ISSN={1999-4907}, ...@@ -70,7 +70,7 @@ ISSN={1999-4907},
number = {1}, number = {1},
pages = {5-19}, pages = {5-19},
keywords = {Habitat suitability models, LiDAR, object-oriented metrics, point-cloud area-based metrics}, keywords = {Habitat suitability models, LiDAR, object-oriented metrics, point-cloud area-based metrics},
doi = {https://doi.org/10.1002/rse2.117}, doi = {10.1002/rse2.117},
url = {https://zslpublications.onlinelibrary.wiley.com/doi/abs/10.1002/rse2.117}, url = {https://zslpublications.onlinelibrary.wiley.com/doi/abs/10.1002/rse2.117},
eprint = {https://zslpublications.onlinelibrary.wiley.com/doi/pdf/10.1002/rse2.117}, eprint = {https://zslpublications.onlinelibrary.wiley.com/doi/pdf/10.1002/rse2.117},
abstract = {Abstract Habitat suitability models (HSMs) are widely used to plan actions for species of conservation interest. Models that will be turned into conservation actions need predictors that are both ecologically pertinent and fit managers’ conceptual view of ecosystems. Remote sensing technologies such as light detection and ranging (LiDAR) can describe landscapes at high resolution over large spatial areas and have already given promising results for modeling forest species distributions. The point-cloud (PC) area-based LiDAR variables are often used as environmental variables in HSMs and have more recently been complemented by object-oriented (OO) metrics. However, the efficiency of each type of variable to capture structural information on forest bird habitat has not yet been compared. We tested two hypotheses: (1) the use of OO variables in HSMs will give similar performance as PC area-based models; and (2) OO variables will improve model robustness to LiDAR datasets acquired at different times for the same area. Using the case of a locally endangered forest bird, the capercaillie (Tetrao urogallus), model performance and predictions were compared between the two variable types. Models using OO variables showed slightly lower discriminatory performance than PC area-based models (average ΔAUC = −0.032 and −0.01 for females and males, respectively). OO-based models were as robust (absolute difference in Spearman rank correlation of predictions ≤ 0.21) or more robust than PC area-based models. In sum, LiDAR-derived PC area-based metrics and OO metrics showed similar performance for modeling the distribution of the capercaillie. We encourage the further exploration of OO metrics for creating reliable HSMs, and in particular testing whether they might help improve the scientist–stakeholder interface through better interpretability.}, abstract = {Abstract Habitat suitability models (HSMs) are widely used to plan actions for species of conservation interest. Models that will be turned into conservation actions need predictors that are both ecologically pertinent and fit managers’ conceptual view of ecosystems. Remote sensing technologies such as light detection and ranging (LiDAR) can describe landscapes at high resolution over large spatial areas and have already given promising results for modeling forest species distributions. The point-cloud (PC) area-based LiDAR variables are often used as environmental variables in HSMs and have more recently been complemented by object-oriented (OO) metrics. However, the efficiency of each type of variable to capture structural information on forest bird habitat has not yet been compared. We tested two hypotheses: (1) the use of OO variables in HSMs will give similar performance as PC area-based models; and (2) OO variables will improve model robustness to LiDAR datasets acquired at different times for the same area. Using the case of a locally endangered forest bird, the capercaillie (Tetrao urogallus), model performance and predictions were compared between the two variable types. Models using OO variables showed slightly lower discriminatory performance than PC area-based models (average ΔAUC = −0.032 and −0.01 for females and males, respectively). OO-based models were as robust (absolute difference in Spearman rank correlation of predictions ≤ 0.21) or more robust than PC area-based models. In sum, LiDAR-derived PC area-based metrics and OO metrics showed similar performance for modeling the distribution of the capercaillie. We encourage the further exploration of OO metrics for creating reliable HSMs, and in particular testing whether they might help improve the scientist–stakeholder interface through better interpretability.},
......
This diff is collapsed.
...@@ -351,7 +351,7 @@ raster::plot(edgeOutside, main = "Edges (detection by dilation)", lege ...@@ -351,7 +351,7 @@ raster::plot(edgeOutside, main = "Edges (detection by dilation)", lege
<h2 class="unnumbered">References</h2> <h2 class="unnumbered">References</h2>
<div id="refs" class="references csl-bib-body hanging-indent"> <div id="refs" class="references csl-bib-body hanging-indent">
<div id="ref-Glad20" class="csl-entry"> <div id="ref-Glad20" class="csl-entry">
Glad, Anouk, Björn Reineking, Marc Montadert, Alexandra Depraz, and Jean-Matthieu Monnet. 2020. <span>“Assessing the Performance of Object-Oriented LiDAR Predictors for Forest Bird Habitat Suitability Modeling.”</span> <em>Remote Sensing in Ecology and Conservation</em> 6 (1): 5–19. https://doi.org/<a href="https://doi.org/10.1002/rse2.117">https://doi.org/10.1002/rse2.117</a>. Glad, Anouk, Björn Reineking, Marc Montadert, Alexandra Depraz, and Jean-Matthieu Monnet. 2020. <span>“Assessing the Performance of Object-Oriented LiDAR Predictors for Forest Bird Habitat Suitability Modeling.”</span> <em>Remote Sensing in Ecology and Conservation</em> 6 (1): 5–19. <a href="https://doi.org/10.1002/rse2.117">https://doi.org/10.1002/rse2.117</a>.
</div> </div>
</div> </div>
</div> </div>
...@@ -17,15 +17,15 @@ knitr::opts_chunk$set(fig.align = "center") ...@@ -17,15 +17,15 @@ knitr::opts_chunk$set(fig.align = "center")
``` ```
--- ---
Licence: CC-BY / Source page Licence: CC-BY / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/forest.structure.metrics.Rmd)
The R code below presents a forest structure metrics computation workflow from Airborne Laser Scanning (ALS) data. Workflow is based on functions from packages `lidaRtRee` and `lidR`, packages `vegan` and `foreach` are also required. Metrics are computed for each cell of a grid defined by a resolution. Those metrics are designed to describe the 3D structure of forest. The R code below presents a forest structure metrics computation workflow from Airborne Laser Scanning (ALS) data. Workflow is based on functions from packages `lidaRtRee` and `lidR`, packages `vegan` and `foreach` are also required. Metrics are computed for each cell of a grid defined by a resolution. Those metrics are designed to describe the 3D structure of forest.
Different types of metrics are computed: Different types of metrics are computed:
* 1D height metrics * 1D height metrics
+ 2D metrics of the canopy height model (CHM) + 2D metrics of the canopy height model (CHM)
+ tree segmentation metrics + tree segmentation metrics (see also [tree segmentation tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/tree.detection.Rmd))
+ forest gaps and edges metrics + forest gaps and edges metrics (see also [gaps and edges detection tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/gaps.edges.detection.Rmd))
The forest structure metrics derived from airborne laser scanning can be used for habitat suitability modelling and mapping. This workflow has been applied to compute the metrics used in the modeling and mapping of the habitat of Capercaillie (*Tetrao urogallus*)(@Glad20). For more information about tree segmentation and gaps detection, please refer to the corresponding tutorials. The forest structure metrics derived from airborne laser scanning can be used for habitat suitability modelling and mapping. This workflow has been applied to compute the metrics used in the modeling and mapping of the habitat of Capercaillie (*Tetrao urogallus*)(@Glad20). For more information about tree segmentation and gaps detection, please refer to the corresponding tutorials.
...@@ -126,7 +126,7 @@ The next step is to compute the canopy height model (CHM). It will be used to de ...@@ -126,7 +126,7 @@ The next step is to compute the canopy height model (CHM). It will be used to de
* 2D canopy height metrics related to multi-scale vertical heterogeneity (mean and standard deviation of CHM, smoothed at different scales) * 2D canopy height metrics related to multi-scale vertical heterogeneity (mean and standard deviation of CHM, smoothed at different scales)
+ tree metrics from tree top segmentation + tree metrics from tree top segmentation
+ gaps and edges metrics from gap extraction + gaps and edges metrics from gap segmentation
The CHM is computed and NA values are replaced by 0. A check is performed to make sure low or high points are not present. The CHM is computed and NA values are replaced by 0. A check is performed to make sure low or high points are not present.
...@@ -351,8 +351,10 @@ Some outputs are displayed hereafter. ...@@ -351,8 +351,10 @@ Some outputs are displayed hereafter.
par(mfrow=c(1,3)) 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) 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,
raster::plot(metrics.trees$Tree.meanCrownVolume, main="Mean crown volume of detected trees (m3)") main="Density of (detected) trees > 20m (/ha)")
raster::plot(metrics.trees$Tree.meanCrownVolume,
main="Mean crown volume of detected trees (m3)")
``` ```
### Point cloud (1D) metrics ### Point cloud (1D) metrics
...@@ -394,8 +396,6 @@ Based on those metrics, additional metrics are computed (Simpson index, relative ...@@ -394,8 +396,6 @@ Based on those metrics, additional metrics are computed (Simpson index, relative
metrics.1d$H.simpson <- raster::stackApply(metrics.1d[[n.breaksH[c(-1,-length(n.breaksH))]]], metrics.1d$H.simpson <- raster::stackApply(metrics.1d[[n.breaksH[c(-1,-length(n.breaksH))]]],
1, 1,
function(x, ...){vegan::diversity(x, index="simpson")}) function(x, ...){vegan::diversity(x, index="simpson")})
#
# convert to points HERE (easier-> avoid for loops)
# Relative density # Relative density
for (i in n.breaksH[c(-1,-length(n.breaksH))]) for (i in n.breaksH[c(-1,-length(n.breaksH))])
{ {
...@@ -404,23 +404,22 @@ for (i in n.breaksH[c(-1,-length(n.breaksH))]) ...@@ -404,23 +404,22 @@ for (i in n.breaksH[c(-1,-length(n.breaksH))])
} }
# Penetration ratio # Penetration ratio
# cumulative sum # cumulative sum
# cumS <- t(apply(metrics.1d[,c("nbSol",n.breaksH[-length(n.breaksH)])],1,cumsum)) cumu.sum <- metrics.1d[["nbSol"]]
cum.sum <- metrics.1d[["nbSol"]]
for (i in n.breaksH) for (i in n.breaksH)
{ {
cum.sum <- raster::addLayer(cum.sum, cum.sum[[dim(cum.sum)[3]]] + metrics.1d[[i]]) cumu.sum <- raster::addLayer(cumu.sum, cumu.sum[[dim(cumu.sum)[3]]] + metrics.1d[[i]])
} }
names(cum.sum) <- c("nbSol", n.breaksH) names(cumu.sum) <- c("nbSol", n.breaksH)
# penetration ratio WHY USE 1-ratio -> should be interception ratio ? # compute interception ratio of each layer
pen.ratio <- cum.sum[[-1]] intercep.ratio <- cumu.sum[[-1]]
for (i in 2:dim(cum.sum)[3]) for (i in 2:dim(cumu.sum)[3])
{ {
pen.ratio[[i-1]] <- 1 - cum.sum[[i-1]]/cum.sum[[i]] intercep.ratio[[i-1]] <- 1 - cumu.sum[[i-1]]/cumu.sum[[i]]
} }
names(pen.ratio) <- paste0(names(cum.sum)[-1], "ratio") names(intercep.ratio) <- paste0(names(cumu.sum)[-1], "ratio")
# merge rasterstacks # merge rasterstacks
metrics.1d <- raster::addLayer(metrics.1d, pen.ratio) metrics.1d <- raster::addLayer(metrics.1d, intercep.ratio)
rm(cum.sum, pen.ratio) # rm(cumu.sum, intercep.ratio)
# #
metrics.1d metrics.1d
``` ```
...@@ -432,12 +431,12 @@ Some outputs are displayed hereafter. ...@@ -432,12 +431,12 @@ Some outputs are displayed hereafter.
par(mfrow=c(1,3)) par(mfrow=c(1,3))
raster::plot(metrics.1d$I.mean.1stpulse, main="Mean intensity (1st return)") 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.p50, main="Percentile 50 of heights")
raster::plot(metrics.1d$H.nb2_5ratio, main="Penetration ratio 2-5 m") raster::plot(metrics.1d$H.nb2_5ratio, main="Interception ratio 2-5 m")
``` ```
## Batch processing ## Batch processing
The following code uses parallel processing to handle multiple files and arranges all metrics in a single dataframe. Code in the "parameters" section has to be run beforehand. The following code uses parallel processing to handle multiple files and arranges all metrics in a raster stack. Code in the "parameters" section has to be run beforehand.
```{r batchProcessing, include = TRUE, eval=TRUE, warning=FALSE, message=FALSE, fig.width = 9, fig.height = 2.6} ```{r batchProcessing, include = TRUE, eval=TRUE, warning=FALSE, message=FALSE, fig.width = 9, fig.height = 2.6}
# processing by tile # processing by tile
...@@ -487,9 +486,29 @@ metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar% ...@@ -487,9 +486,29 @@ metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar%
# multiplying factor to compute proportion # multiplying factor to compute proportion
mf <- 100/(resolution/resCHM)^2 mf <- 100/(resolution/resCHM)^2
# compute raster metrics # 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) metrics.2dchm <-
}, lidaRtRee::rasterMetrics(
output="raster") 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 # compute gap metrics
...@@ -623,7 +642,8 @@ metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar% ...@@ -623,7 +642,8 @@ metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar%
metrics.edges <- raster::crop(metrics.edges, metrics.2dchm) metrics.edges <- raster::crop(metrics.edges, metrics.2dchm)
metrics.trees <- raster::crop(metrics.trees, metrics.2dchm) metrics.trees <- raster::crop(metrics.trees, metrics.2dchm)
# #
temp <- raster::addLayer(metrics.1d, metrics.2dchm, metrics.gaps, metrics.edges, metrics.trees) temp <- raster::addLayer(metrics.1d, metrics.2dchm, metrics.gaps,
metrics.edges, metrics.trees)
temp[is.na(temp)] <- 0 temp[is.na(temp)] <- 0
} # end of vegetation-only points presence check } # end of vegetation-only points presence check
} # end of buffer-less points presence check } # end of buffer-less points presence check
...@@ -655,22 +675,22 @@ for (i in n.breaksH) ...@@ -655,22 +675,22 @@ for (i in n.breaksH)
cum.sum <- raster::addLayer(cum.sum, cum.sum[[dim(cum.sum)[3]]] + metrics[[i]]) cum.sum <- raster::addLayer(cum.sum, cum.sum[[dim(cum.sum)[3]]] + metrics[[i]])
} }
names(cum.sum) <- c("nbSol", n.breaksH) names(cum.sum) <- c("nbSol", n.breaksH)
# absorption ratio # interception ratio
pen.ratio <- cum.sum[[-1]] intercep.ratio <- cum.sum[[-1]]
for (i in 2:dim(cum.sum)[3]) for (i in 2:dim(cum.sum)[3])
{ {
pen.ratio[[i-1]] <- 1 - cum.sum[[i-1]]/cum.sum[[i]] intercep.ratio[[i-1]] <- 1 - cum.sum[[i-1]]/cum.sum[[i]]
} }
names(pen.ratio) <- paste0(names(cum.sum)[-1], "ratio") names(intercep.ratio) <- paste0(names(cum.sum)[-1], "ratio")
# merge rasterstacks # merge rasterstacks
metrics <- raster::addLayer(metrics, pen.ratio) metrics <- raster::addLayer(metrics, intercep.ratio)
# #
# ########################### # ###########################
# export as raster files # export as raster files
# for (i in 1:names(metrics)) # for (i in 1:names(metrics))
# { # {
# print(i) # print(i)
# writeRaster(metrics[[i]],file=paste("output/raster_",i,"_",resolution,"m.tif",sep=""),overwrite=T) # writeRaster(metrics[[i]],file=paste("output/raster_",i,"_",resolution,"m.tif",sep="")
# } # }
# ########################### # ###########################
# display # display
...@@ -678,3 +698,4 @@ raster::plot(metrics[[c("I.mean.1stpulse","H.simpson")]], ...@@ -678,3 +698,4 @@ raster::plot(metrics[[c("I.mean.1stpulse","H.simpson")]],
main=c("Mean intensity of 1st pulse","Simpson indice of point heights")) main=c("Mean intensity of 1st pulse","Simpson indice of point heights"))
``` ```
## References
Markdown is supported
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