Aberfoyle Species Classification — AlphaEarth Embeddings

Random Forest classification with spatial cross-validation

This document reproduces the random forest workflow with a single predictor source: 64-band Google/DeepMind AlphaEarth annual embedding raster (year 2023). Each pixel carries a 64-dimensional learned representation of the full year of satellite observations.

Predictor stack: AlphaEarth embeddings

The AlphaEarth raster has 64 layers (A00–A63), each an abstract embedding dimension rather than a physical spectral band, at 10 m.

layer_stack <- rast(file.path(dir_embeddings, "alphaearth_aberfoyle_2023.tif"))

message("Number of layers: ", nlyr(layer_stack))


predictor_names <- names(layer_stack)

Embedding dimensions have no direct visual meaning, so a true-colour composite isn’t possible. As a purely illustrative check that the data varies spatially in a sensible way, the plot below maps three arbitrary embedding dimensions to red/green/blue.

plotRGB(
  layer_stack, r = 47, g = 1, b = 41,
  stretch = "lin", smooth = TRUE, axes = FALSE,
  main = "AlphaEarth embeddings (A46 / A00 / A40) — not true colour"
)

Field data: response variable

The response variable is tree species. The species included in the analysis were selected because they adequately represent the vegetation components present in the study area that are suitable for classification.

Training samples were selected from sub-compartments dominated by a single species. This approach ensures that the AlphaEarth extracted information is representative of only one species and helps prevent contamination of the training data caused by the presence of multiple species within the same sampling unit.

Several less abundant species, such as Douglas fir, Western hemlock…, were grouped into the Other Conifers class, while species such as beech…, alder were grouped into the Other Broadleaves class. As a result, these classes contain a mixture of species with potentially different spectral and structural characteristics, which increases within-class variability and can lead to lower classification accuracy compared with classes represented by a single species.

samples <- read_sf(file.path(dir_vectors, "V7_samples_aberfoyle.shp")) |>
  st_transform(crs(layer_stack)) |>
  mutate(species = as.factor(species), ID = row_number())

count(st_drop_geometry(samples), species)

species	n
Birch	50
Larch	72
Lodgepole pine	20
Norway spruce	38
Oak	50
Open	21
Other broadleaves	19
Other conifers	29
Scots pine	80
Sitka spruce	180

library(mapview)
mapview(
  samples, zcol = "species",
  col.regions = rainbow(n_distinct(samples$species)),
  layer.name = "Samples", map.types = "Esri.WorldImagery"
)

Extract predictor values at sample locations

#samples_vect <- vect(samples)
#
#extracted <- terra::extract(layer_stack, samples_vect, xy = TRUE, ID = TRUE) |>
#  as_tibble() |>
#  drop_na()
#
#saveRDS(extracted, file.path(dir_outputs, "extracted.rds"))

extracted <- readRDS(file.path(dir_outputs, "extracted.rds"))

message("Total number of 10 × 10 m sample locations extracted: ", nrow(extracted))

Data by species:

label_lookup <- samples |> st_drop_geometry() |> select(ID, species)

trainDat_raw <- extracted |> left_join(label_lookup, by = "ID")

count(trainDat_raw, species)

species	n
Birch	579
Larch	577
Lodgepole pine	102
Norway spruce	396
Oak	480
Open	345
Other broadleaves	102
Other conifers	86
Scots pine	865
Sitka spruce	2788

Recall that pixels within one polygon are highly spatially autocorrelated, which is exactly why the folds below are built at the polygon (spatial block) level rather than by randomly shuffling pixels.

Species distribution in AlphaEarth embedding space

Principal Component Analysis (PCA) is used to project the 64-dimensional embedding space onto a small number of orthogonal axes, making it easier to visualise how observations from different species are distributed within the AlphaEarth feature space.

emb_pca <- trainDat_raw |>
  select(all_of(predictor_names)) |>
  prcomp(scale. = TRUE)

scores <- emb_pca$x |>
  as_tibble() |>
  bind_cols(species = trainDat_raw$species)

ggplot(scores, aes(PC1, PC2, colour = species)) +
  geom_point(alpha = 0.6) +
  stat_ellipse(linewidth = 0.8) +
  theme_minimal() +
  labs(
    title = "Species separation in AlphaEarth embedding space"
  )

