Scaling spatiotemporal transformers for regional food security: a tri-stream latent-dynamic approach to pre-harvest yield and price forecasting

1 Introduction

Food insecurity affects over 820 million people globally, with the most severe impacts concentrated in regions experiencing the compound effects of climate variability, conflict, and economic instability (FAO, 2023). The Lake Chad Basin, spanning Cameroon, Chad, Nigeria, and Niger, exemplifies these challenges: the region has experienced a 90% reduction in lake surface area since the 1960s, persistent armed conflict, and recurrent drought conditions that threaten the livelihoods of over 30 million people dependent on rain-fed agriculture (Mahmood et al., 2019; Nagarajan et al., 2018). In such contexts, accurate pre-harvest yield forecasting - ideally 90 days before harvest - can enable timely interventions including food aid prepositioning, market stabilization measures, and adaptive agricultural support (Brown and Funk, 2008).

Traditional yield forecasting approaches rely on crop simulation models calibrated with ground-truth data, but these methods face significant limitations in data-sparse regions where meteorological stations are scarce, soil surveys incomplete, and administrative yield statistics often delayed by months or years (Jones et al., 2003; van Ittersum et al., 2013). Satellite remote sensing offers a promising alternative, providing consistent, near-real-time observations of vegetation conditions across large spatial extents (Atzberger, 2013). The Sentinel-1 and Sentinel-2 missions, in particular, deliver complementary synthetic aperture radar (SAR) and multispectral optical imagery at 10-m resolution with revisit times of 5–12 days, enabling detailed monitoring of crop phenology even in cloud-prone tropical regions (Drusch et al., 2012; Torres et al., 2012).

Recent advances in deep learning have demonstrated substantial improvements over traditional statistical approaches for yield forecasting (Kamilaris and Prenafeta-Boldú, 2018; van Klompenburg et al., 2020; Chlingaryan et al., 2018). Convolutional neural networks (CNNs) effectively extract spatial patterns from satellite imagery (Kussul et al., 2017), while recurrent architectures including Long Short-Term Memory (LSTM) networks capture temporal dynamics in vegetation index time series (Sun et al., 2019; Hochreiter and Schmidhuber, 1997). The Transformer architecture, originally developed for natural language processing (Va et al., 2017), has shown particular promise for time series forecasting due to its ability to model long-range dependencies through self-attention mechanisms (Zhou et al., 2021; Wu et al., 2021). However, existing satellite-based yield prediction systems typically treat remote sensing observations in isolation, failing to integrate the economic and climatic context that substantially influences agricultural outcomes in smallholder-dominated systems.

The relationship between predicted yields and actual food security outcomes is mediated by complex market dynamics, particularly in regions where infrastructure limitations and conflict create fragmented markets with high price volatility (Minot, 2014; Bellemare, 2015). The Soudure period - the pre-harvest lean season when food stores are depleted but new crops are not yet available - represents a critical vulnerability window during which food prices can spike dramatically, particularly following poor harvests (Devereux et al., 2014). Predicting not only yields but also price stability conditions therefore has direct operational value for humanitarian response planning.

This paper introduces the Tri-Stream Latent-Dynamic Transformer (TLDT), a novel architecture designed to address these challenges through three key innovations. First, we propose a tri-stream encoder architecture that processes spatial-spectral features (from Sentinel-1/2 imagery), temporal phenological dynamics (monthly vegetation index sequences), and exogenous signals (rainfall anomalies, price volatility, Soudure stress indicators) through specialized encoder modules. Second, we introduce a Tri-Stream Cross-Attention (TSCA) mechanism that employs spatial features as queries, economic context as keys, and temporal phenology as values, enabling the model to learn how market conditions modulate the relationship between observed crop growth and predicted yields. Third, we develop an Adaptive Positional Encoding (APE) scheme that handles the irregular temporal sampling characteristic of cloud-affected optical satellite data by learning continuous positional embeddings from observation time deltas.

We evaluate TLDT on a dataset covering the Lake Chad Basin, integrating satellite imagery from Google Earth Engine, historical yield statistics from FAOSTAT, rainfall data from CHIRPS, and market prices from the World Food Programme. Because yield is recorded only at national granularity, the resulting 288 location-year-crop samples reduce to just 18 independent (country, year) units - far fewer than would typically be needed to train and evaluate a multi-head Transformer with confidence. Under a leave-one-year-out cross-validation protocol appropriate to this sample size, TLDT performs comparably to rather than clearly better than a standard Transformer baseline on yield prediction (RMSE = 0.360 ± 0.103 vs. 0.350 ± 0.057 t/ha; mean R ² = −0.10), and a leakage diagnostic shows that its price-stability classification accuracy is statistically indistinguishable from a trivial logistic regression trained on three input features alone. We report these results in full, together with the data-pipeline corrections that produced them, as this paper’s central contribution: a transparent, cautionary case study of what a tri-stream fusion architecture does and does not demonstrably add in a genuinely data-sparse agricultural forecasting setting, rather than a claim of architectural superiority the data cannot support. The price-stability task is presented as exploratory only, given the leakage diagnostic above, and the Harvest-Expectation Index is presented as a conceptual illustration of how such outputs could eventually be operationalized, not as a validated early-warning tool.

The remainder of this paper is organized as follows. Section 2 reviews related work on yield forecasting, satellite remote sensing, and deep learning approaches. Section 3 details the TLDT architecture and training methodology. Section 4 describes our study area and data sources. Section 5 presents experimental results and analysis. Section 6 discusses implications and limitations. Section 7 concludes with directions for future work.

2 Related work

2.1 Satellite-based crop yield estimation

Satellite remote sensing has been employed for crop monitoring and yield estimation since the 1970s, with early approaches focusing on empirical relationships between vegetation indices and crop productivity (Tucker, 1979; Rouse et al., 1974). The Normalized Difference Vegetation Index (NDVI), computed from red and near-infrared reflectance, remains widely used due to its correlation with leaf area index and photosynthetic capacity (Pettorelli et al., 2005). More sophisticated indices including the Enhanced Vegetation Index (EVI), which incorporates blue band reflectance to reduce atmospheric effects, and the Normalized Difference Water Index (NDWI), sensitive to crop water content, provide complementary information about crop condition (Huete et al., 2002; Gao, 1996).

