A Low-Complexity, County-Scale Yield Prediction Method for Maize and Sunflower Using Harmonized Landsat–Sentinel (HLS) Data

1. Introduction

As the world population increases and environmental stresses intensify, the reliable prediction of crop production has become both an economic necessity and a strategic priority [1,2]. Timely and accurate crop yield estimation at global, national, and regional scales is essential for modern agricultural management, ensuring global food security, balancing international grain trade, and promoting sustainable agricultural development [3,4,5].

The application of remote sensing technologies in agriculture began with the launch of the Landsat satellite, equipped with the Multispectral Scanner (MSS) in 1972 [6]. Since then, remote sensing has become a widely used technique in agriculture [7] and has played a significant role in assessing crop production over the years [8]. Satellite-based remote sensing technology has rapidly developed into a widely applied tool for crop yield estimation [9,10]. For large-scale yield monitoring, it is necessary that remotely sensed data be available reliably and consistently over several years, and that the temporal resolution of the data be sufficiently dense. The combination of high temporal frequency, wide spatial coverage, and low cost makes remote sensing an ideal choice for national- and regional-level crop yield estimation applications [11].

In previous decades, this requirement was met, for example, by data from the AVHRR (Advanced Very High Resolution Radiometer, ~1 km spatial resolution) instrument aboard the NOAA (National Oceanic and Atmospheric Administration) satellites, or from the MODIS (Moderate Resolution Imaging Spectroradiometer, 250–500 m spatial resolution) sensors on the Terra and Aqua satellites. With the advent of the new-generation multispectral Sentinel-2 satellite constellation (Sentinel-2A, launched in 2015; Sentinel-2B, launched in 2017; and Sentinel-2C, launched in 2024) by the European Space Agency new opportunities have emerged for long-term, frequent monitoring of agricultural landscapes [12]. The Multispectral Imager (MSI) onboard all three satellites provides high-resolution optical observations in 13 spectral bands with a spatial resolution between 10 m and 60 m, and a revisit frequency of 2–5 days when two satellites are operational [13]. Around the same time, in 2013, the National Aeronautics and Space Administration (NASA) launched Landsat-8, the new member of the more than 40-year-old Landsat series [14], followed by Landsat-9 in 2021 [15]. Both satellites carry the Operational Land Imager (OLI) instrument, with a spatial resolution of 30 m (15 m for panchromatic) and a revisit time of 16 days [16]. Landsat-9 has the important mission of working together with Landsat-8 to reduce the revisit period of Landsat Earth observations to eight days [17].

It has been well established that combining data from different satellite sensors provides greater opportunities for cloud-free surface observation [18,19,20]. Given that Landsat-OLI and Sentinel-MSI make similar measurements in terms of spectral, spatial, and angular characteristics [12,21], it was a natural step to integrate data from the two instruments into a unified system. This project, known as the Harmonized Landsat/Sentinel-2 (HLS), provides high-frequency global imagery products (excluding Antarctica) with a spatial resolution of 30 m every two to three days since the launch of Landsat-9 [22]. Such dense and consistent data are critical for various applications, including continuous crop monitoring and yield estimation [23].

Remote sensing–based yield estimation models are generally categorized into three main types: data assimilation models, light-use efficiency models, and empirical statistical models [24]. Empirical statistical models do not require numerous input parameters [25]; therefore, their use is widespread due to their relative simplicity, rapidity, and computational ease [26,27]. In recent years, these conventional models have produced remarkable results in yield estimation [28,29]. Such models mainly rely on linear relationships between crop yield and characteristic spectral bands or vegetation indices , sometimes supplemented by basic meteorological or agronomic data [30,31,32,33,34,35]. Here, observations from the past are used to inform what is happening in the present, without a strong need for understanding of the causality [36]. In these studies, the most commonly used vegetation index is the Normalized Difference Vegetation Index (NDVI) [37,38], which is strongly correlated with photosynthetic capacity, and thus yield [36].

Building on a previously developed robust yield estimation method based on a regression method, which has been successfully applied to several crop species [39,40], the present study aims to further refine and extend this approach using a simple empirical regression yield estimation method for maize and sunflower based on HLS NDVI data. To achieve this, we addressed the following main goals: (i) improve the previously developed robust method based on 30-m data, (ii) test the model on input data from areas with different vegetation cover, (iii) estimate the predictive power of the model, and (iv) test the accuracy of the model for forecasts at different dates.

