Multi-Layer Perceptron classification with kNNDM spatial cross-validation
Author
Forest Research
Published
17 September 2026
Show code
#| include: false# Spatial datalibrary(terra)library(sf)# Modellinglibrary(tidymodels)# Spatial cross-validation by nearest-neighbour distance matchinglibrary(CAST)# MLP engine (torch backend)library(brulee)# Remote sensinglibrary(RStoolbox)# Data wrangling / visualisationlibrary(tidyverse)library(mapview)library(uwot)library(rasterVis)library(plotly)#Access performinglibrary(tictoc)library(future)tidymodels_prefer()set.seed(4127)# --- Project paths ---------------------------------------------------------dir_raster <-"Planet"dir_vectors <-"Shapefiles"dir_vi <-"VI"# Artifacts shared with the companion random-forest workflow (raw predictor# stack, pixel extraction, kNNDM folds) are re-used as-is — none of them# depend on which classifier is fitted downstream.dir_outputs_rf <-"OutputsRF"# New artifacts specific to this MLP workflow are kept separate so the two# analyses never overwrite each other's cached objects.dir_outputs <-"OutputsMLP_PlanetVI"dir.create(dir_outputs, showWarnings =FALSE)
This document reproduces the Planet-imagery species-classification workflow, but with a Multi-Layer Perceptron (MLP) in place of the random forest used in the companion analysis. An MLP requires feature scaling and benefits from explicit regularisation (dropout + early stopping), which shapes the recipe and tuning choices below; the predictor generation, sampling and spatial cross-validation design are otherwise unchanged from the random forest workflow.
Generate predictor variables from the Planet images
Planet images
Having images from different seasons provides a greater diversity of data, which can help train a more robust model. Bands are loaded as SpatRaster objects and renamed per season.
Show code
band_names <-c("coastal_blue", "blue", "green_I", "green","yellow", "red", "rededge", "nir")season_dirs <-c(winter ="10 March 2023_psscene_analytic_8b_sr_udm2",spring ="30 May 23_psscene_analytic_8b_sr_udm2",summer ="6sep24_ConvexHull_psscene_analytic_8b_sr_udm2",autumn ="11 November 2023_psscene_analytic_8b_sr_udm2")# Read one season's 8-band composite and give each band a season-tagged nameload_season <-function(season) { r <-rast(file.path(dir_raster, season_dirs[[season]], "composite.tif"))names(r) <-str_c(band_names, season, sep ="_") r}seasons <-names(season_dirs)season_rasters <-map(seasons, load_season) |>set_names(seasons)
The following is a false-colour (red / NIR / blue) composite of Aberfoyle Forest for each season.
Show code
par(mfrow =c(2, 2), mar =c(2, 2, 2, 2))season_labels <-c(winter ="Winter - 10th March 2023",spring ="Spring - 30th May 2023",summer ="Summer - 06th September 2024",autumn ="Autumn - 11th November 2023")iwalk(season_rasters, \(r, season) {plotRGB(r, r =6, g =8, b =2, stretch ="lin",main = season_labels[[season]], smooth =TRUE, axes =FALSE)})
Create vegetation indices (Summer)
Vegetation indices from the summer scene can improve the model by providing features that capture the peak of vegetative growth, when foliage is densest and most uniform — this allows for a clearer distinction between species. Vegetation indices also help reduce noise from raw spectral data. See the index reference for details.1
Since the closed canopy limits soil influence, indices sensitive to chlorophyll content, canopy structure, and vegetative vigour are prioritised over indices designed to correct for soil or water effects:
MCARI is highly sensitive to variations in chlorophyll and canopy structure, useful for distinguishing species with subtle differences in greenness.
GNDVI uses the green band and is more responsive to chlorophyll levels than NDVI, helping differentiate species with similar biomass but differing health.
NDVI remains a reliable general-purpose baseline, complementing MCARI and GNDVI.
MSAVI2 adds structural contrast, valuable when species vary in crown density or canopy cover.
Show code
#tic("time to create indexes: ")#VI <- spectralIndices(# season_rasters$summer,# redEdge1 = "rededge_summer",# green = "green_summer",# red = "red_summer",# nir = "nir_summer",# blue = "blue_summer",# indices = c("NDVI", "GNDVI", "MCARI", "MSAVI2"),# scaleFactor = 1,# skipRefCheck = TRUE#)##writeRaster(VI, file.path(dir_vi, "Indices_summer.tif"), overwrite = TRUE)#toc()VI <-rast(file.path(dir_vi, "Indices_summer.tif"))names(VI)
The predictor stack combines the original Planet bands for each season with the summer vegetation indices. This is identical to, and reused from, the random-forest workflow — the predictor stack does not depend on the downstream classifier.
Show code
#layer_stack <- c(# season_rasters$winter,# season_rasters$spring,# season_rasters$summer,# season_rasters$autumn,# VI#)### Save the predictor stack to disk#writeRaster(# layer_stack,# file.path(dir_outputs_rf, "layer_stack.tif"),# overwrite = TRUE#)### Save the layer names separately (GeoTIFF does not preserve them)#saveRDS(# names(layer_stack),# file.path(dir_outputs_rf, "layer_stack_names.rds")#)layer_stack <-rast(file.path(dir_outputs_rf, "layer_stack.tif"))names(layer_stack) <-readRDS(file.path(dir_outputs_rf, "layer_stack_names.rds"))nlyr(layer_stack)
[1] 36
Field data: response variable
The response variable is tree species from the Forester sub-compartment database, restricted to sub-compartments covered by a single species so that spectral signatures are as clean as possible.
Show code
# One polygon per sub-compartment; ID is the join key back to extracted pixelssamples <-read_sf(file.path(dir_vectors, "V7_samples_aberfoyle.shp")) |>mutate(species =as.factor(species), ID =row_number())count(st_drop_geometry(samples), species)
The sample polygons are unevenly distributed across species. This matters a great deal for what follows: with few independent spatial units per class, any cross-validation fold that happens to contain all polygons of a rare species will leave the model with no training data for it, and vice versa. Metrics below should be read with this limitation in mind, not as a fully resolved, high-confidence classification.
Extract predictor values at sample locations
Extraction does not depend on the classifier, so the cached pixel table from the random-forest workflow is reused directly.
Show code
tic("Time to extract data: ")samples_vect <-vect(samples)## Pull every predictor value for every pixel falling inside a sample polygon.## xy = TRUE keeps pixel coordinates; ID links each pixel to its polygon.#extracted <- terra::extract(layer_stack, samples_vect, xy = TRUE, ID = TRUE) |># as_tibble() |># drop_na()#toc()#saveRDS(extracted, file.path(dir_outputs_rf, "extracted.rds"))extracted <-readRDS(file.path(dir_outputs_rf,"extracted.rds"))nrow(extracted)
At the pixel level, class sizes are much healthier than the polygon counts above — but recall that pixels within one polygon are highly spatially autocorrelated, which is exactly why the folds below are built at the polygon (spatial unit) level rather than by randomly shuffling pixels.
Species distribution in the Planet + VI predictor space (PCA)
Principal Component Analysis is used to project the 36-dimensional predictor space (32 seasonal spectral bands + 4 summer vegetation indices) onto a small number of orthogonal axes, making it easier to visualise how observations from different species are distributed within that predictor 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 36 predictors are reduced to two UMAP axes here, with species labels overlaid. Unlike the AlphaEarth embedding workflow (where a cosine metric is preferred because only the direction of the learned embedding vector is meaningful), these predictors are physically-scaled reflectance values and vegetation indices, so the ordinary Euclidean metric is used instead.
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. kNNDM (k-fold nearest-neighbour distance matching) builds the folds so that the geometry of the validation problem resembles the geometry of the real prediction problem: it arranges the folds so that the distances from held-out points to their nearest training point follow the same distribution as the distances from the prediction locations to their nearest sample point.2 The same spatial cross-validation scheme is used as in the companion random-forest workflow — the classifier changes below, not the resampling.
Prediction domain and sample representatives
kNNDM needs two sets of locations:
the training points — the labelled samples. We use one representative point per sub-compartment polygon (its centroid). These polygons are the genuinely independent spatial units: pixels within a polygon are near-duplicates, so folding at the polygon level keeps a whole sub-compartment together and reflects the true amount of independent data.
the prediction domain — the area the final map is produced for. The map is clipped to the forest sub-compartments, so those polygons define the domain, and kNNDM samples prediction points from within it.
Show code
# Prediction domain: the sub-compartments the final map is clipped tosubcompartments <-read_sf(file.path(dir_vectors, "subcompartmentCLIP.shp")) |>st_transform(st_crs(samples))# One representative point per independent spatial unit (sub-compartment)sample_points <-st_centroid(samples)
Generate the folds
knndm() clusters the sample points and searches for the fold assignment whose held-out-to-training distance distribution best matches the prediction distance distribution. k = 5 is used: with the class imbalance present here, more folds would leave several folds missing entire classes. The fold design is identical to the random-forest workflow — it depends only on the sample geometry and prediction domain, not on the classifier — so the cached folds are reused directly.
Show code
# kNNDM is cheap but cached to match the rest of the workflow. Uncomment to# regenerate the fold design.#tic("time to build kNNDM folds: ")#set.seed(4127)#knndm_folds <- knndm(# tpoints = sample_points, # labelled sample representatives# modeldomain = subcompartments, # area we predict to# dist_space = "geographical", # match distances in geographic space# k = 5 # number of cross-validation folds#)##saveRDS(knndm_folds, file = file.path(dir_outputs_rf, "knndm_folds_Planet.rds"))#toc()knndm_folds <-readRDS(file.path(dir_outputs_rf, "knndm_folds_Planet.rds"))knndm_folds
knndm object
Space:
Clustering algorithm: hierarchical
Intermediate clusters (q): 346
W statistic: 85.1676
Number of folds: 5
Observations in each fold: 122 112 99 114 112
Inspect the quality of the distance matching
kNNDM reports a W statistic and produces three empirical cumulative distribution functions (ECDFs) that let us judge whether the folds represent the prediction problem well.
Show code
# The plot method returns a ggplot comparing the three distance ECDFs.plot(knndm_folds)
The three curves are:
Ĝij (prediction) — distances from prediction points to their nearest sample point. This is the target distribution: the distances the model must actually predict at.
**Ĝ*j (kNNDM CV) — distances from held-out points to their nearest training point under the kNNDM folds. We want this curve to lie on top of** the prediction curve.
Ĝj (LOO) — the leave-one-out reference, shown as the optimistic extreme (held-out points sit very close to a training point).
Show code
# W is the Wasserstein distance (area between the ECDFs) separating the kNNDM CV# curve from the prediction curve, in the units of the distance axis (metres).# Smaller W means the CV folds reproduce prediction conditions more faithfully.knndm_folds$W
[1] 85.16762
Assign folds to samples
Each sample point corresponds to one polygon, so knndm_folds$clusters is that polygon’s fold. Because sample_points was built in the same row order as samples, we attach the fold directly. The column is named spatial_fold so the rest of the workflow can consume it unchanged.
Show code
samples$spatial_fold <-factor(knndm_folds$clusters)# Class representation per fold — confirm no fold is missing a class entirelycount(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
13
14
7
9
7
3
2
12
22
33
2
9
13
2
14
15
1
8
4
8
38
3
8
7
6
4
8
7
1
4
17
37
4
10
18
2
7
10
4
6
6
18
33
5
10
20
3
4
10
6
2
3
15
39
Show code
# Spatial view of the fold assignment across the forestmapview( samples, zcol ="spatial_fold",col.regions =rainbow(n_distinct(samples$spatial_fold)),layer.name ="kNNDM fold", map.types ="Esri.WorldImagery")
Build the predictor table and fold assignment
Each pixel inherits the spatial fold of the polygon it was extracted from, so a whole sub-compartment is always held out together.
# Summary of the minimum number of samples per class and foldcount(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
40
531
228.8
Other conifers
56
367
174.6
Open
84
1940
773.8
Other broadleaves
84
483
215.0
Larch
446
2069
1266.6
Norway spruce
530
1543
865.6
Oak
635
1956
1039.2
Scots pine
647
2562
1926.0
Birch
861
1866
1292.4
Sitka spruce
5691
6693
6200.0
Show code
# Turn the fold labels into an rsample object: for each fold, everything else is# the analysis (training) set and the fold itself is the assessment set.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. 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 singletune() 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. The single-hidden-layer version is therefore the natively supported tuning path taken here.
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 36-dimensional Planet + VI 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.
Show code
mlp_rec <-recipe(species ~ ., data = trainDat) |>update_role(ID, spatial_fold, new_role ="id") |>step_zv(all_predictors()) |># drop zero-variance predictorsstep_normalize(all_numeric_predictors()) # MLPs need scaled inputsmlp_spec <-mlp(hidden_units =tune(), # width of the single hidden layerpenalty =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 layerepochs =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 stoppingstop_iter =15# patience (epochs without improvement) ) |>set_mode("classification")mlp_wflow <-workflow() |>add_recipe(mlp_rec) |>add_model(mlp_spec)mlp_wflow
# 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.## These ranges are deliberately more conservative than a first, unconstrained# guess would suggest. The training table here has ~70,000 pixel rows (one# row per 3 m pixel across all sample polygons) — an order of magnitude more# rows than a typical tabular tuning example — and a scratch timing test on# this machine showed each fit costs roughly 1 second per training epoch on# the CPU-only torch backend. An initial attempt at a wider grid (hidden# units up to 128, epochs up to 300) triggered out-of-memory / allocation# errors partway through fitting on this machine. The ranges below were# chosen to keep every candidate architecture inexpensive enough to fit# reliably across all 5 spatial folds without exhausting memory, while still# spanning small-to-moderate network capacity.mlp_params <-extract_parameter_set_dials(mlp_wflow) |>update(hidden_units =hidden_units(range =c(8L, 128L)),dropout =dropout(range =c(0, 0.6)),epochs =epochs(range =c(50L, 250L)) )mlp_params
name
id
source
component
component_id
object
hidden_units
hidden_units
model_spec
mlp
main
integer , 8 , 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 , 250 , TRUE , TRUE , # Epochs
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, Lodgepole pine) 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.
Reflecting the compute-budget reasoning above, the initial design and the number of Bayesian iterations are both kept small (6 space-filling points and up to 6 further Bayesian iterations, with early stopping if 6 iterations pass without improvement) — enough to explore the reduced hyperparameter space without the memory pressure seen with a larger grid.
Bayesian search did not improve on the best point from the initial space-filling design within the (small) iteration budget used here — this is plausible given how few distinct architectures a budget this size can explore, and is reported plainly rather than treated as evidence that the network hyperparameters are unimportant.
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. 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.
# Accuracy, kappa, macro F1 and weighted F1 on the pooled out-of-fold predictionsfinal_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.7605493
kap
multiclass
0.6591993
f_meas
macro
0.6900935
f_meas_weighted
macro_weighted
0.7846361
Show code
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))
Show code
species_levels <-levels(oof_preds$species)# Producer's accuracy = recall; user's accuracy = precision, computed one class# at a time in a one-vs-rest fashion.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
6462
0.785
0.684
Larch
6333
0.867
0.888
Lodgepole pine
1144
0.000
NA
Norway spruce
4328
0.000
NA
Oak
5196
0.802
0.857
Open
3869
0.924
0.991
Other broadleaves
1075
0.006
0.128
Other conifers
873
0.000
NA
Scots pine
9630
0.519
0.691
Sitka spruce
31000
0.963
0.736
Producer’s/user’s accuracy for the rarest species should be treated cautiously: with few sample polygons, their out-of-fold predictions come from very few folds, so these numbers are based on 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 36-predictor Planet + VI set — the dropout and internal early-stopping validation split, rather than a manual importance-based feature filter, are what control overfitting here — then applied to the full raster stack and clipped to the forest sub-compartments.
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.
Why the tuning ranges and search budget are conservative: the ~70,000-row pixel table and CPU-only torch backend made wider hidden-unit/epoch ranges and a larger grid computationally expensive and, at the widest settings tried, prone to out-of-memory errors on this machine. Narrower ranges (8–48 hidden units, 15–40 epochs) and a small initial design (6 points) plus limited Bayesian refinement (up to 6 further iterations) were chosen to keep every candidate fit affordable across all 5 spatial folds.
Why no importance-based feature reduction step: impurity importance is specific to tree ensembles (ranger) and is not used here. For the MLP, dropout, the fixed weight decay of zero, and the internal early-stopping validation split serve the equivalent purpose of controlling overfitting on the 36-dimensional Planet + VI predictor 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. Lodgepole pine/Other conifers with few polygons), accuracy alone would favour hyperparameters that just predict the majority class well; macro F1 weights all species equally regardless of prevalence.
Notes and caveats
Sample size. The training set is small and uneven across species (see the class counts above). Treat the accuracy figures as indicative, not as a rigorous validation — more independent polygons per class would materially improve confidence in the results.
Spatial CV, not random CV. Folds are built by kNNDM at the sub-compartment level and distance-matched to the prediction domain, so accuracy reflects prediction to new areas, not just new pixels within already-seen sub-compartments. The W statistic and ECDF plot above document how closely the CV represents those prediction conditions.
