Needleleaf Index — Conifer Detection at Aberfoyle

Threshold calibrated from local training samples

Author

Miguel Ibañez Alvarez

Published

September 2, 2026

Forest Research

Background

The Needleleaf Index (NI) uses the fact that conifer canopies reflect less energy in the near-infrared and shortwave-infrared range than broadleaf vegetation, but more than water. The index is simply the sum of those three infrared reflectance values:

\[ \text{NI} = \rho_{\text{NIR}} + \rho_{\text{SWIR1}} + \rho_{\text{SWIR2}} \]

Lower NI → darker, denser canopy (conifers); higher NI → more open or broadleaf canopy.

For Sentinel-2 surface reflectance the bands are:

Role Sentinel-2 band Wavelength
NIR B8 842 nm
SWIR1 B11 1610 nm
SWIR2 B12 2190 nm

Rather than applying thresholds from the literature (calibrated for North American boreal forests on Landsat), we derive the classification threshold directly from field-collected training samples at Aberfoyle, ensuring it reflects the local spectral conditions.

Setup

Show code 
library(terra)   # raster operations
library(sf)      # vector operations
library(dplyr)   # data manipulation
library(tidyr)   # reshaping
library(purrr)   # iteration
library(ggplot2) # visualisation

Load data

Show code 
# Sentinel-2 surface reflectance image (10 m, July 2026)
s2 <- rast("Sentinel/20260715_S2.tif")

# Training polygons — field-surveyed species labels
samples <- st_read("Shapefiles/V7_samples_aberfoyle.shp", quiet = TRUE)
Show code 
samples |>
  st_drop_geometry() |>
  count(species, name = "polygons") |>
  arrange(desc(polygons))
Table 1: Training polygons by species.
species polygons
Sitka spruce 180
Scots pine 80
Larch 72
Birch 50
Oak 50
Norway spruce 38
Other conifers 29
Open 21
Lodgepole pine 20
Other broadleaves 19

The training data cover 559 polygons across 10 classes, including both coniferous and broadleaf species as well as open ground.

Compute the Needleleaf Index

Show code 
# NI is the sum of three infrared bands
ni <- s2[["B8"]] + s2[["B11"]] + s2[["B12"]]
names(ni) <- "NI"

Extract NI values at the training locations

The training polygons are in a different coordinate system (UTM zone 30N) from the raster (British National Grid), so we reproject them first.

Show code 
# Reproject training polygons to match the raster
samples_bng <- st_transform(samples, crs = crs(ni))

# Extract NI values for every pixel inside each training polygon
ni_pixels <- terra::extract(ni, vect(samples_bng), ID = TRUE) |>
  as_tibble() |>
  # Join species label back using the polygon row index
  left_join(
    samples_bng |>
      st_drop_geometry() |>
      mutate(ID = row_number()) |>
      select(ID, species),
    by = "ID"
  ) |>
  filter(!is.na(NI))

cat("Total training pixels extracted:", nrow(ni_pixels), "\n")
Total training pixels extracted: 6320

Explore NI distributions by species

Show code 
# Order species by their median NI value so the plot is easier to read
species_order <- ni_pixels |>
  group_by(species) |>
  summarise(median_NI = median(NI)) |>
  arrange(median_NI) |>
  pull(species)

# Label each species as conifer or non-conifer for colouring
conifer_species <- c(
  "Sitka spruce", "Scots pine", "Norway spruce",
  "Lodgepole pine", "Other conifers", "Larch"
)

ni_pixels <- ni_pixels |>
  mutate(
    type    = if_else(species %in% conifer_species, "Conifer", "Non-conifer"),
    species = factor(species, levels = species_order)
  )

ggplot(ni_pixels, aes(x = NI, y = species, fill = type)) +
  geom_violin(alpha = 0.6, colour = "grey40", linewidth = 0.3) +
  geom_boxplot(width = 0.15, outlier.size = 0.5, colour = "grey20") +
  scale_fill_manual(values = c(Conifer = "#2d7a2d", `Non-conifer` = "#cc6600")) +
  labs(
    x    = "Needleleaf Index",
    y    = NULL,
    fill = NULL
  ) +
  theme_minimal() +
  theme(legend.position = "top")
Figure 1: NI distribution for each species, ordered by median value. Conifers (green) cluster at lower NI than broadleaves and open ground (orange/grey). Larch, being deciduous, sits between the two groups in the July image.

The plot shows a reasonably clear NI gap between conifers (≲ 0.50) and broadleaves / open ground (≳ 0.50), with some overlap. Larch — a deciduous conifer — has a higher summer NI than evergreen conifers because its fully expanded leaves behave spectrally more like broadleaves in July.

Derive the classification thresholds

