# Remote sensing based yield monitoring: Application to winter wheat in United States and Ukraine

B. Franch, E.F. Vermote, S. Skakun, Jean-Claude Roger, I. Becker-Reshef, E.

Murphy, C. Justice

► To cite this version:

B. Franch, E.F. Vermote, S. Skakun, Jean-Claude Roger, I. Becker-Reshef, et al., Remote sensing based yield monitoring: Application to winter wheat in United States and Ukraine. International Journal of Applied Earth Observation and Geoinformation, 2019, 76, pp.112-127. 10.1016/j.jag.2018.11.012. hal-03469851

HAL Id: hal-03469851

https://uca.hal.science/hal-03469851v1

Submitted on 11 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 **Remote sensing based yield monitoring: application to winter wheat in**

2 **United States and Ukraine**

3 B. Franch¹,², E.F. Vermote², S. Skakun¹,², J.C. Roger¹,², I. Becker-Reshef¹, E. Murphy¹,²,

4 C. Justice¹

5

6 ¹University of Maryland, Dept. of Geographical Sciences

7 ²Terrestrial Information Systems Laboratory, NASA Goddard Space Flight Center

8

9 **ABSTRACT**

10 Accurate and timely crop yield forecasts are critical for making informed agricultural
11 policies and investments, as well as increasing market efficiency and stability. Earth
12 observation data from space can contribute to agricultural monitoring, including crop
13 yield assessment and forecasting. In this study, we present a new crop yield model based
14 on the Difference Vegetation Index (DVI) extracted from Moderate Resolution Imaging
15 Spectroradiometer (MODIS) data at 1km resolution and the un-mixing of DVI at coarse
16 resolution to a pure wheat signal (100% of wheat within the pixel). The model was
17 applied to estimate the national and subnational winter wheat yield in the United States
18 and Ukraine from 2001 to 2017. The model at the subnational level shows very good
19 performance for both countries with a coefficient of determination higher than 0.7 and a
20 root mean square error (RMSE) of lower than 0.6 t/ha (15-18%). At the national level for
21 the United States (US) and Ukraine the model provides a strong coefficient of
22 determination of 0.81 and 0.86, respectively, which demonstrates good performance at
23 this scale. The model was also able to capture low winter wheat yields during years with
extreme weather events, for example 2002 in US and 2003 in Ukraine. The RMSE of the model for the US at the national scale is 0.11 t/ha (3.7%) while for Ukraine it is 0.27 t/ha (8.4%).

27

## 1. INTRODUCTION

Agriculture faces major challenges in this century. According to the Food and Agriculture Organization of the United Nations (FAO), global food production will need to increase by 70% by 2050 in order to meet the growing global demand [1]. However, despite tremendous improvements in technology and crop yield potential, food production remains highly dependent on climate [2]. Since Earth Observation (EO) satellites became available, it was recognized that the technology held considerable promise for agricultural monitoring [3]. The NASA Large Area Crop Inventory Experiment (LACIE) and Agriculture and Resources Inventory Surveys Through Aerospace Remote Sensing (AgriSTARS) programs made significant advances in crop monitoring, but were seriously constrained by satellite data availability. With recent advances in sensing technology, observation frequency, data products, computational analytics models and information technology, we are now better able to forecast and estimate crop production. Within-season crop yield assessment and forecasting is one of the major components of crop production monitoring. Timely and accurate information on crop yields at global, national, and regional scales is extremely important for informing a range of agricultural and food security decisions. For example, at the global and national level, timely information about crop development and production prospects
46 in food exporting countries can help to enable actions to mitigate the effects of any food
47 supply crisis (URL 1).

48 Observations from the NASA Earth Observing System (EOS) Moderate Resolution
49 Imaging Spectroradiometer (MODIS) sensors have several of the key qualities needed for
50 global crop yield monitoring such as global coverage at coarse spatial resolution (250m),
51 an excellent temporal resolution of two images every day that combined, minimize the
52 data gaps associated with clouds, and a suite of validated products. A good temporal
53 resolution is required, when working on agriculture applications given that the crop
54 phenology and conditions (water supply, pests, environmental...) can change very
55 quickly. Additionally, the MODIS dataset of 17+ years allows us to build a strong
56 empirical model based on a large number of statistics and covering a wide range of crop
57 conditions relative to other available satellite data. Generally, quantitative crop yield
58 forecast regression models are based on EO-derived vegetation indices (e.g., maximum
59 Normalized Difference Vegetation Index (NDVI) for the season [4] or temporal NDVI
60 integration [5]–[9], senescence rate [10] or combining land surface temperature (LST)
61 and NDVI [11]).

62 However, the main disadvantage of coarse to moderate resolution sensors is the spatial
63 resolution that often mixes, in a given pixel, signals from different land cover types and
64 crops. Stratifying a region into different crop types (commonly termed as crop masking)
65 is an important step in developing EO-based yield models [12]. Such masks enable the
66 isolation of the remotely sensed, crop-specific signal throughout the growing season,
67 reducing the noise on the signal from other land cover or crop types. In the United States
68 (US), the US Department of Agriculture (USDA) generates a yearly national Cropland
69 Data Layer (CDL) since 2007 [13] and Canada provides yearly national Annual Crop 70 Inventory Maps (ACIM) since 2009 [14]. However, these masks are provided at the end 71 of the growing season and there are no available crop type masks for other countries. As 72 a consequence, generally EO-based crop yield models use static cropland or crop type 73 masks. Next, we briefly describe some examples of studies found in the literature that use 74 coarse resolution EO data to monitor crop yields, how they deal with the signal mixing 75 and which masks they use as reference. 76

76 Becker-Reshef et al. (2010) [4] developed an empirical, but generalized crop yield model 77 based on MODIS Climate Modeling Grid (CMG, 0.05 degree spatial resolution) data and 78 applied it to Kansas and Ukraine. The model, that was calibrated in Kansas, is based on 79 the relationship between the NDVI value at the peak for the season and the final yield 80 value, corrected for the average purity or the percentage of crop of the 5% purest pixels 81 within the area studied. They used static crop type masks by aggregating the USDA CDL 82 of 2007 to 0.05 degrees spatial resolution and developed their own mask for Ukraine. 83

83 Later, Franch et al. (2015) [15] improved the temporal aspect of the model by including 84 Growing Degree Days (GDD) information. The model was applied satisfactorily (errors 85 lower than 10%) to monitor wheat yield at the national level over the US, Ukraine and 86 China. It also showed good performance, when applied to the Advanced Very High 87 Resolution Radiometer (AVHRR) Long Term Data Record (LTDR) [16]. However, this 88 model was designed to work primarily at the national level. 89

89 Kowalik et al. (2017) [17] developed an empirical model to forecast wheat yield over 90 Europe. They used SPOT-VEGETATION products of NDVI and Fraction of Absorbed 91 Photosynthetically Active Radiation (FAPAR) at 1 km spatial resolution with errors
ranging from 10-17% depending on the location. They divided Europe into different agro-climatic zones and aggregated those pixels using as a reference the arable land fraction from the Corine Land Cover (CLC, [18]) and estimated the weighted average of pixels with arable land fraction land higher than 50%. Later, López-Lozano *et al.* (2015) [19] analyzed the correlation between EO indicators (FAPAR from SPOT-VEGETATION at 1km spatial resolution) and wheat, barley and grain maize yields in Europe. In this case they divided Europe into different sub-national areas depending on the available statistics for each country. In this work, they also analyzed the weighted average of pixels within the 90th percentile arable lands fraction for each administrative unit, using the CLC as a reference. They found that crops growing under moderately or strongly water-limited conditions present a high inter-annual variability of leaf area development and senescence that influences the final yield, which is observable by EO systems. However, in agro-climatic areas, where crop growth is not frequently water-limited (e.g. irrigated crops), yield variations are explained by a complex interaction of different factors not necessarily related to the crop leaf area dynamics, and, therefore, weather indicators are needed to complement the FAPAR time series.

Kogan *et al.* (2013) [20] assessed the feasibility and relative efficiency of using NDVI MODIS data to forecast winter wheat yield at oblast (sub-national) level in Ukraine. They spatially aggregated each oblast using a cropland map as a reference (Rainfed croplands class) that was extracted from the ESA GlobCover map at the 300 m resolution. They concluded that NDVI values 2-3 months prior to harvest allowed reliable yield forecasts at the oblast level with RMSE ranging from 0.5 to 0.8 t/ha.
Anderson *et al.* (2016) [21] analyzed the correlation between the Evaporative Stress Index (ESI) from MODIS at 1km spatial resolution and yield for major crops grown in Brazil. They aggregated those pixels, classified as majority agricultural land, using a land use map of Brazil, generated by the Brazilian Geographical and Statistical Institute (IBGE). Their results showed better correlation than other indices such as LAI for most crops and regions. Additionally, in the case of soybeans, the yields showed earlier response (10 to 25 days) to extreme events than LAI.

