Impact of STARFM on Crop Yield Predictions: Fusing MODIS with Landsat 5, 7, and 8 NDVIs in Bavaria Germany

: Rapid and accurate yield estimates at both ﬁeld and regional levels remain the goal of sustainable agriculture and food security. Hereby, the identiﬁcation of consistent and reliable methodologies providing accurate yield predictions is one of the hot topics in agricultural research. This study investigated the relationship of spatiotemporal fusion modelling using STRAFM on crop yield prediction for winter wheat (WW) and oil-seed rape (OSR) using a semi-empirical light use efﬁciency (LUE) model for the Free State of Bavaria (70,550 km 2 ), Germany, from 2001 to 2019. A synthetic normalised difference vegetation index (NDVI) time series was generated and validated by fusing the high spatial resolution (30 m, 16 days) Landsat 5 Thematic Mapper (TM) (2001 to 2012), Landsat 7 Enhanced Thematic Mapper Plus (ETM+) (2012), and Landsat 8 Operational Land Imager (OLI) (2013 to 2019) with the coarse resolution of MOD13Q1 (250 m, 16 days) from 2001 to 2019. Except for some temporal periods (i.e., 2001, 2002, and 2012), the study obtained an R 2 of more than 0.65 and a RMSE of less than 0.11, which proves that the Landsat 8 OLI fused products are of higher accuracy than the Landsat 5 TM products. Moreover, the accuracies of the NDVI fusion data have been found to correlate with the total number of available Landsat scenes every year (N), with a correlation coefﬁcient (R) of +0.83 (between R 2 of yearly synthetic NDVIs and N) and − 0.84 (between RMSEs and N). For crop yield prediction, the synthetic NDVI time series and climate elements (such as minimum temperature, maximum temperature, relative humidity, evaporation, transpiration, and solar radiation) are inputted to the LUE model, resulting in an average R 2 of 0.75 (WW) and 0.73 (OSR), and RMSEs of 4.33 dt/ha and 2.19 dt/ha. The yield prediction results prove the consistency and stability of the LUE model for yield estimation. Using the LUE model, accurate crop yield predictions were obtained for WW (R 2 = 0.88) and OSR (R 2 = 0.74). Lastly, the study observed a high positive correlation of R = 0.81 and R = 0.77 between the yearly R 2 of synthetic accuracy and modelled yield accuracy for WW and OSR, respectively.


