Within-Field Relationships between Satellite-Derived Vegetation Indices, Grain Yield and Spike Number of Winter Wheat and Triticale

The aims of this study were to: (i) evaluate the relationships between vegetation indices (VIs) derived from Sentinel-2 imagery and grain yield (GY) and the number of spikes per square meter (SN) of winter wheat and triticale; (ii) determine the dates and plant growth stages when the above relationships were the strongest at individual field scale, thus allowing for accurate yield prediction. Observations of GY and SN were performed at harvest on six fields (three locations in two seasons: 2017 and 2018) in three regions of Poland, i.e., northeastern (A—Brożówka), central (B—Zdziechów) and southeastern Poland (C—Kryłów). Vegetation indices (Normalized Difference Vegetation Index (NDVI), Soil-Adjusted Vegetation Index (SAVI), modified SAVI (mSAVI), modified SAVI 2 (mSAVI2), Infrared Percentage Vegetation Index (IPVI), Global Environmental Monitoring Index (GEMI), and Ratio Vegetation Index (RVI)) calculated for sampling points from mid-March until mid-July, covering within-field soil and topographical variability, were included in the analysis. Depending on the location, the highest correlation coefficients (of about 0.6–0.9) for most of VIs with GY and SN were obtained about 4–6 weeks before harvest (from the beginning of shooting to milk maturity). Therefore, satellite-derived VIs are useful for the prediction of within-field cereal GY as well as SN variability. Information on GY, predicted together with the results for soil nutrient availability, is the basis for the formulation of variable fertilize rates in precision agriculture. All examined VIs were similarly correlated with GY and SN via the commonly used NDVI. The increase in NDVI by 0.1 unit was related to an average increase in GY by about 2 t ha−1.