In this study, we present a new empirical model based on remote sensing data at coarse spatial resolution (MODIS aggregated to 1km). This model is based on the wheat signal un-mixing from the signal of other surfaces, surrounding the wheat fields by using yearly crop type masks. Additionally, based on the analysis of the un-mixed time series of the wheat, we use three different regressors that we calibrate for each administrative unit against the historical yield statistics. The main goal of developing this model is to make it sensitive to extreme weather events and applicable at the subnational level. This model is based on the Difference Vegetation Index (DVI) [22], that is estimated by the simple difference between surface reflectance in the red and the Near Infrared (NIR) spectral bands. The following aspects are incorporated into the proposed model:

1. We use BRDF-corrected daily 1 km MODIS data derived from both Terra and Aqua satellites.

2. We un-mix wheat surface reflectance's to estimate the signal as if the MODIS pixels were 100% wheat. We use a satellite-derived amplitude of the peak DVI along with the time period that the peak remains (Length) and Evaporative
Fraction (EF) to correlate with winter wheat yields. The EF is a parameter related to the evapotranspiration and is further detailed in section 3.2.

In the next section (2.), we describe the study area and the data considered. Then, we describe the method proposed and how we derive the different input parameters (Section 3.). Finally, we show the results both at the subnational and the national levels for the two countries considered (Section 4.).

142

## 2. STUDY AREA AND DATA DESCRIPTION

### 2.1. Study area

The study is carried out over the US and Ukraine from 2001 until 2017. The US is one of the main producers and exporters of wheat globally. In 2016 the US was the #1 wheat exporting country, shipping 14.8% of global wheat exports (URL 2). Wheat is produced in almost every state in the United States and winter wheat varieties dominate US production, representing between 70-80% of the total wheat production. The winter wheat is planted in the Fall and harvested during June-July. Generally, wheat is rain-fed and just 7% of the national production is irrigated. The main wheat class is Hard Red Winter Wheat, which is grown primarily in the Great Plains, with Kansas being the largest producing state.

Ukraine is another critical player in the global wheat market. In 2016, Ukraine ranked #6 by exporting 7.2% of global wheat exports (URL 2). Wheat is grown all across the country, although the central and southern regions are the key growing areas (Forest-Steppe and Steppe zones in Figure 1). About 95% of Ukraine's wheat production is winter wheat, planted in the Fall and harvested during June-July of the following year.
Generally, wheat is rain-fed and is not irrigated. Ukraine produces mostly Hard Red Winter Wheat.

Tables 1 and 2 show the primary characteristics of the main climatic zones in the US and Ukraine, where most of the winter wheat fields are planted and the average yields for each area using as a reference 2016.

<table><caption>Table 1. Main climate zones in the US (classification from Vocke & Ali (2013) and yield values from NASS statistics).</caption><thead><tr><td>Zone</td><td>States</td><td>Climate conditions</td><td>Major Type of winter wheat</td><td>Average yield (2016), t/ha</td></tr></thead><tbody><tr><td>North Central</td><td>MT, ND, SD, MN</td><td>Humid continental</td><td>Hard Red</td><td>3.6</td></tr><tr><td>Southern Plains</td><td>OK, TX</td><td>Temperate climate</td><td>Hard Red</td><td>2.4</td></tr><tr><td>Central Plains</td><td>KS, CO</td><td>Temperate climate</td><td>Hard Red</td><td>3.5</td></tr><tr><td>Northern Plains</td><td>NE, WY</td><td>Humid continental</td><td>Hard Red</td><td>3.0</td></tr><tr><td>Northwest</td><td>WA, OR, ID</td><td>Oceanic</td><td>Soft White</td><td>5.0</td></tr></tbody></table>

<table><caption>Table 2. Main climate zones in Ukraine.</caption><thead><tr><td>Zone</td><td>Conditions</td><td>Average yield (2016), t/ha</td></tr></thead><tbody><tr><td>Plane-Polissya</td><td>Mixed forest zone in the north, 26% of the entire Ukrainian territory</td><td>4.36</td></tr><tr><td>Forest-Steppe</td><td>Mainly cultivated area in the center</td><td>4.88</td></tr></tbody></table>
<table><tr><td></td><td>(34%)</td><td></td></tr><tr><td>Steppe</td><td>Intensive cultivated area (40%)</td><td>3.64</td></tr></table>

<figure><img src="image_2.png" /><figcaption>Figure 1. Agro-climatic zones of Ukraine [20]</figcaption></figure>

## 2.2. Crop statistics on winter wheat

In each country, we run the wheat yield model for each Administrative Unit (AU). Thus, we work at different spatial scales depending on the data availability. In the case of the US, we work at the county level (average area of 258,000 ha), while in the case of Ukraine we work at the oblast level (average area of 2,414,000 ha).

