coregistration.Rmd 14.2 KB
Newer Older
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
1
---
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
2
title: "R workflow for field plot coregistration with ALS data"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
3
author: "Jean-Matthieu Monnet"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
4
date: "`r Sys.Date()`"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
5
6
7
output:
  pdf_document: default
  html_document: default
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
8
bibliography: "./bib/bibliography.bib"
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
9
10
11
12
13
14
15
16
17
18
19
---

```{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)
knitr::opts_chunk$set(fig.align = "center")
```


---
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
20
This workflow compares a canopy height model (CHM) derived from airborne laser scanning (ALS) data with tree positions inventoried in the field, and then proposes a translation in plot position for better matching. The method is described in @Monnet2014. Here it is exemplified with circular plots, but it can be applied to any shape of field plots. The workflow is based on functions from packages `lidaRtRee` and `lidR`. Example data were acquired in the forest of Lac des Rouges Truites (Jura, France).
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
21

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
22
Licence: CC BY / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/coregistration.Rmd)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38


## Material

### Field data

The study area is a part of the forest of Lac des Rouges Truites. 44 plots have been inventoried, 15 are available for testing. Plots have a 14.10 m radius. A data.frame `p` contains the positions of the center of plots. Attributes are:

* `placette`: plot id
- `Xtheo` and `Ytheo`: XY coordinates initially sampled when preparing the field inventory
- `XGPS` and `YGPS`: XY coordinates recorded in the field with a GNSS receiver during the field inventory
- `Prec`: GNSS precision in meter specified by the receiver
- `dist`: horizontal distance between the sample and recorded coordinates.

```{r plots, include = TRUE}
# load plot coordinates (data.frame with lines corresponding to the las objects)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
39
load(file="./data/coregistration/plotsCoregistration.rda")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
head(p, n=3L)
```

On each plot, five trees which were considered  suitable for coregistration (vertical trunk, high and peak-shaped crown, well separated from neighboring trees) have been positioned relatively to the plot center. From the XY coordinates recorded by the GNSS receiver and the relative coordinates (slope, distance, azimuth), the XY coordinates of those trees are computed. Data.frame `ap` contains the following attributes:

* `plac`: plot id
- `n`: tree id
- `dia`: tree diameter in cm
- `distR`: slope distance from the plot center to tree center, in m
- `azimutG`: azimuth from the plot center to the tree center, in grades
- `pente.`: slope from plot center to tree center, in grades
- `XYZGPS`: coordinates of the plot center
- `xyz`: coordinates of the tree center
- `d`: horizontal distance between plot center and tree in m


```{r trees, include = TRUE}
# load inventoried trees (data.frame with plot id info )
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
58
load(file="./data/coregistration/treesCoregistration.rda")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
59
60
61
62
63
64
65
head(ap, n=3L)
```

### ALS data

Airborne laser scanning data on the study area is part of a campaign acquired in 2016 with an airborne RIEGL LMS Q680i sensor. Acquisition was funded by the Région Franche-Comté.

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
66
ALS data over the plots is provided as a list of LAS objects in `rda` file. 
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
67
68
69

```{r las, include = TRUE}
# load point cloud over reference plots (list of las objects)
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
70
load(file="./data/coregistration/lasCoregistration.rda")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
```

Display point cloud of plot 1.

```{r lasplot, include=TRUE, eval=FALSE, webgl=TRUE, warning=FALSE}
# plot point cloud
plot(las[[1]])
```

The code to extract LAS objects from a directory containing the LAS files is (code corresponding to parameters in the next paragraph has to be run beforehand):

```{r extractlas, include = TRUE, eval=FALSE}
# create catalog of LAS files
cata <- lidR::catalog("/directory_with_classified_laz_files/")
# set coordinate system
sp::proj4string(cata) <- sp::CRS(SRS_string = "EPSG:2154")
# extract LAS data on plot extent
las <- lidR::clip_circle(cata, p$XGPS, p$YGPS, p.radius + b.size + 5)
# normalize heights if point cloud are not already normalized
las <- lapply(las, function(x) {lidR::normalize_height(x, lidR::tin())})
# save as rda file for later use
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
92
save(las, file="./data/coregistration/lasCoregistration.rda")
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
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
```

