1. Introduction

Timely and accurate information on crop yields and production at global, national, and regional scales is extremely important for many agriculture applications ^[1]. At national/regional scale, it can be an input to local authorities to make decisions on food security issues or deciding on subsidies in case of extreme weather conditions such as droughts. At field scale, spatial variability of yields can help to obtain objective information, for example, for farmers to improve management practices and identify yield gaps ^[2], or for insurance companies to feed this information into insurance models ^[3,4].

Owing to its coverage, temporal and spatial resolution, remote sensing images from space has always been a powerful tool to develop empirical models for predicting and assessing yields at regional and national scales ^{[5,6,7,8,9,10,11]}, or assimilating biophysical parameters into crop growth models ^[12,13,14]. In particular, coarse resolution sensors, e.g. Moderate Resolution Imaging Spectroradiometer (MODIS), Advanced Very High Resolution Radiometer (AVHRR), SPOT-VEGETATION, thanks to its daily coverage and availability of historical datasets dating back to 1980s–1990s, have extensively been used for building empirical models for crop yield forecasting and assessment. These models connect satellite-derived features, for example vegetation indices (VIs) such as Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Vegetation Health Index (VHI) and/or biophysical parameters such as Leaf Area Index (LAI), Fraction of Photosynthetically Active Radiation (FPAR), with reference yield data. For example, Johnson (2016) ^[5] analyzed efficiency of multiple MODIS land products including NDVI, EVI, LAI, FPAR, and Gross Primary Production (GPP) to assess crop yield at county level in US for ten major agriculture commodities. He found positive correlations of vegetation products against yield for all crops, except rice, and that finer spatial resolution improved the correlations. López-Lozano et al. (2015) ^[6] investigated the use of the Fraction of Absorbed Photosynthetically Active Radiation (fAPAR) derived from SPOT-VEGETATION at 1 km spatial resolution to assess crop yields (wheat, barley and maize) at province level in Europe. They found high correlations (R² > 0.6) in water-stressed regions; however, lower correlations (R² < 0.5) were observed for regions with high yields. Salazar et al. (2007) applied AVHRR-derived VHI to estimate winter wheat yield in Kansas, US, and found high correlations with official statistics for 1982–2004. NDVI, as well as biophysical parameters LAI and fAPAR, also proved to be efficient in predicting winter wheat yields at different scales in Ukraine ^[15]. In order to overcome some limitations of empirical models in terms of robustness, Becker-Reshef et al. ^[10] developed a generalized winter wheat yield forecasting model that was calibrated for one region (Kansas, US) and successfully applied for another one (Ukraine) to provide an error of less than 10% that can be suitable for operational context. Adding meteorological data, in particular temperature, has usually had a positive effect on crop yield models reducing the error and improving timeliness ^[5,6,7]. Though these models are empirical and based on relative simple equations, they perform at the same level, or even better, than more comprehensive crop models that are based on crop growth simulations ^[8,16]. The reasons for that are: complexity of accounting multiple factors influencing the yield, lack of high-quality data required to calibrate and run such models, and difficulties of upscaling "point" estimates to higher spatial scale ^[17].

The use of Landsat and Landsat-like 30 m data for crop yield forecasting and assessment has been limited due to its infrequent revisit rate of only 16 days. There have been works fusing Landsat with higher temporal frequency, but spatially coarser, MODIS sensor ^[18,19], or combining Landsat with biophysical models ^[20,21]. For example, Lobell et al. (2015) proposed a scalable satellite-based crop yield mapper (SCYM) that is based on training a statistical model from crop simulations and applying the statistical model to Landsat-5/7 images and meteorological data ^[21]. They achieved on average R² values of 0.35 and 0.32 for corn and soybeans, respectively, for the large areas in the Midwestern United States. However, these approaches showed varying results in terms of errors and still have limitations constrained by lower frequency of moderate resolution images. With the combined use of Landsat-8 and Sentinel-2 remote sensing satellites that will enable acquisition of an image every 3–5 days globally, it becomes possible to implement approaches similar to those used for MODIS/AVHRR to develop next generation agriculture products at higher spatial resolution (30 m).

