[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Next Article in Journal
Spatio-Temporal Assessment of Tuz Gölü, Turkey as a Potential Radiometric Vicarious Calibration Site
Previous Article in Journal
Response to Sheppard, C. Atoll Rim Expansion or Erosion in Diego Garcia Atoll, Indian Ocean? Comment on Hamylton, S.; East, H. A Geospatial Appraisal of Ecological and Geomorphic Change on Diego Garcia Atoll, Chagos Islands (British Indian Ocean Territory). Remote Sens. 2012, 4, 3444–3461
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 6
Issue 3
10.3390/rs6032473
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mapping Crop Cycles in China Using MODIS-EVI Time Series

by
Le Li
1,2,*,
Mark A. Friedl
2,*,
Qinchuan Xin
3,
Josh Gray
2,
Yaozhong Pan
4 and
Steve Frolking
5
1
State Key Laboratory of Vegetation and Environmental Change, Institute of Botany, Chinese Academy of Sciences, Beijing 100093, China
2
Department of Earth and Environment, Boston University, Boston, MA 02215, USA
3
Ministry of Education Key Laboratory for Earth System Modeling, Tsinghua University, Beijing 100084, China
4
State Key Laboratory of Earth Processes and Resource Ecology, Beijing Normal University, Beijing 100875, China
5
Institute for the Study of Earth, Oceans and Space, University of New Hampshire, Durham, NH 03824, USA
*
Authors to whom correspondence should be addressed.
Remote Sens. 2014, 6(3), 2473-2493; https://doi.org/10.3390/rs6032473
Submission received: 30 December 2013 / Revised: 6 March 2014 / Accepted: 10 March 2014 / Published: 20 March 2014
Browse Figures

<p>The climate zones in China as derived from Peel <span class="html-italic">et al.</span> [<a href="#b42-remotesensing-06-02473" class="html-bibr">42</a>]. The areas of South China Sea Islands are not displayed; there is little cropland in these islands.</p> ">
<p>Flowchart showing the algorithm used for mapping agricultural intensity.</p> ">
<p>15-month time series fits for MODIS EVI data shown for pixels with (<b>A</b>) single cropping, (<b>B</b>) double cropping, (<b>C</b>) triple cropping, and (<b>D</b>) frequent clouds. Savitzky-Golay (SG), asymmetric Gaussian (AG), and double logistic (DL) smoothing functions were applied in TIMESAT, and the adaptive Savitzky-Golay filter was chosen for all subsequent analyses.</p> ">
<p>The agricultural mask for mainland China as derived from the IGBP cropland layer in MODIS Land Cover Type product. Both cropland (class 12) and cropland/natural vegetation mosaic (class 14) are used to define the agricultural areas in this study.</p> ">
<p>Agricultural intensity maps for China shown for (<b>A</b>) 2006 and (<b>B</b>) 2007.</p> ">
<p>Regional-scale views of agricultural intensity in China in 2006 for major agricultural regions: (<b>A</b>) Northeast China Plain, (<b>B</b>) North China Plain and Middle-lower Yangzte Plain, (<b>C</b>) Sichuan Basin, and (<b>D</b>) Pearl River Delta (see annotations in <a href="#f1-remotesensing-06-02473" class="html-fig">Figure 1</a>). Heilongjiang province (<a href="#f6-remotesensing-06-02473" class="html-fig">Figure 6A</a>) has the most arable land in China and Henan (<a href="#f6-remotesensing-06-02473" class="html-fig">Figure 6B</a>) is China’s largest grain-producing province.</p> ">
<p>Comparison of (<b>A</b>) arable areas and (<b>B</b>) gross sown areas between estimates from the MODIS Land Cover Type product (MCD12Q1) and national survey data reported by National Bureau of Statistics of China (NBSC) at the province level in 2006. Values for R<sup>2</sup> (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.</p> ">
<p>Comparisons of gross sown areas in (<b>A</b>) 2006 and (<b>B</b>) 2007 for provinces in China estimated from the agricultural intensity map with national survey data from China Statistical Yearbook. Values for R<sup>2</sup> (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.</p> ">
<p>Comparisons of the gross-sown areas for prefectures in Henan province between estimates from the agricultural intensity map and survey data from Henan Statistical Yearbook. Values for R<sup>2</sup> (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.</p> ">
Versions Notes

Abstract

:
As the Earth’s population continues to grow and demand for food increases, the need for improved and timely information related to the properties and dynamics of global agricultural systems is becoming increasingly important. Global land cover maps derived from satellite data provide indispensable information regarding the geographic distribution and areal extent of global croplands. However, land use information, such as cropping intensity (defined here as the number of cropping cycles per year), is not routinely available over large areas because mapping this information from remote sensing is challenging. In this study, we present a simple but efficient algorithm for automated mapping of cropping intensity based on data from NASA’s (NASA: The National Aeronautics and Space Administration) MODerate Resolution Imaging Spectroradiometer (MODIS). The proposed algorithm first applies an adaptive Savitzky-Golay filter to smooth Enhanced Vegetation Index (EVI) time series derived from MODIS surface reflectance data. It then uses an iterative moving-window methodology to identify cropping cycles from the smoothed EVI time series. Comparison of results from our algorithm with national survey data at both the provincial and prefectural level in China show that the algorithm provides estimates of gross sown area that agree well with inventory data. Accuracy assessment comparing visually interpreted time series with algorithm results for a random sample of agricultural areas in China indicates an overall accuracy of 91.0% for three classes defined based on the number of cycles observed in EVI time series. The algorithm therefore appears to provide a straightforward and efficient method for mapping cropping intensity from MODIS time series data.

1. Introduction

Information regarding the spatial distribution of land cover types is essential for a wide range of applications including vegetation monitoring, natural resource management, environmental protection, and climate change studies [1–4]. Remotely sensed data from the MODerate Resolution Imaging Spectroradiometer (MODIS) onboard NASA’s Terra and Aqua satellites are well-suited for global land-cover and land-use monitoring and mapping activities because they provide data with high temporal and moderate spatial resolution [5–8]. The global MODIS Land Cover Type product (MCD12) [9,10], for example, is a critical input to many global climate, vegetation, and atmospheric transport models [11–13].
In addition to its importance for modeling, regional-to-global land cover and land use information is essential for studies focused on interactions between humans and natural systems [14–17]. Cropping intensity, which we define here as the number of cropping cycles per year, is an important dimension of land use that is strongly influences water demand and agricultural production [18–21], but has received relatively little attention. In areas, such as Asia, which have limited lands available for arable expansion, crop production is commonly increased by planting crops more than once a year in the same field [22–25]. As such, the gross sown area (the total area sown with crops in a year, accounting for multiple crop cycles) in many parts of Asia is frequently larger than the corresponding amount of arable land. To date, however, information regarding how the amount of arable land differs from total harvested or sown area arising from multi-cropping and fallow rotations has not been provided by global land cover products from instruments such as MODIS. As agricultural intensification has significant societal and environmental implications and is poorly characterized and understood [26], improved information regarding the geographic distribution and dynamics of cropping intensity is essential in regions of the world where multi-cropping is prevalent.
The primary technical challenge in mapping cropping intensity from remotely sensed data is to identify crop phenological cycles from time series of satellite-derived vegetation indices, such as the Enhanced Vegetation Index (EVI) or Normalized Difference Vegetation Index (NDVI) [3,27–29]. Specifically, areas with multiple cropping should exhibit cyclic patterns in time series of MODIS vegetation indices that reflect cropping cycles on the ground and distinguish them from areas with single cropping (crops planted once in a year) or natural vegetation [30–32]. For this work, we use MODIS EVI time series because the EVI has a greater dynamic range than NDVI and is therefore better suited to capturing phenological variation in croplands with large seasonal amplitudes in leaf area [33].
A number of previous studies have attempted to map cropping intensity in China [24,34,35]. However, operational mapping of agricultural intensity based on MODIS data remains a challenge for a variety of reasons. First, trajectories of MODIS EVI time series in croplands are highly variable over large areas [36,37]. For example, double cropping systems of wheat-maize or wheat-soybean dominate the North China Plain, while triple cropping systems, such as rice-rice-rice, wheat-rice-rice, or rape-rice-rice, are common in Southern China [22]. Many of these areas are managed in rotation systems, and planting dates are highly variable from year to year. Second, MODIS pixels frequently include a mixture of land-cover types that introduce scale-dependent errors in land cover maps [38]. Hence, separating multiple cropping from single cropping is especially difficult when both cropping strategies are located in the same pixel. Third, noise from cloud contamination, variable viewing and illumination angles, resampling procedures [39,40], or persistent missing data due to clouds in humid climates often introduce large amounts of noise and missing data in MODIS EVI time series. Finally, areas planted with winter wheat usually manifest double peaks (one before dormancy with a second in the spring) each year [1,41]. For all of the above-mentioned reasons, conventional remote sensing classification algorithms may not provide good results.
The objective of the research described in this paper is to use MODIS EVI time series to derive multi-cropping information that complements and augments information related to agricultural extent provided by the MODIS land cover product [10]. To achieve this goal we developed a straightforward and efficient algorithm for mapping multi-cropping practices from MODIS time series data, and evaluated the results using a random sample of MODIS pixels, national survey data, and field-scale samples.

2. Study Area and Materials

2.1. Description of the Study Area

Mainland China (Figure 1), our study area, presents a challenging environment for mapping agricultural intensity because many croplands in China have small field sizes. Further, the geographic scale and diversity of climatic regimes leads to different cultivation habits throughout China. According to the Koppen-Geiger climate classification scheme, most croplands in China are located in continental hot summer, temperate dry winter hot summer, and temperate humid with hot summer climate zones (Figure 1) [42]. Cropping intensity increases gradually from one to three harvests per year from Northern China to Southern China [2].

2.2. MODIS Data

MODIS land products available from the USGS EROS Data Center (USGS: United States Geological Survey; EROS: Earth Resources Observation and Science) [43] are organized in a tiling system using a Sinusoidal projection. Eighteen MODIS tiles are required to cover our study area encompassing mainland China. For this work we used eight years of data (2005–2012) for multiple MODIS land products to map agricultural intensity. The primary inputs to our algorithm were the Terra and Aqua MODIS eight-day 500 m surface reflectance products (MOD09A1 and MYD09A1). Atmospheric corrections are implemented in the production of MO[Y]D09A1, but bidirectional reflectance distribution function (BRDF) effectsare not fully accounted for [44]. We also used the accompanying quality control flags to screen for low quality observations.
Additionally, the MODIS Land Surface Temperature (LST) product (MOD11A2) [45] and the MODIS Collection 5 Land Cover Type product (MCD12Q1) [10] were used to refine our mapping results. The eight-day composite MOD11A2 product provides daytime/nighttime LST at 1 km resolution (resampled to 500 m to match other products in our study). The MODIS land cover maps are produced for each calendar year for five different classification schemes; for this work we used MODIS land cover mapped according to the International Geosphere Biosphere Program (IGBP) classfication scheme [46]. More details about the MODIS products are available on the USGS website [47].

3. Methods

The workflow to map agricultural intensity is presented schematically in Figure 2. The algorithm consists of three main steps: (1) Time-series smoothing of MODIS EVI data, (2) identification of phenology cycles from the fitted EVI time series, and (3) incorporating ancillary MODIS data to produce a final map of agricultural intensity.

3.1. Preprocessing of MODIS Surface Reflectance Data

MODIS surface reflectance data from the Terra and Aqua satellites are composited for 8 day periods based on the “minimum blue compositing method” [40,48]. As land surface reflectance data are sensitive to aerosol contamination in the blue waveband, this compositing method is designed to select the clearest conditions over the compositing period [44,49]. The use of both Terra and Aqua observations minimizes the influence of non-ideal observations such as clouds, which are frequent in Southern China. The cloud mask defined in the quality control data is used to eliminate poor-quality pixels; specifically, a cloud flag is assigned if all Terra and Aqua observations in the same 8-day compositing period are cloudy.
The composited MODIS surface reflectance data are used to compute MODIS EVI time series:
EVI = G × ρ nir − ρ red ρ nir + C 1 × ρ red − C 2 × ρ blue + L
where ρnir, ρred, and ρblue are surface reflectance in the near-infrared, red, and blue bands, respectively, and L = 1, C1 = 6, C2 = 7.5, and G (gain factor) = 2.5 [33].
We chose to derive MODIS EVI time series from MODIS surface reflectance products instead of using the available standard VI products. Specifically, the MODIS Vegetation Indices product (MOD13) [33] is sub-optimal for mapping agricultural intensity because the compositing period (16-days) is insufficient to capture rapid temporal variations in crop phenology.

3.2. Time-Series Smoothing of MODIS EVI Data

Time-series smoothing of MODIS EVI data is required to reduce noise unrelated to vegetation development or management practices [30,50,51]. For this work, we used the open-source software TIMESAT [52] to fit MODIS EVI time-series data because it offers flexibility for this application. In TIMESAT, three kinds of smoothing functions are available (asymmetric Gaussian, double logistic, and adaptive Savitzky-Golay) with user-defined weighting schemes. To minimize the influence of clouds, observations with cloud flags derived from MODIS quality control data were assigned lower weights in our time-series smoothing.
Figure 3 shows examples of MODIS EVI time series with three TIMESAT smoothing functions fitted. While the three smoothing methods have similar performance for single-cropped pixels, results vary for multi-cropped pixels (Figure 3B,C). Specifically, the EVI trajectory associated with multicropping activities tends to be over-smoothed when asymmetric Gaussian or double logistic smoothing functions are applied. The adaptive Savitzky-Golayfilter provides a better fit, particularly in pixels where multicropping is present. In the adaptive Savitzky-Golay filter, the window size controls the number of successive observations and the polynomial degree controls the order of functions for time-series smoothing. We set a window size of 4 for 8-day composite data (i.e., 32-days) and used second order polynomials to preserve the subtle information related to variation in crop phenology.
Our analysis suggests that the MODIS cloud product is conservative because some clear-sky observations are flagged as cloudy (Figure 3B,C). For cloudy areas, such as Southern China, high frequency and persistent missing data present substantial difficulties for capturing phenological variation in crops (Figure 3D). Successfully smoothed time series is therefore both important and challenging when continuous observations during the growing season are missing due to clouds.

3.3. Identifying Phenological Cycles Based on Smoothed MODIS EVI Time Series

We employed TIMESAT to smooth the EVI time series using an adaptive Savitzky-Golay filter and developed our own approach to identifying cropping cycles from the resulting EVI time series. TIMESAT can also identify distinct growth cycles at each pixel based on the amplitude between the secondary maximum and the primary maximum of time-series data [53]. However, because the current version of TIMESAT only allows a maximum of two growth cycles each year, it does not accommodate triple cropped areas.
Our approach to identify phenological cycles from the temporal trajectories of the fitted MODIS EVI time series includes three steps. First, a moving window is applied to smoothed MODIS EVI time series. Potential peaks (local maxima) and troughs (local minima) are flagged when the central MODIS observation in the window has the highest or lowest EVI value. Selecting a proper window size is important because selecting a window that is too small can result in too few peaks and troughs, and vice versa. Based on tests applied to a wide range of time series, we chose a window size of nine 8-day composites, meaning that peaks/troughs reflect maxima/minima within 72-day moving windows. Second, because pixels with sparse vegetation may also have peaks and troughs that are unrelated to vegetation development and EVI values in croplands generally exceed 0.4 at peak growth [30], spurious peaks were discarded if the corresponding EVI values were less than 0.35. Third, potential peaks and troughs were checked to make sure that only one trough is present between peaks; where present, successive potential peaks with no intervening trough were merged. This process reduces the influence of occasional undetected clouds in MODIS EVI time series.

3.4. Mapping Agricultural Intensity by Incorporating Ancillary MODIS Data

Refinement of detected phenological cycles is necessary to avoid errors caused by a variety of factors. In particular, EVI values over winter wheat usually exhibit an additional peak before cold acclimation [54], so simply counting detected phenological cycles can result in overestimation of crop cycles. Moreover, the fitted curve may not be accurate when there are large data gaps in the time series. Gap filling long series of consecutive missing values can therefore lead to under-prediction when entire crop cycles are missed.
To reduce overestimation of multi-cropping in regions with winter wheat, we use the 8-day MODIS LST product to refine potential peaks based on assumptions regarding ecologically realistic temperatures. Other phenology studies have employed similar techniques [3,27,29], because vegetation processes are sensitive to low temperatures and vegetation photosynthesis increases rapidly once nighttime surface temperature exceeds a minimum threshold (usually −2∼5 °C) [55,56]. Unlike natural vegetation, annual crops such as maize and soybean are often planted when soil temperatures are warm enough for germination (above 8–11 °C) and are harvested four to twelve weeks after peak canopy development [57]. The EVI of paddy rice fields increases after planting, peaks around the heading period, and decreases as leaves wither before harvesting [58]. Winter wheat has a unique phenology cycle: it requires biological activity for three to four weeks in autumn, enters a state of dormancy during the winter, and emerges from dormancy in spring when temperatures exceed about 9 °C [1,54]. Here, we apply a conservative minimum threshold of 5 °C to nighttime LST to define the thermal growing season of crops based on analysis of the four major crops in China [59]. To avoid counting the autumn growth of winter wheat as a separate crop cycle, the period for peak EVI was constrained to occur between the start and three weeks before the end of the calendar-year thermal growing season. Finally, acknowledging the impacts of persistent cloud cover, we also identify pixels with large data gaps in their EVI time series. If four or more successive MODIS observations during the growing season are missed due to clouds, the pixels are flagged accordingly.
Using this approach, the number of phenology cycles can be determined by counting the number of seasonal peaks each year. Pixels with more than three seasons detected in a single year are assigned to have three cycles, since four or more cropping cycles per year is rare [28]. Finally, to produce the final map of agricultural intensity, our algorithm uses an agricultural land cover mask derived from the MODIS Land Cover product (Figure 4), which includes all pixels belonging to “Croplands” (class 12) and “Croplands/Natural vegetation mosaic” (class 14) in the IGBP classification scheme (Type 1 in MCD12Q1).

3.5. Accuracy Assessment

Accuracy assessment was performed at three scales of spatial aggregation. At the provincial and the prefectural level, gross sown areas derived from MODIS were compared with national survey data. To do this, we used nationwide survey data that was collected using a stratified sampling method that included 130,000 plots from 20,000 villages (note, however, that while provincial data is available for all of China, prefectural data is only available for some provinces because of data restriction policies). These data provide the arable area and gross sown area for provinces, and are derived from China Statistical Yearbooks for 2006–2012 (China Statistical Yearbook, 2006; English version is available on the website of National Bureau of Statistics of China [60]). Although there are well-known problems with these data [61,62], the national survey data represent the best available source of land use reference data in China. The gross sown areas for prefectures in Henan (China’s largest crop-producing province) were derived from the Henan Statistical Year books for 2007–2008 (Chinese version with English annotations).
To derive the amount arable area within prefectures and provinces from MODIS, we tabulated areas for agriculture and agriculture mosaic classes from the MODIS land cover product and then scaled them to account for the fact that the agriculture and agricultural mosaic classes in the MODIS Land Cover Type product are not defined to include 100% agriculture. Specifically, croplands (IGBP class 12) are defined to be between 60% and 100% agriculture by area and cropland/natural vegetation mosaics (IGBP class 14) are defined to be 30%–60% agriculture by area. To account for this, we used fractional weights of 0.7 and 0.4 to adjust the final area of agriculture assigned to pixels in class 12 and class 14, respectively. While these proportions clearly vary from pixel to pixel, comparison of estimated arable area with national survey data at the provincial level show that at the scale of provinces, these adjustments successfully account for sub-pixel mixtures of agriculture with other land cover classes. In a similar fashion, gross sown areas at the prefectural and provincial level from MODIS were counted twice (three times) for double (triple) cropping in our 500-m maps of agricultural intensity.
At the pixel level, we assessed our mapping accuracy for results from 2006 using 4500 individual sample points allocated throughout croplands in mainland China (Appendix). These sample points were randomly selected based on a stratified random sampling method with 1,500 samples allocated to each of three strata: single, double, and triple cropping. Independent professionals, with no knowledge of the mapping results, manually labeled the number of crop cycles for each selected pixel based on MODIS EVI time series (see Figure 3 for examples) and high-resolution imagery from Google Earth. Interpreters also provided commentary flags to indicate whether the EVI trajectory at each sample pixel was too noisy to determine the number cropping cycles. To assess the final map, we computed the producer’s accuracy, user’s accuracy, and overall accuracy, which are widely used accuracy metrics [9,63,64].

4. Results

We mapped annual agricultural intensity in mainland China at 500 m spatial resolution for 2005–2012. Figure 5 shows maps of agricultural intensity in 2006 and 2007 as examples and illustrates two important results. First, the overall pattern of multi-cropping is consistent with expectations and national inventory data (Table 1). Second, the pattern of cropping systems is stable from year to year. Some areas in the North China Plain or the Sichuan Basin may switch between on-year-one-harvest and one-year-two-harvest due to crop rotations.
Regional-scale views of agricultural intensity for 2006 are shown in Figure 6 for several major agricultural regions in China, and a summary of cropland area from China Statistical Yearbook in China in 2006 is shown in Table 1. The spatial distribution of agricultural intensity is consistent with local knowledge and varies with climate and geographic features. In the Northeast China Plain (Figure 6A), where the temperatures are relatively cold in winter, annual crops, such as corn or rice, are planted once per year. By comparison, double cropping systems are widely adopted in the North China Plain (Figure 6B) by planting winter wheat in combination with maize or crops with similar phenology cycles, and which are planted directly following winter wheat harvest in late spring or early summer. The Sichuan Basin (Figure 6C) is a major agricultural region with complex cropping systems that depend on freshwater availability and temperature regimes. In the Middle-lower Yangtze Plain (Figure 6B) and the Pearl River Delta (Figure 6D), where the climate is hot and there is ample freshwater, paddy rice is planted two or three times a year. This region also includes areas with extensive data gaps because of persistent clouds (e.g., some areas in the Pearl River Delta or Sichuan Basin).
At the provincial level, we estimated the total agricultural area in each province from the MODIS Land Cover Type product and compared the results to national survey data in 2006. While the derived agricultural areas match well with the reported arable areas (Figure 7), the MODIS land cover product provides incomplete information related to agricultural land use (gross sown area), which incorporates multi-cropping. In Inner Mongolia, for example, gross sown areas are overestimated because off allow rotations; conversely, gross sown area is under estimated in top agricultural provinces, such as Henan, due to multiple croppings.
Multi-cropping information from MODIS successfully captures patterns in gross sown area. At the provincial level gross sown areas derived from MODIS agree well with national survey data (Figure 8). Specifically, agricultural intensity maps from MODIS capture variability in inventory-based measures of gross sown areas that are not captured in the land cover map alone (R2 = 0.92 compared to 0.73 in 2006). In this context, it is worth noting that very high agricultural intensity makes Henan province the largest crop-producing province in China, even though Heilongjiang province has much more arable area. As a result, underestimation of gross sown area in Henan is reduced when double cropping is taken into account. Similarly, overestimation of gross sown area in Inner Mongolia based on MCD12Q1 is significantly reduced because fallow rotations and sparse agricultural areas with low EVI values are excluded. Extending this analysis to include five years of data, we find that the relationship between our results and national survery data is consistent across years (Table 2), with R2 values that are consistently high (0.89–0.92) and low (but positive) mean errors (ME) that suggest a small but consistent overestimation of 300–600 kha per province.
Assessment of results at finer spatial scale (prefectural level) for Henan province show similar agreement to results obtained with national survey data, with R2 values exceeding 0.97 for 2006 and 2007 (Figure 9). Unfortunately, data policies in China prevent us from obtaining prefectural-level data for other provinces. However, because agriculture is both extensive and intensive in Henan, the province provides an excellent test region for this analysis and the prefectoral-level results we obtained provide confidence that our results are robust and extensible to other provinces.
Finally, to assess the effectiveness of our multi-cropping algorithm, we compared our results with visual assessments based on manual interpretation of MODIS time series. Table 3 presents an error matrix showing results from this analysis, and indicates that the overall accuracy of multi-cropping classes is 91.0% for the 2006 version of China’s agricultural intensity map. Producer and user accuracies are higher than 86% for the three agricultural classes (single-, double-, and triple-cropping). Only one non-cropping sample point was erroneously classified as single cropping, apparently because of over fitting by our algorithm to the EVI time series. Confusion between single- and double-cropping, and between double- and triple-cropping cause the most errors of commission. Overall, the results are very encouraging; however, further evaluation for our algorithm requires field-level ground reference data for cropping intensity.

5. Discussion

5.1. Factors That Influence the Mapping Accuracy

Several factors influence the accuracy of agricultural intensity maps produced using our method. First, while our algorithm screens for good-quality observations, data gaps due to clouds create difficulties identifying crop cycles. Even though the smoothing process can reduce the influence of intermittent cloud-cover, cropping intensity is probably under estimated in areas with persistent cloud-cover, such as Southern China. To address this, we implemented two quality controls measures: (1) The time-series smoothing process incorporates MODIS quality control data, and (2) a quality control layer is generated to indicate if large data gaps were present during the growing season for the year of interest. Similarly, snow cover at high latitudes reduces the number of good observations available from MODIS in the winter. However, this constraint is much less significant than cloud cover in the south because wintertime EVI values are excluded based on the thermal growing season derived from the MODIS LST product.
The MODIS Land Cover Type product is used in this study to distinguish croplands from natural vegetation, and, thus, significantly determines the accuracy of our maps. Hence, errors in the MCD12Q1 product obviously propagate into our results. However, the results we obtain through this work suggest that the MCD12Q1 provides a realistic representation of agriculture in this area. Indeed, a key goal of this work is to demonstrate that cropping intensity maps derived from MODIS EVI time series provide unique and complementary information to conventional land cover products such as MCD12Q1.
A related source of uncertainty in our results is subpixel spatial heterogeneity. To exclude areas that do not show strong seasonal cycles (and hence are unlikely to include crops), our algorithm requires that the magnitude of EVI peak for any detected cycle exceed 0.35. However, this method may bias our results because field sizes in China vary widely and many (if not most) cropland areas are mixed at the spatial resolution of MODIS. Previous studies have reported field sizes of 0.01 to 2 ha in Central and Eastern China [38,66,67], which are well below the spatial resolution of MODIS 500 m pixels (6.25 ha). Our algorithm does not account for this scale of sub-pixel variation, which means that results will be most reliable for areas with large field sizes. The effects of spatial heterogeneity is a common challenge for MODIS-based studies of land use and land cover change because most land use and land cover conversions occur below the spatial resolution of MODIS [68–70].

5.2. Potential Refinements

The results from our study demonstrate the potential of our proposed algorithm for mapping cropping intensity over large geographic regions from MODIS data. However, several refinements might improve our algorithm and make it more robust for operational applications outside of China. From a practical perspective, the availability of reliable ground reference data is a primary constraint on algorithm development and accuracy assessment for this problem domain. National survey data in China and other countries with high agricultural intensity only provide coarse-scale information that cannot be used to evaluate results at the scale of individual pixels. While field-based studies that collect data in-situ may mitigate this challenge, scaling issues (e.g., whether the footprints of field samples and MODIS pixels match each other) can complicate assessment efforts [71,72]. Further, in-situ data are difficult and expensive to obtain, and cannot be collected over large areas without very high costs. For this work, we interpreted MODIS EVI time-series manually for individual pixels. However, it is often difficult to reliably identify cropping cycles in this way, especially when pixels are mixed or when the trajectories are contaminated with noise. Therefore, there is a critical need to develop reliable reference datasets for future research efforts.
Finally, and as we discussed above, subpixel spatial heterogeneity is a significant challenge. Some studies have attempted to address this using time-series decomposition [37] or machine-learning algorithms [73,74]. However, when different cropping system are mixed in the same pixel, MODIS EVI time series are often quite noisy, making it difficult to identify cropping cycles, even by manual interpretation. Time-series decomposition and machine learning methods may be most useful for deriving subpixel information in areas planted with annual crops, but both methods are relatively untested in areas with multi-cropping.

6. Conclusions

We present a straight forward and efficient algorithm to map agricultural intensity over large geographic regions using eight-day MODIS EVI time series data derived from Terra and Aqua MODIS surface reflectance products. Time series are gap-filled and smoothed using an adaptive Savitzky-Golay filter. To detect annual cropping cycles for individual pixels, the algorithm iteratively applies two moving windows to the fitted EVI time series. Results from this process are then refined using the MODIS LST product to define the thermal growing season at each pixel, and using the MODIS land cover product to determine agricultural areas.
Results from applying this method provide agricultural intensity maps over multiple years that are consistent with known geographic patterns of multi-cropping throughout China. Comparison of our results against national survey data shows substantial complementary information related to estimated gross sown areas at provincial and prefectural levels compared to using MODIS Land Cover data alone. Accuracy assessment using data generated by expert classification of randomly selected pixel samples indicates an overall accuracy of 91.0% for three agricultural intensity classes. More generally, the algorithm shows significant potential to automatically estimate reliable cropping intensity information in support of large-scale studies of agricultural land use and land cover dynamics.

Acknowledgments

The authors gratefully acknowledge funding support for this research provided by the National Aeronautics and Space Administration (grant number NNX11AE75G) and by the National Science Foundation (grant number EAR-1038907& EAR-1038818). This study was partially funded by the China Scholarship Council and China Postdoctoral Science Foundation (grant number 2013M540087). We thank Professor Lars Eklundh at Lund University for suggestions on TIMESAT. We also thank Damien Sulla-Menashe and Xiaoman Huang at Boston University for their data processing on MODIS land cover maps and helpful discussions.

Author Contributions

Le Li, Mark A Friedl, and Qinchuan Xin designed the research and performed data analysis. Josh Gray and Steve Frolking made major contributions to their interpretations. Yaozhong Pan provided the Chinese data compilations and helped collecting statistic data. All authors contributed with ideas, writing and discussions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Pan, Y.; Li, L.; Zhang, J.; Liang, S.; Zhu, X.; Sulla-Menashe, D. Winter wheat area estimation from MODIS-EVI time series data using the Crop Proportion Phenology Index. Remote Sens. Environ 2012, 119, 232–242. [Google Scholar]
  2. Piao, S.; Ciais, P.; Huang, Y.; Shen, Z.; Peng, S.; Li, J.; Zhou, L.; Liu, H.; Ma, Y.; Ding, Y. The impacts of climate change on water resources and agriculture in China. Nature 2010, 467, 43–51. [Google Scholar]
  3. Zhang, X.; Friedl, M.A.; Schaaf, C.B.; Strahler, A.H.; Hodges, J.C.; Gao, F.; Reed, B.C.; Huete, A. Monitoring vegetation phenology using MODIS. Remote Sens. Environ 2003, 84, 471–475. [Google Scholar]
  4. Atzberger, C. Advances in remote sensing of agriculture: Context description, existing operational monitoring systems and major information needs. Remote Sens 2013, 5, 949–981. [Google Scholar]
  5. Justice, C.; Townshend, J.; Vermote, E.; Masuoka, E.; Wolfe, R.; Saleous, N.; Roy, D.; Morisette, J. An overview of MODIS land data processing and product status. Remote Sens. Environ 2002, 83, 3–15. [Google Scholar]
  6. Muchoney, D.; Borak, J.; Chi, H.; Friedl, M.; Gopal, S.; Hodges, J.; Morrow, N.; Strahler, A. Application of the MODIS global supervised classification model to vegetation and land cover mapping of central America. Int. J. Remote Sens 2000, 21, 1115–1138. [Google Scholar]
  7. Myneni, R.; Hoffman, S.; Knyazikhin, Y.; Privette, J.; Glassy, J.; Tian, Y.; Wang, Y.; Song, X.; Zhang, Y.; Smith, G. Global products of vegetation leaf area and fraction absorbed par from year one of MODIS data. Remote Sens. Environ 2002, 83, 214–231. [Google Scholar]
  8. Pittman, K.; Hansen, M.C.; Becker-Reshef, I.; Potapov, P.V.; Justice, C.O. Estimating global cropland extent with multi-year MODIS data. Remote Sens 2010, 2, 1844–1863. [Google Scholar]
  9. Friedl, M.A.; McIver, D.K.; Hodges, J.C.; Zhang, X.; Muchoney, D.; Strahler, A.H.; Woodcock, C.E.; Gopal, S.; Schneider, A.; Cooper, A. Global land cover mapping from MODIS: Algorithms and early results. Remote Sens. Environ 2002, 83, 287–302. [Google Scholar]
  10. Friedl, M.A.; Sulla-Menashe, D.; Tan, B.; Schneider, A.; Ramankutty, N.; Sibley, A.; Huang, X. MODIS collection 5 global land cover: Algorithm refinements and characterization of new datasets. Remote Sens. Environ 2010, 114, 168–182. [Google Scholar]
  11. Lim, Y.K.; Cai, M.; Kalnay, E.; Zhou, L. Observational evidence of sensitivity of surface climate changes to land types and urbanization. Geophys. Res. Lett 2005, 32. [Google Scholar] [CrossRef]
  12. Running, S.W.; Nemani, R.R.; Heinsch, F.A.; Zhao, M.; Reeves, M.; Hashimoto, H. A continuous satellite-derived measure of global terrestrial primary production. Bioscience 2004, 54, 547–560. [Google Scholar]
  13. Turner, D.P.; Ritts, W.D.; Cohen, W.B.; Gower, S.T.; Zhao, M.; Running, S.W.; Wofsy, S.C.; Urbanski, S.; Dunn, A.L.; Munger, J. Scaling gross primary production (GPP) over boreal and deciduous forest landscapes in support of MODIS GPP product validation. Remote Sens. Environ 2003, 88, 256–270. [Google Scholar]
  14. Monfreda, C.; Ramankutty, N.; Foley, J.A. Farming the planet: 2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000. Global Biogeochem. Cycle 2008, 22. [Google Scholar] [CrossRef]
  15. Ramankutty, N.; Evan, A.T.; Monfreda, C.; Foley, J.A. Farming the planet: 1. Geographic distribution of global agricultural lands in the year 2000. Global Biogeochem. Cycles 2008, 22. [Google Scholar] [CrossRef]
  16. Pervez, M.S.; Brown, J.F. Mapping irrigated lands at 250-m scale by merging MODIS data and national agricultural statistics. Remote Sens 2010, 2, 2388–2412. [Google Scholar]
  17. You, X.; Meng, J.; Zhang, M.; Dong, T. Remote sensing based detection of crop phenology for agricultural zones in China using a new threshold method. Remote Sens 2013, 5, 3190–3211. [Google Scholar]
  18. Doll, P.; Siebert, S. Global modeling of irrigation water requirements. Water Resour. Res 2002, 38. [Google Scholar] [CrossRef]
  19. Foley, J.A.; DeFries, R.; Asner, G.P.; Barford, C.; Bonan, G.; Carpenter, S.R.; Chapin, F.S.; Coe, M.T.; Daily, G.C.; Gibbs, H.K. Global consequences of land use. Science 2005, 309, 570–574. [Google Scholar]
  20. Tilman, D.; Fargione, J.; Wolff, B.; D’Antonio, C.; Dobson, A.; Howarth, R.; Schindler, D.; Schlesinger, W.H.; Simberloff, D.; Swackhamer, D. Forecasting agriculturally driven global environmental change. Science 2001, 292, 281–284. [Google Scholar]
  21. Xin, Q.; Gong, P.; Yu, C.; Yu, L.; Broich, M.; Suyker, A.E.; Myneni, R.B. A production efficiency model-based method for satellite estimates of corn and soybean yields in the midwestern US. Remote Sens 2013, 5, 5926–5943. [Google Scholar]
  22. Frolking, S.; Qiu, J.; Boles, S.; Xiao, X.; Liu, J.; Zhuang, Y.; Li, C.; Qin, X. Combining remote sensing and ground census data to develop new maps of the distribution of rice agriculture in China. Global Biogeochem. Cycle 2002, 16. [Google Scholar] [CrossRef]
  23. Siebert, S.; Portmann, F.T.; Doll, P. Global patterns of cropland use intensity. Remote Sens 2010, 2, 1625–1643. [Google Scholar]
  24. Xiao, X.; Boles, S.; Frolking, S.; Li, C.; Babu, J.Y.; Salas, W.; Moore, B., III. Mapping paddy rice agriculture in south and southeast Asia using multi-temporal MODIS images. Remote Sens. Environ 2006, 100, 95–113. [Google Scholar]
  25. Yu, L.; Wang, J.; Clinton, N.; Xin, Q.; Zhong, L.; Chen, Y.; Gong, P. FROM-GC: 30 m global cropland extent derived through multi-source data integration. Int. J. Digit. Earth 2013, 6. doi: 0.1080/17538947.2013.822574. [Google Scholar]
  26. Hertel, T.W.; Burke, M.B.; Lobell, D.B. The poverty implications of climate-induced crop yield changes by 2030. Global Environ. Change 2010, 20, 577–585. [Google Scholar]
  27. Ganguly, S.; Friedl, M.A.; Tan, B.; Zhang, X.; Verma, M. Land surface phenology from MODIS: Characterization of the collection 5 global land cover dynamics product. Remote Sens. Environ 2010, 114, 1805–1816. [Google Scholar]
  28. Liu, J.; Zhu, W.; Cui, X. A shape-matching cropping index (CI) mapping method to determine agricultural cropland intensities in China using MODIS time-series data. Photogramm. Eng. Remote Sens 2012, 78, 829–837. [Google Scholar]
  29. Tan, B.; Morisette, J.T.; Wolfe, R.E.; Gao, F.; Ederer, G.A.; Nightingale, J.; Pedelty, J.A. An enhanced timesat algorithm for estimating vegetation phenology metrics from MODIS data. IEEE J. STARS 2011, 4, 361–371. [Google Scholar]
  30. Galford, G.L.; Mustard, J.F.; Melillo, J.; Gendrin, A.; Cerri, C.C.; Cerri, C.E. Wavelet analysis of MODIS time series to detect expansion and intensification of row-crop agriculture in Brazil. Remote Sens. Environ 2008, 112, 576–587. [Google Scholar]
  31. Potgieter, A.; Apan, A.; Hammer, G.; Dunn, P. Early-season crop area estimates for winter crops in NE Australia using MODIS satellite imagery. ISPRS J. Photogramm. Remote Sens 2010, 65, 380–387. [Google Scholar]
  32. Wardlow, B.D.; Egbert, S.L. Large-area crop mapping using time-series MODIS 250 m NDVI data: An assessment for the US central great plains. Remote Sens. Environ 2008, 112, 1096–1116. [Google Scholar]
  33. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ 2002, 83, 195–213. [Google Scholar]
  34. Sakamoto, T.; van Phung, C.; Kotera, A.; Nguyen, K.D.; Yokozawa, M. Analysis of rapid expansion of inland aquaculture and triple rice-cropping areas in a coastal area of the Vietnamese Mekong Delta using MODIS time-series imagery. Landscape Urban Plan 2009, 92, 34–46. [Google Scholar]
  35. Li, X.; Wang, X. Changes in agricultural land use in China: 1981–2000. Asian Geogr 2003, 22, 27–42. [Google Scholar]
  36. Lobell, D.B.; Asner, G.P. Cropland distributions from temporal unmixing of MODIS data. Remote Sens. Environ 2004, 93, 412–422. [Google Scholar]
  37. Ozdogan, M. The spatial distribution of crop types from MODIS data: Temporal unmixing using independent component analysis. Remote Sens. Environ 2010, 114, 1190–1204. [Google Scholar]
  38. Ozdogan, M.; Woodcock, C.E. Resolution dependent errors in remote sensing of cultivated areas. Remote Sens. Environ 2006, 103, 203–217. [Google Scholar]
  39. Roy, D.; Lewis, P.; Schaaf, C.; Devadiga, S.; Boschetti, L. The global impact of clouds on the production of MODIS bidirectional reflectance model-based composites for terrestrial monitoring. IEEE Geosci. Remote Sens. Lett 2006, 3, 452–456. [Google Scholar]
  40. Tan, B.; Woodcock, C.; Hu, J.; Zhang, P.; Ozdogan, M.; Huang, D.; Yang, W.; Knyazikhin, Y.; Myneni, R. The impact of gridding artifacts on the local spatial properties of MODIS data: Implications for validation, compositing, and band-to-band registration across resolutions. Remote Sens. Environ 2006, 105, 98–114. [Google Scholar]
  41. Ren, J.; Chen, Z.; Zhou, Q.; Tang, H. Regional yield estimation for winter wheat with MODIS-NDVI data in Shandong, China. Int. J. Appl. Earth Obs. Geoinf 2008, 10, 403–413. [Google Scholar]
  42. Peel, M.C.; Finlayson, B.L.; McMahon, T.A. Updated world map of the Koppen-Geiger climate classification. Hydrol. Earth Syst. Sci. Discuss 2007, 4, 439–473. [Google Scholar]
  43. The USGS EROS Data Center. Available online: http://eros.usgs.gov/ (accessed on 30 December 2013).
  44. Vermote, E.F.; El Saleous, N.Z.; Justice, C.O. Atmospheric correction of MODIS data in the visible to middle infrared: First results. Remote Sens. Environ 2002, 83, 97–111. [Google Scholar]
  45. Wan, Z.; Zhang, Y.; Zhang, Q.; Li, Z.L. Validation of the land-surface temperature products retrieved from Terra Moderate Resolution Imaging Spectroradiometer data. Remote Sens. Environ 2002, 83, 163–180. [Google Scholar]
  46. Loveland, T.; Belward, A. The IGBP-DIS global 1km land cover data set, DISCover: First results. Int. J. Remote Sens 1997, 18, 3289–3295. [Google Scholar]
  47. The USGS Website. Available online: https://lpdaac.usgs.gov/ (accessed on 30 December 2013).
  48. Roy, D.P.; Borak, J.S.; Devadiga, S.; Wolfe, R.E.; Zheng, M.; Descloitres, J. The MODIS land product quality assessment approach. Remote Sens. Environ 2002, 83, 62–76. [Google Scholar]
  49. Grogan, K.; Fensholt, R. Exploring patterns and effects of aerosol quantity flag anomalies in MODIS surface reflectance products in the tropics. Remote Sens 2013, 5, 3495–3515. [Google Scholar]
  50. Chen, J.; Jonsson, P.; Tamura, M.; Gu, Z.; Matsushita, B.; Eklundh, L. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sens. Environ 2004, 91, 332–344. [Google Scholar]
  51. Sakamoto, T.; Yokozawa, M.; Toritani, H.; Shibayama, M.; Ishitsuka, N.; Ohno, H. A crop phenology detection method using time-series MODIS data. Remote Sens. Environ 2005, 96, 366–374. [Google Scholar]
  52. Jonsson, P.; Eklundh, L. Timesat—A program for analyzing time-series of satellite sensor data. Comput. Geosci 2004, 30, 833–845. [Google Scholar]
  53. Jonsson, P.; Eklundh, L. Seasonality extraction by function fitting to time-series of satellite sensor data. IEEE Trans. Geosci. Remote Sens 2002, 40, 1824–1832. [Google Scholar]
  54. Becker-Reshef, I.; Vermote, E.; Lindeman, M.; Justice, C. A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data. Remote Sens. Environ 2010, 114, 1312–1323. [Google Scholar]
  55. Jolly, W.M.; Nemani, R.; Running, S.W. A generalized, bioclimatic index to predict foliar phenology in response to climate. Glob. Chang. Biol 2005, 11, 619–632. [Google Scholar]
  56. Larcher, W.; Bauer, H. Ecological Significance of Resistance to Low Temperature. In Physiological Plant Ecology I; Lange, O.L., Nobel, P.S., Osmond, C.B., Ziegler, H., Eds.; Springer: Heidelberg, Germany, 1981; pp. 403–437. [Google Scholar]
  57. Li, A.; Liang, S.; Wang, A.; Qin, J. Estimating crop yield from multi-temporal satellite data using multivariate regression and neural network techniques. Photogramm. Eng. Remote Sens 2007, 73, 1149–1157. [Google Scholar]
  58. Sari, D.K.; Ismullah, I.H.; Sulasdi, W.N.; Harto, A.B. Detecting rice phenology in paddy fields with complex cropping pattern using time series MODIS data. ITB J. Sci 2010, 42, 91–106. [Google Scholar]
  59. Lobell, D.B.; Schlenker, W.; Costa-Roberts, J. Climate trends and global crop production since 1980. Science 2011, 333, 616–620. [Google Scholar]
  60. National Bureau of Statistics of China. Available online: http://www.stats.gov.cn/tjsj/ndsj/ (accessed on 30 December 2013).
  61. Seto, K.C.; Kaufmann, R.K.; Woodcock, C.E. Landsat reveals China’s farmland reserves, but they’re vanishing fast. Nature 2000, 406. [Google Scholar] [CrossRef]
  62. Smil, V. China’s agricultural land. China Q 1999, 158, 414–429. [Google Scholar]
  63. Congalton, R.G.; Green, K. A practical look at the sources of confusion in error matrix generation. Photogramm. Eng. Remote Sens 1993, 59, 641–644. [Google Scholar]
  64. Gong, P.; Wang, J.; Yu, L.; Zhao, Y.; Zhao, Y.; Liang, L.; Niu, Z.; Huang, X.; Fu, H.; Liu, S. Finer resolution observation and monitoring of global land cover: First mapping results with Landsat TM and ETM+ data. Int. J. Remote Sens 2013, 34, 2607–2654. [Google Scholar]
  65. World Agriculture: Towards 2010. Available online: https://www.mpl.ird.fr/crea/taller-colombia/FAO/AGLL/pdfdocs/wat2010.pdf (accessed on 30 December 2013).
  66. Ellis, E.C. Long-Term Ecological Changes in the Densely Populated Rural Landscapes of China. In Ecosystems and Land Use Change; Defries, S.R., Asner, G.P., Houghton, R.A., Eds.; American Geophysical Union: Washington, DC, USA, 2013. [Google Scholar]
  67. Podwysocki, M.H. Analysis of Field Size Distributions, Lacie Test Sites 5029, 5033, and 5039, Anhwei Province, People’s Republic of China; Report No. N76-27652/6; U.S. Department of Commerce: Virgunua, VI, USA, 1976. [Google Scholar]
  68. Hilker, T.; Wulder, M.A.; Coops, N.C.; Linke, J.; McDermid, G.; Masek, J.G.; Gao, F.; White, J.C. A new data fusion model for high spatial-and temporal-resolution mapping of forest disturbance based on Landsat and MODIS. Remote Sens. Environ 2009, 113, 1613–1627. [Google Scholar]
  69. Xin, Q.; Olofsson, P.; Zhu, Z.; Tan, B.; Woodcock, C.E. Toward near real-time monitoring of forest disturbance by fusion of MODIS and Landsat data. Remote Sens. Environ. 2013, 135, 234–247. [Google Scholar]
  70. Ozdogan, M.; Yang, Y.; Allez, G.; Cervantes, C. Remote sensing of irrigated agriculture: Opportunities and challenges. Remote Sens 2010, 2, 2274–2304. [Google Scholar]
  71. Turner, D.P.; Gower, S.T.; Cohen, W.B.; Gregory, M.; Maiersperger, T.K. Effects of spatial variability in light use efficiency on satellite-based NPP monitoring. Remote Sens. Environ 2002, 80, 397–405. [Google Scholar]
  72. Woodcock, C.E.; Strahler, A.H. The factor of scale in remote sensing. Remote Sens. Environ 1987, 21, 311–332. [Google Scholar]
  73. Gislason, P.O.; Benediktsson, J.A.; Sveinsson, J.R. Random forests for land cover classification. Pattern Recognit. Lett 2006, 27, 294–300. [Google Scholar]
  74. Gopal, S.; Woodcock, C.E.; Strahler, A.H. Fuzzy neural network classification of global land cover from a 1 AVHRR data set. Remote Sens. Environ 1999, 67, 230–243. [Google Scholar]

Appendix Sample Points for Accuracy Assessment

Remotesensing 06 02473f10 1024
Figure A1. The spatial distribution of the sample points used for accuracy assessment. The sample points (4500 in total) were chosen based on an equalized random sampling method. Each of the three defined classes (single, double, and triple cropping) contains 1500 sample points.
Figure A1. The spatial distribution of the sample points used for accuracy assessment. The sample points (4500 in total) were chosen based on an equalized random sampling method. Each of the three defined classes (single, double, and triple cropping) contains 1500 sample points.
Remotesensing 06 02473f10
Remotesensing 06 02473f1 1024
Figure 1. The climate zones in China as derived from Peel et al. [42]. The areas of South China Sea Islands are not displayed; there is little cropland in these islands.
Figure 1. The climate zones in China as derived from Peel et al. [42]. The areas of South China Sea Islands are not displayed; there is little cropland in these islands.
Remotesensing 06 02473f1
Remotesensing 06 02473f2 1024
Figure 2. Flowchart showing the algorithm used for mapping agricultural intensity.
Figure 2. Flowchart showing the algorithm used for mapping agricultural intensity.
Remotesensing 06 02473f2
Remotesensing 06 02473f3 1024
Figure 3. 15-month time series fits for MODIS EVI data shown for pixels with (A) single cropping, (B) double cropping, (C) triple cropping, and (D) frequent clouds. Savitzky-Golay (SG), asymmetric Gaussian (AG), and double logistic (DL) smoothing functions were applied in TIMESAT, and the adaptive Savitzky-Golay filter was chosen for all subsequent analyses.
Figure 3. 15-month time series fits for MODIS EVI data shown for pixels with (A) single cropping, (B) double cropping, (C) triple cropping, and (D) frequent clouds. Savitzky-Golay (SG), asymmetric Gaussian (AG), and double logistic (DL) smoothing functions were applied in TIMESAT, and the adaptive Savitzky-Golay filter was chosen for all subsequent analyses.
Remotesensing 06 02473f3
Remotesensing 06 02473f4 1024
Figure 4. The agricultural mask for mainland China as derived from the IGBP cropland layer in MODIS Land Cover Type product. Both cropland (class 12) and cropland/natural vegetation mosaic (class 14) are used to define the agricultural areas in this study.
Figure 4. The agricultural mask for mainland China as derived from the IGBP cropland layer in MODIS Land Cover Type product. Both cropland (class 12) and cropland/natural vegetation mosaic (class 14) are used to define the agricultural areas in this study.
Remotesensing 06 02473f4
Remotesensing 06 02473f5 1024
Figure 5. Agricultural intensity maps for China shown for (A) 2006 and (B) 2007.
Figure 5. Agricultural intensity maps for China shown for (A) 2006 and (B) 2007.
Remotesensing 06 02473f5
Remotesensing 06 02473f6 1024
Figure 6. Regional-scale views of agricultural intensity in China in 2006 for major agricultural regions: (A) Northeast China Plain, (B) North China Plain and Middle-lower Yangzte Plain, (C) Sichuan Basin, and (D) Pearl River Delta (see annotations in Figure 1). Heilongjiang province (Figure 6A) has the most arable land in China and Henan (Figure 6B) is China’s largest grain-producing province.
Figure 6. Regional-scale views of agricultural intensity in China in 2006 for major agricultural regions: (A) Northeast China Plain, (B) North China Plain and Middle-lower Yangzte Plain, (C) Sichuan Basin, and (D) Pearl River Delta (see annotations in Figure 1). Heilongjiang province (Figure 6A) has the most arable land in China and Henan (Figure 6B) is China’s largest grain-producing province.
Remotesensing 06 02473f6
Remotesensing 06 02473f7 1024
Figure 7. Comparison of (A) arable areas and (B) gross sown areas between estimates from the MODIS Land Cover Type product (MCD12Q1) and national survey data reported by National Bureau of Statistics of China (NBSC) at the province level in 2006. Values for R2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.
Figure 7. Comparison of (A) arable areas and (B) gross sown areas between estimates from the MODIS Land Cover Type product (MCD12Q1) and national survey data reported by National Bureau of Statistics of China (NBSC) at the province level in 2006. Values for R2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.
Remotesensing 06 02473f7
Remotesensing 06 02473f8 1024
Figure 8. Comparisons of gross sown areas in (A) 2006 and (B) 2007 for provinces in China estimated from the agricultural intensity map with national survey data from China Statistical Yearbook. Values for R2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.
Figure 8. Comparisons of gross sown areas in (A) 2006 and (B) 2007 for provinces in China estimated from the agricultural intensity map with national survey data from China Statistical Yearbook. Values for R2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.
Remotesensing 06 02473f8
Remotesensing 06 02473f9 1024
Figure 9. Comparisons of the gross-sown areas for prefectures in Henan province between estimates from the agricultural intensity map and survey data from Henan Statistical Yearbook. Values for R2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.
Figure 9. Comparisons of the gross-sown areas for prefectures in Henan province between estimates from the agricultural intensity map and survey data from Henan Statistical Yearbook. Values for R2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.
Remotesensing 06 02473f9
Table
Table 1. Summary of Cropland Areas in Chinaa in 2006.
Table 1. Summary of Cropland Areas in Chinaa in 2006.
ProvinceArable Land (kha)Gross Sown Area (kha)Cropping Indexb (100%)
Beijing343.9318.092.5%
Tianjin485.6499.4102.8%
Hebei6883.38785.5127.6%
Shanxi4588.63795.482.7%
Inner Mongolia8201.06215.775.8%
Liaoning4174.83796.790.9%
Jilin5578.44954.188.8%
Heilongjiang11,773.010,083.785.7%
Shanghai315.1403.6128.1%
Jiangsu5061.77641.2151.0%
Zhejiang2125.32837.9133.5%
Anhui5971.79172.5153.6%
Fujian1434.72481.3172.9%
Jiangxi2993.45251.4175.4%
Shandong7689.310,736.1139.6%
Henan8110.313,922.7171.7%
Hubei4949.57279.4147.1%
Hunan3953.07977.6201.8%
Guangdong3272.24815.4147.2%
Guangxi4407.96489.2147.2%
Hainan762.1778.1102.1%
Chongqing2067.63487.7168.7%
Sichuan9169.19480.2103.4%
Guizhou4903.54804.198.0%
Yunnan6421.66053.894.3%
Tibet362.6235.064.8%
Shaanxi5140.54201.881.7%
Gansu5024.73726.074.2%
Qinghai688.0476.769.3%
Ningxia1268.81099.386.6%
Xinjiang3985.73731.293.6%
Notes:
aExcludes Taiwan, Hong Kong, and Macao.
bCropping Index here is calculated as arable land area divided by gross sown area [65].
Table
Table 2. Statistical summary for the comparisons between gross sown areas as derived from our maps of agricultural intensity and those obtained from national survey data. Values for R2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided.
Table 2. Statistical summary for the comparisons between gross sown areas as derived from our maps of agricultural intensity and those obtained from national survey data. Values for R2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided.
YearR2RMSE (1000 kha)ME (1000 kha)
20060.9211.170.47
20070.8901.390.60
20080.8991.380.59
20090.8591.480.38
20100.8971.240.29
20110.8861.310.31
Table
Table 3. The error matrix for the China’s agricultural intensity map in 2006.
Table 3. The error matrix for the China’s agricultural intensity map in 2006.
Land Use ClassesVisual Interpretation of MODIS EVI Time Series
User’s Accuracy
Non-CroppingSingle-CroppingDouble-CroppingTriple-Cropping
non-cropping0000
single-cropping11392100792.8%
double-cropping010113594090.6%
triple-cropping035120134589.7%
Producer’s accuracy91.1%86.1%96.6%
overall accuracy = 91.0%

Share and Cite

MDPI and ACS Style

Li, L.; Friedl, M.A.; Xin, Q.; Gray, J.; Pan, Y.; Frolking, S. Mapping Crop Cycles in China Using MODIS-EVI Time Series. Remote Sens. 2014, 6, 2473-2493. https://doi.org/10.3390/rs6032473

AMA Style

Li L, Friedl MA, Xin Q, Gray J, Pan Y, Frolking S. Mapping Crop Cycles in China Using MODIS-EVI Time Series. Remote Sensing. 2014; 6(3):2473-2493. https://doi.org/10.3390/rs6032473

Chicago/Turabian Style

Li, Le, Mark A. Friedl, Qinchuan Xin, Josh Gray, Yaozhong Pan, and Steve Frolking. 2014. "Mapping Crop Cycles in China Using MODIS-EVI Time Series" Remote Sensing 6, no. 3: 2473-2493. https://doi.org/10.3390/rs6032473

APA Style

Li, L., Friedl, M. A., Xin, Q., Gray, J., Pan, Y., & Frolking, S. (2014). Mapping Crop Cycles in China Using MODIS-EVI Time Series. Remote Sensing, 6(3), 2473-2493. https://doi.org/10.3390/rs6032473

Article Metrics

Back to TopTop