tree.detection.Rmd 29.2 KB
Newer Older
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
1
---
2
title: "R workflow for tree segmentation from ALS data"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
3
4
5
6
author: "Jean-Matthieu Monnet"
date: "`r Sys.Date()`"
output:
  html_document: default
7
8
  pdf_document: default
papersize: a4
9
bibliography: "./bib/bibliography.bib"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
10
11
12
13
14
15
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# Set so that long lines in R will be wrapped:
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=80),tidy=TRUE)
16
17
18
19
20
knitr::opts_chunk$set(fig.align = "center")
# for display of rgl in html
knitr::knit_hooks$set(webgl = rgl::hook_webgl)
# output to html
html <- TRUE
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
21
22
23
```

---
24
The code below presents a tree segmentation workflow from Airborne Laser Scanning (lidar remote sensing) data. The workflow is based on functions from R packages `lidaRtRee` and `lidR`, and it includes the following steps:
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
25

26
27
28
* treetop detection, 
+ crown segmentation,
+ accuracy assessment with field inventory, 
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
29
30
+ species classification

31
Steps 1 and 3 are documented in [@Monnet10; @Monnet11c]. The detection performance of this algorithm was evaluated in a benchmark [@Eysn15].
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
32

33
Licence: GNU GPLv3 / [source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/tree.detection.Rmd) 
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
34
35
36
37
38
39
40
41
42
43
44
45
46

## 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 loadTreeInventory, 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 prepareTreeInventory, eval=FALSE}
# import field inventory
47
fichier <- "chablais3_listeR.csv"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
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` diameter at breast height (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 displaytreeinventory, echo=FALSE}
head(treeinventorychablais3, n=3L)
```

Function `plotTreeInventory` is designed to plot forest inventory data.

```{r plotTreeInventory, include = TRUE, out.width = '50%', fig.dim=c(5.5, 5.5)}
# display inventoried trees
lidaRtRee::plotTreeInventory(
  treeinventorychablais3[, c("x", "y")],
  treeinventorychablais3$h,
  species = as.character(treeinventorychablais3$s)
)
```
77
78
79
80
81
82
83
84
85
86
87
88
89
90
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"]
library(ggplot2)
ggplot(treeinventorychablais3, aes(x = x, y = y, group = s)) +
  geom_point(aes(color = s, size = d)) +
  coord_sf(datum = 2154) +
  scale_color_manual(values = plot.species) +
  scale_radius(name="Diameter") + 
  geom_text(aes(label=n, size=20), hjust=0, vjust=1) +
  labs(color = "Species") # titre de la légende
```
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
91
92
93
94
95
96
97
98
99
100
101
102

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 roi, include = TRUE}
# buffer to apply around ROI (meters)
ROI.buffer <- 10
# ROI limits
ROI.range <- data.frame(round(apply(treeinventorychablais3[,c("x","y")],2,range)))
```

### Airborne Laser Scanning data

103
In this tutorial, ALS data available in the `lidaRtRee`` package is used. 
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
104

105
```{r loadALS, include = TRUE, message = FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
106
107
# load data in package lidaRtRee (default)
data(laschablais3, package="lidaRtRee")
108
laschablais3
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
109
110
```

111
112
113
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}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
# directory for laz files
lazdir <- "/directory_of_las_files/"
# build catalog of files
cata <- lidR::readLAScatalog(lazdir)
# set coordinate system
sp::proj4string(cata) <- sp::CRS(SRS_string = "EPSG:2154")
# set sensor type
lidR::sensor(cata) <- "ALS"
# extract points in ROI plus additional buffer
laschablais3 <- 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)
# save as rda for easier access:
# save(laschablais3, file="laschablais3.rda", compress = "bzip2")
```

## 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 subtracting the DTM to the DSM.

```{r computeDEMs, include=TRUE}
# 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)
# surface model
dsm <- lidaRtRee::points2DSM(laschablais3,res=0.5, dtm@extent@xmin,
                             dtm@extent@xmax, dtm@extent@ymin, dtm@extent@ymax)
