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:
  html_document: default
  pdf_document: default
papersize: a4
bibliography: "../bib/bibliography.bib"

# erase all
cat("\014")
rm(list = ls())
# knit options
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) 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 (tested with version r packageVersion("lidaRtRee")) and lidR (tested with version r packageVersion("lidR")).

Licence: GNU GPLv3 / Source page

Many thanks to Pascal Obstétar for checking code and improvement suggestions.

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.

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 borders 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.

Files are too large to be hosted on the gitlab repository, they can be downloaded as a zip file from Google drive, and should be extracted in the folder "data/aba.model/ALS/" before proceeding with the processing. Files can be automatically downloaded thanks to the googledrive package with the following code, this requires authenticating yourself and authorizing the package to deal on your behalf with Google Drive.

# set temporary file
temp <- tempfile(fileext = ".zip")
# download file from google drive
dl <- googledrive::drive_download(googledrive::as_id("1ripo-PLZ8_IjE7rAQ2RECj-fjg1UpC5i"), 
                     path = temp, 
                     overwrite = TRUE)
# unzip to folder
out <- unzip(temp, exdir = "../data/aba.model/ALS/")
# remove temporary file
unlink(temp)

Filters are directly specified in the catalog options, in order to load only relevant points. In particular, ground and vegetation points are loaded for height catalog, whereas only ground and water points are loaded for the altitude catalog. Points with heights larger than a threshold are alos filtered out.

# build catalogs
# option to read only xyzirnc attributes (coordinates, intensity, echo order and classification) from altitude files
# build filter to keep certain classes and discard high points
filter_cata_height <- paste(c("-keep_class", class_points, "-drop_z_above", h_points), collapse = " ")
cata_height <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.norm.laz/",
                                       progress = FALSE,
                                       select = "xyzirnc",
                                       filter = filter_cata_height)
# option to read only xyzc attributes (coordinates and classification) from height files
# option to read only ground and water points from altitude files
cata_altitude <- lidR::readALSLAScatalog("../data/aba.model/ALS/tiles.laz/",
                                         progress = FALSE,
                                         select = "xyzc",
                                         filter = "-keep_class 2 9")
# set projection info
lidR::projection(cata_height) <- lidR::projection(cata_altitude) <- 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 of IGN. Both layers are available under the open licence Etalab 2.0.

# 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. Metrics are computed in each cell of 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.

# resolution for metrics map
resolution <- 25

The following paragraph show how to compute metrics for a 300 x 300 m^2^ subsample of the first ALS tile.

# load ALS tile
point_cloud_height <- lidR::clip_rectangle(cata_height, 899600, 6447600, 899900, 6447900)
# only ground points are read thanks to catalog option on cata_altitude
point_cloud_ground <- lidR::clip_rectangle(cata_altitude, 899600, 6447600, 899900, 6447900)

Point cloud metrics

For computation of point cloud metrics, the same function used in lidaRtRee::clouds_metrics is now supplied to lidR::pixel_metrics. Some metrics are displayed hereafter.

# compute point metrics
metrics_points <- lidR::pixel_metrics(point_cloud_height, aba_point_metrics_fun, res = resolution, pkg = "terra")

par(mfrow = c(1, 3))
for (i in c("zmean", "zsd", "imean"))
{
  terra::plot(metrics_points[[i]], main = i)
}

Tree metrics

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).

It is very important that tree segmentation parameters are the same as in the calibration steps, and that the function for aggregating tree attributes into raster metrics is the also the same as in the calibration step (except for the area which might differ slightly between the field plots and target cells).

