Forest Growth Over Time

Aberfoyle Study Area — LiDAR Increment Analysis

Author

Miguel Ibañez Alvarez

Published

28 August 2026

Forest Research

Show code 
library(tidyverse)
library(sf)

Overview

This document visualises how two key forest growth metrics evolved across the Aberfoyle study area between consecutive LiDAR acquisitions. For each of the five increment periods, per-cell differences are summarised as period means ± one standard deviation, providing an indication of both the typical growth signal and the spatial variability across the stand.

The two metrics plotted are:

Metric Field Units
Stand Volume increment dVol95_NA m³
Top Height increment dTopHeight95_NA m

Data source: grid/LiDAR_increments.gpkg

Data Loading

Each increment period is stored as a separate layer in the GeoPackage. We read all five layers, drop the geometry (we only need the attribute table for these plots), and combine them into a single tidy data frame.

Show code 
gpkg_file <- "grid/LiDAR_increments.gpkg"

# Read each layer and drop geometry — we only need the attribute table
inc_2002_2006 <- st_read(gpkg_file, layer = "Increment_2002_2006", quiet = TRUE) |>
  st_drop_geometry()

inc_2006_2008 <- st_read(gpkg_file, layer = "Increment_2006_2008", quiet = TRUE) |>
  st_drop_geometry()

inc_2008_2012 <- st_read(gpkg_file, layer = "Increment_2008_2012", quiet = TRUE) |>
  st_drop_geometry()

inc_2012_2021 <- st_read(gpkg_file, layer = "Increment_2012_2021", quiet = TRUE) |>
  st_drop_geometry()

inc_2021_2023 <- st_read(gpkg_file, layer = "Increment_2021_2023", quiet = TRUE) |>
  st_drop_geometry()

# Stack all periods into one tidy data frame
all_increments <- bind_rows(
  inc_2002_2006,
  inc_2006_2008,
  inc_2008_2012,
  inc_2012_2021,
  inc_2021_2023
)

glimpse(all_increments)
Rows: 267,626
Columns: 5
$ identification  <int> 232755, 366861, 351248, 351203, 369003, 153198, 351140…
$ period          <chr> "2002–2006", "2002–2006", "2002–2006", "2002–2006", "2…
$ dTopHeight95_NA <dbl> 1.785153, 1.060787, 0.831206, 0.385152, 1.154493, -5.6…
$ dVol95_NA       <dbl> 3.1417811, 0.9711802, 6.5506892, 1.7553394, -4.5116837…
$ dGYCspec95_NA   <dbl> 0.63405507, 0.08853695, -0.09590797, -0.45251957, 0.16…

Data Preparation

We add two columns used for plotting:

  • year_start and year_end: extracted from the period string, used to compute the midpoint year for the x-axis.
  • year_mid: the temporal midpoint of each period, placing each observation at the centre of its measurement window.
  • period_label: a clean label for axis ticks and legends.
Show code 
increments_tidy <- all_increments |>
  # Extract the two years from the period string (e.g. "2002–2006")
  mutate(
    year_start = str_extract(period, "^\\d{4}") |> as.integer(),
    year_end   = str_extract(period, "\\d{4}$") |> as.integer()
  ) |>
  # Midpoint year positions each period at the centre of its window
  mutate(
    year_mid     = (year_start + year_end) / 2,
    period_label = paste0(year_start, "–", year_end)
  ) |>
  # Keep period_label ordered chronologically for legend and facets
  mutate(
    period_label = fct_reorder(period_label, year_mid)
  )

# Check the period structure
increments_tidy |>
  distinct(period_label, year_start, year_end, year_mid) |>
  arrange(year_mid)
period_label year_start year_end year_mid
2002–2006 2002 2006 2004.0
2006–2008 2006 2008 2007.0
2008–2012 2008 2012 2010.0
2012–2021 2012 2021 2016.5
2021–2023 2021 2023 2022.0

Summarise by period

For each period we compute the mean and standard deviation of each metric, plus the cell count (n). The SD ribbon will span mean ± 1 SD.

Show code 
growth_summary <- increments_tidy |>
  group_by(period_label, year_mid, year_start, year_end) |>
  summarise(
    # Stand Volume
    vol_mean  = mean(dVol95_NA,       na.rm = TRUE),
    vol_sd    = sd(dVol95_NA,         na.rm = TRUE),
    vol_n     = sum(!is.na(dVol95_NA)),

    # Top Height
    th_mean   = mean(dTopHeight95_NA, na.rm = TRUE),
    th_sd     = sd(dTopHeight95_NA,   na.rm = TRUE),
    th_n      = sum(!is.na(dTopHeight95_NA)),

    .groups = "drop"
  ) |>
  # Pre-compute ribbon bounds to keep plot code readable
  mutate(
    vol_lo = vol_mean - vol_sd,
    vol_hi = vol_mean + vol_sd,
    th_lo  = th_mean  - th_sd,
    th_hi  = th_mean  + th_sd
  )