UMAP representation of the embedding space

UMAP (Uniform Manifold Approximation and Projection) is a non-linear dimensionality reduction technique that can reveal cluster structure that PCA — being constrained to linear projections — may miss. The same 64 embedding dimensions are reduced to two UMAP axes here, with species labels overlaid.

set.seed(123)

umap_coords <- uwot::umap(
  trainDat_raw |> dplyr::select(all_of(predictor_names)),
  n_neighbors = 15,
  min_dist    = 0.1,
  metric      = "cosine", #For learned embeddings (AlphaEarth, transformer embeddings, etc.), it's usually the direction of the vector that matters, not so much its magnitude. That is why it is usually recommended ‘cosine’, vs ‘euclidean’
  scale       = TRUE
)

umap_df <- tibble(
  UMAP1   = umap_coords[, 1],
  UMAP2   = umap_coords[, 2],
  species = trainDat_raw$species
)

ggplot(umap_df, aes(UMAP1, UMAP2, colour = species)) +
  geom_point(alpha = 0.5, size = 0.8) +
  stat_ellipse(linewidth = 0.8) +
  theme_minimal() +
  labs(
    title    = "Species separation in AlphaEarth embedding space (UMAP)",
    subtitle = "n_neighbours = 15, min_dist = 0.1, Cosine metric",
    colour   = "Species"
  ) +
  guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1)))

Random Forest classification

Spatial cross-validation with blockCV

Random shuffling of pixels for cross-validation would place near-duplicate, spatially autocorrelated pixels from the same sub-compartment in both the training and assessment folds, giving an overly optimistic estimate of accuracy. blockCV estimates the spatial autocorrelation range of the predictor stack, and folds are then built so that whole sub-compartments — not individual pixels — are held out together.¹

#set.seed(4127)
#sac <- cv_spatial_autocor(r = layer_stack, num_sample = nrow(extracted), plot = TRUE)
#range_m <- round(sac$range, 0)
#range_m
#
#saveRDS(sac, file = file.path(dir_outputs, "sac_object.rds"))

sac <- readRDS(file.path(dir_outputs, "sac_object.rds"))
sac[["plots"]][["barchart"]]

range_m <- round(sac$range, 0)

message("The suggested block size is ", range_m, " m")

The suggested block size (1273 m) is the median autocorrelation range across all 64 embedding layers. Blocks of at least this size are then assigned to k cross-validation folds. With 529 sample polygons spread across 10 classes k = 5 is used.

The output below shows the number of training and testing samples per species in each of the five spatial cross-validation folds.

set.seed(4127)
blocks <- cv_spatial(
  x = samples,
  column = "species",
  r = layer_stack,
  size = range_m,
  k =5,
  selection = "random",
  progress = FALSE,
  report = TRUE,
  plot = TRUE
)

  train_Birch train_Larch train_Lodgepole pine train_Norway spruce train_Oak
1          33          51                   16                  30        44
2          43          59                   17                  31        32
3          45          59                   16                  33        39
4          42          61                   19                  32        43
5          37          58                   12                  26        42
  train_Open train_Other broadleaves train_Other conifers train_Scots pine
1         15                      18                   22               60
2         17                      16                   20               67
3         15                      13                   23               69
4         19                      14                   23               64
5         18                      15                   28               60
  train_Sitka spruce test_Birch test_Larch test_Lodgepole pine
1                140         17         21                   4
2                153          7         13                   3
3                139          5         13                   4
4                142          8         11                   1
5                146         13         14                   8
  test_Norway spruce test_Oak test_Open test_Other broadleaves
1                  8        6         6                      1
2                  7       18         4                      3
3                  5       11         6                      6
4                  6        7         2                      5
5                 12        8         3                      4
  test_Other conifers test_Scots pine test_Sitka spruce
1                   7              20                40
2                   9              13                27
3                   6              11                41
4                   6              16                38
5                   1              20                34

The table below shows the distribution of samples by species across the five spatial cross-validation folds, indicating the number of observations assigned to the test set in each fold.

