Winter wheat mapping combining variations before and after estimated heading dates

ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

Contents lists available at ScienceDirect

ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs

Winter wheat mapping combining variations before and after estimated heading dates Bingwen Qiu a,⇑, Yuhan Luo a, Zhenghong Tang b, Chongcheng Chen a, Difei Lu a, Hongyu Huang a, Yunzhi Chen a, Nan Chen a, Weiming Xu a a Key Laboratory of Spatial Data Mining & Information Sharing of Ministry of Education, National Engineering Research Centre of Geospatial Information Technology, Fuzhou University, Fuzhou 350002, Fujian, China b Community and Regional Planning Program, University of Nebraska-Lincoln, Lincoln 68558, NE, USA

a r t i c l e

i n f o

Article history: Received 1 September 2015 Received in revised form 27 August 2016 Accepted 14 September 2016

Keywords: Time series analysis Winter wheat MODIS EVI2 Phenological stages Intra-class variability

a b s t r a c t Accurate and updated information on winter wheat distribution is vital for food security. The intra-class variability of the temporal profiles of vegetation indices presents substantial challenges to current time series-based approaches. This study developed a new method to identify winter wheat over large regions through a transformation and metric-based approach. First, the trend surfaces were established to identify key phenological parameters of winter wheat based on altitude and latitude with references to crop calendar data from the agro-meteorological stations. Second, two phenology-based indicators were developed based on the EVI2 differences between estimated heading and seedling/harvesting dates and the change amplitudes. These two phenology-based indicators revealed variations during the estimated early and late growth stages. Finally, winter wheat data were extracted based on these two metrics. The winter wheat mapping method was applied to China based on the 250 m 8-day composite Moderate Resolution Imaging Spectroradiometer (MODIS) 2-band Enhanced Vegetation Index (EVI2) time series datasets. Accuracy was validated with field survey data, agricultural census data, and Landsat-interpreted results in test regions. When evaluated with 653 field survey sites and Landsat image interpreted data, the overall accuracy of MODIS-derived images in 2012–2013 was 92.19% and 88.86%, respectively. The MODIS-derived winter wheat areas accounted for over 82% of the variability at the municipal level when compared with agricultural census data. The winter wheat mapping method developed in this study demonstrates great adaptability to intra-class variability of the vegetation temporal profiles and has great potential for further applications to broader regions and other types of agricultural crop mapping. Ó 2016 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.

1. Introduction Wheat is one of the major grain crops in China and in the world (Franch et al., 2015; Ren et al., 2008; Upadhyay et al., 2015). Timely and accurate winter wheat mapping is very important for both food security and environmental sustainability. Global land cover products are available in China, e.g. the 1 km Global Land Cover 2000 (GLC2000) (Bartholomé and Belward, 2005) and 500 m Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Collections (MAD12C5) (Friedl et al., 2010). However, more specific information on crop type is still limited in existing land

⇑ Corresponding author at: Science Building, Floor 13th, Gongye Road 523, Fuzhou University, Fuzhou 350002, Fujian, China. Fax: +86 591 87892536. E-mail address: [email protected] (B. Qiu).

cover products (Dong et al., 2015). Therefore, new methods should be explored to provide additional information on land-cover data (Sun et al., 2012). Remote sensing has proved to be an effective tool for individual crop mapping (Pan et al., 2012). Given that every existing crop has a unique phenology during the growth season (Pan et al., 2012), multiple-temporal remote sensing data are efficient in characterizing and extracting individual crops (Lunetta et al., 2010; Qiu et al., 2015, 2014; Upadhyay et al., 2015). Considerable numbers of supervised algorithms have been developed through calculating similarities with the temporal profile of standard vegetation indices (VIs) from different agricultural crops (Gumma et al., 2015; Sun et al., 2012). Sun et al. (2012) developed a curveshape based algorithm to map winter wheat areas and demonstrated the usefulness of the temporal signatures in time-series

http://dx.doi.org/10.1016/j.isprsjprs.2016.09.016 0924-2716/Ó 2016 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.

36

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

MODIS vegetation index data for winter wheat classification. This kind of time-series based approach measures similarity directly without any transformation. However, the biggest challenge is the intra-class variability of VIs temporal profiles (Lunetta et al., 2010; Qin et al., 2015). The remote-sensing vegetation indices time series curves are closely related to crop phenology (Sakamoto et al., 2010). And crop phenology is significantly influenced by many factors such as latitude and altitude (Zhang et al., 2014). Therefore, crop profiles have been found to vary across different agroregions in different years (Foerster et al., 2012; Siebert and Ewert, 2012). For one specific agricultural crop, the temporal profiles of vegetation indices generally display considerable intra-class variability due to variations of crop calendars, natural conditions, and tillage management (Qiu et al., 2015; Wardlow et al., 2007). Besides the original time-series based methods, there are other two categories: transformation-based approach and the metricbased approach (Lhermitte et al., 2011). The latter two approaches assess vegetation similarity based on specific components or statistical characteristics of a time series (Lhermitte et al., 2011). Phenological metrics derived from VI time series have been successfully applied to enrich the information for vegetation mapping and monitoring (Wardlow and Egbert, 2010; Zhong et al., 2011). Nevertheless, phenological metrics were not selected for the top ten features for crop identification, because of phenological character variations in crop development schedules (Hao et al., 2015). Besides phenological metrics, another winter wheat mapping method specifically calculates the winter bump to its growing season (Pan et al., 2012). Under the assumption of the winter bump method, winter wheat can be discriminated from others by two distinctive peaks observed from temporal profiles of the vegetation indices during its growing period (Wang et al., 2015). However, the uniqueness and recognizability of the winter bump over large areas has not been fully exploited. The winter bump of winter wheat might be less obvious due to weak growth during winter season in some regions with considerably colder weather conditions (Wang et al., 2015). Most of the previous studies were conducted on relatively small areas (Wang et al., 2015). Little research has been paid attention to the intra-class variability over large regions. Few crop mapping studies have accounted for the variations of crop profiles resulting from specific weather and management conditions (Foerster et al., 2012). In order to fill this gap, this study developed a new approach to identify winter wheat over large regions through combining variations before and after estimated heading dates. The winter wheat mapping algorithm and its applications in China were illustrated in the following sections.

