1 Objective

The primary aim of this hands-on exercise is to generate wall-to-wall raster products describing forest structure for a selected area of interest (AOI) using airborne laser scanning (ALS) data. Specifically, participants will compute the Digital Terrain Model (DTM), Digital Surface Model (DSM), Canopy Height Model (CHM), and canopy cover maps at multiple height thresholds. The entire processing workflow is implemented in R, providing a fully reproducible and open-source environment for processing and analysing airborne LiDAR data for forest and ecological applications.

2 Methodology

The exercise follows a sequential, grid-based workflow. First, raw ALS point cloud data are read, inspected, and projected to a common coordinate reference system. Ground points are then identified and used to generate a DTM, while all points are used to compute a DSM. The CHM is derived by normalising elevations relative to the terrain surface. Finally, canopy cover is calculated on a regular grid at multiple height thresholds using both all returns and first returns to characterise vegetation structure and canopy closure across the study area.

3 Input data

Airborne laser scanning (.laz / .las)

4 Expected output

Digital Terrain Model (DTM) raster (GeoTIFF)
Digital Surface Model (DSM) raster (GeoTIFF)
Canopy Height Model (CHM) raster (GeoTIFF)
Canopy cover rasters (≥ 1 m, ≥ 2 m, ≥ 3 m; all returns and first returns) (GeoTIFF)

5 R packages

lidR: Offers tools for reading, processing, visualising, and extracting metrics from airborne LiDAR (ALS) point cloud data, especially for forestry and ecological uses (Roussel et al., 2015)
RCSF – Cloth Simulation Filter (CSF) is an airborne LiDAR (Light Detection and Ranging) ground points filtering algorithm which is based on cloth simulation (Zhang et al., 2016).
raster – Supports reading, writing, and manipulating raster spatial data used in legacy geospatial workflows (Hijmans et al., 2015)
terra – A modern framework for handling large raster and vector datasets, successor to raster (Hijmans et al., 2022)
sf – Handles vector spatial data using simple features, fully integrated with tidy R.
dplyr – Facilitates fast, readable data manipulation like filtering, grouping, and summarising.
ggplot2 – Creates layered, high-quality, customizable data visualisations.
plot3D – Provides functions for static 3D visualisations, helpful for exploring LiDAR point clouds.

6 Hands-on Steps

6.1 Loading R packages

packages <- c(
  "lidR",
  "RCSF",
  "terra",
  "sf",
  "dplyr",
  "ggplot2",
  "plot3D"
)

# install only missing packages
missing_pkgs <- packages[!packages %in% installed.packages()[, "Package"]]

if (length(missing_pkgs) > 0) {
  install.packages(missing_pkgs, dependencies = TRUE)
}

# load packages quietly
suppressMessages(
  lapply(packages, library, character.only = TRUE)
)

6.2 Loading raw LiDAR data

LiDAR point clouds are commonly distributed in several formats, including LAS, LAZ, PLY, and others. Here, users can upload LiDAR datasets by selecting the desired file path. Once loaded, key LiDAR features such as coordinate system, point density, number of points, memory usage, and acquisition type are printed to support informed processing decisions.

# ============================================================
# 1. LiDAR File Selection
# ============================================================
las_path <- file.choose()  # Select LAS/LAZ file interactively

# ============================================================
# 2. LiDAR Data Import with Fallback Reader
# ============================================================
las <- tryCatch(
  readALSLAS(las_path),
  error = function(e) {
    message("readALSLAS() failed. Using lidR::readLAS() instead.")
    lidR::readLAS(las_path)
  }
)  # Read LiDAR file with robust error handling

# ============================================================
# 3. Input Data Validation
# ============================================================
if (is.empty(las)) stop("LAS file is empty")  # Stop if no points are present

# ============================================================
# 4. Coordinate Reference System Assignment
# ============================================================
#lidR::projection(las) <- sf::st_crs(32633)  # Assign projected CRS (UTM zone 33N)