Introduction
Accurate predictions of grain yield at both field and regional scales remain a goal for sustainable agriculture and food security [1,2].The delivery of timely crop monitoring and accurate crop yield estimates is of great value for the formulation of food policies, the regulation of food prices, and agricultural management and is urgently needed for the development of sustainable agriculture [3,4].Among different crop types, oil-seed rape (OSR) (Brassica napus) and winter wheat (WW) (Triticum aestivum) are major crops with high economic value for animal feed, biodiesel production, pollination, biodiversity, and human consumption in the European Union [5,6].In Germany, WW (total production in 2016 was 24.6 million tons) and OSR (4.9 million tons) are crops of significant importance, generally cultivated as high input and conventionally managed monocultures [6][7][8][9].The future climatic changes and increasing climatic variability have diverted the increasing grain yield trend of these crops towards maintaining yield stability [7].Therefore, the accurate yield estimates of WW and OSR could contribute positively to agricultural management practises and optimise resource use to stabilise yields in the future.
Remote sensing (RS) technology can be used to determine and monitor the features of the earth's surface by providing synoptic, timely, and cost-effective information about the earth's surface [10,11].Many studies have implemented RS-based methodologies to estimate the crop production of different crop types at different geographical locations [12][13][14][15][16][17].Landsat (L), SPOT, World View, and Sentinal-2 (S) satellite data with a medium spatial resolution of 10-100 m were utilised to assess and estimate agricultural production at regional and local scales [1,10].The availability of historical RS data since 1972 has also increased the potential of science to invest, design, and implement accurate and reliable methodologies by validating the methods with old yield data sets [18][19][20].Until now, various studies have implemented different methodologies (such as interpolation [21,22], extrapolation [23,24], vegetation indices [25,26], linear regression models [27], crop growth models (CGMs) [17,28], machine learning (ML) [29][30][31], and deep learning (DL) [32,33] using the RS data and accurately predicted crop yields in almost every corner of the world.However, to adequately justify their methods' reliability, stability, and preciseness, very few studies have consistently tested their methodologies for yield prediction for more than five years.
CGMs using the RS data as input parameters successfully attempted to estimate crop yields by covering vast spatial scales and updating the information temporally [17,[34][35][36][37].Many CGMs have been used in crop monitoring for different design purposes, regional environments, and crop types [35].Some very famous models driven by various factors such as radiation, water, or soil are named as AquaCrop [38], soil-water-atmosphere-plant (SWAP) [39], agricultural production systems simulator (APSIM) [40], simple and universal crop growth simulator (SUROS) [41], semi-empiric light use efficiency (LUE) model [42], world food study model (WOFOST) [43], Carnegie-Ames-Stanford Approach (CASA) [44], and the simple algorithm for yield estimate (SAFY) model [45].However, most CGMs are complicated and time-consuming and require many input parameters that could be difficult to obtain or substitute through RS data.LUE and AquaCrop are proven to be more precise, accurate, and reliable by the previous literature [17].However, their performance stability is not determined, as no study has analysed their performance for more than two years at the same study site.
Crop yield prediction at regional, national, and global scales has been conducted based on both climate data and RS data [46].Temperature, solar radiation, and precipitation, as well as the normalised difference vegetation index (NDVI) and leaf area index (LAI), are generally considered the primary climatic and satellite-based input variables used in CGMs [47,48].Therefore, the quality of RS input to CGMs might impact the accuracy of the predicted yield.Even though the RS has broadened the spatial and temporal range of CGMs, the cloud and shadow gaps in the optical satellite data can hinder or limit CGMs from producing accurate yield results [49,50].Many studies have successfully used multitemporal data fusion, combining the data obtained from two different sensors with different spatial and temporal scales, to fill the data gaps [17,[51][52][53].Due to its public availability of code and simplicity of usage, the spatial and temporal adaptive reflectance fusion model (STARFM) [54] is widely used to combine L/S with the moderate resolution imaging spectroradiometer (MODIS) for the application of crop monitoring [55][56][57][58].In a previous study, we tested blending different high (L (30 m, 16 days) and S (10 m, 5-6 days)) and coarse (MODIS: MCD43A4, MOD13Q1, MOD09GQ, and MOD09Q1) spatial resolution products for different land use classes using the STARFM.The study found that both L-MOD13Q1 (30 m, 16 days) (R 2 = 0.62 and RMSE = 0.11) and S-MOD13Q1 (10 m, 16 days) (R 2 = 0.68 and RMSE = 0.13) are suitable for the application of agricultural monitoring, with the former having the upper hand due to its fast and easy processing with lesser storage requirements [59].
Thus, the present study uses the L-MOD13Q1 NDVI product and high-resolution climate parameters (2 km, eight days) as inputs to the LUE model (considered the most accurate, precise, and reliable [17]) for predicting crop yields of WW and OSR at a regional scale for Bavaria from 2001 to 2019.This long-term yield prediction of both crop types would investigate the stability and preciseness of the LUE model by validating the modelled yield with district level Bayerisches Landesamt für Statistik (LfStat) data of Bavaria with a 95% confidence interval.The specific research objectives include: (i) finding the potential of STARFM for blending the long-term NDVI time series; (ii) investigating the preciseness and stability of the LUE model by validating the modelled yield at district level in Bavaria from 2001 to 2019; and (iii) exploring the impact of the fused NDVI input time series on the accuracy of the modelled yields.

Materials and Methods
The general workflow of the study is shown in Figure 1.The flow diagram is divided into three parts: (1) data fusion; (2) generation and validation of L-MOD13Q1 NDVI time series from 2001 to 2019; and (3) comparative analysis between fused (L-MOD13Q1) and non-fused (L-MOD13Q1) products in crop yield modelling 2019; and then, modelling crop yields using L-MOD13Q1 NDVI for WW and OSR from 2001 to 2019.The first part was a testing phase that investigated the suitable synthetic NDVI product (which is L-MOD13Q1) for the agricultural class of Bavaria for the year 2019 (completed in the preceding work [59]).The second section is an extension of the first section, and it generates and validates the NDVI time series of L-MOD13Q1 for eighteen more years (i.e., from 2001 to 2018) using the same methodology as the previous section (as used for 2019).In the third section, the output NDVI time series of part 2 and the climate elements are used as inputs to the LUE model, which estimates the crop yields of WW and OSR from 2001 to 2019 in Bavaria.The satellite NDVI and the climate data are selected for the respective starts and ends of the seasons for WW and OSR from 2001 to 2019.Both inputs are masked for WW and OSR using the InVeKos data that was available from 2005 to 2019 (source: www.ec.europa.eu/info/index_en,accessed on 21 June 2021).
As crop field information was unavailable from 2001 to 2004, InVeKos field data from 2005 to 2009 was used to classify the WW and OSR fields in their respective years.Finally, the obtained crop yield is validated using the LfStat data at the regional level in Bavaria (the regional map is shown in Figure 2).Because the validation data is available at a regional scale, the field outputs of every region were converted to a single regional value by totalling the pixel values of every field.The satellite data sets were downloaded and preprocessed in Google Earth Engine (GEE), and the fusion analysis is performed in R (version 4.0.3)using R-Studio.Overview of the study region.The LC map of Bavaria is obtained by combining multiple inputs of landcover maps, such as the Amtliche Topographisch-Kartographische Informations System, Integrated Administration Control System (which provides the crop field information), and the Corine LC, into one map.Agriculture (peach green) dominates mainly in the northwest and southeast of Bavaria, while forest and grassland classes (dark green and yellow, respectively) dominate in the northeast and south.The LC map is overlayed by the district map of Bavaria.The enlargement (displayed with a dark red box on the top right map) shows the urban area of the city of Würzburg, with the oil-seed rape (OSR) fields (dark orange) and the winter wheat (WW) fields (dark green) in 2019.A brief description of the regions of Bavaria is shown in Figure A1.

Study Area
The study area is Bavaria, which is one of the federal states of Germany located between 47 • N and 50.5 • N and between 9 • E and 14 • E (Figure 2).As the largest state of Germany, Bavaria covers an area of approx.70,550 km 2 , covering almost one-fifth of Germany.The diverse topography of the region, with higher elevations in the south (Bavarian Alps) and east (Bavarian Forest and Fichtel Mountains), impacts the climate of the state.The mean annual temperature ranges from −3.3 • C to 11 • C, and the mean annual precipitation sums range from approx.500 to above 3100 mm.In 2019, about 36.91% of the area of the state was covered by forest, and 31.67% by agriculture [59].More than half of the arable land is used to grow cereals, where WW predominates with 37%, followed by winter barley (25%), summer barley (12%), and grain maize (8%) [60].While OSR predominates in the oil-producing crops in the state.The federal state is divided into 71 Landkreise (rural districts) and 26 Kreisfreie Städte (city districts).A brief description of the regions of Bavaria is shown in Figure A1.

Data
The study collected satellite data (with different spatial and temporal resolutions), climate data, and vector data for the period of 2001 to 2019.A brief description of the data used in the present study, with their spatial and temporal resolutions and references, is shown in Table 1.
Table 1.A summary of the collected datasets for fusion modelling and winter wheat's (WW) and oilseed rape's (OSR) crop modeling.The satellite data used for fusion and crop modelling are Landsat 5, 7, and 8 and Moderate Resolution Imaging Spectroradiometer (MODIS) MOD13Q1; the climate parameters are minimum temperature ( • C) (Tmin), maximum temperature ( • C) (Tmax), dewpoint temperature ( • C) (Tdew), relative humidity (%) (RH), evaporation (mm) (Ep), transpiration (mm) (Tp), and solar radiation (MJm −2 day −1 ) (Rs); Shuttle Radar Topography Mission (SRTM) elevation data of Bavaria; InVeKos data provides the fields of WW and OSR for Bavaria from 2005 to 2019; the Bayerisches Landesamt für Statistik (LfStat) data provides the crop yield information (dt/ha) of WW and OSR at district level in Bavaria from 2001 to 2019.synthetic time series using the STARFM algorithm.With the aim of generating a continuous cloud-free and shadow-free time series (that covers the time frame of 2001 to 2019), highpair data sets such as Landsat 5 Thematic Mapper (TM) (1984 (launched)-2013 (ended)), Landsat 7 Enhanced Thematic Mapper Plus (ETM+) (1999-2003 (stripes in the data after this date due to scan line corrector failure)), and Landsat 8 Operational Land Imager (OLI) (2013-present) were used.The Landsat data arrived with different spectral bands, i.e., coastal/aerosol, blue, green, red, near-infrared (NIR), shortwave infrared (SWIR) 1, and SWIR 2. The snow, shadow, and cloud cover were removed from the Landsat data using the "pixel_qa" quality assessment band generated using the C function of the mask (CFMask) algorithm.The number of cloud-free scenes (0% cloud cover based on CFMask) available every year (N) is shown in Table 2. Due to the difference in surface reflectance and atmospheric conditions, there is a considerable variation between the spectral values of Landsat sensors, which may have significant influences depending on the Landsat data application [61].Therefore, the study performed the inter-sensing harmonisation of the NDVI bands (calculated using NIR and red bands) of Landsat sensors, applying the coefficients proposed by [61] and derived using ordinary least squares (OLS) regression.The pre-processing steps were performed using the platform Google Earth Engine.The Landsat products were generated using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS), which applies atmospheric correction, geometric correction, and calibration procedures to the raw data.During the atmospheric correction step, the impact of atmospheric scattering and absorption is removed from the raw data, and a surface reflectance product is generated independent of atmospheric effects.The geometric correction corrected the viewing angles to remove the effects of the satellite's position and attitude at the time of image acquisition.This correction ensures that the pixels are accurately georeferenced and aligned with each other.Lastly, the calibration procedures applied during the LEDAPS processing correct for spectral band configuration, ensuring that the reflectance values across different spectral bands are consistent and accurate.
In addition, for the low pair, the study selected the MOD13Q1 V6 product, which provides an NDVI value per pixel with 250 m spatial and 16-day temporal resolution.Based on the quality information, pixels with noise (NDVI values <−1 and >+1) were masked out.Both the day of acquisition and quality information were considered while generating the NDVI values from the product.For crop modelling, this study input the eight-day satellite datasets from the stem elongation phases till the flowering stages of both WW and OSR.The parameters for the growth season of OSR were taken from a literature review that accurately monitored the growth timing and condition of the crop based on NDVI and the normalised difference yellowness index (NDYI) (calculated using the green and blue bands [62]) using the unmanned aerial vehicles (UAVs) in Germany [63].
The phenological stages for WW were referenced from the literature that detected the phonological development of the crop using the time series of Sentinel-1 and Sentinel-2 in Germany [64].The study compared the phenology results with the BBCH scale (Biologische Bundesanstalt, Bundessortenamt, and CHemische Industrie), which is a system used worldwide by research and administration to standardise phenologically similar growth stages of multiple plant species [64,65].Therefore, the start (the stem elongation phase) and end (the flowering stage) of the seasons of OSR and WW were taken as 15 February to 20 April [63] and 15 April to 30 June from 2001 to 2019 [64], respectively.

Climate Data
For this study, the climate data from 2001 to 2019 with one-day temporal resolution were obtained by dynamically downscaling the ECWMF reanalysis 5th generation (ERA5) dataset to a horizontal grid resolution of 2000 m using the hydrologically enhanced weather research and forecasting model [66][67][68].The ERA5 data were provided by the European Centre for medium-range weather forecasts.A detailed analysis of the downscaling approach is provided by [69,70].The climate data were used as one of the inputs to the LUE model, which requires temperature, solar radiation, evapotranspiration, and relative humidity (Figure 1).Prior to input to the model, all climate elements were synchronised with the LUE model by aggregating them into eight days of temporal periods.Similar to the satellite data, the present study considered the eight-day climate data for the same start and end of the seasons for WW and OSR as described in the Section 2.2.1.

Elevation Data
The study made use of the shuttle radar topography mission (SRTM) digital elevation data for Bavaria [71].The data had a spatial resolution of 30 m.For this study, the SRTM was used to correlate modelled crop yields with the elevation above sea level.The visualisation of the data is shown in Figure A2.

InVeKos Data
The present study made use of the InVeKos data to obtain the field base information of WW and OSR from 2005 to 2019 for Bavaria.The InVeKos data were collected through the integrated administration control system (www.ec.europa.eu/info/index,accessed on 21 June 2021), which was available for all agricultural plots in European Union (EU) countries by allowing farmers to graphically indicate their agricultural area.

LfStat Data
The Bayerisches Landesamt für Statistik (LfStat) provided the crop yield information for 29 crop categories, including WW and OSR, in Bavaria on a district level from 2001 to 2019 (source: www.statistikdaten.bayern.de/genesis/online/,accessed on 21 June 2021, Statistics Code: 41241).The LfStat data were used to validate the modelled yield information of the LUE model.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 regions, and the city districts were excluded (Figure A1).

Method 2.3.1. STARFM
The STARFM method [54] was used to fuse Landsat and MOD13Q1 to generate the synthetic NDVI time series with high spatial and temporal resolution from 2001 to 2019.As this paper is an extension of our previous paper, the detailed methodology of STARFM's generation of L-MOD13Q1 time series was explained in [17,59].

LUE Model
The LUE model was based on the light use efficiency principle [72,73], and it was coupled with the RS data by using a similar methodology as [17,42].The model was based on a semi-empirical approach and calculated the FPAR [74] and aboveground biomass (at an eight-day temporal resolution) as follows: where PAR is photosynthetically active radiation (MJ m −2 d −1 ), FPAR is the fraction of PAR absorbed by the canopy, SOS and EOS are the start and end of seasons of WW and OSR, and є is the actual light-use efficiency (g C M J −1 ).The total aboveground biomass calculated by the LUE model is equivalent to the net primary productivity (NPP) (kg ha −1 yr −1 ).A brief explanation of the model with a flow diagram was described in our previous study [17].The specific model was not only selected for its performance but also for its high processing speed and low requirement of input parameters as compared to the other CGMs.The model was calibrated by using values from the previous literature, as follows: The study used a minimum lethal temperature value of −2 • C for both WW and OSR [75][76][77].In the other studies, the optimal minimum values of temperature for WW and OSR at growth stages were 10 • C and 12 • C, respectively [75][76][77].For the vapour pressure deficit (VPD), the present study followed [78], which had analysed the environmental impact on leaf gas exchange in WW with minimum and maximum values of 1.5 and 4.0 kPa, respectively.The value for optimal light use efficiency was used as 3 gC/MJ [79].

Sensitivity Analysis
The study performed a sensitivity analysis of the LUE model for both WW and OSR in Bavaria from 2001 to 2019.During the analysis, the impact of climate stress factors was nullified, and the biomass was calculated by replacing the actual light use efficiency (ε) values with the optimal (ε o ) values (Equation ( 4)).

Statistical Analysis
Both the STARFM NDVI and the LUE-modelled crop yield of WW and OSR were validated using the observed NDVI and LfStat crop yield (with 95% confidence intervals) from 2001 to 2019, respectively.The quality (R 2 ) and the precision (root mean square error (RMSE)) of the obtained results were calculated using the linear regression model (LRM), which aimed to establish a linear relationship between the referenced NDVI/or measured yield (an independent variable) and the synthetic NDVI/or modelled yield (a dependent variable).The correlation plots between the number of Landsat scenes and the synthetic NDVI accuracy from 2001 to 2019 were generated by calculating the correlation coefficient (R) (Equation ( 5)).R values lie between −1 (strong negative correlation between two variables) and 1 (strong positive correlation between two variables).The statistical parameters used to validate the accuracy of modelled yield and synthetic NDVI are R 2 (Equation ( 6)), mean error (ME) (Equation ( 7)), and RMSE (Equation ( 8)).The Equation ( 9) calculates the yield percent difference (%), which was calculated for every region of Bavaria.The yield percent difference was analysed in six categories: less than −4, −4 to −2, −2 to 0, 0 to 2, 2 to 4, and more than 4.
Yield Percent Difference = Mean referenced yield y − modelled yield y referenced yield y * 100 (9) where P i is the predicted value, O i is the observed value, P' is the predicted mean, O is the observed mean value, and n is the total number of observations, referenced yield y is the LfStat yield of every district from 2001 to 2019, and modelled yield y is the LUE generated yield of every district from 2001 to 2019.The significance of the obtained results was obtained by observing the probability value (p-value), which was calculated using the LRM with a H 0 that there is no correlation between the referenced and the modelled or synthetic values and an H 1 that the correlation exists.The test was performed at a significance level (or alpha (α)) of 0.05.A p-value lower than 0.05 indicated that the model is significant, and it rejected the H 0 that there was no correlation.The correlation was calculated between the accuracy of synthetic NDVI and crop yield on a yearly basis using Equation ( 3).This calculated the relationship of data fusion with crop yield prediction results by the LUE model.

Validation of Synthetic Remote Sensing Time Series from 2001 to 2019
For nineteen years (i.e., from 2001 to 2019), the STARFM performed significantly for yearly synthetic output (having a p-value < 0.05); this rejected the H 0 of the linear regression model that there was no correlation between the synthetic and referenced NDVI (Figure 3a-s).After generating the yearly scatter plots, the synthetic products' R 2 and RMSE values were analysed.Among all years, the highest accuracy and precision were obtained for 2016 and 2018, with an average R 2 of 0.75 and RMSE of 0.09.For 2005For , 2006For , 2007For , 2009For , 2011For , 2013For , 2014For , 2017, and 2019, the R 2 values were higher than 0.60 and the RMSE values were lower than 0.12.In other years, such as 2003,2004,2008, and 2010, the R 2 and RMSE values lied within 0.60 to 0.62 and 0.10 to 0.14, respectively.However, the rest of the temporal period (i.e., 2001, 2002, 2012, and 2015) resulted in lower R 2 (<0.60) and RMSE (>0.13) values.
The results proved that the yearly accuracy assessment of the synthetic products is impacted by the total number of Landsat scenes (N) available every year (Figure 4a,b).A high positive (R = +0.83)and negative (R = −0.84)correlation was seen between the yearly quality (R 2 ) and preciseness (RMSE) of the synthetic NDVI products with N, respectively.For example, 2011, 2016, and 2018 were the most accurate years (R 2 > 0.68 and RMSE = 0.09) with a total N of more than 7. Similarly, 2001 and 2002 had the least R 2 (< 0.50) and highest RMSE (> 0.15) with the fewest available Landsat scenes (N = 2/3) in both years.The overall accuracy of L-MOD13Q1 for nineteen years was R 2 of 0.62 and RMSE of 0.12, with an average of 5 N every year.(p) (q) (r) (s)

Comparative Analysis between Crop Yield Accuracy of MOD13Q1 and L-MOD13Q1 Using the Light Use Efficiency Model in 2019
Figure 6a-c displayed the crop yield accuracies between the modelled and referenced crop yields of WW and OSR obtained with different satellite products using the LUE model in 2019.The figures show that the fused product (L-MOD13Q1) obtained a higher R 2 (0.

Statistical Analysis between Reference and Modelled Crop Yields of WW and OSR from 2001 to 2019 using the Light Use Efficiency Model
For both WW and OSR, the LUE model performed significantly for every year (having a p-value < 0.05); this rejected the H0 of the linear regression model that there was no correlation between the referenced and modelled crop yield from 2001 to 2019 (Figure A3a-s and Figure A4a-s).After generating the scatter plots, all crop yield products' R 2 , RMSE, and ME values were analysed.For both WW and OSR, the years 2007 through 2018 and 2019 were the most accurate years where the estimated crop yield resulted in high R 2 values (>0.79).However, both 2018 and 2019 for WW resulted in higher RMSE (4.74 and 4.98 dt/ha) and ME (3.46 and 3.71 dt/ha) values, respectively (Figure A3).The remaining years for WW showed a similar trend in R 2   A4).A mostly, similar trend in R 2 values was observed in the OSR, with values ranging from 0.63 to 0.80.The overall accuracies of both WW and OSR for 19 years were recorded as R 2 of 0.79 and 0.86 and RMSE of 4.51 dt/ha and 2.47 dt/ha, respectively (Figure 7a,b).Negative correlations were seen between the regional mean elevations and the modelled yields of WW (−0.30) and OSR (−0.38), respectively (Figure 8a,b).

Statistical Analysis between Reference and Modelled Crop Yields of WW and OSR from 2001 to 2019 Using the Light Use Efficiency Model
For both WW and OSR, the LUE model performed significantly for every year (having a p-value < 0.05); this rejected the H 0 of the linear regression model that there was no correlation between the referenced and modelled crop yield from 2001 to 2019 (Figures A3a-s and A4a-s).After generating the scatter plots, all crop yield products' R 2 , RMSE, and ME values were analysed.For both WW and OSR, the years 2007 through 2018 and 2019 were the most accurate years where the estimated crop yield resulted in high R 2 values (>0.79).However, both 2018 and 2019 for WW resulted in higher RMSE (4.74 and 4.98 dt/ha) and ME (3.46 and 3.71 dt/ha) values, respectively (Figure A3).The remaining years for WW showed a similar trend in R 2   A4).A mostly, similar trend in R 2 values was observed in the OSR, with values ranging from 0.63 to 0.80.The overall accuracies of both WW and OSR for 19 years were recorded as R 2 of 0.79 and 0.86 and RMSE of 4.51 dt/ha and 2.47 dt/ha, respectively (Figure 7a,b).Negative correlations were seen between the regional mean elevations and the modelled yields of WW (−0.

Sensitivity Analysis
The sensitivity analysis compared the model's performance by excluding the effect of climate stress factors from 2001 to 2019 for both WW and OSR in Bavaria.The LUEmodelled yield showed a higher correlation with the referenced yield when the climate stress factors were included, and vice versa.The model showed higher R 2 and lower RMSE values when compared with the yield values obtained during the sensitivity analysis (Figure 9a,b).The overall accuracies obtained during the sensitivity analysis of both WW and OSR for 19 years were recorded as R 2 of 0.68 and 0.78 and RMSE of 5.88 dt/ha and 3.41 dt/ha, respectively (Figure 9c,d).

Sensitivity Analysis
The sensitivity analysis compared the model's performance by excluding the effect of climate stress factors from 2001 to 2019 for both WW and OSR in Bavaria.The LUEmodelled yield showed a higher correlation with the referenced yield when the climate stress factors were included, and vice versa.The model showed higher R 2 and lower RMSE values when compared with the yield values obtained during the sensitivity analysis (Figure 9a,b).The overall accuracies obtained during the sensitivity analysis of both WW and OSR for 19 years were recorded as R 2 of 0.68 and 0.78 and RMSE of 5.88 dt/ha and 3.41 dt/ha, respectively (Figure 9c,d).

Sensitivity Analysis
The sensitivity analysis compared the model's performance by excluding the effect of climate stress factors from 2001 to 2019 for both WW and OSR in Bavaria.The LUEmodelled yield showed a higher correlation with the referenced yield when the climate stress factors were included, and vice versa.The model showed higher R 2 and lower RMSE values when compared with the yield values obtained during the sensitivity analysis (Figure 9a,b).The overall accuracies obtained during the sensitivity analysis of both WW and OSR for 19 years were recorded as R 2 of 0.68 and 0.78 and RMSE of 5.88 dt/ha and 3.41 dt/ha, respectively (Figure 9c,d).

Statistical Analysis between Reference and Modelled Crop Yields of WW and OSR from 2001 to 2019 using the Light Use Efficiency Model at Regional Level
On comparing the long-term crop yield at the regional level, the yearly spatial change from the mean referenced and modelled yield was displayed for both WW and OSR (Figures 10 and 11).For WW, most of the regional yield lied between 65 and 80 dt/ha (Figure 10).Districts such as, Regen, Freyung-Grafenau, Bad Tölz-Wolfratshausen, and Garmisch-Partenkirchen, the average percent difference was calculated as −25.10% (LUE: ~75 dt/ha), −18.68% (~60 dt/ha), −8.08% (~62 dt/ha), and −5.58% (~65 dt/ha), which showed that the model highly overestimated the crop yield values as compared to the referenced yield (Figures 12a and A5a).The positive yield percent difference (where the model underestimated the crop yield) between 0 and 4% had an accuracy greater than 0.80 as compared to the negative yield percent difference between −4 and 0% with an accuracy less than 0.70 (Figure 13).Similarly, the model underestimated the crop yield of Oberallgäu, Miltenberg, Deggendorf, and Dachau with 4.65% (~78 dt/ha), 3.91% (~68

Statistical Analysis between Reference and Modelled Crop Yields of WW and OSR from 2001 to 2019 Using the Light Use Efficiency Model at Regional Level
On comparing the long-term crop yield at the regional level, the yearly spatial change from the mean referenced and modelled yield was displayed for both WW and OSR (Figures 10 and 11).For WW, most of the regional yield lied between 65 and 80 dt/ha (Figure 10).Districts such as, Regen, Freyung-Grafenau, Bad Tölz-Wolfratshausen, and Garmisch-Partenkirchen, the average percent difference was calculated as −25.10% (LUE: ~75 dt/ha), −18.68% (~60 dt/ha), −8.08% (~62 dt/ha), and −5.58% (~65 dt/ha), which showed that the model highly overestimated the crop yield values as compared to the referenced yield (Figures 12a and A5a).The positive yield percent difference (where the model underestimated the crop yield) between 0 and 4% had an accuracy greater than 0.80 as compared to the negative yield percent difference between −4 and 0% with an accuracy less than 0.70 (Figure 13).Similarly, the model underestimated the crop yield of Oberallgäu, Miltenberg, Deggendorf, and Dachau with 4.65% (~78 dt/ha), 3.91%

Correlation Analysis between the Accuracy Assessments of the Input Synthetic Products and the Crop Yield Modelling
The bar and scatter plots compared and linked the yearly accuracies of the input synthetic time series with the crop yield modelling for WW and OSR from 2001 to 2019, respectively (Figures 14 and 15).For WW, the correlation coefficient showed a higher positive correlation of 0.81 between the R 2 of synthetic accuracy and the modelled yield accuracy (Figure 15a).Except 2015 (yield R 2 : 0.77, synthetic R 2 : 0.53) and 2013 (yield R 2 : 0.71, synthetic R 2 : 0.65), where the fusion accuracies were negatively correlated with crop yield accuracy (Figure 14a).Similarly, for OSR, the correlation coefficient was found to be 0.77 (Figure 15b).For 2001 and 2002, the fusion accuracy was lower (R 2 < 0.50); however, the crop yield accuracy for the same years resulted in an R 2 of more than 0.65 (Figure 14b).

Visualisation of the Modelled Crop Biomass and the NDVI of Different Years at a Field Level
The side-by-side spatial visualisation of the input synthetic NDVI product (DOY 169: 18 June) and the WW-modelled biomass for selected years (2005, 2013, and 2019) is shown in Figure 16, respectively.For every year, the spatial trend of crop biomass and NDVI in every field was seen differently.Likewise, NDVI values were rising from 2005 to 2019 from 0.4 to 0.8, and crop biomass had been observed rising from less than 550 g/m 2 in 2005 to more than 850 g/m 2 in 2019.In most of the fields, the crop biomass was dependent on the higher NDVI values.The NDVI values higher than 0.8 impacted higher crop biomass of more than 850 g/m 2 in almost every year.In 2005, the average field crop biomass resulted in less than 650 g/m 2 ; however, in 2019, the crop biomass resulted in more than 650 g/m 2 .

Correlation Analysis between the Accuracy Assessments of the Input Synthetic Products and the Crop Yield Modelling
The bar and scatter plots compared and linked the yearly accuracies of the input synthetic time series with the crop yield modelling for WW and OSR from 2001 to 2019, respectively (Figures 14 and 15).For WW, the correlation coefficient showed a higher positive correlation of 0.81 between the R 2 of synthetic accuracy and the modelled yield  15a). Except 2015 (yield R 2 : 0.77, synthetic R 2 : 0.53) and 2013 (yield R 2 : 0.71, synthetic R 2 : 0.65), where the fusion accuracies were negatively correlated with crop yield accuracy (Figure 14a).Similarly, for OSR, the correlation coefficient was found to be 0.77 (Figure 15b).For 2001 and 2002, the fusion accuracy was lower (R 2 < 0.50); however, the crop yield accuracy for the same years resulted in an R 2 of more than 0.65 (Figure 14b).5)).

Visualisation of the Modelled Crop Biomass and the NDVI of Different Years at a Field Level
The side-by-side spatial visualisation of the input synthetic NDVI product (DOY 169: 18 June) and the WW-modelled biomass for selected years (2005, 2013, and 2019) is shown in Figure 16, respectively.For every year, the spatial trend of crop biomass and NDVI in every field was seen differently.Likewise, NDVI values were rising from 2005 to 2019 from 0.4 to 0.8, and crop biomass had been observed rising from less than 550 g/m 2 in 2005 to  15a).Except 2015 (yield R 2 : 0.77, synthetic R 2 : 0.53) and 2013 (yield R 2 : 0.71, synthetic R 2 : 0.65), where the fusion accuracies were negatively correlated with crop yield accuracy (Figure 14a).Similarly, for OSR, the correlation coefficient was found to be 0.77 (Figure 15b).For 2001 and 2002, the fusion accuracy was lower (R 2 < 0.50); however, the crop yield accuracy for the same years resulted in an R 2 of more than 0.65 (Figure 14b).