2. Materials and Methods

2.1. Study Area and Crop Information

The study area was Fejér County (corresponding to the NUTS (Nomenclature of Territorial Units for Statistics)-3 level), one of the 19 counties of Hungary (Figure 1), located in Central Europe. The county is situated in central Hungary and has a continental climate characterized by relatively warm summers and moderately cold winters. During the study period 2017–2024, the mean annual temperature was 12.3 °C, while the mean summer temperature exceeded 16 °C. The average annual precipitation ranged between 500 and 600 mm [41]. Fejér County covers 4373 km², of which approximately 3000 km² is agricultural land.

In the present study, the two spring-sown row crops with the largest sown areas—maize and sunflower—were investigated. In the 2017–2024 period, the average sown area of maize was 749 km² and that of sunflower was 385 km². The yield data used to validate the yield estimation method (Figure 2) were obtained from the official reports of the HCSO (Hungarian Central Statistical Office) [42].

2.2. Land Cover and Crop Information

The crop yield estimation analysis was performed for agricultural pixels identified using species and habitat type information from the Ecosystem Map of Hungary developed within the NÖSZTÉP (Nemzeti Ökoszisztéma Szolgáltatás-Térképezés és Értékelés) project [43,44]. This map (hereafter referred to as the NÖSZTÉP dataset) distinguishes 56 ecosystem categories (Level 3) at a spatial resolution of 20 m and covers the entire territory of Hungary. The dataset was resampled to the grid of the remote sensing data by summing up the individual shares, resulting in fractional cover values (in percent) for each ecosystem category in every pixel.

To refine our study area and make our method crop-specific, we incorporated crop maps from 2017 to 2024. These maps were provided by the Lechner Knowledge Center [45], and were produced using machine learning methods applied to Sentinel-2 and Sentinel-1 time series with a resolution of 20 meters. They provide detailed information about crop types. Based on these maps, we identified and selected only maize and sunflower pixels for the study. These crops had a precision of around 95% during the study period [46]. Despite using different methods, the areas of the main crops showed good agreement with the values reported by the HCSO [47] (Figure 3). To minimize error, we applied a 20-meter negative buffer to exclude border pixels that could affect our analysis with undesirable edge effects. Then, we resampled the maps for the HLS grid (30 m) and used only those with at least 90% coverage.

2.3. Processing of the HLS Dataset / Calculating County Mean NDVI from HLS Data

The Harmonized Landsat-Sentinel (HLS v2.0) dataset [48,49] was used to calculate NDVI at 30 m spatial resolution for the period 2017–2024. This dataset combines measurements from the OLI on Landsat-8 and Landsat-9 and the MSI on Sentinel-2A and 2B, providing Nadir BRDF (bidirectional reflectance distribution function)-Adjusted Reflectance (NBAR) data at 30 m resolution. Both HLSL30 (Landsat-based) and HLSS30 (Sentinel-based) datasets were downloaded from NASA EarthData Search for the tiles T33TYN and T33TYM in UTM (Universal Transverse Mercator) zone 33, resulting in a total of 2980 granules for the eight years [50].

Pixels were classified as agricultural when the share of agricultural area exceeded 90%. From these data, daily mean NDVI values for Fejér County were calculated separately for Landsat (HLSL30) and Sentinel (HLSS30) datasets, using only the best-quality observations (pixels with low or moderate aerosol content). The calculations were performed for all NÖSZTÉP-based agricultural pixels and, separately, for maize and sunflower fields according to the crop-type information available on the same grid.

2.4. Crop Yield Estimation

