Assessing within-Field Corn and Soybean Yield Variability from WorldView-3, Planet, Sentinel-2, and Landsat 8 Satellite Imagery

Crop yield monitoring is an important component in agricultural assessment. Multi-spectral remote sensing instruments onboard space-borne platforms such as Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer (MODIS), and Visible Infrared Imaging Radiometer Suite (VIIRS) have shown to be useful for efficiently generating timely and synoptic information on the yield status of crops across regional levels. However, the coarse spatial resolution data inherent to these sensors provides little utility at the management level. Recent satellite imagery collection advances toward finer spatial resolution (down to 1 m) alongside increased observational cadence (near daily) implies information on crops obtainable at field and within-field scales to support farming needs is now possible. To test this premise, we focus on assessing the efficiency of multiple satellite sensors, namely WorldView-3, Planet/Dove-Classic, Sentinel-2, and Landsat 8 (through Harmonized Landsat Sentinel-2 (HLS)), and investigate their spatial, spectral (surface reflectance (SR) and vegetation indices (VIs)), and temporal characteristics to estimate corn and soybean yields at sub-field scales within study sites in the US state of Iowa. Precision yield data as referenced to combine harvesters’ GPS systems were used for validation. We show that imagery spatial resolution of 3 m is critical to explaining 100% of the within-field yield variability for corn and soybean. Our simulation results show that moving to coarser resolution data of 10 m, 20 m, and 30 m reduced the explained variability to 86%, 72%, and 59%, respectively. We show that the most important spectral bands explaining yield variability were green (0.560 μm), red-edge (0.726 μm), and near-infrared (NIR − 0.865 μm). Furthermore, the high temporal frequency of Planet and a combination of Sentinel-2/Landsat 8 (HLS) data allowed for optimal date selection for yield map generation. Overall, we observed mixed performance of satellite-derived models with the coefficient of determination (R2) varying from 0.21 to 0.88 (averaging 0.56) for the 30 m HLS and from 0.09 to 0.77 (averaging 0.30) for 3 m Planet. R2 was lower for fields with higher yields, suggesting saturation of the satellite-collected reflectance features in those cases. Therefore, other biophysical variables, such as soil moisture and evapotranspiration, at similar fine spatial resolutions are likely needed alongside the optical imagery to fully explain the yields.


