Mapping Crop Cycles in China Using MODIS-EVI Time Series

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 OPEN ACCESS Remote Sens. 2014, 6 2474 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.


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][2][3][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][6][7][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][12][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][15][16][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][19][20][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][23][24][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][28][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][31][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.

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].[42].The areas of South China Sea Islands are not displayed; there is little cropland in these islands.

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].

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.

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: where nir  , red  , and blue  are surface reflectance in the near-infrared, red, and blue bands, respectively, and L = 1, 1 C = 6, 2 C = 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.

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.

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.

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).

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].

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.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 (R 2 = 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 R 2 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.1).Heilongjiang province (Figure 6A) has the most arable land in China and Henan (Figure 6B) is China's largest grain-producing province.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.

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][69][70].

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.

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.

Figure 1 .
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 2 .
Figure 2. Flowchart showing the algorithm used for mapping agricultural intensity.

Figure 3 .
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 4 .
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 5 .
Figure 5. Agricultural intensity maps for China shown for (A) 2006 and (B) 2007.

Figure 6 .
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 Figure1).Heilongjiang province (Figure6A) has the most arable land in China and Henan (Figure6B) is China's largest grain-producing province.

Figure 7 .
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 R 2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.

Figure 8 .
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 R 2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.

Figure 9 .
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 R 2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided in figures.

Table 1 .
Summary of Cropland Areas in China a in 2006.
[65]cludes Taiwan, Hong Kong, and Macao.bCroppingIndex here is calculated as arable land area divided by gross sown area[65].

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 R 2 (coefficient of determination), RMSE (root mean square error), and ME (mean error) are provided.

Table 3 .
The error matrix for the China's agricultural intensity map in 2006.