# canopy height model
chm <- dsm - dtm
```

```{r plotDEMs, echo=FALSE, fig.width = 12, fig.height = 4.5, out.width='100%'}
151
par(mfrow = c(1, 3))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
152
# display DTM
153
raster::plot(dtm, main = "DTM")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
154
# display DSM
155
raster::plot(dsm, main = "DSM")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
156
# display CHM
157
raster::plot(chm, main = "CHM")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
```

### 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 computeMask, 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.

180
```{r plotPlot, include = TRUE, out.width = '70%', fig.dim=c(6.5, 4.5), warnings=FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
181
# display CHM
182
raster::plot(chm, col=gray(seq(0,1,1/255)),
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
             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 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)
# 
par(mfrow=c(1,3)) 
# display pre-processed chm
207
raster::plot(segms$smoothed.dem, main="Pre-processed CHM")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
208
# display selected local maxima
209
raster::plot(segms$local.maxima, main="Selected local maxima")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
210
211
212
# display segments, except ground segment
dummy <- segms$segments.id
dummy[dummy==0] <- NA
213
raster::plot(dummy %% 8, main="Segments (random colors)", col=rainbow(8), legend=FALSE)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
```

### 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 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)
head(trees, n=3L)
# convert segments from raster to polygons
v.segments <- raster::rasterToPolygons(segms[[2]], dissolve=T)
#
# display initial image
240
raster::plot(chm, col=gray(seq(0,1,1/255)), main ="CHM and detected positions")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
241
242
243
244
245
# 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
246
graphics::points(trees$x, trees$y, col="blue", cex=trees$h/20, pch = 2)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
247
248
249
250
251
```

## Detection evaluation
### Tree matching

252
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 [@Monnet10]. 
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
253

254
```{r plotMached, include=TRUE, out.width = '70%', fig.dim=c(6.5, 4.5)}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
# match detected treetops with field trees based on relative distance of apices
matched <- lidaRtRee::treeMatching(treeinventorychablais3[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)
```

### Detection accuracy

```{r detectionStats, include=FALSE}
# height histogram of detections
detection.stats <- lidaRtRee::histDetection(treeinventorychablais3[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.

```{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")],
                                            cbind(trees@coords, trees$h), matched)
```

### Height estimation accuracy

```{r heighRegression, include=FALSE}
heightReg <- lidaRtRee::heightRegression(treeinventorychablais3[selec,c("x","y","h")],
                                         cbind(trees@coords, trees$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 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")],
                                         cbind(trees@coords, trees$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 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())
# 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
309
310
311
312
lasl <- lapply(lasl, function(x){lidR::LAS(x,lasn@header)})
# set coordinate system
dummy <- sp::CRS(SRS_string = "EPSG:2154")
for (i in 1:length(lasl)) {sp::proj4string(lasl[[i]]) <- dummy}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
```

### Metrics computation

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(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)
head(metrics, n=3L)
```

### Merge with reference and detected trees

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")])
# 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.
.
345
```{r Intensity, include=TRUE, out.width = '70%', fig.dim=c(6.5, 4.5)}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
# 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 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 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),]
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 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")],
        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",
        xlab="Specie", main="Highest inventoried tree in segment", las=2,
         varwidth = TRUE)
```

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)
plot(afd)
```

## Display point cloud

The point cloud can be displayed colored by segment, with poles at the location of inventoried trees.

401
402
```{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
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
403
404
405
406
407
408
409
410
411
412
413
# 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$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)
#
# add inventoried trees
treeinventorychablais3$z <- 0
for (i in 1:nrow(treeinventorychablais3))
{
414
415
416
417
  rgl::lines3d(rbind(treeinventorychablais3$x[i]-974300,
                     treeinventorychablais3$x[i]-974300),
               rbind(treeinventorychablais3$y[i]-6581600,
                     treeinventorychablais3$y[i]-6581600),
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
               rbind(treeinventorychablais3$z[i],
                     treeinventorychablais3$z[i]+treeinventorychablais3$h[i]))
}
```

```{r plotTreeModel1, 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 plotTreeModel2, 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"
    )
  }
}
```

## Batch processing

The following code exemplifies how to process numerous LAS files and extract trees for the whole area with parallel processing. The processing runs faster if data is provided as chunks to the segmentation algorithm, and results are then aggregated, rather than running on the full coverage of the data. 
In order to avoid border effects, chunks are provided to the algorithm as overlapping tiles. Results are cropped to prevent the same tree from appearing twice in the final results. Tile size and buffer size are important parameters :

- tile size is a trade-off between the number of chunks to process and the amount of RAM required to process a single tile ;
- buffer size is a trade-off between redundant processing of the overlap area, and assuring that a tree is which treetop is located at the border of a tile has its full crown within the buffer size.

Tiles can be processed with parallel computing, within limits of the cluster's RAM and number of cores.

The steps for processing a batch of las/laz files are :

- build catalog of files
- provide tiling parameters
- provide segmentation parameters
- provide output parameters
- set cluster options for parallel processing
- compute the X and Y coordinates of tiles
- parallelize the processing with the buil-in package `parallel`, by sending to the clusters the coordinates of tiles to process, and a function with the instructions to proceed :
    - load tile corresponding to the coordinates
    - compute CHM
    - segment and extract trees
- aggregate list of results

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)}
rm(list=ls())
# BUILD CATALOG OF FILES
# folder containing the files
538
lazdir <- "./data/forest.structure.metrics"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
539
540
541
542
543
544
545
546
547
548
# build catalog
cata <- lidR::catalog(lazdir)
# set coordinate system
sp::proj4string(cata) <- sp::CRS(SRS_string = "EPSG:2154")
# set sensor type
lidR::sensor(cata) <- "ALS"
#
# BATCH PROCESSING PARAMETERS
# tile size to split area into chunks to process
# trade-off between RAM capacity VS total number of tiles to process
549
tile.size <- 70 # here 70 for example purpose with small area
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
550
# buffer size: around each tile to avoid border effects on segmentation results
551
552
# trade-off between making sure a whole tree crown is processed in case its top 
# is on the border VS duplicate processing
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
buffer.size <- 10 # 5 m is minimum, 10 is probably better depending on tree size
#
# TREE SEGMENTATION PARAMETERS
# set raster resolution
resolution <- 1
#
# OUTPUT PARAMETERS
# option to vectorize tree crowns (set to FALSE if too slow)
vectorize.trees <- TRUE
# output canopy height models ? (set to FALSE if too much RAM used)
output.chm <- TRUE
# save chms on disk
save.chm <- FALSE
#
# CLUSTER PARAMETERS
# working directory
wd <- getwd()
# run with two cores
571
572
clust <- parallel::makeCluster(getOption("cl.cores", 2))
#, outfile = "/home/jean-matthieu/Bureau/output.txt")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
# export global variables to cluster because they will be called inside the function
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
  )#[1:2]
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
)#[1:2]
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))
# 
# PARALLEL PROCESSING
# function takes coordinates of tile as arguments
resultats <- parallel::mcmapply(
  # i and j are the coordinated of the low left corner of the tile to process
  FUN = function(i, j)
  {
    # set working directory
    setwd(wd)
    # initialize output
    output <- list()
    output$name <- paste0(i, "_", j)
    # extract tile plus buffer
    las.tile <- lidR::clip_rectangle(
      cata,
      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)
    {
      # normalization if required
      # 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)
      # check all chm is not NA
      if (!all(is.na(raster::values(chm))))
      {
        # output chm if asked for
625
626
        if (output.chm) output$chm <- raster::crop(chm, raster::extent(i, i+tile.size,
                                                                       j, j+tile.size))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
627
        # save on disk
628
629
630
631
632
        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)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
633
634
635
636
637
638
639
640
        #
        # tree detection
        segms <- lidaRtRee::treeSegmentation(chm)
        # tree extraction
        trees <- lidaRtRee::treeExtraction(segms$filled.dem,
                                           segms$local.maxima,
                                           segms$segments.id)
        # remove trees in buffer area
641
642
        trees <- trees[trees$x >= i & trees$x < i + tile.size &
                         trees$y >= j & trees$y < j + tile.size, ]
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
643
644
645
646
647
648
649
650
651
652
653
654
        # 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)
        {
          # vectorize
          v.trees <- 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"
          # add attributes
          # errors when using sp::merge so using sp::over even if it is probably slower
655
656
          # merge(v.trees@data, trees, all.x = TRUE)
          v.trees@data <- cbind(v.trees@data, sp::over(v.trees, trees))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
          # save in list
          output$v.trees <- v.trees
        }
        # save trees in list
        output$trees <- trees
      } # end of raster is not all NAs check
    } # end of nrow LAS check
    output
  }, x, y, SIMPLIFY = FALSE) # function applied to the lists of coordinates (x and y)
parallel::stopCluster(cl = clust)
#
# RESULTS AGGREGATION
# extract results from nested list into separate lists and then bind data
id <- unlist(lapply(resultats, function(x) x[["name"]]))
#
# trees
trees <- lapply(resultats, function(x) x[["trees"]])
# remove NULL elements
trees[sapply(trees, is.null)] <- NULL
# bind remaining elements
trees <- do.call(rbind, trees)
# 
# chm
if (output.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)
  names(chm) <- id
}
# v.trees
if (vectorize.trees)
{
  v.trees <- lapply(resultats, function(x) x[["v.trees"]])
  # remove NULL elements
  v.trees[sapply(v.trees, is.null)] <- NULL
  v.trees <- do.call(rbind, v.trees)
695
696
  # 1-pixel overlapping in v.trees might be present because image segmentation
  # is not fully identical in overlap areas of adjacent tiles.
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
697
698
699
700
701
}
```

The following image displays the results for the whole area.

702
```{r batch.plot, include=TRUE, out.width = '90%', fig.dim=c(8.5, 5.5)}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
703
704
705
# threshold outsiders in chm
chm.all[chm.all > 40] <- 40
chm.all[chm.all < 0] <- 0
706
707
708
709
710
# display chm
raster::plot(chm.all,
             main ="Canopy Height Model and segmented trees")
# display segments border
sp::plot(v.trees, border = "white", add = T)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
711
# add trees
712
sp::plot(trees, cex = trees$h/40, add = TRUE, pch = 2)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
713
714
715
716
717
718
719
720
721
722
723
724
725
726
```

The following lines save outputs to files.

```{r batch.export, eval = FALSE}
# merged chm
raster::writeRaster(chm.all, file= "chm.tif")
#trees
raster::shapefile(trees, file="trees.points.shp")
#vectorized trees
if (vectorize.trees) raster::shapefile(v.trees, file="trees.polygons.shp")
```

## References