Introduction

Precision management at the subfield scale using remote sensing and big data technology (i.e. field sensors on combines, planters, and applicators) holds promise for sustaining and increasing crop yields whilst helping to reduce negative environmental impacts of intensive agricultural production systems (Liu et al., 2022; Veldhuizen et al., 2020; Weiss, 2001; Shuai & Basso, 2022). Use of yield stability maps and remote sensing vegetation indices to observe in season dynamics is needed to better understand the mechanistic drivers of yield heterogeneity. Such real time information equips producers to improve yield and production efficiency while reducing nutrient pollution and agricultural greenhouse gas (GHG) emissions (Basso et al., 2019). The concurrent availability of high-resolution remote sensing imagery, advanced computing, and high quality, ground truth yield data opens the door to reliable assessment of subfield scale crop growth and performance (Deines et al., 2021; Hatfield et al., 2020). Immediate, major questions and objectives are focused on the agronomic challenges of applying optimum inputs, in particular nitrogen (N) applications to maximize maize yields (Hunt et al., 2019; Ji et al., 2021; Jin et al., 2019; Kayad et al., 2019) while mitigating N losses through leaching and GHG emissions (Basso et al., 2019; Liu et al., 2022). Use of these approaches and tools can also enable a systemic analysis of food security, landscape biodiversity, and land use diversification informed by the intrinsic productivity capacity of the landscape (Hernández-Ochoa et al., 2022).

The heterogenous response of the landscape to agricultural management can be separated into homogeneous yield stability zones. Characterization of common yield zones with similar soil properties, topographic features, and historical yield has been conducted for decades (Nawar et al., 2017; Plant, 2001). The yield stability zones (YSZ) segment the field relative to stable and unstable yields based on temporal variations of crop yield in specific areas. The Stable zones can be further classified into high and low yield levels. This analysis produces high and stable (HS), low and stable (LS), and unstable (UN) yield stability zones, but does not provide distributed absolute yield. In any specific year, there are high-yield and low-yield areas in the unstable zone, caused by interactions of weather, location, and management. Identification of yield stable zones helps improve nitrogen use efficiency (NUE) and farm profitability by enabling farmers to adopt more precise N management strategies that match the yield potential of these stable zones, increasing resource use efficiency (water, nutrients) (Basso et al., 2019). Unstable zones require in-season real time remote sensing imagery to inform tactical management, as yields in these zones are driven by seasonally variable factors such as interactions of hillslope position and precipitation (Maestrini & Basso, 2018; Martinez-Feria & Basso, 2020). Despite this solid finding, we still do not know if the reflectance of a pixel in a low and stable zones would lead to a final low yield or could change to produce high yields, if located in stable zones.

The estimation of post-hoc crop yield for maize has been extensively studied and validated in the literature, including on smallholder fields and at the landscape scale (Deines et al., 2021; Jin et al., 2019; Padilla et al., 2018). However, to make timely in-season management decisions, real-time spatial information about crop conditions at the sub-field scale is crucial (Shuai & Basso, 2022). Several commonly used vegetation indices, such as the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Ratio Vegetation Index (RVI), are widely used to represent crop biophysical parameters like leaf area index (LAI) and chlorophyll content (Clevers & Gitelson, 2013; Hassan et al., 2019; Xue & Su, 2017). However, these indices have limitations, particularly in high LAI settings typically found in maize fields, as they tend to saturate and lose sensitivity. In contrast, the Green Chlorophyll Vegetation Index (GCVI) exhibits less saturation and provides a more reliable representation of crop condition (Gitelson et al., 2003).

Deines et al. (2021) successfully utilized GCVI to model corn yields relative to onboard yield monitor data. Their approach involved incorporating weather data, smoothing the data with harmonic best fit regression, and utilizing two images, specifically at the peak of the harmonic GCVI and 30 days thereafter. Building upon this foundation, in our study, we employed GCVI as an indicator of crop condition. By leveraging GCVI, we aimed to capture the dynamic nature of crop growth and monitor its progress throughout the season. The use of GCVI provides valuable insights into crop health and vigor, enabling more accurate and timely assessment of crop conditions.

