To Blend or Not to Blend? A Framework for Nationwide Landsat–MODIS Data Selection for Crop Yield Prediction

The onus for monitoring crop growth from space is its ability to be applied anytime and anywhere, to produce crop yield estimates that are consistent at both the subfield scale for farming management strategies and the country level for national crop yield assessment. Historically, the requirements for satellites to successfully monitor crop growth and yield differed depending on the extent of the area being monitored. Diverging imaging capabilities can be reconciled by blending images from high-temporal-frequency (HTF) and high-spatial-resolution (HSR) sensors to produce images that possess both HTF and HSR characteristics across large areas. We evaluated the relative performance of Moderate Resolution Imaging Spectroradiometer (MODIS), Landsat, and blended imagery for crop yield estimates (2009–2015) using a carbon-turnover yield model deployed across the Australian cropping area. Based on the fraction of missing Landsat observations, we further developed a parsimonious framework to inform when and where blending is beneficial for nationwide crop yield prediction at a finer scale (i.e., the 25-m pixel resolution). Landsat provided the best yield predictions when no observations were missing, which occurred in 17% of the cropping area of Australia. Blending was preferred when <42% of Landsat observations were missing, which occurred in 33% of the cropping area of Australia. MODIS produced a lower prediction error when ≥42% of the Landsat images were missing (~50% of the cropping area). By identifying when and where blending outperforms predictions from either Landsat or MODIS, the proposed framework enables more accurate monitoring of biophysical processes and yields, while keeping computational costs low.


Introduction
The world's human population is projected to increase by more than 35% by 2050 [1]. To contribute to improved global food security, the next generation of crop models and agricultural decision support tools needs to efficiently and consistently operate across various scales [2]. Accurate nationwide crop yield forecasts may ensure food security to the citizens. More accurate crop yield prediction at the subfield scale can provide farmers with more detailed information for guiding, within the growing season, in-field variable rate applications of fertilizer, herbicides, and pesticides. An efficient Table 1. Review of studies that used data fusion approaches to estimate crop yields and other variables highly related to crop yield (e.g., growing season evapotranspiration and crop phenology). Records are sorted chronologically then alphabetically. When more than one site is reported in a single paper, these are identified with an (A) and (B) in the relevant columns as required. R 2 : the coefficient of determination; RMSE: root-mean-square error; MAE: mean absolute error; RMAE: relative MAE; NR: not reported. In Australia, dryland agriculture is practiced across the region known colloquially as the "wheatbelt". STARFM-Spatial and Temporal Adaptive Reflectance Fusion Model; ESTARFM-enhanced STARFM; MODIS-Moderate Resolution Imaging Spectroradiometer; NIR-near infrared; US-United States; NE-New England; NASA-National Aeronautics and Space Administration; L-M-Landsat/MODIS blend; HTF-high temporal frequency; HSR-high spatial resolution.

RS-Bases Model to Estimate Crop Yields
Key Results and Accuracy [28] Spatial  The daily ET derived from the blended data produced the RMAE of 19% with the observed ET (mm/day). The spatial pattern of cumulative ET corresponded to the measured yield.  The county-level correlation between observed and predicted maize yields improved from 0.47 to 0.93 when aligning the ratio of actual-to-reference ET by emergence date rather than calendar date.

Study Area
The southern Australian agricultural system is dominated by dryland agriculture where cereals (e.g., wheat, barley, and oats), oilseeds (e.g., canola), and legumes (e.g., lupins, chickpeas, field peas, and soybeans) are planted, often in rotation with annual pastures and fallows. Australian wheatbelt ( Figure 1) spans an extremely variable agroecological environment with respect to the climate across the continent, leading to the high spatial variability in Australian grain production. The precipitation varies enormously across the country ( Figure S1, Supplementary Materials); winter (June-August) precipitations are dominant in Western Australia and South Australia, while summer (December-February) precipitations are dominant in Queensland and northern New South Wales. In southern New South Wales and Victoria, total precipitation is more uniformly distributed throughout the seasons where summer precipitation is more intense indicating that winter precipitation is more frequent [36]. Long-term average monthly accumulated precipitation records strongly correlate with the average number of rainy days ( Figure S1, Supplementary Materials).

