On-farm soybean yield estimation using earth observation data and machine learning models

Highlights

• This study utilizes an extreme boosting model and vegetation indices for soybean yield prediction using the “leave a field out method”.
• Compare spatial generalization through maps.
• This study also utilizes SHAP-based analysis to understand input features in predicting yield.

1 Introduction

Soybeans yield estimate can help farm managers at make multiple decisions on their farms, such as when and where to apply fertilizer, fungicides, and irrigation, as well as planning for harvest, marketing, and storage (Joshi et al., 2023). This early yield information can also be useful to government organizations for policy development related to trade, market price fixation, and food security.

Multiple methods of yield estimation currently exist, such as a labor- and time-intensive method of counting plants per acre, then pods per acre, then seed per acre, and weight per acre, which is then extrapolated to the entire field (Engledow and Wadham, 1923; Wei and Molin, 2020; Vogel et al., 2021). This method works well for small areas within a field, but traditional methods are labor-intensive and takes a lot of time. Because of this, traditional methods cannot be used for the entire field. Additionally, traditional methods cannot create high-resolution yield prediction maps. Another method that seeks to overcome the limitations of the previous method is the application of Earth observations (EO), such as optical (spectral reflectance from satellite) and non-optical (digital elevation model, DEM) data, into machine learning (ML) models to estimate soybean yield (Stepanov et al., 2020; Mohite et al., 2023; Schwalbert et al., 2020; Joshi et al., 2023). Using optical satellite data, various vegetative indices (VIs) and soil indices (SIs) derived from spectral reflectance have been identified as accurate indicators of soil and plant characteristics (Birth and McVey, 1968). Various VIs that can be derived are Normalized Difference Vegetation Index (NDVI), Green Normalized Difference Vegetation Index (GNDVI), Soil-Adjusted Vegetation Index (SAVI), Green Chlorophyll Index (GLI), Difference Vegetation Index (DVI), and Modified Soil-Adjusted Vegetation Index (MSAVI), which can track different physiological traits of the crop, in-season crop growth, health, and vigor patterns of the crops over time across various growth stages (Chen, 1996; Gitelson et al., 2002; Gitelson et al., 2003; Haboudane et al., 2004; Gitelson et al., 2005; Matiza et al., 2024). For example, GNDVI is effective at detecting variations in chlorophyll content and is less prone to saturation than NDVI, making it particularly useful during full canopy closure later in the growing season (Hunt et al., 2011). Additionally, VIs such as SAVI and MSAVI minimize soil background effects by incorporating red and near-infrared reflectance, thereby enhancing accuracy in areas with sparse vegetation or exposed soil (Huete, 1988). Similarly, soil indices derived from bare-soil imagery can be used to characterize soil properties (Al-Kaisi et al., 2017). For instance, the brightness index (BI) represents the overall reflectance intensity of the soil surface (Ibáñez-Asensio et al., 2013), the saturation index (SI) describe the spectral slope, the hue index (HI) captures dominant color characteristics, the coloration index (CI) relates to soil color variation induced due to different organic matter and moisture content, and the redness index (RI) is associated with the presence of iron oxides such as hematite (Mathieu et al., 1998; Sunoj et al., 2024). Similarly, using non-optical DEM data provides topographic information that reveals landscape differences within a field, such as land slope, which affects water flow and, in turn, the amount of water at different landscape positions (Battaglia et al., 2024). Aspect influences the orientation of sunlight, affecting the temperature and microclimate conditions. Hill-shade provides a graphical representation of topography, highlighting regions with varying sunlight intensity. Similarly, changes in elevation directly impact the distribution of temperature and precipitation (Leuthold et al., 2022; Gxasheka et al., 2023; Kumhálová et al., 2011).