The launch of the Sentinel missions has substantially enhanced capabilities for agricultural monitoring. Sentinel-2’s Multispectral Instrument (MSI) provides 13 spectral bands at 10–60 m resolution with a 5-day revisit time at the equator, enabling detailed phenological monitoring (Defourny et al., 2019). Sentinel-1’s C-band SAR provides all-weather, day-night imaging capability particularly valuable in cloud-prone tropical regions, with backscatter signatures sensitive to crop structure and water content (Mandal et al., 2020; Veloso et al., 2017). Studies have demonstrated the complementary value of combining optical and SAR observations for crop classification and yield estimation (Orynbaikyzy et al., 2019; Steele-Dunne et al., 2017).

Regional yield forecasting systems typically aggregate pixel-level observations to administrative units for comparison with official statistics. Approaches range from simple regression models relating peak-season vegetation indices to reported yields (Bolton and Friedl, 2013) to more sophisticated time series analysis methods that exploit the full temporal profile of crop development (Cai et al., 2019; Sakamoto et al., 2013). However, the coarse temporal resolution of many satellite products and the irregular sampling imposed by cloud contamination in optical imagery remain significant challenges, particularly for identifying critical phenological stages (Whitcraft et al., 2015).

2.2 Deep learning for agricultural forecasting

Deep learning methods have achieved substantial improvements over traditional approaches for crop yield prediction. Convolutional neural networks effectively extract spatial features from satellite imagery, learning hierarchical representations that capture field-level patterns (Zhong et al., 2019; Ji et al., 2018). You et al. (You et al., 2017) demonstrated that CNNs applied to MODIS imagery could predict county-level corn yields in the United States with lower error than traditional regression approaches. Subsequent work has explored deeper architectures including ResNet and DenseNet variants (Nevavuori et al., 2019; Khaki and Wang, 2019).

Recurrent neural networks, particularly LSTM architectures, have proven effective for modeling the temporal dynamics of crop growth (Jiang et al., 2020; Luo et al., 2022). Satellite time series naturally represent sequential observations, and LSTM’s gating mechanisms enable learning of both short-term and long-term dependencies in vegetation index trajectories (Wolanin et al., 2020). Hybrid CNN-LSTM architectures combine spatial and temporal processing, with CNNs extracting features from individual images and LSTMs modeling their evolution over the growing season (Tian et al., 2021; Khaki et al., 2020).

The Transformer architecture has recently emerged as a powerful alternative to recurrent models for sequence modeling (Va et al., 2017). Self-attention mechanisms enable direct modeling of dependencies between any pair of time steps, avoiding the sequential processing bottleneck of RNNs. For time series forecasting, Transformers have demonstrated strong performance across diverse domains including energy demand, traffic flow, and weather prediction (Zhou et al., 2021; Lim et al., 2021; Wu et al., 2020). Applications to satellite time series analysis include crop classification (Rußwurm and Körner, 2020), land cover mapping (Garnot et al., 2020), and phenology extraction (Yuan and Lin, 2021).

However, most existing approaches treat satellite observations as the primary or sole input, with limited integration of the broader context that influences agricultural outcomes. Some studies have incorporated weather data as auxiliary inputs (Peng et al., 2020; Lobell et al., 2015), but systematic integration of economic indicators, market dynamics, and vulnerability factors remains largely unexplored.

2.3 Food security early warning systems

Operational food security monitoring systems, including FEWS NET, the Global Information and Early Warning System (GIEWS), and the Integrated Food Security Phase Classification (IPC), combine multiple data sources to assess food security conditions and project future trends (Funk et al., 2010; Global Partners, 2021). These systems typically integrate satellite-derived indicators of crop and rangeland conditions with market price monitoring, conflict tracking, and household vulnerability assessments (Brown et al., 2014).

Market price dynamics play a crucial role in food security outcomes, particularly in regions where a large share of households are net food purchasers (Headey and Fan, 2008). Price volatility can undermine food access even when aggregate production is adequate, and price spikes frequently coincide with the Soudure period when household food stocks are depleted (Devereux et al., 2014). The relationship between production shortfalls and price increases is mediated by market integration, trade policies, and speculative behavior, creating complex dynamics that are difficult to predict from production indicators alone (Gilbert and Morgan, 2010; Tadesse et al., 2014).

Machine learning approaches have been applied to food security prediction, including classification of IPC phases from household survey data (Lentz et al., 2019) and prediction of acute food insecurity from satellite and conflict indicators (Martini et al., 2022). However, these approaches typically treat yield and price outcomes as separate prediction targets rather than modeling their joint dynamics. Our work addresses this gap through multi-task learning that simultaneously predicts yield quantities and price stability classifications.

2.4 Multi-modal fusion and cross-attention

Integrating information from heterogeneous data sources is a fundamental challenge in machine learning, with approaches ranging from early fusion (concatenating inputs) to late fusion (combining model outputs) to more sophisticated attention-based mechanisms (Baltrušaitis et al., 2019). Cross-attention, which computes attention weights between queries from one modality and keys/values from another, has proven particularly effective for multi-modal learning (Lu et al., 2019; Tsai et al., 2019).

In remote sensing, fusion of optical and SAR imagery has been extensively studied, with deep learning approaches including shared-weight networks, adversarial training, and attention-based selection mechanisms (Hong et al., 2020; Schmitt and Zhu, 2016). Temporal fusion strategies for satellite time series have explored attention mechanisms that weight observations by quality, cloud cover, or phenological relevance (Nyborg et al., 2022; Interdonato et al., 2019).

Our Tri-Stream Cross-Attention mechanism extends these approaches by explicitly modeling the interaction between spatial observations, temporal dynamics, and economic context. By employing spatial features as queries, economic context as keys, and temporal phenology as values, we enable the model to learn attention patterns that reflect how market conditions modulate the relationship between observed vegetation conditions and agricultural outcomes.

3 Methods

3.1 Problem formulation

We formulate pre-harvest yield and price prediction as a multi-task learning problem. Given observations up to time t (corresponding to approximately 90 days before harvest), we aim to simultaneously predict the final harvest yield y∈R⁺ (tons per hectare) and classify price stability conditions c∈{0,1,2} representing Stable, Vulnerable, and Crisis states respectively.

For each location-year sample, we define three input streams. The spatial-spectral stream Xs ∈ R^C×H×W represents multi-band satellite imagery with C = 10 channels (including Sentinel-2 bands B2, B3, B4, B8, B11, B12 and derived indices NDVI, EVI, SAVI, NDWI), with spatial dimensions H=W=32. The temporal stream X_t​ ∈R^T×D represents monthly observations over T = 12 months with D = 8 spectral-temporal features (NDVI, EVI, SAVI, NDWI from Sentinel-2; VV, VH, VH/VV ratio from Sentinel-1; precipitation). The exogenous stream X_e​ ∈ R³ comprises three scalar features: rainfall anomaly during the growing season, price volatility coefficient, and Soudure stress indicator.

