2.1.2. Data

MODIS MOD13Q1 Enhanced Vegetation Index (EVI) products in 2018 with 250 m resolution were used and were synthesized over 16 days based on the maximum value composite principle. The image tiles including h27v05, h27v06, h28v05, and h28v06 (h: horizontal, v: vertical) were downloaded from the National Aeronautics and Space

Administration website (https://modis.gsfc.nasa.gov/, accessed on 1 March 2023). Image preprocessing, including mosaic, resampling and reprojection, was conducted using the MODIS Reprojection Tool. To eliminate the interference of clouds, snow, shadows and other factors, Savitzky–Golay filtering was adopted to reconstruct the original time series. Data quality information was extracted based on the pixel reliability layer of MOD13Q1 products.

Sentinel-2 data with a spatial resolution of 10 m were used to generate validation data. Sentinel-2A (launched 2015) and Sentinel-2B (launched 2017) sensors together offer 5-day revisit with global coverage [38]. NDVI was calculated for the Sentinel-2 data by using near infrared and red band.

Sentinel-2 NDVI composite from three key phenological phases was used to prepare the reference cropping intensity data. The first phase was in mid March, which represents the peak of the first growing season of double-cropping (termed GS1). The second phase was from late May to early June, which is the transition period between the two growing seasons (termed TGS). Both single-cropping and double-cropping had low EVI values at this stage. The third phase was in mid July, which is the peak of the second growing season of double-cropping (termed GS2). Double-cropping is shown in magenta, and single-cropping in blue in the false-color composite image (Figure 2a,b).

**Figure 2.** Unsupervised classification of Sentinel-2 images: (**a**,**b**) are false-color composites of Sentinel-2 NDVI from three phenological phases (GS1, TGS, GS2); (**c**,**d**) are ISODATA classification results for the two regions (A) and (B).

The classification was conducted based on these NDVI composites using ISODATA, in which the classification number was between 10 and 20, and the iteration time was set to 20. The classification results were merged into three categories: double-cropping, single-cropping and non-cropland (Figure 2c,d). Finally, they were aggregated to 250 m fractional images to match the spatial resolution of MODIS data.

The crop pattern sample data were used to evaluate the accuracy of the Sentinel-2 derived cropping intensity data. The producers' accuracies (PA) of single-cropping and double-cropping were 91.06% and 91.16%, respectively, their users' accuracies (UA) were 92.44% and 90.58%, respectively (Table 1), and the overall accuracy (OA) and Kappa coefficient were 92.6% and 0.888, respectively, suggesting the reliability of the Sentinel-2 derived cropping intensity data.


**Table 1.** Accuracies of Sentinel-2 derived cropping intensity data using crop pattern sample data.

To evaluate the accuracy of the Sentinel-2 derived cropping intensity data, we used 1779 crop samples (537 single-cropping samples, 622 double-cropping samples and 620 other samples) as reference data. These crop patterns were transferred into cropping intensity: single cropping, double cropping and other (including non-crop cultivation and non-cropland). These samples were from a filed survey in 2018 as ground truth (Figure 1), and augmented by visual interpreting high-resolution images on the Google Earth Engine platform. The ground truth data were collected using a mobile application, GPSTool 4.0. The augmentation were delineated manually by overlaying ground truth data with highresolution images. The sample data were roughly distributed evenly throughout the validation area and were independent and identically distributed.

In addition to the Sentinel-2 derived reference data, three cropping intensity products were further used for the validation. The first was the MCD12Q2 V6 Land Cover Dynamics product, which provides global estimates of the timing of vegetation phenology at 500 m resolution [39]. The NumCycles layer in MCD12Q2 provides the total number of valid vegetation cycles in a year. The annual average of the NumCycles was calculated and then used as reference data. The second product was the Global Cropping Intensity (GCI) dataset, which is an annual global multi-cropping index distribution map covering the period from 2001 to 2019 at 250 m resolution [40]. The third was a global cropping intensity map dataset at 30 m resolution (GCI30) from 2016 to 2018 [41].