Many studies have explored the application of VIs to ML models for soybean yield estimation because ML models enable the analysis of greater volumes of complex data and support the in-depth understanding of non-linear relationships between predictors and crop yields (Amankulova et al., 2024; Bolton and Friedl, 2013; Egbueri et al., 2023; Carneiro et al., 2023; He et al., 2020). For example, Joshi et al. (2023) used 10 different VIs derived from high-resolution satellite imagery at various soybean growth stages as input to a multilayer perceptron deep learning model to predict the field for each 100 m² (10 x 10 m) area. Various other studies have used ML models for yield estimation, including linear regression (LR), multiple linear regression (MLR), random forest (RF), support vector regression (SVR), k-nearest neighbors (KNN), deep neural networks (DNN), and convolutional neural networks (CNN) (Belgiu and Drăgut, 2016; Jodas et al., 2023; Li et al., 2023; Sarkar et al., 2023; Gauriau et al., 2024; Raju et al., 2024). Most of these studies have used VIs to understand physiological and growth patterns and to inform ML-based yield predictions. With advancements in data analytics and computational capabilities, ML models have significantly improved crop yield prediction. However, the black-box nature and the lack of model transferability and validation in unseen fields remain limiting factors of ML models in yield prediction. This is happening because proper validation methods are not being used. For example, Joshi et al. (2023) used the same field data for both training and testing, which may have overestimated model performance because the training and test data originated from the same spatial domain. Such yield predictions may not perform well when applied to new fields with different soil types, landscapes, and microclimates, thereby reducing model transferability. In addition, due to the “black box” nature of the ML model, this study was unable to interpret how, to what extent, and in what direction different input variables affected soybean yield predictions.

While many studies have extensively used optical remote sensing-based Vegetation Indices (VIs) for predicting soybean yields (Amankulova et al., 2024), the role of topographic and soil information (SI) variables is often overlooked. These variables provide insights into within-field landscape variation, including soil characteristics, water flow, soil color, and mineral content, which are critical for understanding yield variability. Despite their importance, the application of these factors in soybean yield estimation remains inadequately studied, highlighting a gap in current research. Moreover, the leave-one-out approach, which evaluates the generalization and transferability of ML models, and new model interpretability techniques are needed for soybean yield estimation. Therefore, the objectives of this study were: 1) to assess the robustness of a single ML model, i.e., XGB, in predicting soybean yield, 2. To evaluate model generalizability, by implementing a novel approach of a spatio-temporal holdout validation to enable model transferability across both spatial (field) and temporal (year) variability, and 3) To propose a new approach of ML model interpretability technique, i.e., Shapley Additive Explanations (SHAP) analysis, to interpret the relative contribution of different predictors.

2 Materials and methods

2.1 Study sites

Three on-farm agricultural fields located in Edmunds, Hamlin, and Miner counties, South Dakota (SD), in the West Northcentral U.S., were studied during the 2019 and 2021 growing seasons (Figure 1) (Joshi et al., 2023). These counties were located at the border of the Dfa (hot summer humid continental climate), Dfb (warm summer humid continental climate), and Bsk (cold semi-arid) Koppen climate zones (Joshi et.al, 2023). Each field followed a corn-soybean rotation, and was managed using standard regional tillage and planning practices. Weather data (temperature and precipitation) were obtained from the nearest NOAA weather stations in Edmunds (45.441975°, -98.751588°), Hamlin (44.72683°, -97.02859°) and Miner Counties (44.004956°, -97.629446°) located 28.16 km, 24.14 km, and 10.46 km away from the fields, and soil texture information was obtained from USDA-NRCS, (2024) county soils surveys. Key soybean planting, harvesting, and field characteristics are provided in Table 1. Additional information about the study sites is also provided in Joshi et al. (2023).

Table 1. Summary of the study sites, including location, area, study year, planting, harvest time, and soil texture.

2.2 Soybean yield data

This study used yield data collected by Joshi et al. (2023). Yield data were collected using a John Deere combine harvester equipped with a yield monitor system with a 9.1-m header width and calibrated according to protocols outlined in Kharel et al. (2018). Raw yield data were processed using SMS Advanced software (Ag Leader Technology, Ames, IA) to remove sensor noise, correct start and end pass delays, and adjust values to 13% moisture (Joshi et al., 2023; Ag Leader Technology, Ames, IA, USA).

2.3 Soil and vegetation indices using satellite imagery

In this study, PlanetScope multispectral imagery (3.12×3.12-m resolution) and georeferenced yield data were provided by Joshi et al. (2023). PlanetScope satellite provides multispectral bands, including Blue (455–515 nm), Green (500–590 nm), Red (590–670 nm), and Near Infrared (780–860 nm). The multispectral bands from the PlanetScope satellite aid in various applications: 1. Blue (455–515 nm) band is useful for studying water quality and detecting phytoplankton. 2. Green (500–590 nm) band is useful for sensitive to healthy vegetation, helping assess plant health and map vegetation types. 3. Red (590–670 nm) band is essential for evaluating vegetation health and calculating indices like NDVI to monitor land cover changes. 4. The near-infrared (780–860 nm) band is critical for differentiating land cover and assessing plant health, especially in agriculture. Together, these bands enable detailed analysis of land use, agriculture, and ecological changes.

