Evaluation of MODIS, Landsat 8 and Sentinel-2 Data for Accurate Crop Yield Predictions: A Case Study Using STARFM NDVI in Bavaria, Germany

Abstract

The increasing availability and variety of global satellite products and the rapid development of new algorithms has provided great potential to generate a new level of data with different spatial, temporal, and spectral resolutions. However, the ability of these synthetic spatiotemporal datasets to accurately map and monitor our planet on a field or regional scale remains underexplored. This study aimed to support future research efforts in estimating crop yields by identifying the optimal spatial (10 m, 30 m, or 250 m) and temporal (8 or 16 days) resolutions on a regional scale. The current study explored and discussed the suitability of four different synthetic (Landsat (L)-MOD13Q1 (30 m, 8 and 16 days) and Sentinel-2 (S)-MOD13Q1 (10 m, 8 and 16 days)) and two real (MOD13Q1 (250 m, 8 and 16 days)) NDVI products combined separately to two widely used crop growth models (CGMs) (World Food Studies (WOFOST), and the semi-empiric Light Use Efficiency approach (LUE)) for winter wheat (WW) and oil seed rape (OSR) yield forecasts in Bavaria (70,550 km²) for the year 2019. For WW and OSR, the synthetic products’ high spatial and temporal resolution resulted in higher yield accuracies using LUE and WOFOST. The observations of high temporal resolution (8-day) products of both S-MOD13Q1 and L-MOD13Q1 played a significant role in accurately measuring the yield of WW and OSR. For example, L- and S-MOD13Q1 resulted in an R² = 0.82 and 0.85, RMSE = 5.46 and 5.01 dt/ha for WW, R² = 0.89 and 0.82, and RMSE = 2.23 and 2.11 dt/ha for OSR using the LUE model, respectively. Similarly, for the 8- and 16-day products, the simple LUE model (R² = 0.77 and relative RMSE (RRMSE) = 8.17%) required fewer input parameters to simulate crop yield and was highly accurate, reliable, and more precise than the complex WOFOST model (R² = 0.66 and RRMSE = 11.35%) with higher input parameters. Conclusively, both S-MOD13Q1 and L-MOD13Q1, in combination with LUE, were more prominent for predicting crop yields on a regional scale than the 16-day products; however, L-MOD13Q1 was advantageous for generating and exploring the long-term yield time series due to the availability of Landsat data since 1982, with a maximum resolution of 30 m. In addition, this study recommended the further use of its findings for implementing and validating the long-term crop yield time series in different regions of the world.

1. Introduction

Crop yields play a significant role in the world’s agricultural development; however, the combined effects of climate change, increase in global population, and degradation of soil and water resources requires main methods that provide a timely and accurate assessment of crop production and contribute towards increasing the sustainability of agricultural food production [1,2,3]. Over the past few years, the growth in publicly available satellite data and the emergence of new technologies has provided the potential to generate and explore a new level of data with different spatial, temporal, and spectral resolutions [4,5,6,7,8]. However, the fundamental requirements of newly generated synthetic data, i.e., their optimal spatial or temporal resolutions in accurately predicting crop yields, still need to be explored [4,9,10].

To ensure the accurate monitoring of crop yields, many studies in the past two decades have started to examine the relationship between plants and their growing environment and proposed crop models to simulate the crop growth status [11,12,13,14,15,16,17,18]. Since then, crop models have advanced in monitoring crop growth from the qualitative to the quantitative level and have been modified from the simulation of the growth process at a plant level to the field or regional level. Over time, crop growth models (CGMs) such as World Food Studies (WOFOST), Agricultural Production Systems Simulator (APSIM), AquaCrop, Cropping Systems Simulation Model (CropSyst), and Light Use Efficiency (LUE) have been refined and updated to simulate better crop growth status and yield [5,11,15,17,18,19,20,21,22]. However, when crop yields are examined at field scales, CGMs need to account for the spatial variation by providing the spatial distribution of climate variables (temperature, precipitation, soil moisture) and biophysical parameters (leaf area index (LAI), biomass, fraction of absorbed photosynthetic active radiation (FPAR)) [23]. The unavailability of spatial information in crop modelling causes uncertainties that affect the whole model’s physiological growth simulation process and leads to more significant errors in crop yield estimation [19,24].

As an alternative, the remote sensing (RS) approach can fill the spatial gap of CGMs by providing timely, ubiquitous, and frequent observations of the land surface at a range of spatial scales [24]. However, having the marked advantages, the RS approach also has important disadvantages. Optical RS data can suffer from significant gaps in the data record due to the cloud and shadow cover that causes uncertainties in the retrieved set of parameters [4,25,26]. Moreover, RS data is limited to only retrieving a few variables of interest, limiting its potential for accurately monitoring agricultural applications [27,28,29]. On the other hand, because of the heterogeneity in agricultural landscapes (various crop and soil types), CGMs have trouble parameterising for large-scale applications. Thus, to accurately monitor crop phenology and improve crop models’ results, combining CGMs and RS data is desirable to obtain the suitability of both realms [5,30]. Many studies have successfully utilised RS-based LAI or FPAR derived from vegetation indices, e.g., the Normalised Difference Vegetation Index (NDVI), in combination with CGMs to estimate crop biomass or yield at different study regions worldwide [5,31,32,33,34,35,36,37,38,39].

For five decades, the availability of RS data has grown historically, globally, and technically in terms of different spatial, temporal, and spectral resolutions, which has created new possibilities for generating accurate datasets for agricultural monitoring [4,5,40]. However, the significant cloud- and shadow-generated gaps in freely available satellite products (such as Landsat and Sentinel-2) hinders RS applications’ accurate and timely-dense monitoring. Therefore, filling the data gaps in the RS input data is more realistic before implementing the synergistic approach (where CGMs are linked with the RS data) for crop monitoring.

To fill the observation gaps in the RS data, spatial-temporal data fusion, where a high spatial resolution product is synchronised with a coarse or low-resolution product, is considered the most effective solution recommended by many studies on detecting vegetation changes [41,42,43,44]. The Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM), which blends the coarse spatial resolution of MODIS and high spatial resolution of Landsat data, was the first initiative in fusion modelling. Since then, many spatiotemporal models have been developed with a successful validation of new synthetic data [6,7,8,45,46,47,48,49]. Moreover, generating new-resolution synthetic products provides geoscience applications with multi-spatial and multi-temporal resolution data. It then outputs different spatial and temporal data of the ground objects [4,5]. However, the potential of newly generated synthetic data obtained from fusion modelling in crop yield predictions using crop modelling still needs to be explored. Inputting RS data with high spatial and temporal resolution could be further used to improve the time series simulation of crop models and increase the models’ simulation accuracy. In addition, the high spatial resolution of RS data could be used to reduce the problem of mixed pixels and then increase the accuracy of different spatial properties at the field scale [19].

In the current study, the STARFM-based synthetic NDVI time series for the application of agriculture is selected from [4,5], where the fusion of MOD13Q1 (20 m, 16 days) is individually achieved with Landsat (30 m, 16 days; L) and Sentinel-2 (10 m, 5–6 days; S). Therefore, intending to investigate the importance of synthetic and real NDVI products with a different spatial and temporal resolutions, this research paper compares other output products which calculate the crop yield of winter wheat (WW) and oil seed rape (OSR) for Bavaria in 2019. The crop yield output of six different RS products (real: MOD13Q1 (250 m, 8 and 16 days); synthetic: L-MOD13Q1 (30 m, 8 and 16 days) and S-MOD13Q1 (10 m, 8 and 16 days)) with two widely used CGMs (WOFOST and LUE), for the respective crops, is tested. Eventually, for accurate crop yield modelling of WW and OSR in 2019, this study answers three research questions:

• What is the optimal spatial resolution (10 m, 30 m, or 250 m)?
• What is the optimal temporal resolution (8 or 16 days)?
• Which is the suitable CGM (LUE or WOFOST)?

Investigating RS products’ optimal spatial and temporal resolutions for accurate crop yield predictions using CGMs requires heavy preprocessing of multiple synthetic and non-synthetic remote sensing datasets. Therefore, knowing the suitable data inputs for crop modelling would save time and computation power for future crop yield prediction and precision farming studies.

2. Materials and Methods

The general workflow of this study is shown in Figure 1. The flow diagram is divided into (1) Data fusion and (2) Crop yield modelling for 2019. The first part was a testing phase that investigated the suitable synthetic NDVI products (which were L-MOD13Q1 and S-MDO13Q1) for the agricultural land cover class of Bavaria for the year 2019 (completed in the preceding work [4]). The “index-then-blend” (IB) technique is used in the previous study to first produce the NDVI from the high pair (Landsat or Sentinel-2) and low pair (MOD13Q1) images before blending them for the data fusion [50]. The IB technique combines only one band, the NDVI. Therefore, it was faster and less expensive to compute. In the second section, the selected output NDVI time series of part 1 (two real: MOD13Q1 (250 m, 8 and 16 days) and four synthetic: L-MOD13Q1 (30 m, 8 and 16 days) and S-MOD13Q1 (10 m, 8 and 16 days)) and the climate elements were used as an input to the LUE and WOFOST models estimating the crop yield of WW and OSR for 2019 in Bavaria. The satellite NDVI and the climate data were selected for the respective start and end of the season for WW and OSR for 2019. Both inputs are masked for WW and OSR using the InVeKos data (source: www.ec.europa.eu/info/index_en, accessed on 21 June 2021).

In the last steps, Bavaria’s obtained crop yield is validated using the Bayerisches Landesamt für Statistik (LfStat) data at the regional level (with a 95% confidence interval). The satellite datasets were downloaded and preprocessed in Google Earth Engine (GEE), and the fusion analysis was performed in R (version 4.0.3) using RStudio.

2.1. Study Area

The federal state of Bavaria is located between 47°N and 50.5°N and between 9°E and 14°E in the southeastern part of Germany (Figure 2). The topography mainly influences the region’s climate, with higher elevations in the south (northern edge of the Alps) and east (Bavarian Forest and Fichtel Mountains). The mean annual temperature ranges from −3.3 to 11 °C, but in most of the territory, the temperature ranges between 8 and 10 °C [4]. The mean annual precipitation sums range from approximately 500 to above 3100 mm, with wetter conditions in the southern part of Bavaria. In 2019, the land cover was highly dominated by forest (36.91%) and agriculture (31.67%) (based on the LC map of Bavaria, 2019). The agricultural areas are mainly found in the northwest and southwest of Bavaria, while forest cover dominates towards the Alps and the east. The other land cover classes include grassland, urban, natural-semi, and water cover of approx. 19.16%, 8.97%, 1.84%, and 1.44%, respectively, for the territory (estimates based on the LC map of Bavaria, 2019) [4]. With an area of approx. 70,500 km², Bavaria covers almost one-fifth of Germany. The federal state is divided into 96 counties with 71 rural districts (so-called “Landkreise”) and 25 city districts (so-called “Kreisfreie Städte”). A brief description of the regions of Bavaria is shown in Figure A1.

2.2. Data