Introduction
Satellite remote sensing (RS) helps in the mapping of current crop status and the assessment of biophysical parameters. Currently, RS data are publicly accessible due to the availability of an unprecedented amount of Free Sentinel data from the Copernicus Program, established by the European Space Agency [1]. Constellations of the Sentinel-2 satellites (2A and 2B) can be also used for precision farming applications [2]. Thanks to the high spatial resolution (pixel size of 10 m), relatively short revisit time (about 5 days) and multispectral sensors, it is possible to observe and analyze the crop    According to the Köppen classification [19], the majority of Poland lies in the warm temperate climate zone with fully humid continental climate and hot summers (Cfb). The locations of Kryłów (C) and Zdziechów (B) are situated within this zone, while the location of Brożówka (A) is placed in the transitional area between the Cfb zone and the snowy climate zone, with fully humid and warm summers (Dfb).
The sum of precipitation in two studied seasons (from the beginning of September to the end of August) ranged from about 520 mm (location C, season 2016/2017) to 740 mm (location A, season 2016/2017). The lowest air temperatures were observed for location A, while locations B Agronomy 2020, 10, 1842 4 of 18 and C were characterized by similar air temperatures. Growing degree days calculated at a base air temperature of 5 • C for the period from the sowing date to the end of July were much lower in 2017 (for locations: A-1364, B-1531, and C-1492) in comparison to 2018 (for locations: A-1795, B-1975, and C-1930). Monthly sums of precipitation and mean air temperatures for the nearest meteorological stations from the Institute of Meteorology and Water Management National Research Institute IMGW-PIB (https://danepubliczne.imgw.pl/) are presented in the Figure 2. (C) and Zdziechów (B) are situated within this zone, while the location of Brożówka (A) is placed in the transitional area between the Cfb zone and the snowy climate zone, with fully humid and warm summers (Dfb).
The sum of precipitation in two studied seasons (from the beginning of September to the end of August) ranged from about 520 mm (location C, season 2016/2017) to 740 mm (location A, season 2016/2017). The lowest air temperatures were observed for location A, while locations B and C were characterized by similar air temperatures. Growing degree days calculated at a base air temperature of 5 °C for the period from the sowing date to the end of July were much lower in 2017 (for locations: A-1364, B-1531, and C-1492) in comparison to 2018 (for locations: A-1795, B-1975, and C-1930). Monthly sums of precipitation and mean air temperatures for the nearest meteorological stations from the Institute of Meteorology and Water Management National Research Institute IMGW-PIB (https://danepubliczne.imgw.pl/) are presented in the Figure 2. Monthly sums of precipitation, mean air temperature and growing degree days (GDD) for meteorological stations located closest to the research sites (for location A-Olecko, for location B-Puławy, for location C-Strzyżów) from the IMGW-PIB national institute (source of data: https://danepubliczne.imgw.pl/).

Soil and Plant Sampling
Within each of the fields, sampling points of soil and plants were georeferenced using a Global Positioning System (GPS) receiver. The locations of the sampling points were selected to obtain large variability in the studied traits, i.e., soil properties as well grain yield and ear number ( Figure 3). The characteristics of each of the three research sites are presented in Table 1.
Soil samples were collected from a layer of 0-90 cm at the end of winter or in early spring before the starting of the vegetation. Grain yield and number of ears per square meter were determined at the full maturity of plants in plots of 2 m 2 (2 × 1 m 2 subplots, each located very closely to the soil and plant sampling points). The following statistical parameters were calculated: average, standard deviation, minimum, and maximum for each field in a given year.

Soil and Plant Sampling
Within each of the fields, sampling points of soil and plants were georeferenced using a Global Positioning System (GPS) receiver. The locations of the sampling points were selected to obtain large variability in the studied traits, i.e., soil properties as well grain yield and ear number ( Figure 3). The characteristics of each of the three research sites are presented in Table 1.
Soil samples were collected from a layer of 0-90 cm at the end of winter or in early spring before the starting of the vegetation. Grain yield and number of ears per square meter were determined at the full maturity of plants in plots of 2 m 2 (2 × 1 m 2 subplots, each located very closely to the soil and plant sampling points). The following statistical parameters were calculated: average, standard deviation, minimum, and maximum for each field in a given year.

Satellite Data
Satellite Sentinel-2 (2A and 2B) images of Level-1C with a spatial resolution of 10 to 20 m in the range of visible and near-infrared (NIR-band 8 at central wavelength 833 nm) light were downloaded from the Copernicus Open Access Hub [20]. Level-1C processing includes radiometric and geometric corrections, namely ortho-rectification and spatial registration, on a global reference system with a sub-pixel accuracy. Then, the acquired Level 1C products were subjected to atmospheric correction using the Dark Object Subtraction (DOS1) method in Semi-Automatic Classification Plugin documentation (SCP), developed by Luca Congedo (2016) [21].

Satellite Data
Satellite Sentinel-2 (2A and 2B) images of Level-1C with a spatial resolution of 10 to 20 m in the range of visible and near-infrared (NIR-band 8 at central wavelength 833 nm) light were downloaded from the Copernicus Open Access Hub [20]. Level-1C processing includes radiometric and geometric corrections, namely ortho-rectification and spatial registration, on a global reference system with a sub-pixel accuracy. Then, the acquired Level 1C products were subjected to atmospheric correction using the Dark Object Subtraction (DOS1) method in Semi-Automatic Classification Plugin documentation (SCP), developed by Luca Congedo (2016) [21].
Vegetation indices were derived for the representative soil and plant sampling points established within fields and imported together with the coordinates of these points. Then, they were compiled Agronomy 2020, 10, 1842 6 of 18 into tables and subjected to regression analyses. Depending on the location, despite the revisit time of about 5 days, it was only possible to acquire cloudless satellite images during the intensive growth of cereals in spring and early summer for 4-6 dates ( Table 2).

Spectral Vegetation Indices
The selection of the vegetation indices (VIs) used for the winter wheat and triticale yield analysis depended on several factors. First of all, we only considered the same satellite data spatial resolution. Another factor determining our choices was the prevalence in the literature and the comparability with the yield of cereal crops. The VIs were calculated on the basis of the reflectance of wavelength bands registered by Sentinel-2 at spatial resolution 10 m (squared pixels 10 × 10 m). The Normalized Differential Vegetation Index (NDVI) allows us to determine the indirect absorption of photosynthetic radiation on a landscape scale. NDVI calculation is based on the contrast between the largest reflection in the near-infrared band and the absorption in the red band. In the case of this index, the difference in reflection in the near-infrared and red bands is divided by their sum. This approach compensates for the differences in the radiation amounts in both bands. It is believed that NDVI is more sensitive to small amounts of vegetation [22,23].
Another VI used in this study was the Soil-Adjusted Vegetation Index (SAVI). This indicator was designed to minimize the effect of soil reflection on red and near-infrared radiation by adding an estimated background correction factor [24].
where L is a canopy background adjustment factor. An L value of 0.5 was adopted to minimize soil brightness variations and eliminate the need for additional calibration for different soils.
In the literature, modifications of the SAVI index, known as modified SAVI (mSAVI) and modified SAVI 2 (mSAVI2), are often used. These two VIs also include the soil background correction factor [25]: where L is: and: Agronomy 2020, 10, 1842 7 of 18 Crippen (1990) proposed the Infrared Percentage Vegetation Index (IPVI), which simplifies the calculation by eliminating the subtraction of the red radiation value in the NDVI indicator counter. This simplification of a new VI calculation proved to be important for fast image processing [26].
In 1992, Pinty and Verstraete proposed the Global Environmental Monitoring Index (GEMI), a new nonlinear crop status indicator, which includes a correction factor adjusted by the effect of soil reflectance and atmospheric background [27].
The last VI used in this study for a comparison with NVDI was the Ratio Vegetation Index (RVI), sometimes also referred to as the Simple Ratio (SR). This VI is a ratio of NIR and red reflectance [28].
The investigated Vis, except RVI, are normalized and placed on a comparative scale. These indices have values from −1 to 1, where values close to −1 indicate water, or any inanimate matter, values from 0 to 0.20 relate to barren rock, sand, barren soil and plants at an early stage of growth, values of 0.20-0.50 relate to sparse vegetation such as shrubs and grasslands, while values above 0.60 indicate crops at their peak growth stages.

Statistical Analysis
The Orfeo Toolbox (OTB Development Team, 2018) [29] was used for the calculation of vegetation indices for pixels located nearest to the soil and plant sampling points in 2017 and 2018. The indices was subjected to statistical analysis and the mean, standard deviation, coefficient of variation (CV), minimum, maximum were determined (Tables 3 and 4). Relationships between VI values and grain yield as well as the number of spikes per square meter were evaluated using Pearson's correlation and a linear regression analysis. These analyses were conducted for the dates from mid-March to mid-July ( Table 2). The choice of dates remained strictly dependent on the availability of cloudless Sentinel-2 satellite images. Depending on location, there were 4 to 6 dates (obtained, on average, every 3-4 weeks) in which the images were available without clouds.   The correlation coefficient values for each date are presented over time in a graphical form for better visualization of the relationship change between NDVI and grain yield (Tables 5 and 6). For periods when maximal values of the correlation coefficients were observed, the diagrams of linear regression are presented ( Figures 5 and 6), together with regression equations and coefficients of determination (R 2 ).
Spatial analyses were conducted, and maps were prepared using QGIS 2.18 software (QGIS Development Team, Gossau, Switzerland) [30], while the statistical analyses were performed using GNU R and TIBCO Statistica software.

Grain Yield and Spike Number
The highest grain yield as well the highest number of spikes per m 2 was observed for location C (Kryłów) in both years (Table 3). This was mainly due to the very favorable soil conditions where the research sites were located in both years in (Table 1). In 2017, soil conditions were very good and uniform within the very flat field (denivelation of less than 2 m). This caused very low yield variability (SD = 0.63) within this field. In 2018, soil conditions were slightly less favorable and more variable, but still quite uniform in comparison to soil conditions in the other two locations (A and B). In location A (Brożówka), the average grain yield of winter wheat was much higher (8.33 t·ha −1 ) in 2017 than the yield of winter triticale in 2018 (5.10 t·ha −1 ). The main reason for the much lower grain yield as well as the lower spike number per m 2 (458 spikes in 2017 versus 317 in 2018) was the shorter tillering time in 2017/2018, related to the later sowing date of winter triticale when compared to winter wheat (2016/2017). Moreover, we must mention that winter wheat is usually sown later than winter triticale. In both years in location A, yield variability within this field was quite high, mainly due to the undulated surface of the field (denivelation of about 20 m), which determined the soil variability. In location B, in both years, a grain yield of about 5 t ha −1 and a spike density of about 400 spikes per m 2 were relatively low and very variable within the field. The main reason for such high within-field variability of the grain yield traits was the variable soil texture (Table 1), which caused varied water availability for plants. A shortage of water usually occurred in later growth stages, especially during heading and grain filling. This caused very low-weight grains and, despite quite large numbers of spikes, grain yield was very low. Such a situation was especially visible in the sandy parts of the fields in location B.

Changes in NDVI over the Vegetation Season and in Research Locations
Due to cloud cover, the availability of useful satellite images from Sentinel-2 was limited to 4-6 scenes per season during intensive growth of winter cereals, i.e., from half of March to half of July. Despite this limitation, it was possible to evaluate changes in VI values during the growing seasons. The most substantial differences between the average and a range of NDVI values for the three locations were observed from the end of March to the beginning of April. For location A, mean values of NDVI in both years were very low (0.3-0.4) until the beginning of April. The main cause of these low NDVI values in location A were lower air temperatures and consequently growing degree days (from sowing date to the end of April in 2017: location A-283, B-298, and C-451; in 2018: A-415; B-501, and C-605) during early spring than in the other two locations (Figure 2). The highest values of NDVI (of about 0.7 at the beginning of April) were observed in location C. This was caused by warmer conditions in both autumn and early spring as well as by more favorable soil conditions, which caused more intensive plant growth. The very high values of NDVI (0.7 or higher) in location C were observed for a much longer time, even until the end of June, in comparison to the other two locations. Consequently, this longer and more intensive crop growth in location C allowed us to obtain higher grain yields. In location B, NDVI values in early spring were at a medium level (higher than in location A and lower in comparison to location C) and very variable depending on the season at the later growth stages. In all locations, in 2017, a maximum NDVI value was achieved in very late June and much earlier in 2018, at the end of April (Figure 4, Table 4).

Relationships between NDVI and Grain Yield and Spike Number
One of the aims of the study was to determine the dates and plant growth stages when the relationships between VI values and grain yield, as well as number of spikes, were the strongest. The strongest correlations were achieved at various growth stages (Tables 5 and 6).
For fields located in northeastern Poland, the highest correlation coefficients between grain yield and as well as number of spikes and NDVI were observed at the end of May (approximately shortly before or during the heading stage) in both years. On the basis of the regression slope values, it can be concluded that an increase in NDVI by 0.1 unit corresponded to the increase in grain yield by about 1.5-1.6 t•ha −1 in both years (Figures 5 and 6).

Relationships between NDVI and Grain Yield and Spike Number
One of the aims of the study was to determine the dates and plant growth stages when the relationships between VI values and grain yield, as well as number of spikes, were the strongest. The strongest correlations were achieved at various growth stages (Tables 5 and 6).
For fields located in northeastern Poland, the highest correlation coefficients between grain yield and as well as number of spikes and NDVI were observed at the end of May (approximately shortly before or during the heading stage) in both years. On the basis of the regression slope values, it can be concluded that an increase in NDVI by 0.1 unit corresponded to the increase in grain yield by about 1.5-1.6 t·ha −1 in both years (Figures 5 and 6). In location B, the strongest relationship between NDVI and grain yield was registered in late June of 2017 (end of milk maturity of winter wheat) and in the second half of May 2018 (shooting and heading of winter triticale). In this research site, the relationships were the strongest among the three examined locations, i.e., the R 2 value was about 0.70-0.75 (Figures 6 and 7). Moreover, the relationship between NDVI and the number of spikes was slightly weaker. The increase in NDVI by 0.1 unit was related to an increase in grain yield of about 8.0 t•ha −1 in 2017 and about 2.5 t•ha −1 in 2018.
In location C in 2017, the strongest correlation between NDVI and both grain yield and number of spikes was observed, respectively, on 21 June and 11 July (milk and dough maturity). In 2018, the strongest relationships between NDVI and both grain yield and spike number per square meter were obtained at the end of April and June. Due to a lack of data (clouds), it was not possible to verify this relationship at the end of May. The relationship was much stronger in 2018 than in 2017, mainly due to the considerably higher within-field variability of soil and consequently grain yield variability (Coefficient of Variation-CV) for grain yield was 6% and 16%, respectively, for 2017 and 2018). Based on the regression equations ( Figures 5 and 6), it was evaluated that an increase in NDVI by 0.1 unit was related to the increase in grain yield by about 1.0 t•ha −1 in 2017 and about 2.1 t•ha −1 in 2018. For all six site years, the strongest correlation between NDVI and grain yield and number of spikes was  observed on similar dates and at similar growth stages. This is because the number of spikes per square meter usually very strongly affects grain yield.  In location B, the strongest relationship between NDVI and grain yield was registered in late June of 2017 (end of milk maturity of winter wheat) and in the second half of May 2018 (shooting and heading of winter triticale). In this research site, the relationships were the strongest among the three examined locations, i.e., the R 2 value was about 0.70-0.75 (Figures 6 and 7). Moreover, the relationship between NDVI and the number of spikes was slightly weaker. The increase in NDVI by 0.   In location C in 2017, the strongest correlation between NDVI and both grain yield and number of spikes was observed, respectively, on 21 June and 11 July (milk and dough maturity). In 2018, the strongest relationships between NDVI and both grain yield and spike number per square meter were obtained at the end of April and June. Due to a lack of data (clouds), it was not possible to verify this relationship at the end of May. The relationship was much stronger in 2018 than in 2017, mainly due to the considerably higher within-field variability of soil and consequently grain yield variability (Coefficient of Variation-CV) for grain yield was 6% and 16%, respectively, for 2017 and 2018). Based on the regression equations ( Figures 5 and 6), it was evaluated that an increase in NDVI by 0.1 unit was related to the increase in grain yield by about 1.0 t·ha −1 in 2017 and about 2.1 t·ha −1 in 2018. For all six site years, the strongest correlation between NDVI and grain yield and number of spikes was observed on similar dates and at similar growth stages. This is because the number of spikes per square meter usually very strongly affects grain yield.