Introduction
Crop yield assessment and forecasting is one of the major components in agriculture production monitoring [1]. Timely, accurate, and spatially explicit information on crop yields at global, national, and regional scales is extremely important [2] for understanding and allocating the food supply. Furthermore, spatial variability of crop yields at the field scale can help provide objective information to farmers to improve management practices and identify yield gaps [3][4][5] or to improve insurance companies' risk models [6].
The increased number of very high spatial resolution (VHR, <5 m) sensors onboard space-borne platforms, including micro and nanosatellites, are providing opportunities for a new generation of spatially detailed geo-information products for agricultural monitoring. High temporal resolution (near daily cadence) is additionally required for agricultural applications to monitor and capture the dynamics of crop growth [7][8][9][10]. Images acquired by the WorldView-3 (WV-3,~1.25 m) satellite or the PlanetScope (Dove-Classic,~3 m) constellation of satellites have finer spatial resolution compared to images provided by moderate spatial resolution satellites such as Landsat 8 (30 m) [11] and Sentinel-2 (10-20 m) [12]. While WV-3 provides better spatial resolution and richer spectral information in visible and near-infrared (NIR) spectrum (eight bands vs. four) compared to Planet/Dove, only the latter provides a near-daily revisit cycle.
The combined spatial and temporal attributes from WorldView-3 and Planet can provide new insight on the variability of crop parameters at field scale, including yields [13]. Satellite data availability during the growing season can be used to predict before harvest within-the-field crop yields or analyze areas of higher and lower yields for precision farming. It can also be used for stratification purposes within the sampling-based estimation of yields [14,15]. The latter is especially important, since if satellite data could capture within-field crop variability during the vegetation season, they would allow one to guide, and potentially improve, the sampling strategies for ground-based crop yield estimation and validation. This is similar to the approach of using satellite-derived classification maps for guiding sampling strategies to estimate areas of land cover land use change (LCLUC) [16,17]. However, in the case of LCLUC, one must deal with categorical information, while in the case of yields, it is continuous information. To address that aspect, dedicated frameworks such as VALidation of European Remote sensing Instruments (VALERI) [18] protocol have been setup to validate satellite-derived maps of biophysical variables from coarse to moderate spatial resolutions. The VALERI protocol introduces an elementary sampling unit (ESU) corresponding to the spatial resolution of a satellite-derived map (e.g., a 30 × 30 m pixel) and a number of point-based measurements (e.g., 10-15 samples per ESU). Thus, the variability of the biophysical parameter inside the ESU can be captured.
Previous studies have attempted to connect remote sensing data (both reflectances and vegetation indices (VIs)) acquired from space-borne VHR sensors with crop yields at field scale and showed promising results [13,[19][20][21][22][23]. Most studies used IKONOS, Quickbird, SPOT5, or the RapidEye satellites. Even through VHR, each would only provide one or a few image acquisitions throughout the season, which would bring large uncertainties in terms of optimal timing for correlating satellite-derived features with yields.
Therefore, space-borne sensors at very high spatial and temporal resolution can provide new insights and opportunities for precision agriculture [13,24] and further expand applications of remote sensing sensors in agricultural domains at various geographical scales. The present study takes advantage of VHR satellite data made available within the National Aeronautics and Space Administration's (NASA) program to explore the strengths and the limitations of multi-spectral sensors to explain crop yield variability at field scale. Specifically, the present study aims to address the following research questions: • Q3: What is the optimal timing of satellite image acquisitions and optimal spectral bands for explaining within-field variability?
The study was performed for 2018 and 2019 using VHR data (1-3 m), acquired by PlanetScope (Dove-Classic) and WorldView-3 satellites. Additionally, use comparisons were made to moderate spatial resolution data (30 m) as acquired by Landsat-8 and Sentinel-2 satellites. The crops of corn and soybean were both analyzed with the study sites being within Iowa, US. The within-field ground refence data were precision collected by harvesting machines.

Study Area and Reference Yield Data
Corn and soybean yields were obtained for 30 fields in 2018 and 2019 from Iowa State University. Fields were scattered throughout 7 counties in the state of Iowa, namely Boone, Buena Vista, Hamilton, Jasper, Lucas, Pottawatomie, and Story ( Figure 1). The detailed characteristics of each field are described in Table A1. The yield data were collected by combine harvesting machines, which provide dry volumetric yields roughly every three meters along track. These yields were pre-processed before correlating them to the satellite imagery. Headlands, row ends, and samples outside ±3 standard deviations from the field mean were removed. In cases when several data records with different yield values were within 3 m of each other, only the data record with the maximum yield was retained, while others were filtered out [20,25,26].
The study was performed for 2018 and 2019 using VHR data (1-3 m), a PlanetScope (Dove-Classic) and WorldView-3 satellites. Additionally, use co were made to moderate spatial resolution data (30 m) as acquired by Landsattinel-2 satellites. The crops of corn and soybean were both analyzed with the being within Iowa, US. The within-field ground refence data were precision c harvesting machines.

Study Area and Reference Yield Data
Corn and soybean yields were obtained for 30 fields in 2018 and 2019 from University. Fields were scattered throughout 7 counties in the state of Iow Boone, Buena Vista, Hamilton, Jasper, Lucas, Pottawatomie, and Story (Figure tailed characteristics of each field are described in Table A1. The yield data wer by combine harvesting machines, which provide dry volumetric yields roug three meters along track. These yields were pre-processed before correlating t satellite imagery. Headlands, row ends, and samples outside ±3 standard devia the field mean were removed. In cases when several data records with differen ues were within 3 m of each other, only the data record with the maximum retained, while others were filtered out [20,25,26]. Average per-field yield values as well as planting and harvest dates are Figure 2. The average yield for corn was 12.1 t/ha and varied from 9.6 t/ha to 13 for soybean was 3.9 t/ha and varied from 2.4 t/ha to 4  Average per-field yield values as well as planting and harvest dates are shown in Figure 2. The average yield for corn was 12.1 t/ha and varied from 9.6 t/ha to 13.6 t/ha and for soybean was 3.9 t/ha and varied from 2.4 t/ha to 4.5 t/ha. On average, corn and soybean were planted on 4 May (day of year (DOY) 124) and 23 May (DOY 143). They

Satellite Data at Very High Spatial Resolution (VHR): WorldView-3 and Planet
WorldView-3 (WV-3) data were provided by the NASA's Commercial Archive Data for NASA investigators (cad4nasa.gsfc.nasa.gov) [27]. WV-3 acquires data of the Earth's surface in eight spectral bands in the visible and infrared range of electromagnetic spectrum: coastal blue (0.426 μm), blue (0.479 μm), green (0.552 μm), yellow (0.610 μm), red (0.662 μm), red-edge (0.726 μm), and two near-infrared (NIR1 0.831 μm and NIR2 0.910 μm) channels. At-nadir spatial resolution is 1.25 m. WV-3 also acquires data in eight short-     [27]. WV-3 acquires data of the Earth's surface in eight spectral bands in the visible and infrared range of electromagnetic spectrum: coastal blue (0.426 µm), blue (0.479 µm), green (0.552 µm), yellow (0.610 µm), red (0.662 µm), red-edge (0.726 µm), and two near-infrared (NIR1 0.831 µm and NIR2 0.910 µm) channels. At-nadir spatial resolution is 1.25 m. WV-3 also acquires data in eight short-wave infrared bands (ranging from 1195 nm to 2365 nm), but those were unavailable and, hence, not used. Since only archival images were selected and no new acquisitions were available, we were able to obtain two WV-3 images over the Coles fields, for which within-fields yields were available ( Figure 4). Though the timing of the WV-3 images (2 and 21 July 2018) might not have been optimal, those data were still useful in analyzing WV-3 spectral bands at 1.25 m spatial resolution to explain within-field corn and soybean yield variability early in the season (planting dates were 16 May and 6 June 2018 for corn and soybean, respectively). Images' digital numbers were converted to the topof-atmosphere (TOA) reflectance values and further processed for atmospheric correction. We used an adapted version of the Land Surface Reflectance Code (LaSRC), previously developed for the moderate spatial resolution satellites such as Landsat 8 [28] and Sentinel-2 [29]. LaSRC is a generic atmospheric correction algorithm for estimating land surface reflectance (SR), which takes into account absorption by atmospheric gases and scattering by molecules and aerosols. It is based on the inversion of the "6SV" radiative transfer code. The result of the atmospheric correction is the estimate of SR in the corresponding spectral bands of WV-3 satellite. wave infrared bands (ranging from 1195 nm to 2365 nm), but those were unavailable and, hence, not used. Since only archival images were selected and no new acquisitions were available, we were able to obtain two WV-3 images over the Coles fields, for which withinfields yields were available ( Figure 4). Though the timing of the WV-3 images (2 and 21 July 2018) might not have been optimal, those data were still useful in analyzing WV-3 spectral bands at 1.25 m spatial resolution to explain within-field corn and soybean yield variability early in the season (planting dates were 16 May and 6 June 2018 for corn and soybean, respectively). Images' digital numbers were converted to the top-of-atmosphere (TOA) reflectance values and further processed for atmospheric correction. We used an adapted version of the Land Surface Reflectance Code (LaSRC), previously developed for the moderate spatial resolution satellites such as Landsat 8 [28] and Sentinel-2 [29]. LaSRC is a generic atmospheric correction algorithm for estimating land surface reflectance (SR), which takes into account absorption by atmospheric gases and scattering by molecules and aerosols. It is based on the inversion of the "6SV" radiative transfer code. The result of the atmospheric correction is the estimate of SR in the corresponding spectral bands of WV-3 satellite. Planet data, acquired by Dove-Classic sensors, were acquired through the NASA program [30]. PlanetScope is composed of a constellation of satellites that provide a neardaily global coverage of the Earth surface [31]. Dove-Classic sensors acquire data in the four spectral bands: blue (0.485 µ m), green (0.545 µ m), red (0.630 µ m), and near-infrared (0.820 µ m) at 3.25 m spatial resolution. We acquired all cloud-free Planet data over the fields for the whole growing season in 2018 and 2019 and an SR product provided by Planet. There were 1167 and 1095 unique scenes, respectively, for those two years.

Satellite Data at Moderate Spatial Resolution: Landsat 8 and Sentinel-2
We used NASA's Harmonized Landsat Sentinel-2 (HLS, ver. 1.4) product at 30 m [32] to correlate satellite-derived features with crop yields. HLS is an analysis ready data product, which is derived from original TOA reflectance data by applying a series of pre-processing steps, including co-registration, cloud screening, atmospheric and angular corrections, and spectral adjustments. This ensures that data acquired by Landsat 8 and Sentinel-2 satellites are geometrically and radiometrically harmonized and can be used together. Combination of Landsat 8 and Sentinel-2A/B satellites provides a global median average revisit interval of 2.9 days [33] and already proved to be beneficial for crop yield assessment [9]. HLS is organized using a Sentinel-2 tiling system with each tile covering the area Planet data, acquired by Dove-Classic sensors, were acquired through the NASA program [30]. PlanetScope is composed of a constellation of satellites that provide a neardaily global coverage of the Earth surface [31]. Dove-Classic sensors acquire data in the four spectral bands: blue (0.485 µm), green (0.545 µm), red (0.630 µm), and near-infrared (0.820 µm) at 3.25 m spatial resolution. We acquired all cloud-free Planet data over the fields for the whole growing season in 2018 and 2019 and an SR product provided by Planet. There were 1167 and 1095 unique scenes, respectively, for those two years.

Satellite Data at Moderate Spatial Resolution: Landsat 8 and Sentinel-2
We used NASA's Harmonized Landsat Sentinel-2 (HLS, ver. 1.4) product at 30 m [32] to correlate satellite-derived features with crop yields. HLS is an analysis ready data product, which is derived from original TOA reflectance data by applying a series of pre-processing steps, including co-registration, cloud screening, atmospheric and angular corrections, and spectral adjustments. This ensures that data acquired by Landsat 8 and Sentinel-2 satellites are geometrically and radiometrically harmonized and can be used together. Combination of Landsat 8 and Sentinel-2A/B satellites provides a global median average revisit interval of 2.9 days [33] and already proved to be beneficial for crop yield assessment [9]. HLS is organized using a Sentinel-2 tiling system with each tile covering the area of approximately 110 km × 110 km. We used the following tiles, which covered corn and soybean fields in the study region: 15TVH, 15TVG, 15TVF, 15TUH, and 15TUF. We used the following spectral bands: green (wavelength~0.560 µm), red (~0.660 µm), NIR (~0.865 µm from the 8A band of Sentinel-2), and two short-wavelength infrared (SWIR) (~1.5 and~2.2 µm). Overall, 407 and 400 unique HLS scenes were used in 2018 and 2019, respectively.

Methods
First, we explored the effect spatial resolution has on capturing field-scale yield variability. For this, the harvest data from machinery were used to simulate how much yield variability is reduced when spatial resolution decreases. Yield data from the combine machinery were again available approximately every 3 m. Yield maps were generated at spatial resolutions of 10 m, 20 m, and 30 m by averaging all the point yield data within the coarser pixel. Relative efficiency, which is the ratio between yield variance at coarse resolution (10 m, 20 m, or 30 m) to yield variance at 3 m resolution, was used as a numerical metric for yield variability evaluation: where R(x) is the relative efficiency at spatial resolution x m, V y (x) is the variance of yield at x m, and x 0 = 3 m. This metric is similar to the metric used in land area estimation and in quantifying the efficiency of satellite data in reducing the variance of the area estimates [16,34]. We also used the coefficient of variation (CV) to measure a relative variability of yields within the pixel or a field: where σ y and µ y are standard deviation and average of yields y.
In order to explore the yield variability at field scale, the reference within-field corn and soybean yield values were correlated with the satellite-derived features of SR or VIs using a multivariable linear regression (Table 1). For 1.25 m WV-3 and 3 m Planet data, those features were derived for each yield sample separately, while for 30 m HLS data, those yield values were averaged within each 30 m pixel. Specifically, for each cloud-free satellite acquisition, a linear regression model was built between yields and SRs or VIs, and an adjusted coefficient of determination R 2 adj and a root mean square error (RMSE) were calculated: where y est and y re f are estimate and reference yield values, respectively; y re f = 1/n ∑ n i=1 y re f is the average of the reference values, n is the number of samples, and p is the number of variables (without the constant term). Table 1. Vegetation indices used in the study.

Vegetation Index Formula
Normalized difference vegetation index (NDVI) For Planet and HLS data, the best date in terms of maximum R 2 adj was determined to correlate with crop yields. For WV-3, no time series was available, thus we focused on exploring various spectral bands using a "feature importance" metric derived from the random forest regression [35]. This was done through multiple random permutations of input features, and regression was evaluated in terms of information gain. Features that provided higher informativeness had a higher importance level. Figure 5 shows the effect of spatial resolution on capturing within-field yield variability using a corn field as an example. These results show the difference between high spatial resolution data compared to moderate spatial resolution of 10-30 m. For example, only 59% of yield variability can be explained with 30 m, while the rest is lost because of coarsening and mixture effects inside the 30 m pixel.

Effect of Spatial Resolution on Capturing Field Scale Yield Variability
At present, one of the popular approaches for collecting reference crop yields is through the crop cuts [14,15]. At fine spatial scales (usually at less than 1 m), destructive harvesting of the crops is made, and the yield is estimated from the collected grain. Those yields are then used as reference and compared to satellite-derived yields at moderate spatial resolution. However, sampling within the moderate resolution pixel should be statistically rigorous, thus the selected samples would capture yield variability across the moderate resolution pixel, and therefore the estimates would be unbiased. Figure 6 shows a histogram of the CV (Equation (2)) for the yields at 10 m, 20 m, and 30 m resolutions. CV exceeds 20% and 10% at 30 m and 10 m, respectively. Our results are in line with the previously developed VALERI protocol [18]. Therefore, a similar protocol shall be implemented for any sample-based ground measurements of crop yields, when multiple samples should be selected inside a 10 or 30 m pixel to capture yield variability. The VALERI protocol recommends 10-15 samples of in-situ ground measurements of the biophysical variable inside the pixel. One cannot make a single crop cut inside a 30 m pixel and compare to the satellite-derived yields at moderate spatial resolution.       Table 2 shows values of the coefficient of determination, RMSE values obtained for the linear models correlating corn and soybean yields, and SR values derived from WV-3 images. A model based on WV-3 image acquired on 21 July 2018 was able to explain 48% of corn (RMSE = 1.07 t/ha or 9.2%) and 36% of soybean (RMSE = 0.47 t/ha or 13.2%) yield variability. This was 66 and 45 days after corn and soybean planting, respectively. The SRbased model outperformed a single-regressor model based on the VI. For corn, the best one was with the red edge chlorophyll index (Cire) and gave R 2 = 0.44 and RMSE = 1.11 t/ha (9.5%) (Figure 7). For soybean, the best VI-based model was with the difference vegetation index (DVI) and gave R 2 = 0.25 and RMSE = 0.50 t/ha (14.3%). Overall, VIs based on the red spectral band were not able to explain corn yield variability (the best was Wide dynamic range vegetation index (WDRVI) at 19%), whereas the green-based index Green chlorophyll vegetation index (GCVI) and the red-edge-based index CIre were able to explain 37% and 44% of yield variability, respectively.  (5)); RRMSE: relative RMSE (Equation (6)).

Analysis of Corn and Soybean Yield Varibility with WV-3 Data
Sens. 2021, 13, x FOR PEER REVIEW 9 of 18 Table 2 shows values of the coefficient of determination, RMSE values obtained for the linear models correlating corn and soybean yields, and SR values derived from WV-3 images. A model based on WV-3 image acquired on 21 July 2018 was able to explain 48% of corn (RMSE = 1.07 t/ha or 9.2%) and 36% of soybean (RMSE = 0.47 t/ha or 13.2%) yield variability. This was 66 and 45 days after corn and soybean planting, respectively. The SRbased model outperformed a single-regressor model based on the VI. For corn, the best one was with the red edge chlorophyll index (Cire) and gave R 2 = 0.44 and RMSE = 1.11 t/ha (9.5%) (Figure 7). For soybean, the best VI-based model was with the difference vegetation index (DVI) and gave R 2 = 0.25 and RMSE = 0.50 t/ha (14.3%). Overall, VIs based on the red spectral band were not able to explain corn yield variability (the best was Wide dynamic range vegetation index (WDRVI) at 19%), whereas the green-based index Green chlorophyll vegetation index (GCVI) and the red-edge-based index CIre were able to explain 37% and 44% of yield variability, respectively.  (5)); RRMSE: relative RMSE (Equation (6)).  Not all spectral bands were equally important in explaining corn and soybean variability. Figure 8 shows performance of various models and relative importance of spectral bands. Results suggest that, during the stage of corn development when a WV-3 image was taken, in this case, 66 days after planting, the red-edge (0.726 µm), the green (0.552 µm), and the NIR (0.831 µm) bands were the most important in explaining yield variance at field scale. These are in accordance with results obtained for VIs, when red-edge and green-based indices outperformed red-based indices.

