library(lidR)
library(terra)
library(dplyr)
library(sf)
library(here)
library(XML)
library(stats)
library(tibble)
library(data.table)
library(rprojroot)
<- 40
zmax <- 1
grain.size <- here::here()
project_root
# Choose LAD method: "linear" or "beer"
# Beer–Lambert conversion Notes:
# - Avoids log(0) and 1 by clipping near-extreme values
# - Use when cumulative light absorption or occlusion is relevant
# - Suitable if extinction coefficient is known or estimated from prior studies
<- "beer" # Set to "linear" or "beer"
lad_method
# Optional: extinction coefficient (used only for Beer–Lambert conversion)
<- 0.25
k_extinction
<- file.path(project_root, "data/TLS/LAD_voxDF.rds")
output_voxels <- file.path(project_root, "data/TLS/lad_array_m2m3.rds")
output_array <- file.path(project_root, "data/TLS/lad_vertical_profile.pdf")
output_profile_plot <- file.path(project_root, "data/envimet/tls_envimet_trees.pld")
output_envimet_tls_3d <- file.path(project_root, "data/envimet/als_envimet_trees.pld") output_envimet_als_3d
Using leaf area density (LAD) from TLS data in ENVI-met for 3D plants
Using leaf area density (LAD) from TLS data in ENVI-met for 3D plants.
Background and Method
This section explains the theoretical principles of leaf area density (LAD) and describes how it can be determined using terrestrial laser scanning (TLS). Leaf area density is an important parameter in environmental modeling, for example for radiation balance and microclimate simulations. It indicates the leaf area per volume (m²/m³) and is therefore a decisive factor for microclimate simulations, radiation models, and energy flows in vegetation stands.
Approach Type | Name / Description | Nature |
---|---|---|
Pulse-count based | Simple linear normalization of return counts or voxel hits | Empirical, direct |
Linear normalization | Straightforward normalization of pulse counts by voxel volume or max LAD | Empirical, basic |
Pulse-density normalization | Adjusts for occlusion and scan geometry | Semi-empirical |
Gap fraction models | Estimate LAD/LAI from canopy openness statistics | Semi-empirical |
Beer–Lambert conversion conversion | Uses exponential light attenuation to infer LAD | Physically-based |
Voxel-based inverse modeling | Optimizes 3D LAD to match observed light attenuation or reflectance | Physically-based |
Allometric / geometric reconstruction | Reconstructs crown volume and distributes LAD using QSM or shape fitting | Geometric, structural |
- Linear normalization is a practical baseline: simple, fast, and reproducible.
- Beer–Lambert conversion introduces realism via physical light attenuation.
More advanced models (e.g. voxel inverse or QSM-based) aim for higher biophysical fidelity at the cost of complexity.
The present analysis is based on TLS with a medium-range RIEGL scanner (e.g., VZ-400). This captures millions of 3D points of the vegetation structure with high angular resolution. The point cloud is divided into uniform voxels, from which the leaf area density is estimated in two ways.
Linear normalization (straightforwad)
\[
\text{LAD}_i = \frac{N_i}{N_{\max}} \cdot \text{LAD}_{\max}
\] - \(N_i\): Number of laser points in voxel \(i\)
- \(N_{\max}\): Maximum across all voxels
- \(\text{LAD}_{\max}\): Maximum LAD value from the literature (e.g., 5 m²/m³)
Beer–Lambert conversion
\[ \text{LAD}_i = -\frac{\ln\left(1 - \frac{N_i}{N_{\max}}\right)}{k \cdot \Delta z} \]
- \(k\): Extinction coefficient (typically 0.3–0.5)
- \(\Delta z\): vertical voxel height
Overall Workflow
What happens in the script?
Step | Description | Relevant Code |
---|---|---|
1. Read & Filter LAS | Load TLS data, optionally crop and clean it | readLAS() and las = filter_poi(...) |
2. Voxel Grid Setup | Set up 3D grid at defined grain.size |
passed to pixel_metrics(..., res = grain.size) |
3. Count Pulses | Count returns in each voxel height bin | pointsByZSlice() function |
4. Normalise Pulse Counts | Divide by global max (relative LAD) | in convert_to_LAD() : lad = (count / max) * LADmax |
5. Export Raster | Convert metrics to raster stack | terra::rast() from voxel_df |
6. Visualization | Plot LAD profiles | see plotting section |
7. Export to Plant3D | Exports the LAD to ENVI-met | see export section |
Implemetation
To use this ENVI-met tree modeling workflow in R, follow these steps to load and initialize the project correctly:
Clone GitHub Repo in RStudio
Option 1: RStudio GUI
- Go to File → New Project → Version Control → Git
- Enter the repository URL:
https://github.com/gisma-courses/tls-tree-climate.git - Choose a project directory and name
- Click Create Project
Option 2: Terminal
git clone https://github.com/gisma-courses/tls-tree-climate.git
Then open the cloned folder in RStudio (via .Rproj
file or “Open Project”).
Note: Make sure Git is installed and configured in
RStudio → Tools → Global Options → Git/SVN
The use of the
{here}
package depends on having a valid RStudio project. Without this, file paths may not resolve correctly.
Data Input Parameters and Paths
Set global parameters for the workflow, such as file paths, voxel resolution, and maximum LAD value for normalization.
Voxelization of TLS data
Voxelisation turns a 3D TLS point cloud into a grid of cubes (voxels), where each voxel holds structural information. The number of points per voxel is used to estimate Leaf Area Density (LAD), typically normalized relative to the voxel with the most returns.
- Each voxel = a 1×1×1 m³ cube
- Count the laser hits per voxel
- Normalize to maximum
- Multiply by a literature-based LAD_max (e.g. 5 m²/m³)
This gives a spatially distributed LAD profile suitable for further analysis or models like ENVI-met.
Conversion to LAD (m²/m³)
The conversion to LAD (Leaf Area Density, in m²/m³) from TLS-based voxel pulse counts is done using a relative normalization heuristic which is adopted as a practical approximation in voxel-based canopy structure analysis using TLS (Terrestrial Laser Scanning) data.:
For each voxel layer (e.g. pulses_2_3m
), the LAD is calculated as:
\[ \text{LAD}_{\text{voxel}} = \left( \frac{\text{pulse count in voxel}}{\text{maximum pulse count over all voxels}} \right) \times \text{LAD}_{\text{max}} \]
Where:
pulse count in voxel
= number of returns in this voxel layer (from TLS)max_pulse
= the maximum pulse count found in any voxel (used for normalization)LAD_max
= a fixed normalization constant (e.g. 5.0 m²/m³) chosen from literature or calibration
Species / Structure Type | LADₘₐₓ (m²/m³) | Source / Notes |
---|---|---|
Fagus sylvatica (European beech) | 3.5–5.5 | Calders et al. (2015), Chen et al. (2018) |
Quercus robur (English oak) | 3.0–6.0 | Hosoi & Omasa (2006), field studies with TLS voxelization |
Coniferous trees (e.g. pine) | 4.0–7.0 | Wilkes et al. (2017), higher LAD due to needle density |
Mixed broadleaf forest | 3.0–6.0 | Flynn et al. (2023), canopy averaged estimates |
Shrubs / understorey | 1.5–3.0 | Chen et al. (2018),lower vertical structure density |
Urban street trees | 2.0–4.0 | Simon et al. (2020), depending on pruning and species |
LAD values refer to maximum expected per 1 m vertical voxel. Values depend on species, seasonality, and scanning conditions.
What this means conceptually
You’re not measuring absolute LAD, but instead:
- Using the number of TLS returns per voxel as a proxy for leaf density
- Then normalization all voxels relatively to the most “leaf-dense” voxel
- The
LAD_max
defines what value the “densest” voxel should reach in terms of LAD
This is fast, simple, and works well when:
- You want relative structure across the canopy
- You don’t have absolute calibration (e.g. with destructive sampling or hemispheric photos)
Caveats and assumptions
- This approach assumes the TLS beam returns are proportional to leaf area, which is a simplification
- It’s sensitive to occlusion and TLS positioning
- The choice of
LAD_max
is crucial—common values from literature range from 3–7 m²/m³ for dense canopies
The LAD conversion in the following code is a relative, normalized mapping of TLS pulse counts to LAD values, normalized by the highest voxel return and normalized using a fixed LAD_max
. This gives a plausible LAD field usable for analysis, visualization, or simulation input (e.g. for ENVI-met).
::datatable(head(df_lad, 5)) DT
Raster Stack Representation of 3D Vegetation (Voxel-Based)
We represent 3D vegetation using a voxel-based raster stack:
- Space is divided into cubic voxels (e.g. 1 × 1 × 1 m).
- Each raster layer represents a height slice (e.g. 0–1 m, 1–2 m, …).
- Voxels store values like pulse counts or Leaf Area Density (LAD).
This 2D stack structure enables:
- Vertical profiling of vegetation per XY column.
- Layer-wise analysis (e.g. median, entropy).
- Integration with raster data like topography or irradiance.
- Use in raster-based ecological and microclimate models.
It supports both analysis and visualization of vertical structure with standard geospatial tools.
ENVI-met supports custom vegetation input via the SimplePlant method, which requires a vertical LAD profile per grid column. A raster stack derived from TLS data provides exactly this: each layer represents LAD in a specific height slice, and each XY cell corresponds to one vertical profile. This structure can be exported as CSV, ASCII rasters, or custom profile files.
For 3D vegetation parameterization in ENVI-met 5.8+, the raster stack enables preprocessing of spatially explicit LAD or LAI profiles, even if some reformatting is needed.
The raster stack also supports canopy clustering and prototyping. It allows classification of structural types, simplification of complex vegetation, and the creation of representative profiles for simulation.
Visualization
library(terra)
# In SpatRasterStack umwandeln
<- df_lad[, c("X", "Y")]
xy <- df_lad[, grep("^lad_", names(df_lad), value = TRUE)]
lad_vals
<- rast(cbind(xy, lad_vals), type = "xyz")
lad_raster plot(lad_raster)
LAD Profile Visualizations from TLS Data
The plot_lad_profiles()
function visualizes vertical leaf area density (LAD) profiles derived from voxelized TLS (terrestrial laser scanning) data. LAD represents leaf surface area per unit volume (m²/m³). The function provides three main plot styles:
1. XY Matrix Plot (plotstyle = "each_median"
)
Displays a grid of mini-profiles, each representing a 0.5 × 0.5 m (x/y) ground column.
Within each cell, a normalized vertical LAD profile is plotted:
- Y-axis (height) is normalized from 0 to 1 per column.
- X-axis shows LAD values normalized relative to the global LAD maximum.
Useful for comparing structural patterns across space.
2. Overall Median Profile (plotstyle = "all_median"
)
- Aggregates LAD values across all (x/y) locations by height bin.
- Produces a typical vertical profile using the median and smoothed with a moving average.
- Height is shown in absolute units (e.g. meters).
- Captures the dominant vertical canopy structure.
3. Single Profile (plotstyle = "single_profile"
)
- Extracts and plots the LAD profile at a specific (x, y) coordinate.
- Both LAD and height are shown in absolute units.
- Plots the true vertical structure at one location.
The matrix plot shows multiple vertical LAD profiles arranged in a grid, with each small plot corresponding to a specific spatial location. This allows the vertical vegetation structure to be viewed in relation to its position on the ground. To make the individual profiles comparable, both height and LAD values are normalized within the plot. A reference profile on the side shows the overall median LAD distribution by height, which helps interpret the scale and shape of the individual profiles.
Plot the profiles
# Option 1: Profile in each column
plot_lad_profiles(lad_df, plotstyle = "each_median")
# Option 2: Overall vertical LAD profile (median of all)
plot_lad_profiles(lad_df, plotstyle = "all_median")
# Option 3: Single profile at specified coordinates
plot_lad_profiles(lad_df, plotstyle = "single_profile", single_coords = c(57.5, -94.5))