Probabilistic NDVI Forecasting from Sparse Satellite Time Series and Weather Covariates

AI summary

Overview

• Research area: Machine learning for remote sensing and precision agriculture, specifically probabilistic time-series forecasting of vegetation indices from satellite imagery combined with meteorological data.
• Technical level: Advanced. The paper assumes familiarity with transformer architectures, quantile/pinball loss, attention masking, and forecasting evaluation metrics such as CRPS and MASE.
• Scope in one sentence: The paper proposes and evaluates a transformer-based quantile forecasting framework that predicts field-level NDVI up to three clear-sky acquisitions (approximately 14 days) ahead from sparse, irregular Sentinel-2 observations plus past and future weather covariates across European ecozones.

What This Paper Is About

Crop monitoring from satellites is hampered because clouds block many Sentinel-2 overpasses, leaving NDVI measurements sparse and irregularly spaced in time rather than arranged on a clean daily or weekly grid. The authors ask whether a model can still produce useful short-term vegetation forecasts under such conditions, and whether it can output an uncertainty range rather than only a single number. Their goal is a field-level NDVI forecaster that combines past vegetation dynamics with weather information, predicts quantiles rather than point values, and explicitly accounts for the fact that "three steps ahead" can mean a different number of elapsed days for different samples.

Key Contributions

1. A transformer-based quantile model for sparse NDVI forecasting. The architecture decouples a history encoder (past NDVI plus past weather) from a future encoder (known future weather covariates), fusing the two representations to predict NDVI quantiles at q ∈ {0.1, 0.5, 0.9} for multiple future acquisition timestamps in parallel, with binary masks preventing attention from attending to cloudy, unobserved positions.
2. A temporal-distance weighted quantile loss. Because clear-sky observations are irregularly spaced, each future target is down-weighted by w_k = 1 / (1 + α·(τ_{t+k} − τ_t)) with α = 0.1, so the training objective prioritizes nearer, less uncertain targets while retaining supervision at longer horizons.
3. Cumulative and extreme-weather feature engineering. Nine derived features are computed from precipitation and temperature: three "between-target" aggregates over the variable-length interval between consecutive observations (cumulative rainfall, days below 10 °C, days above 30 °C) and six rolling-window aggregates of the same quantities over 7- and 14-day windows.
4. Validation across European ecozones and growing seasons. The model is compared against statistical, recurrent, convolutional, LLM-based, and transformer baselines, with ablation studies isolating the contribution of the two encoding branches, target history, and each engineered component.

Main Findings

• Best results across every reported metric: The proposed model reaches RMSE 0.096 ± 0.155, MAE 0.062 ± 0.074, WMAPE 0.130 ± 0.156, MASE 0.828 ± 0.989, CRPS 0.038 ± 0.050, and Pinball loss 0.021 ± 0.030, outperforming AutoARIMA, LSTM, RNN, DeepAR, InceptionTime, TimeLLM, PatchTST, and Chronos-2 on the clear-sky NDVI task.
• Lightweight relative to strong baselines: The model uses 2.16 M parameters and 111.97 MFLOPs. For comparison, TimeLLM uses 124.44 M parameters and 510.61 MFLOPs, Chronos-2 uses 28 M parameters and 316.87 MFLOPs, and PatchTST uses 3.18 M parameters and 722.42 MFLOPs. TimeLLM is deterministic, so its CRPS and Pinball loss are not reported.
• Gains are not merely capacity-driven: The model outperforms smaller recurrent baselines such as LSTM (RMSE 0.122, MAE 0.088), RNN (RMSE 0.134, MAE 0.099), and DeepAR (RMSE 0.116, MAE 0.075), while remaining lighter than the large LLM-based and transformer baselines.
• Differences are statistically significant: All pairwise comparisons against baselines using the Diebold–Mariano test with a Newey–West heteroskedasticity and autocorrelation consistent variance correction, with the truncation lag set to the three-step horizon, yield p < 0.001.
• Temporal loss weighting matters most among the two proposed components: With both temporal weighting (w_k) and feature engineering enabled, the model reaches RMSE 0.096, MAE 0.062, WMAPE 0.130, MASE 0.828, CRPS 0.038, Pinball 0.021. The weakest configuration in this ablation reports RMSE 0.099, MAE 0.064, WMAPE 0.134, MASE 0.854, CRPS 0.040, Pinball 0.022, with the two mixed configurations falling in between (RMSE 0.097 and 0.098). The paper states temporal weighting provides the largest gains, while feature engineering yields smaller but systematic improvements.
• Target history dominates; weather adds incremental value: In the input-modality ablation, removing target history causes the largest degradation. The full multimodal configuration reports RMSE 0.096, MAE 0.062, WMAPE 0.130, MASE 0.828, CRPS 0.038, Pinball 0.021, whereas the worst reported configuration in that table reaches RMSE 0.165, MAE 0.129, WMAPE 0.273, MASE 1.734, CRPS 0.078, Pinball 0.044. Future covariates alone are described as insufficient for competitive performance but contribute incrementally in the full setting.
• Climate-dependent performance degradation: Stratifying by aggregated Köppen–Geiger groups, MAE rises from 0.042 (Semi-arid, BSk) to 0.073 (Continental, Dfa+Dfb+Dfc+Dsb) and R² falls from 0.791 to 0.557. Mean bias is small and positive across all groups (0.003–0.011), indicating a secondary overestimation tendency. C-Med, C-Temp, and Continental climates show broader spread at intermediate NDVI levels, attributed to higher phenological variability and observation noise from persistent cloud cover and seasonal effects.
• Robustness to imperfect weather forecasts: When future covariates are perturbed from 0% to 20% in 5% increments, degradation is moderate up to 10% and increases more sharply beyond it, with RMSE rising by 2% and CRPS by 1.6% at 20% perturbation. This supports the choice of 10% as the base training perturbation.