For the US, we use the official archive of county-level statistics on yield, area harvested, and production that are available from the USDA National Agricultural Statistics Service (NASS) Quick Stats database (http://www.nass.usda.gov/Quick_Stats/). The NASS crop statistics are based on data obtained from multiple frame-based sample surveys of farm operators, objective yield surveys, agribusinesses, shippers, processors and commercial
184 storage firms.

185 For Ukraine, oblast-level crop statistics are obtained from the State Statistics Service of
186 Ukraine (SSS) (http://www.ukrstat.gov.ua) for winter wheat area harvested and yield.

187 These official statistics are based on farm surveys collected from all the agricultural
188 enterprises (large-scale farms that produce commodities exclusively for sale).

189

190 **2.3. Crop type masks**

191 For the US, a national winter wheat mask is available from the CDL produced by NASS
192 from 2008 to 2016. The CDL is a rasterized land cover map using field level training data
193 from extensive ground surveys, farmer reports provided to the US Farm Service Agency
194 (FSA), and remotely sensed data from Landsat Thematic Mapper (TM), Landsat
195 Enhanced Thematic Mapper (ETM+) and Advanced Wide Field Sensor (AWiFS). These
196 data are used in a decision tree classifier in order to produce a land cover classification
197 that distinguishes between different crop types, including winter wheat [13], [24]. In this
198 study, yearly CDL layers are used to identify winter wheat growing areas by estimating
199 the percentage of 30m CDL pixels within each MODIS 1km pixel. At 1km spatial
200 resolution MODIS pixels mix the signal from wheat with the signal from other surfaces
201 surrounding the wheat fields. In order to isolate the wheat signal through the un-mixing
202 process (see section 3.3), we used the CDL masks to identify not only the wheat but also
203 the other surfaces and how much percentage of each class is observed by each pixel.
204 Therefore, looking at the main land cover classes surrounding the wheat fields in the US,
205 we consider a total of eight different cover classes that are included in the un-mixing
206 method: winter wheat, spring wheat, corn, soybean, alfalfa, potato, grassland and forest.
Before 2008 the CDL was not produced at the national scale and has gaps for some states. To generate yearly winter wheat masks for US from 2001 to 2007 and for Ukraine for the 2001–2016 period, we used MODIS-derived 8-day composited NDVI values at 250 m spatial resolution [25] and an automatic phenology-based approach for winter wheat mapping [26]. This approach showed an accuracy of >90% for US and Ukraine, and good correspondence to official statistics with an average coefficient of determination >0.85.

Additionally for Ukraine, we extracted static grassland and forest masks from the 30 m land cover map for Ukraine for 2010 [27], [28]. The reported classification accuracies are 97% and 90% for forest and grassland, respectively. Note that these classes are just considered when they are surrounding the wheat fields, that is, when their signal is mixed with wheat in order to un-mix it. Thus, by using static masks of grassland and forests in Ukraine we assume that these surfaces surrounding the wheat fields have not changed considerably during the time period analyzed.

## 2.4. MODIS surface reflectance time-series

This study is based on MODIS daily surface reflectance Collection 6 data (M{OY}D09GQ) distributed by the Land Processes Distributed Active Archive Center (LP DAAC, https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table), which are gridded in the sinusoidal projection at 250m resolution. Additionally, we used the product M{OY}D09GA to extract the geometry of observation/illumination of each image. Since the nominal 250 m MODIS resolution decreases for the off-nadir
229 observations and due to inaccurate registration [29], we re-scaled the 250m surface
230 reflectance to 1km spatial resolution to mitigate that effect by aggregating 4x4 pixels.

231

232 **3. METHODS**

233 **3.1. General workflow**

234 The diagram in Figure 2 shows the workflow followed in this study. The model is based
235 on three different regressors (blue boxes in the diagram): the amplitude of the peak, the
236 time period that the peak remains for and the average of the Evaporative Fraction (EF)
237 during the thirty-day period following the peak. In order to determine these parameters,
238 the inputs needed (yellow boxes in the diagram) are the MODIS Surface Reflectance
239 product, the MODIS Land Surface Temperature product, the CDL crop mask in the US
240 and official statistics for model calibration and validation. In the following subsections
241 we describe more detail about the methods (green boxes in the diagram).

242
<figure><img src="image_3.png" /><figcaption>Figure 2. The workflow followed in this study.</figcaption></figure>

244

245

246

247

## 3.2. Surface reflectance normalization and albedo estimation

248

The wide swath MODIS sensor allows for near global coverage of the Earth every day.

249

However, it has a 16-day repeat cycle, which means that every day the geometry of
observation is different and can include View Zenith Angles (θv) of up to 65 degrees. As

250

a consequence, the surface reflectance that is defined for a given geometry of

251
observation-illumination has different values every day. When working with normalized vegetation indices such as the NDVI, the geometry effects are minimized, but that is not the case, when working with non-normalized parameters such as DVI. In order to normalize the BRDF effects on the surface reflectance, we used the VJB method [30], [31]. This method uses longer compositing periods (5 years in [30]), than the MCD43 product (16 days) [32], which reduces the noise in the normalized reflectance time series [29]. In this study, we derive the nadir BRDF parameters at 1km spatial resolution using the most recent five years (2012 to 2016). By using the daily surface reflectance (from both Aqua and Terra) and its angular conditions during the five-year period considered, we derive the variables that define the BRDF shape (V and R in Eq. 1) using the approach used by [31]. Then we apply Eq. 1 for deriving the normalized surface reflectance ($\rho^N$) in MODIS band 1 and 2, and we also use the BRDF coefficients for estimating the surface albedo:

$$ \rho^N(45,0,0) = \rho(\theta_s, \theta_v, \phi) \frac{1+VF_1(45,0,0)+RF_2(45,0,0)}{1+VF_1(\theta_s, \theta_v, \phi)+RF_2(\theta_s, \theta_v, \phi)}, \quad (1) $$

where $\rho$ is the directional surface reflectance, $\theta_s$ is the sun zenith angle, $\theta_v$ is the view zenith angle, $\phi$ is the relative azimuth angle, $F_1$ is the volume scattering kernel, based on the Ross-Thick function derived by [33], but corrected for the Hot-Spot process proposed by [34], $F_2$ is the geometric kernel, based on the Li-sparse model [35], but considering the reciprocal form given by [36], $V$ represents the volume BRDF parameter since it is linked to the Volume kernel and $R$ represents the roughness BRDF parameter since it is
linked to the geometric kernel. These parameters (V and R) represent the shape of the BRDF.

In order to estimate the surface albedo, that is needed to estimate the EF, we follow the same methodology to derive the surface albedo as the official MCD43 MODIS product [32].

Downwelling flux may be written as the sum of a direct component and a diffuse component.

Black-sky albedo ($\alpha_{bs}$) or directional-hemispherical reflectance (DHR) is defined as albedo in the absence of a diffuse component and is a function of solar zenith angle. White-sky albedo ($\alpha_{ws}$), or bihemispherical reflectance, is defined as albedo in the absence of a direct component, when the diffuse component is isotropic. It is independent of the geometry of illumination-observation. Therefore, the albedo can be written by integrals of the BRDF model through the black-sky albedo and the white-sky albedo.

$$ \alpha_{bs}(\theta_i, \lambda) = \sum_k f_k(\lambda)h_k(\theta_i), \qquad (2) $$

$$ \alpha_{ws}(\lambda) = \sum_k f_k(\lambda)H_k, \qquad (3) $$

where $f_k(\lambda)$ are the BRDF model parameters and

$$ h_k(\theta_i) = \int_0^{2\pi} \int_0^{\pi/2} k_k(\theta_i, \theta_v, \phi) \sin\theta_v \cos\theta_v d\theta_v d\phi, \qquad (4) $$

$$ H_k = 2 \int_0^{\pi/2} h_k(\theta_i) \sin\theta_i \cos\theta_i d\theta_i, \qquad (5) $$

where $k_k(\theta_i, \theta_v, \phi)$ are the BRDF model functions of the observation geometry ($F_1$ and $F_2$ in Eq. (1)).
The albedo, $\alpha(\theta, \lambda)$, under actual atmospheric conditions can also be modeled quite accurately as an interpolation between the black-sky (direct beam) albedo and white-sky (completely diffuse) albedo as a function of the fraction of diffuse skylight $S(\theta, \tau(\lambda))$, which, in turn, is a function of optical depth $\tau$ [37], [38].

$$ \alpha(\theta, \lambda) = [1 - S(\theta, \tau(\lambda))]\alpha_{bs}(\theta, \lambda) + S(\theta, \tau(\lambda))\alpha_{ws}(\lambda) \quad (6) $$

We simulate the diffuse skylight with 6S [39] using as an input the MODIS Aerosol Optical Thickness (AOT) at 550nm product, which is included in the Climate Modeling Grid (CMG) product.

### 3.3. Un-mixing of the wheat signal

When working with medium to coarse spatial resolution sensors, the signal retrieved by each pixel generally mixes different land cover classes. Therefore, crop yield modeling requires availability of crop-specific masks, so that only the remotely sensed signal from the crops in question will be considered. In order to work with the purest wheat signal possible, in this work we present a method based on the un-mixing of the signal at 100% wheat purity. For each AU and for a given date, the total DVI signal from each pixel, $i$, can be written as the sum of the DVI signal from the wheat ($DVI_{wheat}$) multiplied by the percentage of wheat within the pixel or wheat purity ($Wpct$) and the DVI from other surfaces within the pixel ($DVI_{others}$) multiplied by the remaining percentage ($1 - Wpct$):.

$$ DVI_{i} = DVI_{i,wheat} \cdot Wpct_{i} + DVI_{i,others} \cdot (1 - Wpct_{i}) \quad (7) $$

Eq. (7) has the structure of a linear regression function and can be solved with the Ordinary Least Squares (OLS) method. Figure 3 shows an example of the linear
regression between the DVI<sub>i</sub> and the wheat purity, Wpct<sub>i</sub>, (black dots) for the Day of the Year (DOY) 79 of 2016 for Harper County (Kansas).

<figure><img src="image_4.png" /><figcaption>Figure 3. Linear regression of the DVI<sub>i</sub> versus the wheat purity (Wpct<sub>i</sub>) through the Harper County for DOY 79 of 2016.</figcaption></figure>

In the case of the US, where the CDL product provides masks for other crops, we can
take advantage of such information to improve the un-mixing of the signal. As described
in the previous section, after studying the main classes that generally surround the US wheat fields, we consider a total of 8 different classes plus the “others” contribution. The eight classes considered are: winter wheat (#1), spring wheat (#2), corn (#3), soybean (#4), alfalfa (#5), potato (#6), grassland (#7) and forest (#8) with crop purities (percentage of a given crop area within the pixel) Cpct (corn). Then, equation 7 is written as:

$DVI_i = \sum_{Y=1}^{8} DVI_{i,Y} \cdot Ypct_i + DVI_{i,others} \cdot (1 - \sum_{Y=1}^{8} Ypct_i), \quad (8)$

where Y is the class number, Ypct_i is the purity of class Y (or percentage of class Y area within the pixel). This equation also has the structure of a linear regression function.
This un-mixing method is applied to each AU and day, isolating in each case the DVI signal coming from wheat. Based on the time series of the DVI_wheat, we estimate the amplitude of the peak value of this variable and the time period that this peak remains (hereafter referred as length), defined as the number of days when the DVI_wheat is higher than 80% of the amplitude value.

## 3.4. Evaporative Fraction estimation

Any change in the water resources of the crop takes some time to show an effect on the biomass, that is directly linked to the vegetation indices. Therefore, any water stress condition happening during or after the peak might not be accounted for, when using just the amplitude and length. However, these changes can be detected earlier, when using thermal indices [21]. One of the key factors that drive the land surface energy process is
345 evapotranspiration (ET). ET is defined as the flux of water evaporated at the earth- 346 atmosphere interface (from soil, water bodies and interception) and transpired by 347 vegetation through stomata in their leaves as a consequence of the photosynthetic 348 processes. There are several ET models and products [40]–[46]. However, they generally 349 rely on the estimation of several parameters such as net radiation, the sensible and latent 350 heat flux and the ground flux. In this study, we examine the use of the Evaporative 351 Fraction (EF), as one of the inputs to the proposed yield model. The EF describes the 352 ratio of evapotranspiration to available total energy (equation 9).

$$EF = \frac{LET}{LET + H} = \frac{LET}{R_n - G}, \qquad (9)$$

354 where LET is the latent heat flux, H is the sensible heat flux, R_n is the net radiation and G 355 is the soil heat flux. In this study, EF is derived based on the Simplified Surface Energy 356 Balance Index (S-SEBI) method developed by [44]. The surface temperature and 357 reflectance of areas with constant atmospheric forcing are correlated and the relationships 358 can be applied to determine the effective land surface properties [40], [47]. The S-SEBI 359 method uses as reference the LST vs albedo diagrams to determine the dry and the wet 360 conditions of the area considered. The dry conditions are defined by a baseline 361 representative of the maximum sensible heat with a reflectance dependent maximum 362 temperature. Meanwhile, the wet conditions are defined by a baseline of the maximum 363 potential evapotranspiration with a reflectance dependent minimum temperature, hence 364 the sensible and latent heat fluxes are partitioned according to the actual surface 365 temperature. The assumption of this method is that global radiation and atmospheric 366 conditions tend to be constant and that there are a sufficient number of wet and dry pixels
to determine and right wet and dry edge [44]. A detailed analysis on these baselines
estimation uncertainty can be found in [48]. The EF is derived using an approximation
according to [44], [49], [50]:

$$
371 \quad EF = \frac{T_H - T_S}{T_H - T_{LE}} = \frac{a_H \alpha + b_H - T_S}{(a_H - a_{LE})\alpha + (b_H - b_{LE})}, \quad (10)
$$

372
373 where $T_H$ is land surface temperature for dry pixels; $T_{LE}$ is land surface temperature for
374 wet pixels; $T_S$ is land surface temperature; $\alpha$ is surface albedo; $a_H$, $a_{LE}$ are the slopes of
375 the line of, respectively, the high and low temperatures as a function of surface albedo;
376 $b_H$, $b_{LE}$ are the intercepts of the line of, respectively, the high and low temperatures as a
377 function of surface albedo.
378

378 In this study, the LST is extracted from the official MODIS (MYD11) LST product [51].
379 Note that the EF is estimated using the Aqua satellite, whose overpass time is during the
380 afternoon, when higher water stress can be detected. Section 3.2 describes the method
381 followed to derive the surface albedo. We derive the LST vs albedo diagram parameters'
382 for each day at the state level in the US and at the oblast level for Ukraine. Note that in
383 the case of the US, the LST vs albedo diagrams are built at state level (instead of county
384 level) to ensure having enough pixels with wide variability of moisture conditions. Figure
385 4 shows an example of the LST-albedo diagram in Kansas. Then, for each pixel, Eq. 10 is
386 implemented, providing a pixel-based EF estimation. Finally, at the county level we
387 apply the same un-mixing procedure as the DVI in order to isolate the EF for the wheat.
388 This parameter is considered in the model by estimating its average during the 30 days
389 after the seasonal peak (grain filling period). We consider this period to account for any
water stress conditions happening during or after the peak and that will have a negative impact on the final yield [52], [53].

392

<figure><img src="image_5.png" /></figure>

393

Figure 4. Example of a LST-albedo diagram in Kansas on April 15th, 2016. The red lines represent the maximum sensible heat (dry boundary, top line) and the maximum potential evapotranspiration (wet boundary, bottom line). Red dots represent the data considered to derive the red lines.

398

399 **3.5. Evaluation of the un-mixing method**
We test the performance of the un-mixing of the wheat signal by applying the un-mixing method to all Landsat 8 OLI cloud-free data available from 2013 to 2016 in three different counties in the US: Harper (Kansas, 17 images), Walla Walla (Washington, 23 images) and Bingham (Idaho, 14 images). The validation of the method is just focused on Landsat data to avoid any impact of spectral adjustment between two sensors or any error associated with angular effects correction. Thus, first Landsat data are transformed to MODIS-like images by upscaling the spatial resolution of Landsat to 1km and re-projecting them to the sinusoidal projection. Then, using the same CDL crop masks that we use for MODIS, we apply the un-mixing method. Finally, we validate the un-mixing of the wheat signal against the average throughout the county of all the 100% purity wheat pixels at the original Landsat spatial resolution.

411

## 3.6. Calibration of the yield model

For each AU, a separate linear regression was built with different combinations of regressors (Equations 11-16).

$$415 \quad yield = a_{10} + a_{11} \cdot A \quad (11)$$

$$416 \quad yield = a_{20} + a_{21} \cdot L \quad (12)$$

$$417 \quad yield = a_{30} + a_{31} \cdot EF \quad (13)$$

$$418 \quad yield = a_{40} + a_{41} \cdot A + a_{42} \cdot L \quad (14)$$

$$419 \quad yield = a_{50} + a_{51} \cdot A + a_{52} \cdot EF \quad (15)$$

$$420 \quad yield = a_{60} + a_{61} \cdot A + a_{62} \cdot L + a_{63} \cdot EF \quad (16)$$

A p-value for each regressor was calculated which provides a test of the null hypothesis that a coefficient for this regressor is equal to zero, i.e. the regressor has no effect on the
dependent variable (crop yield in our case). A small *p*-value indicates that the null hypothesis can be rejected, meaning that there is a statistically significant relationship between the regressor and dependent variable. We used this statistic index to decide in each AU which variables were significant (by using a threshold of 0.2 in *p*-value) and which equation to consider (lowest *p*-value for each regressor).

There are several winter wheat varieties planted in the US (over 30,000) distributed over different regions with diverse climates, soils, treated with different fertilizers and susceptible to different diseases. Thus, specific modelling of such characteristics would be complex, data demanding, and leading to highly parametrized models. Therefore, our approach was to empirically calibrate the model at the county level assuming that homogeneity in the above-mentioned parameters. Additionally, given that we use 1km spatial resolution data whose wheat signal is un-mixed, we need to ensure that the wheat purity is sufficient to invert the model. At low wheat purities, the wheat signal is low and it is competing with the signals from other surfaces. In these circumstances, the un-mixing can lead to high errors. Therefore, the model is applied to every county with wheat purity values higher than 40%. After testing other purities, this 40% threshold showed the best performance or balance between number of counties considered (using a higher threshold reduces the number of counties where the model can be applied) and noise in the extrapolated data (when using lower thresholds). During the calibration process, none of the three regressors considered in the proposed model (amplitude, length and EF) was able to capture the late frost that mostly affected Kansas and Oklahoma in 2007. This late frost happened during early April and is the only event with these characteristics that occurred during the 17 years analyzed. As a consequence, there are
446 not enough data to study the inclusion of another variable in the model capable to
447 respond to these extreme conditions. Following the approach of [4], every county in these
448 two states during 2007 is excluded from the calibration and the validation processes.

449 In Ukraine we follow the same methodology as in the US, by calibrating the crop yield
450 model individually each AU (oblast).

451

452 **3.7. Validation of the yield model**

453 Each country is validated both at national and subnational level from 2001 to 2017. Note
454 that the period 2001 to 2016 is used both for calibration and validation, while 2017 is just
455 used for validation. Additionally, we applied at subnational level a leave-one-out cross-
456 validation (LOOCV) procedure, in which iteratively a year is retained for validation,
457 while other years are used for model's calibration.

458 For the Ukraine, the model can be applied to every oblast. Therefore, the national number
459 match the addition of the oblast-level statistics. However, this is not the case for the US,
460 where there are some counties excluded from the study because the maximum purity of
461 the wheat signal is lower than 40% in any 1 km pixel due to small field size.
462 Additionally, in the official county-level statistics there are some counties with no
463 statistics or some temporal gaps in other counties whose data is not publicly available due
464 to privacy issues but are accounted at national level.

465

466 **4. RESULTS**

467 **4.1. Evaluation of the un-mixing method**
Figure 5 shows the evaluation of the un-mixing method using Landsat 8 as a reference.

The results show a good agreement of the extrapolated signal with a high coefficient of
determination (>0.96) and errors of around 0.02. As an example of the impact and utility
of the un-mixing approach, Figure 6 plots the average of every Landsat 1km pixel with
any wheat purity against the same x-axis as in Figure 5 (just averaging those pixels with
wheat purity equal to 100%). Given that these coarse resolution pixels mix different
surfaces' DVI, this variable does not show any variability throughout the year not
allowing any phenological analysis.

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5. Evaluation of the un-mixing method using Landsat 8 as a reference for Harper County, Kansas (top left); Walla Walla County, Washington (top right); and Bingham County, Idaho (bottom). The error bars represent the standard deviation of the wheat pixels' average at the Landsat spatial resolution.</figcaption></figure>

<figure><img src="image_8.png" /></figure>

482
483 Figure 6. Example of the impact of not using the un-mixing method for Harper County, Kansas.

484 The error bars represent the standard deviation of the wheat pixels' average at the Landsat spatial

485 resolution.

486

487 **4.2. Calibration of the yield model**

488 **4.2.1. The US**

489 Figure 7 shows the p-values associated with each input statistical significance in the

490 calibration equations of the model and the determination coefficient for each county.

491

<figure><img src="image_9.png" /></figure>

a)

<figure><img src="image_10.png" /></figure>

b)

492

493

494

495
<figure><img src="image_11.png" /></figure>

c)

<figure><img src="image_12.png" /></figure>

d)

496

497

498

499

500

501

502

Figure 7. Geographical distribution of *p*-values derived for a) the amplitude of the peak, b) length of the peak and c) EF. Corresponding coefficients of determination for each county is shown in subplot (d).

503

504

505

506

507

508

The amplitude is used in most counties and shows generally p-values lower than 0.05 in the most important producing counties along the Great Plains, the North-West, North Central regions and Michigan. However, it shows values higher than 0.10 in California.

The length is the second most used regressor and generally shows values lower than 0.05 in Kansas and South Dakota. Finally, the EF is the least used regressor and it is mostly
used in the Great Plains area where most counties show low p-values. The coefficient of
determination shows that the model works very well ($R^2>0.6$) along the Great Plains and
some norther counties. However, the coefficient of determination is low in California
(2% of national production) and the counties around the Mississippi (3% of national
production). California winter wheat is one of the few states that plants most classes of
winter wheat such as hard red, hard white, soft white or durum. This heterogeneity of
classes maybe the reason of the low coefficients of determination in this area. Regarding
Mississippi and Arkansas counties, they generally show low purities (slightly above
40%). This impacts the accuracy of the un-mixing given that the wheat signal is
competing with other surfaces whose signal will be stronger.

4.2.2. Ukraine

