Predicting vegetation across 133 global ecoregions with 93% accuracy

Overview

Desertflow Global Bundle is a noncommercial research project studying how vegetation responds to climate, soil, fire, and human-driven stressors across the world's terrestrial ecoregions. The core question: can per-region machine-learning specialists, combined with a principled cross-region router, beat traditional global vegetation models — and do so without secretly training on the answer?

The system fetches per-cell, per-year aggregates for 133 World Wildlife Fund (WWF) terrestrial ecoregions across 22 years (2002–2024), trains a CatBoost regression specialist for each region in a 89-region training split, then combines those specialists at prediction time via a per-cell Bayesian-ridge router. Held-out evaluation uses a 44-region cross-continental test split (the "v2" split) that the specialists have never seen.

The framework, all driver modules, the validation harness, and trained model artifacts are being prepared for open-source release alongside an EarthArXiv preprint.

Why this matters

Desertification and land degradation affect over 3.2 billion people and cost the global economy an estimated $490 billion per year (UNCCD figures). Accurate, regionally-tuned vegetation forecasting is foundational for:

• Land Degradation Neutrality (SDG 15.3): tracking SDG indicator 15.3.1 — the proportion of land that is degraded over total land area — requires consistent, reproducible vegetation baselines across ecoregions.
• Climate adaptation planning: identifying which regions are most vulnerable to vegetation collapse under projected 2030+ climate scenarios.
• Conservation prioritization: directing limited resources toward ecosystems with the highest combined risk and recovery potential.
• Policy verification: independent, methodologically-transparent baselines for verifying claims under the UN Convention to Combat Desertification (UNCCD) and REDD+ frameworks.

This project does not aim to replace existing operational systems (MODIS-based monitoring, Copernicus land services, etc.). It aims to provide a methodologically clean reference baseline that other studies can compare against — one that has been explicitly audited for the leakage failure modes that quietly inflate accuracy in much of the published NDVI forecasting literature.

The problem: cross-region NDVI forecasting

Predicting NDVI (Normalized Difference Vegetation Index) is the standard benchmark task for vegetation modeling. The challenge is that vegetation response to climate varies enormously between ecoregions:

• Sahara: NDVI is near-zero year-round; a single rainfall event causes brief, sharp greening.
• Amazon basin: NDVI is high and seasonally damped; the dominant signal is cloud cover, not climate.
• Yakutian larch forest: a strong but short summer growing season; winter NDVI is meaningless under snow.
• Indo-Gangetic plain: NDVI is dominated by human agriculture, not natural climate response.

Global vegetation models that train one set of weights across all biomes typically achieve cross-region R² in the 0.6–0.7 range — the global average smooths over what's actually happening in any one place. Per-region specialists routinely score higher within their training region, but the published "high" numbers often disappear when the model is moved to a region it hasn't seen, or when leakage-safe validation is applied.

Desertflow targets the middle ground: per-region specialists combined with a principled router that handles regions the specialists weren't trained on, evaluated on a held-out 44-region cross-continental test set with no NDVI-derived inputs.

System architecture

The system uses a two-layer architecture, internally referred to as the "B2" router pattern:

Key design choices:

• Specialists are isolated. Specialist N never sees data from specialist M's region. This eliminates the most common form of cross-region leakage.
• The router is the only place cross-region information is combined, and it does so via per-cell weights — not by sharing training data.
• Router calibration is also leakage-safe: it uses a calibration year window separate from the prediction year, so the router never sees future labels at calibration time.

Project details

• Spatial coverage: 133 WWF terrestrial ecoregions spanning approximately 45 million km² of land surface — including the Sahara, Sahel, Amazon basin, Congo basin, Indo-Gangetic plain, Siberian taiga, Patagonian steppe, Australian outback, Mongolian steppe, and 124 other ecoregions across all inhabited continents.
• Temporal coverage: 22 continuous years (2002–2024), monthly cadence, enabling decadal-scale trend analysis.
• Spatial resolution: ~4 km native grid (MODIS-derived), aggregated per ecoregion for training.
• Target variable: V11_mean_ndvi — region-mean MODIS NDVI per cell-year, computed from MOD13Q1.
• Driver candidates: ~250 leakage-safe candidate features per cell, of which ~33 final drivers were selected after rigorous leakage auditing.
• Train/test split: 89 training regions / 44 cross-continental held-out test regions (the "v2" split).
• Validation rule: the "2030 rule" — no model is allowed any feature that wouldn't be available for a 2030 forward prediction (rules out NDVI-derived features like LAI, fPAR, GPP, and most phenology metrics).