This study investigated relevant satellite data, with different spatial and temporal resolutions used to predict the crop yields of Bavaria on a regional level. Several climate parameters were inputted into the crop models along with the satellite data. Further, the updated InVeKos data of 2019 (https://ec.europa.eu/info/index_en, accessed on 16 February 2023) are used to obtain the reference field information of WW and OSR for every district of Bavaria.

Table 1. A summary of the collected datasets for crop modelling of winter wheat (WW) and oil seed rape (OSR) in 2019. The satellite data used for crop yield modelling are synthetic L-MOD13Q1, S-MOD13Q1, and real Moderate Resolution Imaging Spectroradiometer (MODIS) MOD13Q1. The climate parameters are minimum temperature (°C) (Tmin), maximum temperature (°C) (Tmax), dewpoint temperature (°C) (Tdew), solar radiation (MJ m⁻² day⁻¹) (Rs), sunshine duration (hours) (N), evaporation (mm) (Ep), transpiration (mm) (Tp), runoff (mm) (Roff), and precipitation (mm) (P). InVeKos data provide the fields of WW and OSR for Bavaria in 2019; the Bayerisches Landesamt für Statistik (LfStat) data provide the crop yield information (dt/ha) of WW and OSR at the district level of Bavaria in 2019.

Data	Product Name	Resolution (Spatial-Temporal)	References
Climate data	Tmin, Tmax, Tdew, Rs, N, Ep, Tp, Roff, P	2000 m, 8 and 16 days	https://www.uni-augsburg.de/de/fakultaet/fai/geo/ (accessed on 21 June 2021)
Satellite data	L-MOD13Q1	30 m, 8 and 16 days	[4]
	S-MOD13Q1	10 m, 8 and 16 days	[4]
	MODIS (MOD13Q1)	250 m, 8 and 16 days	https://lpdaac.usgs.gov/ (accessed on 21 June 2021)
Vector data	InVeKos	2019	www.ec.europa.eu/info/index_en (accessed on 21 June 2021)
	LfStat	2019	https://www.statistikdaten.bayern.de/genesis/online/ (accessed on 21 June 2021)

briefly describes the used data and indicates the spatial and temporal resolutions.

2.2.1. Satellite Data

This study employed two freely available spatially high-resolution products obtained from the Sentinel-2 Copernicus program and Landsat 8 Land Surface Reflectance Code (LASRC). The LASRC Tier 1 offers seven spectral bands (coastal/aerosol, blue, green, red, near-infrared (NIR), shortwave infrared (SWIR) 1, SWIR 2) with a spatial resolution of 30 m on a Universal Transverse Mercator (UTM) projection. Using the snow, shadow, and cloud masks, the created C function of the mask (CFMask) method removed snow (Bit 4), clouds (Bit 5), and cloud shadows (Bit 3) using the “pixel_qa” band. After preprocessing, the available snow-free, cloud-free, and shadow-free Landsat images were acquired in 2019 for the state of Bavaria on the following day-of-year (DOY), respectively: 49 (18 February), 81 (22 March), 145 (25 May), and 177 (26 June) (Figure 3).

Moreover, this study also used Sentinel-2 data, which enabled global coverage, five-day return frequency, and multi-spectral imaging with 12 spectral bands at spatial resolutions of 10–20 m. Sentinel-2’s surface reflectance data were processed using the Google Earth Engine after being acquired from the Copernicus Open Access Hub (accessed on 2 August 2021) [51]. The data were computed using sen2cor, which used three quality assessment (QA) bands to create cloud-free images with a QA60 bitmask band containing cloud mask information. After preprocessing, the available Sentinel-2 images were acquired in 2019 for the state of Bavaria at the following DOY, respectively: 49 (18 February), 81 (22 March), 97 (7 April), 113 (23 April), 145 (25 May), and 177 (26 June) (Figure 3).

For data fusion, the coarse resolution MOD13Q1 V6 product was used in this study to generate L-MOD13Q1 and S-MOD13Q1 by fusing it with the preprocessed Landsat and Sentinel-2 data. MOD13Q1 provided an NDVI value per pixel with 250 m spatial and 16-day temporal resolution. In the composed product of MOD13Q1, the NDVI value of a pixel value is assigned with the minor rules and best viewing geometry to the first date of a 16-day time frame. Pixels with constraints (e.g., shadows, clouds) were masked using the quality information (QA) provided along with the NDVI band. Considering the day of acquisition and the QA, linear interpolation of all NDVI values was performed on the product [52] to generate a time series without gaps.

The present study used the synthetic L-MOD13Q1 (30 m, 16 days) and S-MOD13Q1 (10 m, 16 days) NDVI time series generated by [4] as input to the two CGMs obtaining crop yields. Both synthetic products (16 days), L-MOD13Q1 and S-MOD13Q1, were generated using the STARFM and the 8-day products were further developed by applying the linear interpolation approach on 16-day products. The 8- and 16-day time series for RS products were obtained for DOYs from the stem elongation phases until the flowering stages of both WW and OSR. For OSR, the start of the season was 15 February and the end of the season was 20 April 2019 [53]. Moreover, for WW, the start- and end-of-the-season period lay between 15 April and 30 June 2019 [54]. In addition, MOD13Q1 (i.e., just the MODIS NDVI time series without image fusion) was also chosen as an input to the CGMs to allow a comparison between the synthetic and the raw RS time series for crop yield estimation. Based on our previous study, the accuracy assessments of STARFM-generated L-MOD13Q1 and S-MOD13Q1 NDVI products (further used as input for the two CGMs) with the real Landsat and Sentinel-2 NDVI for the agricultural LC class are shown in

Table 2. Statistical analysis between the NDVI values obtained from raw Landsat (L) and Sentinel-2 (S) images with synthetic images (L- and S-MOD13Q1) on different days of the year (DOY). The analysis shows the R² and mean RMSE obtained between the synthetic and reference NDVI in Bavaria for the land cover class of agriculture (which covers 31.67% of Bavaria) and for the overall classes (including agriculture, forest, urban, water, natural, semi-natural, and grassland) in 2019.

NDVI Product	LC Class	DOY
		OSR					WW
		49	81	97	Mean R2	Mean RMSE	113	145	177	Mean R2	Mean RMSE
L-MOD13Q1	Agriculture	0.41	0.49	-	0.45	0.11	-	0.66	0.65	0.65	0.10
	Overall	0.43	0.50	-	0.47	0.11	-	0.61	0.62	0.62	0.11
S-MOD13Q1	Agriculture	0.49	0.74	0.85	0.69	0.10	0.76	0.50	0.60	0.62	0.12
	Overall	0.48	0.67	0.80	0.65	0.13	0.81	0.64	0.65	0.70	0.13

[4]. Our previous study briefly discussed the accuracy assessment of different spatial, temporal products [4]. However, the present study also evaluated synthetic NDVI products’ performance by comparing them with the real NDVI products of Landsat, Sentinel-2, and MOD13Q1. This study compared the mean NDVI values for all RS products used by taking 10,000 random points in Bavaria.

2.2.2. Climate Data

The climate data for 2019 with a daily temporal resolution includes variables such as maximum, minimum, and dewpoint temperature, solar radiation, sunshine hours, precipitation, soil moisture, solar radiation, evaporation, and transpiration. This information was obtained by dynamical downscaling the ERA5 reanalysis dataset, provided by the European Centre for Medium-Range Weather Forecasts [55], to a horizontal grid resolution of 2 km using the hydrologically enhanced Weather Research and Forecasting model [56,57]. A more detailed description of the chosen downscaling approach is provided by [58] and [59]. For this study, the daily climate data are aggregated into 8- and 16-day temporal periods and adapted to the CGMs. Like the satellite data, the present study considers the 8-day climate data for the same SOS and EOS for WW and OSR as described in Section 2.2.1.

2.2.3. InVeKos Data

The field-based InVeKos data are used to identify the fields of WW and OSR in 2019 country-wide. InVeKos data are collected through the Integrated Administration Control System (IACS) (www.ec.europa.eu/info/index_en, accessed on 21 June 2021), which is available for all agricultural plots in European Union (EU) countries by allowing farmers to indicate their agricultural areas graphically. In the IACS, European Union countries are responsible for the administration and control of payments to farmers through a principle called shared management.

2.2.4. LfStat Data

The Bayerisches Landesamt für Statistik (LfStat) is the crop yield validation data providing the crop yield information of 29 crop categories, including WW and OSR in Bavaria on a district level in 2019 (source: www.statistikdaten.bayern.de/genesis/online/, Statistics Code: 41241, accessed on 21 June 2021). The LfStat data were used to validate the LUE-modelled yields of WW and OSR. The validation results were used to check the model’s accuracy, consistency, and stability in generating the yield results in the region. The validation was limited to the rural areas, with the city districts excluded (Figure A1).

2.3. Method

2.3.1. WOFOST

The WOFOST model is a process-based mechanistic model that describes crop parameters such as crop biomass and yield by considering crop genetic properties and climatic and management parameters [60]. The significant processes are CO₂ assimilation, transpiration, respiration, phenological development, and dry matter formation (crop biomass). The model states crop biomass as a function of solar radiation, temperature, and daily crop characteristics. LAI is a vital state variable in the WOFOST model to analyse dynamic growth processes. The model simulates the daily crop growth rate, the gross CO₂ assimilation rate that depends on the LAI, and incoming radiation. In addition, LAI is an essential parameter for the calculation of the potential transpiration in the model. The inputted LAI is calculated as a function of NDVI for WW and OSR (

Table 3. A summary of equations used to calculate LAI from NDVI satellite products for WW and OSR. These equations are derived from the growth stage of a crop until the flowering stage. The LAI product is used as an input to the WOFOST model.

Crop Type	Equation	R2	References
WW	$\mathrm{L}\mathrm{A}\mathrm{I}=-16.606+20.612\ast \mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$	0.70	[61]
OSR	$\mathrm{L}\mathrm{A}\mathrm{I}=-0.393+4.769\ast \mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$	0.81	[62]

).

The daily gross assimilation rate of the crop is calculated by the daily absorbed radiation and the photosynthetic characteristics of each leaf, and it further calculates the total carbohydrates (CH₂O) produced. Some fractions of the CH₂O produced are used to provide energy for respiration (maintenance respiration), with the remaining energy converted into dry matter. The model calculates the growth rate as:

$$∆\mathrm{G}={\mathrm{C}}_{\mathrm{e}}(\mathrm{A}-{\mathrm{R}}_{\mathrm{m}})$$

(1)

where ∆G is the growth rate (kg dry matter ha⁻¹ d⁻¹), A is the gross assimilation rate (kg CH₂O ha⁻¹ d⁻¹), R_m is the maintenance respiration rate (kg CH₂O ha⁻¹ d⁻¹), and C_e is the conversion efficiency (kg dry matter kg⁻¹ CH₂O). Based on Monteith’s Principle of Light Use Efficiency, the calculation of total dry matter (kg dry matter ha⁻¹ yr⁻¹) in the WOFOST model is equivalent to the net primary production (NPP) (kg ha⁻¹ yr⁻¹) [20,63].

This study adopted the complete working methodology of this model from the WOFOST 6.0 documentation prepared by [64], which is also available online (www.wur.nl). The complete calculation, with input parameters to calculate the growth of a WW and a OSR biomass, makes the WOFOST model complex in design. The conceptual diagram of the complete WOFOST model is shown in Figure A2. The WOFOST model is calibrated using the values given in Table 4.

2.3.2. LUE

The LUE model is based on a Light Use Efficiency principle [20,21] and is coupled with the RS data using a similar methodology [5,65]. The model is based on a semi-empirical approach and calculates the FPAR [66] and the daily aboveground biomass as:

$$\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}=\sum \left(\mathrm{P}\mathrm{A}\mathrm{R}\ast \mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}\right)\ast \in$$