For this analysis, we selected two satellite images: one before planting, representing early bare-soil imagery, and another at the R4/R5 growth stage, characterized by full pods and seed initiation. The early imagery was used to derive SIs, including Saturation Index (SI), Coloration Index (CI), Hue Index (HI), and Redness Index (RI), to capture soil brightness, mineralogy, and moisture-related properties before canopy closure. These indices help characterize field-level soil variability influencing water and nutrient availability. Moreover, imagery at the R4/R5 growth stage was selected for its superior predictive performance, as concluded by Joshi et al. (2023). Different ten VIs were derived from selected satellite imagery to characterize various physiological aspects of crop growth. Detailed information on the VIs, including their mathematical equations, is provided in Supplementary Table 1S.

2.4 Topographic variables

A Digital Elevation Model (DEM) with a 1-meter-per-pixel resolution was obtained from the USGS website (U.S. Geological Survey, 2025). Collected DEM was then processed in QGIS (v3.x) to derive topographic variables. Elevation was obtained directly from the DEM. Slope and aspect were computed using the Raster- Terrain Analysis toolset, which first computes first-order terrain derivatives from local elevation gradients. Contour lines were generated from the DEM using the Contour tool at a specified interval to represent elevation changes across the landscape. All derived layers were generated in the same spatial reference system and resolution as the original DEM to maintain consistency across subsequent analyses. Each topographic variable provided different information about the sites. For example, the slope of the land affects water flow, which, in turn, affects the amount of water at different landscape positions. Aspect influences the orientation of sunlight, affecting the temperature and microclimate conditions. Hill-shade provides a graphical representation of topography, highlighting regions with varying sunlight intensity.

2.5 Database development, feature preprocessing

The overall steps and workflow used in the analysis are presented in Figure 2. Following the collection of yield data, SIs, topographic, and VIs, all three study sites were divided into 10x10 m grids of interest. The number of grid cells varied by field size, with 3, 281, 5, 246, and 3, 614 at the Edmunds, Hamlin, and Miner sites, respectively. Across the two study years and three sites, the total number of grid cells was 24, 282. For these, all variables in each grid, including topographic, SI, VIs, and yield-monitoring data, were aggregated using the “Rasterio” library in Python (Bosch, 2019). On average, across all six site years, each field grid cell consisted of approximately 5 to 8 yield data points.

To ensure comparability among predictor variables all feature values were scaled between 0 and 1 using the Max-Min Scaler normalization. Vegetation indices, soil features from Table 2, and topographic attributes (e.g., elevation, curvature, aspects, slope, and hillshade) were selected to estimate soybean yield using an ML model. Following feature selection, we applied a combination of variable metrics to estimate soybean yield. Combinations of variables have been designed to provide unique perspectives on insights into factors influencing crop performance (Genuer et al., 2010). Yield (kg ha^-1) was standardized (z score) for model training, and back-transformed for interpretation. Two scenarios were constructed for estimating yield:

Table 2. Soil indices derived from pre-planting PlanetScope multispectral imagery, including mathematical equations, spectral bands used, and literature references.

These indices were calculated from early-season bare-soil imagery to characterize field-scale soil properties potentially associated with soybean yield variability.

BI, Brightness Index; SI, Saturation Index; HI, Hue Index; CI, Coloration Index; RI, Redness Index. Blue, Green, and Red represent PlanetScope surface reflectance values for the blue (455–515 nm), green (500–590 nm), and red (590–670 nm) spectral bands, respectively. Soil indices were calculated from pre-planting bare-soil imagery before canopy closure to minimize vegetation effects.

• VI-only, in this scenario, only 10 different VIs that were derived using satellite imagery at R4/R5 growth stage were used for yield prediction and
• topography + soil indices + VIs: In this scenario, a model was developed by using a comprehensive dataset that includes topographic variables, SIs, and VIs.

2.6 Model implementation and validation