Quality Assessment of Synthetic Remote Sensing Time Series from 2001 to 2019
The present study investigates the potential of the STARFM over the Bavarian state of Germany to generate the synthetic NDVI time series from 2001 to 2019 by selecting the best-performing high (Landsat) and low (MODIS) pair obtained for the agricultural class from the previous literature.Many studies prefer using ESTARFM (Enhanced STARFM) for better fusion accuracy [80,81]; however, some studies found STARFM performing significantly better than ESTARFM [82,83].Simple in its design, faster to implement, and capable of fusing the entire state of Bavaria (which covers almost one-fifth of the area of Germany) for two decades, the study finds STARFM to be more advantageous over ESTARFM.ESTARM was complex, time-consuming, and computationally  The present study investigates the potential of the STARFM over the Bavarian state of Germany to generate the synthetic NDVI time series from 2001 to 2019 by selecting the best-performing high (Landsat) and low (MODIS) pair obtained for the agricultural class from the previous literature.Many studies prefer using ESTARFM (Enhanced STARFM) for better fusion accuracy [80,81]; however, some studies found STARFM performing significantly better than ESTARFM [82,83].Simple in its design, faster to implement, and capable of fusing the entire state of Bavaria (which covers almost one-fifth of the area of Germany) for two decades, the study finds STARFM to be more advantageous over ESTARFM.ESTARM was complex, time-consuming, and computationally expensive for covering extensive data for extended periods [84,85].One of the strengths of ESTARFM is that it incorporates additional information, such as a land-cover map, to improve the accuracy of the fusion [80].The study incorporates Bavaria's accurate and updated land cover map into the STARFM to balance its input requirements with the ESTARFM.It provided homogeneity to the STARFM and increased its fusion accuracy (as discussed briefly in our previous study [59]).
In our previous study, we found that L-MOD13Q1 (30 m, 16 days) (R 2 = 0.62 and RMSE = 0.11) was suitable for the application of agricultural monitoring due to its fast and easy processing with lesser storage requirements [59].Moreover, the present study focuses on two decades (2001 to 2019); therefore, the paper generates and validates a Landsat-based synthetic NDVI time series (L-MOD13Q1) due to its continuous availability since 1982 with a maximum resolution of 30 m.As NDVI is among the most effective and widely used vegetation indices, many spatiotemporal fusion-based studies have used it as their primary input [5,30,36,37].However, many spatiotemporal fusion algorithms are based on reflectance fusion, which requires more processing time and storage than NDVI (or one-band blending) fusion [86,87].Having high computation power with fewer storage problems for the long-term time series of 2001 to 2019 for complete Bavaria (70,550 km 2 ), the research uses the strategy "index-then-blend" (IB), which generates the NDVI from Landsat and MOD13Q1 before they are blended for fusion [88].The IB strategy is used in multiple works of the literature with highly accurate and precise fusion outputs [17,59,88].
The analysis found that the accuracies of the fusion products are dependent on the available number of Landsat scenes per year (N) [59], such that the higher N, the higher the fusion accuracy of the synthetic NDVI product in a respective year.For instance, the positive R (+0.75) shows the positive correlation between R 2 of yearly synthetic NDVIs and N (representing the higher quality of the fused product), and the negative R (−0.73) shows the negative correlation between RMSEs and N (representing the higher precision).However, as the research made use of Landsat 8 Operational Land Imager (OLI) (from 2013 to 2019) and Landsat 5 Thematic Mapper (TM) (from 2001 to 2013), it was found that Landsat OLI-based fusion with MOD13Q1 resulted in higher accuracy as compared to Landsat TM [89].For example, the years 2001,2002,2004,2005, and 2012 (using Landsat 5 and 7) have a lower R 2 (<0.60) and a higher RMSE (>0.12) than the remaining years (using Landsat 8).The reason could be that Landsat 8 has improved upon the quality of Landsat 5 and 7, offering improved data accuracy.Moreover, the accuracy of the year 2012 is affected due to the gaps generated by the scan line corrector (SLC) failure in Landsat 7.
On comparing the fusing results on a DOY basis, the study finds that the few cloudfree DOYs could create large gaps between the available Landsat scenes that might affect the accuracy of the fusion product [17,59].For example, the DOYs 33 to 97 (N = ~6) result in a low R 2 (0.54) and a high RMSE (0.16) as compared to the DOYs 113 to 193 (N = ~8), which have a high R 2 (0.64) and a low RMSE (0.10).