Relationships between the Other Vegetation Indices and Grain Yield and Spike Number
Regression analyses performed to evaluate the strength of the relationship between the vegetation indices (SAVI, mSAVI, mSAVI2, GEMI, IPVI, RVI) and grain for each crop in 2017 and 2018 (Figure 7) indicated that the strongest relationships were observed for most of the VIs on the same dates. In location A, the highest correlations were observed at the end of May (end of shooting stage) in 2017, and from late May to mid-June (heading stage) in 2018. In location B, in 2017, the value of the correlation coefficients differed for individual VIs on several dates, while, in 2018, the strongest relationships were achieved in mid-May (beginning of heading stage) and at the end of June (milk maturity stage) for all VIs. The highest differentiation in the correlation coefficient values for the relationship between various VIs and yields was observed at the end of June and beginning of July (milk maturity stage) of 2017 in location C. In 2018, in all locations, the correlation coefficients had similar values for all VIs across the whole season.
The vegetation indices (VI) are ranked by order of decreasing correlation (data not shown) with yield as follows: SAVI, MSAVI, NDVI, MSAVI2, IPVI, RVI, GEMI. This ranking was carried out by averaging the R 2 values from all of the measurement dates and research areas. The strength of the correlation of VIs with the number of ears was evaluated in the same way as above, and the order of decreasing correlation was as follows: NDVI, SAVI, MSAVI, MSAVI2, IPVI, RVI, GEMI.