samples$spatial_fold <- factor(blocks$folds_ids)
count(st_drop_geometry(samples), spatial_fold, species) |>
  pivot_wider(names_from = species, values_from = n, values_fill = 0)

spatial_fold	Birch	Larch	Lodgepole pine	Norway spruce	Oak	Open	Other broadleaves	Other conifers	Scots pine	Sitka spruce
1	17	21	4	8	6	6	1	7	20	40
2	7	13	3	7	18	4	3	9	13	27
3	5	13	4	5	11	6	6	6	11	41
4	8	11	1	6	7	2	5	6	16	38
5	13	14	8	12	8	3	4	1	20	34

Build the predictor table and fold assignment

The predictor dataset was created by joining the extracted pixel values with the corresponding species labels and spatial fold assignments. The resulting table contains the predictor variables, the target species, and the spatial fold identifier for each sample. A summary of sample counts per species and fold was then generated to assess the distribution of observations across the spatial cross-validation folds.

label_lookup <- samples |> st_drop_geometry() |> select(ID, species, spatial_fold)

trainDat <- extracted |>
  left_join(label_lookup, by = "ID") |>
  select(-x, -y)

predictor_names <- setdiff(names(trainDat), c("ID", "species", "spatial_fold"))
count(trainDat, spatial_fold, species) |>
  pivot_wider(names_from = species, values_from = n, values_fill = 0)

spatial_fold	Birch	Larch	Lodgepole pine	Norway spruce	Oak	Open	Other broadleaves	Other conifers	Scots pine	Sitka spruce
1	182	156	10	69	52	67	4	19	207	619
2	98	102	45	95	179	84	31	22	169	422
3	57	134	24	51	133	152	22	15	122	644
4	79	95	2	61	64	14	31	24	160	590
5	163	90	21	120	52	28	14	6	207	513

# Summary of the minimum number of samples per class and fold
count(trainDat, spatial_fold, species) |>
  group_by(species) |>
  summarise(
    min_fold = min(n),
    max_fold = max(n),
    mean_fold = mean(n)
  ) |>
  arrange(min_fold)

species	min_fold	max_fold	mean_fold
Lodgepole pine	2	45	20.4
Other broadleaves	4	31	20.4
Other conifers	6	24	17.2
Open	14	152	69.0
Norway spruce	51	120	79.2
Oak	52	179	96.0
Birch	57	182	115.8
Larch	90	156	115.4
Scots pine	122	207	173.0
Sitka spruce	422	644	557.6

make_fold <- function(f) {
  make_splits(
    list(
      analysis   = which(trainDat$spatial_fold != f),
      assessment = which(trainDat$spatial_fold == f)
    ),
    data = trainDat
  )
}

fold_ids <- sort(unique(trainDat$spatial_fold))
spatial_folds <- manual_rset(
  map(fold_ids, make_fold),
  ids = str_c("Fold", fold_ids)
)
spatial_folds

Baseline model: all 64 predictors

A random forest (ranger engine) is tuned over mtry using the spatial folds above, optimising Cohen’s kappa (more informative than accuracy here, given the class imbalance).

rec_full <- recipe(species ~ ., data = trainDat) |>
  update_role(ID, spatial_fold, new_role = "id") |> # ID and Spatial_fold should not be used to train the model. It marks them as identification variables.
  step_zv(all_predictors()) #Eliminates variables with zero variance.

rf_spec <- rand_forest(mtry = tune(), trees = 300) |>
  set_engine("ranger", importance = "impurity", num.threads = 1) |>
  set_mode("classification")

wf_full <- workflow() |> add_recipe(rec_full) |> add_model(rf_spec)

## Usar todos los cores menos uno
#workers <- parallel::detectCores() - 1
#
#plan(multisession, workers = workers)
#
#tic("time to tune")
#
#set.seed(4127)
#
#tune_full <- tune_grid(
#  wf_full,
#  resamples = spatial_folds,
#  grid = tibble(mtry = c(4, 8, 12, 16, 24, 32, 48)),
#  metrics = metric_set(accuracy, kap),
#  control = control_grid(
#    save_pred = TRUE,
#    verbose = TRUE,
#    allow_par = TRUE
#  )
#)
#
#toc()
#
## Volver a ejecución secuencial
#plan(sequential)
#
## Ver resultados ordenados por Kappa
#collect_metrics(tune_full) |>
#  filter(.metric == "kap") |>
#  arrange(desc(mean))
#
## Guardar resultado
#saveRDS(
#  tune_full,
#  file.path(dir_outputs, "tune_full_spatial_rf.rds")
#)