(2)

$$\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}=1.222\ast \mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}+0.1914$$

(3)

$$\in =\sum \left({{\mathrm{T}\mathrm{m}\mathrm{i}\mathrm{n}}_{min}}^{\prime }\ast {\mathrm{V}\mathrm{P}\mathrm{D}}^{\prime }\ast \mathrm{K}\mathrm{s}\right)\ast {\in }_o$$

(4)

where PAR is photosynthetically active radiation (MJ m⁻² d⁻¹), FPAR is the fraction of PAR absorbed by the canopy, ∈ is the actual Light Use Efficiency (g C M J⁻¹), ∈_o is the actual Light Use Efficiency (g C M J⁻¹), Tmin_min′ is the minimum of the minimum temperature (°C) index, VPD^’ is the vapour pressure deficit (kPa) index, and Ks is the soil moisture stress index. The temperature and vapour pressure indexes are calculated using the minimum and maximum values for the study region. The total aboveground biomass calculated by the LUE model is equivalent to the net primary productivity (NPP) (kg ha⁻¹ yr⁻¹). A brief explanation of the model with a flow diagram is described in our previous study [5]. The specific model is not only selected for its performance but also its high processing speed and low requirement of input parameters compared to the other CGMs. The linear regression equations used to calculate crop yields of WW and OSR for different satellite biomass products using LUE are shown in

Table A1. Summary of linear regression equations used to calculate crop yield from biomass obtained from different satellite products (MOD13Q1, Landsat (L)-MOD13Q1, and Sentinel-2 (S)-MOD13Q1) for WW and OSR using LUE and WOFOST models. The yield obtained is in dt/ha.

Crop Type	Crop Model	Equation	R2
WW	LUE	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{MOD13Q1}=-56.549+0.2231\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{MOD13Q1}$	0.73
WW	LUE	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{L-MOD13Q1}=22.278+0.0743\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{L-MOD13Q1}$	0.82
WW	LUE	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{S-MOD13Q1}=-14.377+0.147\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{S-MOD13Q1}$	0.85
WW	WOFOST	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{MOD13Q1}=52.533+0.0599\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{MOD13Q1}$	0.69
WW	WOFOST	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{L-MOD13Q1}=58.027+0.0537\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{L-MOD13Q1}$	0.75
WW	WOFOST	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{S-MOD13Q1}=58.670+0.0573\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{S-MOD13Q1}$	0.78
OSR	LUE	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{MOD13Q1}=-25.192+0.1158\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{MOD13Q1}$	0.67
OSR	LUE	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{L-MOD13Q1}=-5.823+0.0807\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{L-MOD13Q1}$	0.80
OSR	LUE	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{S-MOD13Q1}=-6.035+0.0816\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{S-MOD13Q1}$	0.82
OSR	WOFOST	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{MOD13Q1}=-4.4375+0.0721\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{MOD13Q1}$	0.62
OSR	WOFOST	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{L-MOD13Q1}=-16.345+0.089\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{L-MOD13Q1}$	0.63
OSR	WOFOST	${\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}_{S-MOD13Q1}=-8.592+0.815\ast {\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}_{S-MOD13Q1}$	0.64

.

Both models (LUE and WOFOST) were calibrated by using values shown in Table 4. This study used a minimum lethal temperature of −2 °C for WW and OSR [67,68,69]. In the other studies, the optimal minimum values of temperature of WW and OSR at growth stages were 10 °C and 12 °C, respectively [67,68,69]. For the Vapor Pressure Deficit (VPD), the present study followed [70], which analysed the environmental impact on leaf gas exchange of WW with minimum and maximum values of 1.5 and 4.0 kPa, respectively. The value used for optimal Light Use Efficiency was 3 gC/MJ [71].

Table 4. Description of model calibration values taken from the related literature for the WOFOST and LUE models, plus, the climate thresholds used to calculate the climate stress indexes used in the design of a model.

Parameter	Description	Model(s)	Value	Units	Reference
ξ	Scattering coefficient	WOFOST	0.2	-	[60]
kdf	Diffusion coefficient	WOFOST	0.72	-	[72]
Am	Gross assimilation rate	WOFOST	4	g/m2	[73]
Ce	Conversion coefficient	WOFOST	0.0399	-	[74]
∈o	Light use efficiency	WOFOST&LUE	3	gC/MJ	[71]
Tmin min	Minimum of minimum temperature	WOFOST&LUE	−2	°C	[67]
Tmin max	Maximum of minimum temperature	WOFOST&LUE	12	°C	[70]
VPD min	Minimum VPD	LUE	1.3–1.5	k Pa	[75,76]
VPD max	Maximum VPD	LUE	3.6–4	k Pa	[75,76]
Zr	Maximum root depth	WOFOST&LUE	1.5–1.8	m	[77]
P	Average fraction of TAW	WOFOST&LUE	0.55	-	[77]

Table 4. Description of model calibration values taken from the related literature for the WOFOST and LUE models, plus, the climate thresholds used to calculate the climate stress indexes used in the design of a model.

Parameter	Description	Model(s)	Value	Units	Reference
ξ	Scattering coefficient	WOFOST	0.2	-	[60]
kdf	Diffusion coefficient	WOFOST	0.72	-	[72]
Am	Gross assimilation rate	WOFOST	4	g/m2	[73]
Ce	Conversion coefficient	WOFOST	0.0399	-	[74]
∈o	Light use efficiency	WOFOST&LUE	3	gC/MJ	[71]
Tmin min	Minimum of minimum temperature	WOFOST&LUE	−2	°C	[67]
Tmin max	Maximum of minimum temperature	WOFOST&LUE	12	°C	[70]
VPD min	Minimum VPD	LUE	1.3–1.5	k Pa	[75,76]
VPD max	Maximum VPD	LUE	3.6–4	k Pa	[75,76]
Zr	Maximum root depth	WOFOST&LUE	1.5–1.8	m	[77]
P	Average fraction of TAW	WOFOST&LUE	0.55	-	[77]

2.3.3. Sensitivity Analysis

This study performed a sensitivity analysis of the LUE and WOFOST models for both WW and OSR in Bavaria in 2019. The values of climate variables were optimised in the design of every model. During the analysis, the impact of climate stress factors was nullified, and the biomass calculation replaced the actual Light Use Efficiency (ε) values with the optimal (ε_o) values.

2.3.4. Statistical Analysis

Both the referenced and the modelled (LUE and WOFOST) crop yields of WW and OSR were validated using LfStat crop yield (with a 95% confidence interval) for 2019, respectively. The quality (R²) and the precision (root mean square error (RMSE)) of the obtained results were calculated using a linear regression model (LRM), which aimed to establish a linear relationship between the referenced (independent variable) and modelled yields (dependent variable) of WW and OSR at different spatial (10, 30, and 250 m) and temporal (8 and 16 days) scales. The statistical parameters used to validate and compare the accuracies of the LUE- and WOFOST-modelled yields with the referenced yield are R² (Equation (5)), Mean Error (ME) (Equation (6)), RMSE (Equation (7), and relative RMSE (RRMSE) (Equation (8)). To compare the yield outputs of both models, this study considered RRMSE < 15% as good agreement, 15–30% as moderate agreement, and > 30% as poor agreement [78]. The lower the value of ME, RMSE, and RRMSE, the better the model performed.

$${\mathrm{R}}^2=\frac{\left(\right({\sum {\mathrm{P}}_{\mathrm{i}}-{\mathrm{P}}^{\mathrm{\prime }}\left)\right({\mathrm{O}}_{\mathrm{i}}-{\mathrm{O}}^{\mathrm{\prime }}\left)\right)}^2}{{\left(\sum {\mathrm{P}}_{\mathrm{i}}-{\mathrm{P}}^{\mathrm{\prime }})\right)}^2{\left(\sum {\mathrm{O}}_{\mathrm{i}}-{\mathrm{O}}^{\mathrm{\prime }})\right)}^2},$$

