forest.structure.metrics.Rmd 30 KB
Newer Older
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
1
---
2
title: "R workflow for forest structure metrics computation from ALS data"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
3
author: "Jean-Matthieu Monnet, with contributions from B. Reineking and A. Glad"
4
date: "`r Sys.Date()`"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
5
output:
6
  html_document: default
7
  pdf_document: default
8
papersize: a4
9
bibliography: "../bib/bibliography.bib"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
10
11
12
---

```{r setup, include=FALSE}
13
14
15
16
# erase all
cat("\014")
rm(list = ls())
# knit options
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
17
18
knitr::opts_chunk$set(echo = TRUE)
# Set so that long lines in R will be wrapped:
19
knitr::opts_chunk$set(tidy.opts = list(width.cutoff = 80), tidy = TRUE)
20
knitr::opts_chunk$set(fig.align = "center")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
21
22
23
```

---
24
Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/forest.structure.metrics.Rmd)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
25

26
27
Many thanks to Pascal Obstétar for checking code and improvement suggestions.

28
The R code below presents a forest structure metrics computation workflow from Airborne Laser Scanning (ALS) data. Workflow is based on functions from `R` packages [lidaRtRee](https://cran.r-project.org/package=lidaRtRee) (tested with version `r packageVersion("lidaRtRee")`) and [lidR](https://cran.r-project.org/package=lidR) (tested with version `r packageVersion("lidR")`). Package `vegan` is 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.
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
29
30
31
32
Different types of metrics are computed:  

* 1D height metrics
+ 2D metrics of the canopy height model (CHM)
33
34
+ tree segmentation metrics (see also [tree segmentation tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/tree.detection.Rmd))
+ forest gaps and edges metrics (see also [gaps and edges detection tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/gaps.edges.detection.Rmd))
35
36
37

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.

38
The workflow processes normalized point clouds provided as las/laz tiles of rectangular extent. Parallelization is used for faster processing, packages `future` and `future.apply`  are used. A buffer is loaded around each tile to prevent border effects in tree segmentation and CHM processing.
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
39
40
41

## Parameters

42
Set the number of cores to use for parallel computing.
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
43
```{r parallelSetup, include = TRUE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
44
45
# create parallel frontend, specify to use two parallel sessions
future::plan("multisession", workers = 2L)
46
47
# remove warning when using random numbers in parallel sessions
options(future.rng.onMisuse = "ignore")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
48
49
50
51
52
53
54
```

Numerous parameters have to be set for processing.
```{r parameters, include = TRUE}
# output metrics map resolution (m)
resolution <- 10
# canopy height model resolution (m)
55
res_chm <- 0.5
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
56
57
# buffer size (m) for tile processing (20 m is better for gaps metrics, 10 m is enough for
# tree metrics)
58
buffer_size <- 20
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
59
60
# height threshold (m) for high points removal (points above this threshold are considered
# as outliers)
61
points_max_h <- 60
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
62
# points classes to retain for analysis (vegetation + ground)
63
class_points <- c(0, 1, 2, 3, 4, 5)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
64
# ground class
65
class_ground <- 2
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
66
# Gaussian filter sigma values for multi-scale smoothing
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
67
68
sigma_l <- list(0, 1, 2, 4, 8, 16)
# fonction to computed raster statistics from multi-scale smoothing
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
smoothed_raster_stats <- function(x) {
  data.frame(
    CHM0_sd = sd(x$smoothed_image_0),
    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$smoothed_image_0),
    CHM_PercInf0_5 = sum(x$smoothed_image_0 < 0.5) * mf,
    CHM_PercInf1 = sum(x$smoothed_image_0 < 1) * mf,
    CHM_PercSup5 = sum(x$smoothed_image_0 > 5) * mf,
    CHM_PercSup10 = sum(x$smoothed_image_0 > 10) * mf,
    CHM_PercSup20 = sum(x$smoothed_image_0 > 20) * mf,
    CHM_Perc1_5 = (sum(x$smoothed_image_0 < 5) - sum(x$smoothed_image_0 < 1)) * mf,
    CHM_Perc2_5 = (sum(x$smoothed_image_0 < 5) - sum(x$smoothed_image_0 < 2)) * mf
  )
}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
87
# height breaks for penetration ratio and density
88
breaks_h <- c(-Inf, 0, 0.5, 1, 2, 5, 10, 20, 30, 60, Inf)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
89
90
91
# percentiles of height distribution
percent <- c(0.10, 0.25, 0.5, 0.75, 0.9)
# surface breaks for gap size (m2)
92
breaks_gap_surface <- c(4, 16, 64, 256, 1024, 4096, Inf)
93
#
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
94
# gap surface names
95
96
n_breaks_gap <- gsub("-", "", paste0("G_s", breaks_gap_surface[c(-length(breaks_gap_surface))], "to",
                                     breaks_gap_surface[c(-1)]
97
))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
98
# height bin names
99
n_breaks_h <- gsub("-", "", paste0("nb_H", breaks_h[c(-length(breaks_h))], "to", breaks_h[c(-1)]
100
))
101
# 
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
102
103
104
```
The first step is to create a catalog of LAS files (should be normalized, non-overlapping rectangular tiles). Preferably, tiles should be aligned on a multiple of resolution, and points should not lie on the northern or eastern border when such borders are common with adjacent tiles.

105
```{r lascatalog, include = TRUE, fig.dim = c(3.5, 2.5), out.width='40%', warning=FALSE}
106
107
# create catalog of LAS files
cata <- lidR::readALSLAScatalog("../data/forest.structure.metrics")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
108
# set coordinate system
109
lidR::projection(cata) <- 2154
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
110
111
# disable display of catalog processing
lidR::opt_progress(cata) <- FALSE
112
113
# option to read only xyzc attributes (coordinates, intensity, echo order and classification) from height files (saves memory)
lidR::opt_select(cata) <- "xyzirnc"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
114
# display
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
115
lidR::plot(cata, main = "Bounding box of las files")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
116
117
118
119
120
121
122
123
124
125
126
```

## Workflow exemplified on one tile

### Load point cloud
The tile is loaded, plus a buffer on adjacent tiles. Buffer points have a column "buffer" filled with 1. 

```{r loadPointCloud, include = TRUE, message=FALSE, warnings=FALSE, fig.dim = c(4.5, 3.5), out.width='60%', }
# tile to process
i <- 1
# tile extent
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
127
b_box <- sf::st_bbox(cata[i,])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
128
# load tile extent plus buffer
129
a <- lidR::clip_rectangle(
130
131
  cata, b_box[1] - buffer_size, b_box[2] - buffer_size, b_box[3] + buffer_size,
  b_box[4] + buffer_size
132
)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
133
# add 'buffer' flag to points in buffer with TRUE value in this new attribute
134
a <- lidR::add_attribute(
135
  a, a$X < b_box[1] | a$Y < b_box[2] | a$X >= b_box[3] | a$Y >= b_box[4],
136
137
  "buffer"
)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
138
139
140
141
142
143
a
```

Only points from desired classes are retained. In case some mis-classified, high or low points remain in the data set, the point cloud is filtered and negative heights are replaced by 0.
```{r filterPointCloud, include = TRUE, message=FALSE, warnings=FALSE}
# remove unwanted point classes, and points higher than height threshold
144
a <- lidR::filter_poi(a, is.element(Classification, class_points) & Z <= points_max_h)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
145
# set negative heights to 0
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
146
a$Z[a$Z < 0] <- 0
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
147
148
149
150
151
152
153
154
#
summary(a@data)
```

### Canopy height model

The next step is to compute the canopy height model (CHM). It will be used to derive:

155
156
* 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
157
+ gaps and edges metrics from gap segmentation
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
158
159
160

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.

161
```{r computeCHM, include = TRUE, fig.width = 6, fig.height = 4.3,  message=FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
162
# compute chm
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
163
chm <- lidR::rasterize_canopy(a, res = res_chm, algorithm = lidR::p2r(), pkg = "terra")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
164
# replace NA, low and high values
165
chm[is.na(chm) | chm < 0 | chm > points_max_h] <- 0
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
166
# display CHM
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
167
terra::plot(chm, asp = 1, main = "Canopy height model")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
168
169
170
171
172
173
```

## Metrics computation

### 2D CHM metrics

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
174
The CHM is smoothed with a Gaussian filter with different `sigma` values. Smoothed results are stored in a list and then integrated ito a single  raster
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
175

176
```{r 2dchmMetrics, include = TRUE, fig.width = 12, fig.height = 6.5, message=FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
177
# for each value in list of sigma, apply filtering to chm and store result in list
178
179
st <- lapply(sigma_l, FUN = function(x) {
  lidaRtRee::dem_filtering(chm, nl_filter = "Closing", nl_size = 5, sigmap = x)$smoothed_image
180
})
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
181
# convert to raster
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
182
183
184
st <- terra::rast(st)
# modify layer names
names(st) <- paste0(names(st), "_", unlist(sigma_l))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
185
# display
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
186
terra::plot(st, range = range(terra::values(st[[1]])))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
187
188
```

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
189
The raster is then converted to points. Points outside the tile extent are removed and summary statistics (mean and standard deviation) of smoothed CHM heights are computed for each pixel at the output metrics resolution. Canopy covers in different height layers are also computed for the initial CHM after applying a non-linear filter designed to fill holes. 
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
190
191
192

```{r 2dmetrics2points, include = TRUE, message=FALSE, warning=FALSE}
# crop to bbox
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
193
st <- terra::crop(st, terra::ext(b_box[1], b_box[3], b_box[2], b_box[4]))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
194
# multiplying factor to compute proportion
195
mf <- 100 / (resolution / res_chm)^2
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
196

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
197
# compute raster metrics
198
199
metrics_2dchm <- lidaRtRee::raster_metrics(st,
                                           res = resolution,
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
200
                                           fun = smoothed_raster_stats,
201
                                           output = "raster"
202
)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
203
#
204
metrics_2dchm
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
205
206
207
208
```

Some outputs are displayed hereafter. On the first line are the standard deviation of CHM smoothed with different `sigma` values. On the second line are the CHM mean and percentages of CHM values below 1 m and above 10 m.

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
209
```{r 2dmetricsOutputs, include = TRUE, fig.width = 12, fig.height = 7.6, warning=FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
210
terra::plot(metrics_2dchm[[c(
211
212
  "CHM0_sd", "CHM2_sd", "CHM8_sd", "CHM_mean",
  "CHM_PercInf1", "CHM_PercSup10"
213
)]])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
214
215
216
217
```

### Gaps and edges metrics

218
Gaps are computed with the function `gap_detection`. 
219
```{r gapMetrics, include = TRUE, fig.width = 12, fig.height = 2.9, warning=FALSE, message=FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
220
# compute gaps
221
222
223
224
gaps <- lidaRtRee::gap_detection(chm,
                                 ratio = 2, gap_max_height = 1,
                                 min_gap_surface = min(breaks_gap_surface), closing_height_bin = 1,
                                 nl_filter = "Median", nl_size = 3, gap_reconstruct = TRUE
225
)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
226
# display maps
227
par(mfrow = c(1, 3))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
228
229
terra::plot(chm, main = "CHM")
terra::plot(log(gaps$gap_surface), main = "Log of gap surface")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
230
# display gaps
231
dummy <- gaps$gap_id
232
dummy[dummy == 0] <- NA
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
233
terra::plot(dummy %% 8, asp = 1, main = "Gap id (random colors)", col = rainbow(8))#, legend = FALSE)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
234
235
236
237
238
```
Summary statistics are then computed: pixel surface occupied by gaps in different size classes. Results are first cropped to the tile extent.

```{r gapStats, include = TRUE, warning=FALSE, message=FALSE}
# crop results to tile size
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
239
gaps_surface <- terra::crop(gaps$gap_surface, terra::ext(
240
241
  b_box[1], b_box[3],
  b_box[2], b_box[4]
242
))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
243
# compute gap statistics at final metrics resolution, in case gaps exist
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
244
if (!all(is.na(terra::values(gaps_surface)))) {
245
246
247
  metrics_gaps <- lidaRtRee::raster_metrics(gaps_surface,
                                            res = resolution,
                                            fun = function(x) {
248
                                              hist(x$gap_surface,
249
250
251
252
                                                   breaks = breaks_gap_surface,
                                                   plot = F
                                              )$counts * (res_chm / resolution)^2
                                            }, output = "raster"
253
  )
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
254
  # compute total gap proportion
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
255
  metrics_gaps$sum <- terra::tapp(metrics_gaps, rep(1, length(names(metrics_gaps))),
256
                                         fun = sum
257
  )
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
258
  # if gaps are present
259
  names(metrics_gaps) <- c(n_breaks_gap, paste("G_s", min(breaks_gap_surface), "toInf", sep = ""))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
260
261
  # replace possible NA values by 0 (NA values are present in lines of the target raster
  # where no values >0 are present in the input raster)
262
  metrics_gaps[is.na(metrics_gaps)] <- 0
263
} else {
264
  metrics_gaps <- NULL
265
}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
266
#
267
metrics_gaps
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
268
269
```

270
Edges are detected  with the function `edge_detection` as the outer envelope of the previously delineated gaps. 
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
271

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
272
```{r edgeDetection, include = TRUE, fig.width = 12, fig.height = 3.4, warning=FALSE, message=FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
273
# Perform edge detection by erosion
274
edges <- lidaRtRee::edge_detection(gaps$gap_id > 0)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
275
# display maps
276
par(mfrow = c(1, 3))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
277
278
279
terra::plot(chm, main = "CHM")
terra::plot(gaps$gap_id > 0, main = "Gaps", legend = FALSE)
terra::plot(edges, main = "Edges", legend = FALSE)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
280
281
282
283
284
```
The percentage of surface occupied by edges is then computed as summary statistic. Results are first cropped to the tile extent.

```{r edgeStats, include = TRUE, warning=FALSE, message=FALSE}
# crop results to tile size
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
285
edges <- terra::crop(
286
  edges,
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
287
  terra::ext(b_box[1], b_box[3], b_box[2], b_box[4])
288
)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
289
# compute gap statistics at final metrics resolution, in case gaps exist
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
290
if (!all(is.na(terra::values(edges)))) {
291
292
293
  metrics_edges <- lidaRtRee::raster_metrics(edges,
                                             res = resolution,
                                             fun = function(x) {
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
294
                                               sum(x$lyr.1) * (res_chm / resolution)^2
295
                                             }, output = "raster"
296
  )
297
  names(metrics_edges) <- "edges.proportion"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
298
299
  # replace possible NA values by 0 (NA values are present in lines of the target raster
  # where no values >0 are present in the input raster)
300
  metrics_edges[is.na(metrics_edges)] <- 0
301
} else {
302
  metrics_edges <- NULL
303
}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
304
#
305
metrics_edges
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
306
307
308
309
```

Some outputs are displayed hereafter. The proportion of surface occupied by gaps of various sizes is depicted as well as the proportion of surface occupied by edges.

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
310
```{r gapDisplay, include = TRUE, fig.width = 12, fig.height = 3.4, warning=FALSE, message=FALSE}
311
par(mfrow = c(1, 3))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
312
313
314
315
# terra::plot(metrics_gaps[["G_s4to16"]], main="Prop of gaps 4-16m2")
terra::plot(metrics_gaps[["G_s64to256"]], main = "Prop of gaps 64-256m2")
terra::plot(metrics_gaps[["G_s4toInf"]], main = "Prop of gaps >4m2")
terra::plot(metrics_edges, main = "Proportion of edges")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
316
317
318
319
320
```

### Tree metrics

Tree tops are detected with the function `treeSegmentation` and then extracted with `treeExtraction`. 
321
```{r treeMetrics, include = TRUE, fig.width = 9, fig.height = 3.1, warning=FALSE, message=FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
322
# tree top detection (default parameters)
323
segms <- lidaRtRee::tree_segmentation(chm, hmin = 5)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
324
# extraction to data.frame
325
trees <- lidaRtRee::tree_extraction(segms)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
326
# display maps
327
par(mfrow = c(1, 2))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
328
terra::plot(chm, main = "CHM and tree tops")
329
points(trees$x, trees$y, pch = 4, cex = 0.3)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
330
# display segments, except ground segment
331
dummy <- segms$segments_id
332
dummy[dummy == 0] <- NA
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
333
terra::plot(dummy %% 8, asp = 1, main = "Segments (random colors)", col = rainbow(8))
334
points(trees$x, trees$y, pch = 4, cex = 0.3)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
335
336
```

337
Results are cropped to the tile extent. Summary statistics from `std_tree_metrics` are then computed for each pixel. Canopy cover in trees and mean canopy height in trees are computed in a second step because they can not be computed from the data of the extracted trees which crowns may span several output pixels. They are calculated from the images obtained during the previous tree segmentation.
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
338
339
340

```{r cropTreeMetrics, include = TRUE, warning=FALSE, message=FALSE}
# remove trees outside of tile
341
trees <- trees[trees$x >= b_box[1] & trees$x < b_box[3] & trees$y >= b_box[2] & trees$y < b_box[4], ]
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
342
# compute raster metrics
343
344
345
346
347
metrics_trees <- lidaRtRee::raster_metrics(trees[, -1],
                                           res = resolution,
                                           fun = function(x) {
                                             lidaRtRee::std_tree_metrics(x, resolution^2 / 10000)
                                           }, output = "raster"
348
)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
349
# remove NAs and NaNs
350
metrics_trees[!is.finite(metrics_trees)] <- 0
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
351
#
352
353
# compute canopy cover in trees and canopy mean height in trees
# in region of interest, because it is not in previous step.
354
r_tree_chm <- segms$filled_dem
355
# set chm to NA in non segment area
356
r_tree_chm[segms$segments_id == 0] <- NA
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
357
# compute raster metrics
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
358
dummy <- lidaRtRee::raster_metrics(terra::crop(
359
  r_tree_chm,
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
360
  terra::ext(
361
362
    b_box[1], b_box[3], b_box[2],
    b_box[4]
363
364
365
366
367
  )
),
res = resolution,
fun = function(x) {
  c(
368
369
    sum(!is.na(x$filled_dem)) / (resolution / res_chm)^2,
    mean(x$filled_dem, na.rm = T)
370
371
372
373
  )
},
output = "raster"
)
374
names(dummy) <- c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
375
#
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
376
377
dummy <- terra::extend(dummy, metrics_trees)
metrics_trees <- c(metrics_trees, dummy)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
378
#
379
metrics_trees
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
380
381
382
```
Some outputs are displayed hereafter. 

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
383
```{r displayTreeMetrics, include = TRUE, fig.width = 12, fig.height = 3.4, warning=FALSE, message=FALSE}
384
par(mfrow = c(1, 3))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
385
terra::plot(metrics_trees$Tree_meanH, main = "Mean (detected) tree height (m)")
386
points(trees$x, trees$y, cex = trees$h / 40)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
387
terra::plot(metrics_trees$TreeSup20_density,
388
             main = "Density of (detected) trees > 20m (/ha)"
389
)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
390
terra::plot(metrics_trees$Tree_meanCrownVolume,
391
             main = "Mean crown volume of detected trees (m3)"
392
)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
393
394
395
396
397
398
399
400
```

### Point cloud (1D) metrics

Buffer points are first removed from the point cloud, then metrics are computed from all points and from vegetation points only.

```{r 1dmetrics, include = TRUE}
# remove buffer points
401
a <- lidR::filter_poi(a, buffer == 0)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
402
403
#
# all points metrics
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
404
metrics_1d <- lidR::pixel_metrics(a, as.list(c(
405
406
407
408
  as.vector(quantile(Z, probs = percent)),
  sum(Classification == 2),
  mean(Intensity[ReturnNumber == 1])
)), res = resolution)
409
names(metrics_1d) <- c(paste("Hp", percent * 100, sep = ""), "nb_ground", "Imean_1stpulse")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
410
411
#
# vegetation-only metrics
412
a <- lidR::filter_poi(a, Classification != class_ground)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
413
dummy <- lidR::pixel_metrics(a, as.list(c(
414
415
416
417
  mean(Z),
  max(Z),
  sd(Z),
  length(Z),
418
  hist(Z, breaks = breaks_h, right = F, plot = F)$counts
419
)), res = resolution)
420
names(dummy) <- c("Hmean", "Hmax", "Hsd", "nb_veg", n_breaks_h)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
421
# merge rasterstacks
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
422
423
dummy <- terra::extend(dummy, metrics_1d)
metrics_1d <- c(metrics_1d, dummy)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
424
425
rm(dummy)
# crop to tile extent
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
426
metrics_1d <- terra::crop(metrics_1d, terra::ext(b_box[1], b_box[3], b_box[2], b_box[4]))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
427
428
429
430
```

Based on those metrics, additional metrics are computed (Simpson index, relative density in height bins and penetration ratio in height bins)

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
431
```{r 1dmetricsAdditional, include = TRUE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
432
# Simpson index
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
433
metrics_1d$Hsimpson <- terra::app(
434
  metrics_1d[[n_breaks_h[c(-1, -length(n_breaks_h))]]],
435
436
437
438
  function(x, ...) {
    vegan::diversity(x, index = "simpson")
  }
)
439
440
# Relative density 
for (i in n_breaks_h[c(-1, -length(n_breaks_h))])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
441
{
442
443
  metrics_1d[[gsub("nb", "relativeDensity", i)]] <-
    metrics_1d[[i]] / (metrics_1d$nb_veg + metrics_1d$nb_ground)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
444
445
446
}
# Penetration ratio
# cumulative sum
447
448
cumu_sum <- metrics_1d[["nb_ground"]]
for (i in n_breaks_h)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
449
{
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
450
  cumu_sum <- c(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics_1d[[i]])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
451
}
452
names(cumu_sum) <- c("nb_ground", n_breaks_h)
453
# compute interception ratio of each layer
454
455
intercep_ratio <- cumu_sum[[-1]]
for (i in 2:dim(cumu_sum)[3])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
456
{
457
  intercep_ratio[[i - 1]] <- 1 - cumu_sum[[i - 1]] / cumu_sum[[i]]
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
458
}
459
names(intercep_ratio) <- gsub("nb", "intercepRatio", names(cumu_sum)[-1])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
460
# merge rasterstacks
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
461
metrics_1d <- c(metrics_1d, intercep_ratio)
462
# rm(cumu_sum, intercep_ratio)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
463
#
464
metrics_1d
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
465
466
467
468
```

Some outputs are displayed hereafter. 

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
469
```{r display1dmetrics, include = TRUE, fig.width = 12, fig.height = 3.4, warning=FALSE, message=FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
470
# display results
471
par(mfrow = c(1, 3))
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
472
473
474
terra::plot(metrics_1d$Imean_1stpulse, main = "Mean intensity (1st return)")
terra::plot(metrics_1d$Hp50, main = "Percentile 50 of heights")
terra::plot(metrics_1d$intercepRatio_H2to5, main = "Interception ratio 2-5 m")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
475
476
477
478
```

## Batch processing

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
479
The following code uses parallel processing to handle multiple files of a catalog, and arranges all metrics in a raster. Code in the "parameters" section has to be run beforehand.
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
480

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
481
```{r batchProcessing, include = TRUE, eval=TRUE, warning=FALSE, message=FALSE, fig.width = 12, fig.height = 3.1}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
482
# processing by tile
483
484
485
486
metrics <- future.apply::future_lapply(
  as.list(1:nrow(cata)),
  FUN = function(i) {
    # tile extent
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
487
    b_box <- sf::st_bbox(cata[i,])
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
    # load tile extent plus buffer
    a <- try(lidR::clip_rectangle(
      cata,
      b_box[1] - buffer_size,
      b_box[2] - buffer_size,
      b_box[3] + buffer_size,
      b_box[4] + buffer_size
    ))
    #
    # 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 <= points_max_h)
    # check number of remaining points
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
511
    if (nrow(a) == 0)
512
513
      return(NULL)
    # set negative heights to 0
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
514
    a$Z[a$Z < 0] <- 0
515
516
    #
    # compute chm
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
517
    chm <- lidR::rasterize_canopy(a, res = res_chm, algorithm = lidR::p2r(), pkg = "terra")
518
519
520
521
522
523
524
525
526
527
528
    # replace NA, low and high values
    chm[is.na(chm) | chm < 0 | chm > points_max_h] <- 0
    #-----------------------
    # compute 2D CHM metrics
    # for each value in list of sigma, apply filtering to chm and store result in list
    st <- lapply(
      sigma_l,
      FUN = function(x) {
        lidaRtRee::dem_filtering(chm,
                                 nl_filter = "Closing",
                                 nl_size = 5,
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
529
                                 sigmap = x)$smoothed_image
530
531
532
      }
    )
    # convert to raster
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
533
534
535
    st <- terra::rast(st)
    # modify layer names
    names(st) <- paste0(names(st), "_", unlist(sigma_l))
536
537
    # crop to bbox
    st <-
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
538
      terra::crop(st, terra::ext(b_box[1], b_box[3], b_box[2], b_box[4]))
539
540
541
542
543
544
545
    # multiplying factor to compute proportion
    mf <- 100 / (resolution / res_chm) ^ 2
    # compute raster metrics
    metrics_2dchm <-
      lidaRtRee::raster_metrics(
        st,
        res = resolution,
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
546
        fun = smoothed_raster_stats,
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
        output = "raster"
      )
    #
    # -------------------
    # compute gap metrics
    gaps <- lidaRtRee::gap_detection(
      chm,
      ratio = 2,
      gap_max_height = 1,
      min_gap_surface = min(breaks_gap_surface),
      closing_height_bin = 1,
      nl_filter = "Median",
      nl_size = 3,
      gap_reconstruct = TRUE
    )
    # crop results to tile size
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
563
564
    gaps_surface <- terra::crop(gaps$gap_surface,
                                 terra::ext(b_box[1], b_box[3], b_box[2],
565
566
                                                b_box[4]))
    # compute gap statistics at final metrics resolution, in case gaps exist
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
567
    if (!all(is.na(terra::values(gaps_surface)))) {
568
569
570
571
      metrics_gaps <- lidaRtRee::raster_metrics(
        gaps_surface,
        res = resolution,
        fun = function(x) {
572
          hist(x$gap_surface,
573
574
575
576
577
578
579
               breaks = breaks_gap_surface,
               plot = F)$counts * (res_chm / resolution) ^
            2
        },
        output = "raster"
      )
      # compute total gap proportion
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
580
581
582
583
      metrics_gaps$sum <- terra::tapp(metrics_gaps,
                                      rep(1, length(names(metrics_gaps))),
                                      fun = sum
      )
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
      # if gaps are present
      names(metrics_gaps) <-
        c(n_breaks_gap, paste("G_s", min(breaks_gap_surface), "toInf", sep = ""))
      # replace possible NA values by 0 (NA values are present in lines of the target raster
      # where no values >0 are present in the input raster)
      metrics_gaps[is.na(metrics_gaps)] <- 0
    } else {
      metrics_gaps <- NULL
    }
    #
    # --------------------
    # compute edge metrics
    # Perform edge detection by erosion
    edges <- lidaRtRee::edge_detection(gaps$gap_id > 0)
    # crop results to tile size
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
599
    edges <- terra::crop(edges, terra::ext(b_box[1], b_box[3],
600
601
                                                b_box[2], b_box[4]))
    # compute edges statistics at final metrics resolution, in case edges exist
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
602
    if (!all(is.na(terra::values(edges)))) {
603
604
605
606
      metrics_edges <- lidaRtRee::raster_metrics(
        edges,
        res = resolution,
        fun = function(x) {
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
607
          sum(x$lyr.1) * (res_chm / resolution) ^ 2
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
        },
        output = "raster"
      )
      names(metrics_edges) <- "edges.proportion"
      # replace possible NA values by 0 (NA values are present in lines of the target raster
      # where no values >0 are present in the input raster)
      metrics_edges[is.na(metrics_edges)] <- 0
    } else {
      metrics_edges <- NULL
    }
    #
    # --------------------
    # compute tree metrics
    # tree top detection (default parameters)
    segms <- lidaRtRee::tree_segmentation(chm, hmin = 5)
    # extraction to data.frame
    trees <-
      lidaRtRee::tree_extraction(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::raster_metrics(
      trees[,-1],
633
634
      res = resolution,
      fun = function(x) {
635
        lidaRtRee::std_tree_metrics(x, resolution ^ 2 / 10000)
636
637
638
      },
      output = "raster"
    )
639
640
641
642
643
    # remove NAs and NaNs
    metrics_trees[!is.finite(metrics_trees)] <- 0
    #
    # extend layer in case there are areas without trees
    metrics_trees <-
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
644
645
646
      terra::extend(metrics_trees, metrics_2dchm)
    # set NA values to 0
    metrics_trees[is.na(metrics_trees)] <- 0
647
648
649
650
651
652
653
    # compute canopy cover in trees and canopy mean height in trees
    # in region of interest, because it is not in previous step.
    r_tree_chm <- segms$filled_dem
    # set chm to NA in non segment area
    r_tree_chm[segms$segments_id == 0] <- NA
    # compute raster metrics
    dummy <- lidaRtRee::raster_metrics(
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
654
      terra::crop(r_tree_chm, terra::ext(b_box[1], b_box[3], b_box[2], b_box[4])),
655
656
657
658
659
660
      res = resolution,
      fun = function(x) {
        c(sum(!is.na(x$filled_dem)) / (resolution / res_chm) ^ 2,
          mean(x$filled_dem, na.rm = T))
      },
      output = "raster"
661
    )
662
663
    names(dummy) <-
      c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
664
665
666
    dummy <- terra::extend(dummy, metrics_2dchm)
    # set NA values to 0
    dummy[is.na(dummy)] <- 0
667
    #
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
668
    metrics_trees <-c(metrics_trees, dummy)
669
    #
670
671
672
673
674
    # -------------------------
    # compute 1D height metrics
    # remove buffer points
    a <- lidR::filter_poi(a, buffer == 0)
    #
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
675
    if (nrow(a) == 0) {
676
677
678
      metrics_1d <- NULL
    } else {
      # all points metrics
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
679
      metrics_1d <- lidR::pixel_metrics(a, as.list(c(
680
681
682
        as.vector(quantile(Z, probs = percent)),
        sum(Classification == 2),
        mean(Intensity[ReturnNumber == 1])
683
      )), res = resolution)
684
685
686
687
      names(metrics_1d) <-
        c(paste("Hp", percent * 100, sep = ""),
          "nb_ground",
          "Imean_1stpulse")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
688
      #
689
690
      # vegetation-only metrics
      a <- lidR::filter_poi(a, Classification != class_ground)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
691
692
      if (nrow(a) != 0) {
        dummy <- lidR::pixel_metrics(a, as.list(c(
693
694
695
696
697
698
699
700
701
702
703
704
705
          mean(Z),
          max(Z),
          sd(Z),
          length(Z),
          hist(
            Z,
            breaks = breaks_h,
            right = F,
            plot = F
          )$counts
        )), res = resolution)
        names(dummy) <-
          c("Hmean", "Hmax", "Hsd", "nb_veg", n_breaks_h)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
706
707
708
        # merge raster
        dummy <- terra::extend(dummy, metrics_1d)
        metrics_1d <- c(metrics_1d, dummy)
709
710
711
712
        rm(dummy)
        # -------------
        # merge metrics
        metrics_1d <-
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
713
714
715
          terra::extend(metrics_1d, metrics_2dchm)
        # metrics_1d[is.na(metrics_1d)] <- 0
        metrics_1d <- terra::crop(metrics_1d, metrics_2dchm)
716
        metrics_gaps <-
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
717
718
719
          terra::extend(metrics_gaps, metrics_2dchm)
        # metrics_gaps[is.na(metrics_gaps)] <- 0
        metrics_gaps <- terra::crop(metrics_gaps, metrics_2dchm)
720
        metrics_edges <-
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
721
722
723
724
          terra::extend(metrics_edges, metrics_2dchm)
        # metrics_edges[is.na(metrics_edges)] <- 0
        metrics_edges <- terra::crop(metrics_edges, metrics_2dchm)
        metrics_trees <- terra::crop(metrics_trees, metrics_2dchm)
725
        #
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
726
        temp <- c(metrics_1d,
727
728
729
730
731
732
733
                                 metrics_2dchm,
                                 metrics_gaps,
                                 metrics_edges,
                                 metrics_trees)
        temp[is.na(temp)] <- 0
      } # end of vegetation-only points presence check
    } # end of buffer-less points presence check
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
734
    return(terra::wrap(temp))
735
736
  }
)
737
# --------------
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
738
# mosaic rasters
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
739
740
741
742
743
# unpack rasters
metrics <- lapply(metrics, terra::rast)
# names_metrics <- names(metrics[[1]])
metrics <- do.call(terra::merge, metrics)
# names(metrics) <- names_metrics
744
# --------------------------
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
745
746
# compute additional metrics
# Simpson index
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
747
metrics$Hsimpson <- terra::app(metrics[[n_breaks_h[c(-1,-length(n_breaks_h))]]],
748
749
750
                                       function(x, ...) {
                                         vegan::diversity(x, index = "simpson")
                                       })
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
751
752
#
# Relative density
753
for (i in n_breaks_h[c(-1,-length(n_breaks_h))])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
754
{
755
756
  metrics[[gsub("nb", "relativeDensity", i)]] <-
    metrics[[i]] / (metrics$nb_veg + metrics$nb_ground)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
757
758
759
}
# Penetration ratio
# compute cumulative sum
760
761
cumu_sum <- metrics[["nb_ground"]]
for (i in n_breaks_h)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
762
{
763
  cumu_sum <-
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
764
    c(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics[[i]])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
765
}
766
names(cumu_sum) <- c("nb_ground", n_breaks_h)
767
# interception ratio
768
769
intercep_ratio <- cumu_sum[[-1]]
for (i in 2:dim(cumu_sum)[3])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
770
{
771
  intercep_ratio[[i - 1]] <- 1 - cumu_sum[[i - 1]] / cumu_sum[[i]]
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
772
}
773
names(intercep_ratio) <- gsub("nb", "intercepRatio", names(cumu_sum)[-1])
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
774
# merge rasterstacks
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
775
metrics <- c(metrics, intercep_ratio)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
776
#
777
# ----------------------
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
778
779
780
781
# export as raster files
# for (i in 1:names(metrics))
# {
#  print(i)
782
#  writeRaster(metrics[[i]],file=paste("output/raster_",i,"_",resolution,"m.tif",sep="")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
783
# }
784
# -------
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
785
# display
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
786
terra::plot(metrics[[c("Imean_1stpulse", "Hsimpson")]],
787
788
             main = c("Mean intensity of 1st pulse",
                      "Simpson indice of point heights")
789
)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
790
791
```

792
## References