In this study, two yield estimation models were developed for Fejér County, covering eight years (2017–2024), based on HLS satellite data. In both cases, the calculations were derived from daily NDVI values averaged over a given area. The two models differ in that, in one case, the data originated from the entire agricultural area of the county, whereas in the other case, a crop-specific mask was applied. The advantage of the first method is that it does not require detailed information on the sown area of each crop, which may vary from year to year. In this approach, it is assumed that the general conditions of vegetation-covered surfaces also represent the conditions of individual crops with similar vegetation cycles [40,51]. The second model, by contrast, requires detailed classification each year to determine the exact sown area of each crop. The first model was previously referred to as the robust method [39], and this terminology is retained in the present study. The second model is hereafter referred to as the crop-specific (CS) model. In the robust procedure, the entire agricultural area of the county was considered, while clearly non-agricultural areas (settlements, forests, large water bodies, etc.) were masked out. The robust model thus minimises data requirements, whereas the CS model incorporates more detailed spatial information but requires accurate annual classification. Subsequently, annual NDVI time series were generated for each crop from the averaged daily NDVI data. As an illustration, Figure 4 presents the NDVI time series for 2024.

As the next step, a double-Gaussian curve was fitted to the annual NDVI time series between appropriately selected starting (t_b) and ending (t_e) days, depending on the crop type, to eliminate stochastic noise caused by atmospheric and other effects [30]:

$$G\left(t\right)=A_1·e^{{\displaystyle \frac{{-(t-t_1)}^2}{h_1^2}} }{+ A}_2·e^{{\displaystyle \frac{{-(t-t_2)}^2}{h_2^2}} }+S,$$

(1)

where G(t) is the fitted double-Gaussian curve, S is the predefined value of the so-called soil line, A₁ and A₂ are the amplitudes, t₁ and t₂ are the locations of the maxima of the Gaussian curves, and h₁ and h₂ are their respective half-widths. The value of t_e corresponds to the date of yield estimation.

The values of the resulting double-Gaussian curve were then integrated between a suitably selected starting (d_b) and ending (d_e) day, thus obtaining an index that can be directly compared with the yield data published by the HCSO. It is important to note that these dates (d_b and d_e) are distinct from the curve-fitting parameters (t_b and t_e). While the latter are defined based on the beginning and end of the vegetation period of the given crop, the d_b and d_e values were determined to maximize the correlation between the estimated and actual yield data (provided by the HCSO). The indices obtained using the CS method were named GYURI (General Yield Reference Index), and those obtained using the robust procedure were named GYURRI (General Yield Robust Reference Index) [39]. For example, in the case of GYURI, the value for a given year y is defined as follows:

$${ GYURI}_y={\sum }_{t=d_b}^{d_e}G\left(t\right),$$

(2)

where G(t) is the double-Gaussian function fitted to the daily NDVI data measured in year y for the given crop, and d_b and d_e are the starting and ending dates of integration. Then, a straight line was fitted to the annual GYURI_y (or GYURRI_y) – HCSO point pairs, from whose parameters the yield for the given year can be simply calculated.

2.5. Statistical Analysis, Diagnostic and Predictive Approach

In this study, we distinguished between diagnostic and predictive applications of the two above-mentioned yield estimation models. The diagnostic model aims to describe the relationship between vegetation index–based indicators (e.g., GYURI) and the observed yields in the eight years within this study as accurately as possible. The predictive model, in contrast, aims to estimate future yields under operational conditions, where the model performance is evaluated by cross-validation (leave-one-year-out) [23,35,52]. This distinction separates the interpretive value of the method (how yield fluctuations relate to vegetation dynamics) from its forecasting capability (how well future yields can be predicted). In the following, we refer to the diagnostic method as R-D for robust model and CS-D for crop-specific model, while the predictive method as R-P and CS-P.

The efficiency of the models was tested with three statistical indicators: square of the Pearson linear correlation coefficient (R²), mean absolute error (MAE), and mean absolute percentage error (MAPE). Calculating the correlation coefficient is one of the most commonly used methods in statistical studies [53], MAE is a widely used metric for evaluating models, similarly to root mean square error (RMSE) [54,55], while MAPE is a measure of quality for regression models [56]:

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(y_i-{\mathrm{ŷ}}_i)}^2}{{\sum }_{i=1}^n{(y_i-{\mathrm{y}}_{\mathrm{a}})}^2} ,}$$

(3)

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{(y_i-{\mathrm{ŷ}}_i)}^2} ,$$

(4)

