R files moved to folder R

bc48ab2c · Monnet Jean-Matthieu · fc1d9997 · bc48ab2c · bc48ab2c · bc48ab2c
Commit bc48ab2c authored 4 years ago by Monnet Jean-Matthieu
Hide whitespace changes
Inline Side-by-side

Showing

with 256 additions and 43 deletions
+256 -43
--- a/area-based.1.data.preparation.Rmd
+++ b/area-based.1.data.preparation.Rmd
@@ -6,7 +6,7 @@ output:
  html_document: default
  pdf_document: default
 papersize: a4
-bibliography: "./bib/bibliography.bib"
+bibliography: "../bib/bibliography.bib"
 ---
 ```{r setup, include=FALSE}
@@ -18,7 +18,7 @@ knitr::opts_chunk$set(fig.align = "center")
 The code below presents a workflow to prepare inventory data for the calibration of area-based models with airborne laser scanning data and field measurements.
-Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.1.data.preparation.Rmd)
+Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.1.data.preparation.Rmd)
 Required `R` libraries : `ggplot2`, `sf`, `ggmap`
@@ -40,7 +40,7 @@ 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/",
+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
@@ -49,7 +49,7 @@ 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))))
+  cbind(dummy,data.frame(plotId=rep(substr(x,31,34),nrow(dummy))))
 })
 # bind elements of list into a single data.frame
 trees <- do.call(rbind, trees)
@@ -97,7 +97,7 @@ The next step is to import plot coordinates.
 ```{r plotLocation, include=TRUE, fig.width = 4, fig.height = 6}
 # list location files by pattern matching
