# install.packages(c("terra", "sf", "dplyr", "ggplot2", "usdm", "caret", "randomForest", "CAST"))
library(terra) # raster handling
library(sf) # vector/point handling
library(dplyr)
library(ggplot2)
library(usdm) # VIF
library(caret) # model training wrapper
library(randomForest)
library(CAST) # Area of ApplicabilityGeospatial Machine Learning in R
Classification (Land cover) and Regression (Canopy height) Tasks
1 Background
In the GEE, we:
- Built a Sentinel-2 composite for the year 2026 within our roi
- Sampled a classification training set (forest, tea, built, water, and planted_forest) and a canopy-height regression training set (from the Meta canopy-height product (Tolan et al. 2024)), and
- Exported both the raster and the point tables (with a random 70/30 train/validation split already applied) to Google Drive.
This script picks up from those exports and does in R (R Core Team 2026) what is not trivial in GEE (Gorelick et al. 2017); checking predictor multi-collinearity, and mapping the Area of Applicability (AOA), i.e. where the trained model is safe to extrapolate to and where it isn’t. We also have the flexibility to use more algorithms and unlimited computation time. The AoA point matters a lot once these maps are used for economic analyses (e.g. land valuation, agricultural suitability, ecosystem service accounting): a prediction outside the AOA is a guess, not a measurement.
Tip
Why do we need this as agricultural/environmental economists?
To create new variables from satellite images.
To quantify how much of a given land cover exists in an area.
To assess change in land cover over time, among others.
2 Setup
Create an R Project and keep all the data from the GEE in a data folder within it. Specifically, the raster data composites for classification (landcover) and regression (canopy height) and their corresponding GeoJSON training.
We rely on terra (Hijmans et al. 2026) for raster data handling, sf (Pebesma and Bivand 2023) for vector data, usdm (Naimi et al. 2014) for multi-collinearity checks, caret (Kuhn 2007) for model setting, randomForest [Liaw and Wiener (2002)] for the model, and CAST (Meyer et al. 2026) for area of applicability checks. Both dplyr (Wickham et al. 2026) and ggplot2 (Wickham 2016) were for wrangling and visualization.
3 Part A — Land cover classification
3.1 Load the composite raster, training, and validation points
composite <- rast("data/composite_image.tif")
train_pts <- st_read("data/training_points.geojson")Reading layer `training_points' from data source
`C:\Users\Wyclife\R_PROJECTS\ml_milan_2026\data\training_points.geojson'
using driver `GeoJSON'
Simple feature collection with 354 features and 14 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 35.26485 ymin: -0.4331227 xmax: 35.35225 ymax: -0.3518252
Geodetic CRS: WGS 84valid_pts <- st_read("data/validation_points.geojson")Reading layer `validation_points' from data source
`C:\Users\Wyclife\R_PROJECTS\ml_milan_2026\data\validation_points.geojson'
using driver `GeoJSON'
Simple feature collection with 146 features and 14 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 35.26602 ymin: -0.4344702 xmax: 35.35279 ymax: -0.3510167
Geodetic CRS: WGS 84train_pts$label <- factor(train_pts$label)
valid_pts$label <- factor(valid_pts$label)3.2 Multi-collinearity among spectral bands
Some Sentinel-2 bands are highly correlated by nature (e.g. the visible bands, or the red-edge/NIR bands). RandomForest is fairly robust to collinear predictors for prediction, but collinearity still muddies interpretation (e.g. variable importance), so it’s worth checking.
bands <- names(composite)
band_vals <- st_drop_geometry(train_pts)[, bands]
vif_result <- vifstep(band_vals, th = 10) # flags predictors with VIF > 10
vif_result7 variables from the 11 input variables have collinearity problem:
B7 B8A B6 B11 B3 B4 B8
After excluding the collinear variables, the linear correlation coefficients ranges between:
min correlation ( B9 ~ B2 ): -0.09111023
max correlation ( B12 ~ B2 ): 0.6899148
---------- VIFs of the remained variables --------
Variables VIF
1 B2 4.495327
2 B5 4.496314
3 B9 2.523310
4 B12 2.078630If several bands are flagged, options are: drop the redundant ones. You can also replace the raw bands with a smaller set of indices (NDVI, NDWI, etc.), or run a PCA on the band stack and classify on the components instead.
3.3 Train a random forest classifier
set.seed(248)
ctrl <- trainControl(method = "cv", number = 5)
rf_class <- train(
# x = band_vals[, c("B2", "B5", "B9", "B12")],
x = band_vals,
y = train_pts$label,
method = "rf",
trControl = ctrl,
tuneLength = 3,
ntree = 500
)
rf_classRandom Forest
354 samples
11 predictor
5 classes: '0', '1', '2', '3', '4'
No pre-processing
Resampling: Cross-Validated (5 fold)
Summary of sample sizes: 283, 283, 284, 282, 284
Resampling results across tuning parameters:
mtry Accuracy Kappa
2 0.9746457 0.9682893
6 0.9746065 0.9682468
11 0.9661949 0.9577256
Accuracy was used to select the optimal model using the largest value.
The final value used for the model was mtry = 2.
Warning
Whether the model runs regression or classification depends on the class of y.
3.4 Validate
valid_vals <- st_drop_geometry(valid_pts)[, bands]
pred_valid <- predict(rf_class, newdata = valid_vals)
confusionMatrix(pred_valid, valid_pts$label)Confusion Matrix and Statistics
Reference
Prediction 0 1 2 3 4
0 26 1 0 1 0
1 0 35 1 0 0
2 0 0 28 0 1
3 0 0 0 28 0
4 0 0 0 0 25
Overall Statistics
Accuracy : 0.9726
95% CI : (0.9313, 0.9925)
No Information Rate : 0.2466
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.9656
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: 0 Class: 1 Class: 2 Class: 3 Class: 4
Sensitivity 1.0000 0.9722 0.9655 0.9655 0.9615
Specificity 0.9833 0.9909 0.9915 1.0000 1.0000
Pos Pred Value 0.9286 0.9722 0.9655 1.0000 1.0000
Neg Pred Value 1.0000 0.9909 0.9915 0.9915 0.9917
Prevalence 0.1781 0.2466 0.1986 0.1986 0.1781
Detection Rate 0.1781 0.2397 0.1918 0.1918 0.1712
Detection Prevalence 0.1918 0.2466 0.1986 0.1918 0.1712
Balanced Accuracy 0.9917 0.9816 0.9785 0.9828 0.98083.5 Area of Applicability (AoA or MESS)
trainDI_class <- trainDI(train = band_vals, variables = bands, verbose = FALSE)
aoa_class <- aoa(newdata = composite, trainDI = trainDI_class)
plot(aoa_class$AOA, main = "Area of Applicability — Classification")
Pixels where AOA == 0 are spectrally unlike anything in the training set. The classifier is extrapolating there, so those predictions should be treated with caution (or masked out) in any downstream analysis.
3.6 Generate the land cover map
land_cover_r <- predict(composite, rf_class, na.rm = TRUE)
class_colors <- c("darkgreen", "lightgreen", "red", "blue", "green")
class_names <- c("Forest", "Tea", "Built", "Water", "Planted forest")
levels(land_cover_r) <- data.frame(
id = seq_along(class_names),
class = class_names
)
plot(land_cover_r, col = class_colors, main = "Predicted Land Cover")
# plot(land_cover_r, main = "Predicted Land Cover (R)")3.7 Run zonal statistics
Assume we have some disticts/states/ecozones covering the mapped area and we want to know proportion of each of the landcover classes we predicted.
zones <- vect("data/zones.geojson")
plot(land_cover_r, col = class_colors, main = "Predicted Land Cover")
plot(zones, add = T, border = "black", lwd = 6)
text(zones, labels = zones$name, cex = 3)
a <- cellSize(land_cover_r, unit = "ha")
r <- c(land_cover_r, a)
x <- terra::extract(r, zones)
zone_names <- zones$name
x_summary <- x |>
mutate(names = zone_names[ID]) |>
group_by(names, class) |>
summarise(
area_ha = sum(area, na.rm = TRUE),
.groups = "drop"
)
x_summary# A tibble: 25 × 3
names class area_ha
<chr> <fct> <dbl>
1 A Forest 440.
2 A Tea 1081.
3 A Built 664.
4 A Water 0.437
5 A Planted forest 167.
6 B Forest 703.
7 B Tea 804.
8 B Built 116.
9 B Water 15.6
10 B Planted forest 424.
# ℹ 15 more rowsggplot(x_summary, aes(names, area_ha)) +
facet_wrap(~class, ncol = 2) +
geom_col() +
theme(text = element_text(size = 20)) +
labs(x = "Zone", y = "Area (ha)")
4 Part B — Canopy height regression
4.1 Load data
composite_canopy <- rast("data/composite_image_canopy.tif")
names(composite_canopy) <- bands
train_canopy <- st_read("data/training_points_canopy.geojson")Reading layer `training_points_canopy' from data source
`C:\Users\Wyclife\R_PROJECTS\ml_milan_2026\data\training_points_canopy.geojson'
using driver `GeoJSON'
Simple feature collection with 708 features and 14 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 35.26404 ymin: -0.435099 xmax: 35.35333 ymax: -0.3511964
Geodetic CRS: WGS 84valid_canopy <- st_read("data/validation_points_canopy.geojson")Reading layer `validation_points_canopy' from data source
`C:\Users\Wyclife\R_PROJECTS\ml_milan_2026\data\validation_points_canopy.geojson'
using driver `GeoJSON'
Simple feature collection with 292 features and 14 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 35.26395 ymin: -0.4350092 xmax: 35.35306 ymax: -0.3510167
Geodetic CRS: WGS 844.2 Train a random forest regressor. Be careful with column names in train_canopy.
bands <- names(composite_canopy)
train_vals_c <- st_drop_geometry(train_canopy)[, c(bands, "canopy_height")]
rf_reg <- train(
x = train_vals_c[, bands],
y = train_vals_c$canopy_height,
method = "rf",
trControl = ctrl,
tuneLength = 3,
ntree = 500
)
rf_regRandom Forest
708 samples
11 predictor
No pre-processing
Resampling: Cross-Validated (5 fold)
Summary of sample sizes: 567, 567, 566, 566, 566
Resampling results across tuning parameters:
mtry RMSE Rsquared MAE
2 3.831780 0.8257425 2.385577
6 3.820493 0.8269319 2.353346
11 3.858581 0.8237153 2.372213
RMSE was used to select the optimal model using the smallest value.
The final value used for the model was mtry = 6.4.3 Validate: RMSE and R²
valid_vals_c <- st_drop_geometry(valid_canopy)[, bands]
pred_c <- predict(rf_reg, newdata = valid_vals_c)
obs <- valid_canopy$canopy_height
rmse <- sqrt(mean((obs - pred_c)^2))
r2 <- cor(obs, pred_c)^2
cat("RMSE:", round(rmse, 2), "m\nR-squared:", round(r2, 3))RMSE: 3.97 m
R-squared: 0.805ggplot(data.frame(obs, pred_c), aes(obs, pred_c)) +
geom_point(alpha = 0.5) +
geom_abline(slope = 1, intercept = 0, color = "red") +
labs(x = "Observed canopy height (m)",
y = "Predicted canopy height (m)") +
theme_minimal()
4.4 Predict the canopy height map
canopy_pred_r <- predict(composite_canopy, rf_reg, na.rm = TRUE)
plot(canopy_pred_r, main = "Predicted Canopy Height (R)")
4.5 Area of Applicability
trainDI_reg <- trainDI(train = train_vals_c[, bands], variables = bands, verbose = FALSE)note: variables were not weighted either because no weights or model were given,
no variable importance could be retrieved from the given model, or the model has a single feature.
Check caret::varImp(model)note: No model and no CV folds were given. The DI threshold is therefore based on all training dataaoa_reg <- aoa(newdata = composite_canopy, trainDI = trainDI_reg)Computing DI of new data...Computing AOA...Finished!plot(aoa_reg$DI, main = "Dissimilarity Index — Canopy Height")
plot(aoa_reg$AOA, main = "Area of Applicability — Canopy Height")
4.6 Mean canopy height per zone
We could be interested in knowing which zone has more tall trees. Could be important for forest restoration, other than NDVI etc.
zonal_stats <- extract(canopy_pred_r, zones, fun = mean, na.rm = TRUE, bind = TRUE)
zonal_df <- as.data.frame(zonal_stats) |>
rename(mean_cht = lyr1) |>
select(-id) |>
arrange(desc(mean_cht))
zonal_df name mean_cht
1 E 13.165224
2 D 8.617190
3 B 7.135768
4 C 5.785033
5 C 5.785033
6 A 3.3849635 GEE vs. R: what each is good for
| Google Earth Engine | R | |
|---|---|---|
| Strength | Cloud-scale access to imagery, fast compositing/masking over huge areas | Statistical diagnostics, uncertainty quantification, reproducible reporting |
| Classifier/regressor | ee.Classifier.smileRandomForest (fast, but limited diagnostics) |
caret/randomForest (same algorithm family, richer tooling around it) |
| Multicollinearity checks | Not built in | usdm::vifstep, corrplot |
| Area of Applicability | Not available | CAST::aoa |
| Typical role in a pipeline | Data preparation and export | Modeling, diagnostics, and communication (e.g. this Quarto document) |
6 Suggested exercises
- Drop the bands flagged by
vifstepand re-run the classifier, does accuracy change? - Compare the AOA maps for the classification and regression models, do they mask the same areas?
- Try
method = "ranger"instead of"rf"incaret::trainand compare speed and accuracy. - Overlay the AOA mask on the predicted maps and report what fraction of the study area falls outside the AOA. This is a useful caveat when presenting these maps to policymakers.
7 References
Gorelick, Noel, Matt Hancher, Mike Dixon, Simon Ilyushchenko, David Thau, and Rebecca Moore. 2017. “Google Earth Engine: Planetary-Scale Geospatial Analysis for Everyone.” Remote Sensing of Environment, ahead of print. https://doi.org/10.1016/j.rse.2017.06.031.
Hijmans, Robert J., Andrew Brown, and Márcia Barbosa. 2026. Terra: Spatial Data Analysis. https://doi.org/10.32614/CRAN.package.terra.
Kuhn, Max. 2007. Caret: Classification and Regression Training. The R Foundation. https://doi.org/10.32614/cran.package.caret.
Liaw, Andy, and Matthew Wiener. 2002. Classification and Regression by randomForest. 2: 18–22. https://CRAN.R-project.org/doc/Rnews/.
Meyer, Hanna, Carles Milà, Marvin Ludwig, Jan Linnenbrink, and Fabian Schumacher. 2026. CAST: ’Caret’ Applications for Spatial-Temporal Models. https://doi.org/10.32614/CRAN.package.CAST.
Naimi, Babak, Nicholas a.s. Hamm, Thomas A. Groen, Andrew K. Skidmore, and Albertus G. Toxopeus. 2014. Where Is Positional Uncertainty a Problem for Species Distribution Modelling. 37: 191–203. https://doi.org/10.1111/j.1600-0587.2013.00205.x.
Pebesma, Edzer, and Roger Bivand. 2023. Spatial Data Science: With Applications in r. https://doi.org/10.1201/9780429459016.
R Core Team. 2026. R: A Language and Environment for Statistical Computing. https://www.r-project.org/.
Tolan, Jamie, Hung-I Yang, Benjamin Nosarzewski, et al. 2024. “Very High Resolution Canopy Height Maps from RGB Imagery Using Self-Supervised Vision Transformer and Convolutional Decoder Trained on Aerial Lidar.” Remote Sensing of Environment 300 (January): 113888. https://doi.org/10.1016/j.rse.2023.113888.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. https://ggplot2.tidyverse.org.
Wickham, Hadley, Romain François, Lionel Henry, Kirill Müller, and Davis Vaughan. 2026. Dplyr: A Grammar of Data Manipulation. https://doi.org/10.32614/CRAN.package.dplyr.