Skip to content
GitLab

tree detection workflow · Changes

Page history
Update tree detection workflow authored by Monnet Jean-Matthieu's avatar Monnet Jean-Matthieu
Hide whitespace changes
Inline Side-by-side
tree-detection-workflow.md
View page @ 0b4027f5
[workflow.treedetection.Rmd](/uploads/223064d4cfa07260a7e78e206d66b0d9/workflow.treedetection.Rmd) ---
\ No newline at end of file title: "Workflow for tree detection"
author: "Jean-Matthieu Monnet"
date: "April 27, 2018"
output: html_document
bibliography: workflow.treedetection.bib
---
```{r setup, include=FALSE}
library(knitr)
knitr::opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(webgl = hook_webgl)
```
---
The code below presents a tree detection workfow from Airborne Laser Scanning (ALS) data. Workflow is based on functions from packages lidaRtRee and lidR.
```{r include = FALSE}
library(lidR)
```
* treetop detection and crown segmentation
+ accuracy assessment with field inventory
+ species classification
## Material
### Field inventory
The field inventory corresponds to a 50m x 50m plot located in the Chablais mountain (France) [@Monnet11c, pp. 34]. All trees with a diameter at breast height above 7.5 cm are inventoried. The data is available in package lidaRtRee.
```{r include = TRUE}
# load dataset from package (default)
data(treeinventorychablais3, package="lidaRtRee")
```
Otherwise you can load your own data provided positions and heights are measured.
```{r eval=FALSE}
# import field inventory
fichier <- "lidaRtRee/inst/extdata/chablais3_listeR.csv"
tree.inventory <- read.csv(file=fichier,sep=";",header=F)
names(tree.inventory) <- c("x","y","d","h","n","s","e","t")
# save as rda for later access
save(tree.inventory,file="tree.inventory.rda")
```
Attributes are:
* x easting coordinate in Lambert 93
+ y northing coordinate in Lambert 93
+ d dbh (cm)
+ h tree height (m)
+ n tree number
+ s species abreviated as GESP (GEnus SPecies)
+ e appearance (0: missing or lying, 1: normal, 2: broken treetop, 3: dead with branches, 4: snag)
+ t tilted (0: no, 1: yes)
```{r echo=FALSE}
head(treeinventorychablais3)
```
Function plotTreeInventory is designed to plot forest inventory data.
```{r include = TRUE}
# display inventoried trees
lidaRtRee::plotTreeInventory(treeinventorychablais3[,c("x","y")], treeinventorychablais3$h, species=as.character(treeinventorychablais3$s), add=FALSE)
```
We define the region of interest (ROI) to crop ALS data to corresponding extent before further processing. ROI is set on the extent of tree inventory, plus a 10 meter buffer.
```{r include = TRUE}
# buffer to apply around ROI (meters)
ROI.buffer <- 10
# ROI limits
ROI.range <- apply(treeinventorychablais3[,c("x","y")],2,range)
# ROI center
ROI.center <- round((ROI.range[2,]+ROI.range[1,])/2)
# ROI dimensions
ROI.size <- round((ROI.range[2,]-ROI.range[1,])/2)
```
### Airborne Laser Scanning data
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 include = FALSE}
# load data in package lidaRtRee (default)
data(laschablais3, package="lidaRtRee")
```
```{r eval=FALSE}
# directory for laz files
lazdir <- "/media/data/las/chab/laz/"
# build catalog
cata <- lidR::catalog(lazdir)
# extract points in ROI
las <- lidR::catalog_queries(cata,ROI.center[1],ROI.center[2],ROI.size[1]+ROI.buffer+5,ROI.size[2]+ROI.buffer+5,ScanAngle=F, RGB=F, pulseID=F)
laschablais3 <- las[[1]]
# save as rda for easier access
save(las, file="las.rda")
```
```{r echo=FALSE}
laschablais3
```
## Data preparation
### Digital Elevation Models
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 substracting the DTM to the DSM.
```{r include=TRUE}
# terrain model
dtm <- lidaRtRee::points2DTM(laschablais3 %>% lasfilter(Classification==2),res=0.5, ROI.center[1]-ROI.size[1]-ROI.buffer, ROI.center[1]+ROI.size[1]+ROI.buffer, ROI.center[2]-ROI.size[2]-ROI.buffer, ROI.center[2]+ROI.size[2]+ROI.buffer)
# surface model
dsm <- lidaRtRee::points2DSM(laschablais3,res=0.5, ROI.center[1]-ROI.size[1]-ROI.buffer, ROI.center[1]+ROI.size[1]+ROI.buffer, ROI.center[2]-ROI.size[2]-ROI.buffer, ROI.center[2]+ROI.size[2]+ROI.buffer)
# canopy height model
chm <- dsm - dtm
```
```{r echo=FALSE, fig.width = 12, fig.height = 4}
# set coordinate reference
chm@crs <- dtm@crs <- dsm@crs <- sp::CRS("+init=epsg:2154")
#
par(mfrow=c(1,3))
# display DTM
raster::plot(dtm,asp=1, main="DTM")
# display DSM
raster::plot(dsm,asp=1, main="DSM")
# display CHM
raster::plot(chm,asp=1, main="CHM")
```
### Visual comparison of field inventory and ALS data
A plot mask is computed from the inventoried positions, using a height-dependent buffer. Indeed, tree tops are not necessairily located verticaly above the trunk positions. In order to compare detected tree tops with inventoried trunks, buffers have to be applied to trunk positions to account for the non-verticality of trees.
```{r include=TRUE, message=FALSE}
# select trees with height measures
selec <- which(!is.na(treeinventorychablais3$h))
# plot mask computation based on inventoried positions
# convex hull of plot
ChullMask <- lidaRtRee::rasterChullMask(treeinventorychablais3[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)
# union of convexHull and tree buffers
plotMask <- max(ChullMask , TreeMask)
# vectorize plot mask
v.plotMask <- raster::rasterToPolygons(plotMask,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.
```{r include = TRUE, fig.width = 6.5, fig.height = 6, warnings=FALSE}
# display CHM
raster::plot(chm,asp=1, col=gray(seq(0,1,1/255)), 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)
# display plot mask
raster::plot(v.plotMask,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:
* image pre-processing (non-linear filtering and smoothing for noise removal)
+ local maxima filtering and selection for tree top detection
+ image segmentation with a watershed algorithm for tree delineation.
The first two steps are documented in [@Monnet11c, pp. 47-52]
```{r include=TRUE, fig.width = 12, fig.height = 4}
# tree detection (default settings), applied on canopy height model
segms <- lidaRtRee::treeSegmentation(chm)
#
par(mfrow=c(1,3))
# display pre-processed chm
raster::plot(segms$smoothed.dem,asp=1, main="Pre-processed CHM")
# display selected local maxima
raster::plot(segms$local.maxima,asp=1, main="Selected local maxima")
# display segments
raster::plot(segms$segments.id %% 8,asp=1, main="Segments")
```
### 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 :
* 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
+ s: crown surface
+ v: crown volume
+ sp: crown surface inside plot
+ vp: crown volume inside plot
```{r include=TRUE, fig.width = 6.5, fig.height = 6}
## tree extraction only inside plot mask for subsequent comparison
trees <- lidaRtRee::treeExtraction(segms$filled.dem, segms$local.maxima, segms$segments.id, plotMask)
head(trees)
# convert segments from raster to polygons
v.segments <- raster::rasterToPolygons(segms[[2]], dissolve=T)
#
# display initial image
raster::plot(chm,asp=1, col=gray(seq(0,1,1/255)), main ="Canopy Height Model 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)
```
## Detection evaluation
### Tree matching
To assess detection accuracy, reference (field) trees should be linked to detected trees. Despite the possibility of error, automated matching has the advantage of making the comparison of results from different algorithms and settings reproducible and objective. The algorithm presented below is based on the 3D distance between detected treetops and inventory positions and heights [@Monnet2014].
```{r include=TRUE, fig.width = 6.5, fig.height = 6}
# match detected treetops with field trees based on relative distance of apices
matched <- lidaRtRee::treeMatching(treeinventorychablais3[selec,c("x","y","h")],trees[,c("x","y","h")])
# display matching results
lidaRtRee::plot2Dmatched(treeinventorychablais3[selec,c("x","y","h")], trees[,c("x","y","h")], matched, chm, v.plotMask)
```
### Detection accuracy
```{r include=FALSE}
# height histogram of detections
detection.stats <- lidaRtRee::histDetection(treeinventorychablais3[selec,c("x","y","h")], trees[,c("x","y","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.
```{r include=TRUE, fig.width = 6.5, fig.height = 6}
# height histogram of detections
detection.stats <- lidaRtRee::histDetection(treeinventorychablais3[selec,c("x","y","h")], trees[,c("x","y","h")], matched)
```
### Height estimation accuracy
```{r include=FALSE}
heightReg <- lidaRtRee::heightRegression(treeinventorychablais3[selec,c("x","y","h")], trees[,c("x","y","h")], matched, species=treeinventorychablais3$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.
```{r include=TRUE, fig.width = 6.5, fig.height = 6}
# linear regression between reference height and estimated height
heightReg <- lidaRtRee::heightRegression(treeinventorychablais3[selec,c("x","y","h")], trees[,c("x","y","h")], matched, species=treeinventorychablais3$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 labelled with the id of the segment they belong to. A list of las objects corresponding to the segments is then created.
```{r include=TRUE, warning=FALSE, message=FALSE}
# normalize point cloud
lasn <- lidR::lasnormalize(laschablais3, method="delaunay", copy=T)
## add segment id in LAS object
lasn@data$seg.id <- raster::extract(segms[["segments.id"]], lasn@data[,1:2])
# split las object by segment id
lasl <- split(lasn@data,lasn@data$seg.id)
# convert list of data.frames to list of las objects
lasl <- lapply(lasl, function(x){LAS(x,lasn@header)})
```
### Metrics computation
Basic point cloud metrics are computed in each segment (maximum, mean, standard deviation of heights, mean and standard deviation of intensity).
```{r include=TRUE, eval=TRUE}
# compute basic las metrics in each segment
metrics <- as.data.frame(matrix(unlist(lapply(lasl,function(x) {unlist(lidR::lasmetrics(x,list(max(Z),mean(Z), sd(Z), mean(Intensity),sd(Intensity))))})),ncol=5,byrow=T))
row.names(metrics) <- names(lasl)
names(metrics) <- c("maxZ", "meanZ","sdZ","meanI","sdI")
metrics$seg.id <- row.names(metrics)
head(metrics)
```
### Merge with reference and detected trees
Computed metrics are merged with reference and detected tree, thanks to the segment id.
```{r 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")])
# create new data.frame by merging metrics and inventoried trees based on segment id
tree.metrics <- base::merge(treeinventorychablais3,metrics)
# remove non-tree segment
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)
```
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 include=TRUE, fig.width = 6.5, fig.height = 6}
# create raster of segment' mean intensity
r.mean.intensity.segm <- segms[["segments.id"]]
# match segment id with id in metrics data.frame
dummy <- match(raster::values(r.mean.intensity.segm), trees$id)
# replace segment id with corresponding mean intensity
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")
# display tree inventory
lidaRtRee::plotTreeInventory(treeinventorychablais3[,c("x","y")], treeinventorychablais3$h, species=as.character(treeinventorychablais3$s), add=T)
```
### Exploratory analysis
A boxplot of mean intensity in segments per species shows that mean intensity distribution is quite different between species. The analysis can be run by considering only the highest inventoried trees in each segment, as smaller trees are likely to be suppressed and have smaller contribution to the point cloud.
```{r 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),]
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),]} else {
i <- i+1
}
}
tree.metrics.h$s <- factor(tree.metrics.h$s)
```
```{r include=TRUE, fig.width = 12, fig.height = 6}
par(mfrow=c(1,2))
boxplot(meanI~s,data=tree.metrics[,c("s","maxZ", "meanZ", "sdZ", "meanI", "sdI")], ylab="Mean intensity in segment", xlab="Specie", main="All inventoried trees")
boxplot(meanI~s,data=tree.metrics.h, ylab="Mean intensity in segment", xlab="Specie", main="Highest inventoried tree in segment")
```
A linear discriminant analysis shows that main species can be discriminated. In this plot where deciduous trees are in the lower layer, height variables are the ones that best separate deciduous from coniferous. Intensity values somehow allow to distinguish between fir and spruce.
```{r include=TRUE}
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)
```
## Display point cloud
The point cloud can displayed colored by segment, with poles at the location of inventoried trees.
```{r, webgl=TRUE}
# remove points not in tree segments
selecp <- which(lasn@data$seg.id!=0 )
# plot point cloud
rgl::plot3d(lasn@data[selecp,1:3],col=lasn@data$seg.id[selecp]%%10 +1,aspect=FALSE)
#
# add inventoried trees
treeinventorychablais3$z <- 0
for (i in 1:nrow(treeinventorychablais3))
{
rgl::lines3d(rbind(treeinventorychablais3$x[i],treeinventorychablais3$x[i]),rbind(treeinventorychablais3$y[i],treeinventorychablais3$y[i]),rbind(treeinventorychablais3$z[i],treeinventorychablais3$z[i]+treeinventorychablais3$h[i]))
}
```
## References
```{r, include=FALSE, eval=FALSE}
### Inventoried trees
# Using package rLiDAR, a 3d view of the field inventory can be displayed
# shape different between coniferous / deciduous
rgl::rgl.open()
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])) {next}
if (!is.element(as.character(treeinventorychablais3$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,resolution="high",X=treeinventorychablais3$x[i],Y=treeinventorychablais3$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,resolution="high",X=treeinventorychablais3$x[i],Y=treeinventorychablais3$y[i],mesh=F,stemcolor="burlywood4",crowncolor="darkgreen")
}
}
```
```{r include=FALSE, eval=FALSE}
#
# virtual trees from detection
# shape different based on mean lidar intensity value (threshold 55)
library(rLiDAR)
rgl::rgl.open()
rgl::rgl.bg(color="white")
for (i in 1:nrow(trees))
{
if (trees$meanI[i]>55) {
rLiDAR::LiDARForestStand(crownshape="halfellipsoid",CL=0.6*trees$h[i],CW=sqrt(4*trees$s[i]/pi),HCB=0.4*trees$h[i],dbh=trees$h[i]/50,resolution="high",X=trees$x[i],Y=trees$y[i],mesh=F)
} else {
rLiDAR::LiDARForestStand(crownshape="cone",CL=0.5*trees$h[i],CW=sqrt(4*trees$s[i]/pi),HCB=0.5*trees$h[i],dbh=trees$h[i]/50,resolution="high",X=trees$x[i],Y=trees$y[i],mesh=F,stemcolor="burlywood4",crowncolor="darkgreen")
}
}
```
\ No newline at end of file