The model learns a mapping f: (X_s​,X_t​,X_e​)→(y,c) parameterized by neural network weights θ, optimized to minimize a combined loss function balancing yield prediction accuracy and price classification performance.

3.2 Overall architecture

The TLDT architecture, illustrated in Figure 1, comprises four main components: stream-specific encoders for spatial, temporal, and exogenous inputs; the Tri-Stream Cross-Attention fusion module; a feature integration layer; and task-specific prediction heads for yield regression and price classification.

3.3 Spatial-spectral encoder

The Spatial-Spectral Encoder processes satellite imagery through a convolutional neural network backbone designed to extract field-level spatial patterns and spectral signatures indicative of crop type and condition. The architecture comprises three convolutional blocks with progressively increasing channel depth, followed by adaptive average pooling to produce a fixed-length feature vector regardless of input spatial dimensions.

Each convolutional block contains two 3 × 3 convolutional layers with batch normalization and ReLU activation, followed by 2 × 2 max pooling:

h s l = M a x P o o l R e L U B N C o n v 3 x 3 R e L U B N C o n v 3 x 3 h s l − 1

The encoder maps input X_s ∈ ℝ^10×32×32 through channel dimensions [10 → 64 → d], where d is the hidden dimension (128 in our experiments), with final adaptive average pooling producing h_s ∈ ℝ^d.

This architecture captures multi-scale spatial patterns from field boundaries and vegetation texture at fine scales to landscape-level heterogeneity at coarser scales. The choice of relatively shallow architecture reflects the limited spatial extent of our input patches and the need to avoid overfitting given our modest dataset size.

3.4 Temporal phenology encoder with adaptive positional encoding

The Temporal Phenology Encoder processes monthly observation sequences through a Transformer encoder architecture, capturing the temporal dynamics of crop development from planting through maturity. A key innovation is our Adaptive Positional Encoding (APE) scheme, designed to handle the irregular temporal sampling that arises from cloud contamination in optical satellite imagery.

3.4.1 Adaptive positional encoding

Standard positional encodings for Transformers assume regular sampling intervals, using either sinusoidal functions or learned embeddings indexed by position (Va et al., 2017). This assumption is violated in satellite time series, where cloud cover can result in gaps of days to weeks between usable observations. We propose APE, which learns positional embeddings directly from observation time deltas:

A P E δ t = tanh W 2 . G E L U W 1 . δ t 60

Where δ_t is the time delta (in days) since the previous observation, normalized by the maximum expected gap (60 days), and W₁ ∈ ℝ^1×d/2, W₂ ∈ ℝ^d/2×d are learned projection matrices. The GELU activation provides smooth gradients, while the tanh output bounds embeddings to facilitate stable training.

APE enables the model to learn that, for example, a 30-day gap during the dry season has different phenological implications than a 30-day gap during rapid vegetative growth. The continuous formulation also gracefully handles varying gap lengths without requiring discretization.

3.4.2 Transformer encoder

The temporal sequence is first projected to the hidden dimension via a linear layer, then combined with APE embeddings:

Z 0 = X t W p r o j + A P E δ

Where W p r o j ∈ R D × d. The sequence is then processed through L transformer encoder layers, each comprising multi-head self-attention and position-wise feed-forward networks:

Z l ′ = L a y e r N o r m Z l − 1 + M H A Z l − 1 , Z l − 1 , Z l − 1

Z l = L a y e r N o r m Z l ′ + F F N Z l ′

Where MHA denotes multi-head attention and FFN is a two-layer feed-forward network with GELU.

The final temporal representation is obtained by mean pooling over the sequence dimension followed by layer normalization:

h t = L a y e r N o r m 1 T ∑ t − 1 T Z L t

Where L = 2 layers and 4 attention heads in our experiments, balancing model capacity against overfitting risk.

3.5 Exogenous encoder

The Exogenous Encoder processes three scalar features - rainfall anomaly, price volatility, and Soudure stress - through a multi-layer perceptron:

h c = L a y e r N o r m R E L U W 2 . R E L U W 1 . X e

Where W₁ ∈ ℝ^3×d/2 and W₂ ∈ ℝ^d/2×d. Dropout is applied after each ReLU activation for regularization. The rainfall anomaly component is computed as:

a r a i n = P ¯ s e a s o n − μ P σ P

Where P_season is observed seasonal precipitation and μ_P, σ_P are climatological mean and standard deviation. Price volatility is quantified as the coefficient of variation of monthly cereal prices, and the Soudure stress indicator is computed as:

s s o u d u r e = max ⁡ 0 , − a r a i n , I A S

Taking a positive value only when rainfall is below normal during this critical window.

The input exogenous vector X_e is computed using only information available up to the 90-day pre-harvest forecast cutoff (September), while the price-stability label used as the classification target (Section 4.6) is computed from the full realized marketing season, including the October–December post-cutoff period. Price volatility is built from the nearest reporting market to each sampling location, identified by great-circle distance within the WFP/HDX HAPI price records, rather than from a price series pooled across an entire country (Section 4.2.3; Section 4.5). This temporal and spatial separation between the input window and the label construction window creates a genuine forecasting setup rather than a same-window restatement, though it does not eliminate all overlap between input and label; Section 5.3 quantifies the remaining overlap directly via a leakage diagnostic.

3.6 Tri-stream cross-attention fusion

The core architectural innovation of TLDT is the Tri-Stream Cross-Attention (TSCA) module, which integrates information across the three encoded streams through a purpose-designed attention mechanism. Unlike standard cross-attention that operates between two sequences, TSCA employs a three-way configuration:where queries are derived from spatial features, keys from economic context, and values from temporal phenology. This design reflects the intuition that we wish to retrieve relevant temporal patterns (phenological dynamics) based on how spatial observations (current crop conditions) relate to economic context (market conditions and climate stress).

Q = h s W Q , K = h e W K , V = h t W V

Multi-head attention is computed as:

h e a d i = s o f t m a x Q i K i T d k V i

h f = C o n c a t h e a d 1 , … , h e a d H W O

Where H is the number of attention heads and d k = d H is the per-head dimension.

The attention weights A = s o f t m a x Q K T / d k provide interpretable information about how strongly economic context modulates the fusion process for each sample. High attention to economic context indicates samples where market conditions significantly influence the mapping from observed crop conditions to predicted yields - for example, during conflict periods when market disruptions may cause price spikes even with adequate production.

3.7 Feature fusion and prediction heads

