Optional hands-on exercises: Spain examples

Where your models can be trusted: evaluating spatial machine learning reliably

The exercises use the Spain data in data.zip1. The data/spain folder contains a predictor stack, a study-area boundary, a temperature sample, and a PM\(_{2.5}\) sample.

Choose one path for the exercises:

  1. Temperature: a sample with a more random spatial pattern
  2. PM\(_{2.5}\): a more clustered sample

The aim is to practice three ideas from the workshop:

  1. How training sample design affects representativeness
  2. How validation strategy affects estimated performance
  3. How AoA and LPD help interpret where predictions may be more or less trustworthy

1 Setup

The code below attaches the necessary libraries, loads the data, and defines a function to compute RMSE. You may modify it to fit your needs, e.g., by adding more packages or defining a different performance metric.

library(CAST)
library(sf)
library(terra)
library(dplyr)
library(ggplot2)
library(caret)
library(blockCV)

set.seed(44)

if (!dir.exists("data")) {
  data_url <- "https://jakubnowosad.com/ml4eo2026workshop/data.zip"
  download.file(data_url, "data.zip")
  unzip("data.zip", exdir = "data")
}

spain <- read_sf("data/spain/spain.gpkg")
predictors <- rast("data/spain/predictors.tif")
predictor_names <- names(predictors)

samples_temp <- read_sf("data/spain/temp_train.gpkg")
samples_pm25 <- read_sf("data/spain/PM25_train.gpkg")

rmse_fun <- function(obs, pred) {
  sqrt(mean((obs - pred) ^ 2, na.rm = TRUE))
}

2 Exercise 1: Inspect the chosen data

Using your chosen response variable, answer these questions:

  • What are the six predictor variables in the raster stack?
  • How many observations are in the chosen sample?
  • What is the response variable?
  • Does the sample look more random or more clustered in geographic space?

Create one figure that shows the chosen sample on the Spain boundary.

Solution (PM2.5 path) 
predictor_names
nrow(samples_pm25)
# "PM25"
# "more clustered"

plot(st_geometry(spain))
plot(samples_pm25["PM25"], add = TRUE) # check the spatial pattern

3 Exercise 2: Prepare the modeling table

Prepare the predictor table for the chosen sample.

Tasks:

  • If you chose temperature, extract the six predictor values from the raster stack
  • If you chose PM\(_{2.5}\), keep only the same six predictors used in the raster stack
  • Check whether any predictor values are missing
Solution (PM2.5 path) 
pm25_train <- samples_pm25 |>
  st_drop_geometry() |>
  dplyr::select(PM25, all_of(predictor_names))

colSums(is.na(pm25_train))

4 Exercise 3: Compute AoA and LPD for the chosen sample

Compute the Area of Applicability for the chosen training sample using the six predictors.

Tasks:

  • Plot the DI summary output
  • Create an AoA map
  • Create an LPD map
  • Compute the proportion of the study area classified as inside the AoA
  • Try to interpret these results
Solution (PM2.5 path) 
lpd_data <- CAST::aoa(
  newdata = predictors,
  train = pm25_train[, predictor_names],
  variables = predictor_names,
  LPD = TRUE,
  verbose = FALSE
)

plot(lpd_data) +
  ggtitle("DI summary for the PM2.5 sample")
plot(lpd_data$AOA, main = "AoA from the PM2.5 sample")
plot(lpd_data$LPD, main = "LPD from the PM2.5 sample")

# aoa coverage
terra::global(lpd_data$AOA, "mean", na.rm = TRUE)

5 Exercise 4: Evaluate the model with different validation strategies

Fit a model of your choice (e.g., random forest) for the chosen response variable, and evaluate it with different validation strategies.

Compare these validation strategies:

  • random 5-fold CV
  • spatial 5-fold CV with blockCV::cv_spatial()
  • 5-fold kNNDM CV with CAST::knndm()