growth_summary
period_label year_mid year_start year_end vol_mean vol_sd vol_n th_mean th_sd th_n vol_lo vol_hi th_lo th_hi
2002–2006 2004.0 2002 2006 1.6121076 2.116290 11863 1.7573208 1.1418782 11863 -0.5041825 3.728398 0.6154426 2.899199
2006–2008 2007.0 2006 2008 -0.6005422 1.719078 13178 0.8207059 0.8267883 13178 -2.3196203 1.118536 -0.0060824 1.647494
2008–2012 2010.0 2008 2012 3.1264442 3.080925 59367 2.2300231 1.6888870 59367 0.0455193 6.207369 0.5411361 3.918910
2012–2021 2016.5 2012 2021 5.6088798 4.280658 82696 4.6812891 2.7644476 82696 1.3282222 9.889537 1.9168415 7.445737
2021–2023 2022.0 2021 2023 2.2792019 2.100393 100522 1.6381406 0.9371355 100522 0.1788085 4.379595 0.7010051 2.575276

Plots

Show code 
# Shared publication-quality theme applied to both plots
theme_growth <- function() {
  theme_bw(base_size = 13) +
    theme(
      plot.title       = element_text(face = "bold", size = 14),
      plot.subtitle    = element_text(size = 11, colour = "grey40"),
      plot.caption     = element_text(size = 9,  colour = "grey50", hjust = 0),
      axis.title       = element_text(size = 12),
      axis.text        = element_text(size = 11),
      panel.grid.minor = element_blank(),
      legend.position  = "none"
    )
}

Stand Volume increment over time

Show code 
ggplot(growth_summary, aes(x = year_mid)) +

  # SD uncertainty band
  geom_ribbon(
    aes(ymin = vol_lo, ymax = vol_hi),
    fill  = "#2166ac",
    alpha = 0.15
  ) +

  # Horizontal reference line at zero (no change)
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey60", linewidth = 0.6) +

  # Trend line connecting period means
  geom_line(
    aes(y = vol_mean),
    colour    = "#2166ac",
    linewidth = 1.2
  ) +

  # Period mean points
  geom_point(
    aes(y = vol_mean),
    colour = "#2166ac",
    size   = 3.5,
    shape  = 19
  ) +

  # Period labels above each point
  geom_text(
    aes(y = vol_hi + 5, label = as.character(period_label)),
    size   = 3.2,
    colour = "grey30",
    vjust  = 0
  ) +

  scale_x_continuous(
    breaks = growth_summary$year_mid,
    labels = growth_summary$year_mid,
    limits = c(2002, 2025)
  ) +

  labs(
    title    = "Stand Volume Increment Between LiDAR Acquisitions",
    subtitle = "Period mean ± 1 SD across all grid cells  |  Aberfoyle study area",
    x        = "Midpoint year of increment period",
    y        = "Volume increment  (m³)",
    caption  = "Source: LiDAR increments — grid/LiDAR_increments.gpkg\nDashed line = zero change. SD band reflects spatial variability across cells, not measurement error."
  ) +

  theme_growth()

Top Height increment over time

Show code 
ggplot(growth_summary, aes(x = year_mid)) +

  # SD uncertainty band
  geom_ribbon(
    aes(ymin = th_lo, ymax = th_hi),
    fill  = "#1a9850",
    alpha = 0.15
  ) +

  # Horizontal reference line at zero (no change)
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey60", linewidth = 0.6) +

  # Trend line connecting period means
  geom_line(
    aes(y = th_mean),
    colour    = "#1a9850",
    linewidth = 1.2
  ) +

  # Period mean points
  geom_point(
    aes(y = th_mean),
    colour = "#1a9850",
    size   = 3.5,
    shape  = 19
  ) +

  # Period labels above each point
  geom_text(
    aes(y = th_hi + 0.1, label = as.character(period_label)),
    size   = 3.2,
    colour = "grey30",
    vjust  = 0
  ) +

  scale_x_continuous(
    breaks = growth_summary$year_mid,
    labels = growth_summary$year_mid,
    limits = c(2002, 2025)
  ) +

  labs(
    title    = "Top Height Increment Between LiDAR Acquisitions",
    subtitle = "Period mean ± 1 SD across all grid cells  |  Aberfoyle study area",
    x        = "Midpoint year of increment period",
    y        = "Top Height increment  (m)",
    caption  = "Source: LiDAR increments — grid/LiDAR_increments.gpkg\nDashed line = zero change. SD band reflects spatial variability across cells, not measurement error."
  ) +

  theme_growth()