Methodology in Plain English

The authors build on the GreenEarthNet dataset, which contains 24,061 spatio-temporal data cubes over Europe from 2017 to 2022. Each cube holds 30 cloud-masked Sentinel-2 images at 5-day intervals, each 128 × 128 pixels spanning 2.56 × 2.56 km, plus 150 daily meteorological observations. NDVI is computed from the near-infrared (B8A) and red (B04) bands, both at 20 m resolution. For each cube, cloudy pixels are removed using the dataset's cloud mask and the remaining valid pixels are averaged per timestamp; when an overpass is fully cloudy, no NDVI value exists and the timestamp is treated as missing rather than interpolated. The official split is used, with 2017–2019 for training and 2020 as the test set, giving 55,341 training samples and 3,336 test samples.

Forecasting samples are drawn with a sliding window over only the timestamps where clear-sky NDVI actually exists. The model sees the last three clear-sky NDVI values (p = 3) with their past covariates and must predict the next three clear-sky values (h = 3), a shift of four observations that corresponds to roughly 14 days depending on the revisit schedule. Calendar time is encoded with sine and cosine terms of the day of the year up to the third harmonic.

Because future weather is treated as forecast-available, the authors use observed daily meteorology as a proxy for forecasts but deliberately corrupt it during training with horizon-dependent multiplicative noise, g_k = 1 + β·Δt_k, scaled so that g_K = 2 at the final step, meaning perturbation intensity doubles across the window. The base perturbation is 10% of the variable magnitude, chosen by sensitivity analysis. All variables are standardized using training statistics and then passed through an arcsinh transform to limit outlier influence.

The network has two branches. The history branch encodes past NDVI and past covariates, masks missing targets, and compresses the output with masked temporal average pooling. The future branch encodes the whole covariate sequence with self-attention so intermediate timesteps inform the representation, then selects only the embeddings at real Sentinel-2 acquisition times. The pooled history vector is concatenated with those selected future embeddings and passed to a quantile head that outputs all future timesteps at once, without autoregressive decoding. Training uses the pinball loss at three quantile levels, weighted by temporal distance. Both branches use identical transformer encoders with 8 layers, d_model = 128, 8 attention heads, a 512-dimensional feed-forward network, and dropout of 0.1. Training uses PyTorch and Adam for 200 epochs with batch size 128 and an initial learning rate of 10^-4, a 20% validation split of the training data for model selection, and learning-rate reduction by a factor of 0.2 after 20 epochs without validation improvement down to a minimum of 5 × 10^-5. Experiments run on an NVIDIA T4 GPU.

Evaluation covers four point metrics on the median prediction (RMSE, MAE, WMAPE, MASE) and two probabilistic metrics (CRPS and averaged Pinball loss across the three quantile levels). The paper does not report calibration or coverage diagnostics, interval-width statistics, or inference latency.

Why This Matters