## Parameters

Several parameters have to be specified for optimal results.

```{r parameters, include = TRUE}
# vegetation height threshold for removal of high values
hmax <- 50
# field plot radius for computation of pseudo-CHM on the inventory area
p.radius <- 14.105
# raster resolution for image processing
r.res <- 0.5
# maximum distance around initial position to search for coregistration
b.size <- 5
# increment step to search for coregistration (should be a multiple of resolution)
s.step <- 0.5
```

## Coregistration of one plot

### Canopy height model

The first step is to compute the canopy height model from the ALS data, and remove artefacts by thresholding extreme values and applying a median filter.

```{r chm, include = TRUE, out.width='70%', fig.width = 9, fig.height = 4.5}
# Choose a plot as example
i <- 10
# compute canopy height model
chm <- lidaRtRee::points2DSM(las[[i]],res=r.res)
# apply threshold to canopy height model (CHM)
chm[chm>hmax] <- hmax
# fill NA values with 0
chm[is.na(chm)] <- 0
# apply median filter (3x3 window) to CHM
chmfilt <- lidaRtRee::demFiltering(chm,"Median",3, sigmap=0)$non.linear.image
# plot canopy height model
par(mfrow=c(1,2))
raster::plot(chm, main="Raw canopy height model")
raster::plot(chmfilt, main="Filtered canopy height model")
```

### Plot mask from tree inventory

The trees corresponding to the plot are extracted, and a plot mask is computed from the plot center and radius.

Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
139
```{r mask, include = TRUE, out.width = '50%', fig.dim=c(5.5, 4.5), warning=FALSE}
Monnet Jean-Matthieu's avatar
Monnet Jean-Matthieu committed
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
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
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
309
310
311
312
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
345
346
347
348
349
# plot centre
centre <- p[i,c("XGPS","YGPS")]
# extract plot trees
trees <- ap[ap$plac==p$placette[i],]
# create raster with plot extent
r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p.radius, resolution=r.res)
# keep only trees with diameter information
trees <- trees[!is.na(trees[,"dia"]),]
# create plot mask
r.mask <- lidaRtRee::rasterXYMask(rbind(c(centre$XGPS, centre$YGPS),
                                        c(centre$XGPS, centre$YGPS)),
                                  c(p.radius,p.radius), r, binary=T)
# replace 0 by NA
r.mask[r.mask==0] <- NA
# specify projection
raster::crs(r.mask) <- 2154
# display plot mask
raster::plot(r.mask, main="Plot mask and tree positions")
# add tree positions
points(trees[,c("x","y")], cex=trees[,"dia"]/40)
# add plot center
points(centre, pch=3)
```

### Compute correlation and identify optimal translation

First the function computes the correlation for different possible translations of the plot center inside the buffer specified by the user, and outputs an image of correlation values between the ALS CHM and the pseudo CHM. Second this image is analyzed to identify which translation yields the highest correlation. The function outputs a list which first element is the correlation image, and the second one the corresponding statistics.

```{r correlation, include = TRUE}
# compute correlation image 
coreg <- lidaRtRee::coregistration(chmfilt, trees[,c("x", "y", "dia")], mask=r.mask,
                                   buffer=b.size, step=s.step, dm=2, plot=F)
# correlation image statistics
round(coreg$local.max, 2)
```
The maximum of in the correlation image is located at (`r coreg$local.max$dx1`, `r coreg$local.max$dy1`), given by attributes `dx1` and `dy1`.

