Commit b47c3de7 authored by Monnet Jean-Matthieu's avatar Monnet Jean-Matthieu
Browse files

Functions and variables names modifications to comply with lidaRtRee-dev

No related merge requests found
Showing with 1027 additions and 1028 deletions
+1027 -1028
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
...@@ -3,8 +3,8 @@ title: "R workflow for field plot coregistration with ALS data" ...@@ -3,8 +3,8 @@ title: "R workflow for field plot coregistration with ALS data"
author: "Jean-Matthieu Monnet" author: "Jean-Matthieu Monnet"
date: "`r Sys.Date()`" date: "`r Sys.Date()`"
output: output:
pdf_document: default
html_document: default html_document: default
pdf_document: default
papersize: a4 papersize: a4
bibliography: "../bib/bibliography.bib" bibliography: "../bib/bibliography.bib"
--- ---
...@@ -90,7 +90,7 @@ cata <- lidR::readALSLAScatalog("/directory_with_classified_laz_files/") ...@@ -90,7 +90,7 @@ cata <- lidR::readALSLAScatalog("/directory_with_classified_laz_files/")
# set coordinate system # set coordinate system
lidR::projection(cata) <- 2154 lidR::projection(cata) <- 2154
# extract LAS data on plot extent # extract LAS data on plot extent
las <- lidR::clip_circle(cata, p$XGPS, p$YGPS, p.radius + b.size + 5) las <- lidR::clip_circle(cata, p$XGPS, p$YGPS, p_radius + b_size + 5)
# normalize heights if point cloud are not already normalized # normalize heights if point cloud are not already normalized
las <- lapply(las, function(x) { las <- lapply(las, function(x) {
lidR::normalize_height(x, lidR::tin()) lidR::normalize_height(x, lidR::tin())
...@@ -107,13 +107,13 @@ Several parameters have to be specified for optimal results. ...@@ -107,13 +107,13 @@ Several parameters have to be specified for optimal results.
# vegetation height threshold for removal of high values # vegetation height threshold for removal of high values
hmax <- 50 hmax <- 50
# field plot radius for computation of pseudo-CHM on the inventory area # field plot radius for computation of pseudo-CHM on the inventory area
p.radius <- 14.105 p_radius <- 14.105
# raster resolution for image processing # raster resolution for image processing
r.res <- 0.5 r_res <- 0.5
# maximum distance around initial position to search for coregistration # maximum distance around initial position to search for coregistration
b.size <- 5 b_size <- 5
# increment step to search for coregistration (should be a multiple of resolution) # increment step to search for coregistration (should be a multiple of resolution)
s.step <- 0.5 s_step <- 0.5
``` ```
## Coregistration of one plot ## Coregistration of one plot
...@@ -126,13 +126,13 @@ The first step is to compute the canopy height model from the ALS data, and remo ...@@ -126,13 +126,13 @@ The first step is to compute the canopy height model from the ALS data, and remo
# Choose a plot as example # Choose a plot as example
i <- 10 i <- 10
# compute canopy height model # compute canopy height model
chm <- lidaRtRee::points2DSM(las[[i]], res = r.res) chm <- lidaRtRee::points2DSM(las[[i]], res = r_res)
# apply threshold to canopy height model (CHM) # apply threshold to canopy height model (CHM)
chm[chm > hmax] <- hmax chm[chm > hmax] <- hmax
# fill NA values with 0 # fill NA values with 0
chm[is.na(chm)] <- 0 chm[is.na(chm)] <- 0
# apply median filter (3x3 window) to CHM # apply median filter (3x3 window) to CHM
chmfilt <- lidaRtRee::demFiltering(chm, "Median", 3, sigmap = 0)$non.linear.image chmfilt <- lidaRtRee::dem_filtering(chm, "Median", 3, sigmap = 0)$non_linear_image
# plot canopy height model # plot canopy height model
par(mfrow = c(1, 2)) par(mfrow = c(1, 2))
raster::plot(chm, main = "Raw canopy height model") raster::plot(chm, main = "Raw canopy height model")
...@@ -149,23 +149,23 @@ centre <- p[i, c("XGPS", "YGPS")] ...@@ -149,23 +149,23 @@ centre <- p[i, c("XGPS", "YGPS")]
# extract plot trees # extract plot trees
trees <- ap[ap$plac == p$placette[i], ] trees <- ap[ap$plac == p$placette[i], ]
# create raster with plot extent # create raster with plot extent
r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p.radius, resolution = r.res) r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p_radius, resolution = r_res)
# keep only trees with diameter information # keep only trees with diameter information
trees <- trees[!is.na(trees[, "dia"]), ] trees <- trees[!is.na(trees[, "dia"]), ]
# create plot mask # create plot mask
r.mask <- lidaRtRee::rasterXYMask(rbind( r_mask <- lidaRtRee::raster_xy_mask(rbind(
c(centre$XGPS, centre$YGPS), c(centre$XGPS, centre$YGPS),
c(centre$XGPS, centre$YGPS) c(centre$XGPS, centre$YGPS)
), ),
c(p.radius, p.radius), r, c(p_radius, p_radius), r,
binary = T binary = T
) )
# replace 0 by NA # replace 0 by NA
r.mask[r.mask == 0] <- NA r_mask[r_mask == 0] <- NA
# specify projection # specify projection
raster::crs(r.mask) <- 2154 raster::crs(r_mask) <- 2154
# display plot mask # display plot mask
raster::plot(r.mask, main = "Plot mask and tree positions", legend = FALSE) raster::plot(r_mask, main = "Plot mask and tree positions", legend = FALSE)
# add tree positions # add tree positions
points(trees[, c("x", "y")], cex = trees[, "dia"] / 40) points(trees[, c("x", "y")], cex = trees[, "dia"] / 40)
# add plot center # add plot center
...@@ -180,13 +180,13 @@ First the function computes the correlation for different possible translations ...@@ -180,13 +180,13 @@ First the function computes the correlation for different possible translations
```{r correlation, include = TRUE} ```{r correlation, include = TRUE}
# compute correlation image # compute correlation image
coreg <- lidaRtRee::coregistration(chmfilt, trees[, c("x", "y", "dia")], coreg <- lidaRtRee::coregistration(chmfilt, trees[, c("x", "y", "dia")],
mask = r.mask, mask = r_mask,
buffer = b.size, step = s.step, dm = 2, plot = F buffer = b_size, step = s_step, dm = 2, plot = FALSE
) )
# correlation image statistics # correlation image statistics
round(coreg$local.max, 2) round(coreg$local_max, 2)
``` ```
The maximum of the correlation image is located at (`r coreg$local.max$dx1`, `r coreg$local.max$dy1`), given by attributes `dx1` and `dy1`. The maximum of 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} ```{r plotcorrelation, include = TRUE, out.width='100%', fig.width = 12, fig.height = 4.5}
par(mfrow = c(1, 3)) par(mfrow = c(1, 3))
...@@ -194,25 +194,25 @@ raster::plot(chmfilt, main = "Initial tree positions and CHM") ...@@ -194,25 +194,25 @@ raster::plot(chmfilt, main = "Initial tree positions and CHM")
# display initial tree positions # display initial tree positions
graphics::points(trees$x, trees$y, cex = trees$dia / 40) graphics::points(trees$x, trees$y, cex = trees$dia / 40)
# display correlation image # display correlation image
raster::plot(coreg$correlation.raster, raster::plot(coreg$correlation_raster,
main = "Correlation image", main = "Correlation image",
col = cm.colors(16) col = cm.colors(16)
) )
# plot local maximum # plot local maximum
graphics::points(coreg$local.max$dx1, coreg$local.max$dy1, pch = 4) graphics::points(coreg$local_max$dx1, coreg$local_max$dy1, pch = 4)
abline(h = 0, lty = 2) abline(h = 0, lty = 2)
abline(v = 0, lty = 2) abline(v = 0, lty = 2)
legend("topleft", "Maximum", pch = 4) legend("topleft", "Maximum", pch = 4)
# #
raster::plot(chmfilt, main = "Coregistered tree positions and CHM") raster::plot(chmfilt, main = "Coregistered tree positions and CHM")
# display coregistered tree positions # display coregistered tree positions
graphics::points(trees$x + coreg$local.max$dx1, trees$y + coreg$local.max$dy1, graphics::points(trees$x + coreg$local_max$dx1, trees$y + coreg$local_max$dy1,
cex = trees$dia / 40, col = "red" cex = trees$dia / 40, col = "red"
) )
# display initial plot center # display initial plot center
graphics::points(p[i, c("XGPS", "YGPS")], pch = 3) graphics::points(p[i, c("XGPS", "YGPS")], pch = 3)
# display coregistered plot center # display coregistered plot center
graphics::points(p$XGPS[i] + coreg$local.max$dx1, p$YGPS[i] + coreg$local.max$dy1, graphics::points(p$XGPS[i] + coreg$local_max$dx1, p$YGPS[i] + coreg$local_max$dy1,
pch = 3, pch = 3,
col = "red" col = "red"
) )
...@@ -237,15 +237,15 @@ library(foreach) ...@@ -237,15 +237,15 @@ library(foreach)
First CHMs are calculated for each point cloud contained in the list. First CHMs are calculated for each point cloud contained in the list.
```{r batch.chm, include = TRUE, eval=TRUE} ```{r batch.chm, include = TRUE, eval=TRUE}
l.chm <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% { l_chm <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% {
# compute CHM # compute CHM
chm <- lidaRtRee::points2DSM(las[[i]], res = r.res) chm <- lidaRtRee::points2DSM(las[[i]], res = r_res)
# apply threshold to canopy height model (CHM) # apply threshold to canopy height model (CHM)
chm[chm > hmax] <- hmax chm[chm > hmax] <- hmax
# fill NA values with 0 # fill NA values with 0
chm[is.na(chm)] <- 0 chm[is.na(chm)] <- 0
# apply median filter (3x3 window) to CHM # apply median filter (3x3 window) to CHM
chmfilt <- lidaRtRee::demFiltering(chm, "Median", 3, sigmap = 0)[[1]] chmfilt <- lidaRtRee::dem_filtering(chm, "Median", 3, sigmap = 0)[[1]]
return(chmfilt) return(chmfilt)
} }
``` ```
...@@ -253,28 +253,28 @@ l.chm <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% ...@@ -253,28 +253,28 @@ l.chm <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar%
### Trees lists extraction and plot masks computation ### Trees lists extraction and plot masks computation
```{r batch.mask, include = TRUE, eval=TRUE} ```{r batch.mask, include = TRUE, eval=TRUE}
l.field <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% { l_field <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% {
# plot centre # plot centre
centre <- p[i, c("XGPS", "YGPS")] centre <- p[i, c("XGPS", "YGPS")]
# extract plot trees # extract plot trees
trees <- ap[ap$plac == p$placette[i], ] trees <- ap[ap$plac == p$placette[i], ]
# create raster with plot extent # create raster with plot extent
r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p.radius, resolution = r.res) r <- lidaRtRee::circle2Raster(centre$XGPS, centre$YGPS, p_radius, resolution = r_res)
# keep only trees with diameter information # keep only trees with diameter information
trees <- trees[!is.na(trees[, "dia"]), ] trees <- trees[!is.na(trees[, "dia"]), ]
# create plot mask # create plot mask
r.mask <- lidaRtRee::rasterXYMask(rbind( r_mask <- lidaRtRee::raster_xy_mask(rbind(
c(centre$XGPS, centre$YGPS), c(centre$XGPS, centre$YGPS),
c(centre$XGPS, centre$YGPS) c(centre$XGPS, centre$YGPS)
), ),
c(p.radius, p.radius), r, c(p_radius, p_radius), r,
binary = T binary = T
) )
# replace 0 by NA # replace 0 by NA
r.mask[r.mask == 0] <- NA r_mask[r_mask == 0] <- NA
# specify projection # specify projection
raster::crs(r.mask) <- 2154 raster::crs(r_mask) <- 2154
return(list(trees = trees[, c("x", "y", "dia")], r.mask = r.mask)) return(list(trees = trees[, c("x", "y", "dia")], r_mask = r_mask))
} }
``` ```
...@@ -282,10 +282,10 @@ l.field <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar ...@@ -282,10 +282,10 @@ l.field <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar
Then the correlation image and corresponding statistics are computed for each plot. Then the correlation image and corresponding statistics are computed for each plot.
```{r batch.correlation, include = TRUE, eval=TRUE} ```{r batch.correlation, include = TRUE, eval=TRUE}
l.coreg <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% { l_coreg <- foreach::foreach(i = 1:length(las), .errorhandling = "remove") %dopar% {
lidaRtRee::coregistration(l.chm[[i]], l.field[[i]]$trees, lidaRtRee::coregistration(l_chm[[i]], l_field[[i]]$trees,
mask = l.field[[i]]$r.mask, mask = l_field[[i]]$r_mask,
buffer = b.size, step = s.step, dm = 2, plot = F buffer = b_size, step = s_step, dm = 2, plot = FALSE
) )
} }
``` ```
...@@ -297,8 +297,8 @@ translations <- foreach::foreach( ...@@ -297,8 +297,8 @@ translations <- foreach::foreach(
) %dopar% { ) %dopar% {
# create data.frame with coregistration results and new plot coordinates # create data.frame with coregistration results and new plot coordinates
data.frame( data.frame(
plotid = p$placette[i], X.cor = p$XGPS[i] + l.coreg[[i]]$local.max$dx1, 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 Y_cor = p$YGPS[i] + l_coreg[[i]]$local_max$dy1, l_coreg[[i]]$local_max
) )
} }
# #
...@@ -312,21 +312,21 @@ for (i in 1:length(las)) ...@@ -312,21 +312,21 @@ for (i in 1:length(las))
{ {
# CHM # CHM
pdf(file = paste("Coregistration_", p$placette[i], ".pdf", sep = "")) pdf(file = paste("Coregistration_", p$placette[i], ".pdf", sep = ""))
raster::plot(l.chm[[i]], asp = 1) raster::plot(l_chm[[i]], asp = 1)
# display initial tree positions # display initial tree positions
graphics::points(l.field[[i]]$trees$x, l.field[[i]]$trees$y, graphics::points(l_field[[i]]$trees$x, l_field[[i]]$trees$y,
cex = l.field[[i]]$trees$dia / 40 cex = l_field[[i]]$trees$dia / 40
) )
# display coregistered tree positions # display coregistered tree positions
graphics::points(l.field[[i]]$trees$x + l.coreg[[i]]$local.max$dx1, graphics::points(l_field[[i]]$trees$x + l_coreg[[i]]$local_max$dx1,
l.field[[i]]$trees$y + l.coreg[[i]]$local.max$dy1, l_field[[i]]$trees$y + l_coreg[[i]]$local_max$dy1,
cex = l.field[[i]]$trees$dia / 40, col = "red" cex = l_field[[i]]$trees$dia / 40, col = "red"
) )
# display initial plot center # display initial plot center
graphics::points(p[i, c("XGPS", "YGPS")], pch = 3) graphics::points(p[i, c("XGPS", "YGPS")], pch = 3)
# display coregistered plot center # display coregistered plot center
graphics::points(p$XGPS[i] + l.coreg[[i]]$local.max$dx1, graphics::points(p$XGPS[i] + l_coreg[[i]]$local_max$dx1,
p$YGPS[i] + l.coreg[[i]]$local.max$dy1, p$YGPS[i] + l_coreg[[i]]$local_max$dy1,
pch = 3, col = "red" pch = 3, col = "red"
) )
graphics::legend("topleft", c("Initial", "Coregistered"), graphics::legend("topleft", c("Initial", "Coregistered"),
...@@ -342,7 +342,7 @@ for (i in 1:length(las)) ...@@ -342,7 +342,7 @@ for (i in 1:length(las))
Optimized plot center positions are added to the original data. Optimized plot center positions are added to the original data.
```{r batch.merge, include = TRUE, eval=TRUE} ```{r batch.merge, include = TRUE, eval=TRUE}
p <- base::merge(p, translations[, c("plotid", "X.cor", "Y.cor")], p <- base::merge(p, translations[, c("plotid", "X_cor", "Y_cor")],
by.x = "placette", by.x = "placette",
by.y = "plotid" by.y = "plotid"
) )
...@@ -351,15 +351,15 @@ round(head(p, n = 3L), 2) ...@@ -351,15 +351,15 @@ round(head(p, n = 3L), 2)
Tree positions are corrected to account for new center position. Tree positions are corrected to account for new center position.
```{r batch.correct, include = TRUE, eval=TRUE} ```{r batch.correct, include = TRUE, eval=TRUE}
# add plot center coregistered coordinates to trees data.frame # 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") 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 # compute new tree coordinates from coregistered plot center
dummy <- lidaRtRee::polar2Projected( dummy <- lidaRtRee::polar2Projected(
ap$X.cor, ap$Y.cor, ap$ZGPS, ap$azimutG / 200 * pi, ap$X_cor, ap$Y_cor, ap$ZGPS, ap$azimutG / 200 * pi,
ap$distR, atan(ap$pente. / 100), 1.55 / 180 * pi, ap$distR, atan(ap$pente. / 100), 1.55 / 180 * pi,
-2.2 / 180 * pi, 0 -2.2 / 180 * pi, 0
) )
ap$X.cor <- dummy$x ap$X_cor <- dummy$x
ap$Y.cor <- dummy$y ap$Y_cor <- dummy$y
# save new table # save new table
# write.table(round(p.cor,3),file="coregistered_plots.csv", row.names=F,sep=";") # write.table(round(p.cor,3),file="coregistered_plots.csv", row.names=F,sep=";")
``` ```
...@@ -367,26 +367,26 @@ ap$Y.cor <- dummy$y ...@@ -367,26 +367,26 @@ ap$Y.cor <- dummy$y
### Analysis of GPS errors ### Analysis of GPS errors
```{r batch.distance, include = FALSE, eval=TRUE} ```{r batch.distance, include = FALSE, eval=TRUE}
distance <- sqrt((p$X.cor - p$XGPS)^2 + (p$Y.cor - p$YGPS)^2) 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. 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} ```{r batch.distanceplot, include = TRUE, out.width='80%', fig.width = 8, fig.height = 4}
par(mfrow = c(1, 2)) par(mfrow = c(1, 2))
# plot position difference with additionnal jitter to visualize same points # plot position difference with additionnal jitter to visualize same points
plot(jitter(p$X.cor - p$XGPS), jitter(p$Y.cor - p$YGPS), plot(jitter(p$X_cor - p$XGPS), jitter(p$Y_cor - p$YGPS),
asp = 1, col = "black", asp = 1, col = "black",
main = "Corrected-Initial position", xlab = "X difference", ylab = "Y difference" main = "Corrected-Initial position", xlab = "X difference", ylab = "Y difference"
) )
abline(v = 0, lty = 2) abline(v = 0, lty = 2)
abline(h = 0, lty = 2) abline(h = 0, lty = 2)
# Distance between initial and corrected # Distance between initial and corrected
hist(sqrt((p$X.cor - p$XGPS)^2 + (p$Y.cor - p$YGPS)^2), hist(sqrt((p$X_cor - p$XGPS)^2 + (p$Y_cor - p$YGPS)^2),
main = "Histogram of distances", main = "Histogram of distances",
xlab = "Absolute distance corrected-initial", ylab = "Number of plots" xlab = "Absolute distance corrected-initial", ylab = "Number of plots"
) )
``` ```
## References ## References
\ No newline at end of file
This diff is collapsed.
...@@ -85,9 +85,9 @@ The procedure of gap detection is exemplified in @Glad20. ...@@ -85,9 +85,9 @@ The procedure of gap detection is exemplified in @Glad20.
```{r gaps, include = TRUE, fig.width = 6.5, fig.height = 6, warning=FALSE} ```{r gaps, include = TRUE, fig.width = 6.5, fig.height = 6, warning=FALSE}
# Perform gap detection with gap width criterion # Perform gap detection with gap width criterion
gaps1 <- lidaRtRee::gapDetection(chm, ratio = 2, gap.max.height = 1, min.gap.surface = 100) gaps1 <- lidaRtRee::gap_detection(chm, ratio = 2, gap_max_height = 1, min_gap_surface = 100)
# Perform gap detection without gap width nor min gap surface criteria # Perform gap detection without gap width nor min gap surface criteria
gaps2 <- lidaRtRee::gapDetection(chm, ratio = NULL, gap.max.height = 1, min.gap.surface = 0) gaps2 <- lidaRtRee::gap_detection(chm, ratio = NULL, gap_max_height = 1, min_gap_surface = 0)
``` ```
The following figure displays the obtained gaps, colored by size. The following figure displays the obtained gaps, colored by size.
...@@ -95,14 +95,14 @@ The following figure displays the obtained gaps, colored by size. ...@@ -95,14 +95,14 @@ The following figure displays the obtained gaps, colored by size.
# display CHM and areas where vegetation height is lower than 1m # display CHM and areas where vegetation height is lower than 1m
par(mfrow = c(1, 2)) par(mfrow = c(1, 2))
# max value of surface (log) # max value of surface (log)
z.lim <- log(max(raster::values(gaps1$gap.surface))) z_lim <- log(max(raster::values(gaps1$gap_surface)))
# plot gap surface # plot gap surface
raster::plot(log(gaps1$gap.surface), raster::plot(log(gaps1$gap_surface),
main = "log(Gap surface) (1)", zlim = c(0, z.lim), main = "log(Gap surface) (1)", zlim = c(0, z_lim),
col = rainbow(255, end = 5 / 6) col = rainbow(255, end = 5 / 6)
) )
raster::plot(log(gaps2$gap.surface), raster::plot(log(gaps2$gap_surface),
main = "log(Gap surface) (2)", zlim = c(0, z.lim), main = "log(Gap surface) (2)", zlim = c(0, z_lim),
col = rainbow(255, end = 5 / 6) col = rainbow(255, end = 5 / 6)
) )
``` ```
...@@ -111,32 +111,32 @@ Gaps can be vectorized for display over the original CHM. ...@@ -111,32 +111,32 @@ Gaps can be vectorized for display over the original CHM.
```{r gapsvectorize, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5)} ```{r gapsvectorize, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5)}
# convert to vector # convert to vector
gaps1.v <- raster::rasterToPolygons(gaps1$gap.id, dissolve = TRUE) gaps1_v <- raster::rasterToPolygons(gaps1$gap_id, dissolve = TRUE)
# display gaps over original CHM # display gaps over original CHM
raster::plot(chm, main = "Gaps over CHM (1)") raster::plot(chm, main = "Gaps over CHM (1)")
sp::plot(gaps1.v, add = T) sp::plot(gaps1_v, add = T)
``` ```
Statistics on gaps (number, size) can be derived from the raster objects. Example for gaps (2) Statistics on gaps (number, size) can be derived from the raster objects. Example for gaps (2)
```{r gapsstats, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5) } ```{r gapsstats, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5) }
# compute gaps surface # compute gaps surface
gaps.df <- data.frame(t(table(raster::values(gaps2$gap.id))))[, -1] gaps_df <- data.frame(t(table(raster::values(gaps2$gap_id))))[, -1]
# change column names # change column names
names(gaps.df) <- c("id", "surface") names(gaps_df) <- c("id", "surface")
# convert surface in pixels to square meters # convert surface in pixels to square meters
gaps.df$surface <- gaps.df$surface * raster::res(chm)[1]^2 gaps_df$surface <- gaps_df$surface * raster::res(chm)[1]^2
# data.frame # data.frame
head(gaps.df, n = 3) head(gaps_df, n = 3)
# cumulative distribution of gap surface by surface # cumulative distribution of gap surface by surface
Hmisc::Ecdf(log(gaps.df$surface), weights = gaps.df$surface, xlab = "log(surface)") Hmisc::Ecdf(log(gaps_df$surface), weights = gaps_df$surface, xlab = "log(surface)")
``` ```
Histogram can be used to compare the distribution of surface in different classes of gap surface. For better visualization of the differences between the set of criteria, gaps bigger than 500 m^2^ are not taken into account. Histogram can be used to compare the distribution of surface in different classes of gap surface. For better visualization of the differences between the set of criteria, gaps bigger than 500 m^2^ are not taken into account.
```{r gapsstats2, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5) } ```{r gapsstats2, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5) }
# extract vector of total gap surface of pixel # extract vector of total gap surface of pixel
surface1 <- raster::values(gaps1$gap.surface) surface1 <- raster::values(gaps1$gap_surface)
surface2 <- raster::values(gaps2$gap.surface) surface2 <- raster::values(gaps2$gap_surface)
# remove large gaps # remove large gaps
surface1[surface1 > 500] <- NA surface1[surface1 > 500] <- NA
surface2[surface2 > 500] <- NA surface2[surface2 > 500] <- NA
...@@ -152,38 +152,38 @@ In the previous part, a pixel that fulfills the maximum height criterion but doe ...@@ -152,38 +152,38 @@ In the previous part, a pixel that fulfills the maximum height criterion but doe
```{r gapsReconstruct, include = TRUE, out.width='100%', fig.width = 13, fig.height = 4.5, warning = FALSE} ```{r gapsReconstruct, include = TRUE, out.width='100%', fig.width = 13, fig.height = 4.5, warning = FALSE}
# load canopy height model # load canopy height model
chm <- lidaRtRee::chmchablais3 chm <- lidaRtRee::chm_chablais3
chm[is.na(chm)] <- 0 chm[is.na(chm)] <- 0
# Perform gap detection with distance ratio criterion # Perform gap detection with distance ratio criterion
gaps1 <- lidaRtRee::gapDetection(chm, gaps1 <- lidaRtRee::gap_detection(chm,
ratio = 2, gap.max.height = 1, min.gap.surface = 50, ratio = 2, gap_max_height = 1, min_gap_surface = 50,
nlFilter = "None" nl_filter = "None"
) )
# Perform gap detection with distance ratio criterion and gap reconstruction # Perform gap detection with distance ratio criterion and gap reconstruction
gaps1r <- lidaRtRee::gapDetection(chm, gaps1r <- lidaRtRee::gap_detection(chm,
ratio = 2, gap.max.height = 1, min.gap.surface = 50, ratio = 2, gap_max_height = 1, min_gap_surface = 50,
gapReconstruct = TRUE, nlFilter = "None" gap_reconstruct = TRUE, nl_filter = "None"
) )
# display detected gaps # display detected gaps
par(mfrow = c(1, 3)) par(mfrow = c(1, 3))
# plot gaps # plot gaps
raster::plot(chm, main = "Canopy heightmodel") raster::plot(chm, main = "Canopy heightmodel")
raster::plot(gaps1$gap.id > 0, main = "Gaps", legend = FALSE) raster::plot(gaps1$gap_id > 0, main = "Gaps", legend = FALSE)
raster::plot(gaps1r$gap.id > 0, main = "Gaps extended by reconstruction", legend = FALSE) raster::plot(gaps1r$gap_id > 0, main = "Gaps extended by reconstruction", legend = FALSE)
``` ```
Edge detection is performed by extraction of the difference between a binary image of gaps and the result of a morphological erosion or dilation applied to the same image. Edge detection is performed by extraction of the difference between a binary image of gaps and the result of a morphological erosion or dilation applied to the same image.
```{r edgeDetection, include = TRUE, out.width='70%', fig.width = 9, fig.height = 4.5} ```{r edgeDetection, include = TRUE, out.width='70%', fig.width = 9, fig.height = 4.5}
# Perform edge detection by erosion # Perform edge detection by erosion
edgeInside <- lidaRtRee::edgeDetection(gaps1r$gap.id > 0) edge_inside <- lidaRtRee::edge_detection(gaps1r$gap_id > 0)
# Perform edge detection by dilation # Perform edge detection by dilation
edgeOutside <- lidaRtRee::edgeDetection(gaps1r$gap.id > 0, inside = FALSE) edge_outside <- lidaRtRee::edge_detection(gaps1r$gap_id > 0, inside = FALSE)
# display zoom on edges # display zoom on edges
par(mfrow = c(1, 2)) par(mfrow = c(1, 2))
raster::plot(edgeInside, main = "Edges (detection by erosion)", legend = FALSE) raster::plot(edge_inside, main = "Edges (detection by erosion)", legend = FALSE)
raster::plot(edgeOutside, main = "Edges (detection by dilation)", legend = FALSE) raster::plot(edge_outside, main = "Edges (detection by dilation)", legend = FALSE)
``` ```
Percentage of edges can be calculated as the ratio between edge pixels surface and total surface, multiplied by 100. The result is dependent on the raster resolution. Obtained percentage for the "erosion" method is `r round(sum(raster::values(edgeInside))/(nrow(edgeInside)*ncol(edgeInside))*100, 1)`, whereas it is `r round(sum(raster::values(edgeOutside))/(nrow(edgeOutside)*ncol(edgeOutside))*100, 1)` for method "dilation". Percentage of edges can be calculated as the ratio between edge pixels surface and total surface, multiplied by 100. The result is dependent on the raster resolution. Obtained percentage for the "erosion" method is `r round(sum(raster::values(edge_inside))/(nrow(edge_inside)*ncol(edge_inside))*100, 1)`, whereas it is `r round(sum(raster::values(edge_outside))/(nrow(edge_outside)*ncol(edge_outside))*100, 1)` for method "dilation".
## References ## References
\ No newline at end of file
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment