Regional-Scale Cotton Yield Estimate in the Brazilian Cerrado Using Simple Linear Regression Models with Time-Series of Vegetation Indices from EOS-MODIS Satellite Data

1. Introduction

Cotton is an important cash crop in Brazil, which is the world's fourth largest producer and the second largest exporter, with 1.66 million hectares cultivated in the 2022-2023 season and a total lint production of 3.03 million tons [1]. Most of the cotton has been cultivated in the Brazilian Cerrado, a biome typically characterized by flat and arable lands with good rainfall (750 to 2000 mm/year for the whole Cerrado area [2]), under a highly mechanized, rainfed (92% of the total production [3]), intensive double crop production system, generally cultivated from January to July after soybean [4]. The main cotton producing states have been Mato Grosso and Bahia, which together account for around 90% of the total cotton produced in Brazil [5].

Estimating in-season cotton yield over large areas at a regional or national level provides critical information for farmers, policymakers, governments, crop insurers and commodity traders. Yield estimations are generally performed by a combination of field surveys, remote sensing, statistical regression and crop simulation models [6,7]. Particularly for farmers, the mid-season crop forecasting is very important for management, harvesting, storage and transport logistics decisions and also for planning the next crop [8,9]. Forecasting cotton yield in advance (i.e. before 90-100 days after sowing) is a valuable management information for making viable corrective interventions, such as adjustments to the application of planned inputs (top dressing fertilizer, growth regulators and others), comparing predicted to the expected or historical yields of specific fields.

Although crop yield models have shown adequate performances in field-scale estimations, their applications in large areas are difficult due to the huge amount of data and computational processing costs [10]. To overcome this limitation, statistical regression models such as multiple linear regression and machine learning techniques have been applied to combine spectral indices obtained from satellite data and climate variables, at regional or national level [10-12]. In Table 1 is shown several cotton yield estimation studies using statistical regression models, or crop models combined with remote sensing at regional and field scale under different platforms (satellite, airplane, unmanned aerial vehicle (UAV) or ground sensors) in different countries. Three of these studies developed regression models (Random Forest, among other techniques) to predict cotton yield at a regional scale in China [10], India [11] and Australia [12] from satellite-derived vegetation indices (VIs) and climate variables as covariates. They included precipitation, temperature, evapotranspiration, vapor pressure, soil moisture, among other climate variables. Hundreds of data instances from thousands of hectares in different farm fields over many crop years were used to train and test model performances, which resulted in root mean square errors (RMSE) for seed cotton yield of 157 kg ha^-1 [10], 375 kg ha^-1 [11] and 976 kg ha^-1 [12], in China, India and Australia, respectively.

Simple linear regression models have been developed to predict cotton yield using only satellite-derived vegetation indices as independent variables [13-16], unmanned aerial vehicle (UAV) [17-20], airplane [21] and ground sensor [22]. If on the one hand VIs by themselves have shown some limitations for describing the complex relationships among production, plant physiology, climate, soil nutrition, soil water parameters, pest and disease infestation, crop design and management characteristics [10], on the other several studies have indicated that crop yield can be predicted with relatively good accuracies with this approach [7,23,24]. Just to give an idea of the average magnitude assessed from the summary review shown in Table 1, the average of RMSEs for seed cotton yield for the linear regressions was close to that one assessed by using machine learning techniques and crop model approaches, which indicates linear regression models as potentially advantageous since they operate with a smaller set of well-defined (independent) variables, fewer data manipulation and smaller amounts of instances compared to the other approaches. Moreover, the greater mathematical explainability inherent in linear equations for specific satellite products, VI and crop phenological stage can be directly transferred and/or easily adapted to other locations, allowing clearer understandings and comparisons among different models.

Table 1. A review of previous studies for estimation of seed cotton yield using remote sensing.

Reference	Country	Plot	Area	Model	RS source$	Regression	RMSE##	R2
			ha	Aproach&		Model#	kg ha-1
[10]	China	355	-	CV/RS	Modis/Sentinel	LSTM, SVM, RF	375	0.65
[42]	EUA	12	150	CM/RS	Spectroradiometer		468	-
[13]	USA	3	188	RS	Modis/Landsat	LR	673	0.52
[14]	USA	-	-	RS	Modis	LR	-	0.16
[11]	India	-	-	CV, RS	Modis	RF	157	0.69
[12]	Australia	253	-	CV, RS	Landsat	RF	976	-
[17]	USA	1	5	RS	UAV	MLR	261	0.87
[36]	USA	805	0,65	RS	UAV	ANN, RF	-	0.72
[39]	USA	1	57	RS	Modis/Landsat		463	0.84
[43]	USA	-	-	EM/RS	Sentinel		-	-
[18]	USA	2550	6	RS	UAV	LR	550	0.92
[37]	Brazil	1	90	RS	Optical sensor	decision trees	-	0.81
[40]	USA	-	73	RS	Landsat	ANN	375-470	0.71
[38]	USA	2	120	RS	Landsat	exponential	481	0.81
[19]	Australia	90	7	RS	UAV	LR and quadradic	-	0.75
[15]	USA	949	-	RS	Modis	LR	-	0.48
[16]	USA	24	0,2	RS	NASA data	LR	-	0.85
[21]	USA	48	1,5	RS	Airborne	LR	-	0.47
[22]	USA	44	5,3	RS	Spectroradiometer	LR	-	0.89
&CV: Climate Variables, RS: Remote Sensing, CM: Crop Model, EM: Ecosystem Model; %RS: Remote Sensing;
UAV: Unmanned Aerial Vehicle; #LSTM: Long Short-Term Memory; SVM: Support Vector Machine; RF: Random Forest;
LR: Linear Regression, MLR: Multiple Linear Regression, ANN: Artificial Neural, Network;
##Root Mean Square Error of seed cotton yield.

Table 1. A review of previous studies for estimation of seed cotton yield using remote sensing.

Reference	Country	Plot	Area	Model	RS source$	Regression	RMSE##	R2
			ha	Aproach&		Model#	kg ha-1
[10]	China	355	-	CV/RS	Modis/Sentinel	LSTM, SVM, RF	375	0.65
[42]	EUA	12	150	CM/RS	Spectroradiometer		468	-
[13]	USA	3	188	RS	Modis/Landsat	LR	673	0.52
[14]	USA	-	-	RS	Modis	LR	-	0.16
[11]	India	-	-	CV, RS	Modis	RF	157	0.69
[12]	Australia	253	-	CV, RS	Landsat	RF	976	-
[17]	USA	1	5	RS	UAV	MLR	261	0.87
[36]	USA	805	0,65	RS	UAV	ANN, RF	-	0.72
[39]	USA	1	57	RS	Modis/Landsat		463	0.84
[43]	USA	-	-	EM/RS	Sentinel		-	-
[18]	USA	2550	6	RS	UAV	LR	550	0.92
[37]	Brazil	1	90	RS	Optical sensor	decision trees	-	0.81
[40]	USA	-	73	RS	Landsat	ANN	375-470	0.71
[38]	USA	2	120	RS	Landsat	exponential	481	0.81
[19]	Australia	90	7	RS	UAV	LR and quadradic	-	0.75
[15]	USA	949	-	RS	Modis	LR	-	0.48
[16]	USA	24	0,2	RS	NASA data	LR	-	0.85
[21]	USA	48	1,5	RS	Airborne	LR	-	0.47
[22]	USA	44	5,3	RS	Spectroradiometer	LR	-	0.89
&CV: Climate Variables, RS: Remote Sensing, CM: Crop Model, EM: Ecosystem Model; %RS: Remote Sensing;
UAV: Unmanned Aerial Vehicle; #LSTM: Long Short-Term Memory; SVM: Support Vector Machine; RF: Random Forest;
LR: Linear Regression, MLR: Multiple Linear Regression, ANN: Artificial Neural, Network;
##Root Mean Square Error of seed cotton yield.

The present study evaluates the potential of using simple linear models for correlating yield and satellite-derived VIs from EOS-MODIS (Moderate Resolution Imaging Spectroradiometer) to forecast in-season cotton yield in the Brazilian Cerrado. Reflectance time-series from several MODIS spectral bands were extracted for assessing 15-day interval averages (from sowing to harvest) of well-known VIs. These regular interval averages were evaluated as predictors (independent variables) in linear regression models for estimating cotton yield. The on-farm research was carried out in 281 commercial farm fields under different cotton yields, sowing dates and geographically sparse (with enough geographic information for delimitating the cotton farm field contour).

2. Materials and Methods

2.1. Experimental areas and dataset

The regression models were trained and tested using data from 281 commercial cotton farm fields (hereinafter referred to as plots) in the states of Mato Grosso (125 plots), Goiás (83 plots) and Bahia (73 plots) for the growing seasons of 2016 to 2022. Data from 2016, 2017 and 2018 seasons were obtained from surveys carried out by the Mato Grosso Cotton Institute [25,26] and the Brazilian Agricultural Research Corporation [27], and from cotton farmers in seasons 2019 to 2022. In Figure 1 is depicted the spatial distribution of the 281 plots. Most plots in Mato Grosso, Goiás and Bahia were located at the west, south and west parts of the states, respectively.

The dataset includes information of seed cotton yield for the whole plot (kg ha^-1), plot area (ha), cultivar, plant line spacing (m), seed population (seed ha^-1), sowing date and water management (irrigated or rainfed). In Figure 2 are shown histograms for some of these parameters. The average seed cotton yield assessed for the 281 plots was 4,209 kg ha^-1, ranging from 694 kg ha^-1 to 7,361 kg ha^-1 (Figure 2a). Low yield plots were located in Mato Grosso and Goiás under rainfed production in seasons 2016 and 2022, something particularly caused by drought or irregular rainfall distributions during these crop seasons [28]. Most of the high yield plots was located in Bahia under central pivot irrigation and in Mato Grosso and Goiás during years with good rain distributions (2017 and 2018). Only 20% of the plots (56) were cultivated under irrigation, being all of them in Bahia and Goiás, whereas the remaining (80%) were under rainfed production (225 plots). Plots in Mato Grosso (125) were all under rainfed production. Plant line spacing used was mostly 0.9 m with average seed linear density of 10 seed m^-1, resulting in average population of 111,111 seed ha^-1. The average area of the plots was 188 ha (Figure 2c) and average cotton cycle was 196 days (Figure 2d).

Daily accumulated precipitations, maximum and minimum air temperatures data were obtained for each plot by mean of the Brazilian agrometeorological monitoring system (web application) AGRITEMPO [29], by taking the closest climate station to each plot and downloading the daily climate data from 0 to 180 days after sowing (DAS). The data were then averaged for month 1 (0-30 DAS) to 6 (150-180 DAS) and then correlated to the observed yield only to evaluate the effects of the climate variables on yield and to help explain possible bias between observed and predicted yields.

2.2. Satellite data acquisition and preprocessing