The four representations - spatial, temporal, exogenous, and cross-attended - are concatenated and processed through a fusion network:where W f u s e ∈ R 4 d × d. This produces the unified representation used by both prediction heads. The yield prediction head is a two-layer MLP with ReLU activation:

h = L a y e r N o r m R e L U W f u s e h s ; h t ; h e ; h f

y ^ = W y T R e L U W y , 1 h

3.8 Multi-task learning objective

We optimize a combined loss function balancing yield prediction and price classification:

L = λ y L y i e l d + λ c L p r i c e

Where the yield loss is mean squared error:

L y i e l d = 1 N ∑ i − 1 N y ^ i − y i 2

And the price loss is cross-entropy:

L p r i c e = − 1 N ∑ i = 1 N ∑ k = 0 2 1 c i = k log exp ⁡ c ^ i , k ∑ j = 0 2 ⁡ exp ⁡ c ^ i , j

Where λ y = 1.0 a n d λ c = 0.5 based on validation performance, giving higher weight to the primary yield prediction task.

3.9 Training procedure

Models are trained using the AdamW optimizer (Loshchilov and Hutter, 2019) with learning rate 10^–3, weight decay 10^–4, and batch size 32. We employ a reduce-on-plateau learning rate scheduler, halving the learning rate when validation loss plateaus for 5 epochs. Training proceeds for a maximum of 60 epochs with early stopping based on validation loss with patience of 12 epochs.

Gradient clipping with maximum norm 1.0 is applied to prevent exploding gradients. Dropout with probability 0.1 is applied throughout the network for regularization. All models are implemented in PyTorch and trained on a single CPU (experiments can be accelerated with GPU when available).

Data augmentation is limited due to the nature of the task: we apply small Gaussian noise to exogenous features during training but do not augment spatial or temporal features to preserve the integrity of phenological signals.

4 Study area and data

4.1 Study area: Lake Chad basin

Our study focuses on the Lake Chad Basin, a transboundary region encompassing portions of Cameroon, Chad, Nigeria, and Niger. The basin is characterized by a semi-arid climate with a single rainy season from June to October, supporting rain-fed agriculture dominated by sorghum and millet production (Seignobos and Iyébi-Mandjek, 2000). These staple cereals are critical for food security, providing the primary source of calories for the majority of the basin’s population. Figure 2 below shows the spatial distribution of the 24 sampling locations across the three countries, and Figure 3 shows Climatological and Vegetation Profiles of part of the Study locations.

The region faces multiple overlapping challenges including climate variability with frequent droughts, ongoing conflict involving armed groups, population displacement, and limited infrastructure constraining market access (OCHA, 2023). Average sorghum yields in the region range from approximately 0.6–1.2 tonnes per hectare, substantially below potential yields due to limited input use and climatic constraints (Aune et al., 2007).

We selected 24 sampling locations distributed across the three countries, focusing on major agricultural zones and market centers. In Cameroon, locations include Maroua (the regional capital), Garoua, Mokolo, Kaele, Mora, Kousseri, Yagoua, and Guider. In Chad, we include N'Djamena, Moundou, Bongor, Abeche, Sarh, Mongo, Ati, and Pala. In Nigeria, locations comprise Maiduguri, Yola, Bauchi, Gombe, Potiskum, Biu, Damaturu, and Dikwa.

4.2 Satellite data from Google Earth Engine

4.2.1 Sentinel-2 optical imagery

We extracted Sentinel-2 Level-2A surface reflectance imagery from Google Earth Engine for the period 2018–2023 (Gorelick et al., 2017). For each sampling location, we defined a 5 km buffer and computed monthly composites using median aggregation of cloud-free observations. Cloud masking was performed using the Scene Classification Layer (SCL), retaining pixels classified as vegetation, bare soil, or water.

From each monthly composite, we extracted six spectral bands (B2 Blue, B3 Green, B4 Red, B8 NIR, B11 SWIR1, B12 SWIR2) and computed four vegetation indices. The Normalized Difference Vegetation Index (NDVI) was calculated as:

N D V I = ρ N I R − ρ R e d ρ N I R + ρ R e d

The Enhanced Vegetation Index (EVI) incorporates the blue band

E V I = 2.5 × ρ N I R − ρ R e d ρ N I R + ρ R e d − 7.5 ρ B l u e + 1

The soil Adjusted Vegetation Index (SAVI) accounts for soil background effects:

S A V I = ρ N I R − ρ R e d ρ N I R + ρ R e d + 0.5 × 1.5

The Normalized Difference Water Index (NDWI) is sensitive to vegetation water content:

N D W I = ρ N I R − ρ S W I R 1 ρ N I R + ρ S W I R 1

The extraction yielded 1,472 monthly records across all locations, with an average of 61 cloud-free months per location over the 72-month study period (85% coverage).

4.2.2 Sentinel-1 SAR imagery

Sentinel-1 Ground Range Detected (GRD) imagery in Interferometric Wide swath mode was extracted for VV and VH polarizations (Mandal et al., 2020). We filtered for ascending orbit passes and computed monthly median composites within the same 5 km buffers used for Sentinel-2.

SAR backscatter provides complementary information to optical imagery, remaining usable under cloudy conditions and providing sensitivity to vegetation structure and water content (Vreugdenhil et al., 2018). We computed the VH/VV ratio (in decibels), which has been shown to correlate with vegetation biomass and crop growth stage (Khabbazan et al., 2019).

The extraction yielded 1,440 monthly records. Coverage was higher in the western portions of the study area; some locations in eastern Chad (Abeche, Sarh, Mongo, Ati) had limited Sentinel-1 coverage due to orbital characteristics.

4.2.3 Geospatial processing workflow

This subsection documents the geospatial preprocessing pipeline underlying Sections 4.2 and 4.3 precisely as implemented, so that the extraction procedure is fully reproducible. Sampling-location coordinates are stored as point geometries in WGS84 (EPSG:4326), the native CRS of both the Sentinel-1/2 and CHIRPS products used in this study. Around each point we construct a 5 km buffer for Sentinel-1/2 extraction and a 25 km buffer for CHIRPS, reflecting its coarser native pixel size (≈5.5 km at the equator versus 10 m for Sentinel-2). All extraction is executed natively within Google Earth Engine at each product’s native resolution and CRS, with no manual reprojection or resampling performed outside of Earth Engine’s own reducer operations. Monthly composites are built by per-pixel median aggregation of all cloud-free (Sentinel-2) or available (Sentinel-1, CHIRPS) observations within the calendar month, prior to a buffer-mean spatial reduction that yields a single scalar value per band, per month, per location. Missing months arising from residual cloud contamination (Section 4.2.1) are filled by linear interpolation across the temporal sequence; irregular gaps that remain are additionally represented explicitly to the model via the observation time-delta inputs to the Adaptive Positional Encoding (Section 3.4.1), rather than being treated as fully resolved by the interpolation step alone. Cereal market prices are assigned to each sampling location by nearest-market spatial join, described in Section 4.5.

