-
-
Notifications
You must be signed in to change notification settings - Fork 9
/
itd.qmd
235 lines (168 loc) · 14.1 KB
/
itd.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
```{r,echo=FALSE,message=FALSE,warning=FALSE}
r3dDefaults = rgl::r3dDefaults
m = structure(c(0.921, -0.146, 0.362, 0, 0.386, 0.482, -0.787, 0,
-0.06, 0.864, 0.5, 0, 0, 0, 0, 1), .Dim = c(4L, 4L))
rgl::setupKnitr()
r3dDefaults$FOV = 50
r3dDefaults$userMatrix = m
r3dDefaults$zoom = 0.75
rgl::setupKnitr()
knitr::opts_chunk$set(
comment = "#>",
collapse = TRUE,
fig.align = "center")
library(lidR)
options(lidR.progress = F)
```
# Indivitual tree dectection and segmentation {#sec-itd-its}
**Individual Tree Detection (ITD)** is the process of spatially locating trees and extracting height information. **Individual Tree Segmentation (ITS)** is the process of individually delineating detected trees. In `lidR`, detection and segmentation functions are decoupled to maximize flexibility. Tree tops are first detected using the `locate_trees()` function, followed by crown delineation using `segment_trees()`. In the following section, we will use the `MixedConifer.laz` data set, which is included internally within `lidR`, to demonstrate both ITD and ITS with reproducible examples. We will also generate a CHM (@sec-chm) to help visualize the results.
```{r plot-las-mixed-itd, rgl = TRUE, fig.width = 5, fig.height = 4}
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, select = "xyzr", filter = "-drop_z_below 0")
chm <- rasterize_canopy(las, 0.5, pitfree(subcircle = 0.2))
plot(las, bg = "white", size = 4)
```
## Individual Tree Detection (ITD) {#sec-itd}
Tree tops can be detected by applying a Local Maximum Filter (LMF) on the loaded data set. The LMF in `lidR` is point cloud-based, meaning that it finds the tree tops from the point cloud without using any raster. The processing, however is actually very similar. For a given point, the algorithm analyzes neighbourhood points, checking if the processed point is the highest. This algorithm can be applied with the `lmf()` function.
### Local Maximum Filter with fixed windows size {#sec-lmffw}
The LMF can be applied with a constant size windows. Here with a windows size of `ws = 5` meters meaning that for a given point the algorithm looks to the neigbourhood points within a 2.5 radius circle to figure out if the point is the local highest. While the algorithm does not need any CHM to work we chose to display the results on top of a CHM for better visualization.
```{r plot-chm-tree-itd, fig.height=6, fig.width=7}
ttops <- locate_trees(las, lmf(ws = 5))
plot(chm, col = height.colors(50))
plot(sf::st_geometry(ttops), add = TRUE, pch = 3)
```
Tree detection results can also be visualized in 3D!
```{r plot-las-tree-itd, rgl = TRUE, fig.width = 4, fig.height = 3}
x <- plot(las, bg = "white", size = 4)
add_treetops3d(x, ttops)
```
The number of detected trees is correlated to the `ws` argument. Small windows sizes usually gives more trees, while large windows size generally miss smaller trees that are *"hidden"* by big trees that contain the highest points in the neighbourhood. This can be seen in the figure below where too many trees are found with a small window size and only the dominant trees are found with a large windows size.
```{r plot-chm-tree-ws, fig.height=4, fig.width=10}
ttops_3m <- locate_trees(las, lmf(ws = 3))
ttops_11m <- locate_trees(las, lmf(ws = 11))
par(mfrow=c(1,2))
plot(chm, col = height.colors(50))
plot(sf::st_geometry(ttops_3m), add = TRUE, pch = 3)
plot(chm, col = height.colors(50))
plot(sf::st_geometry(ttops_11m), add = TRUE, pch = 3)
```
### Local Maximum Filter with variable windows size {#sec-lmfvw}
The examples above demonstrate the `lmf()` function with a constant window size. A large windows is suitable for large scattered trees while a small windows size if preferable for small and close trees. In reality trees of variable sizes may be present in a single scene leading to sub-optimal outputs. To overcome this issue, a windows size that adapts to the height of pixels or height of the points in our case becomes necessary.
For example, a point at 30 m (a big tree) will be tested with a large window, while a point at 10 m (a smaller tree) will be tested with a smaller window. The window size can therefore be defined as a function of height. Tall trees have comparatively larger crowns, needing larger window sizes to accurately detect their treetops. Variable window sizes are especially suitable for stands of more complex structures, or when tree detection is performed on larger areas, covering heterogeneous stands.
In `lidR` a user can design a function that computes a windows size as a function of point (or pixel) height. When designing a function to define the window size based on point heights we need to determine what the minimum and maximum window size should be related to the minimum and maximum tree heights. In general, the minimum window size should not be smaller than 3 meters. Below we show an example where the windows size is related to the point height with an affine relationship. When a point is at 0 the windows size is 3 meters. At 10 m it is 4 m and so on.
```{r plot-vws, fig.height=4, fig.width=3.4}
f <- function(x) {x * 0.1 + 3}
heights <- seq(0,30,5)
ws <- f(heights)
plot(heights, ws, type = "l", ylim = c(0,6))
```
When applied within `lmf()`, the function yields the following result:
```{r plot-chm-tree-vws, fig.height=4.2, fig.width=5.2}
ttops <- locate_trees(las, lmf(f))
plot(chm, col = height.colors(50))
plot(sf::st_geometry(ttops), add = TRUE, pch = 3)
```
There is no intrinsic limitation to the user-defined function. One does however need to pay attention to special cases. For example, if a point is below 0 m the function above may compute a negative windows size and the detection will fail. Similarly, if an outlier is encountered a massive window may be computed and lead to suspicious results. We therefore recommend using a more robust function with some built in thresholds.
In the next example any points below 2 m will equate to a window size of 3 m, while points above 20 meters equate to a window size of 5 m. Anything between 2 and 20 meter will have a non-linear relationship (for the need of the demo).
```{r plot-vws-2, fig.height=4, fig.width=3.4}
f <- function(x) {
y <- 2.6 * (-(exp(-0.08*(x-2)) - 1)) + 3
y[x < 2] <- 3
y[x > 20] <- 5
return(y)
}
heights <- seq(-5,30,0.5)
ws <- f(heights)
plot(heights, ws, type = "l", ylim = c(0,5))
```
### Local Maximum Filter on a CHM {#sec-lmfchm}
So far tree detection was performed on a point cloud. A CHM can however also be used to find trees instead of the point cloud. This is also a perfectly valid way to process ALS data. Performing the detection on a CHM is faster because there are less data to process, but it is also more complex because the output depends on how the CHM has been built. The spatial resolution (i.e. pixel size), the algorithm used (see section @sec-chm), and any additional post-processing steps will influence tree detection results.
In the examples below we run tree detection on CHMs generated with different algorithms (`p2r`, `pitfree`), different resolutions (0.5 and 1 m), and different post-processing smoothing steps with a median filter applied.
First the different CHMs are generated:
```{r}
# Point-to-raster 2 resolutions
chm_p2r_05 <- rasterize_canopy(las, 0.5, p2r(subcircle = 0.2), pkg = "terra")
chm_p2r_1 <- rasterize_canopy(las, 1, p2r(subcircle = 0.2), pkg = "terra")
# Pitfree with and without subcircle tweak
chm_pitfree_05_1 <- rasterize_canopy(las, 0.5, pitfree(), pkg = "terra")
chm_pitfree_05_2 <- rasterize_canopy(las, 0.5, pitfree(subcircle = 0.2), pkg = "terra")
# Post-processing median filter
kernel <- matrix(1,3,3)
chm_p2r_05_smoothed <- terra::focal(chm_p2r_05, w = kernel, fun = median, na.rm = TRUE)
chm_p2r_1_smoothed <- terra::focal(chm_p2r_1, w = kernel, fun = median, na.rm = TRUE)
```
Then the same tree detection routine with a constant windows size of 5 m is applied to each CHM:
```{r}
ttops_chm_p2r_05 <- locate_trees(chm_p2r_05, lmf(5))
ttops_chm_p2r_1 <- locate_trees(chm_p2r_1, lmf(5))
ttops_chm_pitfree_05_1 <- locate_trees(chm_pitfree_05_1, lmf(5))
ttops_chm_pitfree_05_2 <- locate_trees(chm_pitfree_05_2, lmf(5))
ttops_chm_p2r_05_smoothed <- locate_trees(chm_p2r_05_smoothed, lmf(5))
ttops_chm_p2r_1_smoothed <- locate_trees(chm_p2r_1_smoothed, lmf(5))
```
Finally the detection results are visualized to see the various output found as a function of the CHM building choices:
```{r plot-multichm-lmf, fig.height=12.5, fig.width=9.8}
par(mfrow=c(3,2))
col <- height.colors(50)
plot(chm_p2r_05, main = "CHM P2R 0.5", col = col); plot(sf::st_geometry(ttops_chm_p2r_05), add = T, pch =3)
plot(chm_p2r_1, main = "CHM P2R 1", col = col); plot(sf::st_geometry(ttops_chm_p2r_1), add = T, pch = 3)
plot(chm_p2r_05_smoothed, main = "CHM P2R 0.5 smoothed", col = col); plot(sf::st_geometry(ttops_chm_p2r_05_smoothed), add = T, pch =3)
plot(chm_p2r_1_smoothed, main = "CHM P2R 1 smoothed", col = col); plot(sf::st_geometry(ttops_chm_p2r_1_smoothed), add = T, pch =3)
plot(chm_pitfree_05_1, main = "CHM PITFREE 1", col = col); plot(sf::st_geometry(ttops_chm_pitfree_05_1), add = T, pch =3)
plot(chm_pitfree_05_2, main = "CHM PITFREE 2", col = col); plot(sf::st_geometry(ttops_chm_pitfree_05_2), add = T, pch =3)
```
## Individual Tree Segmentation (ITS) {#sec-its}
While individual tree detection provides very useful information about tree density and size, one may want to go further to segment and extract each tree individually. Several algorithms are available in `lidR` and can be divided in two families.
1. Point cloud-based that perform without a CHM.
2. Raster-based that perform with a CHM.
Each family can be divided into two sub families
1. Algorithms that work in two steps - Individual tree detection followed by segmentation.
2. Algorithms that are all-in-one.
In this section we won't go through all of these possibilities, but its important to recognize that all algorithms are not well suited to every context. In the examples below we used the `dalponte2016()` algorithm, because we found it to perform best using the example data.
### Segmentation of the point-cloud {#sec-its-cloud}
Even when the algorithm is raster-based (which is the case of `dalponte2016()`), `lidR` segments the point cloud and assigns an ID to each point by inserting a new attribute named `treeID` in the `LAS` object. This means that every point is associated with a particular tree. This is because `lidR` is point cloud oriented and we want to provide an immediate way to be able to access a segmented point cloud not the segmented raster.
```{r plot-las-its, rgl = TRUE, fig.width = 5, fig.height = 4, warning = FALSE}
algo <- dalponte2016(chm_p2r_05_smoothed, ttops_chm_p2r_05_smoothed)
las <- segment_trees(las, algo) # segment point cloud
plot(las, bg = "white", size = 4, color = "treeID") # visualize trees
```
We assign an ID to each point because it enables interesting analyses at the point-cloud level. This could involve extracting every tree to derive measurements. In the following we extracted and displayed the tree number 110.
```{r plot-las-tree110, rgl = TRUE, fig.width = 3, fig.height = 4}
tree110 <- filter_poi(las, treeID == 110)
plot(tree110, size = 8, bg = "white")
```
`lidR` provides functions that can be used following `segment_trees()`, like delineation of crown shapes using `crown_metrics()`. This function computes the hulls (either convex or concave) of each tree as well as user-defined metrics (see also sections @sec-metrics, @sec-cba, @sec-aba, @sec-tba, @sec-vba and @sec-pba):
```{r plot-tree-hulls-its, fig.height=4, fig.width=4}
crowns <- crown_metrics(las, func = .stdtreemetrics, geom = "convex")
plot(crowns["convhull_area"], main = "Crown area (convex hull)")
```
### Segmentation of the CHM {#sec-its-chm}
While point cloud segmentation is standard in `lidR`, users may only have access to a CHM. There are many reasons for only using a CHM, and this is why raster-based methods can be run standalone outside `segment_trees()`. Whether using the point cloud or a raster, segmentation results will be exactly the same. The difference will be the data format of the segmentation result. In `lidR`, a `LAS` object will gain a `treeID` attribute, while for rasters, delineated crowns are returned in a raster format. To work outside `segment_trees()` it suffices to call the function standalone like this:
```{r plot-raster-its, fig.height=4.2, fig.width=5.2}
algo <- dalponte2016(chm_p2r_05_smoothed, ttops_chm_p2r_05_smoothed)
crowns <- algo()
plot(crowns, col = pastel.colors(200))
```
The output is a raster, for which `lidR` does not provide any processing tools. At this stage it is up to the user to find the required tools to perform more analysis. The `terra` package is a good place to start!
### Comparaison of tree segmentations {#sec-its-comparaison}
At the point cloud level it's pretty easy to compare tree segmentations by choosing a different attribute names for each method. For example we can compare `dalponte2016` (which is a raster based methods in two steps) and `li2012` (which is an *"all-in-one"* point cloud based method) side by side.
```{r, echo = FALSE}
r3dDefaults <- rgl::r3dDefaults
r3dDefaults$zoom <- 0.7
```
```{r plot-compare-its, rgl = TRUE, warning = FALSE, fig.height=6, fig.width=8}
algo1 <- dalponte2016(chm_p2r_05_smoothed, ttops_chm_p2r_05_smoothed)
algo2 <- li2012()
las <- segment_trees(las, algo1, attribute = "IDdalponte")
las <- segment_trees(las, algo2, attribute = "IDli")
x <- plot(las, bg = "white", size = 4, color = "IDdalponte", colorPalette = pastel.colors(200))
plot(las, add = x + c(100,0), bg = "white", size = 4, color = "IDli", colorPalette = pastel.colors(200))
```
By comparing their crowns with `crown_metrics()` we can more easily see that the `li2012` algorithm doesn't perform well in this example with default parameters - maybe some parameter tuning can lead to better results!
```{r plot-compare-hulls, fig.height=3, fig.width=5.8}
crowns_dalponte <- crown_metrics(las, func = NULL, attribute = "IDdalponte", geom = "concave")
crowns_li <- crown_metrics(las, func = NULL, attribute = "IDli", geom = "concave")
par(mfrow=c(1,2),mar=rep(0,4))
plot(sf::st_geometry(crowns_dalponte), reset = FALSE)
plot(sf::st_geometry(crowns_li), reset = FALSE)
```