```{r plotcorrelation, include = TRUE, out.width='100%', fig.width = 12, fig.height = 4.5}
par(mfrow=c(1,3))
raster::plot(chmfilt, main="Initial tree positions")
# display initial tree positions
graphics::points(trees$x, trees$y, cex=trees$dia/40)
# display correlation image
raster::plot(coreg$correlation.raster, main="Correlation image")
# plot local maximum
graphics::points(coreg$local.max$dx1, coreg$local.max$dy1, pch=3, col="blue")
abline(h=0, lty=2); abline(v=0, lty=2)
#
raster::plot(chmfilt, main="Coregistered tree positions")
# display coregistered tree positions
graphics::points(trees$x+coreg$local.max$dx1, trees$y+coreg$local.max$dy1,
                 cex=trees$dia/40, col="red")
# display initial plot center
graphics::points(p[i, c("XGPS","YGPS")], pch=3)
# display coregistered plot center
graphics::points(p$XGPS[i]+coreg$local.max$dx1, p$YGPS[i]+coreg$local.max$dy1, pch=3,
                 col="red")
graphics::legend("topleft", c("Initial center", "Coregistered"), pch=3,
                 col=c("black", "red"))
# export as pdf
# dev.copy2pdf(file=paste("Coregistration_",p$placette[i],".pdf",sep=""))

```
## Batch processing

The following code processes several plots using multi-core computing. It is based on packages `foreach` and `doParallel` for parallel computing.

```{r batch.cores, include = TRUE, eval=TRUE}
# number of cores to use for parallel processing
doParallel::registerDoParallel(cores=4)
library(foreach)
```

### CHMs computation

First CHMs are calculated for each point cloud contained in the list.
```{r batch.chm, include = TRUE, eval=TRUE}
l.chm <- foreach::foreach (i=1:length(las), .errorhandling="remove") %dopar%
{
  # compute CHM
  chm <- lidaRtRee::points2DSM(las[[i]],res=r.res)
  # apply threshold to canopy height model (CHM)
  chm[chm>hmax] <- hmax
  # fill NA values with 0
  chm[is.na(chm)] <- 0
  # apply median filter (3x3 window) to CHM
  chmfilt <- lidaRtRee::demFiltering(chm,"Median",3, sigmap=0)[[1]]
  return(chmfilt)
}
```

### Trees lists extraction and plot masks computation

```{r batch.mask, include = TRUE, eval=TRUE}
l.field <- foreach::foreach (i=1:length(las), .errorhandling="remove") %dopar%
{
  # plot centre
  centre <- p[i,c("XGPS","YGPS")]
  # extract plot trees
  trees <- ap[ap$plac==p$placette[i],]
  # create raster with plot extent
  r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p.radius, resolution=r.res)
  # keep only trees with diameter information
  trees <- trees[!is.na(trees[,"dia"]),]
  # create plot mask
  r.mask <- lidaRtRee::rasterXYMask(rbind(c(centre$XGPS, centre$YGPS),
                                          c(centre$XGPS, centre$YGPS)),
                                    c(p.radius,p.radius), r, binary=T)
  # replace 0 by NA
  r.mask[r.mask==0] <- NA
  # specify projection
  raster::crs(r.mask) <- 2154
  #
  return(list(trees=trees[,c("x","y","dia")], r.mask=r.mask))
}
```

### Correlation images and corresponding statistics computation

Then the correlation image and corresponding statistics are computed for each plot.
```{r batch.correlation, include = TRUE, eval=TRUE}
l.coreg <- foreach::foreach (i=1:length(las), .errorhandling="remove") %dopar%
{
  lidaRtRee::coregistration(l.chm[[i]], l.field[[i]]$trees, mask=l.field[[i]]$r.mask,
                            buffer=b.size, step=s.step, dm=2, plot=F)
}
```
Finally results for all plots are combined in a single data.frame. 
```{r batch.translations, include = TRUE, eval=TRUE}
translations <- foreach::foreach (i=1:length(las), .combine=rbind,
                                  .errorhandling="remove") %dopar%
{
    # create data.frame with coregistration results and new plot coordinates
  data.frame(plotid=p$placette[i], X.cor=p$XGPS[i]+l.coreg[[i]]$local.max$dx1,
             Y.cor=p$YGPS[i]+l.coreg[[i]]$local.max$dy1, l.coreg[[i]]$local.max)
}
# 
round(head(translations, n=3L),2)
```