Discussion
Constant field monitoring using RS methods is important not only for grain yield prediction, but also for the assessment of the site-specific conditions of plant development during the growing season [31,32]. This may allow us to introduce site-specific management of the field, and thus increase yields or save inputs.

Relationships between NDVI and Grain Yield
Our results for winter cereals, which constitute the majority of crop production in Central Europe, proved moderate and strong relationships (R 2 in range of 0.35-0.81 for all fields in two seasons, Figures 5 and 6) between the NDVI and grain yield of winter wheat and triticale for all of the studied fields. At location A, the correlation coefficient grew from mid-March to the end of May, when it reached the highest value. From the beginning of June to mid-July, we noticed a decrease in this correlation. The highest correlation coefficients are observed from the beginning of June to mid-July due to the growth stages such as heading, milk maturity, or dough maturity for location B. At the site C, we can see high rates of milk maturity at the end of June. A similar tendency, though with lower correlation values, can be observed when comparing the number of spikes per square meter. Similar results were achieved by Benedetti and Rossini (1993) [4] in Italy and by Smith et al. (1995) [5] in Western Australia for wheat. In both studies, AVHRR/NOAA satellite data were used for the calculation of NDVI. In Italy, R 2 coefficient values for the relationship between NDVI and grain yield ranged from 0.24 to 0.749, and for Australia they ranged from 0.46 to 0.72. Ali et al. (2019) [17] observed higher R 2 values of 0.878 and 0.926 for the relationship between the grain yield of winter and bread durum wheat and NDVI when the values of this index were derived from Landsat images. In turn, Labus et al. (2002) [6] achieved low R 2 values ranging from 0.00 to 0.69 when wheat yield prediction was based on NDVI derived from AVHRR/NOAA satellite data. In our study, very low R 2 values were obtained in April and May, and higher R 2 values were mainly close to the heading stage and at the end of the growing season (Tables 5 and 6). Therefore, while comparing many different studies, we can expect that, regardless of the research region, the relationship between NDVI and the grain yield remains at a similar level depending on the growth stage of the studied cereal.