An extreme gradient boosting (XGBoost) model (Chen and Guestrin, 2016) was selected for yield prediction due to strong performance with nonlinear multivariate datasets. It combines both regression and classification approaches. The XGBoost model achieves strong performance by combining multiple weak models (Bentéjac et al., 2021). It iterates in that way, where decision trees are used to reduce the loss function and make attempts to reduce error from previous decision trees (Shafi et al., 2023). To obtain the best hyperparameters (e.g., learning rate, max depth, and subsample) of the XGBoost model, the GridSearch library was used to test several combinations of values. To evaluate generalization across locations, we implemented the “Leave-One-Field-Out Cross-Validation” (LOFO-CV). In general, LOFO-CV cross-validation means that datasets were combined across years and trial fields. LOFO-CV uses the entire dataset for training, except for the field, which is used for testing. LOFO-CV tested the XGBoost model performance of the unseen trial field. LOFO-CV helps to generalize the model and reduce the black-box problem in machine learning. Table 3 summarizes the best-performing ML models and suitable hyperparameters for estimating soybean yield in different trial fields (Edmunds, Hamlin, and Miner) using the LOFO-CV at R4/R5 for each of the two feature matrices identified in section 2.5. The XGBoost model uses 300–500 estimators, a learning rate of 0.05-0.7, colsample_bytree of 0.7-0.9, and subsample of 0.7-0.9. The coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE) were used to assess model performance (i.e., compare estimated and actual yields on the test data).

Table 3. Best model parameters of the XGBoost model for estimating soybean yield from PlanetScope-derived vegetation indices (and topographic/soil features) under Leave-One-Field-Out Cross-Validation (LOFO-CV).

Model performance is reported by comparing satellite-based yield predictions with field-measured yield monitor observations.

Test R², RMSE (kg ha^-1), and MAE (kg ha^-1) were calculated by comparing predicted soybean yields generated from PlanetScope-derived features with observed yields collected using a calibrated combine yield monitor. Lower RMSE and MAE and higher R² indicate better agreement between satellite-based predictions and field measurements. Train and Test refer to the training and unseen test fields used in LOFO-CV.

In this study, Python (version 3.12) and Jupyter Notebook (version 7.3.2) were used along with several libraries, including scikit-learn (version 1.6) for ML, and Pandas, NumPy (version 2.2.1), and matplotlib - for data handling and visualization.

2.7 Model interpretation

Shapley Additive exPlanations (SHAP values) provide an interpretation of model predictions (Lundberg & Lee, 2017). To improve interpretability, SHAP values were computed only for test-set data in each LOFO fold, preventing information leakage. In our study, SHAP analysis was performed on the tested dataset. For each LOFO-CV fold, the best XGBoost model was selected using an inner GroupKfold approach. The top 10 vegetation indices were retained per fold based on the model-derived feature importance. SHAP values were calculated only for test Fold observations with TreeExplainer, which is optimized for tree-based models.

3 Results

3.1 Yield and weather patterns

The effect of spatial agroclimatic zone, agronomical factors, amount of precipitation, and temperature on yield can be observed in Figure 3. Across all trial fields, soybean yield varied substantially between 2019 and 2021. In 2021, Hamlin exhibited the highest mean yield (3, 868 kg ha^-1), followed by Edmunds (2, 897 kg ha^-1), while Miner showed a pronounced yield reduction (556 kg ha^-1). In contrast, 2019 yields were more balanced across sites, with Edmunds highest (2, 987 kg ha^-1) and Hamlin lowest (2, 776 kg ha^-1) (Table 4). Figure 3 shows the distributions, which approximate normality and reflet field-specific variability. More details of the yield can be found in Joshi et al. (2023).

Table 4. Summary statistics of yield across all trial fields for 2019 and 2021.

Bi-monthly average temperature (orange), growing degree days (GDD; black), and precipitation (blue dashed line) for the Edmunds, Hamlin, and Miner study sites during the 2019 and 2021 growing seasons. Overall, the 2021 season was characterized by lower precipitation and higher temperatures across parts of the study region, particularly at Miner, consistent with the reduced soybean yields observed during that season (Figure 3).

3.2 Relationship between yields and Features Variables