4.3 CHIRPS rainfall data

We extracted monthly precipitation totals from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) dataset, available at 0.05° resolution (Funk et al., 2015). For each location, we computed the monthly mean precipitation within a 25 km buffer and derived rainfall anomalies relative to the 1981–2010 climatological baseline.

The growing season (June–October) rainfall anomaly serves as the primary exogenous climate indicator. Negative anomalies during this period indicate drought stress likely to reduce yields, while positive anomalies may indicate favorable conditions or, in extreme cases, flooding risk.

The Soudure stress indicator specifically targets the July–September period when crop water demands are highest. Rainfall deficits during this window have disproportionate impacts on final yields due to sensitivity during flowering and grain filling stages (Sultan et al., 2013).

The extraction yielded 1,728 monthly records with complete coverage across all locations and time periods.

4.4 FAOSTAT yield statistics

Historical crop yield data were obtained from the Food and Agriculture Organization’s Statistical Database (FAOSTAT), Crops and Livestock Products domain (QCL) (FAO, 2024). Sorghum and millet yield statistics for Cameroon, Chad, and Nigeria (2000–2023) were extracted from the FAOSTAT QCL bulk data file (Production_Crops_Livestock_E_All_Data_(Normalized), the same underlying file used by FAOSTAT’s own R package and by reproducible-research pipelines such as Our World in Data’s FAOSTAT ETL), filtered locally to the target countries, items (Sorghum, Millet), and years (2000–2023), yielding 432 FAOSTAT records, cross-referenced in kilograms per hectare (kg/ha).

FAOSTAT provides country-level statistics only; sub-national yield data at the location level are not available for these countries, and every sampling location within a country therefore receives the identical national yield value for a given crop-year. Because sub-national yield statistics are unavailable at location-level granularity within the target area, each micro-sampling node uniformly inherits the corresponding national country-year yield baseline from the FAOSTAT registry. This constraint is treated as a deterministic property of the dataset’s information content, and model evaluation protocols are adjusted accordingly to prevent artificial performance inflation.

Across the corrected, integrated dataset used for model training and evaluation (sorghum and millet combined, 2018–2023), mean observed yields were 1.48 t/ha in Cameroon, 0.71 t/ha in Chad, and 1.10 t/ha in Nigeria, consistent in order of magnitude with published regional statistics (Niang et al., 2014) and with the challenging growing conditions of the Sahel.

4.5 WFP market price data

Cereal market prices were obtained from the World Food Programme’s food-price data, accessed via the Humanitarian Data Exchange (HDX) (WFP, 2024), using the HDX HAPI (Humanitarian API) schema food-price files for Cameroon, Chad, and Nigeria, yielding 14,875 sorghum and millet price records covering 2018–2023.

To establish an explicit spatial relationship between price data and sampling locations, individual market-level coordinates within the HDX HAPI metadata are cross-referenced with each sampling location’s coordinates, and each location is assigned its nearest reporting market by great-circle distance; the price-volatility exogenous feature is constructed from that location-specific market series rather than a pooled country-level series (Section 4.2.3). We note explicitly, and revisit in Section 6.4, that nearest-market assignment is a simplification: it does not capture trade-flow or market-integration effects that can decouple prices between nearby but poorly-connected markets, particularly in conflict-affected areas.

4.6 Dataset construction

The final integrated dataset comprises 288 samples, each representing one location-year-crop combination. For each sample, we constructed:

• Spatial features: 10-channel × 32 × 32 spatial arrays representing growing season composite statistics, synthesized from monthly band means and standard deviations using smooth spatial patterns to approximate field-level heterogeneity.
• Temporal features: 12-month × 8-feature arrays containing monthly values of NDVI, EVI, SAVI, NDWI, VV, VH, VH/VV ratio, and precipitation. Missing values due to cloud cover were filled using linear interpolation.
• Time deltas: 12-element vectors indicating the number of days since the previous monthly observation, with 30-day defaults for consecutive months and 45-day values for months with heavy cloud contamination.
• Exogenous features: 3-element vectors containing rainfall anomaly, price volatility, and Soudure stress.
• Yield targets: Scalar values in tonnes per hectare, derived from FAOSTAT country statistics (Section 4.4). Because yield is recorded at national granularity, all locations within a country share an identical yield value for a given crop-year; see Section 4.7 for how this is handled in evaluation.
• Price labels: Three-class labels (0 = Stable, 1 = Vulnerable, 2 = Crisis) assigned based on a composite risk score incorporating rainfall anomaly, price volatility, and Soudure stress computed over the full realized marketing season (Section 3.5), with yield-based adjustments. Because this composite score is a deterministic function of variables closely related to the model’s own input features, Section 5.3 reports a leakage diagnostic quantifying how much of the label is predictable from the input features alone.

Rather than a single stratified random split, the dataset is evaluated using a group-aware, leave-one-year-out cross-validation protocol, described in Section 4.7.

The 288 location-year-crop samples do not represent 288 independent observations of yield: because FAOSTAT yield is recorded once per country-year and shared across all locations in that country-year, the dataset reduces to only 18 independent (country, year) units (3 countries × 6 years) with respect to the yield-prediction target. Treating the 288 rows as independent, as a simple random train/validation/test split would implicitly do, allows a model to “solve” a test row by having seen a different location from the same country-year during training - a leakage path closed here by the group-aware split described below.

4.7 Group-aware splitting and leave-one-year-out cross-validation

Because only 18 independent (country, year) groups underlie the 288 location-year-crop rows (Section 4.6), a single random train/validation/test split tests on only a handful of these groups: too few, and too volatile, to support a headline comparative claim. To address spatial-temporal dependencies and avoid artificial data leakage from country-year sharing, a group-aware, leave-one-year-out (LOYO) cross-validation scheme was implemented as the sole, authoritative evaluation protocol for this manuscript (Section 5.2).

