Module 4: Classical ML for EO, Done Rigorously

Phase: II, Core analytics
Level: Intermediate to advanced
Estimated time: ~12 hours (about 4 h concepts, 6 h lab, 2 h write-up)
Prerequisites: Phase I (Modules 1 to 3), Phase 0.
Portfolio thread: the foundation of Portfolio Project B. You build a forest/non-forest classifier with honest spatial cross-validation and area estimates with confidence intervals.

Why this module

Most production Earth-observation systems still run on gradient boosting over engineered features, not transformers. The thing that makes such a model credible is not the algorithm, it is the evaluation. Spatial autocorrelation makes random train/test splits lie, often by ten or twenty accuracy points, and a candidate who shows spatially blocked cross-validation and area estimates with confidence intervals is instantly more trustworthy than one quoting a shiny random-split number. This module teaches the classical modeling that quietly powers the field, and the sampling and validation discipline that separate defensible results from optimistic ones.

Learning outcomes

After this module you can:

1. Engineer features from Sentinel-1 and Sentinel-2 for pixel classification (temporal statistics, indices, ratios, texture, terrain).
2. Explain why random k-fold cross-validation leaks under spatial autocorrelation, and implement spatial block cross-validation.
3. Train and tune a gradient-boosted classifier and interpret it with feature importance and SHAP.
4. Design a stratified sample and report area estimates with confidence intervals.
5. Scale inference across an AOI with dask and write the result as a COG.
6. Track experiments reproducibly with MLflow.

Concept lessons (about 4 hours)

Lesson 1: Why gradient boosting still wins for tabular EO