This paper presents one of the first studies to combine Landsat-8 and Sentinel-2A imagery for crop yield mapping by downscaling a generalized empirical model developed for MODIS data ^[7,10]. The model is based on capturing the peak NDVI to correlate with the yield, and growing degree days (GDD) to improve the timeliness of the model. Therefore, the main objectives of the study are: (ⅰ) to assess performance of downscaling the generalized NDVI–based empirical model for winter wheat yield forecasting from coarse spatial resolution to moderate one at 30 m; (ⅱ) to explore the combined use of images acquired by Landsat-8 and Sentinel-2A remote sensing satellites for winter crop mapping and winter wheat yield assessment at regional level.

---

2. Study Area & Materials

---

2.1. Study Area and Reference Data

The study was performed for Kirohohradska oblast in Ukraine for 2016 (Figure 1). Ukraine is a top 10 wheat producer in the world. An oblast is a high-level administrative division of the country (there are 24 oblasts in Ukraine and Autonomous Republic of Crimea), and each oblast is further divided into districts. Kirovhradska oblast is located in the central part of Ukraine and composed of 21 districts with geographical area ranging from 65 to 165 thousand ha and cropland area ranging from 27 to 112 thousand ha. The reason for selecting this region is that it is a top 10 wheat producer in Ukraine and because of availability of reference crop yield and harvested area data at district scale for 2016. Winter wheat is one of the major crops in Kirovhradska oblast accounting for 20% of production of all crops in the region. Winter wheat is mainly rain-fed in the region and usually planted in September-October. After dormancy during the winter, it re-emerges early spring reaching maturity by the end of June. Harvest of winter crops is typically undertaken in July.

