Random Forest modeling with MODIS and ERA5 precipitation
This notebook documents my solution for the Kaggle competition Modeling monthly soil surface water content (Parente et al. 2026). The objective was to predict monthly soil surface water content (sswc_m) for MODIS tile h18v03 from station observations, spatial coordinates, time, and Earth observation covariates.
The final solution used a Random Forest model with MODIS land-surface temperature and vegetation predictors, spatial and seasonal features, ERA5 monthly precipitation, and precipitation from the preceding month. Model performance was evaluated using five-fold spatial cross-validation in which complete monitoring locations, rather than individual rows, were assigned to folds.
The final submission achieved a Kaggle public leaderboard RMSE of 34.36176.
The supplied competition data contained:
sswc_m, representing monthly soil surface water content in \(\mathrm{cm^3/cm^3} \times 1000\).The gridded MODIS predictors for tile h18v03 were obtained from the published MODIS MOD11A2 data archive (Hengl 2026).
The original predictors were:
| Variable | Description |
|---|---|
x, y |
MODIS Sinusoidal coordinates in metres |
dlst |
Monthly mean daytime land-surface temperature |
dlst.95 |
Monthly daytime LST 95th percentile |
nlst |
Monthly mean nighttime land-surface temperature |
nlst.95 |
Monthly nighttime LST 95th percentile |
geom |
Monthly geometric-temperature covariate |
evi |
Monthly MODIS Enhanced Vegetation Index |
YEARMON |
Observation month and year |
The competition data and covariates are described on the competition page (Parente et al. 2026). The target soil water content observations originate from the European station dataset described by Schrier et al. (2026).
Monthly total precipitation was obtained from the ERA5 monthly averaged reanalysis for 2000–2024 (Hersbach et al. 2020). The ERA5 data were supplied in geographic coordinates (WGS 84) at \(0.1^\circ\) resolution. Station coordinates were transformed from the MODIS Sinusoidal projection to the ERA5 coordinate reference system, and precipitation was extracted using bilinear interpolation.
The precipitation variable was converted to monthly millimetres and two features were retained by the final model:
rain_mm: precipitation in the current month;rain_lag1: precipitation in the preceding month.days_per_month <- lubridate::days_in_month(era5_dates)
for (i in seq_len(terra::nlyr(era5_rain))) {
era5_rain_mm[[i]] <-
era5_rain[[i]] * 1000 * days_per_month[i]
}
rain_long <- rain_long |>
arrange(station_id, date) |>
group_by(station_id) |>
mutate(rain_lag1 = lag(rain_mm, 1)) |>
ungroup()The observed response was concentrated around moderate SSWC values, with a smaller number of very high observations.
The observations showed a pronounced annual cycle. SSWC was generally higher during winter and lower during late spring and summer.
Training and test locations were spatially separated. Consequently, a random row-level validation split would overestimate the ability of the model to predict at new locations.
The YEARMON field was converted to a date. Year, month, cyclic seasonal terms, and temperature contrasts were then derived.
prepare_features <- function(data) {
data |>
mutate(
date = lubridate::ymd(paste0(YEARMON, "-01")),
year = lubridate::year(date),
month = lubridate::month(date),
month_sin = sin(2 * pi * month / 12),
month_cos = cos(2 * pi * month / 12),
lst_difference = dlst - nlst,
lst95_difference = dlst.95 - nlst.95,
mean_lst = (dlst + nlst) / 2,
mean_lst95 = (dlst.95 + nlst.95) / 2,
station_id = paste(round(x, 6), round(y, 6), sep = "_")
)
}The final predictor set contained 17 variables:
final_predictors <- c(
"x", "y",
"dlst", "dlst.95", "nlst", "nlst.95",
"geom", "evi", "year",
"month_sin", "month_cos",
"lst_difference", "lst95_difference",
"mean_lst", "mean_lst95",
"rain_mm", "rain_lag1"
)Missing predictor values were replaced by medians calculated from the training partition only. For the final model, medians were recalculated from all training observations and then applied to the test data.
Unique station locations were randomly assigned to five folds using a fixed seed. All monthly observations belonging to the same station remained in the same fold. This design tests transfer to spatially unseen locations and reduces spatial information leakage.
For observed values \(y_i\) and predictions \(\hat{y}_i\), RMSE was calculated as
\[\mathrm{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2}.\]
MAE and \(R^2\) were reported as complementary metrics, although Kaggle ranking was based on RMSE.
The selected Random Forest model used the ranger implementation in R (Wright and Ziegler 2017), with:
mtry = 5;min.node.size = 10;set.seed(2026)
final_era5_model <- ranger::ranger(
formula = sswc_m ~ .,
data = final_training_table,
num.trees = 1000,
mtry = 5,
min.node.size = 10,
splitrule = "variance",
importance = "permutation",
seed = 2026
)oof_metrics |>
mutate(value = round(value, 3)) |>
kable(
col.names = c("Metric", "Spatial OOF value"),
caption = "Pooled spatial out-of-fold performance."
)| Metric | Spatial OOF value |
|---|---|
| RMSE | 36.014 |
| MAE | 26.722 |
| R_squared | 0.453 |
The pooled spatial out-of-fold evaluation produced RMSE = 36.01, MAE = 26.72, and \(R^2\) = 0.453. The model therefore explained approximately 45% of the spatially held-out variation.
fold_metrics |>
mutate(
across(c(RMSE, MAE, R_squared), ~ round(.x, 3))
) |>
kable(
col.names = c(
"Fold", "RMSE", "MAE", "R²", "Validation rows"
),
caption = "Performance in each spatial validation fold."
)| Fold | RMSE | MAE | R² | Validation rows |
|---|---|---|---|---|
| 1 | 38.713 | 27.882 | 0.486 | 30453 |
| 2 | 38.212 | 27.888 | 0.330 | 31317 |
| 3 | 36.665 | 27.347 | 0.457 | 31055 |
| 4 | 35.528 | 27.530 | 0.455 | 30511 |
| 5 | 30.400 | 23.016 | 0.524 | 31258 |
The fold RMSE ranged from approximately 30.4 to 38.7, demonstrating that prediction difficulty varied geographically.
The model captured the main gradient but showed regression toward the mean: low observations were often overestimated, whereas extreme high values were underestimated. This is an important limitation for interpreting extreme wetness conditions.
model_comparison |>
mutate(mean_RMSE = round(mean_RMSE, 2)) |>
kable(
col.names = c("Model", "Mean spatial CV RMSE"),
caption = "Comparison of the best MODIS-only and ERA5-enhanced models."
)| Model | Mean spatial CV RMSE |
|---|---|
| MODIS Random Forest | 36.7 |
| MODIS + ERA5 rainfall | 35.9 |
Adding ERA5 precipitation reduced mean spatial CV RMSE from approximately 36.7 to 35.9, an improvement of about 0.8 RMSE units (approximately 2.2%).
The northing coordinate (y) was the most influential predictor, indicating a strong latitudinal moisture gradient within the tile. geom, upper daytime temperature (dlst.95), easting (x), mean daytime temperature, and current-month rainfall were also influential. The previous-month rainfall contributed less than current rainfall but was retained because the combined rainfall feature set improved spatial validation.
Spatial coordinates are useful predictive covariates, but their high importance also indicates that part of the model performance is based on broad spatial gradients rather than purely process-based relationships.
Continuous maps were generated from full-grid MODIS rasters for tile h18v03. Monthly daytime and nighttime mean and 95th-percentile LST rasters, EVI, and monthly geometric temperature were aligned to the 926.6 m MODIS grid. ERA5 precipitation was projected from WGS 84 to the MODIS Sinusoidal projection using bilinear interpolation.
The final maps used a constant land mask containing 728,588 pixels, ensuring that monthly summary values were calculated for the same spatial support.
The supplied April geometric-temperature raster contained only 32 valid pixels. Missing April geom values were therefore replaced using the training-data median, consistently with the model preprocessing procedure. This limitation should be considered when interpreting the April map.
map_statistics |>
select(YEARMON, minimum, maximum, mean, valid_pixels) |>
mutate(
across(c(minimum, maximum, mean), ~ round(.x, 2))
) |>
kable(
col.names = c(
"Month", "Minimum", "Maximum", "Mean", "Valid pixels"
),
caption = "Summary statistics for the 2024 continuous prediction maps."
)| Month | Minimum | Maximum | Mean | Valid pixels |
|---|---|---|---|---|
| 2024-01 | 94.62 | 476.75 | 158.27 | 728588 |
| 2024-02 | 91.43 | 463.29 | 153.70 | 728588 |
| 2024-03 | 80.09 | 472.85 | 142.93 | 728588 |
| 2024-04 | 87.01 | 465.85 | 139.96 | 728588 |
| 2024-05 | 58.91 | 456.71 | 114.26 | 728588 |
| 2024-06 | 51.19 | 463.46 | 120.03 | 728588 |
| 2024-07 | 64.01 | 453.51 | 127.62 | 728588 |
| 2024-08 | 55.57 | 462.61 | 118.20 | 728588 |
| 2024-09 | 61.74 | 467.40 | 121.13 | 728588 |
| 2024-10 | 63.67 | 471.61 | 140.29 | 728588 |
| 2024-11 | 81.96 | 475.94 | 146.61 | 728588 |
| 2024-12 | 94.70 | 468.31 | 158.02 | 728588 |
The wettest mapped month by spatial mean was January (158.27), closely followed by December (158.02). The lowest spatial mean occurred in May (114.26). The difference between the wettest and driest monthly means was approximately 44.01 units, and the annual mean of monthly spatial means was 136.75.
The final model was trained on all available training observations and used to predict the 73,998 rows of test.csv. Predictions were written in the same order as sample_submission.csv, checked for missing and negative values, and saved as submission_era5_candidate.csv.
submission_era5_candidate <- sample_submission |>
mutate(sswc_m = pmax(final_prediction, 0)) |>
select(ID, sswc_m)
write_csv(
submission_era5_candidate,
"submission_era5_candidate.csv"
)
Several model configurations and predictor sets were tested during the hackathon. The initial baseline was a Random Forest model based on the original MODIS predictors, spatial coordinates and temporal features.
The original predictors included:
x and y spatial coordinates;dlst);nlst);dlst.95 and nlst.95);geom);evi);Additional temporal and temperature-derived features were created:
Different Random Forest configurations were evaluated by changing mtry and min.node.size. The best MODIS-only configuration produced a mean spatial cross-validation RMSE of approximately 36.7.
ERA5 monthly total precipitation was then added to the predictor set. Two rainfall configurations were compared:
Both ERA5 configurations produced a mean spatial cross-validation RMSE of approximately 35.8. The simpler configuration containing current-month rainfall and a one-month lag was selected because it achieved similar performance with fewer predictors.
The final Random Forest configuration used:
mtry = 5;min.node.size = 10;The final model was trained using the complete training dataset and was used to create submission_era5_candidate.csv.
Several additional environmental predictors were considered during model development but were not included in the final Kaggle submission.
Copernicus CLMS soil moisture data were investigated as a potentially important predictor. However, their temporal coverage, raster structure and missing values required additional preprocessing and validation. Because of the limited hackathon time, these data were not included in the final model.
Catchment area from GEDTM30 was also considered as a static topographic predictor. Although the raster was prepared and values were extracted, its contribution was not fully evaluated using the same five-fold spatial cross-validation procedure. Therefore, it was excluded from the final submission to avoid adding an unvalidated predictor.
Atmospheric water vapour was suggested as another potentially useful covariate. It was not tested in the final workflow because of time and data-processing constraints.
The model has several additional limitations:
y, had high variable importance. This suggests that the model partly learned broad geographical gradients rather than only physical relationships.Future improvements could include fully validated CLMS soil moisture, catchment area, atmospheric water vapour and other terrain or soil properties. Gradient boosting models or ensembles could also be compared with Random Forest using the same spatial validation folds.
This project produced a complete spatial-temporal workflow for monthly soil surface water content prediction across MODIS tile h18v03. A Random Forest model combining MODIS temperature and vegetation data, coordinates, seasonal features, and ERA5 rainfall achieved a mean spatial CV RMSE of approximately 35.9 and a Kaggle public RMSE of 34.36176.
ERA5 precipitation produced a modest but consistent improvement over the MODIS-only model. The final continuous maps reproduced a clear seasonal cycle, with higher mean SSWC during winter and lower values in late spring and summer. The results also demonstrate the importance of spatially structured validation and input-raster quality control in Earth observation machine-learning workflows.