Analysis of Corn and Soybean Yield Varibility with WV-3 Data
x FOR PEER REVIEW 10 of 18 Not all spectral bands were equally important in explaining corn and soybean variability. Figure 8 shows performance of various models and relative importance of spectral bands. Results suggest that, during the stage of corn development when a WV-3 image was taken, in this case, 66 days after planting, the red-edge (0.726 µ m), the green (0.552 µ m), and the NIR (0.831 µ m) bands were the most important in explaining yield variance at field scale. These are in accordance with results obtained for VIs, when red-edge and green-based indices outperformed red-based indices.

Analysis of Corn and Soybean Yield Varibility with Planet and HLS Data
Unlike with WV-3, Planet and HLS provided a dense time series of multi-spectral images, which allowed the determination of R 2 optimized timing to correlate satellite-derived features with yields. Figure 9 shows examples of how R 2 changes with time for some of the corn and the soybean fields in 2018 and 2019. Results show that, for some fields, a clear peak was visible, and R 2 reached values up to 0.8, suggesting possibility to explain up to 80% of within-field yield variability, whereas, for other fields, the maximum seasonal R 2 was around 0.1-0.2, suggesting that multi-spectral HLS and Planet data failed to explain within-field yield variability. Table 3 shows DOY values for which maximal was R 2 achieved and its correspondence to the harvesting date. The maximal correlation was observed 90-100 days (3 months) before corn harvest and 70-80 days (2-2.5 months) before soybean harvest.

