Aberfoyle Species Classification — AlphaEarth Embeddings

Multi-Layer Perceptron (MLP) classification with spatial cross-validation

This document reproduces the species-classification 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. Unlike the companion random-forest analysis, the classifier here is a Multi-Layer Perceptron (MLP), which requires feature scaling and benefits from explicit regularisation (weight decay + dropout) given the high predictor dimensionality relative to the number of independent spatial samples.

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

#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_rf, "extracted.rds"))

extracted <- readRDS(file.path(dir_outputs_rf, "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	581
Larch	576
Lodgepole pine	103
Norway spruce	398
Oak	477
Open	343
Other broadleaves	95
Other conifers	88
Scots pine	868
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)))

Multi-Layer Perceptron 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.¹ This spatial-block setup is identical to the random-forest workflow — the classifier changes below, not the resampling scheme.

#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_rf, "sac_object.rds"))

sac <- readRDS(file.path(dir_outputs_rf, "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.

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	183	154	11	66	54	65	2	18	215	622
2	99	99	46	95	175	84	29	26	165	421
3	57	136	25	52	129	152	23	16	121	630
4	80	98	1	63	66	14	29	23	162	596
5	162	89	20	122	53	28	12	5	205	519

# 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	1	46	20.6
Other broadleaves	2	29	19.0
Other conifers	5	26	17.6
Open	14	152	68.6
Norway spruce	52	122	79.6
Oak	53	175	95.4
Birch	57	183	116.2
Larch	89	154	115.2
Scots pine	121	215	173.6
Sitka spruce	421	630	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

Recipe and model specification

Unlike random forest, an MLP is sensitive to the scale of its inputs, so a step_normalize() is added after the zero-variance filter. ID and spatial_fold keep the "id" role so they are carried through the recipe without being treated as predictors.

The architecture is a shallow, single-hidden-layer network fit with the brulee engine (torch backend), using mlp()’s native hidden_units argument as an ordinary tuning parameter (16-128 units). brulee does support deeper networks by passing a vector to hidden_units (e.g. c(64, 32) for two layers), but that argument only accepts a single tune() placeholder — wrapping two separate tune() calls inside c() is not supported by extract_parameter_set_dials() and raises Result must be length 1, not 2. A working two-hidden-layer alternative (a manual grid of fixed architectures, since tune_bayes() cannot search over vector-valued parameters) is shown at the end of this section for reference.

Regularisation combines two mechanisms, both tuned: dropout between layers and an internal validation split with early stopping (validation + stop_iter), which is important given how few independent spatial units (polygons) are available relative to the 64-dimensional input. L2 weight decay (penalty) is not tuned alongside dropout: the brulee engine raises "Both weight decay and dropout should not be specified" if both are non-zero at once, so penalty is fixed at 0 and dropout is used as the sole regulariser.

mlp_rec <- recipe(species ~ ., data = trainDat) |>
  update_role(ID, spatial_fold, new_role = "id") |>
  step_zv(all_predictors()) |>            # drop zero-variance predictors
  step_normalize(all_numeric_predictors()) # MLPs need scaled inputs

mlp_spec <- mlp(
  hidden_units = tune(),  # width of the single hidden layer
  penalty      = 0,       # L2 weight decay disabled: brulee errors if both
                           # penalty and dropout are non-zero, so dropout is
                           # used as the sole regulariser here. Set here (not
                           # via set_engine()) because penalty is a main
                           # parsnip arg, and brulee's own default is 0.001
                           # (not 0), which would silently reintroduce the
                           # same conflict once dropout is tuned.
  dropout      = tune(),  # dropout rate applied after the hidden layer
  epochs       = tune()   # max training epochs (early stopping via stop_iter)
) |>
  set_engine(
    "brulee",
    activation   = "relu",
    learn_rate   = 0.01,
    validation   = 0.15,  # internal holdout used for early stopping
    stop_iter    = 15     # patience (epochs without improvement)
  ) |>
  set_mode("classification")

mlp_wflow <- workflow() |>
  add_recipe(mlp_rec) |>
  add_model(mlp_spec)

mlp_wflow

══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: mlp()

── Preprocessor ────────────────────────────────────────────────────────────────
2 Recipe Steps

• step_zv()
• step_normalize()

── Model ───────────────────────────────────────────────────────────────────────
Single Layer Neural Network Model Specification (classification)

Main Arguments:
  hidden_units = tune()
  penalty = 0
  dropout = tune()
  epochs = tune()

Engine-Specific Arguments:
  activation = relu
  learn_rate = 0.01
  validation = 0.15
  stop_iter = 15

Computational engine: brulee 

# Set sensible tuning ranges for hidden-layer width, dropout, and epochs.
# penalty is fixed at 0 in mlp_spec (not tuned) since brulee disallows
# non-zero penalty and dropout simultaneously.
mlp_params <- extract_parameter_set_dials(mlp_wflow) |>
  update(
    hidden_units = hidden_units(range = c(16L, 128L)),
    dropout       = dropout(range = c(0, 0.6)),
    epochs        = epochs(range = c(50L, 300L))
  )

mlp_params

name	id	source	component	component_id	object
hidden_units	hidden_units	model_spec	mlp	main	integer , 16 , 128 , TRUE , TRUE , # Hidden Units
dropout	dropout	model_spec	mlp	main	double , 0 , 0.6 , TRUE , FALSE , Dropout Rate
epochs	epochs	model_spec	mlp	main	integer , 50 , 300 , TRUE , TRUE , # Epochs

NoteOptional: tuning a two-hidden-layer network instead

tune_bayes()/tune_grid() need a fixed-length numeric or categorical parameter per tunable argument, so a vector-valued hidden_units (multiple layers) can’t be searched the same way as a scalar. The practical workaround is a manual grid over a small, fixed set of candidate architectures, combined with tune_grid() (Bayesian optimisation doesn’t apply to this kind of discrete, structured parameter):

mlp_spec_2layer <- mlp(
  penalty = tune(),
  dropout = tune(),
  epochs  = tune()
) |>
  set_engine(
    "brulee",
    hidden_units = tune("architecture"),  # populated from a list-column below
    activation   = "relu",
    learn_rate   = 0.01,
    validation   = 0.15,
    stop_iter    = 15
  ) |>
  set_mode("classification")

mlp_wflow_2layer <- workflow() |>
  add_recipe(mlp_rec) |>
  add_model(mlp_spec_2layer)

# Candidate architectures: 1 or 2 hidden layers, various widths
architectures <- list(c(32), c(64), c(128), c(64, 32), c(128, 64))

manual_grid <- tidyr::crossing(
  architecture = architectures,
  penalty = 10^seq(-5, -1, length.out = 3),
  dropout = c(0, 0.3, 0.6),
  epochs  = c(100, 300)
)

# tune_grid(mlp_wflow_2layer, resamples = spatial_folds, grid = manual_grid, ...)

Hyperparameter tuning

Tuning is optimised on macro F1 (f_meas, which yardstick defaults to macro-averaging for a multiclass outcome) so that rare species (e.g. Larch, Birch) count as much as abundant ones like Sitka spruce, rather than being swamped by overall accuracy. Accuracy and Cohen’s kappa are also tracked for context. A space-filling initial design seeds a Bayesian optimisation (tune_bayes()) over the 5 spatial folds defined above.

## Parallel back-end
#workers <- max(1, parallel::detectCores() - 1)
#plan(multisession, workers = workers)
#
#tic("time to tune MLP (Bayesian optimisation)")
#
#set.seed(4127)
#initial_grid <- grid_space_filling(mlp_params, size = 20)
#
## Diagnostic: confirm the space-filling design actually has 20 distinct
## rows before spending time fitting it. tune_bayes() needs 2+ (ideally
## more than the number of tuning parameters) usable initial points to
## fit its Gaussian Process surrogate.
#nrow(initial_grid)        # should be 20
#initial_grid               # eyeball for duplicate rows
#
## tune_bayes()'s `initial` arg needs either a positive integer or the
## results of tune_grid() -- a raw grid_space_filling() data frame won't
## work on its own, so we run it through tune_grid() first to get metrics
## attached, then hand those results to tune_bayes() as the initial design.
#initial_res <- tune_grid(
#  mlp_wflow,
#  resamples  = spatial_folds,
#  param_info = mlp_params,
#  grid       = initial_grid,
#  metrics    = metric_set(f_meas, accuracy, kap),
#  control    = control_grid(save_pred = TRUE)   # must match control_bayes(save_pred = TRUE) below
#)
#
## Diagnostic: if tune_bayes() later complains about too few initial grid
## points, it's usually because most candidates in initial_grid errored
## out during fitting rather than initial_grid itself being too small.
## Check for fitting errors/warnings here first.
#collect_notes(initial_res)
#
#mlp_bayes_res <- tune_bayes(
#  mlp_wflow,
#  resamples  = spatial_folds,
#  param_info = mlp_params,
#  initial    = initial_res,
#  iter       = 25,
#  metrics    = metric_set(f_meas, accuracy, kap),
#  control = control_bayes(
#    save_pred   = TRUE,
#    verbose     = TRUE,
#    no_improve  = 10,
#    allow_par   = TRUE,
#    seed        = 4127
#  )
#)
#
#toc()
#plan(sequential)
#
#saveRDS(mlp_bayes_res, file.path(dir_outputs, "mlp_bayes_res.rds"))

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

# Top candidates by macro F1
show_best(mlp_bayes_res, metric = "f_meas", n = 5)

hidden_units	dropout	epochs	.metric	.estimator	mean	n	std_err	.config	.iter
69	0.3157895	181	f_meas	macro	0.7720068	5	0.0208828	pre0_mod10_post0	0
70	0.3060181	170	f_meas	macro	0.7567647	5	0.0323661	iter08	8
68	0.3294730	177	f_meas	macro	0.7555632	5	0.0253399	iter06	6
67	0.3043592	254	f_meas	macro	0.7548116	5	0.0253239	iter10	10
39	0.2842105	286	f_meas	macro	0.7546084	5	0.0164510	pre0_mod05_post0	0

autoplot(mlp_bayes_res, type = "performance", metric = "f_meas") +
  labs(title = "Bayesian optimisation progress — macro F1")

Best hyperparameters

best_mlp <- select_best(mlp_bayes_res, metric = "f_meas")
best_mlp

hidden_units	dropout	epochs	.config
69	0.3157895	181	pre0_mod10_post0

Accuracy assessment

As with the random-forest workflow, accuracy is assessed from pooled, out-of-fold predictions rather than a random hold-out, since a random split would leak spatially autocorrelated pixels between train and test. Because tune_bayes() was run with save_pred = TRUE, the out-of-fold predictions for the best hyperparameter combination can be pulled directly with collect_predictions() — no separate refit over the folds is needed.

oof_preds <- collect_predictions(mlp_bayes_res, parameters = best_mlp)

oof_preds

# Accuracy, kappa, macro F1 and weighted F1 on the pooled out-of-fold predictions
final_metrics <- bind_rows(
  accuracy(oof_preds, truth = species, estimate = .pred_class),
  kap(oof_preds,      truth = species, estimate = .pred_class),
  f_meas(oof_preds,   truth = species, estimate = .pred_class, estimator = "macro"),
  f_meas(oof_preds,   truth = species, estimate = .pred_class, estimator = "macro_weighted") |>
    mutate(.metric = "f_meas_weighted")
)

final_metrics

.metric	.estimator	.estimate
accuracy	multiclass	0.8054456
kap	multiclass	0.7366375
f_meas	macro	0.7789529
f_meas_weighted	macro_weighted	0.8204997

conf_mat_tbl <- conf_mat(oof_preds, truth = species, estimate = .pred_class)

# Row-normalised (i.e. by true class) heatmap, built from the same
# contingency table that autoplot(type = "heatmap") uses internally.
conf_mat_tbl$table |>
  as_tibble() |>
  group_by(Truth) |>
  mutate(prop = n / sum(n)) |>
  ungroup() |>
  ggplot(aes(Prediction, Truth, fill = prop)) +
  geom_tile() +
  geom_text(aes(label = scales::percent(prop, accuracy = 1)), size = 3) +
  scale_fill_gradient(low = "white", high = "steelblue", limits = c(0, 1), name = "Row %") +
  labs(
    title = "Normalised confusion matrix — pooled 5-fold spatial CV predictions (MLP)",
    x = "Predicted species", y = "True species"
  ) +
  theme_minimal() +
  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	581	0.797	0.696
Larch	576	0.800	0.857
Lodgepole pine	103	0.000	NA
Norway spruce	398	0.405	0.603
Oak	477	0.839	0.762
Open	343	0.927	0.964
Other broadleaves	95	0.000	NA
Other conifers	88	0.000	NA
Scots pine	868	0.812	0.724
Sitka spruce	2788	0.925	0.855

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 finalised with the best hyperparameters and refit on all training pixels using the full 64-dimensional predictor set (the strong dropout + weight-decay regularisation, rather than a manual importance-based feature filter, is what controls overfitting here), then applied to the full raster stack and clipped to the forest sub-compartments.

#final_wflow <- finalize_workflow(mlp_wflow, best_mlp)
#
#set.seed(4127)
#final_fit <- fit(final_wflow, data = trainDat)
#
#saveRDS(final_fit, file.path(dir_outputs, "final_fit_mlp.rds"))

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

#prediction <- terra::predict(
#  layer_stack[[predictor_names]],
#  final_fit,
#  fun = function(model, data) {
#
#    data$ID <- 0
#    data$spatial_fold <- NA
#
#    predict(model, new_data = data, type = "class")$.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, "MLP_AberfoyleCLIP_AlphaEarth.tif"), overwrite = TRUE)

prediction_clip <- rast(
  file.path(dir_outputs, "MLP_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 (MLP, AlphaEarth embeddings)"
)

Notes on the modelling choices

• Why a single tuned hidden_units rather than two: extract_parameter_set_dials() expects exactly one tune() placeholder per argument, so a vector like hidden_units = c(tune(), tune()) cannot be extracted as two separate parameters (tune_grid()/tune_bayes() require a fixed-length scalar or categorical value per parameter). The single-layer version tuned above is the natively supported path; the callout in the recipe/spec section shows a manual-grid workaround for searching fixed two-layer architectures.
• Why no importance-based feature reduction step: impurity importance is specific to tree ensembles (ranger). For the MLP, dropout, weight decay and the internal early-stopping validation split serve the equivalent purpose of controlling overfitting on the 64-dimensional embedding space, so all predictors are retained. If a reduced predictor set is preferred, vip::vi_permute() with a custom prediction wrapper is the appropriate permutation-importance route for a brulee model.
• Why macro F1 for tuning: with pronounced class imbalance (e.g. Sitka spruce vs. Larch/Birch with only 1-2 polygons), accuracy alone would favour hyperparameters that just predict the majority class well; macro F1 weights all species equally regardless of prevalence. ```

Footnotes

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