With geographic coordinates of all plots, contour maps (kml format) were generated using the Google Earth Pro (Google, Mountain View, CA), which were then converted in shapefile format files using the QGIS software [30]. The spectral reflectance bands and VIs time-series were obtained from the EOS-MODIS MCD43A4 V6.1 Nadir Bidirectional Reflectance Distribution Function Adjusted Reflectance (NBAR) product, which provides smoothed daily data, with a spatial resolution of 500 meters, corrected based on the 16-day compositions with selection of the best pixels (maximum value composite) in that time period to compose the new interval-mediated time series. The spectral bands in this product are blue (0.459-0.479 m), green (0.545-0.565 m), red (0.62-0.67 m), near infrared-NIR (0.841-0.876 m), short wave infrared-SWIR1 (1.23-1.25 m), SWIR2 (1.628-1.656 m) and SWIR3 (2.105-2.155 m). Google Earth Engine (GEE) platform, by mean of its JavaScript code editor, was programmed to extract the time-series, select the pixels within the plot contour area and generate the average values for each spectral reflectance band for each date from sowing to harvesting. The reflectance image and the time-series graph were visually checked before exporting the representative time-series (comma-separated values - csv format) of each plot. A buffer command of 250 m was applied in the GEE code to avoid pixels located at the plot boundary.

Spectral bands were then combined in an Excel® worksheet and the VIs were calculated, averaged in 15-days intervals from 0 to 240 DAS and correlated to cotton yield for the mutually exclusive intervals. Table 2 lists the VIs evaluated and their definitions. A subset of plots (167) was randomly selected (balanced by yield) and used for generating the regression models (training set) and the remaining dataset (114 plots) was used for the independent model validation (test set), i.e., a well-known cross-validation method set as 60%-40% (train-test). Linear regression models were fitted for each combination of VI (Table 2) and time interval. Model’s accuracies were measured by the Root Mean Square Error (RMSE) and determination coefficient (R²) estimated from the test set.

3. Results and Discussion

Examples of MODIS images (NIR band) and time-series of all spectral reflectance bands are shown in Figure 3 for a plot with very low cotton yield (1,665 kg ha^-1 in season 2016; Figure 3a, left) and high yield (4,965 kg ha^-1 in season 2017; Figure 3a, right). The original smoothed daily time-series are shown in Figure 3b and 15-days averaged time-series in Figure 3c. Reflectance data are plotted in a logarithmic scale for better comparisons among the spectral bands. In general, lower reflectance values (greater absorption of electromagnetic radiation) are observed for the blue, red and green bands, due to the absorption of these light wavelengths by leaves through photosynthesis; higher reflectances in this range can be associated to plant stresses caused by biotic and/or abiotic factors. The spectral range between 0.4 and 0.7 m (visible region) is the photosynthetically active radiation (PAR) that is the most efficient portion of the electromagnetic spectrum used by plant for photosynthesis. Healthy plants present high light absorption in the visible range (0.4-0.7 m) mainly by the leaf pigments like chlorophyll and carotenoids, relatively high reflectance in the NIR range (0.7-1.1 m) due to leaf and cell structure scattering effects, and relatively low reflectance in the SWIR range (1.1-2.5 m) due to water and chemicals into the leaf structure [31,32]. For an actual instance, the low cotton yield (left graph in Figure 3) was mainly caused by low precipitation (water stress) in the 2016 season [28], which is noted by slight increases in red and blue reflectances (lower absorbance by cotton plants), specially between 70 to 170 DAS, and decreasing NIR reflectance, as compared to the high yield plot (right).

In Figure 4a is shown spectral bands time-series averaged for five classes of cotton yield (< 3, 3-3.75, 3.75-4.5, 4.5-5.25, > 5.25 ton ha^-1). Reflectance of NIR and SWIR1 increased with increasing cotton yield, while red, blue, SWIR2, SWIR3 reflectance decreases as yield increased, specially from about 60-80 DAS to harvesting, as can be seen by the correlation coefficients (r) between cotton yield and reflectance shown in Figure 4b. The green band reflectance was higher than red and blue, as expected, but its relationship with yield was not graphically evident (Figure 4a). Best correlations between the spectral reflectance bands and cotton yield were obtained for NIR, SWIR1 (positive), SWIR3, red, SWIR2 and blue (negative) (Figure 4b). These data, obtained using broadband wavelengths satellite sensors are similar to others obtained with narrow band ground measurements by spectroradiometers. For instance, Thenkabail et al. [33] obtained negative correlation between cotton biophysical parameters (leaf area index, wet biomass, plant height and crop yield) and reflectance for wavelengths between 350 nm to about 700 nm (visible light) and positive correlations for wavelengths from about 730 nm to 1050 nm (NIR to part of the SWIR range).