# ============================================================
# 5. LiDAR Object Inspection
# ============================================================
# ============================================================
# 6. Point Density and Point Spacing (Global Metrics)
# ============================================================
# Point density (pts/m²) represents the average number of LiDAR points
# per unit area over the entire dataset. By default, it considers ALL
# returns stored in the LAS file (first, intermediate, and last returns).

# Point spacing (m or cm) is the theoretical average distance between
# neighboring LiDAR points, derived from total area and total point count.
# It is a global summary metric and does NOT describe local point spacing
# or pulse spacing.

las
# cat("------------\n")
# A <- area(las)
# cat("Point spacing (cm):", sqrt(A/npoints(las))*100, "\n")
# cat("Point density (pts/m^2):", npoints(las)/A, "\n")

6.3 CRS transformation and Z conversion from US survey feet to meters

las <- sf::st_transform(las, sf::st_crs(5070))
las$Z <- las$Z * 0.304800609

6.4 Small 3D visualization (didactic example)

This step extracts a small circular subset of the LiDAR point cloud and displays it in 3D to visually illustrate canopy structure, height variation, and point distribution before further processing.

6.5 Ground classification using Cloth Simulation Filter (CSF)

Ground classification separates terrain points from vegetation and objects, essential for terrain models and normalising LiDAR heights. Ground points are identified using the Cloth Simulation Filter (CSF: Zhang et al., 2016), which simulates a virtual cloth over the inverted point cloud. Parameters control cloth resolution, rigidity, and convergence, adapting to terrain roughness. The classified point cloud offers a reliable ground surface for terrain and canopy analysis.

las_g <- classify_ground(las,
                         csf(sloop_smooth=FALSE, class_threshold=0.5,
                             cloth_resolution=0.5, rigidness=1L,
                             iterations=500L, time_step=0.65))

6.6 Canopy height generation

Raster surfaces from classified LiDAR points depict terrain, surface height, and canopy at 0.5 m resolution: • DTM interpolates ground points for terrain elevation without vegetation. • DSM captures maximum elevation, including vegetation and structures. • CHM subtracts DTM from DSM to show vegetation height and canopy structure.

# ============================================================
# 1. Ground Point Extraction
# ============================================================
ground <- filter_ground(las_g)  # Select ground-classified LiDAR points

# ============================================================
# 2. Digital Terrain Model (DTM) Generation
# ============================================================
DTM <- grid_terrain(ground, res=1, algorithm=tin())  # Interpolate ground surface

# ============================================================
# 3. Digital Surface Model (DSM) Generation
# ============================================================
DSM <- grid_metrics(las_g, res=1, ~max(Z))  # Compute canopy top surface

# ============================================================
# 4. Canopy Height Model (CHM) Calculation
# ============================================================
CHM <- DSM - DTM  # Derive canopy height above ground

6.7 Height normalization of LiDAR point clouds

Height normalisation is the process of transforming raw LiDAR elevation values into heights relative to the ground surface. By removing the influence of terrain variability, this step enables vegetation and object heights to be compared consistently across areas with different elevations and slopes.

Several methods can be used for height normalisation, depending on data characteristics and terrain complexity. Common approaches include:

Raster-based normalisation, where point elevations are adjusted using a Digital Terrain Model (DTM)
k-nearest-neighbour interpolation, which estimates ground height locally from nearby ground points
TIN(Triangulated Irregular Network)-based normalisation, which models the terrain as a triangulated irregular network

In this workflow, a TIN-based method is adopted because it provides a balanced and flexible representation of the terrain surface, particularly in areas with variable topography. This approach effectively preserves local ground variations while ensuring robust normalisation for canopy structure analysis.

# ============================================================
# 1. Height Normalization of LiDAR Point Cloud
# ============================================================
nlas <- normalize_height(las, tin())  # Normalize heights using TIN-based ground model

# ============================================================
# 2. Filtering of Non-Physical Points
# ============================================================
nlas <- filter_poi(nlas, Z>=0)  # Remove negative height values

6.8 Canopy cover generation