This work pushes vegetation forecasting away from the assumption of clean, regularly gridded satellite composites and toward the messy reality of cloud-masked, irregularly revisited observations. It shows that a relatively small transformer — 2.16 M parameters and 111.97 MFLOPs — can beat much larger LLM-based and transformer baselines such as TimeLLM (124.44 M parameters) and PatchTST, and can beat a zero-shot Chronos-2, on both pointwise and probabilistic metrics. The temporal-distance weighted loss is a generalizable idea for any sensor with irregular revisit patterns, and the finding that target history dominates while weather covariates add secondary gains gives practitioners a clear prioritization when data is scarce.

Real-world applications:

• Irrigation scheduling: A 14-day-ahead NDVI distribution, with the 10th and 90th percentiles bounding uncertainty, supports decisions about when and how much to water before stress becomes visible.
• Fertilization and nutrient management: Field-level greenness trajectories help time applications to crop phenological stages rather than to fixed calendars.
• Early stress detection and mitigation: Quantile forecasts flag fields where the lower bound of predicted greenness is deteriorating, allowing scouting or intervention before yield loss.
• Regional agricultural monitoring and insurance: Probabilistic forecasts across European ecozones can feed index-based insurance or advisory services where in-field sensors are impractical at scale.

Industry relevance: The model's small size and low computational cost make it a plausible candidate for operational deployment pipelines rather than research-only experiments, and the paper's public code release (https://github.com/arco-group/ndvi-forecasting) lowers the barrier to adaptation. The demonstrated robustness to imperfect weather inputs — only 2% RMSE and 1.6% CRPS degradation at 20% perturbation — matters for real deployments where the weather forecast is never perfect. The explicit caveat that evaluation is restricted to European agricultural areas is important for anyone considering transfer to other regions, where phenology, management practices, and climate regimes may differ substantially.

Future Directions

• Explicit conditioning on climate zone and crop type. The authors name this directly as future work, motivated by the observed performance decline from semi-arid (MAE 0.042, R² 0.791) to continental regimes (MAE 0.073, R² 0.557).
• Generalization beyond Europe. Because the study only covers European agricultural areas, it remains open whether the framework transfers to regions with different crop phenology, management practices, and climate regimes.
• Replacing observed weather with real forecast products. The current setup uses observed daily meteorology as a proxy for weather forecasts, perturbed with a synthetic noise model; testing against actual numerical weather prediction output would quantify real-world degradation and might change the optimal base perturbation level, which was set at 10% via sensitivity analysis.
• Extending the horizon and evaluating longer-range uncertainty. The current configuration predicts three clear-sky acquisitions (roughly 14 days), and the paper provides no evidence about how the temporal weighting and quantile calibration behave at substantially longer horizons where uncertainty accumulates.

Target Audience

This paper is most useful to machine learning researchers working on time-series forecasting under irregular sampling and on multimodal fusion of satellite and meteorological data. It is also relevant to remote sensing and Earth observation scientists interested in operational vegetation monitoring, and to precision agriculture practitioners or agtech engineers who need short-term, uncertainty-aware field-level forecasts and can work with NDVI trajectories rather than dense image products. Readers without background in transformer architectures, quantile regression, or probabilistic forecast evaluation will find the methods section dense, though the problem framing and the results tables are accessible.

Authors’ abstract

Short-term forecasting of vegetation dynamics is a key enabler for data-driven decision support in precision agriculture. Normalized Difference Vegetation Index (NDVI) forecasting from satellite observations, however, remains challenging due to sparse and irregular sampling caused by cloud masking, as well as the heterogeneous climatic conditions under which crops evolve. In this work, we propose a probabilistic forecasting framework for field-level NDVI prediction under sparse, irregular clear-sky acquisitions. The architecture separates the encoding of historical NDVI and meteorological observations from future exogenous covariates, fusing both representations for multi-step quantile prediction. To address irregular revisit patterns and horizon-dependent uncertainty, we introduce a temporal-distance weighted quantile loss that aligns the training objective with the effective forecasting horizon. In addition, we incorporate cumulative and extreme-weather feature engineering to capture delayed meteorological effects relevant to vegetation response. Experiments on European satellite data show that the proposed approach outperforms statistical, deep learning, and time-series baselines on both pointwise and probabilistic evaluation metrics. Ablation studies confirm that target history is the primary driver of performance, with meteorological covariates providing additional gains in the full multimodal setting. The code is available at https://github.com/arco-group/ndvi-forecasting.