The nine vegetation indices (Table 2) were taken as independent variables to fit simple linear regression models to predict cotton yield using the training dataset (167 plots) for each 15-days interval. Figure 5b shows the determination coefficients (R²) obtained for these VIs and DAS intervals. The highest R² were provided by the TVI (0.69), EVI (0.68), SAVI (0.67) and NDVI (0.61) for intervals between 120 and 165 DAS, which are in the cotton late-season stage of development (open boll to defoliation). Time-series for these VIs averaged for different yield classes are shown in Figure 5a. In general, there was a saturation trend for all VIs at the peak of time-series for the highest yield classes, and this effect was more pronounced for NDVI, agreeing with other studies that show NDVI saturation at high amounts of green biomass or high leaf area index [34,35].

For the training dataset, TVI, EVI and SAVI produced the highest determination coefficients for cotton yield forecasts, indicating a good potential for in-season cotton yield estimations from these VIs. Table 3 presents the linear regression models for the four VIs with the highest R² for different 15-days intervals, from 75-90 to 180-195 DAS, and also for 75-195 DAS (average from 75 DAS to about harvesting) and the maximum VI value along the season (denoted as peak). In Figure 6 is shown the relationships between cotton yield and TVI for different DAS intervals. Although highest R² (0.73) were provided by the interval of 75-195 DAS for the TVI (Figure 6) and also for the other VIs (Table 3), such equations would only be useful at the end of the season, since data needs to be accumulated until 195 DAS before averaging them. For in-season forecasts, the best models were obtained from 105-120 DAS to 165-180 DAS; nevertheless 150-165 and 165-180 DAS intervals are also too close to the end of season and they likely do not have any practical for decision making by farmers. Forecasts based on VI peaks have been used in several studies for in-season yield estimates [14,36]. In the present study the peak models presented R² lower than for 105-120 to 165-180 DAS (Table 3), but higher than 90-105 DAS, which mean good potential for earlier in-season forecasts, since the time-series peak values for TVI, EVI, SAVI and NDVI occurred at about 80 to 100 DAS, as can be seen in Figure 5.

Figure 5. MODIS satellite average time-series of NDVI, EVI, SAVI and TVI for five cotton yield classes (a) and determination coefficients (R²) of linear correlations between cotton yield and 9 vegetation indices (Table 2) for the different 15-days DAS (Days After Sowing) classes, for the training dataset (167 plots).

Figure 5. MODIS satellite average time-series of NDVI, EVI, SAVI and TVI for five cotton yield classes (a) and determination coefficients (R²) of linear correlations between cotton yield and 9 vegetation indices (Table 2) for the different 15-days DAS (Days After Sowing) classes, for the training dataset (167 plots).