<figure><img src="image_13.png" /></figure>

Figure 8 shows the p-value of each regressor considered in the model and the
determination coefficient in each oblast.
<figure><img src="image_14.png" /><figcaption>Figure 8. Geographical distribution of *p*-values derived for a) the amplitude of the peak, b) the length of the peak and c) EF. Corresponding coefficients of determination for each oblast is shown in subplot (d).</figcaption></figure>

The amplitude shows low values (<0.05) across the country, and only four oblasts exhibit *p*-values higher than 0.05. Both the length and the EF are found not statistically significant in Ukraine. The length shows low *p*-values (<0.05) in two oblasts (located in the center and north respectively), while the other *oblasts* show higher *p*-values (>0.10). Average EF *p*-values are generally relatively high (>0.10). The coefficient of determination shows very good correlation (R²>0.4) for most oblasts, especially in the southern area, where most of the wheat production is located.

## 4.3. Validation of the yield model

### 4.3.1. The US

Figure 9 (top left) shows the county level calibration and validation of the model in the US. Each point in the plot represents a given county and a given year during the time-period from 2001 to 2016. The statistics show a high correlation with a coefficient of determination of 0.86 and an error of 0.4 t/ha (15.36%). Figure 9 (top right) shows the
542 county level calibration and validation, excluding Kansas and Oklahoma in 2007. Though
543 the global statistics barely improve (14.52%), it removes from the plot the overestimation
544 of the low yields caused by the late frost that the model does not capture. Validation of
545 the model using data from the year of 2017 (the only year not used in the calibration)
546 shows performance consistent with the one derived for 2001 to 2016 validation with an
547 error of 15.67% and coefficient of determination of 0.86. These results show that the
548 model does not exhibit over-fitting. Figure 9 (bottom left) shows the county level
549 LOOCV. Analogously, Figure 9 (bottom left) shows the county level LOOCV, excluding
550 Kansas and Oklahoma in 2007. Figure 10 (left) shows the evaluation of the model by
551 adding up all the county-level statistics, where we can apply the model (all data in Figure
552 9 left) to the national scale and comparing it to the addition of the official county-level
553 statistics of the same counties considered in the model to national-level. Note that at the
554 national level we do not analyze the performance of 2017 separately. The model shows
555 good results for every year except 2007, when the yield is overestimated. Figure 10
556 (right) shows the national validation excluding Kansas and Oklahoma for 2007 since the
557 model is not able to capture the late frost, which improves the validation statistics.
558 If we compare the results of our model to the national statistics released by NASS (Figure
559 11 left) the coefficient of determination remains very similar (0.81 compared to 0.77), but
560 there is a constant negative bias (underestimation of the model) that we can correct by
561 dividing by the slope of the regression (Figure 11 right). Note that in both plots of Figure
562 11, the results of the model are excluded for Kansas and Oklahoma for 2007. Therefore,
563 in order to make 2007 comparable at the national level, we complete the model results
over the rest of the country with the official statistics of these two states. The statistics show a strong coefficient of determination of 0.81 and an error of 0.11 t/ha (3.73%).

566

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

567

Figure 9. County level validation in the US (top) from 2001 to 2016 (left) and excluding Kansas and Oklahoma in 2007 (right). Each black dot represents a county and a year

569
570 from 2001 to 2016 while each red dot represents a county in 2017. Bottom: LOOCV at 571 county scale from 2001 to 2017 (left) and excluding Kansas and Oklahoma in 2007 572 (right).

573

<figure><img src="image_19.png" /></figure>

574

575 Figure 10. National level aggregation of the county level statistics where we can apply 576 the model and when there are county level statistics available (left). Excluding Kansas 577 and Oklahoma in 2007 (right). 2017 is the black triangle.

578
<figure><img src="image_20.png" /><figcaption>Figure 11. National level validation comparing directly the model results to national numbers (left) and correcting for the bias (right). 2017 is the black triangle.</figcaption></figure>

### 4.3.2. Ukraine

Figure 12 (left) shows the oblast level validation of the model. The statistics show very similar numbers to the subnational validation in the US, with a little lower coefficient of determination of 0.73 but a very similar error (to the US) of 0.48 t/ha (14%). Note that in Ukraine there are much fewer oblasts than there are counties in the US and they present a lower variability in yield as compared to the US (the yield at oblast level can get up to 6 t/ha, while in the US counties it can get up to 10 t/ha). Validation of the model using data from the year of 2017 (the only year not used in the calibration) shows a general tendency of underestimation, but performance statistics are similar with the ones derived for 2001 to 2016 validation. Figure 12 (right) shows the oblast level LOOCV. Figure 13 shows the national validation. The results show an excellent coefficient of determination of 0.85,
being able to capture the very low yields in 2003. As in the case of the oblast validation, 594 2017 is underestimated by the model. The RMSE in Ukraine at the national level is 0.27 595 t/ha (8%). 596

597

<figure><img src="image_21.png" /></figure>

598

Figure 12. *Oblast* level calibration validation in Ukraine (left). Each black dot represents 599 an *oblast* and a year from 2001 to 2016 while each red dot represents an *oblast* in 2017. 600
601 LOOCV at oblast level in Ukraine (right).
<figure><img src="image_22.png" /><figcaption>Figure 13. National level validation in Ukraine. 2017 is the black triangle.</figcaption></figure>

## 5. DISCUSSION

Most crop yield models based on remote sensing coarse resolution data use static cropland or crop type masks to aggregate spatially the remotely sensed variables as discussed in the introduction section. When working with static crop type masks, coarsening the spatial resolution of the EO time series data can result in a constant per-pixel proportion over multiple seasons [54]. Besides, cropland static masks are more stable, but do not isolate different crops. The studies working on these masks generally look at a special time of the growing season, when the crop considered shows a particular pattern of the variable studied that allows its isolation from the other crops in the
cropland mask. However, for summer crops, whose phenology is very similar the mix of different crops can lead to wrong conclusions. In this work, we have developed a model to un-mix the wheat signal from other surfaces surrounding the fields by using yearly crop type masks. This method assumes that the DVI<sub>wheat</sub> and the DVI<sub>others</sub> remain constant for an AU for a given date. Therefore, we assume that within the AU the climate is similar, there are no significant elevation differences and the vegetation types in the area are in the same phenological stage. This is a big assumption specially for AUs with very different varieties of wheat and treatments. However, the model works well at sub-national to national scale. Every model has a trade-off and we acknowledge this limitation of coarse spatial resolution data. This will be addressed by using moderate to high spatial resolution data, that will allow the yield estimation at field level. However, the trade-off of such sensors (i.e. Landsat, that dates back to the 70s) is the low temporal resolution with 16 days revisiting time during the historical time series. With the availability of Sentinel 2 data from 2015, its combination with Landsat [55] will address this gap, but by now we can use coarse resolution time series to train the model and learn from it. Looking at the 100% wheat pixels average at the Landsat level for the three different counties considered, the maximum DVI standard deviation is 0.06 which means an error of 10-15% at the peak. Nevertheless, when evaluating the method with Landsat data, the results show good agreement compared to the county 30m wheat pixels average with good correlation and errors of 0.02. One limitation of the method, related to the limitation of working with coarse to moderate spatial resolution data, is that it can only be applied when there is at least one pixel with wheat purity higher than 40% within the AU.
As discussed above, by establishing this threshold, we ensure that there is enough of a wheat signal to be un-mixed from the other surfaces.

Our model is based on the study of different patterns of the wheat DVI evolution through the growing season. In the introduction, we described some yield models that use other vegetation indexes such as the NDVI as a reference. However, the MODIS Surface Reflectance product is now well established [56] and widely used by the science community as a reference. In this context and given the limitations of using the NDVI on dense vegetation (e.g., [57]–[63]), much simpler vegetation indices that were originally discarded as they did not minimize for atmospheric or BRDF effects, can now be used. This is the case of the DVI that we use in this work. The good results obtained using this index demonstrate that it is well correlated to crop yields and is responsive to high yield values. As an example of the performance of this index compared to the NDVI, Figure 14 shows the subnational (top) and national (bottom) validation just using the peak value of the DVI (left) compared to NDVI (right). At subnational level (Figure 14 top), the difference during the calibration period (2001-2016) between using DVI and NDVI is not significant given the high amount of data included in the statistics. However, the validation of a single year (2017) using the peak DVI reduces the RMSE 0.1 t/ha (3%). The difference is also noticed at national level (Figure 14 bottom), where the peak DVI shows better results than the peak NDVI. Note that in this work we choose the DVI for simplicity, but other indexes such as the Enhanced Vegetation Index (EVI) or the Leaf Area Index (LAI) may provide similar results.
<figure><img src="image_23.png" /><figcaption>Figure 14. US county level (top) and national level (bottom) validation when using just the peak DVI (left) and peak NDVI (right).</figcaption></figure>

Besides the un-mixing method and the use of the DVI, the main hypothesis of the model is the correlation of the yield with the DVI seasonal peak amplitude, the time that the peak remains and the average EF 30 days after the peak. The three regressors proposed
are focused on the reproduction stage of the wheat (DVI amplitude and length) and the grain filling process (EF 30-days average). In the introduction, we referenced some EO-based yield models whose main input variables might be the NDVI peak, its temporal integration, its senescence rate, combining the NDVI with the LST or using other indexes such as FAPAR or ESI. While the amplitude of the peak has been used in other studies (but using NDVI instead), the length of the peak has not been analyzed previously (to our best knowledge). Though the peak definition is a singular point, in some cases, the peak can last a few days becoming a short plateau. During this time period, the plant remains in the same phenological stage. Looking at the yield statistics we noticed that the yield is positively correlated with the length of the peak. This is especially the case of 2016 (highest production year) in some counties of Kansas, where the amplitude is very similar to other years in the dataset while the length made a difference. The 30-days average EF, as mentioned in the methods section, accounts for any water stress condition that happen during and after the peak. Regarding the statistical significance of each regressor, the amplitude of the DVI seasonal peak shows the lowest *p*-values of the three regressors both in the US and Ukraine (Figures 5 and 6). Meanwhile, the other two regressors remain significant in some regions of the US as discussed in the results section. However, in Ukraine the length and the EF show less statistical significance throughout the different *oblasts*. Besides, Table 4 shows the number of US counties where each calibration equation is used. The most used equations are 11 (just using the amplitude) and 14 (amplitude and length) that account for over 70% of the counties. However, the equations that include the EF (13, 15 and 16) represent about the 25%. It is also interesting that the three regressors are just considered in 14.2% of counties.
686

687 Table 4. Number of US counties where each equation is used.

<table><thead><tr><td>Equation</td><td>Number of counties</td><td>%</td></tr></thead><tbody><tr><td>11</td><td>167</td><td>38.8</td></tr><tr><td>12</td><td>27</td><td>6.3</td></tr><tr><td>13</td><td>10</td><td>2.3</td></tr><tr><td>14</td><td>147</td><td>34.2</td></tr><tr><td>15</td><td>18</td><td>4.2</td></tr><tr><td>16</td><td>61</td><td>14.2</td></tr></tbody></table>

688

689

690 When applying the method to the US, the proposed model was able to capture most of the extreme events that occurred during the MODIS era. Figure 15 shows the county level validation of two example of extreme events in the US. In 2012, the combination of hot weather and dry conditions during May affected the final yield, especially in the Great Plains (Figure 15 left). Besides, during mid-May 2015 Texas and Oklahoma were hit by a flood event that damaged the wheat fields (Figure 15 right).

696

697

698

699
<figure><img src="image_24.png" /><figcaption>Figure 15. Examples of the model performance in areas hit by drought (left) and flooding (right).</figcaption></figure>

The model was also able to capture high yields, such as 2016, when the rains in late April and mild temperatures in May led to the highest historical yields in the US. However, the model is currently not able to capture late frosts such as the one that impacted Kansas and Oklahoma during early April of 2007. However, it is able to capture the 2003 February-March frost that killed over 50% of the wheat crops in Ukraine. Winter wheat plants are tolerant to temperatures of -20 °C in the vegetative stages, but suffer severe damage at much more moderate temperatures (-5 to -7 °C) during the reproductive stages (Fuller et al., 2007). Thus, the impact of a frost depends on the development stage of the wheat [65]. In Ukraine, the frost occurred in earlier phenological stages of the plant and it affected the biomass by showing much lower peak amplitudes than other “normal” years. Additionally, the failure was result not just of the frost but also of winter kill, low snow cover and a May-June drought. However, the frost in the US happened later in the phenological stage of the plant and did not have any consequence on the amplitude or the
716 length of the peak. Given that the 2007 late frost is the only such occurrence during the
717 MODIS era, we did not have enough data to include a parameter responsive to this event
718 to train the model. In the future, we will explore the application of the model to AVHRR
719 historical data looking for other late frost events happening in the past in the US.

720 The negative bias of the model estimations at the US national level means that the
721 counties not included in the study have generally higher yields than the counties
722 considered. Note that low purity means small field size. By dividing by the slope we
723 assume that the contribution of these counties remains constant through the time series.

724 This assumption is acceptable given that the regression coefficient remains high (0.82
725 compared to 0.88). In Ukraine, the subnational statistics are divided into fewer cases than
726 the US (29 *oblasts* versus about 500 counties). The main advantage of having coarser
727 sub-divisions is that we can apply the model to every *oblast*, that there are statistics for
728 all of them and the *oblast*-level performance is comparable to the county-level (same
729 error of 0.4 t/ha). However, this lower number of subnational statistics leads to a higher
730 error at the national level as compared to the US.

731 Comparing the results of the proposed model to other methods in the literature to monitor
732 wheat in the US and Ukraine, the proposed model improves [4] and [15] errors both in
733 the US (7.5% to 5.5%) and Ukraine (12% to 8%) and increases significantly the
734 coefficient of determination in the US (0.24 to 0.77) and Ukraine (0.71 to 0.88).
735 However, such studies just analyzed data from 2001 to 2012. Compared to [20], Ukraine
736 results at oblast level (0.4 to 0.6 t/ha) are in line with their reported errors (0.5 to 0.8 t/ha)
737 though their study just included the time series from 2001 to 2012.
738 Going forward we will start to run this model in real time to forecast in season yields.
739 This will be based on the processing of in season wheat masks using [26] and the
740 integration of [15] model to forecast the peak amplitude. Regarding the length, we can
741 explore the DVI versus the GDD evolution and adjust it to a function to forecast the
742 length. Finally, the 30-days average EF can be related with different climatology
743 variables that can be forecast by weather models. Additionally, as new data becomes
744 available, it will continue to then improve the model.

745

746 6. CONCLUSIONS

747 This work presents a new EO-based empirical winter wheat yield model. It is based on
748 the un-mixing of the wheat signal by using yearly crop type masks. This wheat signal is
749 used to calibrate the model using as inputs the seasonal amplitude and length of the DVI
750 peak and the average of the EF 30 days after the peak. The model is calibrated in the US
751 at the county level and in Ukraine at the oblast level using historical statistics from 2001
752 to 2016. In each administrative unit, a different calibration coefficient (based on all
753 possible combination of the three regressors) is selected depending on the statistical
754 significance of each variable. Finally, the model is validated not just during 2001 to 2016
755 but also including 2017, which is not used in the calibration process. Overall, the
756 presented model shows a strong performance both at the national (R² higher than 0.7 and
757 RMSE of 0.11 t/ha in the US and 0.27 t/ha in Ukraine) and the subnational level (R²
758 higher than 0.7 and a RMSE lower than 0.6 t/ha). Our goal in the future is to extend the
759 model to forecast the yield and apply it to other crops and other countries.

760
761 ACKNOWLEDGEMENTS

762 This research was funded within the NASA Land-Cover/Land-Use Change (LCLUC)

763 Program, Grant number 80NSSC18K0336, and Cooperative Agreement