2. Study area and datasets 2.1. Study area Our study area is located between latitude 29.08°N–42.67°N and longitude 105.18°E–122.70°E, in China. It covers the North China Plain and the Huanghuai Plain, the most important winter wheat production region (Ren et al., 2008). It comprises 10 provincial-level administrative units (8 provinces plus 2 directcontrolled municipalities) (Fig. 1): Hubei, Anhui, Jiangsu, Shandong, Henan, Shanxi, Shaanxi, Hebei provinces, Beijing and Tianjin municipalities (Lv et al., 2013). There are 456,960 km2 of croplands. The whole terrain is composed of three parts: mountains in the west, hills in the south, and plains in the east and middle areas (Fig. 1). The elevation ranges from around 0 to 3748 m. The study area is separated by the Qin Mountains and the Huai River. The north side of the study area is dominated by a temperate con-

tinental monsoon, and the south side of the study area is characterized by a subtropical continental monsoon. Winter wheat is sown during October-November and harvested during June-July in the following year (Jia et al., 2013). 2.2. Datasets and preprocessing 2.2.1. The MODIS EVI2 time series datasets The 250 m 8 day composite MODIS surface reflectance products from 2001 to 2013 were utilized. The 2-band Enhanced Vegetation Index (EVI2) was calculated by the surface reflectance values from red (620–670 nm) and near infrared red bands (841–876 nm) (Jiang et al., 2008). The EVI2 time series datasets were applied since they eliminated most of the problems associated with sub-pixel and residual clouds (Jiang et al., 2008). 2.2.2. Landsat images and winter wheat extraction Landsat images in six selected regions were applied for winter wheat extraction. Three periods of Landsat 7 Enhanced Thematic Mapper (ETM+) or Landsat 8 Operational Land Imager (OLI) images were collected for each selected region: summer-autumn, winter and late spring seasons. The specific acquisition dates for each region were listed in Table 1. Images acquired in 2002 were Landsat 7 ETM+ images, and the others were Landsat 8 OLI images. Based on the vegetation phenology, forests exhibited distinct features in summer-autumn images. Winter crops including winter wheat displayed obvious characteristics in winter images and winter wheat also demonstrated distinct features in late spring images. Therefore, a three-step procedure was applied to identify winter wheat from these three periods of Landsat images. First, the images in summer-autumn were used to classify forests. The forests were distinguished by high vegetation coverage in the summer-autumn. In remote sensing images, forests were shown large values in the near infrared red (NIR) band, medium values (only larger than non-vegetation) in the short wave infrared band 1 (SWIR1) and low values in the blue band. The non-vegetation was extracted based on the low values in the NIR and the SWIR1. Consequently, the forests and non-vegetation were excluded. Second, the winter crops were classified with images in winter. Winter crops were characterized by considerably high vegetation coverage characteristics in the winter: high values in NIR, medium values in SWIR1 (only larger than water bodies), and low values in the blue band. Water bodies can be extracted by extremely low values in the NIR band. Third, the Landsat images in the late spring were utilized to extract winter wheat. Winter wheat exhibited very high vegetation coverage characteristics during its vigorous growth period in late spring, similar to the spectral characteristics of winter crops in images in the winter. Field survey data together with these above vegetation coverage characteristics were applied to the training data for different vegetation types. The support vector machine (SVM) algorithm was applied for the classification process. The winter wheat areas were finally extracted and further improved through visual interpretation and digitization. A total of 596 points were randomly generated to validate the results. Among them, 551 points were correctly classified based on references interpreted from Google Earth images. Results showed that the Landsat-interpreted images of winter wheat met the precision requirements. We made our best efforts to select images in the corresponding seasons in six selected regions. However, it was limited on the accessibility and quality of Landsat images in those regions. For example, there was no Landsat image in 2001 in region V and VI. Due to the limited Landsat images free from cloud covers and strips during the period of 2012–2013, Landsat images in 2013–2014 were applied for classifications of winter wheat in 2012–2013. Reference sites from Google Earth images was applied again to verify

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

37

Fig. 1. The spatial distribution map of study area and ground survey sites.

Table 1 Acquisition dates of Landsat images. Region

Path/Row

Acquisition dates

Ⅰ Ⅱ III IV V

123034 123036 120037 124038 123033

VI

122038

20131206, 20130613, 20131201, 20131010, 20020319, 20140819, 20020413, 20150425,

20140429, 20140123, 20140424, 20131229, 20020522, 20141202 20020531, 20141005

20140904 20140429 20140611 20140404 20021013, 20140429, 20021209, 20150103,

the changes from 2012–2013 to 2013–2014. It revealed that little changes have taken place in those regions during these two neighboring years, with the rate of change only 1.01%. It was also confirmed from the agricultural census data. Regions I-VI approximately covered 24 municipalities. The total sown areas of winter wheat from agricultural census data was 69,198 km2 in 2012 and 69,490 km2 in 2013, respectively. Landsat images in region VI in 2014–2015 were applied due to heavy cloud cover for most Landsat images in 2014. It revealed that little have taken place in region VI, with a rate of change 0.90%. 2.2.3. Field survey data The field survey data were gathered in early August 2012; early February, late April and late July 2013; mid-January and August 2014, and mid-February, early August and December 2015, respectively. UniStrong G138 or MG858 hand-held GPS receivers were utilized for ground survey in field sites. At each sampling site, the cropping patterns (e.g., winter wheat plus maize) and their cor-