Impact of Synthetic Time Series on Crop Yield Modelling
The objective of the present study is to generate and validate the long-term crop yield time series using the semi-empiric LUE model, which has proven to be more reliable, precise, and simple in the previous literature [17,42].The present study validates the crop yield results of WW and OSR obtained by inputting the synthetic NDVI and climate elements to the LUE model at a regional scale in Bavaria from 2001 to 2019.However, before generating the long-term time series using the synthetic NDVI product, the study finds the potential of fused (L-MOD13Q1) in crop yield prediction by comparing its accuracies with the non-fused (MOD13Q1) product in 2019.The study obtains higher crop yield accuracy with the L-MOD13Q1 (R 2 = 0.81 and RMSE = 3.91 dt/ha) than the MOD13Q1 (R 2 = 0.70 and RMSE = 4.77 dt/ha) irrespective of the crop type (Figure 6a,b).It proves the importance of high-resolution synthetic data for accurate modelling of crop yields.
After generating the long-term crop yield time series, the research finds the significant yearly performance of the model for both WW and OSR; however, some years obtained higher accuracy than the others.For example, 2007, 2018, and 2019 are the most accurate years, with R 2 values of more than 0.79 for both crop types.However, the RMSEs of both 2018 and 2019 are relatively higher (>4.74 dt/ha) than in the other years.Similarly, 2011 and 2016, with a higher number of N (~6), result in lower crop yield accuracy than 2007, 2008, and 2011 (N = ~8).This might be due to the impact of climate variables inputted into the LUE model (discussed briefly in Section 4.3).
The study discusses briefly how the quality of the input data fusion product impacts the accuracy of the CGM.For example, due to the low quality of synthetic NDVI products in 2001, 2002, and 2012, which might impact the accuracy of the input FPAR products generated, the yield prediction accuracy of both WW and OSR is low.The analysis tries to prove that even though synthetic time series would be the preferable solution to input a CGM for yield prediction when the quality of the combined fusion product is low, it could negatively affect the crop yield estimation.In relevance to the above point, high positive correlations have been seen when the accuracies of the synthetic NDVI time series are plotted with the accuracies of modelled crop yield from 2001 to 2019 for WW (R = 0.81) and OSR (0.77).For example, the quality of the NDVI time series for the years 2016 and 2018 is higher with R 2 (>0.73), and the crop yield accuracies are also higher with R 2 of 0.83 (WW)/0.81(OSR), 0.85/0.83,respectively (Figure 3p,r).Similarly, the striped data collected from Landsat 7 in 2012 has deteriorated the quality of the synthetic NDVI product (R 2 = 0.51; RMSE = 0.13), which further negatively affected the crop yield estimations for WW (R 2 = 0.62; RMSE = 5.40 dt/ha) and OSR (R 2 = 0.49; RMSE = 4.13 dt/ha) (Figure 3l).Moreover, the Landsat images were available at different times of the year.This has an impact on the prediction accuracy of both crops.For example, the WW yield results are more accurate than the OSR because the synthetic data in late spring and early summer (DOYs 113 to 193) is usually more precise.
The study compares the long-term crop yield time series by calculating the average percent change from the referenced and modelled yields for both crop types.Previous studies found that the elevation plays a significant role in impacting the regional crop yield [90,91].Most of the studies found lower crop productivity at higher elevations due to complex topography, different climates, and management practices [92,93].Moreover, the cropping intensity at lower elevations is higher as compared to the higher elevations.The survey finds negative correlations between the mean regional elevations and the crop yields of WW (−0.30) and OSR (−0.38).The model is precarious in specific regions, especially the districts at higher elevations in the south (Bavarian Alps) and east (Bavarian Forest and Fichtel Mountains) of Bavaria for both WW and OSR.In regions such as Regen, Freyung-Grafenau, Bad Tölz-Wolfratshausen, and Garmisch-Partenkirchen, the model highly overestimates the crop yield, and for regions such as Oberallgäu, Miltenberg, Deggendorf, and Dachau, it underestimates the yield as compared to the referenced yield for WW.This overestimation of WW yield values has resulted in a decrease in accuracy.The model shows yearly stability in predicting crop yields of WW between 65 and 80 dt/ha for most of the regions.The positive yield percent change (where the model underestimated the crop yield) between 0 and 4% had higher accuracy (R 2 > 0.80) as compared to the percent change between −4 and 0% (R 2 < 0.70).For 48 of the 71 total districts, the model performs relatively well, with a percent change between -2% and +2%.However, unlike WW, both the over-and underestimation-yield values have resulted in a similar increase and decrease in the accuracy of OSR.The positive and negative yield percent change (where the model underand over-estimates the crop yield, respectively) between 0 and +/−4% had an accuracy of more than 0.80.For OSR, the model overestimates the yield for Aichach-Friedberg, Deggendorf, Dingolfing-Lindau, Traunstein, Unterallgäu, Dachau, Rottal-Inn, Miltenberg, and Günzberg and underestimates the yield for Roth, Regen, Kronach, Kitzingen, and Bad Tölz-Wolfratshausen.However, for the 27 districts with OSR, the model performs steadily.Interestingly, the regions where the model's performance went unstable were primarily located in the southern alps, except for Regen, Freyung-Grafenau, Kitzingen, Roth, and Miltenberg.The reason could be the instability of the model at higher elevations or the bad quality of the synthetic NDVI products in specific districts.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.