Compute RMSE for each strategy and compare the results. Interpret the differences.

Solution (PM2.5 path) 
tune_grid <- expand.grid(
  mtry = 2,
  splitrule = "variance",
  min.node.size = 5
)

ctrl_random <- trainControl(
  method = "cv",
  number = 5,
  savePredictions = "final",
  verboseIter = FALSE
)

rf_random <- train(
  x = pm25_train[, predictor_names],
  y = pm25_train$PM25,
  method = "ranger",
  trControl = ctrl_random,
  tuneGrid = tune_grid,
  metric = "RMSE",
  num.trees = 300,
  importance = "impurity"
)

random_rmse <- rmse_fun(rf_random$pred$obs, rf_random$pred$pred)

spatial_folds <- blockCV::cv_spatial(
  x = samples_pm25,
  k = 5,
  seed = 44
)

train_ids <- lapply(spatial_folds$folds_list, function(x) x[[1]])
test_ids <- lapply(spatial_folds$folds_list, function(x) x[[2]])

ctrl_spatial <- trainControl(
  method = "cv",
  index = train_ids,
  indexOut = test_ids,
  savePredictions = "final",
  verboseIter = FALSE
)

rf_spatial <- train(
  x = pm25_train[, predictor_names],
  y = pm25_train$PM25,
  method = "ranger",
  trControl = ctrl_spatial,
  tuneGrid = tune_grid,
  metric = "RMSE",
  num.trees = 300,
  importance = "impurity"
)

spatial_rmse <- rmse_fun(rf_spatial$pred$obs, rf_spatial$pred$pred)

kn <- CAST::knndm(
  tpoints = samples_pm25,
  modeldomain = predictors[[1]],
  k = 5
)

ctrl_knndm <- trainControl(
  method = "cv",
  index = kn$indx_train,
  indexOut = kn$indx_test,
  savePredictions = "final",
  verboseIter = FALSE
)

rf_knndm <- train(
  x = pm25_train[, predictor_names],
  y = pm25_train$PM25,
  method = "ranger",
  trControl = ctrl_knndm,
  tuneGrid = tune_grid,
  metric = "RMSE",
  num.trees = 300,
  importance = "impurity"
)

knndm_rmse <- rmse_fun(rf_knndm$pred$obs, rf_knndm$pred$pred)

data.frame(
  method = c("Random CV", "Spatial CV", "kNNDM CV"),
  rmse = c(random_rmse, spatial_rmse, knndm_rmse)
)

6 Exercise 5: Map AoA, LPD, and expected error from the kNNDM-trained model

Use the model fitted with kNNDM CV to compute:

  • the AoA map
  • the LPD map
  • an error profile based on LPD
  • a map of expected prediction error
  • a map of predicted values masked by the AoA
Solution (PM2.5 path) 
lpd_model_knndm <- CAST::aoa(
  newdata = predictors,
  model = rf_knndm,
  LPD = TRUE,
  verbose = FALSE
)

plot(lpd_model_knndm$AOA, main = "AoA from the kNNDM-trained PM2.5 model")
plot(lpd_model_knndm$LPD, main = "LPD from the kNNDM-trained PM2.5 model")

errormodel_lpd <- errorProfiles(rf_knndm, lpd_model_knndm, variable = "LPD")
plot(errormodel_lpd) +
  ggtitle("Error profile based on LPD")

expected_error_lpd <- terra::predict(lpd_model_knndm$LPD, errormodel_lpd)
plot(expected_error_lpd, main = "Expected prediction error for the PM2.5 model")

predictions <- predict(predictors, rf_knndm, na.rm = TRUE)
predictions_aoa_masked <- mask(predictions, lpd_model_knndm$AOA, maskvalue = 0)
plot(predictions_aoa_masked, main = "Predicted PM2.5 masked by the AoA")

Footnotes

  1. The data is explained and used in Milà et al. (2022)↩︎