We first re-split the data at the (country, year) group level, ensuring no country-year’s rows are ever divided across training and evaluation sets within a fold. LOYO cross-validation then uses each of the 6 years (2018–2023) in turn as a held-out test fold (48 rows/fold: 3 countries × 2 crops × 8 locations), training TLDT and both baselines from scratch on the remaining 5 years for each fold. Within each fold, a second year is rotated from the remaining training years to serve as the early-stopping validation set, so that the test year itself is never used for any training decision. This design directly tests whether the model generalizes to a year it has not seen - the property an operational pre-harvest forecasting system actually needs - and uses every available group as test material at least once, rather than only 3 of 18 groups. The leakage diagnostic reported in Section 5.3 (a logistic regression trained on the input-time exogenous features X_e alone) is evaluated on the identical per-fold splits, so that it is directly comparable to TLDT’s own cross-validated performance.

5 Experiments and results

5.1 Experimental setup

All experiments were conducted using PyTorch 1.12 on the google Colab environment. Models were trained with the hyperparameters described in Section 3.9. We compared TLDT against two baseline architectures:

LSTM Baseline: A hybrid CNN-LSTM architecture with spatial features processed by a simplified CNN (single convolutional layer followed by adaptive pooling), temporal features processed by a 2-layer LSTM with hidden dimension 128, and exogenous features processed by a linear layer. Features are concatenated and passed through a fusion MLP.

Transformer Baseline: Similar architecture to TLDT but without the Tri-Stream Cross-Attention module. Spatial, temporal, and exogenous encoders are identical to TLDT, with features concatenated directly for fusion without cross-modal attention.

Both baselines use the same multi-task loss function and training procedure as TLDT to ensure fair comparison.

5.2 Yield prediction performance

Table 1 presents leave-one-year-out cross-validated performance (mean ± standard deviation across the 6 folds; Section 4.7), which we treat as the authoritative quantitative result of this study. TLDT achieves RMSE = 0.360 ± 0.103 t/ha, MAE = 0.281 ± 0.079 t/ha, and R ² = −0.10 ± 0.56, performing comparably to - not clearly better than - the standard Transformer baseline (RMSE = 0.350 ± 0.057 t/ha), and modestly better than the LSTM baseline (RMSE = 0.390 ± 0.093 t/ha) on this metric. The wide fold-to-fold standard deviations, and the negative mean R ² for both TLDT and LSTM, indicate that no model in this comparison explains yield variance reliably better than predicting the mean; we report this candidly rather than as a limitation to minimize.

TABLE 1. Leave-one-year-out cross-validated performance across models (mean ± SD, 6 folds).

5.3 Price classification performance and leakage diagnostic

Under leave-one-year-out cross-validation, TLDT achieves 56.6% ± 19.2% accuracy and 0.564 ± 0.200 weighted F1-score on three-class price stability classification (Table 1). Because the price-stability label is constructed from a composite formula that shares variables with the model’s own exogenous input features (Section 4.6), we quantify directly how much of this classification signal is attributable to the input features alone via a leakage diagnostic (Table 2): a logistic regression trained on the input-time exogenous features X_e alone (three scalars), evaluated on the identical per-fold splits used for TLDT. This trivial baseline achieves 58.3% ± 19.4% accuracy and 0.551 ± 0.211 weighted F1-score - statistically indistinguishable from TLDT’s own cross-validated performance.

TABLE 2. Leakage diagnostic: exogenous-features-only classifier vs. TLDT, per leave-one-year-out fold.

The price-stability label, as currently defined by a hand-crafted composite risk-score formula (Section 4.6), remains largely a function of a small number of its own contributing variables even after the forecast-cutoff separation described in Section 3.5. We recommend an independent ground-truth label source - for example, Cadre Harmonisé or IPC phase classifications (Global Partners, 2021) - as a necessary next step before this task’s classification accuracy is treated as a genuine forecasting result, and we revisit this point in Section 6.3 and Section 6.6.

5.4 Regional analysis

A statistically defensible approach requires pooling out-of-fold predictions across all six LOYO folds and disaggregating by country, since each fold’s held-out year spans all three countries simultaneously. On closer inspection, however, this pooled breakdown would not resolve the underlying small-sample problem - it would simply relocate it. Because yield is recorded once per country-year and shared across every location within that country-year (Section 4.6), pooling predictions across folds still leaves only 6 independent year-observations behind any single country’s regional metric, one-third of the 18 units already responsible for the wide fold-to-fold standard deviations in Table 1 (e.g., R ² = −0.10 ± 0.56). A country-level RMSE or R ² built from 6 independent points would carry even less statistical meaning than the headline result we already report as inconclusive, and presenting one would sit uneasily against this paper’s own central argument: that small-N metrics in this dataset should not be read as decisive.

A reliable regional breakdown would require sub-national yield data (e.g., agricultural surveys or insurance records) that varies within, not just across, countries; we identify this explicitly as the necessary precondition for future regional analysis in Section 6.4.

5.5 Training dynamics

Figure 4 presents training curves showing loss, RMSE, R ², and accuracy evolution over epochs for all three models on one illustrative training run; exact convergence behaviour (e.g., the epoch at which early stopping triggers) varies across the six leave-one-year-out folds and is not summarized by a single curve.

The LSTM baseline shows the fastest initial convergence but plateaus at higher loss values, suggesting limited model capacity. The Transformer baseline shows more volatile training dynamics, consistent with the known difficulty of training Transformers with limited data (Dosovitskiy et al., 2021).

5.6 Per-fold cross-validation detail

Figure 5 breaks the leave-one-year-out result down by individual fold, for both yield RMSE and price classification accuracy. The fold-to-fold variation is substantial for all three models: no model is uniformly best across years, and the ranking between TLDT, LSTM, and the Transformer baseline changes from fold to fold. This visualizes directly why the aggregate comparison in Table 1 is reported as “comparable performance” rather than a decisive win for any one architecture.

Dashed horizontal lines in Figure 5 mark each model’s mean across all six folds, matching the summary statistics reported in Table 1; no single fold is representative of overall performance on its own, which is why we report cross-validated means and standard deviations throughout this paper rather than a single held-out split.

5.7 Harvest-Expectation Index: A simulated, unvalidated conceptual illustration

Figure 6 illustrates the Harvest-Expectation Index (HEI) concept through a simulated comparison of good harvest, drought, and normal year trajectories. This simulation is a conceptual illustration of how TLDT outputs could be operationalized, not a plot of held-out model predictions.

Under good harvest conditions, HEI rises through the growing season as favorable vegetation development is observed, reaching peak values in August–September before declining slightly to harvest. Drought conditions show declining HEI from May onward as rainfall deficits manifest in reduced vegetation indices. The corresponding price stability index (panel b) shows inverse patterns: stable prices under good harvests, crisis conditions (high prices) under drought.