# compute chm
chm <- lidR::rasterize_canopy(point_cloud_height, res = aba_res_chm, algorithm = lidR::p2r())
# replace NA, low and high values
chm[is.na(chm) | chm < 0] <- 0
# tree top detection (same parameters as used by clouds_tree_metrics in calibration step -> default parameters)
segms <- lidaRtRee::tree_segmentation(chm)
# extraction to data.frame
trees <- lidaRtRee::tree_extraction(segms)
# compute raster metrics
metrics_trees <- lidaRtRee::raster_metrics(trees,
  res = resolution,
  fun = function(x) {
    lidaRtRee::std_tree_metrics(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::raster_metrics(r_treechm,
  res = resolution,
  fun = function(x) {
    c(
      sum(!is.na(x$filled_dem) / (resolution / aba_res_chm)^2),
      mean(x$filled_dem, na.rm = T)
    )
  },
  output = "raster"
)
names(dummy) <- c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
metrics_trees <- c(metrics_trees, dummy)

The detected trees are plotted below, with the CHM as background.

par(mfrow = c(1, 1))
terra::plot(chm, main = "CHM and detected trees")
plot(sf::st_geometry(trees), cex = trees$h / 60, pch = 3, add = TRUE)

Some metrics maps are displayed below.

par(mfrow = c(1, 3))
for (i in c("Tree_meanH", "Tree_sdH", "Tree_density"))
{
  terra::plot(metrics_trees[[i]], main = i)
}

Terrain metrics

Terrain metrics can be computed from the altitude values of ground points. A wrapping function is required to pass lidaRtRee::terrain_points_metrics to lidR::pixel_metrics.

# compute terrain metrics
f <- function(x, y, z) {
  as.list(lidaRtRee::terrain_points_metrics(data.frame(x, y, z)))
}
metrics_terrain <- lidR::pixel_metrics(point_cloud_ground, ~ f(X, Y, Z), res = resolution)

par(mfrow = c(1, 3))
terra::plot(metrics_terrain[["altitude"]], main = "altitude")
terra::plot(metrics_terrain$azimut_gr, main = "azimut (gr)", col = rainbow(128), zlim = c(0,400))
terra::plot(metrics_terrain$slope_gr, main = "Slope (gr)")

Batch processing of tiles

For the batch processing of tiles, the catalog engine of lidR is used. It allows parallel processing with package future.

# specify to use two parallel sessions
future::plan("multisession", workers = 2L)
# remove warning when using random numbers in parallel sessions
options(future.rng.onMisuse = "ignore")

The catalog engine encompasses a buffering to handle the border effects when detecting trees at tile edges. A 10 m buffer is enough for tree metrics. One can also specify points classes to use, or apply a threshold to high points. Some parameters specified in the previous paragraphs are also required for the batch processing (output map resolution, chm resolution for tree segmentation aba_res_chm, functions for metrics computation aba_point_metrics_fun). Points classes to retain (class_points) and height threshold (h_points) have already been set as options for catalog processing.

rm(list = setdiff(ls(), c("resolution", "aba_res_chm", "aba_point_metrics_fun",
                          "cata_altitude", "cata_height", "class_points", "h_points",
                          "model_aba", "model_aba_stratified_mixed", "public", "forest")))
# processing parameters
b_size <- 10
# point filter already in catalog options
lidR::opt_filter(cata_height)
# "-drop_z_above": height threshold (m) for high points removal
# "-keep_class 0 1 2 3 4 5 12" points classes to retain for analysis (vegetation + ground)
# ground should be 2
# here 12 code is also vegetation

This first chunk of code computes tree and point metrics from the normalized ALS tiles.

# SET CATALOG ENGINE OPTIONS
# process by file
lidR::opt_chunk_size(cata_height) <- 0
# buffer size
lidR::opt_chunk_buffer(cata_height) <- b_size
#
# FUNCTION to apply to chunks
routine <-
  function(chunk,
           res_map = resolution,
           points_class = class_points,
           points_h = h_points,
           res_chm = aba_res_chm,
           aba_fun = aba_point_metrics_fun) {
    # extract bounding box (without buffer)
    b_box <- sf::st_bbox(chunk)
    # load chunk (including buffer)
    a <- lidR::readLAS(chunk)
    # return NULL if empty
    if (lidR::is.empty(a))
      return(NULL)
    # add 'buffer' flag to points in buffer with TRUE value in this new attribute
    # the buffer flag is already present in chunk but it seems points lying on 
    # the northern and eastern borders are not considered "buffer"
    a <- lidR::add_attribute(a,
                             a$X < b_box[1] |
                               a$Y < b_box[2] |
                               a$X >= b_box[3] | a$Y >= b_box[4],
                             "buffer")
    # remove unwanted points: not necessary if specified in catalog options
    # a <-
    #   lidR::filter_poi(a, is.element(Classification, points_class) &
    #                      Z <= points_h)
    # check number of remaining points
    if (lidR::is.empty(a)) return(NULL)
    # set negative heights to 0
    a$Z[a$Z < 0] <- 0
    # compute chm
    chm <-
      lidR::rasterize_canopy(a, res = aba_res_chm, algorithm = lidR::p2r())
    # replace NA by 0
    # there should be no values below 0 and above h_points
    chm[is.na(chm)] <- 0
    #
    # compute tree metrics
    # tree top detection (default parameters)
    segms <- lidaRtRee::tree_segmentation(chm)
    # extraction to data.frame
    trees <- lidaRtRee::tree_extraction(segms)
    # if no trees return null
    if (is.null(trees)) return(NULL)
    # remove trees outside of tile
    # no trees should be on a tile border so st_crop can be used
    trees <- sf::st_crop(trees, sf::st_bbox(chunk))
    # if no trees return null
    if (nrow(trees) == 0) return(NULL)
    # compute raster tree metrics
    metrics_trees <- lidaRtRee::raster_metrics(
      trees,
      res = res_map,
      fun = function(x) {
        lidaRtRee::std_tree_metrics(x, res_map ^ 2 / 10000)
      },
      output = "raster"
    )
    # 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 NA in non segment area
    r_treechm[segms$segments_id == 0] <- NA
    # crop to bbox
    r_treechm <- terra::crop(r_treechm, terra::ext(chunk))
    # compute raster metrics in bbox
    metrics_trees2 <- lidaRtRee::raster_metrics(
      r_treechm,
      res = res_map,
      fun = function(x) {
        c(sum(!is.na(x$filled_dem)) / (res_map / aba_res_chm) ^ 2,
          mean(x$filled_dem, na.rm = T))
      },
      output = "raster"
    )
    names(metrics_trees2) <-
      c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
    #
    # compute 1D height metrics
    # remove buffer points
    a <- lidR::filter_poi(a, buffer == 0)
    #
    if (lidR::is.empty(a)) return(NULL)
    # all points metrics
    metrics_points <- lidR::pixel_metrics(a, aba_fun, res = res_map)
    #
    if (is.null(metrics_points)) return(NULL)
    # extend / crop to match bbox
    b_box <- terra::ext(chunk)
    metrics_trees <- terra::extend(metrics_trees, b_box)
    metrics_trees2 <- terra::extend(metrics_trees2, b_box)
    metrics_points <- terra::extend(metrics_points, b_box)
    metrics_trees <- terra::crop(metrics_trees, b_box)
    metrics_trees2 <- terra::crop(metrics_trees2, b_box)
    metrics_points <- terra::crop(metrics_points, b_box)
    # merge rasters
    metrics <- c(metrics_points, metrics_trees, metrics_trees2)
    #
    return(metrics)
  }
#
# APPLY TO CATALOG, specify to merge results
metrics_map <- lidR::catalog_apply(cata_height, routine, .options = list(automerge = TRUE))

The second chunk of code computes terrain metrics from the ALS tiles with altitude values.

# SET CATALOG ENGINE OPTIONS
# process by file
lidR::opt_chunk_size(cata_altitude) <- 0
# no buffer additional buffer
lidR::opt_chunk_buffer(cata_altitude) <- 0
# enable auto-merge of results
lidR::opt_merge(cata_altitude) <- TRUE
# 
# FUNCTION to compute terrain metrics to input to pixel_metrics
f <- function(x, y, z) {
  as.list(lidaRtRee::terrain_points_metrics(data.frame(x, y, z)))
}
# APPLY TO CATALOG
metrics_terrain <- lidR::pixel_metrics(cata_altitude, ~ f(X, Y, Z), res = resolution)
# remove unused layer
metrics_terrain <- terra::subset(metrics_terrain, which(names(metrics_terrain) != "adjR2_plane"))

Finally all metrics are combined in a single object.

metrics_terrain <- terra::extend(metrics_terrain, metrics_map)
metrics_terrain <- terra::crop(metrics_terrain, metrics_map)
# merge rasters
metrics_map <- c(metrics_map, metrics_terrain)

If the metrics maps are to be displayed in external software, the "tif" format can be used to produce one single with all bands. Meanwhile band names are not retained, so it is useful to save them in a separated R archive.

metrics_names <- names(metrics_map)
save(metrics_names, file = "../data/aba.model/output/metrics_names.rda")
terra::writeRaster(metrics_map, file = "../data/aba.model/output/metrics_map.tif", overwrite = TRUE)
save(metrics_map, file = "../data/aba.model/output/metrics_map.rda")

Forest parameters mapping

Single (not stratified) model

Forest parameters are then mapped using the mapped metrics and the previously calibrated model with the function lidaRtRee::aba_predict.

prediction_map <- lidaRtRee::aba_predict(model_aba, metrics_map)
terra::plot(prediction_map, main = "prediction map")

Stratified models

In case strata-specific models have been calibrated, it is necessary to add a layer containing the strata information in the metrics raster. For this, the existing vector layer of public forests is rasterized and added to the metrics raster. The levels in the layer metadata are filled, both in the ID (code: numeric) and stratum (label: character) fields.

# rasterize layer of public forest using value = 1 in polygons
metrics_map$stratum <- terra::rasterize(terra::vect(public), metrics_map, field = 1)
# set to 0 in non-public forests areas
metrics_map$stratum[is.na(metrics_map$stratum)] <- 0
# label levels in the raster metadata
levels(metrics_map$stratum) <- data.frame(id = c(0, 1), stratum = c("private", "public"))
# crop public vector
public_cropped <- sf::st_crop(public, terra::ext(metrics_map))

Then the function lidaRtRee::aba_predict is used to produce the map. If the input is an object containing the stratified model and if the metrics map contains the strata layer, then the model directly maps the output by retaining the corresponding model in each stratum.

# map with stratified model
prediction_map_mixed <- lidaRtRee::aba_predict(model_aba_stratified_mixed, metrics_map, stratum = "stratum")
# for checking purposes
# extracted "private" stratum model from combined model
model_private <- list(model = model_aba_stratified_mixed$model$private, stats = model_aba_stratified_mixed$stats["private", ])
# produce corresponding map
prediction_map_private <- lidaRtRee::aba_predict(model_private, metrics_map)
# extracted "public" stratum model from combined model
model_public <- list(model = model_aba_stratified_mixed$model$public, stats = model_aba_stratified_mixed$stats["public", ])
# produce corresponding map
prediction_map_public <- lidaRtRee::aba_predict(model_public, metrics_map)

# plot all maps
par(mfrow = c(2, 2))
# legend range
limits <- range(c(terra::values(prediction_map_private), terra::values(prediction_map_public)))
# plot raster and vector
terra::plot(prediction_map_private, main = "Private model", range = limits, xaxt = "n", yaxt = "n")
terra::plot(prediction_map_public, main = "Public model", range = limits, xaxt = "n", yaxt = "n")
terra::plot(metrics_map$stratum, main = "Ownership", col = c("green", "blue"), xaxt = "n", yaxt = "n")
# legend("bottom", legend = c("private", "public"), fill = c("green", "blue"))
plot(public_cropped, col = NA, add = TRUE)
terra::plot(prediction_map_mixed, main = "Stratified model", range = limits, xaxt = "n", yaxt = "n")
plot(public_cropped, col = NA, add = TRUE)

Forest mask and thresholds

To avoid applying models in non-forest areas and to remove the extremes values that may have been predicted due to outliers in the ALS point cloud and metrics values, the function lidaRtRee::clean_raster can be applied:

• it applies a lower and upper threshold to map values,
• it sets to 0 the NA values in the map,
• it then multiplies the map with a mask where non-forest cells should have an NA value for those to propagate in the prediction map.

The first step is to rasterize the forest mask. In the current example, all the area is forested so we will used the polygon with ID 27 as mask.

# rasterize the forest mask
forest_mask <- terra::rasterize(terra::vect(forest), metrics_map, field = "FID")
# set values to 1 in polygon with ID 27, 0 outside
terra::values(forest_mask) <- ifelse(terra::values(forest_mask) == 27, 1, NA)
# label levels in the raster metadata
# forest_mask <- terra::as.factor(forest_mask)
levels(forest_mask) <- data.frame(id = c(1), forest = c("forest"))
# crop forest vector
forest_cropped <- sf::st_crop(forest, terra::ext(metrics_map))
# clean raster
# threshold values and set to NA outside of forest
prediction_map_mixed_clean <- lidaRtRee::clean_raster(prediction_map_mixed, c(0, 80), forest_mask)

par(mfrow = c(1, 2))
terra::plot(forest_mask, main = "Forest mask")
plot(sf::st_geometry(forest_cropped), add = TRUE)
# terra::plot(prediction_map_mixed, main = "Stratified predictions map")
terra::plot(prediction_map_mixed_clean, main = "Masked map")
plot(sf::st_geometry(forest_cropped), add = TRUE)

Inference

Work in progress...