Model estimates for cotton yields using TVI equations (Table 3) compared to observed yields with their respective RMSE are shown in Figure 7; Figure 8 displays RMSE (Figure 8a) and R² (Figure 8b) variations with DAS intervals for the four best VIs. Lowest RMSE and highest R² were provided by EVI and TVI for 90-105 to 135-150 DAS (four DAS intervals) allowing in-season yield forecasting with RMSE of about 750 kg ha^-1 (Figure 8a). Averaged VIs for 75-195 DAS and peak models showed similar RMSE values to the best models (EVI, TVI and SAVI) for 15-days intervals (90-105 to 135-150 DAS) (Figure 8c), but for EVI the peak model showed the lowest RMSE (726 kg ha^-1). The average yield predictions for the 114 testing plots are shown in Figure 8d for 105-120 DAS, 75-195 DAS and peaks. The average observed yield (for the 114 plots) was 4,372 kg ha^-1 (dotted horizontal line in Figure 8d) and most of the linear models tend to overestimate yield (about 160 kg ha^-1 in average). In general, there were tendencies of yield overestimation for the lower yield plots and underestimation for the higher yield plots, as can be seen in Figure 7 for TVI, but the 75-195 DAS model seems to reduce such bias. Although NDVI provided the closer average yield forecast, for the 114 plots, to the measured average yield (Figure 8d) its RMSE was much higher than the other VIs (EVI, TVI and SAVI) for all intervals between 75 and 165 DAS. Thus, the best performances were obtained with EVI and TVI for 105-120, 120-35 and 135-150 DAS intervals. If an earlier prediction is required, the 90-105 DAS or peak models are preferred, since RMSE for DAS 75-90 DAS and earlier DAS were too high and R² were inexpressive.

The overestimation trend observed for low yield plots and underestimation for high yields (Figure 7) can be better evaluated correlating the residual error (yield_estimated - yield_observed; estimated by the peak model with TVI) with the observed yield of each plot (Figure 9a), where negative residual errors (underestimation) are verified for yields higher than about 4,000 kg ha^-1 and positive (overestimation) for lower yield values. Several factors may contribute and explain the low and high yield plots, as the climate conditions, soil type, soil fertility, diseases, cultivars and others. Nevertheless, in rainfed cotton production one of the main driven factors that affects yield has been the climatic condition, particularly precipitation and temperature variations along the season. Therefore, the sowing date and cycle duration, established based on the recommended cultivar and local climatic patterns, are impacted by water demands and temperature for each cotton development stages and consequently the expected cotton yield. For the testing dataset, the earliest sowing date was November 19, so observed yield was correlated to sowing after November 1^st as a proxy for sowing data, independent of the year (Figure 9b), showing a negative influence on yield. The cotton cycle duration, which is basically defined by the cultivar used, showed a positive correlation with yield (Figure 9c). These two parameters are intrinsically related to the climatic conditions and both influenced the observed yields. There is a trend of increasing yield for earlier sowing day and for longer cotton cycle duration.

The best correlations between monthly accumulated precipitation and yield occurred in the third and fourth months after sowing (r=0.59 for 60-90 DAS and 0.39 for 90-120 DAS; Figure 9d) that is the cotton mid-season corresponding to the canopy closure until flowering and boll development stages. In the first month, rainfall showed an adverse effect on yield (r=-0.29) since high-intensity rainfall in the very early season can damage the germination process. Monthly averaged maximum and minimum temperatures and their differences (T_max-T_min) also influenced yield mainly after the third month (r=-0.49).

Figure 9. Relationships between the residual error (yieldestimated-yieldobserved; estimated by TVI peak model) and observed cotton yield (a), observed yield and sawing day from November 1st (b), observed yield and cotton cycle (c), and correlation coefficients (r) between observed yield and monthly accumulated precipitation, monthly averaged minimum (Tmin) and maximum (Tmax) temperatures, and Tmax-Tmin (d) for month 1 to 6 after sowing.

Figure 9. Relationships between the residual error (yieldestimated-yieldobserved; estimated by TVI peak model) and observed cotton yield (a), observed yield and sawing day from November 1st (b), observed yield and cotton cycle (c), and correlation coefficients (r) between observed yield and monthly accumulated precipitation, monthly averaged minimum (Tmin) and maximum (Tmax) temperatures, and Tmax-Tmin (d) for month 1 to 6 after sowing.