Crop
Year Harvest HLS Planet Figure 8. Relative importance of WV-3 spectral bands when correlating them with corn yields. WV-3 was acquired on 21 July 2018. Shown also is R 2 for the linear model between SR values and corn yields.

Analysis of Corn and Soybean Yield Varibility with Planet and HLS Data
Unlike with WV-3, Planet and HLS provided a dense time series of multi-spectral images, which allowed the determination of R 2 optimized timing to correlate satellitederived features with yields. Figure 9 shows examples of how R 2 changes with time for some of the corn and the soybean fields in 2018 and 2019. Results show that, for some fields, a clear peak was visible, and R 2 reached values up to 0.8, suggesting possibility to explain up to 80% of within-field yield variability, whereas, for other fields, the maximum seasonal R 2 was around 0.1-0.2, suggesting that multi-spectral HLS and Planet data failed to explain within-field yield variability. Table 3 shows DOY values for which maximal was R 2 achieved and its correspondence to the harvesting date. The maximal correlation was observed 90-100 days (3 months) before corn harvest and 70-80 days (2-2.5 months) before soybean harvest. prominent, 0.32 ± 0.12 (RMSE = 1.92 ± 0.56 t/ha) compared to 0.29 ± 0.13 (RMSE = 1.98 ± 0.57 t/ha) for corn and 0.29 ± 0.20 (RMSE = 0.79 ± 0.41 t/ha) compared to 0.30 ± 0.23 (RMSE = 0.80 ± 0.41 t/ha) for soybean. Overall, R 2 values for corn were higher than those for soybean. Addition of SWIR bands to blue, green, red, and NIR bands in the HLS-based models slightly improved R 2 performance from 0.59 ± 0.17 and 0.62 ± 0.16 for corn and from 0.47 ± 0.20 and 0.50 ± 0.20 for soybean.  Overall, the variable performance in terms of R 2 was observed when using SR and VI satellite-derived features to explain corn and soybean yield variability ( Figure 10). Best results for HLS and Planet were obtained using a linear model incorporating all spectral bands. This outperformed models utilizing vegetation indices. For HLS, the SR-based models gave R 2 = 0.62 ± 0.16 (RMSE = 0.95 ± 0.34 t/ha) and 0.50 ± 0.20 (RMSE = 0.39 ± 0.15 t/ha) for corn and soybean, respectively. Average and standard deviation values of R 2 and RMSE were calculated for all corn and soybean fields in 2018 and 2019. The best VIs models were those based on green and NIR spectral bands such as the WV-3 case. Both GCVI and GNDVI gave R 2 = 0.48 ± 0.20 (RMSE = 1.07 ± 0.47 t/ha) and 0.35 ± 0.22 (RMSE = 0.46 ± 0.17 t/ha) for corn and soybean, respectively. Similar results were achieved for Planet, though the difference between the SR-based model and the best VIs models was less prominent, 0.32 ± 0.12 (RMSE = 1.92 ± 0.56 t/ha) compared to 0.29 ± 0.13 (RMSE = 1.98 ± 0.57 t/ha) for corn and 0.29 ± 0.20 (RMSE = 0.79 ± 0.41 t/ha) compared to 0.30 ± 0.23 (RMSE = 0.80 ± 0.41 t/ha) for soybean. Overall, R 2 values for corn were higher than those for soybean. Addition of SWIR bands to blue, green, red, and NIR bands in the HLS-based models slightly improved R 2 performance from 0.59 ± 0.17 and 0.62 ± 0.16 for corn and from 0.47 ± 0.20 and 0.50 ± 0.20 for soybean. R 2 values obtained for HLS were higher than those obtained for Planet because of different scales employed (3 m vs. 30 m, or an order of 10 × 10 = 100). Averaging 3 m yield values to the 30 m HLS resolution reduces the variance of yield ( Figure 5) that could be explained with satellite data, therefore increasing the coefficient of determination. Overall, a strong correlation between HLS-and Planet-based models was observed (Figure 11), suggesting consistency of results achieved with HLS and Planet data. In other words, the same set of fields gave higher or lower R 2 values for both Planet and HLS.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 18 tations in using multi-spectral satellite imagery in explaining within-field corn and soybean yield variability, especially for fields with large yield values. We did not find strong correlations between R 2 and yield variance or CV. Figure 11. Comparison of R 2 for best models between corn and soybean yields and SR values obtained with HLS and Planet data.
(a) (b) Figure 11. Comparison of R 2 for best models between corn and soybean yields and SR values obtained with HLS and Planet data.
We found a strong correlation between R 2 and field-averaged corn and soybean yields ( Figure 12). As the yield values increased, the performance of HLS-and Planet-based models to explain within-field yield variability decreased. Such dependence was stronger for soybean fields than for corn. R 2 values were less than 0.5 for corn and soybean fields with yields larger than 12.5 t/ha and 3.5 t/ha, respectively. These results show limitations in using multi-spectral satellite imagery in explaining within-field corn and soybean yield variability, especially for fields with large yield values. We did not find strong correlations between R 2 and yield variance or CV.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 18 tations in using multi-spectral satellite imagery in explaining within-field corn and soybean yield variability, especially for fields with large yield values. We did not find strong correlations between R 2 and yield variance or CV.