Study Area
The southern Australian agricultural system is dominated by dryland agriculture where cereals (e.g., wheat, barley, and oats), oilseeds (e.g., canola), and legumes (e.g., lupins, chickpeas, field peas, and soybeans) are planted, often in rotation with annual pastures and fallows. Australian wheatbelt ( Figure 1) spans an extremely variable agroecological environment with respect to the climate across the continent, leading to the high spatial variability in Australian grain production. The precipitation varies enormously across the country ( Figure S1, Supplementary Materials); winter (June-August) precipitations are dominant in Western Australia and South Australia, while summer (December-February) precipitations are dominant in Queensland and northern New South Wales. In southern New South Wales and Victoria, total precipitation is more uniformly distributed throughout the seasons where summer precipitation is more intense indicating that winter precipitation is more frequent [36]. Long-term average monthly accumulated precipitation records strongly correlate with the average number of rainy days ( Figure S1, Supplementary Materials).  [13]. The precipitation data are sourced from Jeffrey et al. [37].   (2009)(2010)(2011)(2012)(2013)(2014)(2015). These images are nadir-corrected, adjusted for bidirectional reflectance distribution function, and topographically corrected followed by Li et al.'s methods [38], resulting in a 25-m pixel resolution. This is as, in Geoscience Australia's processing stream, they oversample the Landsat imagery to 25 m to allow easier integration with other remote sensing data sources. Time series of normalized difference vegetation index (NDVI) are then computed using the red and near-infrared bands [39]. The 16-day composite of MODIS NDVI (MOD13Q1 v006) data are used due to the low clouds, low view angle, and the highest NDVI value at 250-m spatial resolution. The MODIS data are sourced from United States Geological Survey for the corresponding periods and area (16 tiles) and are then downscaled to 25 m by 25 m spatial resolution using nearest neighbor interpolation [40] for input to the blending algorithm and for consistency to enable further analysis.
The inter-annual variability of precipitation ( Figure 2b) illustrates that the probability of rain days is higher during June-August (i.e., tillering to flowering phase), which is a critical period for crop yield prediction [41,42]. Figure 2b also shows the value of 0.58 in July at the 75th quartile, which indicates that the probability of cloud-free Landsat images could be lower than 42% in most areas. Here, 75% of field-scale Landsat missing data occur during the growing season (i.e., 113-289 days of the year (DOY)), and only 25% of the Landsat data are complete sequences ( Figure 3).

Satellite Images and L-M Blended Data
Time series of Landsat-5 Thematic Mapper (TM), Landsat-7 Enhanced Thematic Mapper (ETM)+, and Landsat-8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) standardized surface reflectance data are sourced from Geoscience Australia (http://geoscienceaustralia.github.io/digitalearthau/index.html), which contains 239 scenes covering the Australian wheatbelt across seven years (2009)(2010)(2011)(2012)(2013)(2014)(2015). These images are nadir-corrected, adjusted for bidirectional reflectance distribution function, and topographically corrected followed by Li et al.'s methods [38], resulting in a 25-m pixel resolution. This is as, in Geoscience Australia's processing stream, they oversample the Landsat imagery to 25 m to allow easier integration with other remote sensing data sources. Time series of normalized difference vegetation index (NDVI) are then computed using the red and near-infrared bands [39]. The 16-day composite of MODIS NDVI (MOD13Q1 v006) data are used due to the low clouds, low view angle, and the highest NDVI value at 250-m spatial resolution. The MODIS data are sourced from United States Geological Survey for the corresponding periods and area (16 tiles) and are then downscaled to 25 m by 25 m spatial resolution using nearest neighbor interpolation [40] for input to the blending algorithm and for consistency to enable further analysis.
The inter-annual variability of precipitation ( Figure 2b) illustrates that the probability of rain days is higher during June-August (i.e., tillering to flowering phase), which is a critical period for crop yield prediction [41,42]. Figure 2b also shows the value of 0.58 in July at the 75th quartile, which indicates that the probability of cloud-free Landsat images could be lower than 42% in most areas.
Here, 75% of field-scale Landsat missing data occur during the growing season (i.e., 113-289 days of the year (DOY)), and only 25% of the Landsat data are complete sequences ( Figure 3).   Figure 1) and (b) monthly probability of rain days (bottom). For both parts, the horizontal line represents the median of the data, the box spans from 25th to 75th quartiles of the data, and the circles past the end of the whiskers are outliers, while the rain day threshold is 0 mm/day. Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 20 25th to 75th quartiles of the data, and the circles past the end of the whiskers are outliers, while the rain day threshold is 0 mm/day. To blend MODIS and Landsat NDVI data and create the 16-day time series of gap-filled Landsatlike images, we applied ESTARFM [11], which is superior when the spatial variance is dominant [10]. We follow Jarihani et al.'s [22] recommendation to "index then blend", as it yields more accurate results and provides higher computational efficiency than the "blend then index" approach. Figure  S2 (Supplementary Materials) illustrates the blending process and the resulting Landsat-like image. As a result of data blending, each pixel is described by 22 observations across the calendar year. For more details about ESTARFM see [11,43].