Classifying the image needs two boundaries:

  • An upper threshold separating conifer from broadleaf / open ground (pixels above it are too bright to be conifer).
  • A lower threshold marking the bottom of the conifer NI range (pixels below it are too dark — shadows, water bodies, bare ground — to be reliably identified as conifer).

Upper threshold — conifer vs broadleaf

We find the NI cut-off that best separates conifer from non-conifer pixels in the training data. For every candidate threshold value we calculate:

  • Sensitivity — proportion of conifer pixels correctly identified
  • Specificity — proportion of non-conifer pixels correctly excluded
  • Youden’s J = sensitivity + specificity − 1 (ranges 0–1; higher is better)

The threshold with the highest Youden’s J is selected.

Show code 
# Test candidate thresholds at 0.005 intervals
candidate_thresholds <- seq(0.20, 0.80, by = 0.005)

threshold_performance <- map(candidate_thresholds, function(t) {

  predicted_conifer <- ni_pixels$NI <= t
  actual_conifer    <- ni_pixels$type == "Conifer"

  tp <- sum( predicted_conifer &  actual_conifer)   # true positives
  tn <- sum(!predicted_conifer & !actual_conifer)   # true negatives
  fp <- sum( predicted_conifer & !actual_conifer)   # false positives
  fn <- sum(!predicted_conifer &  actual_conifer)   # false negatives

  tibble(
    threshold   = t,
    sensitivity = tp / (tp + fn),
    specificity = tn / (tn + fp),
    accuracy    = (tp + tn) / nrow(ni_pixels),
    youden_j    = sensitivity + specificity - 1
  )

}) |>
  list_rbind()

# Best threshold = highest Youden's J
best_threshold <- threshold_performance |>
  slice_max(youden_j, n = 1)

best_threshold
threshold sensitivity specificity accuracy youden_j
0.495 0.9137931 0.9601594 0.9248418 0.8739525
Show code 
threshold_performance |>
  pivot_longer(c(sensitivity, specificity, youden_j),
               names_to = "metric", values_to = "value") |>
  mutate(metric = recode(metric,
    sensitivity = "Sensitivity",
    specificity = "Specificity",
    youden_j    = "Youden's J"
  )) |>
  ggplot(aes(x = threshold, y = value, colour = metric)) +
  geom_line(linewidth = 0.9) +
  geom_vline(xintercept = best_threshold$threshold,
             linetype = "dashed", colour = "black") +
  annotate("text",
    x     = best_threshold$threshold + 0.015,
    y     = 0.1,
    label = paste0("NI = ", best_threshold$threshold),
    hjust = 0, size = 3.5
  ) +
  scale_colour_manual(values = c(
    Sensitivity = "#1b7837",
    Specificity = "#762a83",
    "Youden's J" = "#e08214"
  )) +
  labs(x = "NI threshold", y = NULL, colour = NULL) +
  theme_minimal() +
  theme(legend.position = "top")
Figure 2: Sensitivity, specificity and Youden’s J across candidate NI thresholds. The vertical dashed line marks the selected threshold (highest Youden’s J).

The selected upper threshold is NI ≤ 0.495, achieving:

  • Sensitivity 91.4% — proportion of conifer training pixels correctly identified
  • Specificity 96.0% — proportion of non-conifer training pixels correctly excluded
  • Overall accuracy 92.5%

Lower threshold — conifer vs shadow / bare ground

Not everything below the upper threshold is conifer: very dark pixels (cast shadows, water, roads, bare rock) also have low NI. We use the 5th percentile of all conifer training pixels as the lower bound — pixels darker than 95 % of training conifer pixels are excluded from the conifer class.

Show code 
conifer_pixels <- ni_pixels |> filter(type == "Conifer")

# 5th percentile of the conifer NI distribution
lower_threshold <- quantile(conifer_pixels$NI, 0.05) |> round(3)