This visualization illustrates, purely conceptually, how TLDT outputs might eventually be packaged for early-warning use; it is built entirely from simulated trajectories and has not been generated from, or validated against, held-out model predictions or independent food-security outcome data. Given the modest and statistically inconclusive yield and price-classification results reported in Sections 5.2–5.3, we do not present the HEI as an operational tool. Realizing it as one would require the underlying yield and price predictions to be substantially more reliable than demonstrated here, together with independent validation against real early-warning outcomes - both left as future work (Section 6.7).

5.8 Comparison with published benchmarks

Direct comparison with published benchmarks is challenging due to differences in study regions, crops, and evaluation metrics. Our cross-validated RMSE of 0.360 t/ha represents approximately 33% of the mean observed yield (1.10 t/ha across the integrated dataset), and our R ² is close to zero (mean −0.10, consistent with a model that does not reliably outperform predicting the mean). Prior work on sorghum, millet, and other cereal yield prediction in data-sparse, smallholder-dominated African settings has reported comparable or similarly modest performance: Sultan et al. (Sultan et al., 2013) and Burke and Lobell (Burke and Lobell, 2017) both report substantial unexplained variance when relating remote-sensing or climate-based predictors to smallholder cereal yields in West and East Africa, and systematic evaluations of satellite-based yield mapping at national scale (Jin et al., 2019; Cao et al., 2021) similarly report wide error margins in comparable data-sparse regions. We view a modest but honestly evaluated result, obtained end-to-end on real data, as a more useful basis for future work, and note that a full quantitative benchmark against these studies would require matching crops, regions, and evaluation protocols beyond the scope of this research.

6 Discussion

6.1 Value of tri-stream fusion

Under leave-one-year-out cross-validation, TLDT performs comparably to - rather than clearly better than - a standard Transformer baseline without the TSCA module (RMSE 0.360 ± 0.103 vs. 0.350 ± 0.057 t/ha; Table 1), and modestly outperforms an LSTM baseline on RMSE. On price classification, TLDT’s cross-validated accuracy (56.6% ± 19.2%) is statistically indistinguishable from a trivial exogenous-features-only baseline (Section 5.3; Table 2). We do not interpret these results as clear evidence that cross-modal attention enables more effective information integration on this dataset; but that, once evaluated with a design robust to the data’s true information content (18 independent country-year units; Section 4.7), none of the three architectures compared here demonstrates a decisive advantage.

The design of TSCA - using spatial features as queries, economic context as keys, and temporal phenology as values - remains motivated by domain knowledge about agricultural systems: the same vegetation index value may indicate adequate development in a favorable year or stressed crops in a drought year, and using economic context to modulate how spatial observations attend to temporal patterns is, in principle, a reasonable mechanism for capturing this. We see this as an architecturally motivated contribution whose benefit is not yet demonstrated on the present, modestly sized dataset.

6.2 Adaptive positional encoding for irregular sampling

Cloud contamination is a pervasive challenge for optical satellite time series analysis, particularly in tropical regions where cumulus development occurs daily during the rainy season. Our APE scheme addresses this by learning continuous positional embeddings from observation time deltas rather than assuming regular sampling.

While we cannot isolate the contribution of APE from other architectural components in our experiments, the strong performance on temporal sequences with variable gap lengths (averaging 15% missing observations) suggests that APE provides effective handling of irregular sampling. Future ablation studies could compare APE against standard sinusoidal encodings and learned positional embeddings to quantify its contribution.

6.3 Multi-task learning benefits

The joint prediction of yields and price stability classes through multi-task learning retains practical appeal for operational deployment - a single model producing both outputs reduces deployment complexity relative to separate yield and price systems. However, under leave-one-year-out cross-validation, TLDT’s price classification accuracy (56.6% ± 19.2%) is statistically indistinguishable from a trivial logistic regression on the three exogenous input scalars alone (58.3% ± 19.4%; Section 5.3; Table 2), so this operational-simplicity benefit is not accompanied by evidence that the model learns price-relevant signal beyond what a simple baseline already captures. This indicates that, as currently constructed, the price-stability label is largely a restatement of a small number of its own contributing variables, rather than a target that requires the full multi-stream representation to predict.

We do not view this as invalidating the multi-task framing conceptually, but it does mean the current price-classification results should be treated as tentative pending an independent ground-truth label (Section 5.3).

6.4 Regional heterogeneity

As discussed in Section 5.4, a disaggregated, country-level regional breakdown is excluded from our primary performance reporting due to statistical constraints. Because historical yield vectors are reported at country-year intervals and shared across all sub-national sampling nodes, a pooled out-of-fold regional breakdown would rest on only six independent year-observations per country. Compiling metrics from such limited sample slices introduces high variance, making country-specific breakdowns inconclusive.

Several factors plausibly contribute to regional heterogeneity in prediction difficulty, independent of the specific numbers reported. First, the assignment of country-level yields to all locations within a country (Section 4.4) introduces identical targets across locations that likely differ in real yield. Second, conflict impacts including displacement, abandoned fields, and restricted market access affect yields in ways not captured by satellite observations. Third, market fragmentation in conflict-affected areas may decouple local prices from the nearest-market series now used as our price-volatility feature (Section 4.5) - nearest-market assignment by great-circle distance does not capture trade-flow or market-integration effects that can decouple prices between nearby but poorly-connected markets, particularly in conflict-affected areas, and this remains an explicit simplification of our approach.

Addressing these challenges would require sub-national yield data (potentially from agricultural surveys or insurance programs), explicit modeling of conflict exposure, and denser market price observations with trade-flow information. The modular TLDT architecture could accommodate additional input streams capturing these factors.

6.5 Operational considerations

For operational deployment, several considerations would need to be addressed. First, the 90-day forecast horizon assumes that growing season observations up to September are available for October–December harvest prediction. In practice, processing and quality assurance of satellite data introduce latency that may reduce effective lead time.

Second, the model is trained on 2018–2023 data and would need regular retraining or adaptation as conditions evolve. Climate trends, changing agricultural practices, and conflict dynamics all affect the yield-price relationships that the model has learned.

Third, the country-level yield assignments in our training data limit the spatial specificity of predictions. For operational use, integration with sub-national survey data or farmer-reported yields would enable more localized forecasts.

Fourth, uncertainty quantification is not addressed in the current implementation. Bayesian approaches or ensemble methods could provide prediction intervals essential for risk-aware decision making.

6.6 Limitations