Discussion and Conclusions
This study aimed at investigating the efficiency of multi-spectral satellite imagery at various spatial and temporal resolutions to explain within-field yield variability for corn and soybean crops. First, we explored the effect of spatial resolution on yield variability. We found that reducing spatial resolution from 3 m to 30 m reduced yield variance at field scale with relative efficiency (Equation 1) decreasing from 1.0 to 0.6 ( Figure 5), therefore decreasing capabilities of capturing within-field yields with moderate spatial resolution sensors. At the same time, the variance of yields inside 10, 20, or 30 m pixels increased up to 20% ( Figure 6). These results are consistent with previous recommendations established within the VALERI protocol [18] to validate coarse and moderate resolution remote sensing products and have implications on collecting ground-based crop yields, for example, through crop cuts [14,15]. When validating coarse (100 m-1 km) and moderate (10-30 m) resolution crop yield maps, one must ensure that variability inside the pixel is captured through multiple (10-15) point-based samples.
In order to detect the optimal timing of correlating satellite-derived features with corn and soybean yields, a high temporal resolution of satellite data was required. This was achieved through the use of near-daily Planet data and a combination of 3-5 day revisit cycles of Sentinel-2 and Landsat 8 combined [33]. The maximal correlations for corn and soybean were obtained at 3 and 2-2.5 months, respectively, which is consistent with previous studies utilizing coarse resolution MODIS data and county-level yields [36]. We observed a very variable performance of these datasets to explain within-field yield variability: coefficient of determination R 2 (Equations (3) and (4)) varied from 0.21 to 0.88 (average 0.56 ± 0.19) for 30 m HLS and from 0.09 to 0.77 (0.30 ± 0.16) for 3 m Planet. Obtained R 2 values were higher for HLS than Planet because of various scales employed; averaging 3 m yields to 30 m reduced within-field yield variance ( Figure 5), therefore potentially benefiting moderate resolution at the expense of details lost. R 2 values were also higher for corn than soybean, which is consistent with previous works utilizing remote sensing