In this step, we map canopy cover across the entire area of interest using a regular grid (e.g., 3 m). Canopy cover is computed directly from the normalised LiDAR point cloud (heights above ground) by summarising points within each grid cell using pixel_metrics() from the lidR package. We calculate canopy cover at three height thresholds (≥ 1 m, ≥ 2 m, ≥ 3 m) using two complementary definitions:

• All returns canopy cover (CC*_all): the fraction of all LiDAR returns in a pixel with height ≥ the threshold. This reflects overall 3D vegetation “occupancy” and is sensitive to point density and penetration.

• First returns canopy cover (CC*_first): the fraction of first returns only (ReturnNumber == 1) with height ≥ the threshold. This emphasises the top-of-canopy / upper surface (what the laser hits first) and is often closer to how canopy closure is visually perceived from above.e.

cc_metrics <- function(Z, ReturnNumber) {
  # Handle empty pixels
  if (length(Z) == 0) {
    return(list(
      CC1_all = NA_real_, CC2_all = NA_real_, CC3_all = NA_real_,
      CC1_first = NA_real_, CC2_first = NA_real_, CC3_first = NA_real_
    ))
  }

  # ---- All returns canopy cover
  CC1_all <- mean(Z >= 1, na.rm = TRUE)
  CC2_all <- mean(Z >= 2, na.rm = TRUE)
  CC3_all <- mean(Z >= 3, na.rm = TRUE)

  # ---- First returns canopy cover
  first_idx <- (ReturnNumber == 1) & is.finite(Z)
  if (!any(first_idx)) {
    CC1_first <- NA_real_
    CC2_first <- NA_real_
    CC3_first <- NA_real_
  } else {
    Z1 <- Z[first_idx]
    CC1_first <- mean(Z1 >= 1, na.rm = TRUE)
    CC2_first <- mean(Z1 >= 2, na.rm = TRUE)
    CC3_first <- mean(Z1 >= 3, na.rm = TRUE)
  }

  list(
    CC1_all = CC1_all, CC2_all = CC2_all, CC3_all = CC3_all,
    CC1_first = CC1_first, CC2_first = CC2_first, CC3_first = CC3_first
  )
}

cc_stack <- pixel_metrics(nlas, ~cc_metrics(Z, ReturnNumber), res = 3, progress = FALSE)
names(cc_stack)

## [1] "CC1_all"   "CC2_all"   "CC3_all"   "CC1_first" "CC2_first" "CC3_first"

# "CC1_all" "CC2_all" "CC3_all" "CC1_first" "CC2_first" "CC3_first"

6.9 Visualization of raster products

This visualization displays the Digital Terrain Model (DTM), Digital Surface Model (DSM), and Canopy Height Model (CHM) side by side, allowing visual comparison of terrain elevation, surface structure, and vegetation height patterns across the study area.In addition, it displays canopy cover (CC) mapped at three height thresholds (≥ 1 m, ≥ 2 m, ≥ 3 m) using two complementary definitions: CC (all returns): proportion of all LiDAR returns within each pixel above the height threshold (captures overall vegetation occupancy through the canopy). CC (first returns): proportion of first returns only above the threshold (emphasises the upper canopy surface and canopy closure from above).

6.10 Export outputs

This step exports the final raster outputs as GeoTIFFs, including the Digital Terrain Model (DTM), Digital Surface Model (DSM), Canopy Height Model (CHM), and canopy cover maps (≥ 1 m, ≥ 2 m, ≥ 3 m) derived from both all returns and first returns, for further analysis and visualization.