764 #80NSSC18M0039 (NASA Harvest).

765

766 REFERENCES

767 [1] U. FAO, "How to feed the world in 2050," in *Rome: High-Level Expert Forum*,

768 2009.

769 [2] C. Rosenzweig, A. Iglesias, X. B. Yang, P. R. Epstein, and E. Chivian, "Climate

770 Change and Extreme Weather Events; Implications for Food Production, Plant

771 Diseases, and Pests," *Glob. Chang. Hum. Heal.*, vol. 2, no. 2, pp. 90-104, 2001.

772 [3] R. B. Macdonald, "A summary of the history of the development of automated

773 remote sensing for agricultural applications," *IEEE Transactions on Geoscience*

774 and *Remote Sensing*, vol. GE-22, no. 6. pp. 473-482, 1984.

775 [4] I. Becker-Reshef, E. Vermote, M. Lindeman, and C. Justice, "A generalized

776 regression-based model for forecasting winter wheat yields in Kansas and Ukraine

777 using MODIS data," *Remote Sens. Environ.*, vol. 114, no. 6, pp. 1312-1323, 2010.

778 [5] P. J. Pinter, R. D. Jackson, S. B. Idso, and R. J. Reginato, "Multidate spectral

779 reflectance as predictors of yield in water stressed wheat and barley," *Int. J.*

780 *Remote Sens.*, vol. 2, no. 1, pp. 43-48, Jan. 1981.

781 [6] C. J. Tucker, C. Vanpraet, E. Boerwinkel, and A. Gaston, "Satellite remote sensing

782 of total dry matter production in the Senegalese Sahel," *Remote Sens. Environ.*,

783 vol. 13, no. 6, pp. 461-474, 1983.
784 [7] S. N. Goward and D. G. Dye, "Evaluating North American net primary productivity with satellite observations," *Adv. Sp. Res.*, vol. 7, no. 11, pp. 165-174, 1987.

787 [8] P. C. Doraiswamy, T. R. Sinclair, S. Hollinger, B. Akhmedov, A. Stern, and J. Prueger, "Application of MODIS derived parameters for regional crop yield assessment," *Remote Sens. Environ.*, vol. 97, no. 2, pp. 192-202, 2005.

790 [9] M. Meroni, E. Marinho, N. Sghaier, M. M. Verstrate, and O. Leo, "Remote Sensing Based Yield Estimation in a Stochastic Framework — Case Study of Durum Wheat in Tunisia," *Remote Sensing*, vol. 5, no. 2. 2013.

793 [10] F. Baret and G. Guyot, "Monitoring of the ripening period of wheat canopies using visible and near infra red radiometry [reflectance, vegetation index, senescence rate, water plateau]," *Agronomie (France)*, vol. v. 6. 1986.

796 [11] D. M. Johnson, "An assessment of pre- and within-season remotely sensed variables for forecasting corn and soybean yields in the United States," *Remote Sens. Environ.*, vol. 141, pp. 116-128, 2014.

799 [12] J. H. Kastens, T. L. Kastens, D. L. A. Kastens, K. P. Price, E. A. Martinko, and R.-Y. Lee, "Image masking for crop yield forecasting using AVHRR NDVI time series imagery," *Remote Sens. Environ.*, vol. 99, no. 3, pp. 341-356, 2005.

802 [13] D. M. Johnson and R. Mueller, "The 2009 Cropland Data Layer.," *Photogramm. Eng. Remote Sens.*, vol. 76, no. 11, pp. 1201-1205, 2010.

804 [14] T. Fisette et al., "AAFC annual crop inventory," *2013 Second International Conference on Agro-Geoinformatics (Agro-Geoinformatics)*. pp. 270-274, 2013.

806 [15] B. Franch et al., "Improving the timeliness of winter wheat production forecast in
807 the United States of America, Ukraine and China using MODIS data and NCAR
808 Growing Degree Day information," *Remote Sens. Environ.*, vol. 161, pp. 131–148,
809 2015.

810 [16] B. Franch *et al.*, "A 30+ Year AVHRR Land Surface Reflectance Climate Data
811 Record and Its Application to Wheat Yield Monitoring," *Remote Sensing*, vol. 9,
812 no. 3. 2017.

813 [17] W. Kowalik, K. Dabrowska-Zielinska, M. Meroni, T. U. Raczka, and A. de Wit,
814 "Yield estimation using SPOT-VEGETATION products: A case study of wheat in
815 European countries," *Int. J. Appl. Earth Obs. Geoinf.*, vol. 32, pp. 228–239, 2014.

816 [18] G. Büttner, J. Feranec, G. Jaffrain, L. Mari, G. Maucha, and T. Soukup, "The
817 CORINE land cover 2000 project," *EARSEL eProceedings*, vol. 3, no. 3, pp. 331–
818 346, 2004.

819 [19] R. López-Lozano *et al.*, "Towards regional grain yield forecasting with 1km-
820 resolution EO biophysical products: Strengths and limitations at pan-European
821 level," *Agric. For. Meteorol.*, vol. 206, pp. 12–32, 2015.

822 [20] F. Kogan *et al.*, "Winter wheat yield forecasting in Ukraine based on Earth
823 observation, meteorological data and biophysical models," *Int. J. Appl. Earth Obs.
824 Geoinf.*, vol. 23, pp. 192–203, 2013.

825 [21] M. C. Anderson *et al.*, "The Evaporative Stress Index as an indicator of
826 agricultural drought in Brazil: An assessment based on crop yield impacts,"
827 *Remote Sens. Environ.*, vol. 174, pp. 82–99, 2016.

828 [22] A. J. Richardson and C. L. Wiegand, "Distinguishing vegetation from soil
829 background information.," *Photogramm. Eng. Remote Sensing*, vol. 43, pp. 1541–
830 1552, 1977.

831 [23] G. Vocke and M. B. Ali, *US wheat production practices, costs, and yields: Variations across regions*. United States Department of Agriculture, Economic Research Service, 2013.

834 [24] C. Boryan, Z. Yang, R. Mueller, and M. Craig, "Monitoring US agriculture: the US Department of Agriculture, National Agricultural Statistics Service, Cropland Data Layer Program," *Geocarto Int.*, vol. 26, no. 5, pp. 341-358, Aug. 2011.

837 [25] E. Vermote, "MOD09A1 MODIS/Terra Surface Reflectance 8-Day L3 Global 500m SIN Grid V006," *NASA EOSDIS L. Process. DAAC*, vol. 10, 2015.

839 [26] S. Skakun et al., "Early season large-area winter crop mapping using MODIS NDVI data, growing degree days information and a Gaussian mixture model," *Remote Sens. Environ.*, vol. 195, pp. 244-258, 2017.

842 [27] M. Lavreniuk, N. Kussul, S. Skakun, A. Shelestov, and B. Yailymov, "Regional retrospective high resolution land cover for Ukraine: Methodology and results," in *2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS)*, 2015, pp. 3965-3968.

846 [28] M. S. Lavreniuk et al., "Large-Scale Classification of Land Cover Using Retrospective Satellite Data," *Cybern. Syst. Anal.*, vol. 52, no. 1, pp. 127-138, 2016.

849 [29] F. M. Bréon, E. Vermote, E. F. Murphy, and B. Franch, "Measuring the Directional Variations of Land Surface Reflectance From MODIS," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 53, no. 8. pp. 4638-4649, 2015.
853 [30] E. Vermote, C. O. Justice, and F.-M. Bréon, “Towards a generalized approach for correction of the BRDF effect in MODIS directional reflectances,” *IEEE Trans. Geosci. Remote Sens.*, vol. 47, no. 3, pp. 898–908, 2009.

856 [31] B. Franch, E. F. Vermote, J. a. Sobrino, and Y. Julien, “Retrieval of Surface Albedo on a Daily Basis: Application to MODIS Data,” *IEEE Trans. Geosci. Remote Sens.*, vol. 52, no. 12, pp. 7549–7558, 2014.

859 [32] C. B. Schaaf et al., “First operational BRDF, albedo nadir reflectance products from MODIS,” *Remote Sens. Environ.*, vol. 83, no. 1–2, pp. 135–148, 2002.

861 [33] J.-L. Roujean, M. Leroy, and P.-Y. Deschamps, “A bidirectional reflectance model of the Earth’s surface for the correction of remote sensing data,” *J. Geophys. Res.*, vol. 97, no. D18, p. 20455, 1992.

864 [34] F. Maignan, F. M. Bréon, and R. Lacaze, “Bidirectional reflectance of Earth targets: Evaluation of analytical models using a large set of spaceborne measurements with emphasis on the Hot Spot,” *Remote Sens. Environ.*, vol. 90, no. 2, pp. 210–220, 2004.