tune_full <- readRDS(
  file.path(dir_outputs, "tune_full_spatial_rf.rds")
)

best_mtry_full <- select_best(tune_full, metric = "kap")
best_mtry_full

mtry	.config
8	pre0_mod2_post0

Improving the model: importance-based feature selection

Some embedding dimensions are likely more informative than others for separating these species, and dropping the weakest ones can reduce noise. Impurity importance is extracted from the ranger engine fitted in every spatial fold and averaged, giving an importance ranking that is itself cross-validated rather than computed on a single fit.

## Paralelización
#plan(
#  multisession,
#  workers = parallel::detectCores() - 1
#)
#
#cat("Workers disponibles:", future::nbrOfWorkers(), "\n")
#
#tic("time to importance-extraction")
#
## Fijar mejor workflow
#wf_final_full <- finalize_workflow(
#  wf_full,
#  best_mtry_full
#)
#
#set.seed(4127)
#
## Evaluación espacial paralelizada
#res_full <- fit_resamples(
#  wf_final_full,
#  resamples = spatial_folds,
#  metrics = metric_set(accuracy, kap),
#  control = control_resamples(
#    save_pred = TRUE,
#    allow_par = TRUE,
#    extract = \(x) extract_fit_engine(x)$variable.importance
#  )
#)
#
## Guardar inmediatamente por si acaso
#saveRDS(
#  res_full,
#  file.path(dir_outputs, "res_full_spatial_rf.rds")
#)
#
## Extraer importancia de variables
#importance_tbl <- res_full |>
#  select(id, .extracts) |>
#  unnest(.extracts) |>
#  unnest_longer(
#    .extracts,
#    values_to = "importance",
#    indices_to = "predictor"
#  ) |>
#  summarise(
#    mean_importance = mean(importance),
#    .by = predictor
#  ) |>
#  arrange(desc(mean_importance))
#
## Guardar resultados
#saveRDS(
#  importance_tbl,
#  file.path(dir_outputs, "importance_tbl.rds")
#)
#
## Mostrar importancia
#print(importance_tbl)
#
#toc()
#
## Volver a modo secuencial
#plan(sequential)

res_full <- readRDS(
  file.path(dir_outputs, "res_full_spatial_rf.rds")
)

importance_tbl <- readRDS(
  file.path(dir_outputs, "importance_tbl.rds")
)

importance_tbl |>
  slice_max(mean_importance, n = 20) |>
  ggplot(aes(mean_importance, fct_reorder(predictor, mean_importance))) +
  geom_col() +
  labs(
    title = "Mean impurity importance across 5 spatial CV folds",
    x = "Importance", y = NULL
  )

Predictors with above-average importance are kept as the reduced set — a simple, transparent filter rather than a tuned threshold, given how few independent spatial units are available to tune one reliably.

importance_threshold <- mean(importance_tbl$mean_importance)

selected_predictors <- importance_tbl |>
  filter(mean_importance > importance_threshold) |>
  pull(predictor)

length(selected_predictors)

[1] 21

selected_predictors

 [1] "A46" "A00" "A40" "A11" "A22" "A58" "A39" "A25" "A51" "A57" "A21" "A05"
[13] "A63" "A45" "A52" "A62" "A18" "A50" "A07" "A24" "A34"

Refit and compare against the full-predictor model

