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

Minor modifications in plots of tree detection

Updated coregistration
gaps replaced by gaps.edges
parent b75f8427
......@@ -61,3 +61,18 @@ ISSN={1999-4907},
month={May},
pages={1721–1747}
}
@article{Glad20,
author = {Glad, Anouk and Reineking, Björn and Montadert, Marc and Depraz, Alexandra and Monnet, Jean-Matthieu},
title = {Assessing the performance of object-oriented LiDAR predictors for forest bird habitat suitability modeling},
journal = {Remote Sensing in Ecology and Conservation},
volume = {6},
number = {1},
pages = {5-19},
keywords = {Habitat suitability models, LiDAR, object-oriented metrics, point-cloud area-based metrics},
doi = {https://doi.org/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},
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.},
year = {2020}
}
......@@ -5,6 +5,7 @@ date: "`r Sys.Date()`"
output:
pdf_document: default
html_document: default
papersize: a4
bibliography: "./bib/bibliography.bib"
---
......@@ -154,16 +155,17 @@ r.mask[r.mask==0] <- NA
# specify projection
raster::crs(r.mask) <- 2154
# display plot mask
raster::plot(r.mask, main="Plot mask and tree positions")
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
points(centre, pch=3)
legend("topleft", c("Trees", "Plot center"), pch=c(1, 3))
```
### Compute correlation and identify optimal translation
First the function computes the correlation for different possible translations of the plot center inside the buffer specified by the user, and outputs an image of correlation values between the ALS CHM and the pseudo CHM. Second this image is analyzed to identify which translation yields the highest correlation. The function outputs a list which first element is the correlation image, and the second one the corresponding statistics.
First the function computes the correlation for different possible translations of the plot center inside the buffer specified by the user, and outputs an image of correlation values between the ALS CHM and the pseudo CHM. The pseudo CHM is a height model where at the location of inventoried trees pixels are attributed the value corresponding to tree size (e.g. height or diameter). Second the correlation image is analyzed to identify which translation yields the highest correlation. The function outputs a list which first element is the correlation image, and the second one the corresponding statistics.
```{r correlation, include = TRUE}
# compute correlation image
......@@ -172,20 +174,22 @@ coreg <- lidaRtRee::coregistration(chmfilt, trees[,c("x", "y", "dia")], mask=r.m
# correlation image statistics
round(coreg$local.max, 2)
```
The maximum of in 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))
raster::plot(chmfilt, main="Initial tree positions")
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, main="Correlation image")
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=3, col="blue")
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")
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,
cex=trees$dia/40, col="red")
......
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
---
title: "Workflow for forest gaps and edges detection from ALS data"
title: "R workflow for forest gaps and edges detection from ALS data"
author: "Jean-Matthieu Monnet"
date: "Jan 28, 2021"
date: "`r Sys.Date()`"
output:
pdf_document: default
html_document: default
bibliography: workflow.treedetection.bib
papersize: a4
bibliography: "./bib/bibliography.bib"
---
```{r setup, include=FALSE}
......@@ -16,10 +17,9 @@ knitr::opts_chunk$set(fig.align = "center")
```
---
The code below presents a forest gaps and edges detection workfow from Airborne Laser Scanning (ALS) data. Workflow is based on functions from packages `lidaRtRee` and `lidR`.
This tutorial presents R code for forest gaps and edges detection from Airborne Laser Scanning (ALS) data. The workflow is based on functions from packages `lidaRtRee` and `lidR`.
Licence: CC BY
Source page: https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee/-/wikis/Gaps-and-edges-detection-workflow
Licence: CC BY / Source page
Changelog
......
......@@ -3,8 +3,9 @@ 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"
---
......@@ -236,14 +237,13 @@ head(trees, n=3L)
v.segments <- raster::rasterToPolygons(segms[[2]], dissolve=T)
#
# display initial image
raster::plot(chm, col=gray(seq(0,1,1/255)),
main ="CHM and detected positions")
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)
# display plot mask
sp::plot(v.plotMask,border="red",add=T)
# display inventoried trees
graphics::points(trees$x, trees$y, col="blue", cex=trees$h/20)
graphics::points(trees$x, trees$y, col="blue", cex=trees$h/20, pch = 2)
```
## Detection evaluation
......@@ -411,8 +411,10 @@ rgl::plot3d(points.seg[,c("X", "Y", "Z")], col=points.seg$seg.id%%10 +1,aspect=F
treeinventorychablais3$z <- 0
for (i in 1:nrow(treeinventorychablais3))
{
rgl::lines3d(rbind(treeinventorychablais3$x[i]-974300, treeinventorychablais3$x[i]-974300),
rbind(treeinventorychablais3$y[i]-6581600,treeinventorychablais3$y[i]-6581600),
rgl::lines3d(rbind(treeinventorychablais3$x[i]-974300,
treeinventorychablais3$x[i]-974300),
rbind(treeinventorychablais3$y[i]-6581600,
treeinventorychablais3$y[i]-6581600),
rbind(treeinventorychablais3$z[i],
treeinventorychablais3$z[i]+treeinventorychablais3$h[i]))
}
......@@ -546,7 +548,8 @@ lidR::sensor(cata) <- "ALS"
# trade-off between RAM capacity VS total number of tiles to process
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
# 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
#
# TREE SEGMENTATION PARAMETERS
......@@ -565,7 +568,8 @@ save.chm <- FALSE
# working directory
wd <- getwd()
# run with two cores
clust <- parallel::makeCluster(getOption("cl.cores", 2))#, outfile = "/home/jean-matthieu/Bureau/output.txt")
clust <- parallel::makeCluster(getOption("cl.cores", 2))
#, outfile = "/home/jean-matthieu/Bureau/output.txt")
# export global variables to cluster because they will be called inside the function
parallel::clusterExport(cl = clust, ls(), envir = .GlobalEnv)
#
......@@ -618,9 +622,14 @@ resultats <- parallel::mcmapply(
if (!all(is.na(raster::values(chm))))
{
# output chm if asked for
if (output.chm) output$chm <- raster::crop(chm, raster::extent(i, i+tile.size, j, j+tile.size))
if (output.chm) output$chm <- raster::crop(chm, raster::extent(i, i+tile.size,
j, j+tile.size))
# save on disk
if (save.chm) raster::writeRaster(raster::crop(chm, raster::extent(i, i+tile.size, j, j+tile.size)), file= paste0("chm_", i, "_", j, ".tif"), overwrite = TRUE)
if (save.chm) raster::writeRaster(raster::crop(chm,
raster::extent(i, i+tile.size,
j, j+tile.size)),
file= paste0("chm_", i, "_", j, ".tif"),
overwrite = TRUE)
#
# tree detection
segms <- lidaRtRee::treeSegmentation(chm)
......@@ -629,7 +638,8 @@ resultats <- parallel::mcmapply(
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
......@@ -642,7 +652,8 @@ resultats <- parallel::mcmapply(
# names(v.trees) <- "id"
# add attributes
# errors when using sp::merge so using sp::over even if it is probably slower
v.trees@data <- cbind(v.trees@data, sp::over(v.trees, trees))#merge(v.trees@data, trees, all.x = TRUE)
# merge(v.trees@data, trees, all.x = TRUE)
v.trees@data <- cbind(v.trees@data, sp::over(v.trees, trees))
# save in list
output$v.trees <- v.trees
}
......@@ -681,7 +692,8 @@ if (vectorize.trees)
# 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 is not fully identical in overlap areas of adjacent tiles.
# 1-pixel overlapping in v.trees might be present because image segmentation
# is not fully identical in overlap areas of adjacent tiles.
}
```
......@@ -697,7 +709,7 @@ raster::plot(chm.all,
# display segments border
sp::plot(v.trees, border = "white", add = T)
# add trees
sp::plot(trees, cex = trees$h/40, add = TRUE)
sp::plot(trees, cex = trees$h/40, add = TRUE, pch = 2)
```
The following lines save outputs to files.
......
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