Vegetation Indices (VIs): Correlation between VIs and soybean yield varied by field, year, and the R4/R5 and R6/R7 growth stages (Figure 4). Across all location, VIs from the R4/R5 stage consistently showed strongest correlations with soybean yield, with NDVI, RVI, and DVI reaching correlations of 0.5-0.7 at Edmunds and Hamlin in both seasons. At Miner, VI-yield correlations were relatively weak in 2019 but substantially stronger in 2021, likely due to differences in weather, soil moisture, or management practices. Additionally, spatial variability is evident, with Edmunds displaying more consistent VI-yield relationships compared to the trial field of Hamlin and Miner.

Soil Indices (SIs): Figure 5 shows only three soil indices with high correlation; other soil indices were eliminated from the figure based on similar physical and chemical characteristics. Soil indices demonstrated moderate correlation with crop yield (Figure 5), particularly with CI and RI (r ~ 0.59 to 0.60) at Edmunds. Hamlin showed minimal SI-yield correlation with yield, in both 2019 and 2021.

Topographic Variables: Topographic attributes exhibited generally weak with yield (Figure 5). Elevation had the highest correlations (r = 0.10 in 2019; and 0.38 in 2021), while slope and hillshade and aspect, showed minimal SI-yield relationships in both years.

Overall, R4/R5 VIs provide the strongest and most consistent predictors of yield across fields; soil indices contribute meaningful information at specific sites, while topography plays a modest and site-dependent role.

3.3 Machine learning model performance (XGBoost)

The LOFO-CV results for the VI-only and, topography + soil indices + VIs model is summarized in Table 5 for the trial fields of Edmunds, Hamlin, and Miner. Model performance varied across trial fields and years, reflecting the differences in landscape conditions and inter-annual weather patterns.

Table 5. Performance of the XGBoost model for estimating soybean yield using PlanetScope-derived vegetation indices (VI) alone and a combined feature set (Topography + Soil Indices + VI) under Leave-One-Field-Out Cross-Validation (LOFO-CV).

Model performance is reported by comparing satellite-based yield predictions with field-measured soybean yields obtained from a calibrated combine yield monitor across the Edmunds, Hamlin, and Miner trial fields during the 2019 and 2021 growing seasons.

Train R² and Test R² represent the coefficient of determination for the training and unseen test datasets, respectively. RMSE (kg ha^-1) and MAE (kg ha^-1) quantify the difference between predicted yields derived from PlanetScope imagery and observed yields measured by the combine yield monitor. Higher R² values and lower RMSE and MAE values indicate closer agreement between satellite-based predictions and field measurements. LOFO-CV, Leave-One-Field-Out Cross-Validation; VI, Vegetation Indices; SI, Soil Indices.

Best performance: Edmunds (2021) with R² = 0.54, RMSE = 574 kg ha^-1, and MAE = 456 kg ha^-1. Moderate Performance: Hamlin across both years (R² = 0.20-0.39).

Lowest Performance: Miner, where negative R² values were observed due to yield instability and z-score back-transformation effects.

The topography + soil indices + VIs generally did not improve predictive accuracy under the LOFO-CV approach, suggesting that topographic and soil features did not generalize well across locations.

Spatial prediction maps (Figure 6) show that areas of under- and over-estimation correspond to localized field conditions, reinforcing the challenge of field-to-field transferability). Figure 6 shows only the best performance model at Edmunds uisng vegetation indices (VIs) matrix. Spatial maps show the combined years (2019 and 2021) yield prediction and error. Spatial maps show the combined years (2019 and 2021) yield prediction and error. Spatial prediction performance of the XGBoost model is shown in Figure 6. The yield model is applied under spatial generalization and compares observed and model-predicted yields, along with their residuals, across three trial fields (Edmunds, Hamlin, and Miner). Observed yield in the figure represents measured yield and captures within the trial field spatial heterogeneity driven by local field conditions. The predicted maps show the spatial distribution of yield estimated by the XGBoost model using vegetation indices. While residual maps depict prediction errors, where values close to zero indicate accurate predictions, positive and negative values correspond to under- and overestimation of yield, respectively. LOFO-CV reveals that R4/R5 VI-based models transfer reasonably well to fields with similar conditions (e.g., Edmunds) but struggle at locations with extreme variability (e.g., Miner).

3.4 SHAP-based feature importance