868 [35] X. Li and A. H. Strahler, “Geometric-optical bidirectional reflectance modeling of the discrete crown vegetation canopy: effect of crown shape and mutual shadowing,” *IEEE Trans. Geosci. Remote Sens.*, vol. 30, no. 2, pp. 276–292, 1992.

871 [36] W. Lucht, “Expected retrieval accuracies of bidirectional reflectance and albedo from EOS-MODIS and MISR angular sampling,” *J. Geophys. Res.*, vol. 103, no. D8, pp. 8763–8778, 1998.

874 [37] P. Lewis and M. Barnsley, “Influence of the sky radiance distribution on various formulations of the earth surface albedo,” *Proceedings of the Conference on*
876 *Physical Measurements and Signatures in Remote Sensing*. pp. 707–715, 1994.

877 [38] W. Lucht and P. Lewis, “Theoretical noise sensitivity of BRDF and albedo retrieval from the EOS-MODIS and MISR sensors with respect to angular sampling,” *Int. J. Remote Sens.*, vol. 21, no. 1, pp. 81–98, 2000.

880 [39] E. Vermote, D. Tanre, J. L. Deuze, M. Herman, and J. J. Morcrette, “Second Simulation of the Satellite Signal in the Solar Spectrum (6S). 6S User Guide Version 2. Appendix III: Description of the subroutines,” vol. 35, no. 3, 1997.

883 [40] W. G. M. Bastiaanssen, M. Menenti, R. A. Feddes, and A. A. M. Holtslag, “A remote sensing surface energy balance algorithm for land (SEBAL). 1. Formulation,” *J. Hydrol.*, vol. 212, pp. 198–212, 1998.

886 [41] R. Allen, M. Tasumi, and R. Trezza, “Satellite-Based Energy Balance for Mapping Evapotranspiration with Internalized Calibration (METRIC)—Model,” *J. Irrig. Drain. Eng.*, vol. 133, no. 4, pp. 380–394, Aug. 2007.

889 [42] G. Ravazzani, C. Corbari, S. Morella, P. Gianoli, and M. Mancini, “Modified Hargreaves-Samani Equation for the Assessment of Reference Evapotranspiration in Alpine River Basins,” *J. Irrig. Drain. Eng.*, vol. 138, no. 7, pp. 592–599, Jul. 2012.

893 [43] M. Menenti and B. J. Choudhury, “Parameterization of land surface evaporation by means of location dependent potential evaporation and surface temperature range,” *Proc. IAHS Conf. L. Surf. Process.*, no. 212, pp. 561–568, 1993.

896 [44] G. J. Roerink, Z. Su, and M. Menenti, “S-SEBI: A simple remote sensing algorithm to estimate the surface energy balance,” *Phys. Chem. Earth, Part B Hydrol. Ocean. Atmos.*, vol. 25, no. 2, pp. 147–157, 2000.
899 [45] Z. Su, "The Surface Energy Balance System (SEBS) for estimation of turbulent
900 heat fluxes," *Hydrol. Earth Syst. Sci. Discuss.*, vol. 6, no. 1, pp. 85–100, 2002.

901 [46] M. C. Anderson *et al.*, "Evaluation of Drought Indices Based on Thermal Remote
902 Sensing of Evapotranspiration over the Continental United States," *J. Clim.*, vol.
903 24, no. 8, pp. 2025–2044, Apr. 2011.

904 [47] M. Menenti, W. Bastiaanssen, D. van Eick, and M. A. Abd el Karim, "Linear
905 relationships between surface reflectance and temperature and their application to
906 map actual evaporation of groundwater," *Adv. Sp. Res.*, vol. 9, no. 1, pp. 165–176,
907 1989.

908 [48] M. Gómez, A. Olioso, J. A. Sobrino, and F. Jacob, "Retrieval of
909 evapotranspiration over the Alpilles/ReSeDA experimental site using airborne
910 POLDER sensor and a thermal camera," *Remote Sens. Environ.*, vol. 96, no. 3, pp.
911 399–408, 2005.

912 [49] Z. Su, H. Pelgrum, and M. Menenti, "Aggregation effects of surface heterogeneity
913 in land surface processes," *Hydrol. Earth Syst. Sci.*, vol. 3, no. 4, pp. 549–563,
914 Dec. 1999.

915 [50] W. W. Verstraeten, F. Veroustraete, and J. Feyen, "Estimating evapotranspiration
916 of European forests from NOAA-imagery at satellite overpass time: Towards an
917 operational processing chain for integrated optical and thermal sensor data
918 products," *Remote Sens. Environ.*, vol. 96, no. 2, pp. 256–276, 2005.

919 [51] Z. Wan and J. Dozier, "A generalized split-window algorithm for retrieving land-
920 surface temperature from space," *IEEE Trans. Geosci. Remote Sens.*, vol. 34, no.
921 4, pp. 892–905, Jul. 1996.
922 [52] R. A. Fischer and R. Maurer, "Drought resistance in spring wheat cultivars. I. 923 Grain yield responses," *Aust. J. Agric. Res.*, vol. 29, no. 5, pp. 897-912, 1978.

924 [53] F. Giunta, R. Motzo, and M. Deidda, "Effect of drought on yield and yield 925 components of durum wheat and triticale in a Mediterranean environment," *F. 926 Crop. Res.*, vol. 33, no. 4, pp. 399-409, 1993.

927 [54] I. Becker-Reshef *et al.*, "Prior Season Crop Type Masks for Winter Wheat Yield 928 Forecasting: A US Case Study," *Remote Sensing*, vol. 10, no. 10. 2018.

929 [55] M. Claverie *et al.*, "The Harmonized Landsat and Sentinel-2 surface reflectance 930 data set," *Remote Sens. Environ.*, vol. 219, pp. 145-161, 2018.

931 [56] E. F. Vermote, N. Z. El Saleous, and C. O. Justice, "Atmospheric correction of 932 MODIS data in the visible to middle infrared: first results," *Remote Sens. Environ.*, 933 vol. 83, no. 1-2, pp. 97-111, 2002.

934 [57] A. A. Gitelson, "Wide Dynamic Range Vegetation Index for Remote 935 Quantification of Biophysical Characteristics of Vegetation," *J. Plant Physiol.*, 936 vol. 161, no. 2, pp. 165-173, 2004.

937 [58] F. Baret and G. Guyot, "Potentials and limits of vegetation indices for LAI and 938 APAR assessment," *Remote Sens. Environ.*, vol. 35, no. 2, pp. 161-173, 1991.

939 [59] A. A. Gitelson, Y. J. Kaufman, R. Stark, and D. Rundquist, "Novel algorithms for 940 remote estimation of vegetation fraction," *Remote Sens. Environ.*, vol. 80, no. 1, 941 pp. 76-87, 2002.

942 [60] A. A. Gitelson, R. Stark, U. Grits, D. Rundquist, Y. Kaufman, and D. Derry, 943 "Vegetation and soil lines in visible spectral space: A concept and technique for 944 remote estimation of vegetation fraction," *Int. J. Remote Sens.*, vol. 23, no. 13, pp.
945 2537–2562, Jan. 2002.

946 [61] C. Buschmann and E. Nagel, “In vivo spectroscopy and internal optics of leaves as
basis for remote sensing of vegetation,” *Int. J. Remote Sens.*, vol. 14, no. 4, pp.

948 711–722, Mar. 1993.

949 [62] P. J. Sellers *et al.*, “Remote sensing of the land surface for studies of global
change: Models — algorithms — experiments,” *Remote Sens. Environ.*, vol. 51,

951 no. 1, pp. 3–26, 1995.

952 [63] S. Skakun, E. Vermote, J.-C. Roger, and B. Franch, “Combined Use of Landsat-8
and Sentinel-2A Images for Winter Crop Mapping and Winter Wheat Yield

954 Assessment at Regional Scale,” *AIMS Geosci.*, vol. 3, no. 2, pp. 163–186, 2017.

955 [64] M. P. Fuller, A. M. Fuller, S. Kaniouras, J. Christophers, and T. Fredericks, “The
freezing characteristics of wheat at ear emergence,” *Eur. J. Agron.*, vol. 26, no. 4,

957 pp. 435–441, 2007.

958 [65] J. P. Shroyer, M. E. Mikesell, and G. M. Paulsen, *Spring freeze injury to Kansas
wheat*. Cooperative Extension Service, Kansas State University, 1995.

960 Vermote, E., 2015. MOD09Q1 MODIS/Terra Surface Reflectance 8-day L3 Global 250m
SIN Grid V006. NASA EOSDIS Land Processes DAAC.

962 http://dx.doi.org/10.5067/MODIS/MOD09Q1.006 (Accessed 26.02.2018).

963 URL 1. http://www.amis-outlook.org/

964 URL 2. http://www.worldstopexports.com/wheat-exports-country/