responding phenological stages were recorded. A total of 1082 ground truth sites were collected (see locations in Fig. 1). About two-fifths of these collected ground sites were applied to establish the standards for identifying winter wheat, and three-fifths were utilized for accuracy evaluations. Of the 1082 collected ground sites, 475 sites were cultivated with winter wheat plus maize/bean/peanut. Another 289 nonwinter wheat crop sites were planted with single rice, double rice, spring maize and vegetables. A total of 318 field survey sites were covered by natural vegetation, primarily deciduous broadleaf forests and evergreen needle-leaf forests. A large proportion of field survey sites (920) were located on plains. These sites on the plains were generally homogeneous; their average parcel size was normally larger than 500 m
500 m. Other field survey sites in mountains and hills might be more heterogeneous. 2.2.4. Crop calendar data, agricultural census data and other datasets Crop calendar data from 110 agro-meteorological stations were obtained. The crop calendar data recorded the sowing, seedling, tillering, heading, and harvesting dates of winter wheat. The data were provided by the National Meteorological Information Center, China Meteorological Administration. The agricultural census data at the municipal level were provided by the National Bureau of Statistics of China (NBS) (http://www.stats.gov.cn). Other datasets included the 90 m DEM datasets from the Shuttle Radar Topography Mission (SRTM) (Fig. 1). 3. Methodology A new method for mapping winter wheat through Combining variations Before and After estimated Heading dates (CBAH) was

38

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

proposed (Fig. 2). It included the following procedures: data preprocessing, analyzing the EVI2 temporal profiles, dealing with challenges of the intra-class variability, winter wheat extraction and accuracy evaluation. The entire procedure was executed using the Matlab software package (The MathWorks, Natick, Massachusetts, USA).

Second, troughs were examined in late May or June in the next year, revealing the harvesting season. Third, the length of the growing period of winter wheat was longer than that of other crops. Therefore, a winter wheat mapping method was developed through combining variations before and after the estimated heading dates.

3.1. Data preprocessing

3.3. Dealing with challenges of the intra-class variability of EVI2 temporal profiles

A smoothed daily continuous EVI2 time series was developed through the following steps. First, observations with cloud contaminations were identified with references to the Quality Assurance (QA) layers. Second, these data composites with cloud contaminations were discarded. In particular, pixels recognized as serious cloud contaminations were abandoned (e.g., 6 continuous data composites with cloud contaminations during the growing period of winter wheat, indicating that there was no good observation over 48 days). Third, a daily continuous EVI2 time series was established through linear interpolation based on cloud-free observations and then smoothed with the Whittaker Smoother (WS) (Qiu et al., 2015). The Whittaker Smoother was applied due to its good performance in most cases compared with other commonly-utilized noise reduction methods (Atkinson et al., 2012). 3.2. Analyzing the EVI2 temporal profiles of winter wheat and other crops 3.2.1. Distinct intra-class variability of EVI2 temporal profiles of winter wheat The EVI2 temporal profiles of winter wheat from four sites, sites A (111°380 800 E, 31°430 5200 N), B (116°370 1200 E, 33°180 2200 N), C (119°130 2800 E, 36°460 1600 N) and D (117°160 5100 E, 39°420 3300 N) (see locations in Fig. 4), are shown in Fig. 3(a). Distinct intra-class variability was exhibited in the EVI2 temporal profiles of winter wheat. The heading dates of winter wheat were obviously delayed from south to north, and varied from April to May. Additionally, for winter wheat, the curve shape of EVI2 temporal profiles also varied at different locations (from site A to D). The dates on which winter wheat bumps occurred probably changed from November in the former year to February in the next year, depending on specific climatic conditions in different regions. The magnitudes of these bumps/local maximums in the winter season were varied in different studying regions, ranging from 0.15 to 0.3. 3.2.2. Similarities and dissimilarities in EVI2 temporal profiles between winter wheat and others The EVI2 temporal profiles of other primary vegetation types are shown in Fig. 3(b). Two sites for each vegetation type (single rice, double rice, vegetable and forest) were selected, respectively. Similar to winter wheat, there was considerable intra-class variability in their EVI2 temporal profiles: the curve shape and heading/peaking dates of each specific vegetation type varied at different sites. For each specific crop, the growth rates and magnitudes of EVI2 changed at different sites. For example, the heading dates of most agricultural crops were observed in the summer. However, the heading dates of some vegetables could also be examined in spring. Similar to winter wheat, local maximums might also be detected in winter from their EVI2 temporal profiles of other crops (Fig. 3(b)). Therefore, it is not appropriate to identify winter wheat through winter bumps. In spite of this, there were still considerable dissimilarities between winter wheat and others (Fig. 3). The following three features were distinguishable from the EVI2 temporal profiles of winter wheat. First, local minimums were examined during the period from mid-October to late November, indicating the seedling date.

3.3.1. Developing trend surfaces of two key phenological parameters for winter wheat The phenological dates generally increased with altitudinal and latitudinal gradients (Fig. 4). Trend surfaces were developed for two key phenological parameters: the heading date and the length of the early growth stage (from seedling to heading dates). Their trend surfaces were established by altitude and latitude through linear equations. The parameters of the equations were estimated by the crop calendar data from agro-meteorological stations. The equations for the heading date and the length of early growth stage are provided as the following functions (1) and (2), respectively.

y ¼ 2:8027  latitude þ 0:0065  altitude þ 18:5197

ð1Þ

y ¼ 7:3571  latitude þ 0:0222  altitude  76:0539

ð2Þ

For each pixel, these two phenological parameters were derived based on the latitude and altitude. The harvesting date of winter wheat was calculated as 48 days later than the estimated heading dates with references to phenological data from agrometeorological stations. 3.3.2. Designing phenology-based indicators revealing variations during the early and late growth stages For each pixel, two phenology-based indicators were derived based on the following procedures (Fig. 5): First, two growth stages were decided based on the estimated seedling, heading, and harvesting dates of winter wheat (represented in dots in the upper subplot in Fig. 5). The early growth stage was defined as the period from estimated seedling to heading dates (see the central-below subplot in Fig. 5). The late growth stage was determined as the period from the estimated heading to the harvesting dates (see the right-lower subplot in Fig. 5). Second, the EVI2 values were extracted according to the estimated seedling, heading, and harvesting dates for winter wheat (EVI2seedling, EVI2heading, and EVI2harvesting in Fig. 5). Third, the maximum and minimum EVI2 values and changes in amplitude during the early (see the central-lower subplot in Fig. 5) and late growth stages (see the right-lower subplot in Fig. 5) were identified, respectively. Finally, two phenology-based indicators were derived based on EVI2 Variations during the early and late growth stages. One indicator represented the EVI2 Variations during the Early growth stage (EVE for short). The other indicator denoted the EVI2 Variations during the Late growth stage (EVL for short). The first indicator, EVE, was designed through summing up two items: the EVI2 difference between the estimated heading (EVI2heading) and seedling dates (EVI2seedling), and the EVI2 change amplitude (EVI2max1  EVI2min1) during the early growth stage (function (3)). Similarly, the second indicator, EVL, was derived based on another two items: the EVI2 change between the estimated heading (EVI2heading) and harvesting dates (EVI2harvesting), and the EVI2 change amplitude (EVI2max2  EVI2min2) during the early growth stage (function (4)).

EVE ¼ ðEVI2heading  EVI2seedling Þ þ ðEVI2max1  EVI2min1 Þ

ð3Þ

EVL ¼ ðEVI2heading  EVI2harv esting Þ þ ðEVI2max2  EVI2min2 Þ

ð4Þ

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

39

Fig. 2. Overview of methodology.

Fig. 3. EVI2 temporal profiles of winter wheat (a), and others (b) from October to June.

where EVI2seedling, EVI2heading, EVI2harvesting represents the EVI2 value at estimated seedling, heading and harvesting dates of winter wheat, EVI2max1, EVI2min1 denotes the maximum, minimum EVI2 values during the early growth stage of winter wheat, EVI2max2, EVI2min2 indicates the maximum, minimum EVI2 values during the late growth stage of winter wheat, respectively. 3.4. Winter wheat extraction For a non-winter wheat crop pixel, it was almost impossible to have the real seedling, heading and harvesting dates exactly match the estimated phenological dates of winter wheat at that location. For an unknown pixel, either when the observed heading date was

earlier or later than that estimated (Fig. 5), the EVE and EVL cannot gain large values simultaneously. This is due, in part, to the following three reasons: First, EVI2 values at the estimated heading date, EVI2heading, would be much less than expected. Second, either the EVI2max1 or the EVI2max2 was impossible to obtain large values. For example, the EVI2max1 could not acquire a large value when the observed heading date was later than that estimated. Third, the EVI2seedling and EVI2harvesting would not attain very low values concurrently, particularly the EVI2harvesting. The EVI2 values from the estimated harvesting date for winter wheat were much higher than for other non-winter wheat crops. Even when the observed heading date of the non-winter crop corresponded to that estimated by chance, the third reason was still valid.

40

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

Fig. 4. Locations of (a) agro-meteorological sites and (b) observed length of early growth stage and heading dates of winter wheat.

Fig. 5. Schematic diagrams of extracting phenology-based indicators.

Pixels of other land covers could not gain large values of EVE and EVL synchronously. Specifically, the EVL from forested sites could not obtain a large value. Both the EVE and EVL from nonvegetation or sparse vegetation sites were incapable of gaining large values due to the relatively lower EVI2 change amplitude at any stage. Therefore, it is obvious that large values of EVE and EVL were only expected from winter wheat sites. The intra-class variability of winter wheat from different regions was diminished by the proposed CBAH method. First, the intra-class variability introduced by phenological shift was eliminated through developing a trend surface of key phenological parameters of winter wheat. For candidate pixels of winter wheat, the actual phenological dates corresponded substantially to those estimated phenological dates through trend surface. Second, the intra-class variability introduced by various curve shapes of EVI2 temporal profiles was also greatly reduced. Although the curve shapes of EVI2 temporal profiles of winter wheat still differed in various regions, these two phenology-based indicators (EVE and EVL) from the candidate winter wheat sites were expected to obtain large values (Fig. 5). For the winter wheat sites, the estimated heading date from the trend surface corresponded to the maximum and the seedling/harvesting dates matched with the two minimums in the EVI2 temporal profile. The magnitudes of

EVI2 Variations of winter wheat during the early and late growth stages would be larger than those from other land covers. Therefore, two phenology-based indicators revealing the EVI2 Variations during the early and late growth stages, the EVE and EVL, were applied to distinguish them.

If ðEVE > h1 and EVL > h2 Þ; it is winter wheat:

ð5Þ

where h1 and h2 are constant. This algorithm was applied to the whole study area through a per-pixel strategy. No land use/cover data was needed to mask out the non-cultivated areas before classification. 3.5. Accuracy assessment Accuracy evaluation was conducted with the agricultural census data from the NBS, ground truth data and Landsat imageinterpreted maps in some test regions. The total areas of winter wheat estimated from MODIS images were compared to those from agricultural census data at the municipal level. For accuracy assessments of winter wheat map in 2012–2013, the sown areas of winter wheat in 2012 from the NBS were utilized. The growing period of winter wheat actually covered two neighboring years. The harvesting areas of winter wheat in 2013 were sowed in late

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

41

2012. In order to compare with agricultural census data, the number of pixels labeled as winter wheat was calculated at the municipal level from the MODIS-derived map, and then multiplied by the MODIS pixel size. The spatial consistency of MODIS-estimated map of winter wheat was evaluated with 653 ground truth sites and Landsat image-interpreted results on some test regions. From south to north, six test regions were selected (see locations in Fig. 6). Landsat images in regions I-VI were applied for accuracy assessment in 2012–2013. Landsat images in regions V and VI were also utilized for accuracy evaluation in 2001–2002. 4. Results 4.1. MODIS-derived winter wheat map in 2012–2013 by the CBAH method A total of 429 randomly selected ground survey sites were applied for determining the constants in function (5). There were distinct differences in these two indicators, EVE and EVL, between winter wheat and others. The constants for extracting winter wheat in function (5) were decided: h1 was 0.3 and h2 was 0.12. The winter wheat distribution map in 2012–2013 was finally derived by the CBAH method (Fig. 6). Winter wheat dominated in the North China Plain and Huanghuai Plain. It was principally distributed in the cultivated land in the Henan and Jiangsu provinces, the northern Anhui province, the southern Hebei province, and the western Shandong province. 4.2. Evaluation with agricultural statistical data The total MODIS-estimated area of winter wheat in 2012–2013 (Fig. 6) was 198,453 km2. A slightly overestimation was revealed when compared with the total sown areas of winter wheat from the agricultural census data (189,572 km2). The MODIS-derived areas of winter wheat in 2012–2013 were aggregated to the municipal and province levels. At the province level, fairly good agreement was found in most provinces, especially in the Anhui and Shandong provinces. Underestimates were found in several provinces cultivated with limited areas of winter wheat, e.g., the Shanxi and Shaanxi provinces. There were slightly overestimates in some provinces planted with large areas of winter wheat such as the Henan and Jiangsu provinces (Table 2). At the municipal level, the cultivated areas of winter wheat estimated from MODIS images by the CBAH method correlated well with agricultural census data, with a coefficient of determination of 0.8386 and Root Mean Square Error (RMSE) of 1.847 (Fig. 7(c)). 4.3. Spatial agreements with ground truth sites and Landsat imageinterpreted results The MODIS-derived winter wheat map in 2012–2013 was evaluated with 653 ground survey sites. For the 297 winter wheat cultivated sites, 266 sites were correctly classified (89.56% agreement). For the 356 non-winter wheat sites, 336 sites were accurately identified (94.38% agreement). An overall accuracy of 92.19% was obtained when compared with the in-situ observation data (Table 3). For the MODIS-derived winter wheat maps in 2012–2013, accuracy assessments with the Landsat-interpreted maps were conducted on regions I-IV (see locations in Fig. 6). The MODIS-estimated maps of the winter wheat areas were visually consistent with the Landsat-interpreted maps (Fig. 8). In order to implement the evaluation, the Landsat-interpreted maps of winter wheat were aggregated to 250 m with the majority resampling

Fig. 6. Winter wheat distribution map derived from MODIS through CBAH method in 2012–2013.

method. An overall accuracy of 88.86% was obtained (Table 4). The user and producer accuracies of winter wheat were 83.26% and 92.92%, respectively. The kappa index was 0.7793. 4.4. Further applications of the CBAH method during the early 21st century The proposed winter wheat mapping method was further applied during the early 21st century in order to test its robustness and efficiency. Before its applications, the feasibility of the developed trend surfaces of key phenological parameters of winter wheat in 2012–2013 needed to be evaluated. According to observed phenological data from agro-meteorological sites, from 2001 to 2013, the annual absolute difference of observed heading dates was within 1.3 days for two neighboring years. Other phenological parameters such as the sowing, heading, and harvesting dates varied less than three days for one decade during 1981– 2009 in the North China Plain (Xiao et al., 2013). Therefore, the established trend surface of the heading date in 2012 was applied to the period 2008–2012. Another trend surface of heading date of winter wheat in 2006 was developed for the period 2001–2007. The length of the early growth stage of winter wheat almost remained unchanged during the early 21st century with reference to the agro-meteorological sites. Therefore, the trend surface of the length of early growth stage in 2012–2013 was utilized for the study period of 2001–2012. Finally, the spatial distribution maps of winter wheat during these years were achieved through the CBAH method (Fig. 9, results in 2001–2002 and 2005–2006 were provided). The winter wheat maps in 2001–2002, 2005–2006, and 2012– 2013 demonstrated similar spatial distribution patterns (Figs. 6 and 9). During these years, the winter wheat dominated the North China Plain and Huanghuai Plain. The total MODISestimated area of winter wheat in 2001–2002 and 2005–2006 was 196,775 km2 and 182,475 km2, respectively. Similar to 2012– 2013, there was a slight overestimation compared with the agricultural census data from the NBS (177,611 km2 in 2001 and 171,962 km2 in 2005). The coefficient of determination between the MODIS-estimated dataset and the agricultural census data at

42

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

Table 2 Provincial level comparisons of winter wheat-planted areas between NBS reports and MODIS estimates in 2012.

NSB CBAH Curve Curve (mask) Bump Bump (mask)

Beijing

Tianjin

Shanxi

Shaanxi

Hubei

Hebei

Anhui

Jiangsu

Shandong

Henan

Total

522 179 50 35 4381 1373

1051 415 44 37 1713 1553

6885 2187 5963 3794 39,311 11,834

11,276 9297 21,847 15,200 77,936 19,734

10,655 13,795 58,597 41,949 79,699 34,604

24,044 20,971 18,022 14,917 43,930 27,896

24,155 25,793 65,883 48,133 53,808 39,140

21,326 28,798 51,854 41,000 31,366 32,273

36,259 37,897 58,888 48,485 52,363 54,290

53,400 59,119 78,233 64,480 79,429 73,616

189,572 198,453 359,381 278,029 463,937 296,313

Unit: km2. Curve/curve (mask) denoted curve-shaped based method without/with croplands mask. Bump/bump (mask) represented winter bump-based method without/ with croplands mask.

Fig. 7. Municipal-level comparisons of the winter wheat-planted areas between NBS reports and MODIS estimates based on CBAH in (a) 2001–2002, (b) 2005–2006, and (c) 2012–2013, and curve-shape based method (d) and winter bump based method (e) in 2012–2013.

Table 3 Accuracy assessment for three methods in 2012–2013 using ground truth data. Total

Winter wheat Others User accuracy (%) Overall accuracy (%)

297 356

CBAH

Curve-shape based algorithm

Winter bump-based method

Winter wheat

Others

Producer accuracy (%)

Winter wheat

Others

Producer accuracy (%)

Winter wheat

Others

Producer accuracy (%)

266 20 93.01 92.19

31 336 91.55

89.56 94.38

203 139 59.36 64.32

94 217 69.77

68.35 60.96

208 136 60.47 65.54

89 220 71.20

70.03 61.80

Note: Descriptions of curve-shape or winter bump-based methods are provided in related references (Pan et al., 2012; Sun et al., 2012; Wang et al., 2015).

the municipal level was above 0.82. Specifically, it was 0.8283 in 2001–2002 and 0.8304 in 2005–2006, suggesting relatively good agreement (Fig. 7(a) and (b)). The cultivated areas of winter wheat were considerably stable during the period from 2001–2002 to 2012–2013 (Figs. 6 and 9).

However, some local variations were examined. Two regions, one in the north (region V) and another in the south (region VI) (see locations in Fig. 6), were selected for verification. For both regions, the winter wheat areas in the MODIS-derived images decreased from 2001–2002 to 2012–2013. The reductions of winter wheat

43

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

Fig. 8. Winter wheat maps from (a) Landsat and (b) MODIS data in region I, II, III and IV (see locations in Fig. 6).

Table 4 Accuracy assessment for three methods in regions I-IV during 2012–2013 using Landsat data. Landsat images

MODIS classification CBAH

Winter wheat Others User accuracy (%) Overall accuracy (%) Kappa

Curve-shape based algorithm

Winter bump based method

Winter wheat

Others

Producer accuracy (%)

Winter wheat

Others

Producer accuracy (%)

Winter wheat

Others

Producer accuracy (%)

13,613,259 2,740,653 83.24 88.86 0.7793

1,030,870 16,464,118 94.11

92.96 85.73

12,932,607 8,535,662 60.24 69.73 0.4207

1,711,522 10,669,109 86.18

88.31 55.55

11,614,482 7,017,886 62.33 70.32 0.4146

3,029,647 12,186,885 80.09

79.31 63.46

Note: Counted in number of pixels.

Fig. 9. Spatial distribution images of winter wheat in (a) 2001–2002 and (b) 2005–2006.

44

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

Fig. 10. Comparisons of MODIS and Landsat-derived images in regions V and VI in 2001–2002 (a, b, e, f) and 2012–2013 (c, d, g, h).

Table 5 Accuracy assessment for the CBAH method in region V-VI using Landsat data. Landsat images

MODIS classification 2001–2002

Winter wheat Others User accuracy (%) Overall accuracy (%) Kappa

2012–2013

Winter wheat

Others

Producer accuracy (%)

Winter wheat

Others

Producer accuracy (%)

5,247,158 939,342 84.82 88.76 0.7939

1,076,170 10,661,780 90.83

82.98 91.90

4,317,658 980,897 81.49 89.51 0.8208

899,071 11,726,824 92.88

82.77 92.28

Note: Counted in number of pixels.

areas in region V and VI were confirmed by Landsat-interpreted maps (Fig. 10). When compared with Landsat-interpreted data in regions V and VI, the overall accuracies were 88.76% and 89.51% in 2001–2002 and 2012–2013, respectively. In 2001–2002, the user and producer accuracies of winter wheat were 84.82% and 82.98%, respectively (Table 5). The kappa index in 2001–2002 was 0.7939. In 2012– 2013, the user and producer accuracies of winter wheat were 81.49% and 82.77%, respectively (Table 5). The kappa index in 2012–2013 was 0.8208 (Table 5). The efficiency and adaptability of the proposed winter wheat mapping method was verified through its successful application in Northern China during the early 21st century. 5. Discussion 5.1. Comparisons with the curve-shape and winter bump-based algorithms The efficiency of our proposed method was further evaluated through its comparisons with two commonly applied algorithms: the curve shape-based algorithm (Sun et al., 2012) and the winter bump method (Wang et al., 2015). Obvious commission errors occurred for the latter two algorithms. Specifically, for the curve shape-based algorithm, when applied without the mask of croplands data, the overall estimated area of winter wheat in 2012–

2013 was 359,381 km2. Even applied with the mask of croplands data, an overestimate of 46% (278,029 km2) was still found when compared with the agricultural census data (189,572 km2). For the winter bump-based algorithm, overestimation was much more substantial: the total estimated area of winter wheat in 2012–2013 was 296,313 km2 when applied with the croplands mask, and it was 463,937 km2 without the croplands mask. Overestimations were confirmed in most provinces (Table 2). At the municipal level, the cultivated areas of winter wheat estimated by both methods did not correlate well with the agricultural census data. For the winter bump based algorithm with croplands mask, the coefficient of determination was 0.5123 (Fig. 7(d)), with Root Mean Square Error (RMSE) of 3.214. For the winter bump based algorithm with croplands mask, the coefficient of determination was only 0.2513 (Fig. 7(e)), with RMSE of 4.513. Overestimations of the curve shape-based and winter bump methods were also confirmed from the ground survey data (see locations in Fig. 1) and the Landsat-interpreted maps for test regions (I-IV, see locations in Fig. 6). Among the 356 ground truth sites without winter wheat, over 130 sites were classified as winter wheat by these two algorithms (Table 3). The overall accuracy for these two algorithms was around 64–66%. When evaluated with Landsat image-interpreted data (Fig. 8), an overall accuracy of around 70% was obtained when applied with the croplands mask (69.73% for the curve shape-based method and 70.32% for the winter bump-based algorithm). The kappa index was approximately

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

0.41–0.42 for these two methods (Table 3). This indicates that the curve shape-based and winter bump-based algorithms had difficulties in obtaining satisfactory results when applied to larger regions. It is because the signatures for vegetation indices were relatively stable in smaller regions (Sun et al., 2012). The temporal profiles for the vegetation indices varied across different areas and winter bumps might also occur in areas with no winter crops (Fig. 3). Compared with the curve shape-based algorithm and the winter bump-based algorithm, our proposed winter wheat mapping method incorporated the intra-class variability over large regions and thus achieved greater efficiency and accuracy for winter wheat mapping. 5.2. Prospects and future directions to deal with the intra-class variability of vegetation indices The accuracy of estimation from single or multi-temporal images was often limited by the intra-class spectral variability problem in high resolution imagery (Zhang and Weng, 2016). The challenge introduced by the intra-class variability still exists for land cover/crop mapping based on time series images. In order to address this problem, a method of temporal shifted standard MODIS-EVI time series images was developed to map maize area with reference to agro-meteorological data (Zhang et al., 2014). Another study incorporated agro-meteorological information into the classification process (Foerster et al., 2012). One recent study proposed an efficient phenology-based paddy rice mapping method with reference to temperature conditions (Qin et al., 2015). These studies open a new direction for phenology-based mapping approaches with reference to accessorial data in order to account for the intra-class variability of VIs temporal profiles. Supplementary information has the potential to improve land cover classification accuracy and efficiency (Eisavi et al., 2015). Enlightened by these studies, this study developed a spatial distribution map for key phenological dates of winter wheat based on altitude and latitude with reference to the crop calendar data. Besides the phenological shift, the temporal VIs profiles of specific agricultural crops might also differ in curve shape and magnitude. For example, in winter wheat, extreme weather conditions, such as snow or frost, could introduce the intra-class variability of VIs temporal profiles. The extreme weather in winter could hold back the growth of winter wheat, and may cause the small bumps in the winter season. But its duration varied from northern to southern regions. This kind of intra-class variability of VIs temporal profiles from winter wheat cannot be eliminated by incorporating the phenological shifts with references to crop calendar data. One promising solution is to develop metrics based on specific periods that phenological characteristics inherently exhibit (Zhang et al., 2013). Images of several key periods are sufficient for accurate crop mapping (Wardlow et al., 2007). The highest possible accuracy of land cover classification can be achieved by optimal variable combinations of a few bands and a shorter period of time-series data (Zhou et al., 2013). It is possible to propose an integrative approach based on the use of an optimal combination of indices for improving the assessment of crop dynamics (Tornos et al., 2015). In this study, two key phenological periods, the early and late growth stages were determined and integration was applied. Several recent studies have also highlighted the advantages of developing metrics based on variance analysis for land cover/crop classification (Maxwell and Sylvester, 2012; Mosleh and Hassan, 2014; Nuarsa et al., 2012; Qiu et al., 2015). For example, the ratio of change amplitude of LSWI to EVI2 during the period from tillering to heading dates was successfully applied for rice mapping (Qiu et al., 2015). The proposed methodology in this study has been successfully applied to the main winter wheat production areas of China during

45

the early 21st century. It is also hoped that these two phenologybased indicators could be applied to identify other agricultural crops depending on the accuracy of its estimated phenological parameters. This study hoped to contribute to this field through designing phenology-based indicators that are efficient in crop mapping. This study enriched the metric-based approaches for crop extraction. Instead of deriving indicators directly from the EVI2 temporal profiles, this study focused on specific growth stages determined by corresponding the estimated phenological dates of winter wheat from trend surfaces. It provided new thinking into developing indicators which could take account of phenological shifts and variations in the curve shapes of EVI2 temporal profiles. This study opened a new direction for dealing with the intra-class variability of VIs temporal profiles through combining variations during different phenological stages. 6. Conclusions This paper proposed a new winter wheat mapping method through exploring the variations of EVI2 temporal profiles during the estimated early and late growth stages with reference to agro-meteorological data. Two indices were designed based on the differences and variable ranges during the early and late growth stages, respectively. A simple algorithm was developed based on these two indices, since large values would be expected for winter wheat. Its application to the main winter wheat production areas of China testified its efficiency. The overall accuracy in 2012–2013 was 92.19% when compared with in-situ observations. When evaluated with Landsat-image interpreted results, the overall accuracy in 2012–2013 was 88.86%. The proposed algorithm can efficiently deal with the intra-class variability of VIs temporal profiles through exploring the reverse variation of vegetation indices before and after the estimated heading dates. Future work could involve (1) applying the proposed method to other agricultural crops through developing their own phenological date trend models; (2) employing the proposed winter wheat mapping method to time series datasets with high or medium spatial resolution such as Landsat images, Chinese Environmental Disaster Reduction Satellite images, and RapidEye images; and (3) further applications to other regions. Acknowledgements This work was supported by the National Natural Science Foundation of China (grant no. 41471362). We are very grateful for the thorough and constructive comments from the reviewers of the manuscript. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.isprsjprs.2016.09. 016. These data include Google maps of the most important areas described in this article. References Atkinson, P.M., Jeganathan, C., Dash, J., Atzberger, C., 2012. Inter-comparison of four models for smoothing satellite sensor time-series data to estimate vegetation phenology. Remote Sens. Environ. 123, 400–417. Bartholomé, E., Belward, A.S., 2005. GLC2000: a new approach to global land cover mapping from Earth observation data. Int. J. Remote Sens. 26, 1959–1977. Dong, J., Xiao, X., Kou, W., Qin, Y., Zhang, G., Li, L., Jin, C., Zhou, Y., Wang, J., Biradar, C., 2015. Tracking the dynamics of paddy rice planting area in 1986–2010 through time series Landsat images and phenology-based algorithms. Remote Sens. Environ. 160, 99–113.

46

B. Qiu et al. / ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017) 35–46

