Regional-scale monitoring of cropland intensity and productivity with multi-source satellite image time series

Abstract

In the context of growing populations and limited resources, the sustainable intensification of agricultural production is of great importance to achieve food security. As the need to support management at a range of spatial scales grows, decision-support tools appear increasingly important to enable the timely and regular assessment of agricultural production over large areas and identify priorities for improving crop production in low-productivity regions. Understanding productivity patterns requires the timely provision of gapless, spatial information about agricultural productivity. In this study, dense 30-m time series covering the 2004–2014 period were generated from Landsat and MODerate-resolution Imaging Spectroradiometer (MODIS) satellite images over the irrigated cropped area of the Fergana Valley, Central Asia. A light-use efficiency model was combined with machine learning classifiers to assess the crop yield at the field level. The classification accuracy of land cover maps reached 91% on average. Crop yield and acreage estimates were in good agreement (R² = 0.812 and 0.871, respectively) with reported yields and acreages at the district level. Several indicators of cropland intensity and productivity were derived on a per-field basis and used to highlight homogeneous regions in terms of productivity by means of clustering. Results underlined that regions with lower water-use efficiency were not only located further away from irrigation canals and intake points, but also had limited access to markets and roads. The results underline that yield could be increased by roughly 1.0 and 1.4 t/ha for cotton and wheat, respectively, if the access to water would be optimized in some of the regions. The minimum calibration requirement of the method and the fusion of multi-sensor data are keys to cope with the constraints of operational crop monitoring and guarantee a sustained and timely delivery of the agricultural indicators to the user community. The results of this study can form the baseline to support regional land- and water-resource management.

Keywords:

• agricultural management
• cropland use intensity
• crop yield
• crop type classification
• water-use efficiency

Introduction

The most intensive cropping systems, with several harvests per year, occur in dry areas where irrigation is widespread (Salmon et al. 2015; Xiao et al. 2006). Irrigated agriculture covers 16–20% of the total arable area, but contributes to circa 44% of total crop production (Alexandratos and Bruinsma 2012). However, it is also the world´s largest source of freshwater consumption (Bastiaanssen, Molden, and Makin 2000; Sauer et al. 2010). With the soaring demography and the limited potential to expand the suitable cropped areas as backdrop, irrigation is expected to become an increasingly important practice to meet the global food demand (Wichelns and Oster 2006). Irrigation may lead to waterlogging, increased salinization (Abbas et al. 2010; Singh 2015), and land degradation (Dregne 2002; Gibbs and Salmon 2015; United Nations 2015) which in turn affect the productivity as they lower yields below their potential (Aslan and Gundoglu 2007; Bahçeci et al. 2006; Lobell et al. 2009; Singh 2010). Since many irrigation systems perform below their potential, it becomes critical to evaluate their actual productivity to assess potential yield gaps (Lobell et al. 2005).

One prominent example of irrigated agriculture in dry areas is the Fergana Valley in Central Asia where agriculture suffers from low field application efficiencies (Pereira et al., 2009; Reddy et al., 2013), and where percolation and transpiration losses render agricultural productivity ineffective (Horst et al., 2005; Reddy et al., 2013). For instance, the yields of cotton, the main crop since the 1960s (Abdullaev et al., 2009a), decreased from 4.6 t/ha in 1980 to 2.9 t/ha in 2000 (SIC-ICWC, 2016) due to widespread land degradation and soil salinization (Ibrakhimov et al. 2007). Other factors such as climate change or agricultural management practices could further lead to a significant drop in production (Abdullaev et al. 2009b; Bichsel 2009). Therefore, understanding the dynamics of water-use efficiency, crop productivity, and their constraints in both space and time would support the identification of options for increasing productivity and cropping intensity in a sustainable fashion. This is particularly notable in view of recent changes in Uzbekistan’s strategy to significantly cut cotton output by 2020 and devote unproductive lands to alternative crops such as vegetable, potato, and forage production.

The integrated use of geographical information systems (GIS), satellite imagery, crop models, and machine-learning algorithms is an efficient means for tracking the dynamics of agricultural productivity in space and time as well as providing decision-support tools (Roberts et al. 2009; Tayyebi et al. 2016a; Tayyebi et al. 2016b; Johnson 2014; Kussul et al. 2012; Matton et al. 2015). An array of indicators of agricultural productivity has been devised to assess (a) cropland dynamics (De Beurs and Ioffe 2014; Li et al. 2014), (b) cropping intensity (Estel et al. 2016; Kuemmerle et al. 2013; Löw, Conrad, and Michel 2015a; Löw et al. 2015b; Martinez-Casasnovas, Martín-Montero, and Auxiliadora Casterad 2005; Peña-Barragán et al. 2011; Shao et al. 2009), (c) biophysical variables such as biomass or leaf area index which are required for yield estimation (Alganci et al. 2014; Doraiswamy et al. 2005; Duveiller, Baret, and Defourny 2011; Haboudane et al. 2002; Li et al. 2015; Lobell 2013), (d) irrigation and water-use efficiency (Bastiaanssen, Molden, and Makin 2000; Bastiaanssen and Bos 1999; Liu et al. 2010; Sepulcre-Cantó et al., 2013), (e) soil salinity (Metternicht and Zinck 2003), and (f) land degradation (Abbas et al. 2010; Dubovyk et al. 2013; Skakun et al. 2015). In addition, remote sensing can support crop insurance through assessment (De Leeuw et al. 2014) or mitigation of crop pest infestations (Latchininsky 2013; Löw et al. 2016b; Renier et al. 2015; Waldner et al. 2015).