$$MAE={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n\left|y_i-{\mathrm{ŷ}}_i\right| ,$$

(5)

$$MAPE={\displaystyle \frac{100}{n}}{\sum }_{i=1}^n\left|{\displaystyle \frac{y_i-{\mathrm{ŷ}}_i}{y_i}}\right| ,$$

(6)

where

$$y_a={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{y}_i$$

(7)

and n denotes the number of observations, y_i is the observed yield, and ŷ_iis the estimated yield.

3. Results

3.1. The Model Parameters for Maize

One of the most important crops in Hungary, and thus in Fejér county, is maize, whose sown area accounts for approximately 15–20 percent of the total area of the county (Figure 2). The main growing stages for maize in Hungary are the seeding stage (mid-April to early May), the jointing stage (from mid-May to early June), the heading stage (from late June to late July), the milking stage (from late July to mid-August), and the mature stage (from mid-September to early October) [57]. The harvest usually ends around late September–mid-October. The daily NDVI values of the maize fields determined based on classification start to increase around DOY (Day Of Year) 120 (end of April), and then, after reaching maximum values, begin to decrease around DOY 210–220. An exception occurred in 2022, when the extreme late-summer drought caused an abrupt decline. (Figure 5). By the end of the harvest (around DOY 260–280), the NDVI values decline to around 0.3. Based on these patterns, the value of soil-line (S) was set to 0.3, and t_b = 125 and t_e = 260 were chosen as the starting and final days of curve fitting. This also means that the yield estimate was made on DOY 260, which typically precedes the end of the harvest period.

As mentioned in the previous chapter, in this study, four yield forecasting methods were tested: the R-D, CS-D, R-P, and CS-P models. The R and CS methods are completely similar; they only differ in the areas from which the NDVI data originate. Using the robust (R) procedure, these data are averages for the entire agricultural area, while in the CS method, they are averages only for areas classified as maize. Consequently, in the R method, the influence of autumn cereals is obvious in the spring–early summer data (Figure 6), so we examined how the accuracy of the estimate changes if the starting day of curve fitting (t_b) increases. Experience has shown that it gives better results if the value of t_b is chosen as 140, so this value was used in the case of the R method. The starting and ending days of the integration (d_b and d_e) used in the calculation of GYURRI in the R-D procedure were 135 and 250, respectively, while in the CS-D model, the best result was obtained by using d_b = 160 and d_e = 235.

3.2. Results for Maize in the Diagnostic Model

In the eight years studied, maize yields ranged widely, with the lowest yield in 2022 (2.56 t/ha), which was hit by a late summer drought, and the highest in 2018 (10.06 t/ha). This wide yield range improved the robustness of the model calibration, resulting in very high correlation values, R² = 0.986 for the R-D method, and R² = 0.998 for the CS-D method (Figure 7).

The differences between the estimated and actual HCSO yield data (Table 1) show that in the case of the R-D method, the difference exceeded 5% in three of the eight years examined, the maximum error (6.25%) occurring in the drought year of 2022. Applying the CS-D model, the differences remained well below 5% in each year, with the maximum difference being 2.33% in 2024. The CS-D model also performed reliably under extreme weather conditions: in 2022 the error of the estimation was only 1.56%. The MAE, that is the average absolute yield difference for the eight years, is 0.23 t/ha for the R-D model and 0.09 t/ha for the CS-D model, and the MAPE is 3.78% and 1.37%, respectively.

3.3. Results for Maize in the Predictive Model

The results of the diagnostic model are promising, but it should not be overlooked that the estimates were calculated based on only eight data points. Therefore, it is necessary to examine how the results change if one year is always omitted from the determination of the model parameters, and the yield of this omitted year is estimated using parameters calculated based on the remaining seven years. It would be possible to vary the values of d_b and d_e in each case to achieve the best correlation values, but such adjustments would compromise model stability. Therefore, the same fixed values of d_b and d_e were used for each year as in the diagnostic model, and only the parameters of the fitted lines changed from year to year.

The results show (Table 2) that for the R-P model, the MAE was 0.3 t/ha, while the MAPE value increased to over 5%. The drought year of 2022 contributed greatly to this, when the percentage error of the estimate exceeded 16%. In contrast, the CS-P model yielded substantially better performance. Here, the largest difference between the estimated and actual yield was 0.17 t/ha in 2019, and even in the drought year of 2020, the percentage difference did not exceed 5% (4.42%). For the CS-P model, the MAE was 0.11 t/ha, while the MAPE value was 2%. The correlation value (R²) of the R-P and CS-P model was 0.976 and 0.997, respectively (Figure 8). The estimation results for each year, using the four models, are summarized in Figure 9.

As mentioned above, in predictive models, one year was always omitted when determining the model parameters. In the case of the CS-P model, a simulation of an operational yield estimation situation was performed. In this procedure, the 2020 crop yield was first estimated from the GYURI data using parameters determined only based on the first three years (2017–2019). Then, the parameters were determined based on four years (2017–2020), which would be used to estimate the yield for 2021, and so on. Obviously, the estimate for 2024 was the same as the result obtained for the CS-P model.

The results (Table 3) are encouraging; only the year 2022, which had an extremely low yield due to drought, had a significant error exceeding 10%. It can also be seen that after including this year, the parameters of the fitted line (a and b) changed significantly. This also supports the experience that for building a sufficiently accurate model, it is fortunate that the difference between the minimum and maximum yield is large enough. While the difference between the highest and lowest yield in 2017–2021 was 4.18 t/ha, this increased to 7.5 t/ha after including the year 2022. Interestingly, in the CS-P procedure, when the 2023–2024 data were included in addition to the 2017–2021 data for the 2022 estimate, the error, as mentioned above, turned out to be significantly smaller. One reason for this is that although the lines fitted to five or seven years with higher yields differ only slightly, in the low-yield year of 2022, even this small difference causes a significant percentage difference (Figure 10). This also indicates that model stability is expected to increase as longer time series become available.

The above results were obtained with an estimate made on DOY 260. However, examining the evolution of the correlation coefficient (R²) obtained by changing the t_e values in the case of the CS-D model (Figure 11), it appears that it is possible to provide a sufficiently accurate forecast even at an earlier date. The calculations were performed for the R-P and CS-P model by simply bringing the last day of curve fitting (t_e) forward.

The results showed that the CS-P model obtained significantly better results than the R-P model (Figure 12). Using the CS-P model, the MAPE remained below 4% (3.86% to be exact) until DOY 230, about 30–50 days before the end of harvest. Interestingly, in the case of the CS-P model, the predicted results on DOY 250 were even slightly more accurate than the estimates obtained on DOY 260. Table 4c shows the numerical results of the forecast.

3.4. The Model Parameters for Sunflower

Sunflower is the second most important spring crop in Hungary and Fejér county, with a sowing area of about half that of maize (Figure 3). Annual sunflower and maize yields showed a moderate positive correlation (Figure 2), but in the case of sunflower, the fluctuations in yields were much smaller during the period under study. While the largest positive and negative deviations from the eight-year average yield for maize were +46% (2018) and -63% (2022), for sunflower, they were only +17% (2019) and –26% (2022). This indicates that the severe drought in 2022 affected sunflower yields to a lesser extent. The growing season of sunflower is roughly the same as that of maize; the NDVI time series from classified sunflower fields are similar to those observed for maize, although the autumn decreasing phase of the values begins earlier (Figure 13). Sunflower harvest generally begins in late August and lasts until mid-September in Hungary. The timing depends largely on the weather of the given year; in drier and warmer years, the harvest season may be shorter, while in cooler, wetter years, it may be slightly longer depending on the time of sowing and the course of ripening. Based on these characteristics, the expectation was that the yield estimation date could be advanced compared to maize.

3.5. Results for Sunflower in the Diagnostic Model

The starting day for fitting the double-Gaussian curve was chosen as t_b = 125, similar to maize. The final date of the curve-fitting was determined based on the correlation values obtained in the CS-D model. The highest correlation (R² = 0.903) was obtained at t_e = 230 (Figure 14), which also means that the yield estimate can be made with sufficient accuracy already in mid-August (around 18 August).

Given the smaller sown area and lower interannual variability, a slightly weaker correlation was expected than in the case of maize. The results confirmed this expectation, although the correlation was above 0.9 in both the R-D and CS-D cases. Unexpectedly, the R-D model outperformed the CS-D model. The correlation between the estimated and actual yields for the R-D method was R² = 0.954, while for the CS-D model it was R² = 0.903 (Figure 15). The MAE for the R-D method was 0.06 t/ha (MAPE = 2.08%), and for the CS-D method, it was 0.08 t/ha (MAPE = 2.90%) (Table 5). In the case of the R-D method, the estimation error remained below 4% in all years, while in the case of the CS-D model, it was greater than 6% on two occasions.

3.6. Results for Sunflower in the Predictive Model

The predictive model was applied similarly to the one used for maize. As with the diagnostic model, the surprising result was that the R-P model obtained a more accurate estimate than the CS-P model based on purely sunflower pixels (Table 6).

The predictive performance exceeded expectations, as the R-P model had only one error of more than 5% (5.23% in 2019), while the remaining years were below 4%. Moreover, the estimates for 2023 and 2024 proved to be extremely accurate. Using the CS-P model, the estimation error exceeded 5% in only two cases. The MAPE value remained below 4% in both cases, being 2.73% for the R-P model and 3.58% for the CS-P model. The correlation value for the R-P model was R² = 0.928, while for the CS-P model it was R² = 0.864 (Figure 15). The results obtained by the four models are shown in Figure 16.

Similar to maize, the yield estimate simulating the operational situation was carried out, but here the R-P method was used. The results showed that the drought year of 2022 with a very low yield significantly modified the parameters of estimation, although the estimation of the yield of this year was just over 5% (Table 7).

Figure 13 indicates that the forecast horizon could not be extended substantially. In the case of the R-P model, at t_e = 220, the percentage difference between the predicted and actual returns was large in two years (15.3% in 2019 and 9.2% in 2021). However, in the other years, the error remained below 3.5%, resulting in the MAPE value of 4.08%. A similar result was obtained in the case of the CS-P model, in which the yield for the year 2019 was also underestimated by 10.3%, and the MAPE value was 4.4%. Forecast accuracy deteriorated markedly before DOY 220, with the MAPE value at DOY 210 being 9.9% and 7.2% for the R-P and CS-P models, respectively, and at DOY 200 this value increased to over 10% in both cases.

4. Discussion

This study demonstrated the successful application and refinement of two remote-sensing based yield estimation models for maize and sunflower using HLS NDVI time series in Fejér County, Hungary. The analysis covered eight consecutive years (2017–2024), allowing for the evaluation of the models under diverse weather conditions, including the extreme drought of 2022. Four modelling configurations were examined: the robust diagnostic (R-D), crop-specific diagnostic (CS-D), robust predictive (R-P) and crop-specific predictive (CS-P) models. The input data for the robust methods were the daily county-level means of NDVI, calculated over all agricultural pixels, whereas the crop-specific models used NDVI values extracted only from pure maize and sunflower fields. In the latter case, the classification had to be done accurately every year. While the diagnostic models relied on all eight years to establish the NDVI–yield relationships, the predictive models employed leave-one-year-out cross-validation to emulate operational forecasting conditions. This distinction made it possible to assess both the interpretive value and the practical forecasting skill of the modelling framework. From a practical point of view, the predictive model is more important, as this procedure is closer to the real, operational situation, where model parameters calculated based on previous years must be used to estimate the yield of the current year.

For maize, all models achieved very high accuracy, the correlation coefficient (R²) was 0.976 for the R-P model and 0.997 for the CS-P model, indicating excellent agreement in both cases but with a measurable advantage for the crop-specific approach. Using the CS-P model, the MAE value was 0.11 t/ha (RMSE = 0.12 t/ha) and the MAPE was 2.0%.

For the above results, yield estimation was performed approximately at the beginning of the crop harvest, in late September, on DOY 260. Importantly, forecast accuracy remained high even when predictions were made on DOY 230, mid-August, 30–50 days before the end of harvest, when using the CS-P model, the MAPE was 3.86%.

For sunflower, the modelling outcomes showed an unexpected pattern: the robust models outperformed the crop-specific ones. The R-P model reached an R² = 0.928, the MAE was 0.08 t/ha (RMSE = 0.09 t/ha), and the MAPE was 2.73%, while the CS-P model produced somewhat weaker results with R² = 0.866, MAE = 0.10 t/ha (RMSE = 0.13 t/ha), and MAPE = 3.58%. Notably, both models accurately captured the low yields of 2022, with the difference between the estimated and actual yield being 2.5% for the R-P model and only 0.5% for the CS-P model. Simulating a quasi-operational situation using the robust procedure, the results remained consistently reliable. The difference between the estimated and actual yield is 0.03% in 2023 and 0.16% in 2024. However, the robust models performed more consistently for sunflower, despite the crop-specific models using only sunflower pixels, and this difference warrants additional analysis. The statistical results of the four models are summarized in Table 8.

For comparison, Skakun et al. [58] estimated maize yields at sub-field scales within study sites in the US state of Iowa, and observed mixed performance of satellite-derived models with the R² varying from 0.21 to 0.88 (averaging 0.56) for the HLS. Khaki et al. [59] used data from 1132 counties for maize across the United States, and predicted the yield one to four months before the harvest, with the MAPE being 8.74%. In the work of Kumar et al. [60], regression analysis was developed between observed yield from crop cut experimental plots and seasonal maximum NDVI of maize area using Sentinel-2 data in India, and a significant relation (R² = 0.87) was obtained between NDVI and observed yield at respective field plots. Wang et al. [61] used Sentinel-2A and MODIS data to construct a yield forecasting model for summer maize in China, with the MAPE being 5.47%–13.74%. The experiment of Gavilán et al. [62] based on Sentinel-2 data was carried out in maize fields located in the Argentine Pampas Region in summer (from January to March) of 2020, obtaining a high (R² = 0.91) correlation value. Birinyi et al. [63] used a six-category crop condition mapping method based on Sentinel-2 data, which explained 77% of county-level yield variability for maize in July. Adriko et al. [64] compared results based on MODIS and Sentinel-2 data in the Midwestern corn belt in the United States; the correlation value was R² = 0.79 for Sentinel-2 and R² = 0.87 for MODIS. Radocaj et al. [65] analyzed Sentinel-2 vegetation index time-series and phenology metrics for Iowa and Illinois, and reported that the combination of wide dynamic range vegetation index and the maximum-value metric achieved the strongest relationship with maize yield, with a mean correlation of about 0.51. In the case of sunflower, Narin et al. [66] obtained R² = 0.74 approximately three months before the harvest in Turkey using Sentinel-2 data. Amankukova et al. [67] applied multiple linear regression and two different machine learning approaches to predicting the crop yield using Sentinel-2 data. They found that the most suitable time for predicting the crop yield was 86–116 days after sowing, the best machine learning approach for predicting field-scale variability of the crop yield was about R² = 0.6 and RMSE = 0.284–0.473 t/ha.

Despite methodological differences, the county-level performance of our models achieved similar or better explanatory power through a simpler, continuous NDVI-based regression model, thereby reducing model complexity and computational cost.

5. Conclusions

In our current research, we had four main objectives: (i) to improve the previously developed robust method based on 30-m data, (ii) to test the model on input data from areas with different vegetation cover, (iii) to estimate the predictive power of the model, and (iv) to test the accuracy of the model for forecasts at different dates. Our results demonstrate that these objectives were successfully achieved through the combined use of HLS data, and regression-based yield models. Both the robust and the crop-specific approaches provided stable and reliable performance, although with different advantages. In the case of maize, the crop-specific model provided a more accurate estimate; however, due to its simplicity and low data requirements, the robust model can be more suitable where a detailed annual crop map is limited or unavailable. Importantly, the models responded well even under extreme growing season conditions, such as the drought year of 2022, highlighting their resilience to increasing climate variability.

Testing the predictive model, it was found that maize yields can be forecast with high accuracy as early as 30–50 days before the end of harvest, which represents significant value for operational decision-making. For sunflower, the early-season forecasting proved more limited; nevertheless, both models maintained reliable accuracy on DOY 230, around the typical beginning of the harvest. Overall, the presented approach offers a favourable balance between model simplicity, data requirements and predictive ability, outperforming or matching the results of more complex methods reported in the literature while requiring fewer input variables. The results also confirm that dense NDVI time series derived from HLS data are suitable for operational regional-scale yield forecasting. Consequently, the methodology can clearly be integrated into regional or national early warning systems. Future research may extend this work to additional crop types and incorporate meteorological or soil moisture information to increase the accuracy and generalizability of the model. It would also be important to include additional years in the studies so that longer time series enhance the stability of the model.