While efforts have been made to develop empirical yield estimation models using remote sensing images, their generalizability, especially at the subfield scale, remains limited, although certain patterns have been observed. Johnson (2014) investigated the correlations between various satellite products, such as NDVI, surface temperature, and rainfall, with corn and soybean yields in the United States. The strongest correlations were found with NDVI images taken in late July to early August. Lobell et al. (2015) utilized multi-temporal Landsat GCVI composites to generate two in-season images, one during the early season (mid-June to mid-July) and one during the peak season (mid-July to late August), and proposed a scalable method for mapping estimated maize yields across the US. Burke and Lobell (2017) replicated this method to detect yield variations in small fields in Kenya, achieving similar accuracy to survey-based measures. Despite significant efforts, the current level of skill and accuracy in estimating subfield crop yields remains low, with only around half of the yield variation being explained (Asseng et al., 2015; Azzari et al., 2017; Burke & Lobell, 2017; Deines et al., 2021).

Translating remote sensing yield estimates into practical in-season management poses challenges, as the marginal accuracy of remote sensing yield models is influenced by the heterogeneity within and between fields, including variations in soil type, conditions, and landscape location. The selection of the optimal image for yield analysis often assumes homogeneity in crop development within and between fields. However, crop emergence and development timing are heterogeneous due to variations in soil quality, soil moisture, and other factors (Albarenque et al., 2023). Delays in emergence can result in subsequent yield reductions. Moreover, using yield stability zones to predict yield levels may not be suitable for optimizing tactical inputs for crop production when the estimation of crop biophysical parameters at specific times during development is not sufficiently accurate due to non-uniform emergence (Albarenque et al., 2023). These challenges give rise to two central research questions: (1) Do remote sensing images provide additional benefits for in-season tactical management beyond yield stability maps? and (2) Does in-season remote sensing of crop growth dynamics enhance the prediction of harvestable yields compared to yield stability zones?

In this study, our primary focus was to evaluate the effectiveness of remote sensing and yield stability methods in estimating yield levels and understanding the influence of in-season crop growth heterogeneity on final yield outcomes. To achieve this, we aimed to accomplish the following specific objectives:

  1. 1.

    Assess the accuracy and advantages of remote sensing techniques in estimating in-season yield levels, particularly in comparison to the historical yield stability zones.

  2. 2.

    Determine the optimal timing of remote sensing image acquisition within the growing season to enhance the precision of yield estimation.

  3. 3.

    Investigate the impact of spatial-temporal variations in a remote sensing vegetative index on yield levels, providing insights into the relationship between crop growth dynamics and final yield outcomes.

By addressing these objectives, we aimed to enhance our understanding of the capabilities and limitations of remote sensing methods for in-season yield estimation and to gain insights into the significance of within-season crop growth heterogeneity on yield outcomes. This study provides evidence to support future research to improve precision management and crop yield predictions by integrating historical crop yields and remote sensing techniques.

Study area and dataset

Study area

Our study focused on assessing maize yield heterogeneity in the U.S. Midwest region, specifically across 115 fields located in four states: Michigan (MI), Indiana (IN), Iowa (IA), and Minnesota (MN). The selected fields represented a diverse range of agricultural landscapes within the Corn Belt region. During the years 2018, 2019, and 2020, we analyzed a total of 34 fields in Michigan, 29 fields in Indiana, 35 fields in Iowa, and 17 fields in Minnesota (Table 1). These states experience a continental climate characterized by cold winters and hot, humid summers (Modified Köppen classification types Df). In the Midwest, the maize crop is typically cultivated within maize-soybean or maize-soybean-wheat rotations. Our analysis focused specifically on the maize crop season, which begins with planting in late April to mid-May, taking advantage of favorable soil conditions, and harvest in October.

The region’s average monthly rainfall during the crop growing season (May-September) typically ranges from 76 to 127 mm, supplemented by stored soil moisture from the winter and spring months. These moisture conditions generally meet the crop’s water requirements in most years.

By examining fields across these states, we aimed to capture the variability in maize yield patterns within the Midwest region and gain insights into the factors influencing yield heterogeneity at a subfield scale.

Table 1 Number of fields in each state used for yield analysis and field acreage
Full size table

We used various spatial datasets to investigate the accuracy of satellite image and historical YSZ and the impacts of environmental factors on rate of change of the GCVI metric (Table 2).

Table 2 Satellite, yield, and environmental dataset used in this study
Full size table

Remote sensing imagery