From a management point of view, satellite remote sensing has several advantages for operational and reliable crop monitoring. First, it has the ability to bypass or complement in situ measurements, which is specifically important in regions such as the Fergana Valley where access to such information is limited. Second, it can provide estimates at a range of spatial scales, which in turn may reduce the reliance on aggregated statistics. Third, it enables both nowcasting and hindcasting using archive imagery, enabling the understanding inter-annual variations in productivity (Lobell et al. 2005). However, one drawback of using certain satellite systems is that they provide images either too coarse for field-level assessments, e.g., MODerate-resolution Imaging Spectroradiometer (MODIS; Duveiller 2011; Löw and Duveiller 2014), or with insufficient temporal frequency that hampers continuous modelling, e.g., Landsat (Lobell et al. 2015).

Bearing these limitations in mind, this article evaluates a set of methods for crop productivity assessment at the field level and its ability to simultaneously respond to users’ needs and cope with operational constraints. Based on user consultations with regional (Scientific-Information Center of the Interstate Commission for Water Coordination in Central Asia, SIC-ICWC) and international organizations (International Center for Agricultural Research in the Dry Areas, ICARDA), the methods seek to provide end-of-season crop production estimates as well as increase the diagnostic capability and the understanding of the cropland-use intensity in the Fergana Valley. Therefore, descriptors of land-use intensity, crop yield, and water-use efficiency were thus produced on annual and multi-annual bases and used to highlight spatial and temporal patterns of productivity. Because operational monitoring cannot tolerate compromises in the integrity of the product due to the inability to collect remote sensing and in situ data or technical hazards such as sensor failure, data collection cancellation, or high cloud coverage, special care was devoted to: (a) ensure that sufficient satellite data would be available through the use of multi-sensor imagery fusion and (b) integrate methods with minimum calibration data requirements and high generalization capabilities.

Study site and data

Study site

The Fergana Valley is located in the Aral Sea Basin, a vast transboundary river basin at the heart of the Eurasian continent, and surrounded by the Alatau Range in the North, the Tian Shan Mountains in the East, and the Alay Mountains in the South (Figure 1). The valley spans in the upstream region of Central Asia’s most important irrigated region (Abdullaev et al. 2009b; Bichsel 2009). This region became most known for the desiccation of the Aral Sea (Löw et al. 2013b; Micklin 2016), caused by the expansion of irrigation agriculture along the Amu Darya and Syr Darya rivers during the Soviet period (Saiko and Zonn 2000). The Fergana Valley is located around the upper and mid-reach of the Syr Darya river catchment that generates almost 70% of surface water (SIC-ICWC 2016). It virtually represents the sole source of water for irrigation. Climate is continental dry with 100–200 mm average annual precipitation in the flat valley floor (Munoz and Grieser 2006).

Figure 1. Cropland distribution in the Fergana Valley. Rivers, intake points, and irrigation canals are shown in shades of blue. The settlements, in yellow, were derived by means of on-screen digitization of GoogleEarth™ images from 2016. Crop yield and acreage statistics were available for the districts in the eastern and south eastern parts solely (grey polygons).

Cotton is the dominant crop in the region (Abdullaev et al. 2009b), though it has been progressively supplemented by winter wheat after the independence from the Union of Soviet Socialist Republics in 1991 (Abdullaev et al. 2009b). Cotton is sown mid-April and harvested from September to October. Winter wheat is generally sown in autumn, and its harvest spans from June to July, depending on the actual sowing date. After winter wheat, a second crop such as sorghum or maize is usually cultivated. In addition, various vegetables, e.g., carrots or beans, melons, and fruit trees are grown.

Satellite data

Satellite image time series from MODIS and Landsat covering a 11-year period from 2004 to 2014 were at hand. On the one hand, 8-Day MODIS surface reflectance data (250 m) were downloaded (MOD09Q1; tile h23v04). In total, 46 images were available for each year and stacked into 11 annual time series. On the other hand, 319 Landsat TM, ETM+, and OLI images (30 m) were available (two tiles: path 152/row 032 and path 153/row 032). The average Landsat data availability was 12 images per year with a minimum of seven (2007) and maxima of 20 (2008).

Field boundaries and in situ data

Field boundaries were digitized based on GoogleEarth™ satellite images from 2016. The resulting database totalled 87,883 field entries, which covered an area of 6,362.15 km². The average field size is 6.74 ha which was found sufficiently large to identify crop classes using 30-m data (Löw and Duveiller 2014). Fields above 500 m of elevation were not digitized because: (a) they were generally too small to be resolved at Landsat’s resolution and (b) the number of usable images dropped markedly due to the more persistent cloud cover.

