Introduction

This document shows the application of caret for spatial modelling at the example of predicting air temperature in Spain. Hereby, we use measurements of air temperature available only at specific locations in Spain to create a spatially continuous map of air temperature. Therefore, machine-learning models are trained to learn the relationship between spatially continuous predictors and air temperature.

When using machine-learning methods with spatial data, we need to take care of, e.g., spatial autocorrelation, as well as extrapolation when predicting to regions that are far away from the training data. To deal with these issues, several methods have been developed. In this document, we will show how to combine the machine-learning workflow of caret with packages designed to deal with machine-learning with spatial data. Hereby, we use blockCV::cv_spatial() and CAST::knndm() for spatial cross-validation, and CAST::aoa() to mask areas of extrapolation. We use sf and terra for processing vector and raster data, respectively.

The caret package

The caret package contains functions to train machine-learning models, as well as for, e.g., model selection. Its main function is caret::train(), which provides a uniform interface to over 200 machine-learning algorithms. (User-specified-) Cross-Validation methods can be defined via caret::trainControl(). An extensive online tutorial is available at https://topepo.github.io/caret/. Furthermore, a paper (https://doi.org/10.18637/jss.v028.i05), as well as a book (http://appliedpredictivemodeling.com/), describing the use of caret are available.

Install and load required R-packages

# install.packages("caret")
# install.packages("CAST")
# install.packages("blockCV")
# install.packages("sf")
# install.packages("terra")

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

Case study data

Load data needed in this modelling example:

• predictor stack: raster dataset of the environmental predictors used to predict air temperature
• train points: vector dataset of ground measurements of air temperature
• spain: region for which predictions are made

predictor_stack <- terra::rast(
    "https://github.com/LOEK-RS/FOSSGIS2025-examples/raw/refs/heads/main/data/predictors.tif"
)
predictor_names <- names(predictor_stack)
spain <- sf::st_read(
    "https://github.com/LOEK-RS/FOSSGIS2025-examples/raw/refs/heads/main/data/spain.gpkg"
)

Reading layer `spain' from data source
  `https://github.com/LOEK-RS/FOSSGIS2025-examples/raw/refs/heads/main/data/spain.gpkg'
  using driver `GPKG'
Simple feature collection with 1 feature and 0 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -13454.15 ymin: 3988025 xmax: 1020771 ymax: 4859816
Projected CRS: ED50 / UTM zone 30N

train_points <- sf::st_read(
    "https://github.com/LOEK-RS/FOSSGIS2025-examples/raw/refs/heads/main/data/temp_train.gpkg"
)

Reading layer `temp_train' from data source
  `https://github.com/LOEK-RS/FOSSGIS2025-examples/raw/refs/heads/main/data/temp_train.gpkg'
  using driver `GPKG'
Simple feature collection with 195 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 36026.79 ymin: 3988818 xmax: 978160.6 ymax: 4858999
Projected CRS: ED50 / UTM zone 30N

train_data <- terra::extract(
    predictor_stack,
    train_points,
    bind = TRUE,
    ID = FALSE
) |>
    sf::st_as_sf()

plot(sf::st_geometry(spain))
plot(sf::st_geometry(train_points), col = "blue4", add = TRUE)

Standard modelling workflow

Firstly, a simple modelling workflow without feature selection and hyperparameter tuning is shown:

1. split data into training and test data
2. train model using the training data only
3. predict on test data to obtain error metrics
4. predict on predictor stack to obtain spatially continuous prediction of air temperature

# 1. train-test split
set.seed(321)
trainIndex <- caret::createDataPartition(
    train_data$temp,
    p = .8,
    list = FALSE,
    times = 1
)
temperature_train <- train_data[trainIndex, ]
temperature_test <- train_data[-trainIndex, ]

# 2. model training
model <- caret::train(
    temp ~ .,
    data = sf::st_drop_geometry(temperature_train),
    method = "ranger",
    tuneGrid = expand.grid(
        "mtry" = 4,
        "splitrule" = "variance",
        "min.node.size" = 5
    ),
    num.trees = 100
)

# 3. predict on test data
test_df <- temperature_test[, "temp", drop = FALSE]
test_df$prediction <- predict(model, temperature_test)

test_metrics <- caret::postResample(
    pred = test_df$prediction,
    obs = test_df$temp
) |>
    round(3) |>
    t() |>
    as.data.frame()

print(test_metrics)

   RMSE Rsquared   MAE
1 0.913    0.901 0.741

# 4. predict to raster stack
prediction_spatial <- terra::predict(predictor_stack, model, na.rm = TRUE)

Spatial Cross-Validation for model selection

Cross-validation (CV) methods are often employed to obtain optimal hyperparameter values. Therefore, the training data are split into folds, and a model is trained on -1 folds. The fold not used for model training is then used to obtain the test statistic. This is repeated over all folds, and the test metric is averaged over the folds.

In spatial machine-learning, a spatial CV is often needed to prevent very similar data to be in the training and testing fold at the same time, which is often the case if training data are clustered and leads to overly optimistic CV error estimates. R packages that implement spatial CV include, e.g., blockCV and CAST. Here, we will explore the integration of those two with caret.