(5)

$$\mathrm{M}\mathrm{E}=\frac{1}{\mathrm{n}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{i}}-{\mathrm{P}}_{\mathrm{i}}\right)}^2$$

(6)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\mathrm{M}\mathrm{E}},$$

(7)

$$\mathrm{R}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}(\mathrm{\%})=\frac{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}{\frac{1}{\mathrm{n}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{O}}_{\mathrm{i}}}\ast 100,$$

(8)

where P_i is the predicted value, O_i is the observed value, P’ is the predicted mean, O’ is the observed mean value, n is the total number of observations, referenced yield_y is the LfStat yield of every district in 2019, and modelled yield_y is the LUE-generated yield of every district in 2019. The results’ significance was obtained by observing the probability value (p-value) calculated using the LRM with a H₀ that no correlation exists between the referenced and the modelled or synthetic values and an H₁ that the correlation exists. The test was performed at a significance (or alpha (α)) of 0.05. A p-value lower than 0.05 indicates that the model is significant and rejects the H₀ that there is no correlation.

3. Results

3.1. Evaluation of Real (MOD13Q1, Landsat, and Sentinel-2) and Synthetic (L-MOD13Q1 and S-MOD13Q1) Satellite NDVI Products

The spatial visualisation of the products MOD13Q1, Landsat, L-MOD13Q1, Sentinel-2, and S-MOD13Q1 at DOY 145 is shown in Figure 4, respectively. Both synthetic products, L-MOD13Q1 and S-MOD13Q1, had shown higher dependency on their high-resolution products (Landsat and Sentinel-2) than MOD13Q1. Figure 4f shows the spatial location of 10,000 random points that compares real and synthetic NDVI products with their respective low pair (MOD13Q1) and high pair (Landsat or Sentinel-2) products by considering the mean values at different DOYs (Figure 5). Figure 5a,b show the line and box plot comparison of real and synthetic products and their interquartile comparison of NDVI values.

Both synthetic products underestimated the NDVI values compared to their actual NDVI products between DOYs 81 and 145 (Figure 5a). From DOYs 145 to 177, Landsat, L-MOD13Q1, Sentinel-2, and S-MOD13Q1 achieved a mean NDVI of approx. 0.71. The median NDVI values of L-MOD13Q1 and S-MOD12Q1 lie close to their respective high pair product (Figure 5b).

3.2. Statistical Analysis of Crop Yields Obtained from LUE and WOFOST Models for WW and OSR Using Multisource Data in 2019

Both 8- and 16-day NDVI inputs, such as L-MOD13Q1 and S-MOD13Q1 and MOD13Q1, performed significantly for WW and OSR with LUE and WOFOST models (p-value < 0.05); this rejects the H₀ of the LRM that there is no relationship between the modelled and measured crop yield (Figure 6 and Figure 7). The linear regression equations obtained to calculate crop yields of WW and OSR for different satellite biomass products with both models are shown in Table A1. The R² values obtained from the S-MOD13Q1 NDVI (8 and 16 days) products have a higher accuracy compared to the L-MOD13Q1 (8 and 16 days) and MOD13Q1 (8 and 16 days). Based on the R² of the different spatial resolutions of the NDVI products for WW, the models’ descending order is LUE (S-MOD13Q1, 10 m, 8 days), LUE (S-MOD13Q1, 10 m, 16 days), LUE (L-MOD13Q1, 30 m, 8 days), LUE (L-MOD13Q1, 30 m, 16 days), WOFOST (S-MOD13Q1, 10 m, 8 days), WOFOST (L-MOD13Q1, 10 m, 8 days), LUE (MOD13Q1, 250 m, 8 days), WOFOST (S-MOD13Q1, 10 m, 16 days), WOFOST (MOD13Q1, 250 m, 8 days), WOFOST (MOD13Q1, 250 m, 16 days), WOFOST (L-MOD13Q1, 30 m, 16 days), and LUE (MOD13Q1, 250 m, 16 days), with R² values of 0.85, 0.85, 0.82, 0.78, 0.78, 0.75, 0.73, 0.73, 0.69, 0.65, 0.64, and 0.52, respectively. In general, the predicted values by both models with different satellite inputs follow a similar pattern, and none of the models can claim to outclass the others. However, the ME and RMSE values give a complete picture of the model comparisons (8-and 16-day products) and performances (i.e., their quality and precision) with every satellite input. The ME and RMSE of WW from the WOFOST (MOD13Q1 8-day) is slightly lower than the WOFOST (L-MOD13Q1 16-day and S-MOD13Q1 16-day); moreover, the RMSE of the WOFOST (S-MOD13Q1 and L-MOD13Q1 (8-day)) is lower than the WOFOST (MOD13Q1 16-day). The overall results of LUE inputting L-MOD13Q1, S-MOD13Q1, and MOD13Q1 8-to-16-days NDVIs range from 5.46 to 6.32 dt/ha (RMSE), 5.01 to 5.40 dt/ha, and 6.52 to 9.33 dt/ha.

Like WW, the R² of the different spatial resolutions of NDVI satellite products for OSR in descending order are LUE (S-MOD13Q1, 10 m, 8 days), LUE (L-MOD13Q1, 30 m, 8 days), LUE (S-MOD13Q1, 10 m, 16 days), LUE (L-MOD13Q1, 30 m, 16 days), LUE (MOD13Q1, 250 m, 8 days), WOFOST (S-MOD13Q1, 10 m, 8 days), WOFOST (L-MOD13Q1, 10 m, 8 days), LUE (MOD13Q1, 250 m, 16 days), WOFOST (S-MOD13Q1, 10 m, 16 days), WOFOST (L-MOD13Q1, 30 m, 16 days), WOFOST (MOD13Q1, 250 m, 8 days), and WOFOST (MOD13Q1, 250 m, 16 days), with R² values of 0.82, 0.80, 0.80, 0.78, 0.67, 0.64, 0.63, 0.63, 0.63, 0.62, 0.62, and 0.60, respectively. It showed that the LUE model is more accurate at different spatial scales than the WOFOST model. Moreover, the model resulted in higher accuracy for the 8-day products of S-MOD13Q1 and L-MOD13Q1 compared to their 16-day products. The overall results of LUE combining L-MOD13Q1, S-MOD13Q1, and MOD13Q1 8-to-16-days NDVIs range from 2.23 to 2.36 dt/ha (RMSE), 2.11 to 2.39 dt/ha, and 3.02 to 3.40 dt/ha.

For the LUE model for WW, both 8-day products of S-MOD13Q1 (median yield = 71.68 dt/ha) and L-MOD13Q1 (74.65 dt/ha) obtained their yield close to the referenced yield (72.30 dt/ha) (Figure 8). For the LUE model of OSR, the 8-day and 16-day S-MOD13Q1 and L-MOD13Q1 resulted in similar predictions, respectively, where the 8-days (median yield ~33 dt/ha) show the median yield closer to the referenced yield (33.50 dt/ha) than the 16-days (~34 dt/ha) (Figure 8).