### Export graphics as pdf files

```{r batch.plots, include = TRUE, eval=FALSE}
for (i in 1:length(las))
{
    # CHM
    pdf(file=paste("Coregistration_",p$placette[i],".pdf",sep=""))
    raster::plot(l.chm[[i]], asp=1)
    # display initial tree positions
    graphics::points(l.field[[i]]$trees$x, l.field[[i]]$trees$y,
                     cex=l.field[[i]]$trees$dia/40)
    # display coregistered tree positions
    graphics::points(l.field[[i]]$trees$x+l.coreg[[i]]$local.max$dx1,
                     l.field[[i]]$trees$y+l.coreg[[i]]$local.max$dy1,
                     cex=l.field[[i]]$trees$dia/40, col="red")
    # display initial plot center
    graphics::points(p[i, c("XGPS","YGPS")], pch=3)
    # display coregistered plot center
    graphics::points(p$XGPS[i]+l.coreg[[i]]$local.max$dx1,
                     p$YGPS[i]+l.coreg[[i]]$local.max$dy1, pch=3, col="red")
    graphics::legend("topleft", c("Initial", "Coregistered"), pch=15,
                     col=c("black", "red"))
    # export as pdf
    dev.off()
}
```

### Add coregistered plot positions and update tree coordinates

Optimized plot center positions are added to the original data.
```{r batch.merge, include = TRUE, eval=TRUE}
p <- base::merge(p, translations[,c("plotid", "X.cor", "Y.cor")], by.x="placette",
                 by.y="plotid")
round(head(p, n=3L),2)
```
Tree positions are corrected to account for new center position.
```{r batch.correct, include = TRUE, eval=TRUE}
# add plot center coregistered coordinates to trees data.frame
ap <- base::merge(ap, p[,c("placette","X.cor","Y.cor")], by.x="plac", by.y="placette")
# compute new tree coordinates from coregistered plot center
dummy <- lidaRtRee::polar2Projected(ap$X.cor, ap$Y.cor, ap$ZGPS, ap$azimutG/200*pi,
                                    ap$distR, atan(ap$pente./100), 1.55/180*pi,
                                    -2.2/180*pi,0)
ap$X.cor <- dummy$x
ap$Y.cor <- dummy$y
# save new table
# write.table(round(p.cor,3),file="coregistered_plots.csv", row.names=F,sep=";")
```

### Analysis of GPS errors

```{r batch.distance, include = FALSE, eval=TRUE}
distance <- sqrt((p$X.cor- p$XGPS)^2 + (p$Y.cor- p$YGPS)^2)
```

Graphics on the difference between initial and corrected positions are displayed. The mean difference is `r round(mean(p$X.cor- p$XGPS), 2)`m in the X axis and `r round(mean(p$Y.cor- p$YGPS),2)`m in the Y axis. `p.values` of a t.test are respectively `r round(t.test(p$X.cor- p$XGPS)$p.value,2)` and `r round(t.test(p$Y.cor- p$YGPS)$p.value, 2)`. Mean absolute distance is `r round(mean(distance), 2)`m with a standard deviation of `r round(sd(distance),2)`m.


```{r batch.distanceplot, include = TRUE, out.width='80%', fig.width = 8, fig.height = 4}
par(mfrow=c(1,2))
# plot position difference with additionnal jitter to visualize same points
plot(jitter(p$X.cor- p$XGPS), jitter(p$Y.cor-p$YGPS), asp=1, col="black",
     main="Corrected-Initial position", xlab="X difference", ylab="Y difference")
abline(v=0, lty=2); abline(h=0, lty=2)
# Distance between initial and corrected
hist(sqrt((p$X.cor- p$XGPS)^2+ (p$Y.cor-p$YGPS)^2), main="Histogram of distances",
     xlab="Absolute distance corrected-initial", ylab="Number of plots")
```

## References