Sensitivity Analysis
Besides the impact of data fusion, climate variables play an essential role in affecting the accuracy of crop yield predictions [31,94,95].To analyse the impact of climate elements, the study performs sensitivity analysis, where the LUE model calculates the crop yields of WW and OSR without including the climate stress factors from 2001 to 2019.As the referenced yield is already influenced by the climate, the results of the study show that the accuracy of crop yield predictions worsens with the exclusion of climate variables, with a lower R 2 (0.68 (WW)/0.74(OSR)) and a higher RMSE (5.88/3.41dt/ha).However, an increase in R 2 (0.79/0.86) and decrease in RMSE (4.51/2.57dt/ha) have been seen when the climate effect is included in the model.As the relationship between climate and crop yield undergoes significant shifts, it might be the reason that some years (2011 and 2016) with higher N (8) obtained lower crop yield accuracy than years (2007, 2018 and 2019) with comparably lower N (6).Moreover, our previous study, which made use of the machine learning approach with the LUE model, identified the impact of every individual climate element used in crop yield predictions [31].
Furthermore, many studies stated that the availability of coarse climate data negatively affected yield prediction accuracy.In a previous study, the coarse spatial resolution of climate data (ECMWF: ~80 km) used to estimate the biomass resulted in low R 2 and high RMSE using CGMs by inputting coarse synthetic NDVI products [17,42].However, while inputting high spatial resolution NDVI products, the low impact of the high spatial resolution of climate elements is observed.The present study inputs high spatial resolution climate data time series (2 km, daily) to the LUE model, resulting in stable yearly accuracies from 2001 to 2019.Notably, selecting climate thresholds according to the geographical location and crop types is essential in achieving high crop yield accuracy [96][97][98].Different climate thresholds are used for WW and OSR, resulting in accurate and stable yield predictions in Bavaria during the study period.