## Paralelización
#plan(
#  multisession,
#  workers = parallel::detectCores() - 1
#)
#
#cat("Workers disponibles:", future::nbrOfWorkers(), "\n")
#
#tic("time to reduced-model-tuning")
#
## Recipe reducida
#rec_reduced <- recipe(
#  as.formula(
#    str_c(
#      "species ~ ",
#      str_c(selected_predictors, collapse = " + ")
#    )
#  ),
#  data = trainDat
#) |>
#  step_zv(all_predictors())
#
## Workflow reducido
#wf_reduced <- workflow() |>
#  add_recipe(rec_reduced) |>
#  add_model(rf_spec)
#
#set.seed(4127)
#
## Tuning paralelizado
#tune_reduced <- tune_grid(
#  wf_reduced,
#  resamples = spatial_folds,
#  grid = tibble(mtry = c(2, 4, 6, 9, 13, 18)),
#  metrics = metric_set(accuracy, kap),
#  control = control_grid(
#    save_pred = TRUE,
#    allow_par = TRUE,
#    verbose = TRUE
#  )
#)
#
## Guardar inmediatamente
#saveRDS(
#  tune_reduced,
#  file.path(dir_outputs, "tune_reduced_spatial_rf.rds")
#)

tune_reduced <- readRDS(
  file.path(dir_outputs, "tune_reduced_spatial_rf.rds")
)

# Mejor mtry
best_mtry_reduced <- select_best(
  tune_reduced,
  metric = "kap"
)

## Comparación Full vs Reduced
#comparison_tbl <- bind_rows(
#  collect_metrics(tune_full) |>
#    filter(mtry == best_mtry_full$mtry) |>
#    mutate(
#      model = str_c(
#        "Full (",
#        length(predictor_names),
#        " predictors)"
#      )
#    ),
#
#  collect_metrics(tune_reduced) |>
#    filter(mtry == best_mtry_reduced$mtry) |>
#    mutate(
#      model = str_c(
#        "Reduced (",
#        length(selected_predictors),
#        " predictors)"
#      )
#    )
#) |>
#  select(
#    model,
#    mtry,
#    .metric,
#    mean,
#    std_err
#  ) |>
#  pivot_wider(
#    names_from = .metric,
#    values_from = c(mean, std_err)
#  )
#
## Guardar tabla de comparación
#saveRDS(
#  comparison_tbl,
#  file.path(dir_outputs, "comparison_tbl.rds")
#)

comparison_tbl <- readRDS(
  file.path(dir_outputs, "comparison_tbl.rds")
)

# Mostrar tabla
knitr::kable(
  comparison_tbl,
  digits = 3,
  caption = "Best-mtry spatial-CV performance: full vs. reduced predictor set"
)

Best-mtry spatial-CV performance: full vs. reduced predictor set

model	mtry	mean_accuracy	mean_kap	std_err_accuracy	std_err_kap
Full (64 predictors)	8	0.828	0.764	0.020	0.028
Reduced (21 predictors)	2	0.822	0.756	0.024	0.031

The reduced predictor set performs comparably to (and, on this run, no worse than) the full set at its best mtry, so it is carried forward — fewer, less redundant predictors also make the final map faster to produce and somewhat easier to reason about.

Accuracy assessment

Because a random hold-out would leak spatially autocorrelated pixels between train and test, accuracy here is assessed from pooled, out-of-fold predictions of the 5-fold spatial cross-validation on the reduced model — every prediction comes from a fold that never saw that sub-compartment during training.

## Paralelización
#plan(
#  multisession,
#  workers = parallel::detectCores() - 1
#)
#
#cat("Workers disponibles:", future::nbrOfWorkers(), "\n")
#
#tic("time to final-model-evaluation")
#
## Fijar el mejor workflow
#wf_final <- finalize_workflow(
#  wf_reduced,
#  best_mtry_reduced
#)
#
#set.seed(4127)
#
## Evaluación espacial paralelizada
#res_final <- fit_resamples(
#  wf_final,
#  resamples = spatial_folds,
#  metrics = metric_set(accuracy, kap),
#  control = control_resamples(
#    save_pred = TRUE,
#    allow_par = TRUE
#  )
#)
#
#toc()
#
## Guardar resultados por si cierras la sesión
#saveRDS(
#  res_final,
#  file.path(dir_outputs, "res_final.rds")
#)

res_final <- readRDS(
  file.path(dir_outputs, "res_final.rds")
)

# Ver métricas
collect_metrics(res_final)

.metric	.estimator	mean	n	std_err	.config
accuracy	multiclass	0.821625	5	0.0242588	pre0_mod0_post0
kap	multiclass	0.755808	5	0.0309185	pre0_mod0_post0