SHAP analysis identified the vegetation indices DVI, RDVI, GCI, TGI, and GNDVI as the top vegetation indices (VIs) and most important predictors of yield across LOFO folds (Figure 7), with MSAVI2, NDVI, and SAVI ranked as least important. The SHAP values indicate the impact of input features on yield prediction, e.g., vegetation indices. A positive SHAP value means a higher likelihood of increased yield prediction, and a negative SHAP value favors a decreased yield prediction. The widest range of DVI SHAP values among all vegetation indices underscores the importance of DVI SHAP for yield prediction. The SHAP value range varies across vegetation indices. For example, DVI shows the largest SHAP value range (approximately −2 + 1), whereas MSAVI2 and RVI show the smallest positive and negative SHAP value ranges (approximately 0). In Figure 7, DVI consistently dominates model contributions, highlighting its importance for capturing canopy vigor during the critical R4/R5 growth stage.

4 Discussion

Soybean yield is affected by many factors such as weather patterns, management practices, and soil properties. In general, it is assumed that higher VI values in a given year may be indicative of higher crop yields. Similarly, favorable soil conditions tend to produce higher yields. Our studies indicated that soybean yield was affected by extreme weather at the trial field in Miner in 2021, because this year was a drought year in South Dakota. It has been observed that the yield was the lowest among the trial fields. Although correlation values varied with yield, vegetation indices were more clearly delineated in 2021 than in 2019. Similarly, Mishra and Cherkauer (2010) reported that drought conditions impact soybean yield, and higher rainfall regions are typically more favorable for crop production. Guo et al. (2022) reported that lower precipitation and the number of hot days with temperatures above 35 °C had a significant negative impact on the major yield-determining components of the soybean crop: the number of pods and seed weight. These studies also indicating similar pattern to our study.

Additionally, GDD helps in understanding plant structure and leaf formation by providing data for each day’s growth Akyuz et al. 2017). Since GGD is a crucial parameter in determining both vegetation growth and yield production, accumulative GDD is used to estimate crop development. Typically, 1, 500 to 2, 050 GDD (833.3 °C - 1138.9 °C) are required to optimize soybean yield (Santos Hansel et al., 2019); however, this also depends on the soybean variety and climatic conditions. Soybean yield was also examined using GDD. While vegetation, soil, and topographic variables were all evaluated in the prediction models, vegetation indices consistently provided the strongest predictive performance, particularly during the R4–R5 growth stage.

VIs (e.g., NDVI, DVI, GNDVI – things that measure the plant not just the environment) outperform combined features (VIs + soil indices + topography) because they:

• directly capture plant condition and integrated stress signals (chlorophyll content, canopy structure, biomass and vigor)
• temporally aligned with yield formation because yield is the outcome of plant growth. As a result, they provide consistent signals across fields that are transferable. In contrast, the soil terrain variables (e.g., indices and topography) are indirect drivers of yield (influence yield through the plant) and introduce noise ad weak relationships (mismatches). Because of the limited amount of data, their effects are nonlinear and context dependent.

So VIs integrate multiple factors (soil fertility and moisture, management practices, weather impacts) and act like a summary signal of all the underlying processes combined; they encode the net effect.

As a result, the VIs provide the best temporal alignment (i.e., showed strongest correlations at the R4/R5 growth stage) because the plant stress signals are visible in the reflectance. Our R² reflects that our models overfit local conditions (Figure 8).

Our models overfit local conditions (because our soil/topography is highly site-specific, with relationships varying across fields, e.g., slope matters in one field, but not another). These variables dilute the stronger signals of the VIs, confuses the tree models, reduces stability across folds, and increases noise and model complexity. As a result, we have poor transferability to new fields. Thus, VIS are more consistent across locations, providing better generalization.

Adding more features, such as soil indices and topographic variables, did not improve the model performance compared with vegetation indices alone. However, this does not mean that soil and topographic information are unimportant. These variables provide useful background information and could be valuable in process-based crop models or as prior information. The results also showed that prediction accuracy depended strongly on image timing. The highest accuracy was achieved during the R4–R5 growth stages, when crop conditions were most closely related to final yield. Missing these critical growth stages or having fewer satellite images reduced the model’s ability to capture crop stress and development, leading to lower prediction accuracy.