Several studies as those included in Table 1 evaluated the best period for acquiring Satellite and UAV images or ground measurements for best cotton yield predictions. Some works have used a single date image and others have evaluated images at specific phenological stages or complete time-series data. Table 4 summarizes the best periods, in DAS, inferred from the papers listed in Table 1. The best correlations between cotton yield and the VIs varied from the flowering and boll development period [12,18,22,37,38,39] to boll opening, maturation and defoliant application [10,15,16,17,19,21,40,41], with the latter in agreement with the results obtained in the present study. Particularly, the study carried out by Lang et al. [10] in 355 plots of cotton cultivated from 2012 to 2019 in Xinjiang Province, China, showed very similar results to ours, with lowest RMSE and highest R² for the fourth and fifth months after sowing (90-120 to 120-150 DAS).

Several studies as those included in Table 1 evaluated the best period for acquiring Satellite and UAV images or ground measurements for best cotton yield predictions. Some works have used a single date image and others have evaluated images at specific phenological stages or complete time-series data. Table 4 summarizes the best periods, in DAS, inferred from the papers listed in Table 1. The best correlations between cotton yield and the VIs varied from the flowering and boll development period [12,18,22,37,38,39] to boll opening, maturation and defoliant application [10,15,16,17,19,21,40,41], with the latter in agreement with the results obtained in the present study. Particularly, the study carried out by Lang et al. [10] in 355 plots of cotton cultivated from 2012 to 2019 in Xinjiang Province, China, showed very similar results to ours, with lowest RMSE and highest R² for the fourth and fifth months after sowing (90-120 to 120-150 DAS).

The best forecasting models for mid to late-season (90 to 150 DAS) for a regional scale can be useful as advance information for institutions as commodity traders, policymakers and governments, and in some instance, for farmers, mainly regarding their logistics and for planning the next crop season (for example, planning plot fertilizer reposition for the next crop based on the previous crop nutrients exportation). For some in-season interventions, as topdressing fertilizer applications, some forecasting models may be impractical, since most of the topdressing fertilizers (macro and micronutrients) are commonly applied before 100 DAS. Nevertheless, earlier prediction using peak models still provides an opportunity window for topdressing fertilizer intervention, since at about 80 DAS around 45%, 50% and 80% of sulfur, potassium and nitrogen, respectively, have been absorbed by cotton plants [42].

4. Conclusions

The forecasting approach based on simple regression models and time series intervals for cotton yield using the Modis product MCD43A4 V6.1 allowed in-season prediction with RMSE of about 750 kg ha^-1 at a regional scale in three cotton-producing states at the Brazilian Cerrado. This RMSE was higher than the ones obtained by Prasad et al. [11] and Lang et al. [10] for seed cotton yield prediction at regional scales in India and China, respectively, but close to the RMSE obtained by Filippi et al. [12] in Australian fields. The approach is simple and can be easily applied to predict seed cotton yield at farm, region and country levels. One limitation of the approach is the coarse spatial resolution of the Modis product, which restricts its use only for large commercial plots, such those in crop production system in the Brazilian Cerrado.

Among the nine VIs evaluated, EVI and TVI were the best individual predictors. Accuracies, evaluated by mean of RMSEs, were very low up to 75 DAS, which constitute a limitation to the application of such approach to estimate cotton yield during the first stages of cotton development. Best in-season prediction was provided by 15-days intervals between 90-105 to 135-150 DAS that corresponds to part of the mid-season to late season stage of development (boll development, open boll and fiber maturation). Best forecasts for earlier stages (below 90-105 DAS class) were provided by the model fitted for peaks with EVI and TVI, which occur around 80-90 DAS.

Future validation experiments should evaluate the accuracies and usefulness of some models in the next coming cotton seasons at different sub-regions, including their abilities to predict cotton net production in a whole farm, state and at regional and country levels. Model’s usefulness at farm level for some management interventions (within-season topdressing fertilizer application or planning reposition for the next crop based on exported nutrients, for instance) requires due caution and/or other validations.

The proposed approach can be extended to corn, which is widely cultivated as second harvest in the Brazilian Cerrado and likely to soybean (cultivated in the rainy season and with much more cloud interference in the satellite images) that is the main first crop cultivated in Brazil.