Figure 9a,b displayed that the fused products obtained higher R² and lower RMSE values than the non-fused products for WW and OSR. For example, L- and S-MOD13Q1 resulted in an R² = 0.72 and 0.76 and RMSE = 4.91 and 4.49 dt/ha, respectively, and MOD13Q1 resulted in an R² = 0.63 and RMSE = 5.85 dt/ha. Analysing the different temporal resolutions of 8- and 16-day products with LUE and WOFOST models, the 8-day products (median R² = 0.77, RMSE= 6.14 dt/ha) resulted in higher R² and lower RMSE than the 16-day products (median R² = 0.69, RMSE= 8.0 dt/ha) (Figure 9c,d).

3.3. Spatial Analysis of Crop Yields Obtained from LUE and WOFOST Models for WW and OSR Using Multisource Data in 2019

The spatial comparison of crop yield at the regional level from the referenced and modelled yields with multi-source data was displayed for both WW and OSR (Figure 10, Figure 11, Figure 12 and Figure 13). For WW, the LUE model showed consistency in yield prediction for regions such as Straubing-Bogen, Bad Kissingen, Landsberg am Lech, Dillingen a.d. Donau, Fresing, Würzburg, Neuburg-Schrobenhausen, Fürth, Neustadt a.d. Aisch, Bad Windsheim, Rhön-Grabfeld, Oberallgäu, Regensburg, Aschaffenburg, and Ansbach for all satellite inputs. However, the WOFOST model showed stability for regions such as Freising, Tirschenreuth, Neustadt a.d. Waldnaab, Kitzingen, Fürth, Schweinfurt, Weißenburg-Gunzenhausen, Neustadt a.d. Aisch-Bad Windsheim, and Kulmbach. The S-MOD13Q1 8-day showed higher spatial accuracy than other remote sensing inputs used in both models. The S-MOD13Q1 8-day product with LUE predicted a higher yield of more than 85 dt/ha for regions such as Altötting, Passau, Straubing-Bogen, Deggendorf, Fürstenfeldbruck, Donau-Ries, Ebersberg, and Unterallgäu, like the referenced yield (Figure 11a). However, when inputted into the WOFOST model, the exact product underestimated the yield for all regions (except Fürstenfeldbruck and Unterallgäu) (Figure 11b).

Similarly, for the OSR, both models showed consistency in yield prediction in regions such as Ebersberg, Eichstätt, Lichtenfels, Würzburg, Roth, Schweinfurt, Dingolfing-Landau, Neustadt a.d. Waldnaab, Pfaffenhofen a.d. Ilm, Kelheim, and Mühldorf a.Inn for all satellite inputs (Figure 12). The WOFOST model had overestimated the crop yields with MOD13Q1 (8 and 16 days) for nearly 18 regions (>40 dt/ha) compared to the referenced yield (Figure 13a,b). The L-MOD13Q1 8-day resulted in an overestimation of crop yields compared to the L-MOD13Q1 8-day product with both LUE and WOFOST models.

3.4. Sensitivity Analysis

The sensitivity analysis compared the models’ (LUE and WOFOST) performance by excluding the effect of climate stress factors for both WW and OSR in Bavaria in 2019. The LUE- and WOFOST-modelled yields showed a higher correlation with the referenced yield when the climate stress factors were included and vice versa. Both models showed higher R² and lower RMSE values compared with the yield values obtained during the sensitivity analysis (Figure 14). The overall accuracies obtained during the sensitivity analysis of both LUE and WOFOST were recorded as R² of 0.61 and 0.58 and RMSE of 6.13 dt/ha and 6.32 dt/ha, respectively (Figure 14). Including climate parameters improved both models’ performance, reducing the RMSE by −38% (LUE) and −11% (WOFOST) and increasing the R² from 19% to 12%, respectively.

3.5. Suitable Crop Growth Model

The statistical analysis showed the R², RMSE, RRMSE, and ME values of the model’s (LUE and WOFOST) performance, including climate stress factors’ effect on both WW and OSR in Bavaria in 2019 (Figure 15). The LUE model resulted in a higher R² (>0.78) for the 8- and 16-day products of L-MOD13Q1 and S-MOD13Q1 than the WOFOST model (R² < 0.71). Similarly, the RMSE and ME of these products show more accurate results with the LUE model (RMSE <4.5 dt/ha, ME <3.3 dt/ha) than the WOFOST model (RMSE < 7.0 dt/ha, ME < 6.0 dt/ha). MOD13Q1 8-day (R² < 0.66, RMSE < 5.19 dt/ha, ME < 4.03 dt/ha) achieved higher accuracy than MOD13Q1 16-day (R² < 0.62, RMSE < 7.10 dt/ha, ME < 5.86 dt/ha) with both LUE and WOFOST. The RRMSE for both models show better agreement (<15%) between the observed and modelled yields for all satellite products. However, the LUE model (<11.50%) showed an overall lower RRMSE than the WOFOST model (<13.67%) at different spatial and temporal scales. Irrespective of the crop type and satellite spatial scale, the LUE model obtained higher R² and lower RRMSE (average R² = 0.77, RRMSE = 8.17 %) than the WOFOST model (average R² = 0.66, RRMSE = 11.35%) (Figure 16).

3.6. Visualisation of the Modelled Crop Biomass by the LUE Model in 2019

Figure 17 shows the spatial distribution of simulated crop biomass for WW and OSR by the LUE model at 10 and 30 m spatial resolutions with 8- and 16-day temporal resolution for Sentinel-2 (S-MOD13Q1) and Landsat (L-MOD13Q1), respectively. The field-based biomass OSR and WW biomass values vary between 771.36 and 1112.58 g/m², respectively. These values were obtained considering the climate stress factors, such as temperature, VPD, and soil moisture stress. Every figure shows the difference between the 8-day and 16-day biomass products. The difference in 8- and 16-day WW products varies between −72.57 g/m² and 80.50 g/m², respectively. The results indicate that for WW, S-MOD13Q1 had almost similar results at both temporal resolutions; however, a slight variation in L-MOD13Q1 was seen. For OSR, a little difference in the field-based biomass was observed in both 8- and 16-day products of Sentinel-2 and Landsat. The 8-day products in WW and OSR for L-MOD13Q1 and S-MOD13Q1 showed an overestimation in crop biomass compared to the 16-day products.

4. Discussion

This study finds the RS data’s optimal spatial and temporal resolutions combined with CGMs for accurate crop yield predictions for Bavaria in 2019. The results are obtained using WOFOST (complex model) and LUE (simple model) CGMs by individually inputting them with climate variables plus six different remote sensing products (Real: MOD13Q1 (250 m), Synthetic: L-MOD13Q1 (30 m) and S-MOD13Q1 (10 m)) at 8 or 16 days of temporal resolution. This study investigates the significance of the more straightforward, more accurate, and faster LUE model with less input requirement than the complex WOFOST model with a high demand of climatic input variables. The sensitivity analysis is performed to obtain the influence of climate stress factors on crop yield predictions with different satellite inputs. The following sections provide a brief discussion of the points mentioned above.

4.1. Importance of the Synthetic Data in Crop Yield Modelling

