Annual dynamic dataset of global cropping intensity from 2001 to 2019

Data collection

Based on the cropland extent, we first introduced a cropland distribution template, the Self-adapting Statistics Allocation Model of Global Cropland (SASAM-GC)¹⁶, as shown in Fig. S1. The global cropland extent map used herein was a global cropland synergy map with a 500-m spatial resolution representing approximately the year 2010, developed by the Smart Agriculture Innovation Team of the Key Laboratory of Agricultural Remote Sensing (AGRIRS) of the Chinese Academy of Agricultural Sciences in cooperation with the International Food Policy Research Institute (IFPRI) and the International Institute of Applied Systems Analysis (IIASA). The overall accuracy of the SASAM-GC products is 90.8%, which is higher than those of existing global farmland products such as the Climate Change Initiative Land Cover (CCI-LC), GlobeLand30, and Medium Resolution Imaging Spectrometer (MODIS) products.

Vegetation indices are often used to depict crop growth, such as the ratio vegetation index (RVI)¹⁷, normalized difference vegetation index (NDVI), and enhanced vegetation index (EVI)¹⁸. Among them, the EVI is the most sensitive to high-biomass regions and less susceptible to atmospheric and soil interference^19,20. MODIS vegetation index datasets are generated every 8 days or 16 days at spatial resolutions of 250 m, 500 m, and 1000 m. The 250-m spatial resolution is considered the best resolution for detecting crops^21,22. Here, we used EVI time series with a spatial resolution of 250 m reported every 16 days in the MODIS product MOD13Q1 as the primary data to calculate the global cropping intensity; these data can be accessed at https://lpdaac.usgs.gov/products/mod13q1v006/.

In formula 1, ρ_NIR, ρ_Red, and ρ_Blue represent the reflectivity of the near-infrared band, the red band, and the blue band, respectively.

The Food and Agricultural Organization of the United Nations statistical data (FAOSTAT) provides long-time-series cropland-related statistical data at the country level and can be accessed at http://www.fao.org/faostat/en/#data. FAOSTAT cropland data have been widely used in a variety of studies evaluating food security and hindcasting historical land-use changes^23,24. Here, we defined the cropland area as the sum of areas characterizing arable land (land under temporary crops, temporary meadows used for mowing or pasture, market and kitchen gardens, and land that is temporarily fallow; abandoned land resulting from shifting cultivation was excluded) and permanent croplands (land cultivated with long-term crops that do not have to be replanted for several years). Additionally, the harvested area is used; this value refers to the area from which a crop is gathered. If the crop under consideration is harvested more than once during a year due to successive cropping (i.e., the same crop is sown or planted more than once in the same field during the same year), the area is counted as many times as the crop is harvested. Therefore, the sown area is recorded only for the crop reported. The statistical cropping intensity is defined as the harvested area divided by the cropland area and is used to validate the consistency between the cropping intensity obtained herein and that reported in the FAOSTAT product at the country level.

Brief algorithm

In this study, a sixth-order polynomial function was used to reconstruct EVI time series for brief and rapid calculations at the global scale²⁵. As the global cropping intensity ranges from 0 to 3 and the sixth-order polynomial function can have 3 maxima at most, we chose the sixth-order polynomial function (formula 2) in this study:

where t is the time-series length referring to the beginning of the time series, EVI(t) is the fitted EVI time series, and a₀, a₁, … a_n are the fitted parameters of the sixth-order polynomial function. Then, the first derivative EVI′(t) and the second derivative EVI″(t) were calculated to find the real numerical solution of the equation’s maxima when EVI′(t) = 0and EVI″(t) < 0.

The procedures followed herein to obtain the global cropping intensity calculation are outlined in Fig. 1:

We first extracted global MODIS EVI time series from 2000 to 2020 according to the SASAM-GC product with areas with cropland probabilities greater than 10% considered as the cropland extent. Moreover, within the global cropland extent, we used the quality control band (QC) to filter a better MODIS EVI time series and interpolate incorrect time values using time series linear interpolation (see Fig. S2 for an example). Using the EVI period from 2000 to 2020, we expanded the EVI time series of each year from 2001 to 2019 by the three months at the end of the previous year and the three months at the beginning of the next year to obtain a total of 18 months in each annual EVI time series. The first minimum value of each 18-month EVI time series was judged as the starting point of the phenological curve. Inconsistent locations were changed to denote the start of planting as the starting point, and a new, convex 12-month EVI time series was reconstructed for each year. Then, this convex EVI time series could be used with sixth-order polynomial fitting to meet the single-peak to three-peaks conditions, as shown in Fig. 2 for single-peak, double-peak, and triple-peak examples, representing multiple cropping intensities once, twice, and three times, respectively. According to the existing literature, filtering using EVI values with peaks value greater than 0.35, differences between the EVI peaks and valleys greater than 0.01, and durations between the peaks and valleys longer than three months could remove some false peaks caused by other vegetation^14,26. Finally, we deleted the number of fake peaks and overfitted peaks according to the experiential threshold and identified the multiple cropping intensity.

According to the above steps, we calculated global time series of farmland pixels annually from 2001 to 2019 to obtain global multiple cropping index distribution maps (see Fig. S3). Examples of increasing and decreasing cropping intensities are also listed in Fig. S4.

Quality control

To make it easy for users to apply our cropping intensity dataset, we introduced a QC band during the production period. The QC band contains four numbers, with values ranging from 0–3 corresponding to “best”, “good”, “fair” and “poor”. If a pixel is judged as “best”, it should satisfy all of the following conditions; “good” pixels meet two of these conditions, “fair” pixels meet one condition, and “poor” pixels meet none of them.

1. a)

   Half of the annual EVI time series should be acceptable according to the MODIS EVI_QC band;

2. b)

   The interpolation counts should be less than half of the whole EVI time series, for a total of 23 counts;

3. c)

   We introduced the delta value here, an estimate of the standard error in predicting value t by EVI(t). The EVI(t) value ± delta should meet the cropping intensity identification conditions.

Source: Ecology - nature.com