Determination of Dates and Plant Growth Stages When Relationship between NDVI and Grain Yield Was the Strongest
Veloso et al. (2017) [10] found, for southeastern France, that NDVI reached maximum values in the second half of May, while, in Poland, depending on the research location, the highest NDVI values were reached from the end of April to the end of June (Figure 4). The site-related differences in the absolute NDVI values and the date at which they reached their maximum may primarily be influenced by the environmental conditions, mainly weather and soil conditions. This, in turn, closely depends on the climate in a specific region (amount of rainfall, growing degree days) and the amount of nitrogen applied. In the study of Naser et al. (2020) [33], much stronger correlations between NDVI, measured using a ground sensor, and the grain yield of winter wheat were observed for dryland than for irrigated conditions. The reason for such results was the greater range of absolute NDVI values for dryland in comparison to irrigated conditions. In that study, the highest correlations were observed after anthesis for both types of conditions. Very strong relationships between NDVI, measured using a ground spectrometer, and wheat yield were observed in a study where the effect of various nitrogen doses (from 0 to 210 kg/ha) was investigated [34]. The strongest correlations (R 2 up to 0.96) were achieved at the heading stage. According to Satir and Berberoglu (2016) [35], the strength of the relationship between satellite-derived NDVI and the yields of wheat, corn and cotton was modified by soil conditions, e.g., soil salinity. This is because soil surface and soil water content can affect the values of vegetation indices, as well as the fact that water availability for plants is limited in high-salinity conditions. A better prediction of grain yield in the variable soil conditions of south Turkey can be achieved not only using prediction models such as NDVI, but also other spectral indices such as Normalized Difference Water Index (NDWI), Soil-Adjusted Vegetation Index (SAVI) or Tasseled Cap Wetness Index (WETNESS). The highest prediction accuracy for wheat was using the model which includes NDVI, NDWI and WETNESS. Prediction of wheat yield based on satellite-derived NDVI was improved if auxiliary data such as grain yield from previous seasons were included in models [36]. Dempewolf et al. (2013) [12] obtained the highest R 2 values of 0.964 for NDVI and the yield of wheat grown in Pakistan six weeks before harvest. In our studies, the highest R 2 value was observed at the end of April, i.e., about 12 weeks before harvest in the case of the field located in southeastern Poland. These results are similar to those of Kussul et al. [13] (2014), obtained in Ukraine, where a good prediction of winter wheat yields was possible even 2-3 months before harvest. The latest date (beginning of July) when the correlation coefficient for the relationship yield versus NDVI was the highest was observed in location B (Zdziechów) in central Poland in 2017. Similar results were obtained in Hungary by Nagy et al. (2018) [15] and Lopresti et al. (2015) [11] in Argentina, where good yield prediction was possible six weeks before harvest. Ali et al. (2019) [17] observed the highest correlation between NDVI and yield for winter bread wheat from stem elongation to the milk maturity of the grain, i.e., 6-17 weeks before harvest. In summary, the choice of the date of the yield estimation with the best accuracy does not depend only on long-term field and satellite measurements. Environmental, weather, and soil conditions should also be taken into account. As shown above, a yield estimation can be properly performed even between 3 months and 6 weeks before harvest, depending on the research region and the field-specific conditions.