Overall, our results suggest that capturing the plant’s dynamic response at the right time is more important than adding static environmental variables. LOFO-CV provides a novel approach to estimating soybean yield using machine learning models; it works on unseen data and across fields. In our studies, we used LOFO-CV, and the results are promising: when we tried to predict yield on unseen data, many parameters matched those of the trained field. In our case, Edmunds’ trial field performed better than Miner’s trial field. In many cases, results were recorded as overstated and underestimates due to changes in field locations. Our study also observed negative R-squared values during testing, which were due to back-transformation. Desloires et al. (2023) also reported similar outcomes using a “field-out-year” approach using an XGBoost model. This is clearly seen from these studies that two or three years of field data are insufficient to accurately utilize the LOFO-CV approach.

In addition, a two- to three-year dataset can be interpreted in many ways, such as comparing different machine learning models, but the simplest way is to perform SHAP analysis. Because one of the most challenging aspects of machine learning models is the black-box nature of hidden layers, this raises questions about the model’s performance and reliability. SHAP analysis provides answers of important variable and ranks them. DVI was identified as the most influential predictor, suggesting that canopy vigor during the reproductive stage contains the most useful information for estimating soybean yield. Prioritizing these high-impact vegetation indices can simplify future models, improve their interpretability, and enhance operational satellite-based yield monitoring for precision agriculture. In our study, Global SHAP analysis reveals that DVI is a top-ranked feature when using the LOFO-CV approach to predict yield, and years of research support similar findings (Figure 7). SHAP analysis captures the importance of variables in the model and their contribution to yield estimation. It also helps to rank these variables by their negative and positive value ranges. In machine learning, it enhances transparency by demonstrating that high accuracy can be achieved with a few variables.

4.1 Limitations and challenges of the study

The use of climate and topography data at the field level is still not widely studied, even though both factors greatly affect crop growth and yield. Climate conditions like temperature, rainfall, and humidity interact with topographical features such as elevation, slope, and aspect. These interactions influence important factors like water drainage, sunlight exposure, and soil moisture, which are essential for crop productivity. However, yield data collected from fields can be inconsistent due to differences in sampling methods, machine calibration, and challenges during harvest. These inconsistencies add noise to ML models, making it harder to estimate yield accurately. To improve model reliability, it is important to establish standardized methods for collecting yield data.

Frequent satellite data is also crucial for tracking rapid crop growth, especially during key growth stages. If satellite images are taken too far apart, important events like pest attacks, irrigation, or fast growth phases may be missed. This limits the accuracy of crop monitoring models, particularly for short time periods. Additionally, within a single field, variations in soil quality, nutrient levels, and farming practices can impact machine-learning models. If models do not consider these differences, they may struggle to provide accurate estimations across the entire field. To improve accuracy, more detailed field-level data should be collected and included in predictive models.

5 Conclusion and future work

5.1 Conclusion

This study demonstrates that growth-stage-specific modeling substantially improves the accuracy of soybean yield estimation from Earth observation (EO) data. The stronger correlations observed at the R4/R5 stage indicate that earlier-season vegetation indices should be prioritized. NDVI, RVI, and DVI emerged as key indicators for nitrogen application, irrigation scheduling, and harvest planning, whereas NDWI requires cautious interpretation due to occasional negative correlations with yield. Soil indices from the VE/VC stage, when combined with topographic features, contributed meaningfully to yield estimation at R4/R5.

A key novelty of this work is the integration of vegetation indices, soil indices, and topographic metrics within a growth-stage-focused ML framework, calibrated across trial fields in different counties. A “LOFO-CV” approach allowed XGBoost to outperform at the trial field, Hamlin and Miner. Compared with previous studies that focus on single-feature groups or full-season averages, our method explicitly accounts for temporal dynamics in crop development and local variability, thereby providing a more operationally relevant yield-prediction tool. Integration of topography, soil indices, and vegetation indices did not perform well as expected while using the “LOFO-CV” approach with the XGBoost model. It shows limitations due to across varied environmental and management conditions.

5.2 Future work

Future research should integrate high-frequency satellite imagery and UAV data with climatic variables such as temperature and precipitation to improve temporal resolution and stress detection. Incorporating localized soil nutrient profiles, management practices, and finer-scale spatial variability will further enhance prediction accuracy. Expanding the framework to other crops and agroecological zones will enable a broader assessment of its applicability in precision agriculture.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fagro.2026.1923239/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Marija Popovic, Delft University of Technology, Netherlands

Reviewed by: Huong T. X. Nguyen, McGill University, Canada

Nurlan Kurmanov, L.N. Gumilyov Eurasian National University, Kazakhstan