Normalised Annual Rates and Cumulative Increments

Data preparation

We extend growth_summary with two new sets of derived columns:

  • Normalised rate — period mean and SD divided by the number of years in that period, expressing growth as an annual rate (m yr⁻¹ or m³ yr⁻¹). This makes periods of different lengths directly comparable.
  • Cumulative increment — running total of the raw period means, placed at the end-year of each period. Propagated uncertainty is √(Σ SD²), which assumes periods are spatially and temporally independent.
Show code 
# Number of years for each increment period — must match chronological order
n_years_lookup <- tribble(
  ~year_start, ~year_end, ~n_years,
  2002,        2006,      4,
  2006,        2008,      2,
  2008,        2012,      4,
  2012,        2021,      9,
  2021,        2023,      2
)

growth_extended <- growth_summary |>
  # Join the year-length lookup
  left_join(n_years_lookup, by = c("year_start", "year_end")) |>

  # Normalised annual rates: divide mean and SD by the period length
  mutate(
    vol_rate_mean = vol_mean / n_years,
    vol_rate_sd   = vol_sd   / n_years,
    vol_rate_lo   = vol_rate_mean - vol_rate_sd,
    vol_rate_hi   = vol_rate_mean + vol_rate_sd,

    th_rate_mean  = th_mean  / n_years,
    th_rate_sd    = th_sd    / n_years,
    th_rate_lo    = th_rate_mean - th_rate_sd,
    th_rate_hi    = th_rate_mean + th_rate_sd
  ) |>

  # Cumulative totals placed at the end-year of each period
  arrange(year_end) |>
  mutate(
    vol_cumul     = cumsum(vol_mean),
    vol_cumul_sd  = sqrt(cumsum(vol_sd^2)),   # propagated uncertainty
    vol_cumul_lo  = vol_cumul - vol_cumul_sd,
    vol_cumul_hi  = vol_cumul + vol_cumul_sd,

    th_cumul      = cumsum(th_mean),
    th_cumul_sd   = sqrt(cumsum(th_sd^2)),
    th_cumul_lo   = th_cumul - th_cumul_sd,
    th_cumul_hi   = th_cumul + th_cumul_sd
  )

# Quick check
growth_extended |>
  select(period_label, n_years,
         vol_rate_mean, vol_rate_sd,
         vol_cumul,     vol_cumul_sd,
         th_rate_mean,  th_rate_sd,
         th_cumul,      th_cumul_sd)
period_label n_years vol_rate_mean vol_rate_sd vol_cumul vol_cumul_sd th_rate_mean th_rate_sd th_cumul th_cumul_sd
2002–2006 4 0.4030269 0.5290725 1.612108 2.116290 0.4393302 0.2854695 1.757321 1.141878
2006–2008 2 -0.3002711 0.8595390 1.011565 2.726520 0.4103529 0.4133942 2.578027 1.409775
2008–2012 4 0.7816111 0.7702312 4.138010 4.114123 0.5575058 0.4222217 4.808050 2.199955
2012–2021 9 0.6232089 0.4756286 9.746889 5.937175 0.5201432 0.3071608 9.489339 3.532984
2021–2023 2 1.1396009 1.0501967 12.026091 6.297753 0.8190703 0.4685678 11.127479 3.655160

Normalised annual volume rate

Show code 
ggplot(growth_extended, aes(x = year_mid)) +

  # SD ribbon (normalised)
  geom_ribbon(
    aes(ymin = vol_rate_lo, ymax = vol_rate_hi),
    fill  = "#2166ac",
    alpha = 0.15
  ) +

  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey60", linewidth = 0.6) +

  geom_line(
    aes(y = vol_rate_mean),
    colour    = "#2166ac",
    linewidth = 1.2
  ) +

  geom_point(
    aes(y = vol_rate_mean),
    colour = "#2166ac",
    size   = 3.5,
    shape  = 19
  ) +

  # Period labels show the interval and its length
  geom_text(
    aes(
      y     = vol_rate_hi + 0.15,
      label = paste0(period_label, "\n(÷", n_years, " yr)")
    ),
    size       = 3.0,
    colour     = "grey30",
    vjust      = 0,
    lineheight = 0.9
  ) +

  scale_x_continuous(
    breaks = growth_extended$year_mid,
    labels = growth_extended$year_mid,
    limits = c(2002, 2025)
  ) +

  labs(
    title    = "Annual Volume Increment Rate Between LiDAR Acquisitions",
    subtitle = "Period mean ± 1 SD normalised by interval length  |  Aberfoyle study area",
    x        = "Midpoint year of increment period",
    y        = "Volume increment rate  (m³ yr⁻¹)",
    caption  = "Source: LiDAR increments — grid/LiDAR_increments.gpkg\nEach raw period increment divided by the number of years in that period."
  ) +

  theme_growth()