Eisavi, V., Homayouni, S., Yazdi, A., Alimohammadi, A., 2015. Land cover mapping based on random forest classification of multitemporal spectral and thermal images. Environ. Monit. Assess. 187, 1–14. Foerster, S., Kaden, K., Foerster, M., Itzerott, S., 2012. Crop type mapping using spectral–temporal profiles and phenological information. Comput. Electron. Agric. 89, 30–40. Franch, B., Vermote, E.F., Becker-Reshef, I., Claverie, M., Huang, J., Zhang, J., Justice, C., Sobrino, J.A., 2015. Improving the timeliness of winter wheat production forecast in the United States of America, Ukraine and China using MODIS data and NCAR Growing Degree Day information. Remote Sens. Environ. 161, 131– 148. Friedl, M.A., Sulla-Menashe, D., Tan, B., Schneider, A., Ramankutty, N., Sibley, A., Huang, X., 2010. MODIS Collection 5 global land cover: algorithm refinements and characterization of new datasets. Remote Sens. Environ. 114, 168–182. Gumma, M.K., Mohanty, S., Nelson, A., Arnel, R., Mohammed, I.A., Das, S.R., 2015. Remote sensing based change analysis of rice environments in Odisha, India. J. Environ. Manage. 148, 31–41. Hao, P., Zhan, Y., Wang, L., Niu, Z., Shakir, M., 2015. Feature selection of time series MODIS data for early crop classification using random forest: a case study in Kansas, USA. Remote Sens. 7, 5347. Jia, K., Wu, B., Li, Q., 2013. Crop classification using HJ satellite multispectral data in the North China Plain. J. Appl. Remote Sens. 7. 073576-073571-073512. Jiang, Z., Huete, A.R., Didan, K., Miura, T., 2008. Development of a two-band enhanced vegetation index without a blue band. Remote Sens. Environ. 112, 3833–3845. Lhermitte, S., Verbesselt, J., Verstraeten, W., Coppin, P., 2011. A comparison of time series similarity measures for classification and change detection of ecosystem dynamics. Remote Sens. Environ. 115, 3129–3152. Lunetta, R.S., Shao, Y., Ediriwickrema, J., Lyon, J.G., 2010. Monitoring agricultural cropping patterns across the Laurentian Great Lakes Basin using MODIS-NDVI data. Int. J. Appl. Earth Obs. Geoinf. 12, 81–88. Lv, Z., Liu, X., Cao, W., Zhu, Y., 2013. Climate change impacts on regional winter wheat production in main wheat production regions of China. Agric. Forest. Meteorol. 171, 234–248. Maxwell, S.K., Sylvester, K.M., 2012. Identification of ‘‘ever-cropped” land (1984– 2010) using Landsat annual maximum NDVI image composites: Southwestern Kansas case study. Remote Sens. Environ. 121, 186–195. Mosleh, M., Hassan, Q., 2014. Development of a remote sensing-based ‘‘Boro” rice mapping system. Remote Sens. 6, 1938–1953. Nuarsa, I.W., Nishio, F., Hongo, C., Mahardika, I.G., 2012. Using variance analysis of multitemporal MODIS images for rice field mapping in Bali Province, Indonesia. Int. J. Remote Sens. 33, 5402–5417. Pan, Y., Li, L., Zhang, J., Liang, S., Zhu, X., Sulla-Menashe, D., 2012. Winter wheat area estimation from MODIS-EVI time series data using the Crop Proportion Phenology Index. Remote Sens. Environ. 119, 232–242. Qin, Y., Xiao, X., Dong, J., Zhou, Y., Zhu, Z., Zhang, G., Du, G., Jin, C., Kou, W., Wang, J., Li, X., 2015. Mapping paddy rice planting area in cold temperate climate region through analysis of time series Landsat 8 (OLI), Landsat 7 (ETM+) and MODIS imagery. ISPRS J. Photogramm. Remote Sens. 105, 220–233. Qiu, B., Li, W., Tang, Z., Chen, C., Qi, W., 2015. Mapping paddy rice areas based on vegetation phenology and surface moisture conditions. Ecol. Indic. 56, 79–86.