Validation at the District Level
The crop yield validation for the more extended time series of 2001 to 2019 is performed using the LfStat crop yield data (used for validation at a 95% confidence interval) for WW and OSR provided by the Bavarian State Office of Statistics.As the validation data set is provided at a regional scale, pixel-based yield information is converted for both crop types to the regional level.However, transferring the field-based information to the district level could result in some uncertainties in the validation process.For example, in some regions of southern Bavaria (Bad Tölz-Wolfratshausen, Garmisch-Partenkirchen, Traunstein, Unterallgäu, and Oberallgäu), where the model's performance is volatile, this might be due to the uncertainty occurring while transferring the pixel-level information to the district level.The availability of fewer fields of WW and OSR in those regions might be the reason for the model's instability, as the validation data recorded high yield values for the same districts.Therefore, future work should aim to validate crop yield results at the field level, which could help achieve more precise results.Additionally, the availability of field data for FPAR, an important input to the LUE model, would help to validate the FPAR product generated using the NDVI time series.

Conclusions
The present study investigates the relationship of spatiotemporal fusion modelling using STRAFM on crop yield prediction for winter wheat (WW) and oil-seed rape (OSR) using a semi-empirical light use efficiency (LUE) model for Bavaria, Germany, from 2001 to 2019.The research paper concludes the findings as follows: The present study recommends validating crop yields at the field scale, as transferring the pixel-based information to the district level could cause uncertainties in the validation process.The accurate crop yield predictions from the analysis for WW and OSR could be further used for the application of biodiversity, where the impact of land use diversity on crop yields could be estimated.The ease of using spatiotemporal modelling with crop growth models would not be limited to Bavaria.The study's methodology could also be tested by coupling machine/deep learning (ML/DL) approaches with CGMs, which might help to include more climate elements to achieve more precise results.Lastly, the study's two decades of accurate yield estimations could strengthen trust in the decision(/policy) making to achieve sustainability in agriculture.(p) (q) (r) (s)

Figure 1 .
Figure 1.The conceptual framework of the study is divided into three parts: Part 1 states the data fusion for 2019 to investigate the best synthetic NDVI time series product (this section was already completed in our previous study [59]); Part 2 generates and validates the synthetic NDVI time series from 2001 to 2019 for the product L-MOD13Q1; and Part 3 performs the comparative analysis to compare the performance of fused (L-MOD13Q1) and non-fused (MOD13Q1) NDVI time series in crop yield prediction for 2019 and then estimates and validates the crop yield for Bavaria by inputting the L-MOD13Q1 time series and climate elements to a semi-empiric Light Use Efficiency (LUE) model; STARFM = Spatial and Temporal Adaptive Reflectance Fusion Model; NDVI = Normalised Difference Vegetation Index; L-MOD09GQ = Landsat-MOD09GQ; L-MOD09Q1 = Landsat-MOD09Q1; L-MCD43A4 = Landsat-MCD43A4; L-MOD13Q1 = Landsat-MOD13Q1; S-MOD09GQ = Sentinel-2-MOD09GQ; S-MOD09Q1 = Sentinel-2-MOD09Q1; S-MCD43A4 = Sentinel-2-MCD43A4; S-MOD13Q1 = Sentinel-2-MOD13Q1; PAR is photosynthetically active radiation, and FPAR is the fraction of PAR absorbed by the canopy.APAR = Absorbed Photosynthetically Active Radiation.

Figure 1 .
Figure 1.The conceptual framework of the study is divided into three parts: Part 1 states the data fusion for 2019 to investigate the best synthetic NDVI time series product (this section was already completed in our previous study [59]); Part 2 generates and validates the synthetic NDVI time series from 2001 to 2019 for the product L-MOD13Q1; and Part 3 performs the comparative analysis to compare the performance of fused (L-MOD13Q1) and non-fused (MOD13Q1) NDVI time series in crop yield prediction for 2019 and then estimates and validates the crop yield for Bavaria by inputting the L-MOD13Q1 time series and climate elements to a semi-empiric Light Use Efficiency (LUE) model; STARFM = Spatial and Temporal Adaptive Reflectance Fusion Model; NDVI = Normalised Difference Vegetation Index; L-MOD09GQ = Landsat-MOD09GQ; L-MOD09Q1 = Landsat-MOD09Q1; L-MCD43A4 = Landsat-MCD43A4; L-MOD13Q1 = Landsat-MOD13Q1; S-MOD09GQ = Sentinel-2-MOD09GQ; S-MOD09Q1 = Sentinel-2-MOD09Q1; S-MCD43A4 = Sentinel-2-MCD43A4; S-MOD13Q1 = Sentinel-2-MOD13Q1; PAR is photosynthetically active radiation, and FPAR is the fraction of PAR absorbed by the canopy.APAR = Absorbed Photosynthetically Active Radiation.

