Random Forest mapping for MODIS tile h18v03 using MODIS and ERA5 precipitation
This notebook documents my solution for the Kaggle competition Modeling monthly soil surface water content. 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 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 target and covariate data are described on the competition page and in Schrier et al. (2026).
Monthly total precipitation was obtained from ERA5 monthly averaged reanalysis for 2000–2024. 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 model used the ranger implementation 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"
)
x and y indicates strong reliance on spatial gradients.geom layer was incomplete and required median imputation.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.