Terrestrial Laser Scanning (TLS) can be used to provide non-destructive estimates of the above-ground woody volume of trees. Specifically, point clouds of individual trees are “segmented” to separate their points from the points from other trees, then these points are used to fit a Quantitative Structural Model (QSM), essentially a series of connected 3D cylinders, which approximate the volume of the tree. The estimates of volume provided by the QSMs can then be combined with wood density estimates to fit local diameter-biomass allometric models. These TLS-derived biomass allometries can then be applied to all the trees in a forest site to estimate their biomass from field measured diameters.
The trouble is that the process of tree segmentation and in particular manually cleaning the segmented point cloud to remove erroneous points that don’t belong to the tree (e.g. an intersecting canopy from another tree, lianas climbing up the tree, etc.), can take a very long time. It is unfeasible to segment and clean all trees in a TLS plot. The question therefore, is how many trees, and which trees, should be segmented to produce good-enough biomass allometries.
To explore this, I did some simulations in R. First, I generated a synthetic tropical forest stand of 10,000 trees. In real tropical forests, stem size distributions often roughly follow a “reverse-J” curve: thousands of small saplings and only a handful of giant trees. We simulate this using a truncated exponential distribution (minimum diameter of 10 cm, capped at a realistic 200 cm). Using the BIOMASS package, I then assigned wood density values from a realistic range (0.4-0.8) and computed the true “population” AGB to provide a benchmark to model performance against.
library(dplyr)
library(ggplot2)
library(brms)
library(tidybayes)
library(BIOMASS)
# Simulate 10,000 stems in a plot
# using a truncated exponential distribution (reverse-J curve)
# Minimum diameter threshold = 10 cm
# Realistic cutoff < 200 cm
set.seed(123)
dbh <- rexp(10000, rate = 0.05) + 10
dbh <- dbh[dbh < 200]
# Define wood density values. Sample from realistic range
wd <- runif(length(dbh), min = 0.4, max = 0.8)
# Define spatial coordinates. Congo basin rainforest
coords <- c(1, 20)
# Estimate biomass
agb <- computeAGB(
D = dbh,
WD = wd,
coord = coords)
# Create dataframe
pop <- data.frame(dbh, wd, agb)
# Calculate stand-level total AGB
agb_total <- sum(pop$agb)
I tested four sampling strategies, each sampling 50 trees:
- The 50 largest diameter trees - large trees hold the most biomass, so maybe it’s most important to get the biomass of these trees right.
- Stratified sampling - divide the diameter range into 5 even size classes and sample 10 trees from each.
- Random - Because small trees dominate the population, this sample will inevitably be skewed towards smaller diameters.
- No big trees - stratified sampling of 50 trees but excluding the 20% largest trees. Canopy occlusion often makes it difficult to reliably fit QSMs for massive emergent trees.
sample_size <- 50
# Strategy A: Top n largest DBH
sample_big <- pop %>%
slice_max(dbh, n = sample_size) %>%
mutate(method = "Largest trees")
# Strategy B: Stratified - Equal sampling from 5 evenly spaced DBH classes
sample_stratified <- pop %>%
mutate(size_class = cut(dbh, breaks = 5)) %>%
group_by(size_class) %>%
slice_sample(n = sample_size / 5) %>%
ungroup() %>%
mutate(method = "Stratified")
# Strategy C: Randomly sampled trees.
# Will pick smaller trees on average due to their abundance
sample_random <- pop %>%
sample_n(sample_size) %>%
mutate(method = "Random")
# Strategy D: No big trees, otherwise stratified
# Large trees may be difficult to reliably fit QSMs, due to occlusion in upper canopy
sample_nobig <- pop %>%
mutate(size_class = cut(dbh, breaks = 5, labels = 1:5)) %>%
filter(size_class != 5) %>%
group_by(size_class) %>%
slice_sample(n = floor(sample_size / 4)) %>%
ungroup() %>%
mutate(method = "No big trees")
# Combine samples
sample_all <- bind_rows(sample_big, sample_stratified, sample_random, sample_nobig)
# Visualise sampling
ggplot() +
geom_point(data = pop, aes(x = dbh, y = agb), alpha = 0.5) +
geom_point(data = sample_all, aes(x = dbh, y = agb, colour = method)) +
facet_wrap(~method) +
theme_bw() +
theme(legend.position = "none") +
labs(
x = "Stem diameter (cm)",
y = "Above-ground woody biomass (Kg)"
)
I fit Bayesian regression models to predict AGB from log-transformed diameters. I included weak metabolic scaling priors, knowing that the scaling exponent for trees is generally around 2.5. I then used {tidybayes} to propagate parameter uncertainty and individual tree residual variance.
# Define weakly informative priors based on metabolic scaling theory
allo_priors <- c(
prior(normal(0, 10), class = "Intercept"),
prior(normal(2.5, 1), class = "b")
)
# Fit log-normal brms model. Predict biomass from log AGB.
mod_strat_brms <- brm(
bf(agb ~ log(dbh)),
data = sample_stratified,
family = lognormal(),
prior = allo_priors,
chains = 4,
iter = 2000,
cores = 4,
seed = 123
)
# Use update() to fit remaining models without recompiling
mod_big_brms <- update(mod_strat_brms, newdata = sample_big, seed = 123)
mod_random_brms <- update(mod_strat_brms, newdata = sample_random, seed = 123)
mod_nobig_brms <- update(mod_strat_brms, newdata = sample_nobig, seed = 123)
I used the fitted models to predict the biomass of all 10,000 trees in the stand and then summed those predictions to estimate the total stand AGB.
# Function to extract posterior predictions and sum to stand level
estimate_stand_agb_brms <- function(mod, pop_data, method_name, n_draws = 1000) {
pop_data %>%
dplyr::select(dbh) %>%
# Generate predictions incorporating both parameter and residual uncertainty
add_predicted_draws(mod, ndraws = n_draws, seed = 123) %>%
group_by(.draw) %>%
summarise(stand_total_agb = sum(.prediction), .groups = "drop") %>%
mutate(method = method_name)
}
# Generate posterior distributions for all sampling strategies
post_big <- estimate_stand_agb_brms(mod_big_brms, pop, "Largest trees")
post_strat <- estimate_stand_agb_brms(mod_strat_brms, pop, "Stratified")
post_random <- estimate_stand_agb_brms(mod_random_brms, pop, "Random")
post_nobig <- estimate_stand_agb_brms(mod_nobig_brms, pop, "No big trees")
# Combine all posterior draws
post_all <- bind_rows(post_big, post_strat, post_random, post_nobig)
# Summarize 95% CIs for plotting
summary_brms <- post_all %>%
group_by(method) %>%
summarise(
Mean = mean(stand_total_agb),
Lower_95 = quantile(stand_total_agb, 0.025),
Upper_95 = quantile(stand_total_agb, 0.975)
)
When we visualize the 95% credible intervals of the posterior predictions against the true population total, the impact of sampling design becomes clear:
# Plot estimated stand totals with 95% CIs
ggplot(summary_brms, aes(x = method, y = Mean, color = method)) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = Lower_95, ymax = Upper_95), width = 0.2, linewidth = 1) +
# True population AGB
geom_hline(yintercept = agb_total, linetype = "dashed", color = "black", linewidth = 1) +
theme_bw() +
theme(legend.position = "none") +
labs(
title = "Bayesian Stand-Level Total AGB Estimates by Sampling Strategy",
subtitle = "Dashed line indicates true population total AGB",
x = "Sampling Strategy",
y = "Total Estimated Above-Ground Biomass (kg)"
) +
scale_y_continuous(labels = scales::comma)
- Largest trees only: By only training the model on the largest trees, the model mis-understands the scaling relationship of the smaller stems, resulting in a severe underestimation of the total stand AGB and a massive uncertainty.
- Random sampling: Because random sampling captures almost no large trees, it is forced to extrapolate the growth curve out to 200 cm. It often over-estimates the biomass of the largest trees.
- No big trees: Similarly, omitting the largest size class entirely forces the model to extrapolate.
- Stratified sampling: By anchoring the regression line evenly across the structural domain of the forest, stratified sampling provides the most accurate and precise prediction of the true stand total.
When visualizing model fits on a standard arithmetic scale, the models might look surprisingly similar. But, log-transforming the axes shows that models fit on biased samples misrepresent the extreme ends of the distribution.
# Helper function to extract arithmetic space prediction intervals using tidybayes
get_brms_preds <- function(mod, new_data, method_name) {
new_data %>%
# Generates predictions in kilograms, incorporating parameter + residual uncertainty
add_predicted_draws(mod, ndraws = 1000) %>%
# Calculates the mean prediction and 95% intervals (.lower, .upper)
mean_qi(.prediction) %>%
# Rename standard tidybayes columns to match your existing ggplot script
rename(
fit = .prediction,
lwr = .lower,
upr = .upper
) %>%
mutate(method = method_name)
}
# Generate new data for prediction
new_data <- data.frame(
dbh = seq(10, 200, 1)
)
# Generate back-transformed data frames for each Bayesian model
pred_big_brms <- get_brms_preds(mod_big_brms, new_data, "Largest trees")
pred_strat_brms <- get_brms_preds(mod_strat_brms, new_data, "Stratified")
pred_random_brms <- get_brms_preds(mod_random_brms, new_data, "Random")
pred_nobig_brms <- get_brms_preds(mod_nobig_brms, new_data, "No big trees")
pred_brms_all <- bind_rows(pred_big_brms, pred_strat_brms, pred_random_brms, pred_nobig_brms)
# Plot model predictions
pred_plot <- ggplot() +
geom_point(data = pop, aes(x = dbh, y = agb), alpha = 0.3, color = "grey30") +
geom_point(data = sample_all, aes(x = dbh, y = agb, fill = method),
shape = 21, colour = "black") +
# Using the new Bayesian prediction data frame
geom_ribbon(data = pred_brms_all,
aes(x = dbh, ymin = lwr, ymax = upr, fill = method),
alpha = 0.3, color = NA) +
geom_line(data = pred_brms_all, aes(x = dbh, y = fit, colour = method), linewidth = 1) +
facet_wrap(~method) +
labs(
x = "DBH (cm)",
y = "Aboveground Biomass (kg)",
fill = "Sampling Strategy",
colour = "Sampling Strategy"
) +
theme_bw() +
theme(legend.position = "bottom")
pred_plot
pred_plot +
scale_y_continuous(trans = "log10")
Next I investigated how the number of segmented trees is likely to affect model performance. I fit the same model using the same stratified sampling approach as above, but using different numbers of trees in each class.
# Define a sequence of total sample sizes to test (must be divisible by 5 classes)
n_samples <- c(25, 50, 75, 100, 150, 200, 300)
# Initialize an empty list to store results
uncertainty_results <- list()
for (n in n_samples) {
cat("Running simulation for n =", n, "...\n")
# 1. Generate stratified sample for current 'n'
sample_n <- pop %>%
mutate(size_class = cut(dbh, breaks = 5)) %>%
group_by(size_class) %>%
slice_sample(n = n / 5) %>%
ungroup()
# 2. Update the existing brms model with the new sample
# refresh = 0 hides the MCMC sampling progress bars
mod_n <- update(mod_strat_brms, newdata = sample_n, seed = 123, refresh = 0)
# 3. Predict stand-level AGB using our previously defined function
post_n <- estimate_stand_agb_brms(mod_n, pop, method_name = as.character(n))
# 4. Calculate total uncertainty
summary_n <- post_n %>%
summarise(
sample_size = n,
Mean = mean(stand_total_agb),
Lower_95 = quantile(stand_total_agb, 0.025),
Upper_95 = quantile(stand_total_agb, 0.975),
# Calculate the CI width as a percentage of the estimated mean
# This gives us a standardized "uncertainty percentage"
CI_width_perc = ((Upper_95 - Lower_95) / Mean) * 100
)
uncertainty_results[[as.character(n)]] <- summary_n
}
# Combine all results into a single dataframe
uncertainty_df <- bind_rows(uncertainty_results)
# Absolute biomass estimates and uncertainty
ggplot(uncertainty_df, aes(x = sample_size, y = Mean)) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = Lower_95, ymax = Upper_95), width = 5, linewidth = 1) +
geom_hline(yintercept = agb_total, linetype = "dashed", color = "black") +
theme_bw() +
labs(
title = "Stand AGB Estimates by Sample Size",
subtitle = "Dashed line = true population total",
x = "Total Trees Sampled (Stratified)",
y = "Estimated Total AGB (kg)"
) +
scale_y_continuous(labels = scales::comma)
# Relative uncertainty curve
ggplot(uncertainty_df, aes(x = sample_size, y = CI_width_perc)) +
geom_line() +
geom_point(size = 2) +
theme_bw() +
labs(
title = "Relative Uncertainty vs. Sampling Effort",
subtitle = "Width of 95% CI as a percentage of the estimated mean",
x = "Total Trees Sampled (Stratified)",
y = "Uncertainty (95% CI Width as % of Mean)"
)
As expected, with a stratified sample the number of individuals doesn’t affect the mean stand-level AGB estimate that much. However, adding more trees does reduce the uncertainty on these estimates. In this simulated example, 75-100 trees might be the sweet spot beyond which adding new trees doesn’t reduce the uncertainty much more.