Figure 2 .
Figure 2. Overview of the study region.The LC map of Bavaria is obtained by combining multiple inputs of landcover maps, such as the Amtliche Topographisch-Kartographische Informations System, Integrated Administration Control System (which provides the crop field information), and the Corine LC, into one map.Agriculture (peach green) dominates mainly in the northwest and southeast of Bavaria, while forest and grassland classes (dark green and yellow, respectively) dom-

Figure 2 .
Figure 2. Overview of the study region.The LC map of Bavaria is obtained by combining multiple inputs of landcover maps, such as the Amtliche Topographisch-Kartographische Informations System, Integrated Administration Control System (which provides the crop field information), and the Corine LC, into one map.Agriculture (peach green) dominates mainly in the northwest and southeast of

Figure 3 .
Figure 3.The scatter plots (a-s) compare the accuracies of Landsat (referenced NDVI) with L-MOD13Q1 (synthetic NDVI) for 2001 to 2019.values of the statistical parameters such as R 2 and RMSE and the total number of Landsat scenes available every year (N) are displayed at the to each plot.Every plot contains a solid line (1:1 line) that is used to visualise the correlation of pixels between the referenced and synthetic N values.The dashed line represents the regression line.The colour of scatter plots depicts the density of points (yellow: low, blue: high).

Figure 3 .
Figure 3.The scatter plots (a-s) compare the accuracies of Landsat (referenced NDVI) with L-MOD13Q1 (synthetic NDVI) for 2001 to 2019.The values of the statistical parameters such as R 2 and RMSE and the total number of Landsat scenes available every year (N) are displayed at the top of each plot.Every plot contains a solid line (1:1 line) that is used to visualise the correlation of pixels between the referenced and synthetic NDVI values.The dashed line represents the regression line.The colour of scatter plots depicts the density of points (yellow: low, blue: high).

Figure 4 .Figure 5 .
Figure 4.The correlation plots between the total number of Landsat scenes per year (N) and (a) R 2 values and (b) RMSE values obtained during the accuracy assessment of referenced and synthetic NDVI products from 2001 to 2019.The correlation coefficient refers to R (see Equation (5)).

Figure
Figure 6a-c displayed the crop yield accuracies between the modelled and referenced crop yields of WW and OSR obtained with different satellite products using the LUE model in 2019.The figures show that the fused product (L-MOD13Q1) obtained a higher R 2 (0.81) and a lower RMSE (3.91 dt/ha) than the non-fused product (MOD13Q1: R 2 = 0.70 and RMSE = 4.77 dt/ha) for both WW and OSR.Analysing the ME of both products with LUE, the L-MOD13Q1 resulted in a lower ME (3.04 dt/ha) than the MOD13Q1 (3.50 dt/ha) (Figure 6c).

Figure 4 .
Figure 4.The correlation plots between the total number of Landsat scenes per year (N) and (a) R 2 values and (b) RMSE values obtained during the accuracy assessment of referenced and synthetic NDVI products from 2001 to 2019.The correlation coefficient refers to R (see Equation (5)).On comparing the yearly fusion results on a DOY basis, the DOYs 113, 129, and 193 had the highest average accuracy with an R 2 of more than 0.65 and a RMSE lesser than 0.10 (Figure5a,b).The DOYs of 33 to 97 and 145 to 177, with low R 2 (<0.60) and high RMSE (>0.11), were obtained.

Figure 4 .Figure 5 .
Figure 4.The correlation plots between the total number of Landsat scenes per year (N) and (a) R 2 values and (b) RMSE values obtained during the accuracy assessment of referenced and synthetic NDVI products from 2001 to 2019.The correlation coefficient refers to R (see Equation (5)).
Figure 6a-c displayed the crop yield accuracies between the modelled and referenced crop yields of WW and OSR obtained with different satellite products using the LUE model in 2019.The figures show that the fused product (L-MOD13Q1) obtained a higher R 2 (0.81) and a lower RMSE (3.91 dt/ha) than the non-fused product (MOD13Q1: R 2 = 0.70 and RMSE = 4.77 dt/ha) for both WW and OSR.Analysing the ME of both products with LUE, the L-MOD13Q1 resulted in a lower ME (3.04 dt/ha) than the MOD13Q1 (3.50 dt/ha) (Figure 6c).

Figure 5 .
Figure 5.The day of the year (DOY)-based comparison of correlation coefficients between (a) R 2 values and (b) RMSE values obtained during the accuracy assessment of referenced and synthetic NDVI products from 2001 to 2019.The correlation coefficient refers to R (see Equation (5)).

Figure 6a -Figure 6 .
Figure 6a-c displayed the crop yield accuracies between the modelled and referenced crop yields of WW and OSR obtained with different satellite products using the LUE model in 2019.The figures show that the fused product (L-MOD13Q1) obtained a higher R 2 (0.81) and a lower RMSE (3.91 dt/ha) than the non-fused product (MOD13Q1: R 2 = 0.70 and RMSE = 4.77 dt/ha) for both WW and OSR.Analysing the ME of both products with LUE, the L-MOD13Q1 resulted in a lower ME (3.04 dt/ha) than the MOD13Q1 (3.50 dt/ha) (Figure 6c).

Figure 6 .
Figure 6.The dot plots compare the accuracies (a) R 2 , (b) RMSE, and (c) ME of referenced data (at 95% confidence intervals) and modelled yields obtained from multi-source data: MOD13Q1 and L-MOD13Q1 in 2019.

Figure 7 .Figure 8 .
Figure 7.The scatter plots compare the accuracies of the modelled and referenced yields (at a 95% confidence interval) of (a) WW and (b) OSR for 19 years together (i.e., from 2001 to 2019).The values of the statistical parameters such as R 2 , RMSE (dt/ha), ME (dt/ha), and total number of points (n) are displayed at the top of each plot.Every plot contains a solid line (1:1 line) that is used to visualise the correlation of pixels between the modelled and referenced yield values.The dashed line represents the regression line.Different colours of the points display different years.

Figure 7 .Figure 7 .Figure 8 .
Figure 7.The scatter plots compare the accuracies of the modelled and referenced yields (at a 95% confidence interval) of (a) WW and (b) OSR for 19 years together (i.e., from 2001 to 2019).The values of the statistical parameters such as R 2 , RMSE (dt/ha), ME (dt/ha), and total number of points (n) are displayed at the top of each plot.Every plot contains a solid line (1:1 line) that is used to visualise the correlation of pixels between the modelled and referenced yield values.The dashed line represents the regression line.Different colours of the points display different years.

Figure 8 .
Figure 8.The scatter plots correlating the modelled yield and regional mean elevation for (a) WW and (b) OSR.The dashed line represents the regression line.Different colours of the points display different crop types (green for WW and orange for OSR).The correlation coefficient refers to R (see Equation (5)).

Figure 9 .
Figure 9.The bar plots show the yearly comparison of accuracies (a) R 2 values and (b) RMSE values obtained from the referenced yields (at a 95% confidence interval), with LUE-modelled yields including climate stress factors (dark blue) and LUE-modelled yields excluding the climate stress factors (sensitivity analysis) (light blue).The scatter plots compare the accuracies of the modelled and referenced yields (at a 95% confidence interval) of (c) WW and (d) OSR for 19 years together (i.e., from 2001 to 2019).The values of the statistical parameters such as R 2 , RMSE (dt/ha), ME (dt/ha), and total number of points (n) are displayed at the top of each plot.Every plot contains a solid line (1:1 line) that is used to visualise the correlation of pixels between the modelled and referenced yield values.The dashed line represents the regression line.Different colours of the points display different years.