Relationship between Other VIs and Grain Yield
According to Purevdorj et al. (1998) [37], SAVI and its modifications usually show higher correlations with grain at the beginning of the growing season compared to NDVI. According to Ren et al. (2018) [38], SAVI, with its modifications, should be used to estimate vegetation coverage at low vegetation density on arable land. In our study, the value of the correlation coefficient for the relationship between SAVI and grain yield was similar to the relationship with NDVI for both seasons and locations. Moreover, only in locations A and C in 2017 were the correlations between SAVI and grain yield stronger in early spring. For location A, during two vegetative seasons, we saw an increase in the correlation coefficient at the beginning of the growing season and a decrease in the cereals' subsequent growth stages. For location B, the correlation coefficient increases until the beginning of June and slightly decreases until the end of July. In tested area C, an increase in the correlation coefficient value is noticed until mid-June, followed by a sharp decline by the end of July.
Ali et al. (2019) [17] selected NDVI and SR (referred to RVI in our study) as better VIs in comparison to EVI, SAVI, GNDVI, and GCI. Their conclusion was that SAVI never exhibited the strongest correlation with yield in comparison with other VIs, although, its correlation with grain yields was, most frequently, significant and strong. On the contrary, our research results proved that SAVI (and its modifications) and NDVI are ranked as VIs with the highest correlations with yield and number of spikes.
Moreover, Ali et al. (2019) [17] obtained R 2 values of about 0.8 for durum and bread wheat yield and NDVI, EVI, and SR. We achieved Pearson's correlation values of 0.4 between winter wheat and triticale yield and SAVI, but 0.3 with mSAVI, NDVI, IPVI, RVI (SR), and GEMI. These differences most probably arise from a different method being used for the collection of yield data.