Methodology

Training proceeds in five clearly separated phases:

1. Per-region data assembly. For each of the 89 training regions, the Earth Engine pipeline assembles a parquet file with all cell-year combinations and all 19 driver modules. Modules are pre-filtered for leakage at this stage.
2. Per-region specialist training. Each region trains an independent CatBoost regressor. Optuna runs ~30 hyperparameter trials per region. The specialist is then frozen.
3. Router calibration. On a calibration year window inside the training regions, the router learns how to weight specialists for a given cell's environmental profile. Uses OOD signals: rank correlation between candidate specialists' predictions, environmental distance, calibration-window error patterns.
4. Held-out evaluation. On the 44 held-out test regions, the router produces per-cell weights without ever seeing test-region training data. The final R² is computed cell-by-cell, not region-mean.
5. Audit pass. Every input column is checked against the leakage-feature list (leakage_features.yaml) and the audit warning is logged if any leak slips through.

Driver modules (19)

The 19 modules group into five families. Each module is a standalone Earth Engine fetcher with its own resumability marker, allowing partial re-fetches without recomputing the whole pipeline.

Notably not in the driver list (deliberately excluded by the leakage audit): LAI, fPAR, GPP, MODIS land-cover classification, sister vegetation cover percentages, NDVI-derived phenology start-of-season day-of-year, and any other feature derived from NDVI or NDVI-equivalent optical bands.

Validation protocol

Validation is the heart of this work. The protocol enforces four rules in sequence:

1. The 2030 rule. No feature is allowed unless it would be computable for a 2030 forward prediction without seeing 2030 NDVI. This rules out everything NDVI-derived.
2. V/O/C/P checks — four explicit leakage gates:
   • V (Vegetation circularity): blocks features computed from NDVI time series (phenology start-of-season day, GPP, fPAR).
   • O (Optical bands): blocks raw MODIS red/NIR bands and indices derived from them.
   • C (Classifier): blocks features derived from supervised vegetation classifiers (MCD12Q1 land cover, tree cover %).
   • P (Pheno-prior): blocks any feature whose value depends on knowing the target year's NDVI peak.
3. Cross-region split. Specialists are trained on 89 regions; the 44 held-out test regions are never used in training, hyperparameter tuning, or router calibration.
4. Temporal split. Within each region, calibration and prediction year windows are separated, so the router never sees the prediction year's labels at calibration time.

Every column loaded by the framework is checked against leakage_features.yaml at training time. Matched columns are auto-set to "(ignore)" with a FAIL-2030 marker, and a banner alerts the operator that leakage features were detected and skipped.

Technology stack

Leakage-aware validation: a methodological contribution

The distinguishing contribution of this work is a documented leakage-aware validation protocol. Many published NDVI forecasting studies inadvertently leak future information through one of four channels:

1. Sister vegetation features. NDVI-derived land cover, GPP, fPAR, LAI, and phenology metrics are all computed from the same MODIS optical bands as NDVI itself. Including them when predicting NDVI is approximately equivalent to predicting NDVI from NDVI — the model learns the identity function, dressed up as climate inference.
2. Shared cross-region statistics. If feature scalers, gap-fillers, or normalization layers are fit across all regions before splitting train/test, statistics from the test regions leak into the training signal.
3. Improperly designed cross-validation splits. Random k-fold on cells within a region leaks neighboring cells' values; random splits across years leak seasonal patterns; even some "spatial" splits leave too-close test cells.
4. Pheno-prior contamination. Some published features (start-of-season day-of-year, peak NDVI day) are computed from the target NDVI series itself but presented as if they were independent inputs.