Figure 9 .
Figure 9.The bar plots show the yearly comparison of accuracies (a) R 2 values and (b) RMSE values obtained from the referenced yields (at a 95% confidence interval), with LUE-modelled yields including climate stress factors (dark blue) and LUE-modelled yields excluding the climate stress factors (sensitivity analysis) (light blue).The scatter plots compare the accuracies of the modelled and referenced yields (at a 95% confidence interval) of (c) WW and (d) OSR for 19 years together (i.e., from 2001 to 2019).The values of the statistical parameters such as R 2 , RMSE (dt/ha), ME (dt/ha), and total number of points (n) are displayed at the top of each plot.Every plot contains a solid line (1:1 line) that is used to visualise the correlation of pixels between the modelled and referenced yield values.The dashed line represents the regression line.Different colours of the points display different years.

Figure 10 .
Figure 10.Spatial distribution of mean referenced yield (2001-2019) and the year-wise predicted yield for WW from 2001 to 2019 using the LUE model for the state of Bavaria.The white colour represents no available data.A detailed map of the administrative regions of Bavaria is shown in Figure A1.

Figure 10 .
Figure 10.Spatial distribution of mean referenced yield (2001-2019) and the year-wise predicted yield for WW from 2001 to 2019 using the LUE model for the state of Bavaria.The white colour represents no available data.A detailed map of the administrative regions of Bavaria is shown in Figure A1.

Figure 11 .Figure 11 .
Figure 11.Spatial distribution of mean referenced yield (2001-2019) and the year-wise predicted yield for OSR from 2001 to 2019 using the LUE model for the state of Bavaria.The white colour represents no available data.A detailed map of the administrative regions of Bavaria is shown in Figure A1.

Figure 11 .
Figure 11.Spatial distribution of mean referenced yield (2001-2019) and the year-wise predicted yield for OSR from 2001 to 2019 using the LUE model for the state of Bavaria.The white colour represents no available data.A detailed map of the administrative regions of Bavaria is shown in Figure A1.

Figure 12 .
Figure 12.The dot plots show the district-wise distribution of modelled yield for (a) WW and (b) OSR, from 2001 to 2019.The green colour depicts the modelled yield of WW, the orange colour depicts the modelled yield of OSR, and the grey colour depicts both referenced yields of WW and OSR.

Figure 12 .Figure 12 .
Figure 12.The dot plots show the district-wise distribution of modelled yield for (a) WW and (b) OSR, from 2001 to 2019.The green colour depicts the modelled yield of WW, the orange colour depicts the modelled yield of OSR, and the grey colour depicts both referenced yields of WW and OSR.

Figure 13 .
Figure 13.The line plots compare the accuracies with the mean yield percent difference (as calculated in Equation (9)) for WW and OSR for 19 years (i.e., from 2001 to 2019).The accuracies of WW and OSR are analysed in six categories (less than −4, −4 to −2, −2 to 0, 0 to 2, 2 to 4, and more than 4%) of yield percent difference.The negative range shows the overestimation, and the positive range shows the underestimation of the modelled yield values by the LUE compared to the referenced yield values.The green colour depicts WW, and the orange colour depicts OSR.

Figure 13 .
Figure 13.The line plots compare the accuracies with the mean yield percent difference (as calculated in Equation (9)) for WW and OSR for 19 years (i.e., from 2001 to 2019).The accuracies of WW and OSR are analysed in six categories (less than −4, −4 to −2, −2 to 0, 0 to 2, 2 to 4, and more than 4%) of yield percent difference.The negative range shows the overestimation, and the positive range shows the underestimation of the modelled yield values by the LUE compared to the referenced yield values.The green colour depicts WW, and the orange colour depicts OSR.

Figure 14 .Figure 15 .
Figure 14.The bar plots compare the yearly (a) R 2 and (b) RMSE values of estimated OSR yield (orange), WW yield (green), and synthetic NDVI (purple) from 2001 to 2019.The units of the RMSE values of both WW and OSR yields are dt/ha.

Figure 14 .
Figure 14.The bar plots compare the yearly (a) R 2 and (b) RMSE values of estimated OSR yield (orange), WW yield (green), and synthetic NDVI (purple) from 2001 to 2019.The units of the RMSE values of both WW and OSR yields are dt/ha.

Figure 14 .Figure 15 .
Figure 14.The bar plots compare the yearly (a) R 2 and (b) RMSE values of estimated OSR yield (orange), WW yield (green), and synthetic NDVI (purple) from 2001 to 2019.The units of the RMSE values of both WW and OSR yields are dt/ha.

Figure 15 .
Figure 15.The correlation plots between R 2 of synthetic NDVI time series and R 2 of modelled yield time series for (a) WW (green) and (b) OSR (orange), from 2001 to 2019.The correlation coefficient refers to R (see Equation (5)).

Figure 16 .
Figure 16.The side-by-side visualisation of synthetic NDVI products obtained on 18 June 2005, 2013 and 2019 (left) with the WW biomass obtained from the LUE modelled for the years of 2005, 2013 and 2019 (right).

Figure 16 .
Figure 16.The side-by-side visualisation of synthetic NDVI products obtained on 18 June 2005, 2013 and 2019 (left) with the WW biomass obtained from the LUE modelled for the years of 2005, 2013 and 2019 (right).

1 .
Quality Assessment of Synthetic Remote Sensing Time Series from 2001 to 2019 (i) To find the potential of STARFM for long-term time series, the paper generates and validates a synthetic normalised difference vegetation index (NDVI) time series blending the high spatial resolution (30 m, 16 days) of Landsat 5 Thematic Mapper (TM) (2001 to 2012), Landsat 7 Enhanced Thematic Mapper Plus (ETM+) (2012), and Landsat 8 Operational Land Imager (OLI) (2013 to 2019) with the coarse resolution of MOD13Q1 (250 m, 16 days) from 2001 to 2019.Overall, the average accuracy of data fusion for nineteen years has an R 2 of 0.66 and an RMSE of 0.11.The accuracy of data fusion is found to be dependent on the number of Landsat scenes available per year (N).The higher the N, the more accurate is the synthetic NDVI time series per year.(ii) To investigate the stability and precision of the LUE model in crop yield prediction, the paper inputs the synthetic NDVI time series and climate elements to the crop model to estimate and validate yearly crop yields for WW and OSR from 2001 to 2019.The validation of crop yield at regional scale results in an average R 2 of 0.79 (WW)/0.86(OSR) and an RMSE of 4.51 dt/ha/2.46dt/ha, respectively.(iii) Identifying the impact of the input data fusion product on the accuracy assessment of the LUE model, high positive correlations are seen when the accuracies of the synthetic NDVI time series are plotted with the accuracies of modelled crop yield from 2001 to 2019 for WW (R = 0.81) and OSR (0.77).
Figure A2.The digital elevation map of Bavaria.The map is generated from shuttle radar topography mission (SRTM) digital elevation data.The elevation ranges from 93 m to 2943 m.

Figure A2 .
Figure A2.The digital elevation map of Bavaria.The map is generated from shuttle radar topography mission (SRTM) digital elevation data.The elevation ranges from 93 m to 2943 m.

Figure A3 .Figure A3 .Figure A3 .Figure
Figure A3.The scatter plots (a-s) compare the accuracies of modelled and referenced yields of WW for 2001 to 2019.The values of the statistical p rameters such as R 2 , RMSE (dt/ha), and ME (dt/ha) are displayed at the top of each plot.Every plot contains a solid line (1:1 line) that is used to v ualise the correlation of pixels between the modelled and referenced yield values.The green colour of scatter plots represents WW.

Figure A4 .
Figure A4.The scatter plots (a-s) compare the accuracies of modelled and referenced yields of OSR for 2001 to 2019.The values of the statistical parameters such as R 2 , RMSE (dt/ha), and ME (dt/ha) are displayed at the top of each plot.Every plot contains a solid line (1:1 line) that is used to visualise the correlation of pixels between the modelled and referenced yield values.The orange colour of scatter plots represents OSR.

Figure A4 .Figure A5 .
Figure A4.The scatter plots (a-s) compare the accuracies of modelled and referenced yields of OSR for 2001 to 2019.The values of the statistical parameters such as R 2 , RMSE (dt/ha), and ME (dt/ha) are displayed at the top of each plot.Every plot contains a solid line (1:1 line) that is used to visualise the correlation of pixels between the modelled and referenced yield values.The orange colour of scatter plots represents OSR.

Table 2 .
A summary of the collected cloud-and shadow-free Landsat 5, Landsat 7, and Landsat 8 datasets available every year with their day of the years (DOYs) between the start and end of the seasons of WW and OSR from 2001 to 2019.N is the total number of Landsat scenes available per year for WW and OSR.