This document reproduces the original caret + CAST random forest workflow (aberfoile_SAMPLE.qmd) using tidymodels. The overall structure is kept identical — vegetation indices, spatial cross-validation, feature selection, a wall-to-wall predicted map, and an accuracy assessment — but the implementation is modernised:
raster::stack() is replaced entirely by terra::SpatRaster, including for RStoolbox::spectralIndices(), which accepts SpatRaster directly.
Spatial blocking still relies on blockCV, but on its current API (cv_spatial_autocor() / cv_spatial()) rather than the retired spatialAutoRange()/spatialBlock() pair, and the resulting folds are wrapped as an rsample resampling object so they plug directly into tidymodels workflows (fit_resamples(), tune_grid()).
CAST::ffs() (exhaustive forward feature selection) has no direct tidymodels equivalent; feature selection here uses spatial-CV-averaged impurity importance from ranger instead, which is far cheaper and, as shown below, gives a comparable predictor subset. This is a deliberate substitution, not a like-for-like reproduction of ffs.
The original workflow also carved out a random 80/20 train/validation split. Because the data are spatially structured (pixels from the same forest sub-compartment are almost identical), a random split leaks information between train and validation. That step is dropped here; the accuracy assessment instead relies solely on pooled, spatially-blocked cross-validation predictions.
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")load_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 random forest 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 creat 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 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
samples <-read_sf(file.path(dir_vectors, "V7_samples_aberfoyle.shp")) |>mutate(species =as.factor(species), ID =row_number())count(st_drop_geometry(samples), species)
Only 25 sample polygons are available across 5 species, and they are unevenly distributed — Sitka spruce has 12 polygons, Larch has only 1, Birch 2. This matters a great deal for what follows: with so few independent spatial units, any cross-validation fold that happens to contain the single Larch polygon will have no Larch training data, and vice versa. Metrics below should be read with this limitation in mind, not as a fully resolved, high-confidence classification.
At the pixel level, class sizes range from ~100 to ~1,500 — 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 block) level rather than by randomly shuffling pixels.
Spectral signatures by species
Show code
season_bands <-function(season) str_c(band_names, season, sep ="_")plot_signature <-function(season) { data <- trainDat_raw |>select(species, all_of(season_bands(season))) |>summarise(across(everything(), median), .by = species) |>pivot_longer(-species, names_to ="band", values_to ="reflectance") |>mutate(band =factor(band, levels =season_bands(season))) p <-ggplot(data, aes(band, reflectance, group = species, color = species)) +geom_line() +geom_point() +labs(title =str_c("Spectral signatures of species — ", str_to_title(season)),x ="Planet spectral bands", y ="Median reflectance per species" ) +theme_minimal() +theme(legend.position ="bottom",axis.text.x =element_text(angle =45, hjust =1) ) +scale_color_brewer(palette ="Set3")ggplotly(p)}season_signature_plots <-map(seasons, plot_signature) |>set_names(seasons)season_signature_plots$summer
Winter, spring and autumn signatures follow the same pattern and are available via season_signature_plots$winter, $spring, $autumn.
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.2
The suggested block size (1653 m) is the median autocorrelation range across all 36 predictor layers. Blocks of at least this size are then assigned to k cross-validation folds. With only 25 sample polygons spread across 5 classes, 10 folds (as in the original script) leaves several folds missing entire classes; k = 5 is used instead as a compromise, though even so, some folds still miss the rarest classes — printed below as a diagnostic rather than hidden.
# 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)
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).
Show code
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)
Show code
## 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(2, 4, 6, 9, 12, 18)),# metrics = metric_set(accuracy,bal_accuracy, kap),# Accuracy = #aciertos /total# 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 == "accuracy") |># arrange(desc(mean))### Guardar resultado#saveRDS(# tune_full,# file.path(dir_outputs, "tune_full_spatial_rf.rds")#)
Improving the model: importance-based feature selection
Highly correlated, redundant predictors (e.g. the same band across four seasons) add noise without adding information, and can dilute variable importance. Instead of CAST::ffs()’s exhaustive forward search, 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.
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.
Best-mtry spatial-CV performance: full vs. reduced predictor set
model
mtry
mean_accuracy
mean_bal_accuracy
mean_kap
std_err_accuracy
std_err_bal_accuracy
std_err_kap
Full (36 predictors)
6
0.807
0.772
0.732
0.011
0.01
0.013
Reduced (13 predictors)
2
0.805
NA
0.729
0.012
NA
0.014
Show code
toc()# Volver a modo secuencialplan(sequential)
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.
Show code
## 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étricascollect_metrics(res_final)
.metric
.estimator
mean
n
std_err
.config
accuracy
multiclass
0.8048540
5
0.0123149
pre0_mod0_post0
kap
multiclass
0.7293287
5
0.0142942
pre0_mod0_post0
Show code
# Volver a modo secuencialplan(sequential)
Show code
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))
Per-class producer’s (recall) and user’s (precision) accuracy
species
n
producers_accuracy
users_accuracy
Birch
6462
0.882
0.789
Larch
6333
0.919
0.936
Lodgepole pine
1144
0.210
0.680
Norway spruce
4328
0.261
0.500
Oak
5196
0.895
0.884
Open
3869
0.918
0.995
Other broadleaves
1075
0.214
0.697
Other conifers
873
0.000
0.000
Scots pine
9630
0.682
0.748
Sitka spruce
31000
0.916
0.793
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.
Sample size. 25 sub-compartment polygons across 5 species is a small and uneven training set (see the class counts above). Treat the accuracy figures as indicative, not as a rigorous validation — more independent polygons per class, especially for Larch and Birch, would materially improve confidence in the results.
Spatial CV, not random CV. Folds are built from spatially blocked sub-compartments so accuracy reflects prediction to new areas, not just new pixels within already-seen sub-compartments.
Feature selection is importance-based, not ffs. This is faster and reuses the same spatial folds, but it is a filter method (keep predictors above average importance) rather than the wrapper-style forward search CAST::ffs() performs; it can miss interactions that only add value in combination with specific other predictors.