When your data is a table of per-pixel or per-parcel features, gradient-boosted trees (XGBoost, LightGBM, or scikit-learn's HistGradientBoosting) are usually the strongest, fastest, and most interpretable option, and they behave well with limited labels. Deep learning earns its place when spatial context or raw imagery matters (Phase III), but reaching for a transformer when a boosted-tree model on good features would do is a common and expensive mistake. Know both, and choose deliberately.

Lesson 2: Features for Sentinel-1 and Sentinel-2

A good feature stack turns raw scenes into signal:

• Optical (Sentinel-2 or HLS): seasonal composites (dry and wet season medians), plus indices such as NDVI and NBR from Module 2. Vegetation phenology is highly discriminative.
• SAR (Sentinel-1): annual median VV and VH, the VH/VV ratio, the Radar Vegetation Index, and texture measures. SAR adds structure and all-weather signal that optical alone misses.
• Terrain (DEM): elevation, slope, and aspect, which shape both vegetation and radar response.
• Temporal statistics: mean, median, percentiles, and amplitude across a year capture how a pixel changes, which a single date cannot.

Standardize per scene where appropriate, and stack everything into one aligned feature array (reuse cube.py and indices.py). The habit from earlier modules still applies: know what each feature physically represents.

Lesson 3: Spatial cross-validation, the credibility lever

Nearby pixels are similar (spatial autocorrelation), so a random train/test split routinely puts a pixel in training and its near-duplicate neighbor in testing. The model then looks excellent on the test set and fails on genuinely new ground. The fix is to split by space, not at random:

• Spatial block cross-validation: partition the AOI into blocks (for example 0.1 degree tiles) and hold out whole blocks, so training and testing never share a neighborhood.
• Grouped folds: assign each sample a block id and use grouped or stratified-grouped k-fold so a block is never split across folds.
• Area of applicability: after training, estimate where in feature space the model is actually trustworthy, so you do not extrapolate blindly.

Always report both the random and the spatial number, and treat the spatial one as the truth. The gap between them is itself a finding, and closing it (with more geographically diverse labels) is usually where the real gains are.

Lesson 4: Sampling design and area estimation

To make a defensible area statement (for example hectares of forest), you need a probability sample of reference labels, usually stratified by mapped class, and then the Olofsson area-weighted estimators with confidence intervals from Module 2. Handle class imbalance thoughtfully (class weights or careful thresholding rather than naive oversampling that leaks across the spatial split).

Lesson 5: MLOps hygiene from day one

Track every experiment with MLflow (or an equivalent): parameters, metrics, and artifacts such as the model and the feature list. This makes results reproducible, comparisons honest, and your work legible to a team. Starting this habit now pays off through the deep-learning and production phases.

Guided lab (about 6 hours): a spatially validated forest classifier

Step 1: Build the feature stack

Over your tropical AOI (for example a Niger Delta or Cross River tile), assemble: Sentinel-1 annual median VV, VH, and VH/VV ratio; Sentinel-2 or HLS dry- and wet-season composites plus NDVI and NBR; and DEM slope and elevation. Reuse cube.py and indices.py. Stack into an aligned array and sample it at your label locations.

Step 2: Labels with a stratified design

Sample forest and non-forest labels from Hansen Global Forest Change (tree canopy cover in 2000 with a forest threshold, and the loss-year layer to avoid labeling recently cleared pixels as forest), using stratified random sampling so both classes are well represented. Keep the label points as a GeoDataFrame with coordinates.

Step 3: Spatial blocks and folds

"""ml.py: spatial cross-validation helpers."""
from __future__ import annotations
import numpy as np

def assign_blocks(lons, lats, block_deg=0.1):
    """Assign each point to a spatial block id (a grid cell of block_deg degrees)."""
    bx = np.floor(np.asarray(lons, float) / block_deg).astype(int)
    by = np.floor(np.asarray(lats, float) / block_deg).astype(int)
    return np.array([f"{a}_{b}" for a, b in zip(bx, by)])

# tests/test_ml.py
from eo_portfolio.ml import assign_blocks

def test_blocks_group_nearby_and_separate_far():
    b = assign_blocks([7.501, 7.509, 7.62], [5.001, 5.004, 5.10], block_deg=0.1)
    assert b[0] == b[1]      # within the same 0.1 degree block
    assert b[0] != b[2]      # a different block

Use these block ids as groups in scikit-learn's StratifiedGroupKFold, so no block is split across folds and class balance is preserved.

Step 4: Train, compare, interpret

Train a gradient-boosted classifier. Evaluate it two ways on the same data: ordinary StratifiedKFold (random) and StratifiedGroupKFold on the spatial blocks. Report both and show the gap. Then compute feature importance and SHAP values, and run an ablation removing the Sentinel-1 features versus the Sentinel-2 features, explaining the result physically.

Step 5: Area estimate, scaled inference, tracking

Draw a stratified reference sample, build the confusion matrix, and compute area with confidence intervals using the accuracy.py helper from Module 2. Then run inference across the full AOI with dask and write the forest mask as a COG. Log parameters, metrics, and artifacts to MLflow throughout.

Step 6: Write-up and commit

A short report with the random-versus-spatial accuracy comparison, the area estimate with its 95 percent interval, the feature importance and ablation, and an MLflow run link. Commit ml.py, tests, and the notebook.

git add src/eo_portfolio/ml.py tests/test_ml.py notebooks/06_forest_classifier.ipynb
git commit -m "Module 4: forest classifier with spatial block CV and area estimates"
git push

Checkpoint

Self-check (answers below).

1. Why does a random train/test split overstate accuracy on spatially autocorrelated imagery?
2. Your random-fold accuracy is 94 percent and your spatial-block accuracy is 81 percent. Which do you report, and what does the gap tell you?
3. Why report an area estimate with a confidence interval rather than just counting classified pixels?
4. When would you prefer a gradient-boosted model over a deep network for an EO task?

Interview-style questions (practice out loud).

• Design a spatial cross-validation scheme for a country-scale land-cover model and justify it.
• What is the "area of applicability" and why does it matter for a deployed map?
• Your Sentinel-1 features dominate the importance ranking over Sentinel-2. Give two physical reasons that could be legitimate and one that would be a red flag.
• How do you estimate forest area with uncertainty for a report a non-specialist will act on?

Answers. (1) Nearby pixels are similar, so a random split places near-duplicate neighbors in both train and test, and the model is effectively tested on data it has almost seen.
(2) Report the spatial number (81); the gap is the amount your random split was inflated by leakage and signals limited geographic generalization, usually fixed with more spatially diverse labels.
(3) Classified pixel counts are biased by classification error; the Olofsson area-weighted estimate with a confidence interval gives a defensible number and its uncertainty.
(4) When the signal is well captured by engineered tabular features, labels are limited, interpretability matters, or compute is constrained; deep learning wins when spatial context or raw imagery drives the task.

Deliverable

A tested ml.py (spatial block assignment), a notebooks/06_forest_classifier.ipynb that trains a gradient-boosted forest/non-forest classifier, reports random versus spatial-block accuracy, gives an area estimate with a 95 percent confidence interval, includes SHAP and an S1-versus-S2 ablation, and logs to MLflow. This is the foundation of Portfolio Project B.

What a hiring manager sees

A model is easy; honest evaluation is rare. Showing spatial block cross-validation, reporting the spatial number as the truth, quantifying area with confidence intervals, and interpreting the model with SHAP is exactly the competency bar for this field. The random-versus-spatial gap in your write-up, presented plainly, tells a reviewer you understand the single most common way EO models mislead, which is worth more than a higher headline accuracy.

Currency note

Verified August 2026. Hansen Global Forest Change is currently version 1.13 (GFC-2025-v1.13), covering 2000 to 2025 at 30 m, with annual tree-cover loss and the year-of-loss layer, and is available in the Earth Engine catalog (UMD/hansen/global_forest_change_2025_v1_13) and the Hansen download storage; the prior GFC-2024-v1.12 (2000 to 2024) is still published. A version 2.0 reprocessing is planned with no confirmed date, so check for the latest release when you build labels. XGBoost, LightGBM, scikit-learn (StratifiedGroupKFold), SHAP, and MLflow are stable and current; confirm API details and pin versions in your Module 0.1 environment file. For "area of applicability," see the Meyer and Pebesma method and its current implementation.