Reference data on crop yield and harvested area at district level were collected from the Department of Agro-Industry Development of Kirovohrad State Administration (http://apk.kr-admin.gov.ua). The data were made available online as the harvest progressed and were based on farm surveys of all agricultural enterprises (that account of more than 90% of all winter crops production in the region) and samples of household farms the same way as official statistics is collected ^[22]. The number of samples for surveying small household farms is selected in such a way to target a coefficient of variation of 10%. The final estimates for winter crop yields and areas were available at the end of November and were used as reference in this study. Uncertainty of reference data should not exceed 10% ^[23].

---

2.2. Landsat-8/OLI and Sentinel-2A /MSI Datasets

Remote sensing images acquired by the Operational Land Imager (OLI) instrument aboard Landsat-8 satellite and by the Multi-Spectral Instrument (MSI) aboard Sentinel-2A satellite were used in the study. Landsat-8/OLI captures images of the Earth's surface in 9 spectral bands at 30 m spatial resolution (15 m for panchromatic band) ^[24] while Sentinel-2A/MSI captures images of the Earth's surface in 13 spectral bands at 10 m, 20 m and 60 m spatial resolution ^[25]. The main bands that were used in the study are bands 4 (Red) and 5 (NIR) from Landsat-8, and bands 4 (Red) and 8A (NIR) from Sentinel-2A. Band 8A from Sentinel-2A was selected instead of band 8 since spectral response function for band 8A is similar to the Landsat-8's band 5 (Figure 2).

Overall, 51 Landsat-8 and 87 Sentinel-2A scenes were acquired over the study area from March 1, 2016 to July, 31, 2016. Landsat-8 images were downloaded from the USGS's Earth Explorer (Pre-Collection Level 1) and Sentinel-2A images were downloaded from the Copernicus Open Access Hub (SciHub) with baseline processing version ranging from 02.01 to 02.04. Landsat-8 scenes covered the following coordinates (path/row) of the World-wide Reference System (WRS-2): 178/026, 179/026, 179/027, 180/026, 180/027, and 181/026. The tile size of Landsat-8 is approximately 185 km × 180 km. Sentinel-2A scenes covered the following tiles: 35UQQ, 35UQP, 36UUV, 36UUU, 36UVV, 36UVU, 36UWV, and 36UWU. The size of the Sentinel-2A tile is approximately 110 km × 110 km (Figure 3). Dates of Landsat-8 and Sentinel-2A acquisitions are given in Appendix A.

The Landsat-8/OLI and Sentinel-2A/MSI scenes were atmospherically corrected for surface reflectance using the LaSRC algorithm ^[26] (Figure 4 and Figure 5) ensuring consistency between these datasets as well as with MODIS data used for building a generalized crop yield model ^[10,28]. Figure 4 shows an example of true and false color compositions of Landsat-8 and Sentinel-2A acquired on the same date. A quantitative analysis with performance metrics is presented in Figure 5. On the y-axis, the figure shows NDVI, NIR and red values from Landsat-8 images (used as a reference) with accuracy, precision and uncertainty (Eq. 2–4) being calculated between Sentinel-2A and Landsat-8 for each bin.

Cloud and shadow screening for Landsat-8 and Sentinel-2A scenes was performed using the Fmask algorithm ^[27] and inversion residuals from aerosol optical thickness (AOT) estimation ^[26] (Figure 6). The pixels identified as those with high aerosol content were also masked out. Images from Sentinel-2A/MSI were further converted to 30 m to match spatial resolution of Landsat-8/OLI. Since atmospheric correction for Sentinel-2A was performed at 10 m spatial resolution for all spectral bands, conversion to 30 m was carried out by aggregation (averaging).

It was found that Landsat-8/OLI and Sentinel-2A/MSI exhibit misregistration issues ^[29]; therefore, additional co-registration was performed to ensure spatial consistency between the datasets ^[30]. Finally, NDVI was calculated for Landsat-8 scenes using band 5 (near-infra red-NIR) and band 4 (red), and for Sentinel-2A scenes using band 8A (NIR) and band 4 (red) using the following equation ^[31]: NDVI = (NIR-Red)/(NIR + Red).

---

2.3. Meteorological Datasets

We used air temperature derived from the NASA's Modern-Era Retrospective analysis for Research and Applications (MERRA2) reanalysis data set ^[32] to compute growing degree days (GDD) for winter wheat. The data are provided on a regular grid that has 576 points in the longitudinal direction and 361 points in the latitudinal direction, corresponding to a resolution of 0.625 × 0.5. We used the time-averaged, two-dimensional data collection (M2SDNXSLV), to extract daily averaged 2-meter air temperature (T2MMEAN). The data sets were extracted from the netCDF format, transformed to the geo-referenced raster, subset for study areas and linearly interpolated to the Landsat 30 m spatial resolution.

---

3. Methodology

Winter wheat yield mapping and assessment at regional scale consists of the two major steps: (ⅰ) winter crop mapping; (ⅱ) yield assessment at 30 m spatial resolution. Figure 7 illustrates all processing steps along with the input datasets. These steps are described in detail in the following sub-sections.

---

3.1. Winter Crop Mapping

For winter crop mapping, we adopted a previously developed approach for MODIS ^[33] that allows automatic mapping of winter crops using a priori knowledge on crop calendar and without using reference (ground truth) data. The method is based on per-pixel estimation of the peak NDVI (hereafter referred as the metric) during early spring (or early fall depending on the Earth hemisphere), when winter crops have developed biomass, while other crops (spring and summer) have no biomass in that time period. The calculated metric will have high NDVI values for winter crops and low NDVI values for other crops (Figure 8). Then, the metric is fitted using a Gaussian mixture model (GMM) ^[34] to automatically discriminate different crop types (winter versus others). The GMM is a linear combination of Gaussian distributions that can model any continuous distribution:

$ p\left( {\bf{x}} \right) = \sum\nolimits_{k = 1}^K {{\pi _k}} N\left( {{\bf{x}}{\rm{|}}{\mu _k},{\sum _k}} \right), $

(1)

where each Gaussian density $N\left ({{\bf{x}}{\rm{|}}{\mu _k}, {\sum _k}} \right)$ is called a component of the mixture and has its own mean ${{\mu _k}}$ and covariance ${{\sum _k}}$; parameters ${{\pi _k}}$ are weight (mixing) coefficients with $\sum\nolimits_{k = 1}^K {{\pi _k}} = 1$.

Parameters of the GMM model are estimated using an expectation-maximization (EM) algorithm that is run for all pixels identified as cropland. In our study, we used a cropland layer derived from the land cover map generated for Ukraine at 30 m spatial resolution ^[35]. The constraint to utilize cropland pixels only comes from potential confusion with grassland, hay, bulrush that might also have already developed biomass within the indicated time period. The component with the largest mean, i.e. NDVI value, in the obtained GMM model is considered to belong to the winter crop class (Figure 8). Finally, the derived GMM model is applied to all cropland pixels, and a posteriori probability (Eq. 1) of the pixel belonging to the winter crop class is estimated in the final resulting map. Pixels, with the probability larger than 0.5, are considered as winter crops.

---

3.2. Winter Wheat Yield Mapping and Assessment

Peak NDVI estimated on a per-pixel basis from a stack of Landsat-8/OLI and Sentinel-2A/MSI images from March to June was selected as a primary parameter for assessing winter wheat yield. In multiple studies, the seasonal peak NDVI has been shown to be strongly correlated with yields for a variety of crop types ^[5,8,9,10]. Since there are no available historical data for combination of Landsat-8 and Sentinel-2A images to correlate with yield measurements and build a crop yield model at district scale, we used a MODIS-derived winter wheat yield model that was calibrated for US and directly applied for Ukraine ^[7,10]. More specifically, the model takes advantage of daily MODIS data at Climate Modeling Grid (CMG) scale at 0.05 resolution to capture an NDVI peak and correlate with the yield. However, since proportions of winter wheat are variable within the CMG pixels, the model establishes a generalized relationship between the slope of NDVI against yield and pixel purity ^[10]: s = 9.61–0.05*m, where m is the winter wheat proportion at CMG scale (from 0 to 100%), and s is the slope such as yield = s*NDVI.

In case of Landsat-8–Sentinel-2A images, we can assume that purity at 30 m level is 100%, i.e. m = 100. Therefore, we obtain the slope of 4.61 to be applied to an NDVI peak calculated from the combination of Landsat-8 and Sentinel-2A data to map winter wheat yield at 30 m resolution.

Therefore, the MODIS-derived coarse resolution (0.05) winter wheat yield model, that was calibrated for Kansas (US) ^[10], is downscaled using winter wheat purity as a proxy to derive the slope between the peak NDVI and yield at 30 m resolution. This slope (4.61) is directly applied to the peak NDVI calculated from the stack of Landsat-8–Sentinel-2A images to derive a winter wheat yield map at 30 m resolution. These are used to estimate district-level yields by averaging yields at 30 m resolution over winter crop masks (section 3.1) for each district. In addition to the average, a standard deviation and coefficient of variation (CV), defined as a ratio between the standard deviation and the mean, is estimated as well. The estimated district-level yields are validated using independent reference data (section 2.1) collected at district level in Kirvohradska oblast (Ukraine) in 2016.

To improve peak NDVI estimation, we applied the GDD-based approach developed by Franch et al. GDD is used as a proxy to predict an NDVI peak using historical relationship between NDVI and GDD. GDD is calculated as the average daily maximum (T_max) and minimum temperatures (T_min) minus a base temperature (T_base) $GDD = \frac{{{T_{max}} + {T_{min}}}}{2} -{T_{base}}$, where, GDD = 0 if [(T_max + T_min)/2 < T_base], and with T_base= 0. Daily GDD is used to calculate accumulated GDD starting from the biofix date which was set to January, 1. We refer the reader to ^[7] for the details of this approach.

---

3.3. Validation Metrics

For comparison of satellite-derived winter crop areas and winter wheat yields with reference datasets at district level, we used the APU analysis metrics ^[28]:

•    accuracy (A) that shows the average bias of the estimates

$ A = \frac{1}{N}\sum\nolimits_{i = 1}^N {\left( {{P_i} - {O_i}} \right)} , $

(2)

•    precision (P) that shows repeatability of the estimates

$ P = \sqrt {\frac{1}{{N - 1}}\sum\nolimits_{i = 1}^N {{{\left( {{P_i} - {O_i} - A} \right)}^2}} } , $

(3)

•    uncertainty (U) that is the root mean squared error

$ P = \sqrt {\frac{1}{{N - 1}}\sum\nolimits_{i = 1}^N {{{\left( {{P_i} - {O_i} - A} \right)}^2}} } , $

(4)

•    relative uncertainty (rU) normalized by an average of reference values:

$ rU\left( \% \right) = \frac{U}{{\frac{1}{N}\sum\nolimits_{i = 1}^N {{O_i}} }} \times 100\% , $

(5)

where $P$ and ${{O_i}}$ are computed (from satellites) and observed (from reference) values, respectively.

---

4. Results & Discussion

---

4.1. Winter Crop Mapping

The GMM approach to winter crop mapping was applied to the peak NDVI calculated for the time period from March 1 to April 6 using a combination of Landsat-8 and Sentinel-2A, as well as using each of them separately. This was done in order to assess an added value of the combined use of these datasets. The indicated period (March 1 to April 6) was selected in such a way to capture NDVI development of winter crops and avoid confusion with early spring cereals that were planted beginning of March in 2016. Unfortunately, capturing peak NDVI during the emergence in late autumn of the previous year (e.g. during November) usually does not improve mapping of winter crops because of: (ⅰ) considerable cloud cover and unavailability of cloud-free imagery in that time period; (ⅱ) discrepancy of emergence state when much of winter crops have low NDVI. The derived maps were used to calculate the area of winter crops at districts level by pixel-counting. These estimates were compared to reference values and are presented in Table 1 and Figure. 9. The derived winter crop map using Landsat-8 and Sentinel-2A is illustrated in Figure 10.

Combination of Landsat-8 and Sentinel-2A allowed us to achieve R² = 0.9 and relative uncertainty of 11.6% when estimating winter crop areas at district level. When comparing to reference ground measurements, the accuracy of identifying winter wheat fields was 94.1%. It should be noted that these results were achieved in an automatic way utilizing knowledge on crop calendar and without utilizing any ground truth data. The use of Landsat-8 images only did not produce satisfactory results (R² = 0.64 and relative uncertainty of 32.7%) because of unavailability of cloud-free images early spring especially in the eastern districts of the oblast whereas the use of Sentinel-2A yielded R² = 0.88 and relative uncertainty of 13.5%. Overall, these results demonstrate the benefits, in a quantitative way, of the combined use of Landsat-8 and Sentinel-2A satellites comparing to the single-satellite usage.

---

4.2. Winter Wheat Yield Mapping

Comparisons of the estimated winter wheat yields at district level are presented in Table 2 and Figure 11.

As with winter crop areas, the combination of Landsat-8 and Sentinel-2A outperformed the single satellite usage in terms of APU metrics and R². When using either Landsat-8 or Sentinel-2A, the peak NDVI approach underestimated official statistics by -0.48 t/ha and -0.34 t/ha, respectively, while their combination improved accuracy to -0.17 t/ha. In terms of uncertainty, the peak NDVI approach for the Landsat-8–Sentinel-2A combination provided 0.31 t/ha (7.7%) whereas those values were 1.8 times higher for the Landsat-8 usage only (0.57 t/ha, 14.3%) and 1.5 times higher for the Sentinel-2A usage only (0.46 t/ha, 11.5%). These results clearly demonstrate the importance of higher observation frequency achieved with combination of Landsat-8 and Sentinel-2A satellites comparing to the single use. An example of the map showing spatial variability of estimated winter wheat yields at field scale is show in Figure 12.

The results presented in Figure 11(A) were further analyzed for errors. Overall, the points can be divided into 3 groups. The first group is composed of 3 points (shown in orange) representing districts with official statistics yields values close to 4 t/ha and underestimated by the peak NDVI approach. These districts feature relatively large values of CV of 21% whereas the average CV for all other districts is approximately 13%. The reason for that is smaller number of images available for these districts (mainly in the eastern part) which reduces ability to capture the peak NDVI. The second group (shown in red) is composed of districts with official statistics yields larger than 4 t/ha and the peak-NDVI approach underestimating it. The reason for that is saturation of NDVI occurs and the proposed approach fails to discriminate yield values at this level. Figure 13 shows an example of NDVI time-series from Landsat-8 and Sentinel-2A satellites for the districts with reference yields of 4.3 t/ha and 3.4 t/ha, and estimated yields of 4.04 ± 0.40 t/ha and 3.65 ± 0.64 t/ha, respectively, by the peak-NDVI approach. For the district with a higher yield value, NDVI quickly becomes 0.8 on April 29 (day of the year (DOY) 120) and not changing considerably (within 0.8–0.9) during the following days 50 days (until June 18 or DOY = 170). The NDVI values start to decrease when the senescence phase occurs and the crop is eventually harvested. This plot also shows the importance of the integration of both datasets: when using just Sentinel-2A data, we miss the peak.

The third group (shown in green) involves 8 districts with moderate yield values of up to 4 t/ha. The proposed approach is able to explain variations in the winter wheat yield (R² = 0.8) giving a bias of 0.1 t/ha and uncertainty of U = 0.13 t/ha (3.5%).

Additional experiments were performed to explore the effect of using GDD to predict the peak NDVI. In general, adding GDD improved estimates as, for example for the LC8-S2A case, the relative uncertainty decreased from 7.7% to 6.5% and R² increased from 0.45 to 0.5. However, adding GDD to the single satellite did not reach the performance of the combined LC8-S2A use without GDD (Table 2). This once again highlights the importance of the more dense time-series of LC8-S2A.

---

5. Conclusions

This study attempted to explore the combined use of Landsat-8 and Sentinel-2A satellites to winter crop mapping and winter wheat yield assessment at regional level. For both problems, the increased frequency of observations from the Landsat-8 and Sentinel-2A satellites was critical as it allowed us to achieve better performance comparing to the single satellite usage. For winter crop mapping, we adopted a previously developed approach for MODIS that allowed automatic winter crop mapping taking into account a priori knowledge on crop calendar without utilizing ground reference data. When comparing to official statistics on winter crop harvested areas, this approach gave R² = 0.9 and relative error of 11.6%. These results are encouraging as with little data inputs (crop calendar and cropland mask) and high temporal resolution of Landsat-8–Sentinel-2A satellites, it would allow the creation of winter crop maps at global scale at 30 m resolution.

For winter wheat yield mapping, we downscaled the generalized empirical model that is based on peak NDVI approach and implemented using MODIS data, and directly applied this model to the Landsat-8–Sentinel-2A images. Overall, the downscaled peak-NDVI approach with combined use of Landsat-8 and Sentinel-2A images gave uncertainty of 0.31 t/ha (7.7%) and R² = 0.45 substantially outperforming Landsat-8 only (1.8 times in terms of uncertainty) and Sentinel-2A only (1.5 times). The model was efficient in explaining moderate yield values ( < 4 t/ha) with R² = 0.8; however, it failed to capture the variance of high yield values (> 4 t/ha) due to NDVI saturation.

---

Conflict of interest

Both authors declare no conflicts of interest in this paper.

---