Yield Data
Australia's major dryland crops are selected for yield estimates and validation, including 35 fields of canola, 123 fields of wheat, and 52 fields of barley (Table S1, Supplementary Materials). The average area of these fields is ~112 ha with a standard deviation of ~69 ha. Using various yield monitoring systems mounted on grain harvesters operated by farmers, these observed data were collected at these sites across the country from 2009 to 2015 ( Figure 1). The data obtained by these commercially available yield monitors were used to construct a yield data image at 5-m resolution [44], which was upscaled to 25 m and 250 m to match the Landsat and MODIS resolutions, respectively.

Climate Data
Climate data were sourced from Science Information for Land Owners (SILO) which provides nationwide meteorological variables (e.g., maximum and minimum air temperature, and precipitation) at daily temporal frequencies by interpolating observations made by the Australian Bureau of Meteorology onto a 0.05° by 0.05° grid [37].

Yield Prediction
We use a semi-empirical model [C-Crop; 6]) because of its low data requirements for calibration across large areas. C-Crop correlates actual grain weight (t/ha) to end-of-season above-ground plant carbon mass (C), and the estimation of C is based on biophysical principles of plant photosynthesis [6].
is estimated using the carbon mass ( −1 ) from the previous time step and the current period's allocated net assimilation flux (gCm −2 ) minus the dead biomass that enters the litter store, and i (1-22) is the model time step (every 16 days in a calendar year). To blend MODIS and Landsat NDVI data and create the 16-day time series of gap-filled Landsat-like images, we applied ESTARFM [11], which is superior when the spatial variance is dominant [10]. We follow Jarihani et al.'s [22] recommendation to "index then blend", as it yields more accurate results and provides higher computational efficiency than the "blend then index" approach. Figure S2 (Supplementary Materials) illustrates the blending process and the resulting Landsat-like image. As a result of data blending, each pixel is described by 22 observations across the calendar year. For more details about ESTARFM see [11,43].

Yield Data
Australia's major dryland crops are selected for yield estimates and validation, including 35 fields of canola, 123 fields of wheat, and 52 fields of barley (Table S1, Supplementary Materials). The average area of these fields is~112 ha with a standard deviation of~69 ha. Using various yield monitoring systems mounted on grain harvesters operated by farmers, these observed data were collected at these sites across the country from 2009 to 2015 ( Figure 1). The data obtained by these commercially available yield monitors were used to construct a yield data image at 5-m resolution [44], which was upscaled to 25 m and 250 m to match the Landsat and MODIS resolutions, respectively.

Climate Data
Climate data were sourced from Science Information for Land Owners (SILO) which provides nationwide meteorological variables (e.g., maximum and minimum air temperature, and precipitation) at daily temporal frequencies by interpolating observations made by the Australian Bureau of Meteorology onto a 0.05 • by 0.05 • grid [37].

Yield Prediction
We use a semi-empirical model [C-Crop; 6]) because of its low data requirements for calibration across large areas. C-Crop correlates actual grain weight (t/ha) to end-of-season above-ground plant carbon mass (C), and the estimation of C is based on biophysical principles of plant photosynthesis [6]. C i is estimated using the carbon mass (C i−1 ) from the previous time step and the current period's allocated net assimilation flux N i (gCm −2 ) minus the dead biomass that enters the litter store, and i (1-22) is the model time step (every 16 days in a calendar year).
where ρ (periods −1 ) is the reciprocal of carbon longevity (i.e., the turn-over rate of plant live carbon into senesced tissue per time step i). Net primary productivity (NPP) is the rate of carbon assimilation from atmospheric CO 2 to organic material (biomass) for a given area while accounting for the energy loss due to autotrophic respiration [6,45].
where f PAR i 0.95 is the fraction of total assimilation flux G i allocated to the above-ground plant biomass (gCm −2 ) at time step i. The plant maintenance respiration is calculated as a function of leaf nitrogen, air temperature, and previous biomass C i−1 . r 10 is the plant tissue respiration rate at 10 • C; cn stands for plant carbon-to-nitrogen ratio; σ i is a scalar that modifies the respiration rate according to the daily air temperature [46].