Normalised annual top-height rate

Show code 
ggplot(growth_extended, aes(x = year_mid)) +

  geom_ribbon(
    aes(ymin = th_rate_lo, ymax = th_rate_hi),
    fill  = "#1a9850",
    alpha = 0.15
  ) +

  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey60", linewidth = 0.6) +

  geom_line(
    aes(y = th_rate_mean),
    colour    = "#1a9850",
    linewidth = 1.2
  ) +

  geom_point(
    aes(y = th_rate_mean),
    colour = "#1a9850",
    size   = 3.5,
    shape  = 19
  ) +

  geom_text(
    aes(
      y     = th_rate_hi + 0.03,
      label = paste0(period_label, "\n(÷", n_years, " yr)")
    ),
    size       = 3.0,
    colour     = "grey30",
    vjust      = 0,
    lineheight = 0.9
  ) +

  scale_x_continuous(
    breaks = growth_extended$year_mid,
    labels = growth_extended$year_mid,
    limits = c(2002, 2025)
  ) +

  labs(
    title    = "Annual Top-Height Increment Rate Between LiDAR Acquisitions",
    subtitle = "Period mean ± 1 SD normalised by interval length  |  Aberfoyle study area",
    x        = "Midpoint year of increment period",
    y        = "Top-height increment rate  (m yr⁻¹)",
    caption  = "Source: LiDAR increments — grid/LiDAR_increments.gpkg\nEach raw period increment divided by the number of years in that period."
  ) +

  theme_growth()

Cumulative volume increment

Points are placed at the end year of each period — the survey date when that total was first measurable. The ribbon shows propagated uncertainty (√Σ SD²), widening as periods accumulate.

Show code 
ggplot(growth_extended, aes(x = year_end)) +

  geom_ribbon(
    aes(ymin = vol_cumul_lo, ymax = vol_cumul_hi),
    fill  = "#2166ac",
    alpha = 0.15
  ) +

  geom_line(
    aes(y = vol_cumul),
    colour    = "#2166ac",
    linewidth = 1.2
  ) +

  geom_point(
    aes(y = vol_cumul),
    colour = "#2166ac",
    size   = 3.5,
    shape  = 19
  ) +

  # Annotate each point with its cumulative total
  geom_text(
    aes(
      y     = vol_cumul_hi + 1,
      label = paste0(year_end, "\n(+", round(vol_mean, 1), " m³)")
    ),
    size       = 3.0,
    colour     = "grey30",
    vjust      = 0,
    lineheight = 0.9
  ) +

  scale_x_continuous(
    breaks = growth_extended$year_end,
    limits = c(2004, 2025)
  ) +

  labs(
    title    = "Cumulative Stand Volume Increment Since 2002",
    subtitle = "Running total of period means  |  Aberfoyle study area",
    x        = "Survey year",
    y        = "Cumulative volume increment  (m³)",
    caption  = "Source: LiDAR increments — grid/LiDAR_increments.gpkg\nAnnotations show the raw increment added in each period. Ribbon = propagated uncertainty (√Σ SD²)."
  ) +

  theme_growth()

Cumulative top-height increment

Show code 
ggplot(growth_extended, aes(x = year_end)) +

  geom_ribbon(
    aes(ymin = th_cumul_lo, ymax = th_cumul_hi),
    fill  = "#1a9850",
    alpha = 0.15
  ) +

  geom_line(
    aes(y = th_cumul),
    colour    = "#1a9850",
    linewidth = 1.2
  ) +

  geom_point(
    aes(y = th_cumul),
    colour = "#1a9850",
    size   = 3.5,
    shape  = 19
  ) +

  geom_text(
    aes(
      y     = th_cumul_hi + 0.1,
      label = paste0(year_end, "\n(+", round(th_mean, 2), " m)")
    ),
    size       = 3.0,
    colour     = "grey30",
    vjust      = 0,
    lineheight = 0.9
  ) +

  scale_x_continuous(
    breaks = growth_extended$year_end,
    limits = c(2004, 2025)
  ) +

  labs(
    title    = "Cumulative Top-Height Increment Since 2002",
    subtitle = "Running total of period means  |  Aberfoyle study area",
    x        = "Survey year",
    y        = "Cumulative top-height increment  (m)",
    caption  = "Source: LiDAR increments — grid/LiDAR_increments.gpkg\nAnnotations show the raw increment added in each period. Ribbon = propagated uncertainty (√Σ SD²)."
  ) +

  theme_growth()