Many studies employing satellite images aimed to compensate for the gaps in the primary data by fusing data with another source for various remote sensing applications [4,5,10,79]. The data fusion is to increase the spatial resolution of the relatively coarse images collected by satellites with high revisit frequencies. The fused results usually inherit the details of the high spatial resolution images and the temporal revisit of the frequencies of their counterparts [79]. In the past two decades, data fusion techniques such as the STARFM and its variants have been applied to satellite images for different applications, including phenology analysis [80,81,82], yield and drought monitoring [83,84,85], forest mapping [46,86], and biophysical parameter estimation [85,87,88,89]. Landsat and MODIS images dominate data fusion; however, other satellite combinations, such as Sentinel-2, Sentinel-1, or Worldview, are being increasingly adopted. However, despite its advantages, the data fusion technique could have challenges. For example, combining different sensors could result in misalignment and inaccuracy. In addition, lower sensor quality in data fusion can affect the results’ accuracy [90]. Therefore, to analyse the quality of data fusion products, this study evaluated the significance of real and synthetic NDVI products by considering the mean NDVI of 10,000 randomly selected points and comparing their mean values at different DOYs. For both L-MOD13Q1 and S-MOD13Q1, the mean values obtained lie close to their respective high-resolution product. Therefore, the accuracy of these products is higher. It leads this study to investigate the potential of synthetic products in crop yield modelling.

Moreover, to further improve the data quality and reduce the computation cost of data fusion, this research employs the “index-then-blend” (IB) technique, which creates the NDVI from both the high pair and low pair images before blending them for data fusion [50]. The results of a preliminary study [23] also indicate that the STARFM could successfully fuse MODIS (MOD13Q1, 250 m, 16 days) with Landsat (output: L-MOD13Q1, 30 m, 16 days) and Sentinel-2 (output: S-MOD13Q1, 10 m, 5–6 days) imagery using the above approach [4,5,91,92]. The low RMSE and high R² obtained for the agricultural class with both L-MOD13Q1 (R² = 0.60, RMSE = 0.11) and S-MOD13Q1 (R² = 0.68, RMSE = 0.13) through the STARFM are comparable to those obtained by other studies [20,39,40,87]. One of our previous studies stated the high potential of data fusion between Landsat and MCD43A4 MODIS products on achieving an R² of 0.61 and RMSE of 0.10 for WW biomass monitoring at Mecklenburg–West Pomerania in Germany [5]. The higher correlations between the observed and predicted NDVI values indicate the suitability of the algorithm for vegetation monitoring. Other studies with spatiotemporal data fusion have used NDVI as their primary input for applications such as crop biomass and yield monitoring [80,81,82,83,84,85,87,88,89]. The present study highlights the importance of the synthetic NDVI time series in crop yield modelling by analysing the accuracy assessment between the raw satellite imagery MOD13Q1 (without fusion) and L and S-MOD13Q1 (with fusion). The crop yield prediction results conclude the need for data fusion (obtaining high-resolution satellite data) for accurate crop yield prediction. Many studies demonstrated the potential of high spatial and temporal remote sensing data to describe the spatiotemporal variability of crop biophysical parameters [93], where the availability of Landsat and Sentinel-2 images offer new perspectives for crop monitoring and modelling.

This study obtains higher crop yield accuracy with the fused products (L- and S-MOD13Q1: R² = 0.72 and 0.76 and RMSE = 4.91 and 4.49 dt/ha) than the non-fused product (MOD13Q1: R² = 0.63 and RMSE = 5.85 dt/ha) for both WW and OSR irrespective of the crop model (LUE or WOFOST) (Figure 9a,b). While comparing the yield prediction accuracies of both fused products, S-MOD13Q1 results are more accurate than the L-MOD13Q1. Due to its higher temporal frequency, Sentinel-2 (5–6 days) had six cloud-free scenes (DOYs: 49, 81, 97, 113, 145, and 177) than the Landsat (16 days), with only four cloud-free scenes (DOYs: 49, 81, 145, and 177) available for the analysis (Figure 2). Due to this lower gap in Sentinel-2 DOYs, the NDVI-fused product (S-MOD13Q1) results in higher accuracy than the Landsat-based product (L-MOD13Q1) [4], which further improves the crop yield prediction accuracy of the former more than the latter. However, the L-MOD13Q1 product is still advantageous for generating and exploring the long-term yield time series due to the availability of Landsat data since 1982 with a maximum resolution of 30 m [10].

Results from previous studies have also shown that the assimilation of RS with high spatial-temporal resolution can significantly improve the accuracy of the output, e.g., with an R² value of 0.86 for the LAI measurements using Sentinel-2 as shown by [94]. Dhillon et al. [5] measured the accuracy of LUE with MODIS and the STARFM; both proved to be more reliable and significant with high R² (>0.64, >0.82) and low RMSE (<650 g/m², <600 g/m²), where MODIS resulted in lower accuracy due its coarser resolution. Further, Huang et al. [95] found that the low spatial resolution of MODIS degrades the output accuracy in crop modelling up to 60%.

The high temporal resolution data help to improve a crop’s accuracy by covering the complete crop stages and measuring climate variables’ impact. The lower the temporal gaps, the higher the attainable accuracies by the crop models [96]. The present study shows that the 8-day products are more accurate for yield prediction than the 16-day products. The 8-day products are more likely to cover fine details of the crop physiology, resulting in higher accuracy. Analysing the different temporal resolutions of 8- and 16-day products with LUE and WOFOST models, the 8-day products (median R² = 0.77, RMSE = 6.14 dt/ha) show a better relationship between the referenced and modelled yields than the 16-day products (median R² = 0.69, RMSE = 8.0 dt/ha) (Figure 9c,d). Therefore, this study concludes that high spatial and temporal remote sensing products are essential for crop growth monitoring influenced by climatic factors [5,9].

Even though the data fusion products obtained in this study resulted in higher accuracy than the non-fused products, many studies have suggested more improvements in the STARFM algorithm [25,37,40]. For example, [97] discussed the inevitable role of different sensors and image-processing algorithms causing inconsistency in the data.

4.2. Importance of Linking Crop Growth Models with RS in Crop Yield Modelling

Crop yield prediction has been considered significant to food security and sustainable agricultural development [29]. This study merged remotely sensed data with process-oriented crop models, which can yield more accurate estimates of model outputs. It gives our approach an advantage over conventional studies that use CGMs [98,99,100]. The current study used the traditional technique of CGMs to monitor WW and OSR yields of Bavaria by integrating STARFM-generated S-MOD13Q1 (10 m, 8 and 16 days) and L-MOD13Q1 (30 m, 8 and 16 days) and raw MOD13Q1 (240 m, 8 and 16 days) products in the two CGMs: WOFOST and LUE.

The performance of both models is compared based on their complexity in design, processing speed, accuracy, and precision. This study found that the WOFOST model, which requires more input parameters, is complex in its design and needs more computing time to generate the output than the LUE model. Compared to the other CGMs, the LUE model is based on the fundamental principles of photosynthesis, considers each crop’s unique properties, and can be calibrated and validated using RS technology [20,21,101]. The model accounts for the crop’s ability to use solar radiation for photosynthesis by correlating its biomass with the amount of solar radiation it receives [102,103]. By using RS-derived NDVI with the amount of solar radiation the crop is absorbing (i.e., APAR), the LUE model can use these variables for its calibration and validation, which makes it more accurate in predicting crop yields [104,105]. The performance of the LUE model in forecasting crop yields also shows consistency with other studies [5,9,106]. Yuan et al. [103] successfully validated the crop yields using the satellite-based LUE model at 36 crop sites. Similar research effectively used the Light Use Efficiency variable for biomass estimation of WW and maize using the Production Efficiency Model [107]. Comparing the results of LUE obtained by [5], the model resulted in an R² of 0.83 and RMSE of 581.82 g/m², which is very close to the results obtained in the present study (R² = 0.81, RMSE = 5.17 dt/ha). Irrespective of the crop type and satellite spatial scale, the results of this study show that the LUE model (average R² = 0.77, RMSE = 4.45 dt/ha) performed more accurately than the WOFOST model (average R² = 0.66, RMSE = 7.75 dt/ha) (Figure 16).