Qiu, B.W., Fan, Z.L., Zhong, M., Tang, Z.H., Chen, C.C., 2014. A new approach for crop identification with wavelet variance and JM distance. Environ. Monit. Assess. 186, 7929–7940. Ren, J., Chen, Z., Zhou, Q., Tang, H., 2008. Regional yield estimation for winter wheat with MODIS-NDVI data in Shandong, China. Int. J. Appl. Earth Obs. Geoinf. 10, 403–413. Sakamoto, T., Wardlow, B.D., Gitelson, A.A., Verma, S.B., Suyker, A.E., Arkebauer, T.J., 2010. A two-step filtering approach for detecting maize and soybean phenology with time-series MODIS data. Remote Sens. Environ. 114, 2146–2159. Siebert, S., Ewert, F., 2012. Spatio-temporal patterns of phenological development in Germany in relation to temperature and day length. Agric. Forest. Meteorol. 152, 44–57. Sun, H., Xu, A., Lin, H., Zhang, L., Mei, Y., 2012. Winter wheat mapping using temporal signatures of MODIS vegetation index data. Int. J. Remote Sens. 33, 5026–5042. Tornos, L., Huesca, M., Dominguez, J.A., Moyano, M.C., Cicuendez, V., Recuero, L., Palacios-Orueta, A., 2015. Assessment of MODIS spectral indices for determining rice paddy agricultural practices and hydroperiod. ISPRS J. Photogramm. Remote Sens. 101, 110–124. Upadhyay, P., Ghosh, S.K., Kumar, A., 2015. Temporal MODIS data for identification of wheat crop using noise clustering soft classification approach. Geocarto Int. 31, 1–18. Wang, X., Li, X., Tan, M., Xin, L., 2015. Remote sensing monitoring of changes in winter wheat area in North China Plain from 2001 to 2011. Trans. Chin. Soc. Agric. Eng. 31, 190–199. Wardlow, B.D., Egbert, S.L., 2010. A comparison of MODIS 250-m EVI and NDVI data for crop mapping: a case study for southwest Kansas. Int. J. Remote Sens. 31, 805–830. Wardlow, B.D., Egbert, S.L., Kastens, J.H., 2007. Analysis of time-series MODIS 250 m vegetation index data for crop classification in the US Central Great Plains. Remote Sens. Environ. 108, 290–310. Xiao, D., Tao, F., Liu, Y., Shi, W., Wang, M., Liu, F., Zhang, S., Zhu, Z., 2013. Observed changes in winter wheat phenology in the North China Plain for 1981–2009. Int. J. Biometeorol. 57, 275–285. Zhang, J., Feng, L., Yao, F., 2014. Improved maize cultivated area estimation over a large scale combining MODIS–EVI time series data and crop phenological information. ISPRS J. Photogramm. Remote Sens. 94, 102–113. Zhang, L., Weng, Q., 2016. Annual dynamics of impervious surface in the Pearl River Delta, China, from 1988 to 2013, using time series Landsat imagery. ISPRS J. Photogramm. Remote Sens. 113, 86–96. Zhang, M.-Q., Guo, H.-Q., Xie, X., Zhang, T.-T., Ouyang, Z.-T., Zhao, B., 2013. Identification of land-cover characteristics using MODIS time series data: an application in the Yangtze River Estuary. PLoS ONE 8, e70079. Zhong, L., Hawkins, T., Biging, G., Gong, P., 2011. A phenology-based approach to map crop types in the San Joaquin Valley, California. Int. J. Remote Sens. 32, 7777–7804. Zhou, F., Zhang, A., Townley-Smith, L., 2013. A data mining approach for evaluation of optimal time-series of MODIS data for land cover mapping at a regional level. ISPRS J. Photogramm. Remote Sens. 84, 114–129.