Conclusions
Various vegetation indices (SAVI, MSAVI, NDVI, MSAVI2, IPVI, RVI and GEMI), their calculation based on the red and near-infrared radiation derived from the Sentinel−2 imagery, showed similar relationships with grain yield and number of spikes per square meter of winter wheat and triticale as the commonly used NDVI. While comparing the grain yield and the number of spikes from NDVI for all three locations during two seasons, it can be seen that the values of R 2 for the number of spikes per square meter are lower than the values of R 2 for the grain yield. Consequently, the grain yield is a better yield prediction parameter than number of spikes per square meter. In general, the relationship between NDVI values and grain yield was stronger at more advanced growth stages. Depending on the region of Poland, the strongest correlations between NDVI and yield and its main component were obtained for NDVI derived from images obtained at the end of April (beginning of shooting) in southeastern Poland, at the end of May (beginning of heading) in the northeastern part of the country and only at the end of June (at milk maturity) in the central region of Poland. The divergence in these periods may result from the sowing date, but also from different climatic conditions, especially GDD values, in the three micro-climate regions of Poland. The different strengths of the relationships may be caused by the soil water deficiency due to variable weather and soil conditions, depending on the research location, at different growth stages. Within-field water variability was higher, especially on sandy soils (which were common in some parts of the field in location B) at later growth stages when air temperatures were higher and negative balance between rainfall and evapotranspiration was more common than in locations A and C [39]. Moreover, the presence of clouds, which made the registration of useful satellite images impossible on some dates, also contributed to the results obtained in this study. For example, in location C in southeastern Poland, there were no cloudless images available in May in both seasons and, because of this, it was not possible to evaluate the relationships between VIs and yield for some important growth stages of the crop. The proposed estimated time for accurate yield prediction is about 4-6 weeks before harvest (from the beginning of shooting to milk maturity).