cat("Lower threshold (5th percentile of conifer training pixels):", lower_threshold, "\n")
Lower threshold (5th percentile of conifer training pixels): 0.258
Show code 
cat("Upper threshold (best Youden's J):", best_threshold$threshold, "\n")
Upper threshold (best Youden's J): 0.495
Show code 
cat("Conifer NI window: [", lower_threshold, ",", best_threshold$threshold, "]\n")
Conifer NI window: [ 0.258 , 0.495 ]
Show code 
ggplot(ni_pixels, aes(x = NI, fill = type)) +
  # Three-zone shading
  annotate("rect", xmin = -Inf,              xmax = lower_threshold,        ymin = 0, ymax = Inf,
           fill = "grey70",    alpha = 0.25) +
  annotate("rect", xmin = lower_threshold,   xmax = best_threshold$threshold, ymin = 0, ymax = Inf,
           fill = "#2d7a2d",   alpha = 0.12) +
  annotate("rect", xmin = best_threshold$threshold, xmax = Inf,             ymin = 0, ymax = Inf,
           fill = "#cc6600",   alpha = 0.12) +
  geom_density(alpha = 0.55, colour = "grey30", linewidth = 0.4) +
  # Threshold lines
  geom_vline(xintercept = lower_threshold,        linetype = "dashed", colour = "grey40",  linewidth = 0.8) +
  geom_vline(xintercept = best_threshold$threshold, linetype = "dashed", colour = "black", linewidth = 0.8) +
  # Zone labels
  annotate("text", x = lower_threshold / 2,                                    y = Inf,
           label = "Excluded\n(too dark)", vjust = 1.3, size = 3, colour = "grey40") +
  annotate("text", x = (lower_threshold + best_threshold$threshold) / 2,       y = Inf,
           label = "Conifer", vjust = 1.3, size = 3, colour = "#2d7a2d") +
  annotate("text", x = best_threshold$threshold + 0.08,                        y = Inf,
           label = "Broadleaf /\nopen", vjust = 1.3, size = 3, colour = "#cc6600") +
  scale_fill_manual(values = c(Conifer = "#2d7a2d", `Non-conifer` = "#cc6600")) +
  labs(x = "Needleleaf Index", y = "Density", fill = NULL) +
  theme_minimal() +
  theme(legend.position = "top")
Figure 3: NI distributions of conifer and non-conifer training pixels, with the two thresholds marking the conifer detection window. Pixels to the left of the lower threshold are excluded as too dark; pixels to the right of the upper threshold are classified as broadleaf / open.
Note

Both thresholds are derived from training data and represent performance within the training set. They will tend to be optimistic — independent validation with held-out samples would give a more realistic estimate of performance across the full image.

Per-species performance

Show code 
ni_pixels |>
  mutate(
    # Three-class prediction using both thresholds
    predicted = case_when(
      NI < lower_threshold               ~ "Excluded (too dark)",
      NI <= best_threshold$threshold     ~ "Conifer",
      TRUE                               ~ "Non-conifer"
    ),
    correct = case_when(
      type == "Conifer"     ~ predicted == "Conifer",
      type == "Non-conifer" ~ predicted == "Non-conifer"
    )
  ) |>
  group_by(species, type) |>
  summarise(
    n_pixels    = n(),
    pct_correct = round(mean(correct) * 100, 1),
    .groups     = "drop"
  ) |>
  rename(
    Species       = species,
    `True class`  = type,
    Pixels        = n_pixels,
    `% correct`   = pct_correct
  )
Table 2: Proportion of training pixels correctly classified using the two-threshold window. Conifer pixels falling outside [lower, upper] are missed; non-conifer pixels inside the window are false positives.
Species True class Pixels % correct
Lodgepole pine Conifer 102 13.7
Norway spruce Conifer 396 79.8
Sitka spruce Conifer 2788 92.1
Other conifers Conifer 86 88.4
Scots pine Conifer 865 92.0
Larch Conifer 577 67.2
Other broadleaves Non-conifer 102 84.3
Oak Non-conifer 480 99.6
Birch Non-conifer 579 93.8
Open Non-conifer 345 98.3

Larch has the lowest rate among conifers: its high summer NI pushes many pixels above the upper threshold. A small share of all species falls below the lower threshold (counted as “excluded”), mostly pixels on shaded edges of crowns.

Apply the thresholds to the full image

Using both thresholds, each pixel is assigned to one of three classes:

NI range Class
NI < 0.258 Excluded (too dark — shadow, water, bare ground)
0.258 ≤ NI ≤ 0.495 Conifer
NI > 0.495 Broadleaf / open
Show code 
# Conifer: NI within the calibrated window
conifer_raster <- ifel(
  ni >= lower_threshold & ni <= best_threshold$threshold,
  1, NA
)
names(conifer_raster) <- "conifer"

# Broadleaf / open: NI above the upper threshold
broadleaf_raster <- ifel(ni > best_threshold$threshold, 1, NA)
names(broadleaf_raster) <- "broadleaf"
Show code 
# Build a single three-class raster for display: 1 = conifer, 2 = broadleaf, NA = excluded
display_raster <- ifel(
  ni >= lower_threshold & ni <= best_threshold$threshold, 1,
  ifel(ni > best_threshold$threshold, 2, NA)
)

plot(
  display_raster,
  col    = c("#2d7a2d", "#cc6600"),
  type   = "classes",
  levels = c("Conifer", "Broadleaf / open"),
  main   = "NI classification — Aberfoyle"
)
Figure 4: Classification result. Green = conifer; orange = broadleaf / open; white = excluded (too dark to classify).

Area summary

Show code 
pixel_area_ha <- (res(ni)[1] * res(ni)[2]) / 10000   # 1 pixel = 100 m² = 0.01 ha

n_total      <- global(!is.na(ni),        fun = "sum", na.rm = TRUE)[[1]]
n_conifer    <- global(conifer_raster,    fun = "sum", na.rm = TRUE)[[1]]
n_broadleaf  <- global(broadleaf_raster, fun = "sum", na.rm = TRUE)[[1]]
n_excluded   <- n_total - n_conifer - n_broadleaf

tibble(
  Class              = c("Conifer", "Broadleaf / open", "Excluded (too dark)", "Total"),
  `Pixels`           = c(n_conifer, n_broadleaf, n_excluded, n_total),
  `Area (ha)`        = round(c(n_conifer, n_broadleaf, n_excluded, n_total) * pixel_area_ha),
  `% of study area`  = round(c(n_conifer, n_broadleaf, n_excluded, n_total) / n_total * 100, 1)
)
Class Pixels Area (ha) % of study area
Conifer 1023717 10237 24.3
Broadleaf / open 2850741 28507 67.6
Excluded (too dark) 345729 3457 8.2
Total 4220187 42202 100.0

Export outputs

Show code 
# Create output folder
if (!dir.exists("outputs/needleleaf")) {
  dir.create("outputs/needleleaf", recursive = TRUE)
}

# Helper: raster → cleaned polygons
# We aggregate to 20 m before polygonising to reduce single-pixel fragments
raster_to_polygons <- function(r, keep_column) {
  r |>
    aggregate(fact = 2, fun = "modal") |>   # merge neighbouring pixels
    as.polygons()                       |>  # raster → polygons
    st_as_sf()                          |>  # convert to sf
    filter(!is.na(.data[[keep_column]]))    # drop NA (non-classified pixels)
}

# 1. Continuous NI raster
writeRaster(
  ni,
  filename  = "outputs/needleleaf/needleleaf_index.tif",
  datatype  = "FLT4S",
  overwrite = TRUE
)

# 2. Conifer raster + shapefile
writeRaster(
  conifer_raster,
  filename  = "outputs/needleleaf/conifer.tif",
  datatype  = "INT1U",
  overwrite = TRUE
)

conifer_polygons <- raster_to_polygons(conifer_raster, "conifer")

st_write(
  conifer_polygons,
  dsn        = "outputs/needleleaf/conifer.shp",
  delete_dsn = TRUE,
  quiet      = TRUE
)

# 3. Broadleaf / open raster + shapefile
writeRaster(
  broadleaf_raster,
  filename  = "outputs/needleleaf/broadleaf.tif",
  datatype  = "INT1U",
  overwrite = TRUE
)

broadleaf_polygons <- raster_to_polygons(broadleaf_raster, "broadleaf")

st_write(
  broadleaf_polygons,
  dsn        = "outputs/needleleaf/broadleaf.shp",
  delete_dsn = TRUE,
  quiet      = TRUE
)

cat(
  "Outputs written to outputs/needleleaf/\n",
  "  needleleaf_index.tif   — continuous NI values\n",
  "  conifer.tif / .shp     — conifer pixels/polygons\n",
  "  broadleaf.tif / .shp   — broadleaf / open pixels/polygons\n"
)
Outputs written to outputs/needleleaf/
   needleleaf_index.tif   — continuous NI values
   conifer.tif / .shp     — conifer pixels/polygons
   broadleaf.tif / .shp   — broadleaf / open pixels/polygons

Caveats

  • Larch is taxonomically a conifer but drops its needles in winter. In the July image its canopy reflectance is more similar to broadleaves than to evergreen conifers, so the index will miss a portion of larch stands.
  • The lower threshold (5th percentile of conifer training pixels) excludes very dark pixels but is not independently validated. In areas with deep cast shadows or dark water bodies within the forest, some conifer edges may be excluded.
  • Both thresholds were optimised on the training polygons only. Performance may be lower in areas spectrally different from the sampled stands (e.g., young plantations, recently thinned areas, steep north-facing slopes with permanent shadow).
  • The Sentinel-2 B11 and B12 bands have a native resolution of 20 m, resampled to 10 m in this stack. The NI captures structural/spectral differences rather than individual crowns.
  • No independent validation has been performed. Cross-validation or an independent test set would give a more reliable accuracy estimate.

References

Amiri, A., Soltani, K., Gumiere, S. J., & Bonakdari, H. (2025). Forest fires under the lens: needleleaf index — a novel tool for satellite image analysis. npj Natural Hazards, 2, 9. https://doi.org/10.1038/s44304-025-00063-w