**1. 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 approach to monitor crop growth uses satellite observations providing repeated synoptic regional assessments [3–6]. High-temporal-frequency (HTF) observations are required to accurately track crop development [7] and predict yield [8], and high-spatial-resolution (HSR) observations are necessary to capture within-field heterogeneity [9].

There is a trade-off between temporal frequency and spatial resolution [10–12] as no single sensor can regularly image vast areas of the Earth used for nation-wide dryland cropping at a high spatial resolution. Commercial options of products combine HTF and HSR images achieved by increasing the number of satellites in orbit, such as PlanetScope; however, such commercial options are not affordable for national-scale monitoring, especially in a country as large as Australia (7.7 million km2), with the southern Australian cereal-based agricultural system notionally covering 530,000 km2 [13]. The Moderate Resolution Imaging Spectroradiometer (MODIS) provides complete coverage of the globe every day at a 250-m spatial resolution from red and near-infrared bands. This resolution constrains the capacity of describing cropping systems, crop growth, and field heterogeneities, especially when fields are small-to-moderate sized and landscapes are fragmented [14,15]. Sensors with higher spatial resolution, such as Sentinel-2 or Landsat, are more suitable for these smaller fields/management units, but their lower temporal frequencies limit their ability to capture rapidly evolving crop processes, especially when factoring in the potential clouds [16]. Lobell et al. [17] used valid Landsat observations (cloud cover <10%) during the growing season to generate optical-based vegetation indices (Vis) and then fitted a multilinear regression between these VIs and a large number of Agricultural Production Systems Simulator (APSIM) simulations for yield prediction. In a country like Australia, most of the arable non-irrigated land grows winter wheat, barley, oats, and canola, and their growth depends on the in-season rainfall; therefore, totally cloud-free time-series observations for their entire growing season are infrequent.

To overcome these challenges, various methods for spatial filtering [18], temporal gap-filling [19–21], and data fusion [22–24] were devised with a varying degree of success. As spatial filtering and temporal gap-filling disregard the spatial and temporal correlation of a pixel, they are highly sensitive to the choice of size/length of the spatial/temporal window [20,25]. Data fusion techniques, on which this article focuses, were shown to improve the temporal resolution of fine-spatial-resolution data by blending observations from sensors with different spatial and temporal characteristics. Prominent examples are the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM; Gao et al. [12]) and the Enhanced STARFM (ESTARFM; Zhu et al. [11]). However, blending is no "silver bullet", as it often introduces unforeseen spatial and temporal variances [10]. Therefore, it is critical to systematically evaluate the benefits of blending and identify where and when blending helps to improve monitoring.

Blending satellite data with complementary frequency and spatial resolution characteristics, such as MODIS (herein considered as an HTF and low-spatial-resolution (LSR) sensor) and Landsat (herein considered as a low-temporal-frequency (LTF) and HSR sensor), provides a solution of synthetic imagery that is both HTF and HSR [26]. Current literature such as Dong et al. [27] found that using Landsat images provides a higher crop yield prediction accuracy for field scales over MODIS images, and a further improvement can be achieved by combining Landsat and MODIS (L–M) blended data with the incomplete Landsat series. To date, the utility of blended output is not yet tested for regional and national crop yield mapping [28–35] (Table 1); this study fills that niche.

This study quantifies and evaluates the benefits of blending satellite data of different temporal frequencies and spatial resolutions for crop yield prediction. The specific objectives of this study are to (i) estimate yields using MODIS, Landsat, and L–M blended data and then compare the yield prediction at both pixel and field scales, (ii) identify the fraction of missing Landsat data during a growing season considering the 16-day acquisition cycle to determine a threshold where the blended data can improve the prediction, and (iii) quantify the improvements in the yield prediction accuracy based on the threshold.



at various scales based on the threshold.

#### **2. Study Area and Data**

#### *2.1. 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).

**Figure 1.** Study sites and the crop growing season (April–October) average accumulated precipitation (mm for the years 2009–2015) across the Australian wheatbelt (~53 million ha) [13]. The precipitation data are sourced from Jeffrey et al. [37].

#### *2.2. Data*

#### 2.2.1. 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–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 2.** Box–whisker plots of 1901 to 2018 averaged (**a**) monthly accumulated precipitation (mm/month) across the Australian wheatbelt (see 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.

**Figure 3.** The quantile of Landsat missing pixels in the fields for which observed yield data are available (2009 to 2015).

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].