The WOFOST model differs from the LUE model by making the potentially unrealistic assumption that crop growth rates are constant throughout the growing season [108]. For instance, crops may experience periods of stress or damage from pests or diseases, which can affect their growth rate and, ultimately, their yield. It makes the model rely heavily on input data, such as LAI, soil, weather, and management parameters, which may only sometimes be available, and could be the reason for inaccuracies in yield predictions [109]. However, many studies have successfully used the WOFOST model combined with RS-based LAI to predict crop yields accurately and have discussed its potential limitations [110,111]. Similar to this study, a comparison of five different crop growth models was made, where the WOFOST model resulted in an average R² of 0.77 and RMSE of 651 g/m², which matches the results of the present study, where the model for WW resulted in an R² of 0.71 and RMSE of 7.75 dt/ha [5].

Comparing the crop yield at a regional level, the LUE model showed consistency in yield prediction in districts such as Straubing-Bogen, Bad Kissingen, Landsberg am Lech, Dillingen a.d. Donau, Fresing, Würzburg, Neuburg-Schrobenhausen, Fürth, Neustadt a.d. Aisch, Bad Windsheim, Rhön-Grabfeld, Oberallgäu, Regensburg, Aschaffenburg, and Ansbach for all satellite inputs. However, the WOFOST model showed stability for regions such as Freising, Tirschenreuth, Neustadt a.d. Waldnaab, Kitzingen, Fürth, Schweinfurt, Weißenburg-Gunzenhausen, Neustadt a.d. Aisch-Bad Windsheim, and Kulmbach. Both models are uncertain in districts at higher elevations in the south (Bavarian Alps) and east (Bavarian Forest and Fichtel Mountains) of Bavaria for both WW and OSR. The models overestimated the crop yield in regions such as Regen, Freyung-Grafenau, Bad Tölz-Wolfratshausen, and Garmisch-Partenkirchen. The reason could be the complex topography and different climate and management practices of these regions, which have impacted the performance of both models [112,113]. Comparing the two models, the LUE model with S-MOD13Q1 8-day showed higher spatial accuracy than the WOFOST model. Like the referenced yield, the S-MOD13Q1 8-day product with LUE predicted a higher yield of more than 85 dt/ha for regions such as Altötting, Passau, Straubing-Bogen, Deggendorf, Fürstenfeldbruck, Donau-Ries, Ebersberg, and Unterallgäu. However, when inputted into the WOFOST model, the exact product underestimated the yield for all regions (except Fürstenfeldbruck and Unterallgäu). The instability of models at higher elevations could be due to the bad quality of the synthetic NDVI products for specific districts. Like OSR, the WOFOST model overestimated the crop yields with MOD13Q1 (8 and 16 days) for nearly 18 regions by predicting a yield of more than 40dt/ha compared to the referenced yield. S-MOD13Q1 and L-MOD13Q1 8-day performed better when inputted to the LUE model than the WOFOST.

The quality of the synthetic NDVI product might vary for these regions as the districts have no horizontal or vertical overlay of Landsat scenes within the path row, limiting their coverage frequency. Moreover, the continuous cloud cover in some regions of Bavaria could have negatively impacted the yield prediction accuracy of models (Figure 2).

4.3. Sensitivity Analysis

The climate variables have essential contributions impacting the accuracy of crop yield predictions [9,114,115]. This study analyses the impact of climate elements by performing a sensitivity analysis where the LUE and WOFOST models calculate crop yields of WW and OSR without including the climate stress factors in 2019. Having already been influenced by the effect of climate elements, the obtained referenced yield shows poor accuracy with crop model yield results after excluding climate stress factors from both models. This study shows that including climate stress indices improves the performance of both models, reducing the RMSE by −38% (LUE) and −11% (WOFOST) and increasing the R² from 19% to 12%, respectively. Our previous study combined the machine learning approach with crop modelling to identify the impact of every climate element used in crop yield predictions [9]. This study found that solar radiation, soil moisture, and temperature are the most influential variables in increasing the yield accuracy for WW and OSR.

4.4. Outlook

The major outlook is to enhance synthetic NDVI for accurate crop yield predictions of different crop types. Both S-MOD13Q1 and L-MOD13Q1 resulted as reliable input products for the application of crop yield forecasting; therefore, their potential needs to be investigated in different parts of the world. This study validates the crop yield data at a regional level; however, for future studies, validating the CGMs at field level yield data could improve models’ performance and promote sustainable and precision farming. The accurate yield results predicted by this study could be used to investigate the impact of biodiversity or further land use diversity on crop yields at a large scale. As CGMs can only input limited input variables, this study recommends coupling the same methodology with machine or deep learning algorithms to include more climate factors into the analysis for precise results.

5. Conclusions

The present study compares the performance of six different remote sensing products (synthetic: Landsat (L)-MOD13Q1 (30 m, 8 and 16 days) and Sentinel-2 (S)-MOD13Q1 (10 m, 8 and 16 days); real: MOD13Q1 (250 m, 8 and 16 days)) when inputted to crop growth models (CGMs) (WOFOST and LUE) to estimate crop yields of winter wheat (WW) and oil seed rape (OSR) for the entire state of Bavaria in 2019. This study aims to minimise future research efforts by identifying and recommending the most suited synthetic satellite inputs for estimating crop yields by discovering the optimal spatial (10 m, 30 m, or 250 m) and temporal (8- or 16-day) resolutions on a regional scale. Lastly, this study finds the potential of LUE and WOFOST models in generating accurate crop yield results. This research paper concludes the findings as follows:

(i)

To discover the optimal spatial resolution for accurate crop yield predictions, this paper recommends S-MOD13Q1 (10 m) due to its lower uncertainty of mixed pixels information resulting in an increase in the accuracy and precision of the modelled yield. This study obtains higher crop yield accuracy with S-MOD13Q1 (R² = 0.76 and RMSE = 4.49 dt/ha) than L-MOD13Q1 and MOD13Q1 (R² = 0.72 and 0.63 and RMSE = 4.91 and 5.85 dt/ha) for both WW and OSR, respectively. However, the L-MOD13Q1 product is more advantageous for generating and exploring the long-term yield time series due to the availability of Landsat data since 1982, with a maximum resolution of 30 m.

(ii)

To investigate the optimal temporal resolution in yield forecasting, this paper recommends S-MOD13Q1 and L-MOD13Q1 (8-day) as they could improve the accuracy of yield prediction with detailed coverage of crop growth stages and briefly analyse the impact of climate variables simultaneously. The 8-day products (median R² = 0.77, RMSE= 6.14 dt/ha) show a better relationship of referenced yield with the modelled yield than the 16-day products (median R² = 0.69, RMSE= 8.0 dt/ha).

(iii)

To find the suitable crop model with the available input variables, this study finds the LUE model simpler, more reliable, and more accurate than the WOFOST model. Moreover, the LUE model inputs fewer variables, which makes the processing faster than the WOFOST model. Comparably, the LUE model results in a higher mean R² = 0.77 and RMSE = 4.45 dt/ha, while the WOFOST model results in a lower R² = 0.66 and RMSE = 7.75 dt/ha for both WW and OSR yield validations in Bavaria in 2019.

The accurate crop yield measures obtained at the field scale before harvest can contribute to crop yield management decision-making, which could play a crucial role in achieving sustainability in agriculture. However, the availability of field-based yield information in future could be more helpful in testing the potential of high spatial resolution RS products at local scales. The ease of using spatiotemporal modelling with crop growth models would be more comprehensive than one geographical region; therefore, the methodology should be applied globally to obtain food security and maintain biodiversity. For even better accuracy, the synergistic approach of linking RS and CGMs could be linked and tested with machine learning algorithms for various crop fields.