Several limitations of this study should be acknowledged. The sample size (288 rows, only 18 independent country-year units with respect to the yield target; Section 4.6) is small relative to typical deep learning applications and relative to what leave-one-year-out cross-validation with only 6 folds can estimate precisely, as reflected in the wide standard deviations in Table 1. The reliance on country-level yield statistics precludes evaluation of genuine sub-national prediction accuracy. The short study period (2018–2023) limits representation of extreme events including major droughts or exceptionally good harvests.

The spatial feature construction approach - synthesizing patch-level features from point statistics - represents an approximation of true satellite image patches. Access to actual georeferenced image chips would enable learning of genuine spatial patterns at the field level.

The price stability classification, even after the forecast-cutoff separation and nearest-market price assignment (Section 3.5; Section 4.5), remains largely predictable from its own constituent input features (Section 5.3): a trivial logistic-regression baseline on three input scalars performs statistically indistinguishably from the full TLDT model. We recommend an independent ground-truth label source, such as Cadre Harmonisé or IPC phase classifications (Global Partners, 2021), as a necessary next step before price-classification performance on this task is presented as evidence of genuine economic forecasting capability.

Finally, the focus on sorghum and millet in the Lake Chad Basin limits generalizability to other crops and regions, and the baseline cross-validated performance documented in Table 1 should be read as an honest empirical anchor for future architectural and data improvements rather than as a demonstrated advance over simpler models on this task.

6.7 Lessons learned for future researchers

Beyond the specific results reported above, this research surfaced several methodological lessons that generalize to other deep-learning applications in data-sparse, real-world settings, which we summarize here for other researchers working in similar conditions.

Audit data pipelines for silent fallback behaviors. In complex multi-stream architectures, API connectivity drops or extraction failures can cause scripts to silently substitute historical averages or placeholder code without throwing visible errors. We highly recommend that any pipeline combining multiple real-time external data streams implement explicit validation assertions or data checksum steps to enforce hard script terminations rather than silent data substitutions.

Count independent units, not rows, before choosing model complexity. Our dataset contained 288 rows but only 18 independent (country, year) units with respect to the yield target. Model capacity and evaluation design should be matched to the latter number, not the former; a multi-head Transformer is difficult to justify statistically at n = 18.

Use group-aware, leave-one-group-out evaluation whenever rows share a coarser-grained label. A single random split can both overstate and understate performance depending on which groups happen to fall into the test set; we observed both directions across different random seeds. Leave-one-year-out cross-validation, or an equivalent group-aware protocol, should be the default whenever multiple rows share a single ground-truth value.

Check label construction for feature leakage before claiming a genuine forecasting target. A hand-crafted label built from a formula that shares variables with the model’s own input features can look like a learnable target while actually being close to a restatement of the inputs. A trivial baseline trained on the same input features is a cheap and effective diagnostic for this failure mode, and we recommend running it before reporting any classification result derived from a composite or engineered label.

Report single-split results only as volatility illustrations, never as headline numbers. Empirical analysis demonstrates that single-split metrics can present a skewed, overly optimistic representation of model performance compared to robust cross-validated baselines. Consequently, single-split outcomes should be treated strictly as diagnostic evidence of sample volatility rather than as citable predictive capacity, especially within highly constrained smallholder datasets.

Prefer a transparent negative or modest result to an inflated one. A transparent, robustly evaluated modest result offers far greater scientific utility to the remote sensing community than an artificially inflated one. Providing an objective account of architectural limits in data-sparse environments establishes a dependable, verifiable baseline that future researchers can reliably build upon.

7 Conclusion

This paper introduced the Tri-Stream Latent-Dynamic Transformer (TLDT), a deep learning architecture for joint pre-harvest yield forecasting and price stability prediction in data-sparse, climate-vulnerable regions. Through purpose-designed components including Adaptive Positional Encoding for irregular satellite time series and Tri-Stream Cross-Attention for multi-modal fusion, TLDT was evaluated on a dataset covering the Lake Chad Basin, using real Sentinel-1/2, CHIRPS, FAOSTAT, and WFP/HDX data throughout.

Under leave-one-year-out cross-validation - the evaluation protocol this research adopts once the true information content of the dataset (18 independent country-year units) was understood - TLDT achieves RMSE = 0.360 ± 0.103 t/ha for yield prediction, performing comparably to a standard Transformer baseline (0.350 ± 0.057 t/ha) rather than clearly outperforming it, and 56.6% ± 19.2% accuracy on three-class price stability classification, which a leakage diagnostic (Section 5.3) shows is statistically indistinguishable from a trivial baseline using only the three exogenous input features. We report these modest, honestly evaluated results as this paper’s primary contribution: a transparent account of what a tri-stream fusion architecture does and does not demonstrably add on real data in a genuinely data-sparse setting.

The Harvest-Expectation Index framework remains a conceptually useful format for interpretable 90-day pre-harvest forecasts, though its operational value depends on the underlying yield and price predictions being more reliable than the present results demonstrate.

Future work should prioritize: an independent ground-truth label for price stability (e.g., Cadre Harmonisé or IPC phase classifications) to resolve the circularity identified in Section 5.3; access to sub-national yield data to test whether the negative-to-near-zero R ² observed here reflects a genuine ceiling on satellite-based prediction at national-yield granularity or an artifact of the label’s coarseness; sub-national yield data that would make a genuine regional breakdown statistically meaningful, given that country-level yield sharing limits any current regional analysis to just 6 independent year-observations per country (Section 5.4); and expansion to a longer study period and additional crops and regions to test generalizability. The modular TLDT architecture facilitates such extensions while maintaining the core tri-stream fusion approach. A condensed set of methodological takeaways for researchers facing similarly small, data-sparse settings is given in Section 6.7.

A further promising direction, suggested during peer review, is reservoir computing. Its low training cost and demonstrated success on continuous, irregularly-sampled time series including food-price and food-security forecasting in settings closely related to this study (Domingo et al., 2023; Herteux et al., 2024) make it an attractive lower-complexity alternative to benchmark against TLDT and the Transformer/LSTM baselines used here. This comparison is particularly relevant given this paper’s own finding that a more complex attention-based architecture did not demonstrably outperform simpler baselines on the present dataset (Section 6.1): reservoir computing offers a further, substantially cheaper point of comparison for testing whether architectural complexity is warranted for this task at the current data scale.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Satiprasad Sahoo, International Center for Agricultural Research in the Dry Areas Egypt, Egypt

Reviewed by: Bilel Ammouri, University of Carthage, Tunisia

Oday ALRahman, University of Technology, Iraq

Duccio Piovani, World Food Programme, Italy