out_dir <- file.path(dirname(las_path), "ForestSAT2026_01_Exercise")
dir.create(out_dir, showWarnings=FALSE)
writeRaster(DTM, file.path(out_dir,"DTM_1m.tif"), overwrite=TRUE)
writeRaster(DSM, file.path(out_dir,"DSM_1m.tif"), overwrite=TRUE)
writeRaster(CHM, file.path(out_dir,"CHM_1m.tif"), overwrite=TRUE)
# --- Canopy cover (ALL returns)
terra::writeRaster(cc_stack[["CC1_all"]], file.path(out_dir, "CC_ge1m_all.tif"),overwrite = TRUE)
terra::writeRaster(cc_stack[["CC2_all"]], file.path(out_dir, "CC_ge2m_all.tif"), overwrite =TRUE)
terra::writeRaster(cc_stack[["CC3_all"]], file.path(out_dir, "CC_ge3m_all.tif"), overwrite =TRUE)
# --- Canopy cover (FIRST returns)
terra::writeRaster(cc_stack[["CC1_first"]],file.path(out_dir, "CC_ge1m_first.tif"),overwrite = TRUE)
terra::writeRaster(cc_stack[["CC2_first"]],file.path(out_dir, "CC_ge2m_first.tif"),overwrite = TRUE)
terra::writeRaster(cc_stack[["CC3_first"]],file.path(out_dir, "CC_ge3m_first.tif"),overwrite =TRUE)

6.11 Reference

Alvites, C., Marchetti, M., Lasserre, B., & Santopuoli, G. (2022). LiDAR as a tool for assessing timber assortments: A systematic literature review. Remote Sensing, 14(18), 4466.
Alvites, C., Santopuoli, G., Maesano, M., Chirici, G., Moresi, F. V., Tognetti, R., & Lasserre, B. (2021). Unsupervised algorithms to detect single trees in a mixed-species and multilayered Mediterranean forest using LiDAR data. Canadian Journal of Forest Research, 51(12), 1766–1780.
Alvites, C., O’Sullivan, H., Francini, S., Marchetti, M., Santopuoli, G., Chirici, G., & Bazzato, E. (2025). Canopy height Mapper: A Google Earth Engine application for predicting global canopy heights combining GEDI with multi-source data. Environmental Modelling & Software, 183, 106268.
Silva, C. A., Crookston, N. L., Hudak, A. T., Vierling, L. A., Klauberg, C., & Silva, M. C. A. (2017). rLiDAR: An R package for LiDAR data processing. CRAN.
Roussel, J.-R., Auty, D., Coops, N. C., Tompalski, P., Goodbody, T. R. H., Meador, A. S., & Achim, A. (2020). lidR: An R package for analysis of Airborne Laser Scanning (ALS) data. Remote Sensing of Environment, 251, 112061.
Hijmans, R. J., Van Etten, J., Cheng, J., Mattiuzzi, M., Sumner, M., & Greenberg, J. A. (2015). raster: Geographic data analysis and modeling. CRAN.
Hijmans, R. J., Bivand, R., Forner, K., Ooms, J., Pebesma, E., & Sumner, M. D. (2022). terra: Spatial data analysis. CRAN.
Wright, M. N., Wager, S., Probst, P., & Wright, M. M. N. (2019). ranger: A fast implementation of Random Forests. CRAN.
Chen, T., He, T., Benesty, M., & Khotilovich, V. (2019). xgboost: Extreme Gradient Boosting. CRAN.
Zhang, W., Qi, J., Wan, P., Wang, H., Xie, D., Wang, X., & Yan, G. (2016). An easy-to-use airborne LiDAR data filtering method based on cloth simulation. Remote sensing, 8(6), 501.

7 Key takeaways

Canopy height and canopy cover describe complementary aspects of forest structure. Canopy height indicates vertical development, and canopy cover shows horizontal closure and gaps; together, they offer a more complete forest structure view than either metric alone.
Differences between all returns and first returns matter for canopy cover. Canopy cover from all returns is more sensitive to within-canopy density, while first-return cover better indicates canopy closure at the top; thus, both should be examined and compared.
Visual inspection is essential to detect artefacts in canopy height and cover maps.. Edge effects, and unrealistic patterns often stem from LiDAR preprocessing steps like ground classification, normalization, or gridding, emphasizing the need for quality control before analysis.

Exercise 01 – Canopy Height and Canopy Cover Generation

Cesar Alvites - Forest Biometrics, Remote Sensing and Artificial Intelligence Lab (SilvaLab). School of Forest, Fisheries and Geomatics Science, University of Florida, PO Box 110410 Gainesville, FL, USA 32611

2026-02-09