corrected Treecanopy cover metrics. first part of aba3.rmd

R/area-based.1.data.preparation.Rmd

@@ -217,7 +217,7 @@ for (i in 1:length(llas.height)) { llas.height[[i]]@data$Z[llas.height[[i]]@data
 #  {lidR::writeLAS(llas.height[[i]], file=paste0(lazdir, i, ".laz"))}
 ```
 
-Points clouds with altitude values can be used to compute terrain statistics. The folder `plots.laz` contains the point clouds with altitude value extracted on the extent of each field plot. No buffer is added when loading those point clouds. The points not classified as ground are removed.
+Point clouds with altitude values can be used to compute terrain statistics. The folder `plots.laz` contains the point clouds with altitude value extracted on the extent of each field plot. No buffer is added when loading those point clouds. The points not classified as ground are removed.
 
 ```{r pointClouds.Ground, include=TRUE, message=FALSE, warning=FALSE}
 # folder with laz files
@@ -491,7 +491,7 @@ hist(plots$D.mean.cm, main="Mean diameter", xlab="(cm)", ylab="Plot number")
 
 In this area, public forests are generally managed in a different way compared to private forests, resulting in different forest structures. Ownership is also linked to species composition. This classification could be used as a stratification when calibrating relationships between ALS metrics and forest parameters.
 
-Plots will be attributed to strata based on external GIS layers. The first layer is a vector map of public forests (Forêts publiques, ONF Paris, 2019.) which is available under the [Open Licence Version 2.0](https://www.etalab.gouv.fr/licence-ouverte-open-licence).
+Plots will be attributed to strata based on external GIS layers. The first layer is a vector map of public forests (Forêts publiques, ONF Paris, 2019) which is available under the [Open Licence Etalab Version 2.0](https://www.etalab.gouv.fr/licence-ouverte-open-licence).
 
 The following lines load the public forests layer and intersects it with the plot locations.
 

R/area-based.2.model.calibration.Rmd

@@ -3,8 +3,8 @@ title: "R workflow for ABA prediction model calibration"
 author: "Jean-Matthieu Monnet"
 date: "`r Sys.Date()`"
 output:
-  pdf_document: default
   html_document: default
+  pdf_document: default
 papersize: a4
 bibliography: "../bib/bibliography.bib"
 ---
@@ -84,7 +84,10 @@ Two types of metrics can be computed.
 Point cloud metrics are computed with the function `lidaRtRee::cloudMetrics`, which applies the `lidR::cloud_metrics` to all point clouds in the list. Default computed metrics are those proposed by the function [`lidR::stdmetrics`](https://github.com/Jean-Romain/lidR/wiki/stdmetrics). Additional metrics are available with the function `lidaRtRee::ABAmodelMetrics`.
 
 ```{r computeMetrics, include=TRUE}
-point.metrics <- lidaRtRee::cloudMetrics(llas.height, ~lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2))
+# define function for later use
+aba.pointMetricsFUN <- ~lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2)
+# apply function on each point cloud in list
+point.metrics <- lidaRtRee::cloudMetrics(llas.height, aba.pointMetricsFUN)
 round(head(point.metrics[, 1:8], n = 3),2)
 ```
 
@@ -94,12 +97,12 @@ Tree metrics rely on a preliminary detection of trees, which is performed with t
 
 ```{r computeTreeMetrics, include=TRUE}
 # resolution of canopy height model (m)
-resCHM <- 0.5
-# specifiy plot radius to exclude trees located outside plots
+aba.resCHM <- 0.5
+# specify plot radius to exclude trees located outside plots
 plot.radius <- 15
 # compute tree metrics
 tree.metrics <- lidaRtRee::cloudTreeMetrics(llas.height, plots[, c("X", "Y")],
-                                            plot.radius, res = resCHM,
+                                            plot.radius, res = aba.resCHM,
                                             func=function(x)
                                             {
                                               lidaRtRee::stdTreeMetrics(x,
@@ -303,6 +306,6 @@ lidaRtRee::ABAmodelPlot(model.ABA.stratified.mixed, main = paste0(variable, ", s
 
 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
+```{r saveModels, eval=FALSE}
+save(model.ABA.stratified.mixed, model.ABA, aba.pointMetricsFUN, aba.resCHM, file = "../data/aba.model/output/models.rda")
 ```
\ No newline at end of file

R/area-based.3.mapping.and.inference.Rmd

@@ -3,10 +3,10 @@ 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
+  pdf_document: default
 papersize: a4
-bibliography: "./bib/bibliography.bib"
+bibliography: "../bib/bibliography.bib"
 ---
 
 ```{r setup, include=FALSE}
@@ -19,10 +19,130 @@ 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)
+Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/area-based.3.mapping.and.inference.Rmd)
 
 # Load data
 
+## ABA prediction models
+
+ABA prediction models calibrated in the previous tutorial can be loaded from the R archive "data/aba.model/output/models.rda". Three models are available: 
+
+* a single model calibrated with the whole field dataset,
+* a stratified model which is the combination of two models: one calibrated with field plots located in private forests, and the other with plots in public forests.
+
+```{r loadModels, include = TRUE}
+load(file="../data/aba.model/output/models.rda")
+```
+
+## Airborne laser scanning data
+
+For wall-to-wall mapping, the ALS metrics that are included in the models have to be computed for each element of a partition of the ALS acquisition. The area of interest is usually divided in square pixels which surface is similar to the one of field plots used for model calibration.
+
+The workflow does not yet use the catalog processing capabilities of the package `lidR`. The batch processing of an ALS wall-to-wall thus cover relies on parallel processing of ALS point clouds organized in rectangular tiles, non-overlapping tiles. It is best if their border correspond to multiple values of the pixel size.
+For example purpose, two adjacent ALS tiles are provided , in two versions: 
+
+* one with altitude values, in the folder "../data/aba.model/ALS/tiles.laz", these will be used for terrain statistics computations,
+* one with (normalized) height values, in the folder "../data/aba.model/ALS/tiles.norm.laz", these will be used for tree detection and point metrics computation.
+
+```{r loadALS, include = TRUE}
+cata.height <- lidR::catalog("../data/aba.model/ALS/tiles.norm.laz/")
+cata.altitude <- lidR::catalog("../data/aba.model/ALS/tiles.laz/")
+sp::proj4string(cata.height) <- sp::proj4string(cata.altitude) <- sp::CRS(SRS_string = "EPSG:2154")
+cata.height
+```
+
+## GIS data
+
+In order to apply the corresponding models in the case of stratum specific calibration, a GIS layer defining the spatial extent of each category is required. The layer "Public4Montagnes" in the folder "data/aba.model/GIS" corresponds to public forests (Forêts publiques, ONF Paris, 2019). All remaining areas will be considered as private-owned. In order to exclude all non-forest areas from model application, another GIS layer will be used. It is the forest cover as defined by the [BD Forêt](https://inventaire-forestier.ign.fr/spip.php?article646) of IGN. Both layers are available under the [open licence Etalab 2.0](https://www.etalab.gouv.fr/licence-ouverte-open-licence).
+
+```{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)
+# load GIS layer of forest mask
+forest <- sf::st_read("../data/aba.model/GIS/ForestMask.shp", stringsAsFactors = TRUE, quiet = TRUE)
+# set coordinate system
+sf::st_crs(public) <- sf::st_crs(forest) <- 2154
+```
+
+# Forest parameters mapping
+
+## ALS metrics mapping
+
+First the ALS metrics used as independent variables in the ABA prediction models have to be mapped on the whole area. For this the SAME functions used for computing the metrics from the plot-extent cloud metrics should be used. They are used inside the function `lidR::grid_metrics` which performs the computation the area of interest divided in square pixels. Pixel surface should be of the same size as calibration plots. Using round values is advisable in case other remote sensing data are considered for use.
+
+```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
+# resolution for metrics map
+resolution <- 25
+```
+
+### Example on one tile
+
+The following code computes point metrics for a 100 x 100 m^2^ subsample of the first ALS tile.
+
+```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
+# load ALS tile
+point.cloud.height <- lidR::clip_rectangle(cata.height, 900100, 6447100, 900400, 6447400)
+point.cloud.ground <- lidR::filter_ground(lidR::clip_rectangle(cata.altitude, 900100, 6447100, 900400, 6447400))
+# compute point metrics
+point.metrics <- lidR::grid_metrics(point.cloud.height, ~lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2), res=resolution)
+raster::plot(point.metrics[c("")])
+```
+
+For tree metrics, three steps are required:
+
+* computation of the canopy height model, segmentation and extraction of trees;
+* summary statistics of tree features by target pixels (values calculated based on trees which apices are inside the pixel);
+* additional statistics about the percentage of cover and the mean canopy height (values calculated based on segmented tree crowns which are inside the target pixel).
+
+```{r treeMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
+# compute chm
+chm <- lidaRtRee::points2DSM(point.cloud.height, res = aba.resCHM)
+# replace NA, low and high values
+chm[is.na(chm) | chm<0 | chm>60] <- 0
+# tree top detection (same parameters as used by cloudTreeMetrics in calibration step -> default parameters)
+segms <- lidaRtRee::treeSegmentation(chm)
+# extraction to data.frame
+trees <- lidaRtRee::treeExtraction(segms$filled.dem, segms$local.maxima, segms$segments.id)
+# compute raster metrics
+metrics.trees <- lidaRtRee::rasterMetrics(trees[,-1],
+                                          res=resolution,
+                                          fun=function(x)
+                                            {
+                                            lidaRtRee::stdTreeMetrics(x, resolution^2/10000)
+                                            }, output="raster")
+# remove NAs and NaNs
+metrics.trees[!is.finite(metrics.trees)] <- 0
+#
+# compute canopy cover in trees and canopy mean height in trees
+# in region of interest, because it is not in previous step.
+r.treechm <- segms$filled.dem
+# set chm to 0 in non segment area
+r.treechm[segms$segments.id==0] <- NA
+# compute raster metrics
+dummy <- lidaRtRee::rasterMetrics(r.treechm,
+                                  res=resolution,
+                                  fun=function(x)
+                                    {
+                                    c(sum(!is.na(x$filled.dem)/(resolution/aba.resCHM)^2),
+                                      mean(x$filled.dem, na.rm=T))},
+                                  output="raster")
+names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot")
+metrics.trees <- raster::addLayer(metrics.trees, dummy)
+#
+metrics.trees
+```
+
+Terrain statistics can be computed from the altitude values of ground points.
+
+```{r terrainMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6}
+# compute terrain metrics
+f <- function(x, y, z)
+{
+  as.list(lidaRtRee::points2terrainStats(data.frame(x, y, z)))
+}
+terrain.metrics <- lidR::grid_metrics(point.cloud.ground, ~f(X, Y, Z), res=resolution)
+```
+
 ```{r old, include = FALSE}
 # workflow associated with package lidaRtRee
 # Copyright Irstea

R/forest.structure.metrics.Rmd

@@ -320,13 +320,11 @@ metrics.trees <- lidaRtRee::rasterMetrics(trees[,-1],
                                             }, output="raster")
 # remove NAs and NaNs
 metrics.trees[!is.finite(metrics.trees)] <- 0
-# remove missing variables
-metrics.trees <- raster::dropLayer(metrics.trees, c("Tree.CanopyCover", "Tree.meanCanopyH"))
 #
-# compute canopy cover in trees and canopy mean height in trees, because it is not
-# correctly calculated in previous step.
+# compute canopy cover in trees and canopy mean height in trees
+# in region of interest, because it is not in previous step.
 r.treechm <- segms$filled.dem
-# set chm to 0 in non segment area
+# set chm to NA in non segment area
 r.treechm[segms$segments.id==0] <- NA
 # compute raster metrics
 dummy <- lidaRtRee::rasterMetrics(raster::crop(r.treechm,
@@ -335,10 +333,10 @@ dummy <- lidaRtRee::rasterMetrics(raster::crop(r.treechm,
                                   res=resolution,
                                   fun=function(x)
                                     {
-                                    c(sum(!is.na(x$filled.dem)),
-                                      sum(x$filled.dem,na.rm=T))/(resolution/resCHM)^2},
+                                    c(sum(!is.na(x$filled.dem))/(resolution/resCHM)^2,
+                                      mean(x$filled.dem,na.rm=T))},
                                   output="raster")
-names(dummy) <- c("Tree.CanopyCover", "Tree.meanCanopyH")
+names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot")
 #
 dummy <- raster::extend(dummy, metrics.trees)
 metrics.trees <- raster::addLayer(metrics.trees, dummy)
@@ -577,26 +575,24 @@ metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar%
                                                 output="raster")
       # remove NAs and NaNs
       metrics.trees[!is.finite(metrics.trees)] <- 0
-      # remove missing variables
-      metrics.trees <- raster::dropLayer(metrics.trees, c("Tree.CanopyCover", "Tree.meanCanopyH"))
       #
       # extend layer in case there are areas without trees
       metrics.trees <- raster::extend(metrics.trees, metrics.2dchm, values=0)
-      # compute canopy cover in trees and canopy mean height in trees, because it is not correctly
-      # calculated in previous step.
+      # compute canopy cover in trees and canopy mean height in trees
+      # in region of interest, because it is not in previous step.
       r.treechm <- segms$filled.dem
-      # set chm to 0 in non segment area
+      # set chm to NA in non segment area
       r.treechm[segms$segments.id==0] <- NA
       # compute raster metrics
       dummy <- lidaRtRee::rasterMetrics(
         raster::crop(r.treechm, raster::extent(b.box[1], b.box[3],b.box[2], b.box[4])),
         res=resolution,
         fun=function(x){
-          c(sum(!is.na(x$filled.dem)),
-            sum(x$filled.dem,na.rm=T))/(resolution/resCHM)^2
+          c(sum(!is.na(x$filled.dem))/(resolution/resCHM)^2,
+            mean(x$filled.dem,na.rm=T))
         }
         , output="raster")
-      names(dummy) <- c("Tree.CanopyCover", "Tree.meanCanopyH")
+      names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot")
       dummy <- raster::extend(dummy, metrics.2dchm, values=0)
       #
       metrics.trees <- raster::addLayer(metrics.trees, dummy)