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

Initial push

parents
---
title: "Workflow for data preparation for area-based model calibration"
author: "Jean-Matthieu Monnet"
date: "Oct 12, 2020"
output:
pdf_document: default
html_document: default
bibliography: workflow.treedetection.bib
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
The code below presents a workfow to prepare data for the calibration of area-based models with airborne laser scanning data and field measurements.
Licence: LGPL-3
Source page: https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee/wikis/ABA-data-preparation
Changelog:
* Oct, 2020: checked compatibility with lidR 3.0.3
+ Jan, 2020: field names changed to english
+ Oct, 2019: added projection info / updated points clouds in rda archive
+ Aug, 2018: first version
## Field inventory
### Import tree-level inventory
96 plots of 15 m radius have been inventoried in the Quatre Montagnes area (Vercors Mountain, France). Plots are aggregated in clusters of 4. Data are provided in:
+ one file per plot for tree information ("Verc-CLUSTERID-PLOTID-ArbresTerrain.csv)
- one file per cluster for plot center location ("Verc-CLUSTERID-PiquetsTerrain.csv)
The first step is to import all tree information in a single data.frame.
```{r treeInventory, include=TRUE}
# set inventory parameters
# plot radius (m)
p.radius <- 15
# DBH threshold (cm) for inventory, trees smaller than this value are not inventoried
dbh.min <- 7.5
# list tree files by pattern matching
files.t <- as.list(dir(path="./data.ABA.model/field/",
pattern="Verc-[[:alnum:]]{2}-[[:digit:]]{1}_ArbresTerrain.csv",
full.names=TRUE))
# load content of all files with lapply
trees <- lapply(files.t, function(x)
{
# read table
dummy <- read.table(x, sep=";",header=T,stringsAsFactors = F)
# add one column with plotId from file name
cbind(dummy,data.frame(plotId=rep(substr(x,30,33),nrow(dummy))))
})
# bind elements of list into a single data.frame
trees <- do.call(rbind, trees)
# add study area info in plotId
trees$plotId <- paste0("Verc-", trees$plotId)
# change column names to english
names(trees) <- c("treeId", "poleId", "Species", "Azimuth.gr", "Slope.gr", "Diameter.cm",
"Ground.Distance.m", "Height.m", "Appearance", "Tilted", "Remark", "plotId")
head(trees, n=3)
```
Fields are:
* `treeId`: tree id in the plot
+ `poleId`: plot id inside the cluster (1 to 4)
+ `Species`: tree species abreviated as GESP (GEnus SPecies)
+ `Azimuth.gr`: azimuth in grades from the plot center to tree center
+ `Slope.gr`: slope in grades from the plot center to tree center
+ `Diameter.cm`: tree diameter at breast height (1.3 m, upslope), in cm
+ `Ground.Distance.m`: ground distance between the plot center and tree edge at breast height, in m
+ `Height.m`: tree height, in m
+ `Appearance`: 0: lying or missing, 1: live, 2: live with broken treetop, 3: dead with branches, 4: dead without branches (snag)
+ `Tilted`: 0: no, 1:yes
+ `Remark`
+ `plotId`: 0: plot id
### Import plot locations
The next step is to import plot locations.
```{r plotLocation, include=TRUE, fig.width = 4, fig.height = 6}
# list location files by pattern matching
files.p <- dir(path="./data.ABA.model/field/",
pattern="Verc-[[:alnum:]]{2}_PiquetsTerrain.csv",
full.names=TRUE)
# initialize data.frame
plots <- NULL
# load all plot files with a loop
for (i in files.p)
{
# read file
dummy <- read.table(i, sep=";",header=T,stringsAsFactors = F)
# append to data.frame and add plotId info
plots <- rbind(plots, cbind(dummy,data.frame(plotId=rep(substr(i,30,31),nrow(dummy)))))
}
# keep only necessary data in plots data.frame (remove duplicate position measurements)
plots <- plots[is.element(plots$Id, c("p1","p2","p3","p4")),]
# add plotId to clusterId
plots$clusterId <- paste0("Verc-", substr(plots$plotId,1,2))
plots$plotId <- paste(plots$clusterId, substr(plots$Id,2,2),sep="-")
# keep only coordinates and Id in data.frame
plots <- plots[,c("X", "Y", "plotId", "clusterId")]
#
plot(Y~X, data=plots, asp=1)
```
### Import meta data about meridian convergence and declination
In case tree positions have to be precisely calculated in a given projected coordinates system, it is necessary to correct magnetic azimut with magnetic declination at the time of inventory and meridian convergence at the location.
```{r metaData, include=TRUE}
# get meta data about meridian convergence and declination by importing the meta file.
meta <- read.table(file="./data.ABA.model/placettes.csv",sep=";",header=T)
# keep only required attributes
meta <- meta[,c("Id","Convergence_gr","Declinaison_gr")]
# rename attributes
names(meta) <- c("clusterId", "Convergence.gr", "Declination.gr")
# merge to add convergence and declination in plots data.frame
plots <- merge(plots, meta)
head(plots, n=3)
```
### Compute tree coordinates
Tree coordinates can be computed from the plot center coordinates and from the azimut, slope and ground distance measurements. In case ground distance is measured to the tree edge, the tree diameter has to be taken into account to compute the position of the tree center.
```{r treeCoordinates, include=TRUE}
# merge tree and plots data.frames to import plot center
trees <- merge(trees, plots[,c("X","Y","Convergence.gr","Declination.gr","plotId")],
by="plotId")
# compute projected coordinates
dummy <- lidaRtRee::polar2Projected(trees$X, trees$Y, 0, trees$Azimuth.gr/200*pi,
trees$Ground.Distance.m, trees$Slope.gr/200*pi,
trees$Declination.gr/200*pi,
trees$Convergence.gr/200*pi, trees$Diameter.cm/100)
# add coordinates to trees data.frame
trees[, c("X", "Y", "Horiz.Distance.m")] <- dummy[,c("x", "y", "d")]
#
head(trees, n=3)
```
### Check field data
Field data have to be checked to correct potential errors. Regarding all trees :
+ Horizontal distances between the plot center and the tree center have to be checked to make sure the tree is included in the 15 meter radius.
- Diameters can also be checked to avoid that trees below the DBH limit remain in the inventory.
Once values have been checked for potential writing errors, remaining outliers should be removed.
```{r checkData, include=TRUE}
# keep only trees inside plot
trees <- trees[trees$Horiz.Distance.m <=p.radius , ]
# keep trees above the DBH limit
trees <- trees[trees$Diameter.cm >= dbh.min , ]
# keep only live trees
trees <- trees[is.element(trees$Appearance, c(1,2)), ]
```
For each plot, if the slope and azimut are constant on the whole surface, then plotting the tree slope from the plot center as function of the tree azimut from the plot center should draw a sinusoid-shaped point cloud. Outliers are likely to be measurement errors. For trees located close to the plot center, slope measurement precision is lower.
Displaying the histogram of trees, colored by species is also informative to make sure no abnormal measures remain.
```{r plotData, include=TRUE, fig.width = 8, fig.height = 4}
# plot to test
plot.test <- 1
#
par (mfrow=c(1,2))
# plot slope as a function of azimut,
# symbol size proportional to horizontal distance from center
plot(Slope.gr ~ Azimuth.gr, data=trees, subset=which(trees$plotId==plots$plotId[plot.test]),
cex= Horiz.Distance.m/p.radius)
# plot diameter distribution
trees.t <- trees[trees$plotId==plots$plotId[plot.test], ]
dummy <- split(trees.t$Diameter.cm, trees.t$Species)
lidaRtRee::histStack(dummy, col=lidaRtRee::speciesColor()[names(dummy), "col"],
breaks=seq(from=dbh.min, by=5,
to=5*ceiling(max(trees.t$Diameter.cm-2.5)/5)+2.5))
legend("topright", names(dummy), fill=lidaRtRee::speciesColor()[names(dummy), "col"])
```
### Compute stand level parameters
Three stand level parameters are computed by aggregating the tree level information at the plot level:
* basal area (G.ha, m2/ha)
+ mean diameter (D.mean.cm, cm)
+ stem density (in per hectare)N, /ha)
```{r standParameters, include=TRUE, fig.width = 12, fig.height = 4}
# Basal area in m2/ha
dummy <- aggregate(Diameter.cm ~ plotId, trees,
FUN=function(x){sum(pi*(x/200)^2,na.rm=T)*(10000/(pi*p.radius^2))})
names(dummy)[2] <- "G.ha"
plots <- merge(plots, dummy)
# Stem density in m2/ha
dummy <- aggregate(Diameter.cm ~ plotId, trees,
FUN=function(x){length(x)*(10000/(pi*p.radius^2))})
names(dummy)[2] <- "N.ha"
plots <- merge(plots, dummy)
# Mean diameter in cm
dummy <- aggregate(Diameter.cm ~ plotId, trees,
FUN=function(x){mean(x)})
names(dummy)[2] <- "D.mean.cm"
plots <- merge(plots, dummy)
# Summary statistics
summary(plots[, c("G.ha", "N.ha", "D.mean.cm")])
# display histograms
par(mfrow=c(1,3))
hist(plots$G.ha, main="Basal area", xlab="(m2/ha)")
hist(plots$N.ha, main="Stem density", xlab="(/ha)")
hist(plots$D.mean.cm, main="Mean diameter", xlab="(cm)")
```
### Add stratum information
In this area, public forests are generally managed in a different way compared to private forests. This classification could be used as a stratification for subsequent analysis. Plots will be attributed to strata based on external GIS layers.
```{r stratum, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
# load GIS layer of public forests
public <- rgdal::readOGR("./data.ABA.model/", "Public4Montagnes")
# create dummy variable for spatial query
dummy <- plots[,c("X","Y")]
sp::coordinates(dummy) <- ~ X+Y
dummy@proj4string <- public@proj4string
# spatial query
plots$stratum <- sp::over(dummy, public)$FID
# label levels
levels(plots$stratum) <- c("public", "private")
plots$stratum[is.na(plots$stratum)] <- "private"
row.names(plots) <- plots$plotId
# display plots
plot(Y~X, data=plots, col=plots$stratum, asp=1)
sp::plot(public, add=TRUE)
# plots <- plots[, !is.element(names(plots), c("Convergence.gr","Declination.gr", "clusterId"))]
# save(plots, file="./data.ABA.model/plots.rda")
```
## Preparation of ALS data
### Point cloud extraction
Normalized points clouds are extracted over each plot. For the delineation of single trees, a buffer has to be added around the plot border.
```{r pointClouds, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
# code to create subset of whole tiles
# folder with laz files
# lazdir <- "/media/reseau/lessem/ProjetsCommuns/Lidar/data/38_Quatre_Montagnes/norm.laz/"
lazdir <- "./data.ABA.model/plots.norm.laz"
# make catalog of las files
cata <- lidR::catalog(lazdir)
# disable display of catalog processing
lidR::opt_progress(cata) <- FALSE
# extract las point cloud with buffer
llasn <- lidR::lasclipCircle(cata, plots$X, plots$Y, p.radius + 5)
# add "buffer" attribute equal to TRUE to points outside plot
for (i in 1:length(llasn))
{
llasn[[i]] <- lidR::lasadddata(llasn[[i]],
(llasn[[i]]$X-plots$X[i])^2 +
(llasn[[i]]$Y-plots$Y[i])^2 > p.radius^2,
"buffer")
}
# set projection info
dummy <- lapply(llasn, function(x){x@proj4string <- sp::CRS("+init=epsg:2154")})
names(llasn) <- plots$plotId
# set negative height values to 0
for (i in 1:length(llasn)) { llasn[[i]]@data$Z[llasn[[i]]@data$Z<0] <- 0}
# save plot clouds
# for (i in names(llas))
# {lidR::writeLAS(llas[[i]], file=paste0("./data.ABA.model/plots.norm.laz/", i, ".laz"))}
# save(llasn, file="./data.ABA.model/llasn.rda")
```
### Computation of terrain statistics
Terrain statistics can be extracted from the original point cloud (with altitude values).
```{r terrainStats, include=TRUE, warning=FALSE, message=FALSE, fig.width = 4, fig.height = 6}
# folder with laz files
# lazdir <- "/media/reseau/lessem/ProjetsCommuns/Lidar/data/38_Quatre_Montagnes/norm.laz/"
lazdir <- "./data.ABA.model/plots.laz"
# make catalog of las files
cata <- lidR::catalog(lazdir)
cata@processing_options$progress <- FALSE
# extract las point cloud with buffer
llas <- lidR::lasclipCircle(cata, plots$X, plots$Y, p.radius+10)
names(llas) <- plots$plotId
# save plot clouds
# for (i in names(llas))
# {lidR::writeLAS(llas[[i]], file=paste0("./data.ABA.model/plots.laz/", i, ".laz"))}
# compute terrain statistics
terrain.stats <- do.call(rbind,mapply(function(X,Y)
{
lidaRtRee::points2terrainStats(
(lidR::lasfilter(X, Classification==2))@data[,c("X","Y","Z")],
plots[Y,c("X","Y")], p.radius)},
X=llas, Y=as.list(1:length(llas)),
SIMPLIFY=FALSE))
# compute terrain stats without specifying centre
terrain.stats2 <- do.call(rbind,lapply(llas, function(x)
{lidaRtRee::points2terrainStats(
(lidR::lasfilter(x, Classification==2))@data[,c("X","Y","Z")])}))
#
head(cbind(terrain.stats, terrain.stats2), n=3)
# save(terrain.stats, file="./data.ABA.model/terrain.stats.rda")
```
---
title: "Workflow for ABA prediction model calibration"
author: "Jean-Matthieu Monnet"
date: "Oct 14, 2020"
output:
pdf_document: default
html_document: default
bibliography: workflow.treedetection.bib
---
```{r setup, include=FALSE}
library(knitr)
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")
```
---
Licence: CC-BY
Source page: https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee/wikis/ABA-model-calibration
The code below presents a workflow to calibrate prediction models for the estimation of forest parameters from ALS-derived metrics, using the area-based approach (ABA).
## Load required data
The "Quatre Montagnes" dataset from France, prepared as described in https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee/wikis/ABA-data-preparation is loaded from R archives (rda files).
### Field data
Field data should be organised in a data.frame named "plots" with at least two fields: plotID (unique plot identifier) and a forest stand parameter. A factor variable "stratum" is required for stratified models. Each line in the data.frame corresponds to a field plot.
The provided dataset includes three forest stand variables :
* basal area in m^2/ha (G.ha)
+ stem density in /ha (N.ha)
+ mean diameter at breast height in cm (D.mean.cm).
Scatterplots are presented below.
```{r loadFieldData, include=TRUE}
# load plot-level data
load(file="./data.ABA.model/plots.rda")
summary(plots)
# display forest variables
plot(plots[,c("G.ha", "N.ha", "D.mean.cm")])
```
### ALS data
Normalized ALS point clouds are loaded, as well as terrain statistics previously computed from the ALS ground points of each cloud. The dataset is loaded from R archives (rda files) previously prepared in https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee/wikis/ABA-data-preparation. Point clouds corresponding to each field plot are organized in a list of LAS objects from package lidR.
```{r loadFieldData, include=TRUE}
# list of LAS objects normalized point cloud inside plot extent
load("./data.ABA.model/llasn.rda")
# display one point cloud
lidR::plot(llasn[[1]])
#
# terrain statistics previously computed with (non-normalized) ground points inside each plot extent
load("./data.ABA.model/terrain.stats.rda")
head(terrain.stats, n=3)
# check that plots are ordered in the same way in all data.frames
# all(row.names(plots) == row.names(terrain.stats))
# all(row.names(plots) == names(llasn))
```
## Computation of ALS metrics from the point clouds
Two types of metrics can be computed :
* Point cloud metrics are directly computed from the point cloud or from the derived surface models on the whole plot extent.
+ Tree-based metrics are computed from the characteristics of trees detected in the point cloud (or in the derived surface models).
Point cloud metrics are computed with the function `cloudMetrics` of lidaRtRee package, which applies the `cloud_metrics` function of lidR package to all point clouds in the list. Default computed metrics are those proposed by the function `stdmetrics` of the lidR package. Additionnal metrics are available with the `ABAmodelMetrics` of package lidaRtRee.
```{r computeMetrics, include=TRUE}
# compute point cloud metrics
metrics <- lidaRtRee::cloudMetrics(llasn, ~lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2))
head(metrics, n=3)
```
Tree metrics rely on a preliminary detection of trees, which is performed with the `treeSegmentation` function of package lidaRtRee. For more information on tree detection, please refer to https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee/-/wikis/tree%20detection%20workflow. Tree segmentation requires point clouds or canopy height models with an additionnal buffer in order to avoid border effects when computing tree characteristics. Once trees are detected, metrics are derived with the function `stdTreeMetrics` of package lidaRtRee.
```{r computeTreeMetrics, include=TRUE}
# load point clouds with 15m buffer
# load("./data.ABA.model/llasn.buffer.Rdata")
# resolution of canopy height model
resCHM <- 0.5
# specifiy plot radius to exclude trees located outside plots
plot.radius <- 15
# compute tree metrics
tree.metrics <- lidaRtRee::cloudTreeMetrics(llasn, plots[,c("X", "Y")], plot.radius, res=resCHM, func=function(x) {lidaRtRee::stdTreeMetrics(x, area.ha=pi*plot.radius^2/10000)})
head(metrics, n=3)
```
In case terrain metrics have been computed from the cloud of ground points only, they can also be added as independant variables.
```{r bindMetrics, include=TRUE}
metrics <- cbind(metrics[plots$plotId,],tree.metrics[plots$plotId,])
# add terrain stats
# metrics <- cbind(metrics, terrain.stats[placettes$plotID, 1:3])
# print("terrain stats are not used in calibration because function to compute them on grid is not implemented yet")
```
## Model calibration
Once a dependant variable (forest parameter of interest) has been chosen, the function `ABAmodel` of package lidaRtRee is used to select the linear regression model that yields the highest adjusted-R^2, while checking that some model assumptions remain verified. A box-Cox transformation of the dependant variable can be done to normalize its distribution, or a log transformation of all variables. Model details and cross-validation statistics are available from the returned object.
```{r model.calibration, include=TRUE}
###########################
# calibrate one model : for whole dataset or subset
variable <- "G.ha"
# using only a subsample
subsample <- 1:nrow(plots)
# model calibration
modell1 <- lidaRtRee::ABAmodel(plots[subsample,variable], metrics[subsample,], transform="boxcox", nmax=4, test=c("partial_p", "vif"), plots[subsample,c("X", "Y")])
# display linear regression model and additional statistics
modell1$model
modell1$stats
```
```{r model.calibration, include=TRUE}
# plot predicted VS field values
lidaRtRee::ABAmodelPlot(modell1, variable)
# check correlation between residuals and other variables
cor(cbind(placettes[sousech,c("G.ha","N.ha","D.mean.cm")], modell1$values, terrain.stats[sousech, 1:3]))
#
#
##########################
# calibrate stratified models
variable <- "G.ha"
# stratification variable
strat <- "stratum"
#
modells <- list()
for (i in levels(placettes[,strat]))
{
sousech <- which(placettes[,strat]==i)
if (length(sousech)>0)
{
modells[[i]] <- lidaRtRee::ABAmodel(placettes[sousech, variable], metrics[sousech,], transform="log", nmax=4, test=c("partial_p", "vif"), placettes[sousech,c("X", "Y")])
}
}
rm(strat, sousech, i)
# combine list of models into single object
modell <- lidaRtRee::ABAmodelCombineStrata(modells, placettes$plotID)
modell$stats
# example of combination of stratum-specific models obtained with different transformations (requires previous computation of modellsbox for boxcox et modellslog for log)
# modell <- combine.stratified.ABA.models(list(public=modellsbox[[1]], private=modellslog[[2]]), placettes)
#
lidaRtRee::ABAmodelPlot(modell, variable)
```
---
title: "Workflow for plot coregistration with ALS data"
author: "Jean-Matthieu Monnet"
date: "Jan 28, 2021"
output:
pdf_document: default
html_document: default
bibliography: workflow.treedetection.bib
---
```{r setup, include=FALSE}
# library(knitr)
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")
```
---
This workflow compares a canopy height model (CHM) derived from airborne laser scanning (ALS) data and inventoried tree positions, and proposes a translation in plot position for better matching. It works with circular plots. The method is described in @Monnet2014. The workflow uses field and ALS data acquired in the forest of Lac des Rouges Truites. It is based on functions from packages `lidaRtRee` and `lidR`.
Licence: CC BY
Source page: https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee/wikis/forest-plot-coregistration-workflow
Changelog:
* Jan, 2021: checked compatibility with lidR 3.1.0 and lidaRtRee 3.0.0
+ Oct, 2020: checked compatibility with lidR 3.0.3
+ Feb, 2019: checked compatibility with lidR 2.0.0 and lidaRtRee 2.0.0
+ May, 2018: initial version
```{r include = FALSE}
library(lidR)
```
## 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)
load(file="./data.coregistration/plotsCoregistration.rda")
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 )
load(file="./data.coregistration/treesCoregistration.rda")
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é.
ALS data over the plots is provided as a list of LAS objects in rda file.
```{r las, include = TRUE}
# load point cloud over reference plots (list of las objects)
load(file="./data.coregistration/lasCoregistration.rda")
```
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
save(las, file="./data.coregistration/lasCoregistration.rda")
```
## 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.
```{r mask, include = TRUE, out.width = '60%', fig.dim=c(5.5, 4.5), warning="FALSE" }
# 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)
<