Bundle Research documents each failure mode and tests against it before reporting metrics. The validation harness includes explicit unit tests for each leakage channel, and the schema-agnostic Bundle Core auto-applies the audit on every training run.

Results

Headline metrics on the v2 held-out test split (44 regions, not seen during training or router calibration):

The +0.27 lift from the bundle architecture over specialist-only baseline is the empirical justification for the two-layer design. Per-region specialists trained in isolation score essentially no better than a standard global model (0.66 vs ~0.67) — neither approach generalizes cross-continent on its own. The bundle's Bayesian-ridge router is what turns 89 narrow specialists into a system that matches within-region accuracy on regions it has never seen.

Remaining challenges: regions still closest to the floor

In the original baseline B2 ridge configuration, 9 of the 44 held-out regions fell below R² = 0.70 — the "hard regions" problem. After moving to deeper CatBoost specialists (depth-12) combined with the OOD-weighted router, every region now exceeds R² = 0.80. The 9-region failure tier is gone.

The three regions still closest to the floor are below — all at R² ≈ 0.82–0.84. They are not "failures" but they are where the next round of methodological work is targeted:

Improvement from the previous baseline on these three: borneo_rainforest (+0.003), iberian_wildfire (+0.021), mesopotamia (+0.016). Modest absolute gains, but enough to push all three above 0.80. The E26 anomaly modules (VPD extremes, SST anomalies, snow phenology) currently being fetched are targeted at the same kind of regions — anomaly drivers may help where mean climate isn't predictive enough.

Recent progress

Earth Engine compute requirements

This work requires sustained Earth Engine compute for large-scale ecoregion processing. A single full bundle refresh requires approximately 30,000 batch export tasks across 133 regions × 22 years × 19 driver modules. Per-task EECU consumption ranges from ~0.05 (lightweight modules like VPD extremes) to ~0.15 (snow phenology). A full retrain consumes approximately 860 EECU-hours.

The project is registered under Google Earth Engine's noncommercial research program (project ID: celtic-science-493107-c6). The Community Tier monthly EECU budget (150 hours) covers roughly one-sixth of a full refresh, which is why iterative experiments are blocked at the current tier and a Contributor / Partner Tier uplift application is in review.

Roadmap

References & inspirations

Key data sources and methodological references:

• WWF Terrestrial Ecoregions of the World — Olson et al., 2001 — used for the 133-region biogeographic stratification.
• MOD13Q1 Vegetation Indices 16-Day L3 Global 250m — Didan, 2015 — primary NDVI source.
• ERA5 reanalysis — Hersbach et al., 2020 — primary climate driver source.
• CHIRPS precipitation — Funk et al., 2015 — daily precipitation aggregates.
• TerraClimate — Abatzoglou et al., 2018 — monthly climatic water deficit.
• SoilGrids 2.0 — Poggio et al., 2021 — soil texture and organic carbon.
• Google Earth Engine: Planetary-scale geospatial analysis for everyone — Gorelick et al., 2017 — the platform.
• CatBoost: gradient boosting with categorical features support — Prokhorenkova et al., 2018 — the chosen base learner.
• UN Convention to Combat Desertification (UNCCD) — policy framing for the SDG 15.3 connection.

Full citations and additional methodological references will be in the EarthArXiv preprint.

Frequently asked questions

Acknowledgments

Data: NASA (MOD13Q1, MCD64A1, MOD10A1, MCD12Q1, SMAP), ECMWF (ERA5), UCSB CHG (CHIRPS), TerraClimate, ISRIC (SoilGrids), WWF (Ecoregions). Compute: Google Earth Engine noncommercial research program. Software: open-source contributors to CatBoost, scikit-learn, PySide6, Optuna, rclone, and pyarrow.

This is independent research with no institutional affiliation; the work is published openly so others can verify, criticize, and extend the methodology.

Get in touch

For collaboration, citations, mentorship offers, or any questions about the methodology — bagbagig@bundle-research.xyz (click to copy)

Code: github.com/mohtheghost