-files.p <- dir(path="./data/aba.model/field/",
+files.p <- dir(path="../data/aba.model/field/",
               pattern="Verc-[[:alnum:]]{2}_PiquetsTerrain.csv",
               full.names=TRUE)
 # initialize data.frame
@@ -108,7 +108,7 @@ 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)))))
+  plots <- rbind(plots, cbind(dummy, data.frame(plotId=rep(substr(i,31,32),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")),]
@@ -154,7 +154,7 @@ In case tree positions have to be precisely calculated in a given projected coor
 ```{r metaData, include=TRUE}
 # get meta data about meridian convergence and declination by importing the metadata file.
-meta <- read.table(file="./data/aba.model/field/plot.metadata.csv",sep=";",header=T)
+meta <- read.table(file="../data/aba.model/field/plot.metadata.csv",sep=";",header=T)
 # keep only required attributes
 meta <- meta[,c("Id","Convergence_gr","Declinaison_gr")]
 # rename attributes
@@ -191,7 +191,7 @@ Normalized points clouds are extracted over each plot. For the delineation of si
 ```{r pointClouds.Height, include=TRUE, message=FALSE, warning=FALSE}
 # folder with laz files
 # lazdir <- "/media/reseau/lessem/ProjetsCommuns/Lidar/data/38_Quatre_Montagnes/norm.laz/"
-lazdir <- "./data/aba.model/ALS/plots.norm.laz/"
+lazdir <- "../data/aba.model/ALS/plots.norm.laz/"
 # make catalog of las files
 cata <- lidR::catalog(lazdir)
 # set coordinate system
@@ -221,7 +221,7 @@ Points clouds with altitude values can be used to compute terrain statistics. Th
 ```{r pointClouds.Ground, include=TRUE, message=FALSE, warning=FALSE}
 # folder with laz files
-lazdir <- "./data/aba.model/ALS/plots.laz/"
+lazdir <- "../data/aba.model/ALS/plots.laz/"
 # make catalog of las files
 cata <- lidR::catalog(lazdir)
 # set coordinate system
@@ -391,7 +391,7 @@ ggplot(trees.t, aes(x = Diameter.cm, fill = Species)) +
 ### Trees positions and remote sensing data (by plot)
-In case remote sensing data is available, plotting trees positions, sizes and species on high-resolution remote sensing background also helps in identifying errors or incoherence. Airborne laser scanning data is particularly helpful. Please refer to the [field plot coregistration](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/coregistration.Rmd) and [tree segmentation](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/tree.detection.Rmd) tutorials for additional information.
+In case remote sensing data is available, plotting trees positions, sizes and species on high-resolution remote sensing background also helps in identifying errors or incoherence. Airborne laser scanning data is particularly helpful. Please refer to the [field plot coregistration](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/coregistration.Rmd) and [tree segmentation](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/tree.detection.Rmd) tutorials for additional information.
 The tree positions can be plotted with the canopy height model computed with the ALS data.
@@ -497,7 +497,7 @@ The following lines load the public forests layer and intersects it with the plo
 ```{r stratumExtraction, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
 # load GIS layer of public forests
-public <- sf::st_read("./data/aba.model/GIS/Public4Montagnes.shp", stringsAsFactors = TRUE, quiet = TRUE)
+public <- sf::st_read("../data/aba.model/GIS/Public4Montagnes.shp", stringsAsFactors = TRUE, quiet = TRUE)
 # set coordinate system
 sf::st_crs(public) <- 2154
 # spatial query
@@ -524,12 +524,12 @@ ggmap::ggmap(map) +
 ## Save data before next tutorial
-The following lines save the data required for the [area-based model calibration step](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.2.model.calibration.Rmd).
+The following lines save the data required for the [area-based model calibration step](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.2.model.calibration.Rmd).
 ```{r save.plots, include=TRUE}
-# save(plots, file="./data/aba.model/output/plots.rda")
+# save(plots, file="../data/aba.model/output/plots.rda")
-# save(llas.height, file="./data/aba.model/output/llas.height.rda", compress="bzip2")
+# save(llas.height, file="../data/aba.model/output/llas.height.rda", compress="bzip2")
-# save(terrain.stats, file="./data/aba.model/output/terrain.stats.rda")
+# save(terrain.stats, file="../data/aba.model/output/terrain.stats.rda")
 ```
--- a/area-based.2.model.calibration.Rmd
+++ b/area-based.2.model.calibration.Rmd
@@ -6,7 +6,7 @@ output:
  pdf_document: default
  html_document: default
 papersize: a4
-bibliography: "./bib/bibliography.bib"
+bibliography: "../bib/bibliography.bib"
 ---
 ```{r setup, include=FALSE}
@@ -19,11 +19,11 @@ knitr::opts_chunk$set(fig.align = "center")
 ---
 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). The workflow is based on functions from `R` packages `lidaRtRee` and `lidR`.
-Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.2.model.calibration.Rmd)
+Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.2.model.calibration.Rmd)
 # Load data
-The "Quatre Montagnes" dataset from France, prepared as described in the [data preparation tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.1.data.preparation.Rmd) is loaded from the R archive files located in the folder "data/aba.model/output".
+The "Quatre Montagnes" dataset from France, prepared as described in the [data preparation tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.1.data.preparation.Rmd) is loaded from the R archive files located in the folder "data/aba.model/output".
 ## Field data
@@ -39,7 +39,7 @@ Scatterplots of stand parameters are presented below, colored by ownership (gree
 ```{r loadFieldData, include=TRUE, message = FALSE, warning = FALSE}
 # load plot-level data
-load(file="./data/aba.model/output/plots.rda")
+load(file="../data/aba.model/output/plots.rda")
 summary(plots)
 # display forest variables
 plot(plots[,c("G.m2.ha", "N.ha", "D.mean.cm")],
@@ -48,11 +48,11 @@ plot(plots[,c("G.m2.ha", "N.ha", "D.mean.cm")],
 ## ALS data
-Normalized ALS point clouds extracted over each plot, as well as terrain statistics previously computed from the ALS ground points can also be prepared according to the [data preparation tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.1.data.preparation.Rmd). Point clouds corresponding to each field plot are organized in a list of LAS objects. Meta data of one LAS object are displayed below.
+Normalized ALS point clouds extracted over each plot, as well as terrain statistics previously computed from the ALS ground points can also be prepared according to the [data preparation tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.1.data.preparation.Rmd). Point clouds corresponding to each field plot are organized in a list of LAS objects. Meta data of one LAS object are displayed below.
 ```{r loadALSdData, include=TRUE, message = FALSE, warning = FALSE}
 # list of LAS objects: normalized point clouds inside plot extent
-load("./data/aba.model/output/llas.height.rda") 
+load("../data/aba.model/output/llas.height.rda") 
 # display one point cloud # lidR::plot(llasn[[1]])
 llas.height[[1]]
 ```
@@ -61,7 +61,7 @@ The first lines of the terrain statistics are displayed hereafter.
 ```{r loadTerrainStats, echo = FALSE}
 # terrain statistics previously computed with (non-normalized) ground points inside each plot extent
-load("./data/aba.model/output/terrain.stats.rda")
+load("../data/aba.model/output/terrain.stats.rda")
 head(terrain.stats[, 1:3], n=3)
 ```
@@ -90,7 +90,7 @@ round(head(point.metrics[, 1:8], n = 3),2)
 ## Tree metrics
-Tree metrics rely on a preliminary detection of trees, which is performed with the `lidaRtRee::treeSegmentation` function. For more details, please refer to the [tree detection tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/tree.detection.Rmd). Tree segmentation requires point clouds or canopy height models with an additional buffer in order to avoid border effects when computing tree characteristics. Once trees are detected, metrics are derived with the function `lidaRtRee::stdTreeMetrics`. A user-specific function can be specified to compute other metrics from the features of detected trees. Plot radius has to be specified as it is required to exclude trees detected outside of the plot, and to compute the plot surface. Tree segmentation is not relevant when the point cloud density is too low, typically below five points per m^2^. The function first computes a canopy height model which default resolution is 0.5 m, but this should be set to 1 m with low point densities.
+Tree metrics rely on a preliminary detection of trees, which is performed with the `lidaRtRee::treeSegmentation` function. For more details, please refer to the [tree detection tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/tree.detection.Rmd). Tree segmentation requires point clouds or canopy height models with an additional buffer in order to avoid border effects when computing tree characteristics. Once trees are detected, metrics are derived with the function `lidaRtRee::stdTreeMetrics`. A user-specific function can be specified to compute other metrics from the features of detected trees. Plot radius has to be specified as it is required to exclude trees detected outside of the plot, and to compute the plot surface. Tree segmentation is not relevant when the point cloud density is too low, typically below five points per m^2^. The function first computes a canopy height model which default resolution is 0.5 m, but this should be set to 1 m with low point densities.
 ```{r computeTreeMetrics, include=TRUE}
 # resolution of canopy height model (m)
@@ -301,7 +301,7 @@ lidaRtRee::ABAmodelPlot(model.ABA.stratified.mixed, main = paste0(variable, ", s
 # Save data before next tutorial
-The following lines save the data required for the [area-based mapping step](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.3.mapping.Rmd).
+The following lines save the data required for the [area-based mapping step](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.3.mapping.Rmd).
 ```{r saveModels, include=TRUE}
 # IN PROGRESS

--- a/R/area-based.3.mapping.and.inference.Rmd
+++ b/R/area-based.3.mapping.and.inference.Rmd
+---
+title: "R workflow for forest mapping and inference from ABA prediction models"
+author: "Jean-Matthieu Monnet"
+date: "`r Sys.Date()`"
+output:
+  pdf_document: default
+  html_document: default
+papersize: a4
+bibliography: "./bib/bibliography.bib"
+---
+```{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")
+```
+---
+The code below presents a workflow to map forest parameters using calibrated prediction models (see [previous tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.2.model.calibration.Rmd)) and wall-to-wall coverage of the region of interest with airborne laser scanning. This tutorial is the last step of the so-called area-based approach. It is based on functions from `R` packages `lidaRtRee` and `lidR`.
+Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/area-based.3.mapping.and.inference.Rmd)
+# Load data
+```{r old, include = FALSE}
+# workflow associated with package lidaRtRee
+# Copyright Irstea
+# Author(s): Jean-Matthieu Monnet
+# Licence: CeCill 2.1
+#######################"
+# workflow for computation of ALS metrics on LAZ dataset
+# application of prediction model to obtain forest parameter map
+# inference function to compute estimates from whole map
+#
+# resolution for metrics map
+resolution <- 25
+taille.dalle <- 1000
+################
+# metrics computation
+#
+# compute metrics if not already done
+if (1)#!file.exists("./data/map.metrics.rda"))
+{
+  # create catalog of normalized LAZ files : should be rectangular tiles, with edges aligned on multiple values of the resolution
+  cata <- lidR::catalog("/media/reseau/services/LESSEM/ProjetsCommuns/Lidar/data/73-Bauges/norm.laz.pnr")
+  # cata <- lidR::catalog("/media/data/las/4montagnes/laznorm/")
+  library(foreach)
+  doFuture::registerDoFuture()
+  future::plan("multisession", workers = 2L)
+  #
+  map.metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar%
+  {
+    centroid.x <- mean(as.numeric(cata@data[i,c("Min.X", "Max.X")]))
+    centroid.y <- mean(as.numeric(cata@data[i,c("Min.Y", "Max.Y")]))
+    origine.x <- floor(centroid.x/taille.dalle)*taille.dalle
+    origine.y <- floor(centroid.y/taille.dalle)*taille.dalle
+    a <- lidR::readLAS(cata$filename[i])
+    a <- lidR::lasfilter(a, X>=origine.x & X<origine.x+taille.dalle & Y>=origine.y & Y<origine.y+taille.dalle)
+    # set negative heights to 0
+    a@data$Z[a@data$Z<0] <- 0
+    dummy <- lidR::grid_metrics(a, lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2), res=resolution)
+    return(dummy)
+  }
+  # rbind data.frame with option "fill" to handle data.frames with 0 rows and less columns
+  names.metrics <- names(map.metrics[[1]])
+  test <- Reduce(raster::merge, map.metrics)
+  map.metrics <- do.call(raster::merge, map.metrics)
+  names(map.metrics) <- names.metrics
+} else {load("./data/map.metrics.rda")}
+#
+# ADD TERRAIN METRICS <- REQUIRES NON NORMALIZED POINT CLOUD
+# ADD TREE METRICS <- REQUIRES ADDITIONNAL CODE TO COMPUTE THEM ON GRID 
+#################
+# rasterize strata
+public <- rgdal::readOGR("./data/", "Public4Montagnes")
+map.metrics@crs <- public@proj4string
+map.metrics$strata <- raster::rasterize(public, map.metrics)
+map.metrics$strata[is.na(map.metrics$strata)] <- 0
+levels(map.metrics$strata) <- data.frame(ID=c(0,1), propriete=c("private", "public"))
+#
+################
+# requires previous calibration of models (script workflow.prediction.model.calibration.R)
+# apply stratified prediction model
+rast <- lidaRtRee::ABApredict(modell, map.metrics, strata="strata")
+#
+# load forest mask
+ForestMask <- rgdal::readOGR("./data/", "Forest_Mask")
+ForestMask <- raster::rasterize(ForestMask, rast)
+ForestMask[!is.na(ForestMask)] <- 1
+#
+# clean raster (apply lower and upper thresholds to predicted values)
+rast <- lidaRtRee::cleanPredictionRaster(rast, c(0,120), ForestMask)
+raster::plot(rast)
+# display field values
+points(modell$values[,c("x","y")], cex=modell$values$field/quantile(modell$values$field,probs=0.95)*2)
+# add strata
+sp::plot(public,add=T)
+#
+# apply non-stratified prediction model
+rast1 <- lidaRtRee::ABApredict(modell1, raster::dropLayer(map.metrics, which(names(map.metrics)=="strata")))
+# clean raster
+rast1 <- lidaRtRee::cleanPredictionRaster(rast1, c(0,120), ForestMask)
+raster::plot(rast1)
+points(modell1$values[,c("x","y")], cex=modell1$values$field/quantile(modell1$values$field,probs=0.95)*2)
+#
+#################
+# map inference for non-stratified model
+# non-forest areas are already set to NA in the predictions raster
+# r.mask raster is used to specify post-stratification classes (basal area between 0, 25, 40, 55, 120).
+inference.single <- lidaRtRee::ABA.inference(modell1, rast1, r.mask=cut(rast1,quantile(raster::values(rast1),c(0,0.25,0.5,0.75,1),na.rm=T)))
+# map inference for stratified model
+# non-forest areas are already set to NA in the predictions raster
+# r.mask raster is used to specify post-stratification classes (basal area between 0, 25, 40, 55, 120).
+inference.strat <- lidaRtRee::ABA.inference(modell, rast, r.mask=cut(rast,quantile(raster::values(rast),c(0,0.25,0.5,0.75,1),na.rm=T)))
+#
+# Mask of Public forests
+PublicMask <- raster::rasterize(public, ForestMask)
+# set values in private forests to NA
+PublicMask[PublicMask!=1] <- NA
+# Private Mask
+PrivateMask <- PublicMask
+# reverte PublicMask
+values(PrivateMask) <- ifelse(is.na(values(PrivateMask)), 1, NA)
+# set non-forest values to NA
+PrivateMask <- PrivateMask * ForestMask
+#
+# map inference for non-stratified model
+# public forests only
+inference.single.public<- lidaRtRee::ABA.inference(modell1, rast1*PublicMask, r.mask=cut(rast1*PublicMask,quantile(raster::values(rast1*PublicMask),c(0,0.25,0.5,0.75,1),na.rm=T)))
+# private forests only
+# categories set manually because there are no field observations below q0.25
+inference.single.private<- lidaRtRee::ABA.inference(modell1, rast1*PrivateMask, r.mask=cut(rast1*PrivateMask,c(40, 55, 65, 120)))
+#
+# map inference for stratified model
+# public forests only
+inference.strat.public<- lidaRtRee::ABA.inference(modell, rast*PublicMask, r.mask=cut(rast*PublicMask,quantile(raster::values(rast*PublicMask),c(0,0.25,0.5,0.75,1),na.rm=T)))
+# private forests only
+# categories set manually because there are no field observations below q0.25
+inference.strat.private<- lidaRtRee::ABA.inference(modell, rast*PrivateMask, r.mask=cut(rast*PrivateMask,c(40, 55, 65, 120)))
+```
\ No newline at end of file
--- a/coregistration.Rmd
+++ b/coregistration.Rmd
@@ -6,7 +6,7 @@ output:
  pdf_document: default
  html_document: default
 papersize: a4
-bibliography: "./bib/bibliography.bib"
+bibliography: "../bib/bibliography.bib"
 ---
 ```{r setup, include=FALSE}
@@ -20,7 +20,7 @@ knitr::opts_chunk$set(fig.align = "center")
 ---
 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).
-Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/coregistration.Rmd)
+Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/coregistration.Rmd)
 ## Material
@@ -37,7 +37,7 @@ The study area is a part of the forest of Lac des Rouges Truites. 44 plots have
 ```{r plots, include = TRUE}
 # load plot coordinates (data.frame with lines corresponding to the las objects)
-load(file="./data/coregistration/plotsCoregistration.rda")
+load(file="../data/coregistration/plotsCoregistration.rda")
 head(p, n=3L)
 ```
@@ -56,7 +56,7 @@ On each plot, five trees which were considered  suitable for coregistration (ver
 ```{r trees, include = TRUE}
 # load inventoried trees (data.frame with plot id info )
-load(file="./data/coregistration/treesCoregistration.rda")
+load(file="../data/coregistration/treesCoregistration.rda")
 head(ap, n=3L)
 ```
@@ -68,7 +68,7 @@ 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")
+load(file="../data/coregistration/lasCoregistration.rda")
 ```
 Display point cloud of plot 1.
@@ -90,7 +90,7 @@ 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")
+save(las, file="../data/coregistration/lasCoregistration.rda")
 ```
 ## Parameters

--- a/forest.structure.metrics.Rmd
+++ b/forest.structure.metrics.Rmd
@@ -6,7 +6,7 @@ output:
  html_document: default
  pdf_document: default
 papersize: a4
-bibliography: "./bib/bibliography.bib"
+bibliography: "../bib/bibliography.bib"
 ---
 ```{r setup, include=FALSE}
@@ -17,15 +17,15 @@ knitr::opts_chunk$set(fig.align = "center")
 ```
 ---
-Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/forest.structure.metrics.Rmd)
+Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/forest.structure.metrics.Rmd)
 The R code below presents a forest structure metrics computation workflow from Airborne Laser Scanning (ALS) data. Workflow is based on functions from packages `lidaRtRee` and `lidR`, packages `vegan` and `foreach` are 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.
 Different types of metrics are computed:  
 * 1D height metrics
 + 2D metrics of the canopy height model (CHM)
-+ tree segmentation metrics (see also [tree segmentation tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/tree.detection.Rmd))
+ 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/gaps.edges.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))
 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.
@@ -80,7 +80,7 @@ The first step is to create a catalog of LAS files (should be normalized, non-ov
 ```{r lascatalog, include = TRUE, fig.dim = c(3.5, 2.5), out.width='40%', warning=FALSE}
 # create catalog of LAS files 
-cata <- lidR::catalog("./data/forest.structure.metrics")
+cata <- lidR::catalog("../data/forest.structure.metrics")
 # set coordinate system
 sp::proj4string(cata) <- sp::CRS(SRS_string = "EPSG:2154")
 # set sensor type

--- a/gaps.edges.detection.Rmd
+++ b/gaps.edges.detection.Rmd
@@ -6,7 +6,7 @@ output:
  pdf_document: default
  html_document: default
 papersize: a4
-bibliography: "./bib/bibliography.bib"
+bibliography: "../bib/bibliography.bib"
 ---
 ```{r setup, include=FALSE}
@@ -19,7 +19,7 @@ knitr::opts_chunk$set(fig.align = "center")
 ---
 This tutorial presents R code for forest gaps and edges detection from Airborne Laser Scanning (ALS) data. The workflow is based on functions from packages `lidaRtRee` and `lidR`.
-Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/gap.edges.detection.Rmd)
+Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/gap.edges.detection.Rmd)
 ## Study area and ALS data
@@ -47,7 +47,7 @@ chm <- lidaRtRee::points2DSM(las, 1, centre[1]-size, centre[1]+size, centre[2]-s
 If the CHM has been previously calculated.
 ```{r loadals, include = TRUE}
 # load dataset
-load(file="./data/gap.detection/chm.rda")
+load(file="../data/gap.detection/chm.rda")
 ```
 Missing, low and high values are replaced for better visualization and processing.

--- a/tree.detection.Rmd
+++ b/tree.detection.Rmd
@@ -6,7 +6,7 @@ output:
  html_document: default
  pdf_document: default
 papersize: a4
-bibliography: "./bib/bibliography.bib"
+bibliography: "../bib/bibliography.bib"
 ---
 ```{r setup, include=FALSE}
@@ -30,7 +30,7 @@ The code below presents a tree segmentation workflow from Airborne Laser Scannin
 Steps 1 and 3 are documented in [@Monnet10; @Monnet11c]. The detection performance of this algorithm was evaluated in a benchmark [@Eysn15].
-Licence: GNU GPLv3 / [source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/tree.detection.Rmd) 
+Licence: GNU GPLv3 / [source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/tree.detection.Rmd) 
 ## Material
 ### Field inventory
@@ -535,7 +535,7 @@ Be aware that in case tree segments are vectorized, some obtained polygons might
 rm(list=ls())
 # BUILD CATALOG OF FILES
 # folder containing the files
-lazdir <- "./data/forest.structure.metrics"
+lazdir <- "../data/forest.structure.metrics"
 # build catalog
 cata <- lidR::catalog(lazdir)
 # set coordinate system
@@ -568,8 +568,7 @@ save.chm <- FALSE
 # working directory
 wd <- getwd()
 # run with two cores
-clust <- parallel::makeCluster(getOption("cl.cores", 2))
+clust <- parallel::makeCluster(getOption("cl.cores", 2), outfile = "/home/jean-matthieu/Bureau/output.txt")
-#, outfile = "/home/jean-matthieu/Bureau/output.txt")
 # export global variables to cluster because they will be called inside the function
 parallel::clusterExport(cl = clust, ls(), envir = .GlobalEnv)
 # 

--- a/R/workflow.normalisation2021.R
+++ b/R/workflow.normalisation2021.R
+# workflow associated with package lidaRtRee
+# Copyright Irstea
+# Author(s): Jean-Matthieu Monnet
+# Licence: CeCill 2.1
+###################################
+# NORMALISATION D'UN ENSEMBLE DE FICHIER LAS / LAZ AVEC lidR, méthode sans discrétisation (triangulation de Delaunay des points sol)
+#
+library(foreach)
+# number of cores to use
+doParallel::registerDoParallel(cores=4)
+# 
+# faire le catalog des fichiers las : fonctionne avec des dalles rectangulaires alignées sur les coordonnées
+# cata <- lidR::catalog("~/test_normalisation/")
+cata <- lidR::catalog("~/Bureau/NCaledonie/ncaledonie/LAS_retiled/")
+sp::plot(cata)
+#
+# au cas où les dalles ne soient pas bien alignees: refaire le dallage
+reshape_tiles <- FALSE
+tile.size <- 100
+reshape_folder <- "~/test_normalisation/reshape/"
+reshape_prefix <- "CM2015_"
+if (reshape_tiles)
+{
+  # reshape tiles
+  cataR <- catalog_reshape(cata, tile.size, reshape_folder, prefix=reshape_prefix, ext="laz")
+  # replace catalog
+  cata <- cataR
+}
+#
+# taille du tampon (m) pour eviter les effets de bord lors de la triangulation
+bsize <-5
+#
+catad <- cata@data
+#################
+# calcul parallelise du nuage normalise
+# ecriture dans un fichier suffixe par "norm"
+# renvoi un data.frame avec nom de la dalles, nb points avant, nb points après (perte possible en bordure d'acquisition si les points ne sont pas dans le convex hull des points sol)
+stats <- foreach::foreach (i=1:nrow(catad), .combine=rbind) %dopar%
+{
+  # load tile extent plus buffer
+  a <- try(lidR::clip_rectangle(cata, catad$Min.X[i]-bsize, catad$Min.Y[i]-bsize, catad$Max.X[i]+bsize, catad$Max.Y[i]+bsize)) # using only one core to load files because the for loop is already parallelised
+  if (class(a)=="try-error")
+  {
+    # print("no data")
+    temp <- data.frame(id=catad$filename[i], nb.init=NA,nb.norm=NA)
+  } else {
+    if (!is.null(a))
+    {
+      # normalize
+      a <- lidR::normalize_height(a, lidR::tin())
+      # remove buffer points
+      a <- lidR::filter_poi(a, X >= catad$Min.X[i] & Y >= catad$Min.Y[i] & X < catad$Max.X[i] & Y < catad$Max.Y[i])
+      if (nrow(a@data)>0)
+      {
+        # write normalized tile
+        # convert to las if relevant
+        # new name
+        filename <- strsplit(catad$filename[i], "/")
+        filename <- paste0(paste(filename[[1]][-length(filename[[1]])], collapse="/"), "/norm_",gsub(".las",".laz", filename[[1]][length(filename[[1]])]))
+        lidR::writeLAS(a, file=filename)
+        temp <- data.frame(id=catad$filename[i], nb.init=catad$Number.of.point.records[i], nb.norm=nrow(a@data))
+      } else {
+        temp <- data.frame(id=catad$filename[i], nb.init=0,nb.norm=0)
+      }
+    } else {
+      temp <- data.frame(id=catad$filename[i], nb.init=NA,nb.norm=NA)
+    }
+  }
+  return(temp)
+}