diff --git a/area-based.1.data.preparation.Rmd b/R/area-based.1.data.preparation.Rmd similarity index 95% rename from area-based.1.data.preparation.Rmd rename to R/area-based.1.data.preparation.Rmd index 7981a672f42e69c623a82c70cf9fdc90479030bd..cdf676e33a9d66d4c601d5242529ba9cac26dd00 100644 --- a/area-based.1.data.preparation.Rmd +++ b/R/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(llas.height, file="./data/aba.model/output/llas.height.rda", compress="bzip2") -# save(terrain.stats, file="./data/aba.model/output/terrain.stats.rda") +# save(plots, file="../data/aba.model/output/plots.rda") +# save(llas.height, file="../data/aba.model/output/llas.height.rda", compress="bzip2") +# save(terrain.stats, file="../data/aba.model/output/terrain.stats.rda") ``` diff --git a/area-based.2.model.calibration.Rmd b/R/area-based.2.model.calibration.Rmd similarity index 91% rename from area-based.2.model.calibration.Rmd rename to R/area-based.2.model.calibration.Rmd index ada4764542240b0dd4c48942e0d9a92ed7a5a31e..f070c5aac9f79f3452b18cbde3a078c8b4587b04 100644 --- a/area-based.2.model.calibration.Rmd +++ b/R/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 diff --git a/R/area-based.3.mapping.and.inference.Rmd b/R/area-based.3.mapping.and.inference.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..4f0e7bff33f1821078d73875d1e2a27cad883ef9 --- /dev/null +++ b/R/area-based.3.mapping.and.inference.Rmd @@ -0,0 +1,143 @@ +--- +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 diff --git a/coregistration.Rmd b/R/coregistration.Rmd similarity index 97% rename from coregistration.Rmd rename to R/coregistration.Rmd index 019164f4f885f758e47e45b8db92f771cdf1aad7..4cfd263472569974f482f8b8329118589846fab2 100644 --- a/coregistration.Rmd +++ b/R/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 diff --git a/forest.structure.metrics.Rmd b/R/forest.structure.metrics.Rmd similarity index 99% rename from forest.structure.metrics.Rmd rename to R/forest.structure.metrics.Rmd index 105d3fa1376b0bf003bb68153e37feeea78ddcda..1c17244241d1696dc701e5bff013c779014b6228 100644 --- a/forest.structure.metrics.Rmd +++ b/R/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)) -+ 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)) ++ 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/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 diff --git a/gaps.edges.detection.Rmd b/R/gaps.edges.detection.Rmd similarity index 98% rename from gaps.edges.detection.Rmd rename to R/gaps.edges.detection.Rmd index ad9cab03d037235b15bc499237414a96f0511746..28ef1d0b3f7644a773eb6db13aa96c609ff72617 100644 --- a/gaps.edges.detection.Rmd +++ b/R/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. diff --git a/tree.detection.Rmd b/R/tree.detection.Rmd similarity index 99% rename from tree.detection.Rmd rename to R/tree.detection.Rmd index 19aee4c9ed6c2ca3db1d0e2a0a75ea3aeae0e3e5..e278a43f48ed179927e4a1a0cd6645d2f6a55be1 100644 --- a/tree.detection.Rmd +++ b/R/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)) -#, outfile = "/home/jean-matthieu/Bureau/output.txt") +clust <- parallel::makeCluster(getOption("cl.cores", 2), 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) # diff --git a/R/workflow.normalisation2021.R b/R/workflow.normalisation2021.R new file mode 100644 index 0000000000000000000000000000000000000000..a4b1a1b51156d045325bb412e513a3d7529c5ddb --- /dev/null +++ b/R/workflow.normalisation2021.R @@ -0,0 +1,71 @@ +# 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) +}