oof_preds <- collect_predictions(res_final)

conf_mat(oof_preds, truth = species, estimate = .pred_class) |>
  autoplot(type = "heatmap") +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(title = "Confusion matrix — pooled 5-fold spatial CV predictions") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

species_levels <- levels(oof_preds$species)

per_class_accuracy <- map_dfr(species_levels, \(cl) {
  truth_bin <- factor(if_else(oof_preds$species == cl, cl, "other"), levels = c(cl, "other"))
  pred_bin  <- factor(if_else(oof_preds$.pred_class == cl, cl, "other"), levels = c(cl, "other"))
  tibble(
    species = cl,
    n = sum(oof_preds$species == cl),
    producers_accuracy = sens_vec(truth_bin, pred_bin, event_level = "first"),
    users_accuracy = precision_vec(truth_bin, pred_bin, event_level = "first")
  )
})

per_class_accuracy |>
  knitr::kable(digits = 3, caption = "Per-class producer's (recall) and user's (precision) accuracy")

Per-class producer’s (recall) and user’s (precision) accuracy

species	n	producers_accuracy	users_accuracy
Birch	579	0.796	0.799
Larch	577	0.860	0.832
Lodgepole pine	102	0.206	0.600
Norway spruce	396	0.402	0.619
Oak	480	0.871	0.848
Open	345	0.951	0.968
Other broadleaves	102	0.000	0.000
Other conifers	86	0.000	NA
Scots pine	865	0.773	0.771
Sitka spruce	2788	0.951	0.847

Producer’s/user’s accuracy for Larch and Birch should be treated cautiously: with only 1 and 2 sample polygons respectively, their out-of-fold predictions come from at most one or two folds each, so these numbers are based on very little independent evidence rather than a stable estimate.

Spatial model predictions

The final workflow is refit on all training pixels using the reduced predictor set and the selected mtry, then applied to the full raster stack and clipped to the forest sub-compartments.

#set.seed(4127)
#final_fit <- fit(wf_final, data = trainDat)
#
#saveRDS(
#  final_fit,
#  file.path(dir_outputs, "final_fit.rds")
#)

final_fit <- readRDS(
  file.path(dir_outputs, "final_fit.rds")
)

#prediction <- terra::predict(
#  layer_stack[[selected_predictors]],
#  final_fit,
#  fun = \(model, ...) predict(model, ...)$.pred_class,
#  na.rm = TRUE
#)
#
#subcompartments <- st_read(file.path(dir_vectors, "subcompartmentCLIP.shp"), quiet = TRUE) |>
#  st_transform(crs(prediction))
#
#prediction_clip <- prediction |>
#  crop(subcompartments) |>
#  mask(subcompartments)
#
#writeRaster(prediction_clip, file.path(dir_outputs, "RF_AberfoyleCLIP_AlphaEarth.tif"), overwrite = TRUE)

prediction_clip <- rast(
  file.path(dir_outputs, "RF_AberfoyleCLIP_AlphaEarth.tif")
)

classes <- levels(prediction_clip)[[1]]$class

species_colours <- setNames(
  hcl.colors(length(classes), "Dark 3"),
  classes
)

rasterVis::levelplot(
  raster::raster(prediction_clip),
  maxpixel = 1e6,
  col.regions = species_colours[classes],
  scales = list(draw = FALSE),
  main = "Aberfoyle forest — species classification (AlphaEarth embeddings)"
)

Final model: variable importance

rf_engine <- extract_fit_engine(final_fit)

tibble(
  predictor  = names(rf_engine$variable.importance),
  importance = rf_engine$variable.importance
) |>
  slice_max(importance, n = 20) |>
  ggplot(aes(importance, fct_reorder(predictor, importance))) +
  geom_col() +
  labs(
    title = "Final model: predictor importance (Gini impurity)",
    subtitle = str_c(
      "Random forest, ", rf_engine$num.trees, " trees, ",
      length(rf_engine$variable.importance), " predictors, mtry = ", rf_engine$mtry
    ),
    x = "Importance", y = NULL
  )

Footnotes

1. https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.13107