In 2014, in-situ data of the major crop types were randomly collected throughout the study region (Table 1). In addition, reference samples for 2010–2012 were extracted from existing field-based land cover maps (Conrad et al. 2013b; http://geonode.caiag.kg, last access 09/09/2015) and cross-checked by visual inspection of Normalized Difference Vegetation Index (NDVI) signatures. In total, seven classes were of interest: “Cotton,” “Winter wheat-Other,” “Orchards,” “Fallow,” “Rice,” “Other crops,” and “Fishponds.” The class “Winter wheat-Other” is an intra-annual sequence of winter wheat followed by a summer crop (usually rice, sorghum, or maize). Orchards include perennial woody vegetation such as fruit tree plantations or vineyards. “Fallow” fields were neither tilled nor sowed in the survey year. “Other” groups the remaining marginal crops for which insufficient reference data was available.

Table 1. Number of reference fields per year and per crop. In total, there were 3,211 fieldswith known crop type.

Class	Number of reference fields
	2010	2011	2012	2014
Cotton	240	160	262	234
Winter wheat – Other	244	266	190	200
Orchards	124	114	137	102
Fallow	77	116	115	55
Rice	67	95	61	50
Other crops	130	20	56	30
Fishponds	15	18	19	14
Total number	897	789	840	685

Crop area, yield, and ancillary data

Agricultural statistics from seven districts of the Fergana Valley including in-situ crop yield (in tons per hectare for wheat and cotton) and crop acreage, aggregated at the district level, were available (SIC-ICWC 2016) for the years 2010–2012. Several ancillary datasets were also at hand, namely: the locations of irrigation canals and roads were sourced from OpenStreetMap, the location of intake points was sourced from SIC-ICWC (2016), and the location of the settlements was obtained via Google Earth^TM imagery.

The Global Agro-Ecological Zones (GAEZ) dataset, version 3.0., which is based on the Harmonized World Soil Database (HWSD) and climate data for 1961–1990 (Fischer et al. 2002; IIASA and FAO, 2012), was taken from the Food and Agriculture Organization (FAO). It provides a description of the environmental characteristics of a region as well as the crop-suitability index for several crop types.

Methods

The assessment of crop productivity consists in five consecutive steps (Figure 2). First, Landsat and MODIS images were pre-processed and fused to produce 8-day synthetic time series at 30 m (Section 3.1). Second, crop types were identified at the field level based on NDVI signatures (Section 3.2) The third step focused on yield estimation (Section 3.3), which is the prerequisite for assessing the water-use efficiency (Section 3.4). The fifth module sought to derive several indicators of agricultural intensity and productivity (Section 3.5). Finally, these indicators were utilized to identify and map homogeneous regions of similar productivity profile by means of clustering.

Figure 2. Workflow of the method and its associated maps: crop type, crop yield, water-use efficiency, agricultural indicators, and identification of homogeneous areas.

Image data pre-processing

Multi-temporal classification often relies on vegetation indices like NDVI (Rouse et al. 1974) to capture phenological changes of different vegetation classes. Therefore, a time series of Landsat-like NDVI images was generated by blending NDVI images from Landsat and MODIS using the ESTARFM method (Gao et al. 2006; Zhu et al. 2010). An “index-then-blend” approach was implemented as it was reported to be more accurate and less computationally expensive than blending reflectance values before calculating NDVI (Jarihani et al. 2014; Knauer et al. 2016).

Two pre-processing steps were applied to reduce effects of clouds and aerosols and construct consistent time series. First, Landsat-7 ETM+ data acquired after the scan line corrector failure (31 May 2003) were gap-filled using on the algorithm developed by Chen et al. (2011) and clouds were masked based on the results of Fmask (Zhu and Woodcock 2012). Second, MODIS pixels labelled as no-data, snow/ice, or clouds were excluded based on the quality flags. To minimize the impact of outliers and restore missing values, the MODIS time series were smoothed using the approach developed by Chen et al. (2004).

The ESTARFM fusion operates at the pixel level and requires pairs of Landsat and MODIS images. Linear regressions were fitted for a central pixel within a moving window. Within the moving window, all pixels that are spectrally similar to the central pixel are considered. At least one Landsat-MODIS image pair for approximately the same day was required, followed by a MODIS image corresponding to the next day for which the fused Landsat images can be generated (Figure 3). MODIS images were re-projected to the same projection as the Landsat data (Asia North Albers Equal Area Conic, ESRI: 102,025) and resampled to 30 m using the nearest-neighbour algorithm. The available Landsat/MODIS base pairs were then fused individually with a 1,500 m by 1,500 m window size (Gao et al. 2006). The accuracy of the synthetic images was evaluated by comparing the predicted NDVI with the corresponding Landsat NDVI of 2,500 random samples and for each image pair annually. The prediction accuracy was quantified in terms of median Pearson´s correlation coefficient (R²) and median Root Mean Square Error (RMSE) across observation years.

Figure 3. Example of data fusion of NDVI could free maps in 2004. Left: original Landsat-7 image with SLC-off, middle: MODIS image, right: cloud- and gap-free synthetic image.

Figure 4. Generalized scheme of NDVI phenological metrics for the growing season. The phenological metrics are in fact statistical descriptors of the NDVI trajectory of a field.

Crop type classification and validation

Capitalizing on the field boundary dataset, a per-field classification method was implemented. For every year and every image, field-level average and standard deviation NDVI values were computed and served as predictor variables. To improve class discrimination, eight phenology metrics were additionally derived (Figure 4):

• The median NDVI value of the total growing season;

• The 25% and 75% percentile of NDVI value of the total growing season;

• Amplitude of the total growing season;

• Start and end of growing season;

• Length of the total growing season; and

• The area under the curve (AUC) between the start to the end of the growing season.

Such metrics were found to enhance the classification accuracy when used instead or in addition to raw time series data (Alcantara et al. 2012; Jacquin, Sheeren, and Lacombe 2010). The growing season was determined based on results from TIMESAT (Jonsson and Eklundh 2004).

Random Forests (RF; Breiman 2001) and Support Vector Machines (SVM; Cortes and Vapnik 1995) were calibrated as base classifiers and configured according to the decision-fusion approach developed by Löw, Conrad, and Michel (2015a). This method consists in randomly selecting subsets from the full set of predictor variables, calibrate a RF and SVM classifier for each subset, and then fuse the results according to a weighted majority-voting scheme. Weights were defined in proportion to the corresponding classification uncertainty at the pixel level (Policar 2006). In this study, each subset contained half of the predictor variables and 50 random draws were performed. To identify crop types in years without ground truth data, the reference data-sets were pooled together to calibrate one unified classifier algorithm (Hüttich et al. 2011; Löw et al. 2016b; Zhong, Gong, and Biging 2014).

The accuracy of the annual crop maps was assessed based on independent validation samples. Confusion matrices were calculated to estimate overall (Congalton 1991) and class-wise (F-score) accuracies (Van Rijsbergen 1979). In addition, receiver operating characteristic curves (ROC) and the corresponding area under the ROC curves (AUC), which is an increasingly used accuracy metric in machine learning and data mining (Fawcett 2006; Tayyebi et al. 2014; Tayyebi and Pijanowski 2014), have been calculated. By averaging the AUC over the different classes, an overall measure for the quality of the model predictions is obtained (Loosvelt, Peters, and Skriver 2012a):

$AU{C}_{avg}=\sum _{i=1}^{c}AUC\left(i\right)\ast w\left(i\right)$

with $c$ the number of classes (here c = 7), $AUC\left(i\right)$ the area under the ROC curve for land cover class $i$, and $w\left(i\right)$ the weighing factor associated with class $i$. This weighing factor is determined with respect to the contribution of each class in the test dataset.

Crop yield estimation and validation

Methods to derive crop yields range from empirical relationships between satellite-derived vegetation indices and ground-based yield measurements over official statistics (Bolton and Friedl 2013) to more sophisticated crop simulation models but with high input data requirements (Doraiswamy et al. 2004; Doraiswamy et al. 2005). Models based on light use efficiency (LUE; Monteith and Moss 1977) were shown to provide accurate yield estimates (Alganci et al. 2014; Lobell 2013; Lobell et al. 2003) and do not heavily depend on in-situ calibration data (Lobell 2013; Sadras et al. 2015). The LUE model is based on the assumption that biomass production is proportional to total absorption of Photosynthetically Active Radiation (PAR) over the course of the growing season (Vanino et al. 2015). It has several main terms: (a) the fraction of PAR absorbed by the vegetation canopy (fPAR), (b) the LUE term that describes the light-to-biomass conversion efficacy of the investigated crop, and (c) coefficients needed to convert biomass to harvestable yield, including harvest index ($\mathrm{H}\mathrm{I}$). Mathematically, it reads as:

(1) $Y_{crop}=\left(\sum _{SoS}^{EoS}fPAR\ast PAR\ast LU{E}_{act}\right)\ast HI$(1)

where ${\mathrm{Y}}_{\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{p}}$ is final crop yield (in g × m⁻²), SoS and EoS are the start and end of season, $\mathrm{L}\mathrm{U}{\mathrm{E}}_{\mathrm{a}\mathrm{c}\mathrm{t}}$ is the actual light use efficiency (in gC × MJ⁻¹),PAR is the photosynthetic active radiation (in W × m⁻²). In fact, the product in the brackets represents the daily accumulation of aboveground biomass integrated over the growing season. The fraction of photosynthetic active radiation (fPAR) for each field was derived from the synthetic NDVI values using generic NDVI-fPAR equations (Sellers et al. 1996) – an approach that already gave reasonable results and could be applied independently from fPAR field measurements (Alganci et al. 2014; Lex et al. 2015; Lobell, Ortiz-Monasterio, and Lee 2010; Lobell et al. 2003).

Due to the arid climate, the LUE had to incorporate plant stresses by applying certain stress factors, which limit biomass accumulation and thus crop yields:

(2) $LU{E}_{act}=\epsilon \ast T^{\prime }\ast W^{\prime }$(2)

where $\mathrm{L}\mathrm{U}{\mathrm{E}}_{\mathrm{a}\mathrm{c}\mathrm{t}}$ represents the daily actual light-use efficiency, $\epsilon$ is the maximum crop-specific LUE (cotton: 1.9, wheat: 2.0), ${\mathrm{T}}^{\prime }$ represents temperature and ${\mathrm{W}}^{\prime }$ water-induced stresses (both scaled between 0 and 1). The crop stress scalars were calculated from meteorological measurements and the ratio between actual evapotranspiration (ET_act) and crop evapotranspiration (ET_c). ET_act was assessed for each field using the S-SEBI (Simplified Surface Energy Balance Index) algorithm (Löw et al. 2016a; Roerink, Su, and Menenti 2000; Senay et al. 2013) and ET_c was calculated according to Allen et al. (1998). Crop yield is finally derived by applying $HI$ (cotton: 0.25, wheat 0.38) at the end of the season recalculation to t × ha⁻¹. Values for $HI$ and $\epsilon$ were provided by SIC-ICWC (2016), based on field measurements, and were within the range of other studies (Hunink and Droogers 2011; Lobell, Ortiz-Monasterio, and Lee 2010; Steduto et al. 2009; Fritsch 2013; Shi et al. 2007).

The coefficient of determination (R²), and the root mean square error (RMSE) were used for assessing the estimates of crop yields and crop acreage based on district-level statistics.

Water-use efficiency

Crop-specific water-use efficiencies were produced by dividing crop yield with the estimated seasonal water use, based on ET_act (Section 3.3), and using a straightforward approach that relies only on remotely sensed data (Biradar et al. 2008; Platonov et al. 2008; Sadras et al. 2015):

(3) $WU{E}_{crop}=\frac{Y_{crop}}{E{T}_{act}}$(3)

where $WU{E}_{crop}$is the water-use efficiency in kg/m³.

Indicators of agricultural productivity

Four indicators of cropping intensity and two indicators of productivity were computed from the annual crop type, yield, and WUE maps, namely: (a) cropping frequency, (b) fallow cycles, (c) number of annual multi-cropping, (d) major crop class, (e) ratio actual versus attainable crop yield, and (f) ratio actual versus attainable WUE (see Table 2).

Table 2. Multi-year cropland use intensity indicators, their purpose, and selected references.

Description	Formula	Purpose of use and options for management	References
Cropping frequency (CF)
Cropland-use intensity, counts the number of years $i$ a field was used/harvested. In this study, the maximum number of years $Y$ was 11 (2004–2014)	$\mathrm{C}\mathrm{F}=\sum _{\mathrm{i}=1}^{\mathrm{Y}}\mathrm{y}\left(\mathrm{i}\right)$	Facilitates identification of land affected by land degradation, and priority sites for alternative land use options. Supports assessment of potentials for intensification and the identification of abandoned land.	Alcantara et al. (2012); Dubovyk et al. (2013); Estel et al. (2016); Löw et al. (2016c)
Fallow cycles (FCI)
Number of recurring, consecutive sequences of fallow cropland (i.e., not sown, cropped, or ploughed). Recurring periods of 1 (S1), 2 (S2), and 3 years (S3) were considered. $N_S$ gives the maximum number of possible sequences in the observation period (e.g. 6 times S1). The maximum number of years $Y$ was 11 (2004–2014).	$FCI=\sum _{i=1}^{N_{S1}}\left(S1\ast \frac{N_{S1}}{Y}\right)+\left(S2\ast \frac{N_{S2}}{Y}\right)+\left(S3\ast \frac{N_{S3}}{Y}\right)$	Supports assessing management practices (fallow periods, soil regeneration) and intensity of land management.	Estel et al. (2016)
Multi cropping (MC)
Number of harvests over several years. The class wheat-other equals to 2 harvests, the class fallow equal to zero harvest, the reminder classes equal to one harvest.	$\mathrm{M}\mathrm{C}=0.5\ast \sum _{\mathrm{i}=2}^{22}\mathrm{s}\left(\mathrm{i}\right)$	Assessment of cropping intensity, crop practices and input demands (irrigation water, fertilizer)	Christmann et al. (2009); Estel et al. (2016)
Major crop class (MOD)
Most frequently (modus) cultivated crop in the observation period.	$\mathrm{M}\mathrm{O}\mathrm{D}=\bar{x}\left(i\right)$	Enables identification of monoculture, assessment of environmental impact of monoculture (e.g. land degradation, reduction of soil fertility), track diversification of cropping systems toward more integrative systems, for example double cropping with wheat or intercropping (not assessed in this study), interventions to sustain food production	Christmann et al. (2009)
Ratio actual vs. attainable crop yield (RYIELD)
Ratio of actual yield (as median crop yield over years) to attainable yield defined as the 95% percentile $Y_{95}$on the probability density function of crop yield.	$RYIELD=\frac{\overline{Y_{crop}}}{Y_{95}}$	Facilitates identification of local yield gaps and yield growth opportunities. Enables assessment of yield constraining factors.	Lobell (2013); Lobell, Ortiz-Monasterio, and Lee (2010); Sadras et al. (2015)
Ratio actual vs. attainable WUE (RWUE)
Ratio of actual WUE (as median WUE yield over years) to attainable WUE defined as the 95% percentile $WU{E}_{95}$on the probability density function of WUE.	$RWUE=\frac{\overline{WU{E}_{crop}}}{WU{E}_{95}}$	Facilitates identification of opportunities for improving WUE. Supports assessing potentials to save irrigation water, and manageable and agro-ecological system properties that determine spatial pattern of WUE.	Bastiaanssen and Harshadeep (2005); Cai et al. (2011); González-Dugo and Mateos (2008); Molden and Sakthivadivel (1999)

First, the cropping frequency (CF) is defined as the number of years $y$ out of all observed years $Y$ when a field was cropped, i.e., not classified as “fallow.” Next, consecutive sequences of fallow cycles (FCI) defined as a certain number of consecutive years without crop cultivation were flagged. The number of sequences of certain durations, i.e., 1-year (S1), 2-year (S2), or 3-year (S3) sequences across the entire observation period was summarized following a weighting scheme. Weights were defined as the ratio of the total number of years (11) and the number of maximum possible sequences (maximum six S1 sequences possible), whilst there must be at least one active year between two sequences consecutive sequence of fallow land (Estel et al. 2015). Fields under multiple cropping management were considered as fields being harvested at least twice within a single growing season. In the Fergana Valley, usually a double-cropping scheme is applied, i.e., a summer crop is planted directly after the wheat harvest. The number of seasons per year was defined based on the land cover classification, i.e., the class “wheat-other” indicated a double season, the reminder classes single season. Likewise, the multi-cropping (MC) indicator corresponds to the number of annual seasons (s) over the entire observation period. Next, the most frequent crop class (i.e., the mode, MOD) in the observation period was identified as an indicator of crop dominance. Finally, based on the maps of annual crop yield and water-use efficiency, the ratios between actual and attainable crop yield (RYIELD) and WUE (RWUE), respectively, were calculated. The attainable yield and WUE were defined as the 95% percentile on the probability density function of crop yield and WUE, respectively (Lobell 2013). The crop suitability index was derived for cotton and wheat (the prevailing crop and cultivation system in the region) considering medium input levels of irrigated production from the GAEZ dataset. These ratios were computed only for fields within the same agro-ecological zone.

Mapping regions with similar pattern of productivity

Finally, regions with similar land use intensity and productivity characteristics were identified with a K-medoids clustering (Kaufman and Rousseeuw 1990). For this purpose, the agricultural indicators listed in Table 2 were used as input (i.e., six variables per field). Based on the structural Silhouette technique, the optimal number of clusters was automatically determined (Chiang and Mirkin 2010).

Then, several distance measures were summarized per cluster to assess the impact of irrigation management-related factors on spatial variability of productivity (Lobell et al. 2005; Löw et al. 2017a; Ruecker et al. 2012): distances from the fields to canals and intake points (proxies for water supply), and distances to roads and settlements (proxies for market access).

Results

Quality of the methods

There was no significant difference between the actual and synthetic NDVI values. The predicted and observed NDVI values were highly correlated (Median R² ranging from 0.794 to 0.907, median RMSE ranging from 0.013 to 0.068). In addition, there was no apparent difference between areas at the overlap of two Landsat tiles.

The overall accuracy of the crop type maps trained with annual reference data ranged from 0.89–0.91 (Table 3). The classifier method trained with the pooled reference data (Section 3.2) achieved the highest overall classification accuracy (0.93), although this difference was not statistically significant at the 95% confidence level and according to the McNemar test (Foody 2009; McNemar 1947). Cotton and winter wheat could clearly be distinguished from other crops as indicated by high classification accuracies (of the pooled data set), 0.93 (cotton) and 0.98 (winter wheat). Classification accuracy of rice was lower (0.88). Overall, the $AU{C}_{avg}$values (≥0.90) indicate good performance compared to a random classifier.

Table 3. Overall and per-class classification accuracies (F-score) of land cover maps. Pooled refers to the classification based on training data from four years pooled together. Accuracies of individual years with test data are given separately.

Class	2010	2011	2012	2014	Pooled	Average accuracies
Cotton	0.96	0.91	0.88	0.81	0.93	0.90
Fallow	1.00	1.00	1.00	0.61	1.00	0.92
Orchard	0.86	0.85	0.85	0.72	0.91	0.84
Other	0.87	0.79	0.75	0.86	0.81	0.84
Rice	0.90	0.79	0.86	0.59	0.88	0.79
Water/fishponds	1.00	1.00	1.00	0.90	1.00	0.98
Winter wheat – Other	0.92	0.92	1.00	0.87	0.98	0.94
Overall accuracy	0.91	0.91	0.90	0.89	0.93	0.91
Multi-class ROC curves	0.93	0.92	0.93	0.90	0.94	0.92

The crop area estimates derived by remote sensing were then compared with the official statistics for seven districts (Figure 5 (a)). High R² and low RMSEs were observed for both wheat (0.862/568) and cotton (0.822/964). For all crops combined, the coefficients of determination reached 0.871 at the district level.

Figure 5. (a) Estimated crop area through the LUE-model by district and year against reported crop yields for 2010–2012. (b) Same as (a) but for crop yield.

Regarding the yield modelling, the estimates were in good agreement with the official statistics. The LUE model yielded both under- and overestimation compared to the reported yields by approximatly 10% on average (over all districts with statistical data available). Irrespective of the crop type, the relationship between predicted and observed yields reached an R² of 0.812 for all crops. The strength of the relationship decreased when considering the crop types independently: the R² and RMSE were, respectively, 0.557 and 0.536 for wheat, and 0.631 and 0.480 for cotton.

Spatial pattern of cropland use intensity

The analysis of the indicators of cropland use intensity highlighted several spatial patterns. Cotton and wheat were cultivated throughout the landscape, whereas rice covered only minor parts of the valley (2–6%), preferably in its center (Figure 6). Overall, cotton was cultivated in the central and eastern parts of Fergana, and was often in rotation with wheat. Orchards were spread all over the area but occurred more frequently in the alluvial fans and the mountain slopes.

Figure 6. Maps for assessing cropland use intensity and productivity pattern (based on yield and water-use efficiency) at the field level.

Figure 6. (Continued).

Cotton (2,235.38 km²; 34% of the cultivated area), wheat-other (1,876.64 km²; 29%), and orchards (1,407.87 km²; 22%) were the dominant crops in the Fergana Valley, which confirms the conclusions of Löw et al. (2013a). The annual crop maps and statistics indicate that the area of cotton and wheat-other, which covered a major part of Fergana Valley, displayed the most fluctuations in area across years. It is worth noting that these classes were mapped with high accuracy (>0.9; Table 3) which provides confidence in the changes in their corresponding planted area. The area of orchards appeared relatively constant over time, which highlights the reliability of the classification, as orchards are multiannual crops. Around 80% of the fields were cultivated every year. Overall, 3% of the agricultural fields had five–eight cropped years cropped and 92% of the fields even had high cropping frequencies (i.e., 9–11 years). Fields with a low cropping frequency presumably occurred in the central parts of the Fergana Valley and covered less than 5%. The fallow cycle maps showed that 16% of the cropland had at least one out of three fallow cycles. Two percent of the fields had a high fallow periodicity, which occurred mostly in the central parts of Fergana.

Assessing the annual crop maps from 2004–2014 revealed that double cropping was a widespread practice (Table 4). About 84% of the fields were cultivated with wheat-other (i.e., double-cropped) at least once during 2004–2014. Fields with 9–11 double-cropped years totalled to only 2% and were scattered all over the valley.

Table 4. Crop class statistics in Fergana Valley. Cultivated area per crop class are given in square kilometres. Statistics were derived from annual crop maps.

Year	Crop class (km^2)
	Fallow	Cotton	Rice	Other	Orchard	Winter wheat – Other	Water
2004	291.06	2387.71	348.8	183.82	1339.63	1795.51	15.61
2005	232.8	2363.48	457.22	242.07	1343.79	1707.27	15.52
2006	259.54	2523.62	123.24	227.82	1325.08	1887.33	15.52
2007	285.26	1957.95	661.85	175.6	1333.97	1931.91	15.61
2008	454.55	2225.54	181.5	341.5	1355.71	1787.84	15.52
2009	265.57	1495.41	570.48	239.75	1350.17	2425.24	15.52
2010	169.32	2294.58	193.89	344.48	1364.51	1979.86	15.52
2011	336.81	2451.2	125.02	239.09	1378.22	1815.97	15.84
2012	172.53	2563.05	198.99	379.29	1361.25	1671.28	15.78
2013	227.59	1948.27	631.63	289.11	1364.9	1880.27	20.38
2014	253.97	2378.35	239.89	323.12	1385.09	1760.54	21.19
Median	259.54	2363.48	239.89	242.07	1355.71	1815.97	15.61
Average	268.09	2235.38	339.32	271.42	1354.76	1876.64	16.55
Standard deviation	78.97	317.43	206.19	68.22	18.42	204.28	2.11

Spatial patterns of crop productivity

The deviation of crop yields from the estimated regional average yields (Figure 6) displayed two patterns. First, crop yields exhibited a high degree of heterogeneity, often with more than a two-fold difference, even at short distances. Second, lower yields tended to occur in the central parts of the Fergana Valley, where average deviations amounted to 10% and even reached 40%. Cotton displayed a broader distribution than wheat, with fields achieving more than 200% of average yield. The most productive cotton and wheat fields achieved 229% and 190% of the long-year median yields (2004–2014), respectively. However, this applied for less than 1% of all fields. Overall, the spatial patterns of water-use efficiency and crop yield were similar.

Identifying regions with similar productivity characteristics

Three regions of distinct properties regarding cropping intensity and productivity were identified based on clustering (Figure 7). The first region covered 35.4% of the cropped area(excluding fields where cotton and wheat were never cultivated) and was characterized by a relatively high productivity – 104% and 109% cotton and wheat yield of the regional long-term average, respectively. Both the water-use efficiency and the cropping frequency were the highest observed. At the same time, fields belonging to this region tended to be closer to irrigation infrastructure and markets, as indicated by a low median distance (Figure 7).

Figure 7. Homogeneous regions (clusters) with similar cropping intensity and productivity properties, based on the portioning-around-medoids clustering. Clustering was not performed for those fields never planted with cotton or winter wheat.

Whilst productivity in the third region (54.9% of all fields) tended to be in between the other two regions (Figure 8), the second region was characterized by comparatively low productivity (9.7% of observed fields). Average yields of wheat and cotton were only 73% and 61%, respectively, of the regional, long-term average. It occurred presumably in the central parts of the Fergana Valley, where distances to irrigation infrastructure and markets were high (Figure 8).

Figure 8. For three regions (clusters, see Figure 7), (a) distances to irrigation canals, intake points, settlements, and roads; (b) CF = cropping frequency, MC = multi cropping; (c) median values for AYW/AYC = average attainable yield for wheat/cotton, AWW/AWC = average water-use efficiency for wheat/cotton (Table 1), FCI = fallow cycles. Colours of the bars correspond to clusters in Figure 7.

Discussion and perspectives

Generalization potential of the method

This study established that fused time series of different satellite images (MODIS and Landsat) allowed capitalizing on the spatial resolution of Landsat while ensuring a high observation frequency with MODIS (R² ranging from 0.794 to 0.907). Thus, data fusion could be a suitable means to mitigate the inherent trade-off between spatial and temporal resolution in remote sensing-based agricultural monitoring (Atzberger 2013). This is especially critical when implementing an operational system as the data flow has to be sustained. From a crop mapping point of view, Landsat imagery is well-suited for inventory of crop types in agro-ecosystems characterized by fields of 6.7 ha on average, as observed in the Fergana Valley (Löw and Duveiller 2014). Its 16-day revisit-cycle, however, generally lengthened by cloud contamination, can jeopardize the high spatial resolution benefit of Landsat images for crop monitoring, and lead to classification uncertainty (Löw et al. 2013a). In contrast, MODIS has a high revisit frequency, but its coarse spatial resolution can limits its application to monitoring of large fields only (Duveiller, Baret, and Defourny 2011; Löw and Duveiller 2014), especially in fragmented landscapes (Waldner and Defourny 2017).

Compared to previous studies in Fergana (Löw et al. 2017a), the data fusion enabled assessing crop productivity over a long period due to creating a regular and high-resolution time series. This research provides the first continuous assessment at the field level in Fergana. At synoptic scale, the cropping patterns were in agreement with those found in previous studies, even though they mostly focused on 1–4 years only (Conrad et al. 2013b; Löw et al. 2017b; Löw et al. 2013a). Further, the accuracy of the crop maps was comparable or higher than other object-based crop classifications in Central Asia (Basukala et al. 2017; Löw et al. 2015b; Löw, Knöfel, and Conrad 2015c) and other regions (Balaguer et al. 2010; Duro, Franklin, and Dubé 2012). Whilst in the future, errors in the maps are likely to be reduced by incorporating more in-situ data for the training of the classifiers (e.g., Rodriguez-Galiano et al. 2012; Waske et al. 2010) or achieving a better balance in the training class proportion (Waldner, Jacques, and Löw 2017), the current low error rates offered already much prospective for an accurate crop type monitoring, and provide a suitable basis for yield and WUE estimation.

The accuracy of yield estimates suggests that the light-use efficiency model provided realistic information on the spatial variability in the region. Long-term average crop yields of cotton (3.04 t/ha) and wheat (4.65 t/ha) were in the range of other studies in Central Asia (Abdullaev et al. 2009b; Reddy 2012; SIC-ICWC (2016)). Deviations between estimated and reported yield statistics were in the range of previous attempts (Shi et al. 2007; Fritsch 2013; Lobell et al. 2003; Reeves, Zhao, and Running 2005). Further, Doraiswamy et al. (2005) and Alganci et al. (2014) reported a 10% deviation (compared to county-yield statistics) using a crop growth simulation model for corn and soybean, and a LUE model for cotton, respectively.

Overall, a high crop classification accuracy in combination with relatively low uncertainties in the yield estimates demonstrates the suitability of the combined use of satellite image time series and light-use efficiency models for assessing agricultural production. The operational advantage of the approach is the relatively low requirement of input parameters compared to approaches that apply crop simulation models.

Spatial patterns of productivity

Homogeneous regions in terms of crop intensity and productivity were identified based on clustering. Areas closer to intake points (regions 1 and 3, respectively) tended to have persistently higher crop yields and water-use efficiency. Higher crop yields (up to 10% higher than regional average) and cropping intensities presumably occurred in areas with a shorter distance to irrigation canals and were doubly cropping was frequent (especially region 1), while lower yields occur in areas with a higher distance to irrigation canals (especially region 2). Tail-end locations within the canal network in Central Asia usually are less endowed with water than head-end fields, which in turn lowers crop yields (Ruecker et al. 2012). This underlines the importance of differential access to irrigation water at head- and tail-end locations in an irrigation network (Conrad et al. 2013a). Likewise, agricultural production in region 2 appeared to be under water deficits and/or low application efficiencies, as water-use efficiency was lowest and the number of fallow cycles highest. It conforms a study on smaller scale in the south Fergana Valley canal region, where fluctuations in water delivery rates were higher in locations with a higher distance from the water intake points, i.e. in tail end locations (Abdullaev et al. 2009b). Likewise, disparities of water-use efficiency in region 2 was higher (standard deviation across wheat fields: 36%, cotton: 10%) than in region 1 (12% and 7%, respectively), which indicates that infrastructure deterioration and location within the irrigation system are major limiting factors of crop productivity, as reportedly considered in Central Asia (Abdullaev et al. 2009b; Cai, McKinney, and Rosegrant 2003). Small-scale fluctuations of crop yield, which can be in the range of 1 t/ha for neighbouring fields, can be understood in terms of variability within districts (Ruecker et al. 2012; Shi et al. 2007). In “dryer” years (2008, 2009 or 2014), the crop yield for cotton tended to be lower. Such a trend for wheat did not emerge because it generally requires less water to complete its cycle (Bekchanov, Karimov, and Lamers 2010).

Higher distances to settlements and roads are likely to cause another major constraint for maintaining the irrigation network. Factors related to market accessibility were partly more influencing than the status of irrigation infrastructure. This was indicated by the fact that despite fields in region two tended to be closer to canals (1,607 m) and intake points (21,246 m) than in region 1 (1,941 m and 22,575 m, respectively; Figure 7). This is likely to reflect the inability (e.g. due to financial constraints) to maintain the canal infrastructure, specifically in remote areas and also caused by low technical efficiencies of the production technologies (Karimov and Niño-Zarazúa 2014). Options for improving these technical efficiencies are determined by region, crop growing area, location and soil fertility. The results underline that yield could be increased by roughly 1.0 and 1.4 t/ha for cotton and wheat, respectively, if the access to water in region 2 reaches the same level as in region 1.

Limitations and perspectives

Several limitations of the system need to be considered when interpreting its outputs. First, the map uncertainty might be distributed unequally (Loosvelt et al. 2012b; Löw, Knöfel, and Conrad 2015c; Radoux et al. 2010), although the classification accuracy of the classes fallow, cotton, and wheat-other was high. A direct validation of the cropping indicators remains challenging because it would require retrospective reference data at the field-level for every year. Second, crop yield estimates were validated with aggregated statistics and deviations at the field level tend to be higher than those reported at the district level (Biradar et al. 2008; Platonov et al. 2008). In addition, within-field variations of crop types, which in rare occasions can occur, or crop yield are not represented in the maps and remains a challenge for yield gap analysis (Sadras et al. 2015).

This study facilitated the identification yield-limiting factors and might be used as baseline for guiding mitigating them in the future (Lobell 2013; Löw et al. 2017a). This in turn is a prerequisite to target ground efforts and identify opportunities for improving productivity and subsequently farmers’ income. Water-management organizations may now use this information to identify water saving opportunities, enhance productivity, or assess the effectiveness of water distribution (Reddy 2012).

Agricultural indicators still need to find their way into practice. They are at best provided as digital and openly accessible maps for land and water resource managers, or local and regional authorities and agencies. Further, such information must be provided at the field level and over time to continuously trace productivity, identify cold/hot spots productivity or regions threatened by frequent water scarcity, and finally measure target compliance, e.g., not only by policy and management decisions but also according to sustainable development goals (United Nations 2016). To fulfil this purpose, all the maps produced are available through an online tool that is being developed for Central Asia (http://geoagro.icarda.org/en/gis/details/index).

Conclusions

In Central Asia, closing the yield gap while improving the water-use efficiency is one of the chief challenges of the next decades. Better knowledge of the dynamics and spatial patterns of cropland use intensity and productivity as well as their limiting factors would support decision-making. The overarching objective was to develop a method to achieve annual and multi-annual assessment of agricultural production in the Fergana Valley, Uzbekistan. To ensure reliability, timeliness and accuracy, particular attention was paid to implement an approach with minimum calibration data requirements and sustained remotely sensed data flows. Fused time series of MODIS and Landsat images were produced to provide the high frequency gapless 30-m data required for crop classification, yield, and water-use efficiency estimation. Per-field indicators of cropland use efficiency and productivity were then derived. Results exhibited good agreement with reported crop areas and yields across Fergana. As archive satellite data allows hindcasting dynamics of agricultural productivity, the identification of spatial patterns of cropping intensity and productivity became possible. Such maps could not only serve as baseline to identify regions for intensifying agriculture, but also to regions potentially affected by land degradation or irrigation infrastructure deterioration, which entails to equip agricultural managers and decision makers with knowledge of the most important indicators of agricultural productivity, and their constraining factors. Results demonstrated that remotely sensed estimates of agricultural productivity could provide a unique perspective that would enhance the ability to identify management priorities for improving regional production and/or reducing environmental impacts. This work settles the ground for the transition of satellite-driven crop monitoring from research to operations with minimum calibration data requirements. The proposed approach could be implemented in other irrigated regions in a cost-effective fashion, by adapting regional meteorological and crop-specific parameters.

Highlights

• Implemented the first earth observation-based monitoring for tracking dynamics of agricultural productivity in the Fergana Valley at the field level

• Provided maps about land use, crop yield, and water-use efficiency annually and synthetized in multi-year metrics for the first time

• Identified regions with similar land-use intensity and productivity characteristics using clustering techniques

• Found agricultural productivity impacted by the status of irrigation infrastructure and distance of fields to markets

• Provided management options for improving crop production in regions characterized by low productivity

Acknowledgements

We would like to express our sincere gratitude to Prof. Dr Victor Dukhovny and Dr Galina Stulina (SIC-ICWC) for providing data (water inflow and reference data) and profound feedback to this study. Further, we thank the various reviewers for their valuable work on bringing the article into its present form.

Disclosure statement

No potential conflict of interest was reported by the authors.