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 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 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.
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, "layer_stack.tif"),# overwrite = TRUE#)### Save the layer names separately (GeoTIFF does not preserve them)#saveRDS(# names(layer_stack),# file.path(dir_outputs, "layer_stack_names.rds")#)layer_stack <-rast(file.path(dir_outputs, "layer_stack.tif"))names(layer_stack) <-readRDS(file.path(dir_outputs, "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
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, "extracted.rds"))extracted <-readRDS(file.path(dir_outputs,"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.
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
Show code
season_signature_plots$autumn
Show code
season_signature_plots$winter
Show code
season_signature_plots$spring
Random Forest classification
Spatial cross-validation with kNNDM
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
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.
Show code
tic("time to build kNNDM folds: ")set.seed(4127)knndm_folds <-knndm(tpoints = sample_points, # labelled sample representativesmodeldomain = subcompartments, # area we predict todist_space ="geographical", # match distances in geographic spacek =5# number of cross-validation folds)# Persist the CV design so it can be reloaded without recomputingsaveRDS(knndm_folds, file =file.path(dir_outputs, "knndm_folds.rds"))toc()
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
Baseline model: all 36 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).
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()) # Removes 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
## Use all cores but one#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, kap),# control = control_grid(# save_pred = TRUE,# verbose = TRUE,# allow_par = TRUE# )#)##toc()### Back to sequential execution#plan(sequential)## View results ordered by kappa#collect_metrics(tune_full) |># filter(.metric == "kap") |># arrange(desc(mean))### Save result#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.
Show code
## Parallelisation#plan(# multisession,# workers = parallel::detectCores() - 1#)##cat("Workers available:", future::nbrOfWorkers(), "\n")##tic("time to importance-extraction")### Finalise the best workflow#wf_final_full <- finalize_workflow(# wf_full,# best_mtry_full#)##set.seed(4127)### Parallelised spatial evaluation#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# )#)### Save immediately, just in case#saveRDS(# res_full,# file.path(dir_outputs, "res_full_spatial_rf.rds")#)### Extract variable importance#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))### Save results#saveRDS(# importance_tbl,# file.path(dir_outputs, "importance_tbl.rds")#)### Show importance#print(importance_tbl)##toc()### Back to sequential mode#plan(sequential)
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 (14 predictors)
2
0.815
NA
0.745
0.017
NA
0.023
Show code
toc()
Time to extract data: : 73.324 sec elapsed
Show code
# Back to sequential modeplan(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
## Parallelisation#plan(# multisession,# workers = parallel::detectCores() - 1#)##cat("Workers available:", future::nbrOfWorkers(), "\n")##tic("time to final-model-evaluation")### Finalise the best workflow#wf_final <- finalize_workflow(# wf_reduced,# best_mtry_reduced#)##set.seed(4127)### Parallelised spatial evaluation#res_final <- fit_resamples(# wf_final,# resamples = spatial_folds,# metrics = metric_set(accuracy, kap),# control = control_resamples(# save_pred = TRUE,# allow_par = TRUE# )#)##toc()### Save results in case the session closes#saveRDS(# res_final,# file.path(dir_outputs, "res_final.rds")#)res_final <-readRDS(file.path(dir_outputs, "res_final.rds"))# View metricscollect_metrics(res_final)
.metric
.estimator
mean
n
std_err
.config
accuracy
multiclass
0.8152806
5
0.0169227
pre0_mod0_post0
kap
multiclass
0.7447034
5
0.0226680
pre0_mod0_post0
Show code
# Back to sequential modeplan(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))
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.865
0.734
Larch
6333
0.905
0.939
Lodgepole pine
1144
0.631
0.872
Norway spruce
4328
0.284
0.521
Oak
5196
0.905
0.903
Open
3869
0.875
0.997
Other broadleaves
1075
0.096
0.273
Other conifers
873
0.000
0.000
Scots pine
9630
0.731
0.800
Sitka spruce
31000
0.922
0.812
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 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. 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.
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.