Discussion and Conclusions
This study aimed at investigating the efficiency of multi-spectral satellite imagery at various spatial and temporal resolutions to explain within-field yield variability for corn and soybean crops. First, we explored the effect of spatial resolution on yield variability. We found that reducing spatial resolution from 3 m to 30 m reduced yield variance at field scale with relative efficiency (Equation 1) decreasing from 1.0 to 0.6 ( Figure 5), therefore decreasing capabilities of capturing within-field yields with moderate spatial resolution sensors. At the same time, the variance of yields inside 10, 20, or 30 m pixels increased up to 20% ( Figure 6). These results are consistent with previous recommendations established within the VALERI protocol [18] to validate coarse and moderate resolution remote sensing products and have implications on collecting ground-based crop yields, for example, through crop cuts [14,15]. When validating coarse (100 m-1 km) and moderate (10-30 m) resolution crop yield maps, one must ensure that variability inside the pixel is captured through multiple (10-15) point-based samples.
In order to detect the optimal timing of correlating satellite-derived features with corn and soybean yields, a high temporal resolution of satellite data was required. This was achieved through the use of near-daily Planet data and a combination of 3-5 day revisit cycles of Sentinel-2 and Landsat 8 combined [33]. The maximal correlations for corn and soybean were obtained at 3 and 2-2.5 months, respectively, which is consistent with previous studies utilizing coarse resolution MODIS data and county-level yields [36]. We observed a very variable performance of these datasets to explain within-field yield variability: coefficient of determination R 2 (Equations (3) and (4)) varied from 0.21 to 0.88 (average 0.56 ± 0.19) for 30 m HLS and from 0.09 to 0.77 (0.30 ± 0.16) for 3 m Planet. Obtained R 2 values were higher for HLS than Planet because of various scales employed; averaging 3 m yields to 30 m reduced within-field yield variance ( Figure 5), therefore potentially benefiting moderate resolution at the expense of details lost. R 2 values were also higher for corn than soybean, which is consistent with previous works utilizing remote sensing data and usually showing better performance for corn [7,36,37]. The maximal correlation between satellite-derived features and yields was observed 2-3 months before the harvest, therefore potentially enabling yield forecasts at field scale end of July to mid-August. We also found that performance of the satellite-derived models was better for corn and soybean fields with lower average yields: as the yields increase, R 2 values decrease. It means that multi-spectral data saturate over high yields and cannot capture yield variability. Other variables at such 3 to 30 m scale, such as soil moisture and evapotranspiration (ET) [38][39][40], irrigation vs. rainfed, should be considered, as previous studies showed that field-scale yield variability for corn and soybean in the US is mainly attributed to soil differences in the field and the water holding capacity [41].
In terms of spectral bands, green (~0.560 µm) and NIR (~0.865 µm) were the most important for Planet and HLS to explain corn and soybean within-field yield variability. We did not consider red-edge from Sentinel-2, since it is not available in Landsat 8. When analyzing WV-3 data, red-edge (0.726 µm) showed to be the most important in explaining yield variability. Overall, SR-based models outperformed VI-based models highlighting the importance of incorporating SR values directly into the yield models [9]. VIs based on NIR and red-edge or green spectral bands performed better than those based on NIR and red, as the former previously showed a better correspondence to corn and soybean productivity [42].
The further research should focus on the within-field water-stress indicators that can be retrieved using the fusion of satellite data and corresponding models and active remote sensing sensors such as synthetic aperture radar (SAR) and LIDAR.

Acknowledgments:
We would like to thank NASA for providing support in generating the Harmonized Landsat Sentinel-2 (HLS) product. We would like to thank Iowa State University for providing field-scale yield data records. DigitalGlobe/Maxar data were provided by NASA's Commercial Archive Data for NASA investigators (cad4nasa.gsfc.nasa.gov accessed on 29 January 2021) under the National Geospatial-Intelligence Agency's NextView license agreement.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1 provides detailed information on the fields used in the study, including their location, crop type, and registered yields.