set.seed(333)
spatial_blocks <- blockCV::cv_spatial(
    temperature_train,
    k = 5,
    progress = FALSE
)

  train test
1   126   33
2   126   33
3   128   31
4   128   31
5   128   31

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

tr_control_block <- caret::trainControl(
    method = "cv",
    index = train_ids,
    indexOut = test_ids,
    savePredictions = TRUE
)

hyperparameter_grid <- expand.grid(
    "mtry" = c(2, 4, 6, 10, 12),
    "min.node.size" = c(5, 10, 15),
    "splitrule" = "variance"
)

model_spatial <- caret::train(
    temp ~ .,
    data = sf::st_drop_geometry(temperature_train),
    method = "ranger",
    trControl = tr_control_block,
    tuneGrid = hyperparameter_grid,
    num.trees = 100
)

model_spatial$finalModel

Ranger result

Call:
 ranger::ranger(dependent.variable.name = ".outcome", data = x,      mtry = min(param$mtry, ncol(x)), min.node.size = param$min.node.size,      splitrule = as.character(param$splitrule), write.forest = TRUE,      probability = classProbs, ...)

Type:                             Regression
Number of trees:                  100
Sample size:                      159
Number of independent variables:  22
Mtry:                             10
Target node size:                 5
Variable importance mode:         none
Splitrule:                        variance
OOB prediction error (MSE):       0.8557959
R squared (OOB):                  0.8952759 

knndm_folds <- CAST::knndm(
    tpoints = temperature_train,
    modeldomain = spain,
    space = "geographical",
    clustering = "kmeans",
    k = 5
)

tr_control_knndm <- caret::trainControl(
    method = "cv",
    index = knndm_folds$indx_train,
    indexOut = knndm_folds$indx_test,
    savePredictions = TRUE
)

hyperparameter_grid <- expand.grid(
    "mtry" = c(2, 4, 6, 10, 12),
    "min.node.size" = c(5, 10, 15),
    "splitrule" = "variance"
)

model_knndm <- caret::train(
    temp ~ .,
    data = sf::st_drop_geometry(temperature_train),
    method = "ranger",
    trControl = tr_control_knndm,
    tuneGrid = hyperparameter_grid,
    num.trees = 100
)

model_knndm$finalModel

Ranger result

Call:
 ranger::ranger(dependent.variable.name = ".outcome", data = x,      mtry = min(param$mtry, ncol(x)), min.node.size = param$min.node.size,      splitrule = as.character(param$splitrule), write.forest = TRUE,      probability = classProbs, ...)

Type:                             Regression
Number of trees:                  100
Sample size:                      159
Number of independent variables:  22
Mtry:                             12
Target node size:                 10
Variable importance mode:         none
Splitrule:                        variance
OOB prediction error (MSE):       0.8728343
R squared (OOB):                  0.8931909 

selected_hyperparams <- model_knndm$bestTune

model_ffs <- CAST::ffs(
    predictors = sf::st_drop_geometry(temperature_train)[, predictor_names],
    response = sf::st_drop_geometry(temperature_train)$temp,
    method = "ranger",
    num.trees = 100,
    trControl = tr_control_knndm,
    tuneGrid = selected_hyperparams,
    verbose = FALSE
)

plot(model_ffs, plotType = "selected")

# obtain prediction
prediction_ffs <- terra::predict(predictor_stack, model_ffs, na.rm = TRUE)

AOA

Lastly, the area which is too dissimilar from the training data for the models to make reliable predictions (area of applicability, AOA) is delineated using the function CAST::aoa(). The function CAST::aoa() takes as inputs the predictor stack, as well as the trained caret model. The resulting object contains the dissimilarity values, the threshold used to delineate the AOA (every dissimilarity value above this threshold is considered outside the AOA), as well as the final AOA raster. Since our training data are randomly distributed in the study area, most of the area falls within the AOA.

AOA_without_tuning <- CAST::aoa(
    newdata = predictor_stack,
    model = model,
    verbose = FALSE
)
AOA_with_tuning <- CAST::aoa(
    newdata = predictor_stack,
    model = model_ffs,
    verbose = FALSE
)

Compare predictions obtained by the un-tuned model vs the tuned model

par(mfrow = c(2, 2))
plot(prediction_spatial, main = "prediction without tuning")
plot(AOA_without_tuning$AOA, main = "AOA without tuning")
plot(prediction_ffs, main = "prediction with model selection")
plot(AOA_with_tuning$AOA, main = "AOA with model selection")

Conclusion

caret has no functions that explicitly deal with spatial data. However, due to its rather flexible design, caret is compatible with several packages designed for spatial machine-learning. The caret::trainControl() takes a list of CV indices as input, which makes it quite flexible to work with the output of e.g. CAST::knndm() and blockCV::cv_spatial(). Furthermore, caret is easy to use due to its functional programming paradigm. The documentation is extensive, and it’s rather easy to find modelling algorithms and their hyperparameters. Lastly, it should be noted that caret is not actively developed, since its main developer moved to tidymodels.