Gross Primary Productivity
The total assimilation flux G i , = also known as the gross primary productivity (GPP) (gCm −2 ·day −1 ), can be calculated using the remote sensing-based plant light use efficiency (LUE) approach. The chloroplasts use incoming solar radiation with a spectral range between 400 nm and 700 nm in photosynthesis [47]. [48,49], τ ∂ is atmospheric transmissivity calculated using the Bristow-Campbell relationship [50,51] calibrated for Australia, and ρ sw is the ratio of shortwave irradiance at a sloping surface to that at a horizontal surface [52]. fPAR is the portion of PAR that is absorbed by a photosynthetic organism, and it is estimated using a linear relation between fPAR and rescaled NDVI by thresholds (i.e., local minimum and crop-specific maximum NDVIs) [53,54]. LUE is highly linearly related to a diffuse fraction (f D ) and photosynthetic carbon flux [55].
where A x is the maximum photosynthetic capacity (µmolCm −2 ·s −1 ), which is a crop-specific parameter; is 23, 40, and 45 for barley [56], canola [57], and wheat [57], respectively; f D is the ratio of diffused to total solar irradiance varying from 0.2 (under clear skies) to 1.0 (under overcast skies) [48]. For a full description of C-Crop, see Donohue et al. [6].

Validation
Two sets of data are used for validation and further analysis, depending on the pixel-level completeness of time-series Landsat NDVIs at 25-m pixel resolution across the growing season (i: [8][9][10][11][12][13][14][15][16][17][18][19] between April (DOY 113) and October (DOY 304) ( Figure S3, Supplementary Materials). Firstly, the coordinates of complete time-series Landsat pixels are used to obtain the resampled MODIS, L-M blended, and observed yield data for validation at the 25-m resolution as the first dataset. Secondly, these complete time-series Landsat pixels are upscaled at 250 m pixel size to extract the MODIS and the L-M blended data with the same pixel size for the validation at the moderate resolution. Thirdly, the pixel-level yield predictions are aggregated for each respective field by averaging the yield values of the pixels within, for field-level validation. Fourthly, and finally, the predicted yields are evaluated using the model coefficient of determination (R 2 ), the root-mean-square error (RMSE), and the mean bias error (MBE). The R 2 describes the proportion of the response variable that can be explained by the model. RMSE gives more weight to the largest errors, and the MBE indicates the systematic error of the model to under or overestimate. These procedures are then repeated for the second dataset created according to the coordinates of incomplete time-series Landsat NDVIs at the 25-m pixel resolution. As C-Crop requires a time series of NDVI, the results are assessed without the incomplete time-series pixels of Landsat.

Identification of Threshold for When to Blend
The incomplete time-series pixels of Landsat are gap-filled with the L-M blended pixels (L LM ). The threshold to indicate when blending is beneficial is identified by quantifying the impact of the fraction of missing data on the prediction accuracy at 25 m. Firstly, we group MODIS and L LM time series in eight groups, based on the fraction of Landsat missing data (<10%, 10-20%, 20-30%, 30-40%, 40-50%, 50-60%, 60-70%, and 70-80%). Then, the accuracy of each group is analyzed by calculating the R 2 and the RMSE. The performance of C-Crop using MODIS and L LM is compared by the fraction of missing data in Landsat across time. The threshold can be identified, when the model provides the same R 2 and RMSE using L LM as with MODIS ( Figure S3, Supplementary Materials). This threshold determines when, where, or how much L-M blended data improves crop yield prediction when the fraction of missing data in Landsat is lower than the identified threshold.

Evaluation of the Improvement in Yield Prediction Accuracy
The improvement in prediction accuracy using the identified threshold for multiple spatio-temporal data selection is statistically quantified. More specifically, the threshold is applied to the Landsat observations (2000-2018). Firstly, we compute the temporal probability of optimally using MODIS, Landsat, and L LM images for nationwide crop yield predictions during the past two decades, and then map the results to illustrate the spatial variability of multi-sensor data selection for 25-m pixel-level yield prediction across the wheatbelt. We then evaluate the area percentage of the data sources on a yearly basis and analyze their potential correlation to the annual precipitation (mm/year). Finally, the improvement in the accuracy of predicted yields is evaluated on the field level using MODIS and L LM for Western and eastern Australia, against the reported data [58]. The growing season of 2015 is selected due to the availability of a larger quantity of observed yield data. The incomplete 2015 Landsat series are gap-filled with the blended values corresponding to the threshold value.

Yield Prediction
Exactly 66% of the observed fields have a complete time series of Landsat NDVIs at the pixel level across the growing season. For this dataset, C-Crop performs the best with Landsat images at field (R 2 = 0.68; Figure 4a Figure 6a shows that the 25-m pixel-level yield prediction accuracy fluctuates between 0.35 and 0.75 using MODIS and LLM data in time and space, where time-series Landsat observations are incomplete. The R 2 derived from MODIS is steady (between 0.4 and 0.5) where the fraction of missing Landsat data ranged between 0.05 and 0.42, which is lower than the R 2 (0.62-0.5) derived from LLM data. The RMSE derived from both data sources remains approximately 0.9 t/ha for the same fraction of missing Landsat data (Figure 6b). When more incomplete time-series Landsat data are observed (from 0.42 to 0.75 on both Figure 6a,b), the R 2 based on MODIS increases to around 0.7 and the RMSE decreases to 0.4 t/ha (as MODIS is not as cloud-affected as Landsat due to imagery being acquired on more days), whereas the R 2 derived from the LLM fluctuates between 0.6 and 0.4 and the RMSE changes between 1.3 and 0.8 t/ha ( Figure 6). Given this, up to 42% of missing Landsat data in the growing season is defined as the threshold for when L-M blending is optimal to implement. That is, the 25-m pixel-level yield prediction accuracy can be improved using LLM when the fraction of missing  Figure 6a shows that the 25-m pixel-level yield prediction accuracy fluctuates between 0.35 and 0.75 using MODIS and L LM data in time and space, where time-series Landsat observations are incomplete. The R 2 derived from MODIS is steady (between 0.4 and 0.5) where the fraction of missing Landsat data ranged between 0.05 and 0.42, which is lower than the R 2 (0.62-0.5) derived from L LM data. The RMSE derived from both data sources remains approximately 0.9 t/ha for the same fraction of missing Landsat data (Figure 6b). When more incomplete time-series Landsat data are observed (from 0.42 to 0.75 on both Figure 6a,b), the R 2 based on MODIS increases to around 0.7 and the RMSE decreases to 0.4 t/ha (as MODIS is not as cloud-affected as Landsat due to imagery being acquired on more days), whereas the R 2 derived from the L LM fluctuates between 0.6 and 0.4 and the RMSE changes between 1.3 and 0.8 t/ha ( Figure 6). Given this, up to 42% of missing Landsat data in the growing season is defined as the threshold for when L-M blending is optimal to implement. That is, the 25-m pixel-level yield prediction accuracy can be improved using L LM when the fraction of missing Landsat data at the coordinates is below this threshold. For instance, the L LM -based model increases R 2 by up to 20% when the fraction of the missing Landsat data is below 42% (Figure 6a). Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 20 Landsat data at the coordinates is below this threshold. For instance, the LLM-based model increases R 2 by up to 20% when the fraction of the missing Landsat data is below 42% (Figure 6a).

Evaluation of the Improvement in Yield Prediction Accuracy
Given the availability of cloud-free Landsat and MODIS data across the wheatbelt (2000-2018), in concert with the previously determined threshold ( Figure 6) and finding ( Figure 5), we can demonstrate which imagery (i.e., either MODIS, Landsat, and LLM) is best suitable for 25-m pixel-level crop yield prediction (Figure 7). Figure 7 shows the selection when using multiple satellite products for crop yield estimates across the wheatbelt over the past two decades. For the wheatbelt north of 35° south (S), there is a higher probability of obtaining complete cloud-free Landsat observations over the growing season in the east-west overlapping areas of adjacent Landsat Path-Rows, whilst LLM is optimal elsewhere north of 35° S. South of 35° S, MODIS is optimal, with LLM being optimal in the east-west overlapping areas of adjacent Landsat Path-Rows (in Figure 7).

Evaluation of the Improvement in Yield Prediction Accuracy
Given the availability of cloud-free Landsat and MODIS data across the wheatbelt (2000-2018), in concert with the previously determined threshold ( Figure 6) and finding ( Figure 5), we can demonstrate which imagery (i.e., either MODIS, Landsat, and L LM ) is best suitable for 25-m pixel-level crop yield prediction (Figure 7). Figure 7 shows the selection when using multiple satellite products for crop yield estimates across the wheatbelt over the past two decades. For the wheatbelt north of 35 • south (S), there is a higher probability of obtaining complete cloud-free Landsat observations over the growing season in the east-west overlapping areas of adjacent Landsat Path-Rows, whilst L LM is optimal elsewhere north of 35 • S. South of 35 • S, MODIS is optimal, with L LM being optimal in the east-west overlapping areas of adjacent Landsat Path-Rows (in Figure 7). Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 20 Landsat data at the coordinates is below this threshold. For instance, the LLM-based model increases R 2 by up to 20% when the fraction of the missing Landsat data is below 42% (Figure 6a).

Evaluation of the Improvement in Yield Prediction Accuracy
Given the availability of cloud-free Landsat and MODIS data across the wheatbelt (2000-2018), in concert with the previously determined threshold ( Figure 6) and finding ( Figure 5), we can demonstrate which imagery (i.e., either MODIS, Landsat, and LLM) is best suitable for 25-m pixel-level crop yield prediction (Figure 7). Figure 7 shows the selection when using multiple satellite products for crop yield estimates across the wheatbelt over the past two decades. For the wheatbelt north of 35° south (S), there is a higher probability of obtaining complete cloud-free Landsat observations over the growing season in the east-west overlapping areas of adjacent Landsat Path-Rows, whilst LLM is optimal elsewhere north of 35° S. South of 35° S, MODIS is optimal, with LLM being optimal in the east-west overlapping areas of adjacent Landsat Path-Rows (in Figure 7). For 25-m pixel-level yield prediction, blended data should be preferred on average for 33% of the Australian wheatbelt area (Figure 8a). MODIS and Landsat data remain the preferred data sources in 50% and 17% of cases, respectively. Figure 8b shows that the area percentage is positively correlated with annual precipitation for MODIS and negatively correlated for Landsat and the L-M blended data. MODIS has a larger scatter when annual precipitation is greater than 400 mm/year. crop yield prediction, using the probability of MODIS, Landsat, and LLM images. Blue denotes regions where incomplete Landsat series have a fraction of missing data exceeding 42%, thus indicating where MODIS should be used for yield estimates because it provides more frequent observations than Landsat. Green areas show where adjacent Landsat orbits overlap and, thus, where a complete once every 16-day Landsat series over the whole growing season is available. Areas colored red are those where the fraction of Landsat missing data is below the 0.42 threshold identified previously in Figure  6 when L-M blended data improve the yield prediction accuracy.
For 25-m pixel-level yield prediction, blended data should be preferred on average for 33% of the Australian wheatbelt area (Figure 8a). MODIS and Landsat data remain the preferred data sources in 50% and 17% of cases, respectively. Figure 8b shows that the area percentage is positively correlated with annual precipitation for MODIS and negatively correlated for Landsat and the L-M blended data. MODIS has a larger scatter when annual precipitation is greater than 400 mm/year.  Of the 53-million-ha Australia wheatbelt, complete Landsat time series were suitable for 28.6% of that area, MODIS for 31.5%, and L LM for 39.9% during the 2015 growing season (Figure 8a). In this season, there are 104 fields available for assessing the accuracy improvement in yield prediction accuracy for nationwide crop yield prediction (Table S1, Supplementary Materials). Within these observed yield data, 69 fields are located where L-M blended images can improve the prediction accuracy according to the previously defined (see Figure 6) 42% Landsat missing data threshold. Figure 9 shows the predicted yields against the observed values for 63 fields in Western Australia (i.e., WA), and six fields in eastern Australia (i.e., Queensland (QLD), New south Wales (NSW), Victoria (VIC), and South Australia (SA)). season of 2015 falls below 42%, using LLM reduced the field-level yield prediction errors by 0.15 t/ha for Western Australia and 0.28 t/ha for eastern Australia. More specifically, the bias of 598 kilotons in grain production over the cropping area of 39,882 km 2 in Western Australia and of 1672 kilotons in grain production over the cropping area of 59,716 km 2 in eastern Australia can be decreased when using the threshold determined herein to decide the areas where the blended data can improve the prediction accuracy.

Discussion
Globally, wheat-growing regions are distributed within areas that experience a relatively high frequency of clouds [59]. Persistent cloud cover limits the use of HSR sensors over large spatial extents due to the low frequency of the observations, whereas the other sensors that have HTF are constrained by the LSR. Although blending improves the monitoring of rapidly changing processes such as crop growth, it may introduce unforeseen spatial and temporal variances in the blended images when images are often not observed concurrently [10], and implementing the current suite of algorithms across large areas is currently computationally expensive. Therefore, it is important to systematically quantify and evaluate the utility of blending imagery across time and space before embarking on the development of nationwide blending capabilities.
This study developed a parsimonious approach to identify when and where blending is beneficial for crop yield prediction at the 25-m pixel resolution. When incomplete Landsat series falls below 42% of the possible imagery in the growing season, blending is recommended because it improves the accuracy of the yield estimates. When incomplete Landsat series exceeds 42% of data, the use of HTF LSR data (e.g., MODIS) for the 25-m pixel-level yield prediction is recommended, reducing the computational need for blending. LTF HSR data show benefits, for example, in providing detailed information on plant photosynthesis across space; thus, they should be solely used for crop yield prediction when enough cloud-free images are available throughout the growing season. In these areas (see Figure 9), when the fraction of incomplete Landsat series across the growing season of 2015 falls below 42%, using L LM reduced the field-level yield prediction errors by 0.15 t/ha for Western Australia and 0.28 t/ha for eastern Australia. More specifically, the bias of 598 kilotons in grain production over the cropping area of 39,882 km 2 in Western Australia and of 1672 kilotons in grain production over the cropping area of 59,716 km 2 in eastern Australia can be decreased when using the threshold determined herein to decide the areas where the blended data can improve the prediction accuracy.

Discussion
Globally, wheat-growing regions are distributed within areas that experience a relatively high frequency of clouds [59]. Persistent cloud cover limits the use of HSR sensors over large spatial extents due to the low frequency of the observations, whereas the other sensors that have HTF are constrained by the LSR. Although blending improves the monitoring of rapidly changing processes such as crop growth, it may introduce unforeseen spatial and temporal variances in the blended images when images are often not observed concurrently [10], and implementing the current suite of algorithms across large areas is currently computationally expensive. Therefore, it is important to systematically quantify and evaluate the utility of blending imagery across time and space before embarking on the development of nationwide blending capabilities.
This study developed a parsimonious approach to identify when and where blending is beneficial for crop yield prediction at the 25-m pixel resolution. When incomplete Landsat series falls below 42% of the possible imagery in the growing season, blending is recommended because it improves the accuracy of the yield estimates. When incomplete Landsat series exceeds 42% of data, the use of HTF LSR data (e.g., MODIS) for the 25-m pixel-level yield prediction is recommended, reducing the computational need for blending. LTF HSR data show benefits, for example, in providing detailed information on plant photosynthesis across space; thus, they should be solely used for crop yield prediction when enough cloud-free images are available throughout the growing season.
Annual precipitation, which is strongly associated with cloudiness [60,61] and, thus, partly governs the amount of missing Landsat data, indicates which data sources should be preferably used (Figure 8). For instance, MODIS should be used for the prediction for greater than 50% of the wheatbelt during wet years (e.g., 2003, 2010, 2011, and 2016). During dry years (e.g., 2006 and 2018), the proportion of the wheatbelt where Landsat is suitable for yield prediction increases by >10% compared to the average proportion (17%). For crop yield prediction, blended data can be optimally applied across approximately 33% of the wheatbelt during dry to normal years.
Crop production fluctuates from year to year due to the rapidly changing patterns of precipitation and temperature [42,62,63]. Over the last 50 years, the amount of precipitation received by the Australian wheatbelt declined [64]. Based on model projections, it was also suggested that the changing climatic conditions will affect both the frequency and intensity of the extreme events such as floods, heatwaves, and droughts, which will impact agricultural production [65,66]. Within the context of the present study, an increase in the number of clear sky days can be beneficial for the use of RS in crop yield forecasting. Moreover, a decline in the number of rain days in recent decades within the wheatbelt area can also be noted (see Figure S4, Supplementary Materials). Here, this decline in the number of rain days is assumed to correlate with the increase in the number of days with clear skies. However, a decrease in the number of rain days does not necessarily imply an absence of clouds, as there could be non-precipitating clouds. However, studies such as Norris et al. [67], using observations and Coupled Model Intercomparison Project Phase 5 (CMIP5) model outputs, showed that, on average, there is a decline in cloud amount across the region from 60 • S-60 • north (N). This decline was attributed to the poleward shift of mid-latitude storm tracks, among others.
The potential influence of climate change on grain production suggests a decline in world food supply [68] and an increase in the level of exposure of the global population to the risk of hunger [69]. It is, therefore, important to ensure the accuracy and efficiency of yield prediction at anytime, anywhere, and at any scale. A more precise description of grain weight patterns in time and space provides more accurate information for precision agriculture to improve its production and sustainability. The semi-empirical carbon turn-over model used for crop yield prediction is based on historical yield data, and, in the future, relevant drivers (e.g., precipitation and temperature) and their interactions may change under climate change, thus requiring model re-calibration [4,6,70,71]. Our approach was tested with C-Crop but could be extended to other semi-empirical models [72] and models based on machine learning [73].
Development of a national-scale yield prediction system requires time series of HSR and HTF imagery to describe the crop growth; therefore, next-generation HSR and HTF satellite data like Sentinel-2 [74,75] should be tested for nationwide yield estimates post the growing season of 2015. Harmonized Landsat and Sentinel-2 surface reflectance products should be tested for national-scale crop yield prediction [76] when observed yield data for more recent years are available. However, spatio-temporal fusion of MODIS and Landsat data is still necessary for long-term studies that involve historical satellite images collected before 2015 (back to the year 2000, when MODIS was launched). Future studies should also focus on using active RS technologies, such as radio detecting and ranging (Radar) [77,78] and light detection and ranging (LiDAR) [79], for a national yield estimation because of their ability to penetrate clouds and detect a plant's three-dimensional structural characteristics. Data blending between HSR optical imagery and active RS data (e.g., Radar and LiDAR) [80] warrants further study.

Conclusions
Sparse time series of satellite remote sensing, due to LTF and/or cloud contamination, represent one of the main barriers limiting accurate crop yield estimation at regional to national scales. Blending of HSR but LTF images with LSR but HTF images was proposed to increase the temporal resolution and maintain spatial details. In this study, the benefits of blending were tested for crop yield prediction across the Australian wheatbelt. We found that, when time series are gap-free, yield prediction from Landsat is the most accurate on both field and pixel scales. When Landsat time series contain <42% missing values during the growing season, blending is recommended for nationwide crop yield prediction at the 25-m pixel resolution. When Landsat time series contain ≥42% missing observations, MODIS outperforms blending. Across Australia, these recommendations would improve yield estimates by 0.15 t/ha for Western Australia and 0.28 t/ha for eastern Australia on average. By identifying where and when to blend, this work paves the way to more accurate monitoring of biophysical processes and yields, while keeping computational costs low.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/12/10/1653/s1, Table S1: Summary of observed data (the number of field-years) used herein; Figure S1. Long-term (1901Long-term ( -2018 average monthly accumulated precipitation (mm/month) (left) and the long-term average rain days (right) across the wheatbelt for (a, b) Queensland, (c, d) New South Wales, (e, f) Victoria, (g, h) South Australia, (i, j) Western Australia, and (k, l) Tasmania. The precipitation data are sourced from Jeffrey et al. [37], and the wheat belt mask was resampled from native 2-km 2 horizontal resolution onto the precipitation data grid with 5-km 2 horizontal resolution using nearest neighbor [40]. The vertical black lines (a-f) represent one standard deviation of the monthly wheatbelt specific precipitation for each state. The horizontal line (g-l) represents the median of the data, the box spans from 25th to 75th quartiles of the data, and the rain day threshold is 0 mm/day. It should be noted that Northern Territory is not included in the analysis due to its small cropping area; Figure S2. An example of MODIS and Landsat blended NDVI using ESTARFM. The x-axis and y-axis are time (t, DOY) and spatial resolution (m), respectively. Both MODIS and Landsat have valid observations at t 1 and t 2 . The image at the center of the top row is a valid observation of MODIS at t s , and that at the center of the bottom row is Landsat-like blended data at t s created using the 16-day MODIS composite images on the DOYs 257, 273, and 289, as well as the Landsat images acquired on the DOYs 257 and 289; Figure S3: The implementation flowchart of validation and identification of threshold for when to blend. The shaded areas indicate the original data sources; Figure S4