Planet provided high-resolution, real-time images for agricultural monitoring in visible (RGB) and near infrared bands (NIR) (Marta, 2018). We utilized multi-temporal, clear and cloudless PlanetScope time series images from 2018 to 2020 for Indiana, Michigan, and Minnesota, while images for Iowa were obtained only from 2018 to 2019 due to the lack of yield monitor data for some fields in 2020. PlanetScope images are gap-filled, radiometrically calibrated, atmospheric and orthorectification corrected, and harmonized to Sentinel 2, prior to acquisition (Gao & Zhang, 2021). Tile numbers are 19E-190 N (Indiana), 22E-197 N and 22E-198 N (Iowa), 28E-195 N (Michigan), and 20E-201 N and 20E-202 N (Minnesota), with each tile containing 8000 × 8000 pixels at 3.125-meter resolution (Fig. 1). For each year, we acquired images for the first and 15th day of June-September with an additional image from September 30th to cover ~ 75% of the maize growing season. We resampled Planet images into 10-meter resolution to match the standard resolution of our environmental dataset.

Fig. 1
figure 1

The study area included 115 field locations within a single Planet labs imagery tile for each of four states, Minnesota, Iowa, Indiana, and Michigan

Full size image

Agronomic and environmental datasets

In our study, we compiled an extensive set of agronomic and environmental datasets to enhance our analysis. These datasets included elevation data, weather data, and soil spatial data. For elevation data, we obtained a 10-meter digital elevation model (DEM) from the USGS National Elevation Dataset (Arundel et al., 2015). From this DEM, we calculated the topographic position index (TPI) to classify hillslope positions (Weiss, 2001). The TPI compares each pixel in the DEM to the average of its neighboring pixels within a 9 × 9 window. Based on the TPI values, we classified each pixel as a hilltop (TPI > 1 SD), depression (TPI < -1 SD), or other, where SD represents the standard deviation of TPI values across all fields. To characterize weather conditions, we utilized gridMET data at an approximate resolution of 4 km (Abatzoglou, 2013). This dataset provided daily averages, maximums, and minimums of annual temperature and rainfall for each field. For soil properties, we obtained the Gridded Soil Survey Geographic (gSSURGO) soil dataset at a resolution of 10 m (Soil survey staff, 2009, https://gdg.sc.egov.usda.gov/). This dataset provided valuable information on soil characteristics, including soil organic matter (SOM) content, bulk density (BD), soil depth, and soil water availability.

To evaluate yield performance at the subfield scale, we utilized yield stability maps for 115 fields across the four states. We acquired high-resolution yield shapefiles based on observed yield monitor data directly from the farmers for the period from 2008 to 2020. These shapefiles were processed into 2-meter resolution yield rasters. The available data from 2008 to 2017 were used to generate yield stability (YS) maps, while the years 2018 to 2020 were reserved for evaluation.

The yield stability maps classified each field into high and stable (HS), low and stable (LS), and unstable (UN) areas. Yield levels were defined relative to the field and year using a normalized multi-year average, where high areas had a normalized average yield > 0 and low areas had a yield < = 0. Stability was determined based on the multi-year standard deviation, with stable high and low areas having absolute yield standard deviations less than one, and UN areas having absolute yield standard deviations greater than one. The calculation details of the yield stability zones can be found in Shuai and Basso (2022). The stability maps were also resampled to a standard 10-meter raster to match the DEM and gSSURGO datasets.

To illustrate our approach, we provide an example of the input yield maps for one field, which includes historical yield maps from 2008 to 2017, DEM, and TPI layers. We also present the output yield stability (YS) map, which demonstrates the delineation of HS, LS, and UN areas (see Box 1 and Fig. 2). This field was planted with corn for five years and soybeans for four years. Missing values from the 2009 yield map correspond to missing yield data from the farmer. In 2016 and 2017, most of the field area produced high yields, while there was large variation in yields in other years. The southwestern part of this field consistently produced higher yields (blue) and was calculated as high and stable (HS). The northern end of the field consistently produced lower yields compared to the field average and was calculated as low and stable (LS). The DEM map in Fig. 2 shows that the upper northern end of the field is about 20 m higher than the southern edge of the field. The spatial patterns of the depression area (green) in the TPI map often overlayed HS areas in the YS map suggesting the accumulation of moisture as one mechanism driving high yields in this field.

Fig. 2
figure 2

Historical yields from 2008 to 2017 and associated yield stability, digital elevation, and topographical index maps for one example field

Full size image

Methods

Yield level prediction using yield stability maps and remote sensing imagery

We evaluated the performance of the yield level prediction by YSZ compared to ground truth yield monitor data (Fig. 3). We defined the accuracy of the yield level as the fraction of stable pixels that were correctly classified as high or low. Over the years 2018–2020, we compared the historical YS estimation of yield level to the observed yield level. The observed yield level was defined by each pixel’s yield value relative to the field mean yield value; pixels lower than the mean were labeled low and otherwise were labeled high. Each pixel location was referenced against the HS and LS zones of the historic YS and the percentages of low and high in each zone was calculated. The historical YS UN zone provided no predictive information so was not included in the YS analysis, therefore approximately 30% of the total field area were not included in the analysis.

We evaluated the accuracy of remote sensing linear regression and random forest models for the prediction of yield and yield level compared to ground truth yield monitor data. We used GCVI with three modeling approaches to leverage remote sensing imagery to estimate current year yield and yield levels. In the first modeling approach, we used linear regression (LR) of yield by GCVI using the single optimal image from the season, where the single, post-hoc optimal image was defined as the image with maximum correlation with yield, referred to as optimal image in this study. The LR modeled yield as

$$Yeil{d}_{ij}= {\beta }_{1}GCV{I\_opt}_{ij}+{\beta }_{0} 1$$

where yield is the predicted yield, ij are the (x,y) cell index, \({\beta }_{1}\) is the coefficient, GCVI_opt is the pixel GCVI value of the optimal images and \({\beta }_{0}\) is the y intercept. For the second modeling approach we used all the images in each growing season sequentially as inputs to a series of multiple linear regression (MLR) models of yield. The MLR modeled yield as

$$Yeil{d}_{ij}= {\beta }_{1}GCV{I\_t1}_{ij}+\dots +{\beta }_{n}GCV{I\_tn}_{ij}+{\beta }_{0} 2$$

Where GCVI_t1 to GCVI_tn represent the pixel ij of each image from period 1 to n. In the third modeling approach, we applied random forest (RF) regression models to estimate yield using the same, sequential image inputs as the MLR approach. Random Forest is an ensemble learning method constructed from multiple decision trees. In each tree, a random subset of GCVI images at different acquisition dates is extracted for making decisions. The estimated yield value is the mean of the prediction from all trees. RF regression was conducted using the R “randomForest” software package (Liaw & Wiener, 2002). We used the square root of the number of input variables to split the data at each node and set a maximum of 500 trees. For all three models, 20% of pixels were randomly selected as training data. Then, the trained models were used to estimate the absolute yield. The estimated yields were translated into high or low yield levels for comparison to the historic YS approach, following the same logic described above.

Optimal image selection for yield estimation using clustering analysis

We tested the hypothesis that the optimal image date differs between the field and the subfield scale. For each field, we computed the correlation between GCVI and crop yield and obtained the optimal image. To assess the sub field scale, we performed a cluster analysis on the GCVI curve over the growing season for each pixel with the time-series K-means (TK-means) method using the tslearn python package (Tavenard et al., 2020). This time-series clustering was carried out once for each field-year combination. Similar to the K-means algorithm, the TK-means evaluates the weighted distance between m time series \(Z=\{ {z_1},{z_2},{z_3},?,{z_m}\}\) and the centroids of k clusters \(C=\{ {c_1},{c_2},{c_3},?,{c_k}\}\), where weights are assigned to time series values according to the importance of a time span in the clustering process. Details can be found in Huang et al. (2016). We used soft dynamic time warping as the distance metric suggested by Janati et al. (2020) to avoid the susceptibility of Euclidean distance to time shifts in the GCVI curves expected in this data set. We used 8 clusters of GCVI curves with different shapes and did not optimize the segmentation as the objective of clustering was to evaluate if there was variation in yield estimation accuracy between clusters. Similar to field level analysis, we obtained the optimal image for each cluster in each field. To evaluate whether the patch (field or subfield) level optimal date or the field itself influenced the yield prediction, we conducted yield estimates for LR models using optimal image for each field and for each field cluster.

Impacts of the rate of change of crop growth on yield

We evaluated the change in GCVI between successive images during the growing season, also called rate of change (ROC) in this study, to get implications from early season remotely sensed growth dynamics for subsequent in-season management. Each GCVI image was classified as high (H) or low (L) based on the field’s average GCVI value. Then the change between successive GCVI images was determined by associating labels from both dates, producing four temporal classes (HH, HL, LH, and LL). For example, a pixel was labelled as high-high (HH) on July 15th if it was classified as high on July 1st and July 15th images. HL/LH indicated areas with changing GCVI status, and HH/LL indicated areas with no change. We computed the state average percentage of change vs. non change over time and across historical yield stability zones to test if short term dynamics in crop growth related to crop yield outcomes.

Results

Accuracy of yield level prediction using yield stability maps and in-season remote sensing

Yield monitor data defined yield stability zones consistent with past findings. Over all the fields in this study, approximately 53% of field area was defined as HS, 31% as LS, and 16% as UN, consistent with results in our previous research (Basso et al., 2019). Estimated state-level average yields in each stability zone demonstrate the magnitude of temporal and in-field spatial heterogeneity of yields. HS produced around 11.3 Mg/ha, 2.4 Mg/ha higher than the LS and UN areas. The UN areas experienced large temporal fluctuations, i.e. In Indiana, yield was higher in LS than UN in 2019 and vice versa in 2018 and 2020. By comparing crop yield of each pixel to its field-level average, for both years, 58% of field area was classified as high yield and the remaining 42% as low yield.

Table 3 Average maize yield in each stability zone
Full size table

Across all four states for HS and LS zones, YSZ outperformed LR and MLR model approaches but not RF models to predict the yield level (Fig. 3). The YSZ predicted crop yield level with an accuracy of 76% without any in-season images. The YSZ prediction accuracy was similar in both stable zones with 8% higher accuracy in HS than in LS for Iowa and Minnesota fields, and 8% lower in HS than in LS for the other states. The LR method achieved an accuracy of 58% in yield prediction using post-hoc best correlated images, while MLR and RF improved the accuracy to 66% and 86% using all the seasonal images, respectively. The RS modeled methods accuracy capture both pixels remaining in their respective YSZ category (i.e. High pixels in the HS zone and those pixels that were correctly predicted changing YSZ category (i.e. a low pixel correctly identified in the HS zone).

Fig. 3
figure 3

Accuracies of maize yield level prediction by linear regression (LR), multi-linear regression (MLR), random forest (RF), and yield stability map (YS) for each state in the stable high (HS) and stable low (LS) zones, when compared to yield monitor data

Full size image

Iteratively including successive images across a season improved both the MLR and RF models prediction accuracy (Fig. 4). The MLR model prediction accuracy increased dramatically, 6%, with the addition of the Aug 1st image, but only increased moderately, 1–3%, with the addition of other images. The accuracy of the RF model increased ~ 13% between mid-June and early July, then gradually increased with every additional image until September. The YS estimate performance was equivalent to the RF method built with images from June to August and accuracies were similar across the four states in any given year.

Fig. 4
figure 4

Accuracy assessment for maize yield level prediction models for all fields. Each dates represents a model output built using all images from June 1 to the date

Full size image

Variability in optimal imagery dates and impacts on yield prediction –sub field sample

Sub field clustering of GCVI timeseries identifies spatial and temporal in-field variability in crop development. Demonstrating this variability on four randomly selected fields for year 2018, the time series of the mean GCVI by cluster identified different phenological timing (i.e. date of maximum GVCI) and overall crop performance (i.e. magnitude of maximum GVCI) (Fig. 5). For example, in field 3 clusters 1 and 2 achieved max GCVI by 8/1 while peak GCVI was observed by 7/15 for clusters 6 and 8. Both the shape and variability between clusters vary by field (i.e. field 4 had much larger variability between clusters than fields 1–3). Clustering GCVI time series captures in-field variability while avoiding the noise and complexity of pixel-by-pixel comparison. We calculated the dunn index to quantify the compactness and variance of the clustering, which indicates the efficiency of time-series clustering algorithm. A higher Dunn index indicates compact, well-separated clusters while a lower index indicates less well-separated clusters. The dunn index for each field-year ranges from 0.06 to 0.8 with an average index of 0.19, showing that time-series K-means and the number of predefined clusters (8) are capable of separating VI curves with different characteristics.

For cluster level analysis across all fields, the optimal date for predicting yield level with a single image varied both between and within fields. Between fields the optimal image acquisition date range was large, spanning mid-July to early September. Between all fields and clusters this range increased by approximately a month and ranged from early July to mid-September. The standard deviation of the optimal image date for all fields was 3.8 weeks. The standard deviation of subfield cluster optimal images to the field average was 2.4 weeks (Fig. 6). The largest proportion of fields (29%) had an optimal image from August 1st, two weeks before the optimal date of the full dataset. Nevertheless, the optimal images for 79% of fields were between July 15th and August 15th, with an additional 13% of the fields on September 1st. At the subfield scale, 35% of subfield areas had the same optimal image date as the field.

Field and cluster level single-image optimal date improved LR yield level prediction; the single-image optimal date of the full dataset did not. The LR for image dates July 15th, August 1st, and August 15th produced accuracies of 54%, 59%, 60%. The LR models by field and cluster optimal date produced a higher accuracy of 77% and 80%, respectively. This indicated that optimal images are spatially explicit, and determination of optimal image date is important, but the gains from selecting optimum image date by sub field cluster are small.

Fig. 5
figure 5

Time-series of GCVI at cluster level for example fields in 2018

Full size image
Fig. 6
figure 6

Weekly difference between the field level optimal image and that determined by the whole dataset (a) and between the field and subfield level optimal image (b)

Full size image

Rate of change analysis

Across successive image pairs GCVI levels were predominantly stable. On average, 77% of the subfield areas were classified as ‘Non Change’ (HH or LL) across successive image periods. These changes in the spatial pattern of GCVI were more common at the beginning and end of the season. The first month of crop growth contained the largest percentage of changing areas, with 29% from June 1st to June 15th and 31.0% from June 15th to July 1st (Fig. 7). Fields became more stable in the mid-season (July and August) and the percentage of areas labelled as HL or LH during this time decreased to 20%. In September, this percentage remained generally low, but increased to 25% between September 15th to September 30th in MI and September 1st to September 15th in MN.

Fig. 7
figure 7

Percentage of ‘Change’ (HL or LH) and ‘Non change’ (HH or LL) classes during each successive image period for each state

Full size image

Individual periods between successive images did not always provide explanatory value for end of the season yields. Most of the transient HL and LH values, happened in either early or late season, stabilized in the final yield as shown by historical yield stability zones (Fig. 8). We did not observe substantial differences among stability zones, except from 6/15 to 7/1 in MI and MN. In the high yielding area, the percentage of HL increased from less than 10% in the mid-season to 18.6% at the end of the growing season, but the decrease in GCVI did not correlate to reduced yields. In contrast, in the low yielding area, the largest percentage of the LH occurred during the beginning of the growing season. On average, 20% of the ultimately low yielding areas experienced good crop growth during the early period but still produced low yields. In general, trends in ROC dynamics were consistent across yield stability zones and any single period did not provide predictive power for end of season yields.

Fig. 8
figure 8

Percentage of HL class in the high yielding area (a) and LH class in the low yielding area (b) by the stability zone for each state

Full size image

Discussion

In this work, we compare yield stability zones analysis to remote sensing prediction of yield level to assess how well we can observe with-in-season changes in yield level. A single remotely sensing image provided limited additional yield level information, while yield stability zone analysis provides a robust estimation of yield levels for HS and LS areas- matching and often exceeding high-resolution remote sensing. The optimum sign image date occurred relatively late in the growing season and was unstable both within a field and between field. The MLR and RF models built on remote image time series provide stronger predictions of yield but are very dependent on training data. The image-to-image period rate of change of GCVI provided limited additional information. This work demonstrates the challenge of overcoming the noise in remote sensing imagery for yield prediction that is averaged out through long-term yield stability analysis and suggests in season variability due to climate remains a thorny issue.

The Integration of yield stability maps and remote sensing is needed for in-season assessment of yield levels and tactical management. Yield stability zone analysis of yield level explained the majority of yield level variation for each field and each year. These yield levels provide substantial information to inform spatially explicit baseline strategic management (i.e. variable rate nitrogen applications (Basso et al., 2019), cultivar selection, land retirement, adaptive grazing (Peterson et al., 2019), etc.). This result supports the use of yield stability analysis for the delineation precision management zones (Maestrini & Basso, 2018). However, there are limitations of YS; this approach does not provide an in-season estimate of absolute, field average yield values nor the yield level of unstable zones needed for in-season tactical management. Yield history explains about 50% of variability in stable zones but less than 30% of yield history in unstable zones (Martinez-Feria & Basso, 2020).

In-season GCVI images provide valuable information for both of these limitations of YS, but our results show yield predictions is challenged by spatial and temporal variability. The timing of peak GCVI differs among and within fields due to different crop emergence dates or growth conditions. The large range in optimal dates suggests historical best image dates or sufficient image series for YS should be established for each field or region (Fig. 6). We achieved high levels of accuracy with multi-image MLR and RF models similar to other’s findings (Ruß, 2009). Even lower prediction errors may be possible with deep learning of convolutional neural networks (Jiang et al., 2020). These empirical models may not transfer to new regions where the physical mechanisms controlling yield and phenological timing differ (Stock et al., 2023).

Even with optimal images selected by field, remote sensing is challenged to increase yield estimation accuracy above 50% at the sub field scale (Burke & Lobell, 2017). Shuai and Basso (2022) clearly demonstrated the need to account for water stress during the grain filling period to better predict subfield maize yields. Furthermore, given the availability, acquisition, and processing time, optimal dates in late July and August in the Midwest may be too late into the season to enable use of adaptive tactical management; timely in season assessment of UN areas remains a challenge. Uncertainty in the historic yield stability zones captured in unstable zones present the opportunity for powerful synergy with remote sensing integration, particularly when informed or augmented by physically based model output (Jeffries et al. 2020, Shuai & Basso, 2022).

Our analysis of short-term growth dynamic provided an unreliable approach for yield levels (Fig. 8). Better (or worse) growth during early growth periods did not always result in high (or low) crop yield; with zones commonly reverting to their historic yield zones. Alternatively, the images may simply be capturing non-crop, weedy plant growth (Marino, 2023). There are both experimental and biophysical reasons for the limited information from rate of change analysis. Experimentally, the effect of any one period change is likely not great enough to predict the change in yield level at the end of the season. Future analysis using a time series approach to characterize the rate of change and continuous yield predictions is needed. Biophysically, the high yielding areas where there is a large percentage of HL break from the expectation of the green leaf area duration from anthesis to harvest positively correlating with final grain dry weight (Wolfe et al., 1988). The ultimate, mechanistic drivers of final yield at subfield scale are primarily related to water stress (Martinez-Feria & Basso, 2020; Shuai & Basso, 2022). Furthermore, light saturation in images does not allow vegetation indices to capture assimilates translocation into grains (Karthikeyan et al., 2020). High biomass during the season does not automatically translate into higher yield, as shown in this study, these areas are in High and Stable Zones, which are normally characterized by high soil water availability and better soil health/ higher soil organic matter as the result of greater yields returning more crop residue to the soils overtime (Pan et al., 2009).To improve management, automated integration of yield stability zones, remote sensing models, and physically based models to estimate yield and nitrogen use efficiency at the subfield scale is needed (Maestrini et al., 2022). The combination of modeling imagery will require spatially explicit parameterization of the variables driving both long term average yield stability and with-in season yield variability. Future advancement in remote sensing products, with higher return periods and hyperspectral imaging will increase the number of usable features for yield prediction and our ability to relate phenology to yield. Advancing machine learning algorithms will help us identify meaningful patterns from these rich data. Physically based models will conduct the hypothesis testing to relate these statical relationships to physical mechanisms on the landscape. Adapting these tools, here in development, for relatively simple cropping systems, to diverse and highly integrated systems is needed to reduce environmental impact while sustaining exceptional yields (Hatfield et al., 2020). We encourage researchers to further this analysis to consider more crop types or medium-resolution satellite images (such as Landsat) across a broader region to get more generalized conclusions. Learning from simple systems across a continuum of soil and climatic conditions may enable generalizable strategic and tactical approaches to precision management across the agroecological system.

Conclusion

This study evaluated in-season and real time methods to understand spatial crop yield zone variance. We found that:1) Historic yield stability maps provide estimation of corn yield levels in stable zones with accuracy greater than in-season remote sensing, which could be used to support strategic fertilizer management. 2) Use of field and sub-field cluster specific imagery improved remote sensing estimations of yield and yield level. 3) Remotely sensed plant growth dynamics in any given month, while more variable at the beginning and the end of the season, did not have clear impacts on final yield or yield level. The findings of this research provide valuable insights on the need for an integrated geospatial systems approach based on historical yield zones, in-season remote sensing, and physical modeling to better capture the interactions and dynamic feedbacks to better estimate in-season yield predictions and to inform tactical management.