Development of a New Phenology Algorithm for Fine Mapping of Cropping Intensity in Complex Planting Areas Using Sentinel-2 and Google Earth Engine

: Cropping intensity is a key indicator for evaluating grain production and intensive use of cropland. Timely and accurately monitoring of cropping intensity is of great signiﬁcance for ensuring national food security and improving the level of national land management. In this study, we used all Sentinel-2 images on the Google Earth Engine cloud platform, and constructed an improved peak point detection method to extract the cropping intensity of a heterogeneous planting area combined with crop phenology. The crop growth cycle proﬁles were extracted from the multi-temporal normalized difference vegetation index (NDVI) and land surface water index (LSWI) datasets. Results show that by 2020, the area of single cropping, double cropping, and triple cropping in the Henan Province are 52,236.9 km 2 , 74,334.1 km 2 , and 1927.1 km 2 , respectively; the corresponding producer accuracies are 86.12%, 93.72%, and 91.41%, respectively; the corresponding user accuracies are 88.99%, 92.29%, and 71.26%, respectively. The overall accuracy is 90.95%, and the Kappa coefﬁcient is 0.81. Using the sown area in the statistical yearbook data of cities in the Henan Province to verify the extraction results of this paper, the R 2 is 0.9717, and the root mean square error is 1715.9 km 2 . This study shows that using all the Sentinel-2 data, the phenology algorithm, and cloud computing technology has great potential in producing a high spatio-temporal resolution dataset for crop remote sensing monitoring and agricultural policymaking in complex planting areas.


Introduction
The "17 Sustainable Development Goals" (SDGs) issued by the United Nations in 2015 clearly took eliminating hunger and achieving food security as one of its goals and tasks. Therefore, with the rapid growth of the population and economy, how to ensure food security remains the focus of attention of governments and the international community [1][2][3]. The main ways to ensure food security include expanding the area of cropland, increasing the yield per unit area, and increasing the degree of intensification of cropland [4]. Under the background of the decline of cropland area in the Peoples' Republic of China (PRC) [5,6] and the "ceiling effect" of grain yields [7], increasing the degree of cropland intensification has become a key way to increase food production and ensure food security [8]. Cropping intensity is a key indicator to measure the degree of cropland intensification, timely and accurately monitor regional food production, and help achieve SDGs [1,2,9,10].
Remote sensing can fully reflect the spectral differences and phenological characteristics of different crops. At the same time, with the availability of satellite images with different resolutions greatly improved, considerable progress has been made in the production of cropping intensity distribution maps based on remote sensing [11,12]. Over the past few decades, cropping intensity maps have been drawn in different regions using satellite data, such as Brazil [13,14], China [5,15], India [16,17], Arizona [18]. From the selected data sources, the most commonly used satellite data for monitoring cropping intensity are MODIS data and Landsat data. MODIS data have high temporal resolution (1-d), but the spatial resolution is low (250 m to 1000 m). A large number of mixed pixels are difficult to describe the changes in ground crops' details accurately [19][20][21], so it is more suitable for long-term large-scale regional monitoring. Landsat data have high spatial resolution (30-m) and low temporal resolution (8-d to 16-d). Therefore, it is difficult to obtain high-temporal-resolution data of continuous time series [22][23][24][25]. Sentinel-2A/B data not only have high temporal resolution (5-d to 10-d), but also high spatial resolution (10-m to 60-m), which provides more detailed phenological information for monitoring cropping intensity [26]. This provides the possibility to describe many small-scale farmland areas accurately [27]. Thus, it is necessary to evaluate the potential of Sentinel-2A/B data in cropping intensity monitoring and provide new ideas for the fine mapping of cropping intensity in complex planting areas.
From the selected classification approaches, it can be roughly divided into: (1) imagebased spatial statistics and (2) pixel-based temporal statistics. Regarding (1), previous studies extracted cropping intensity based on the differences in the spectrum [28] or vegetation index (VIs) [29,30] in a specific time window, which can best reflect the growth characteristics of the crops. This classification approach needs to collect a wide range of training samples, which not only consumes a lot of time and manpower, but also the quantity and quality of sample points have a great impact on the classification accuracy. Furthermore, due to the influence of atmosphere, clouds, cloud shadow, etc., it is difficult to obtain crop growth characteristics accurately. Regarding (2), this approach generates a cropping intensity map using rule-based algorithms based on constructed time series datasets of VIs [31][32][33][34]. These VIs include the normalized difference vegetation index (NDVI) [14,35], enhanced vegetation index (EVI) [36,37], land surface water index (LSWI) [16], and modified normalized difference water index (mNDWI) [38,39]. This approach is based on the crop phenological characteristics recorded in the VIs' time series. By analyzing the life cycle of the crop, the temporal metrics of crops are obtained, and the classification rules are generated [39][40][41].
Vegetation phenology is a natural phenomenon that reflects the annual cycle of plants under the influence of biological and non-biological factors [25]. Previous studies have proved that crop phenology algorithms have great potential in agricultural remote sensing mapping, such as identifying the start of season (SOS) and the end of season (EOS) [15], identifying the peak growth (PG) and peak drought (PD) periods [42], calculating the length of the growth cycle [43], calculating the growth rate and decline rate of VIs' values [44], and assessing cropping intensity by identifying the number of peaks and troughs in the VIs' time series within a year [5,16,31,45]. A recent study used the NDVI threshold to distinguish between crops and natural vegetation, and used the LSWI threshold to identify bare soil to map the cropping intensity of seven representative grain producing areas in China successfully [46]. In intercropping areas, the next crop is usually sown when the previous crop is not harvested, which results in the bare soil being not being recognized between the harvest and sowing of the two crops. Another study only used the NDVI time series dataset to construct a binary crop phenology profile indicating the growth and non-growth periods, and mapped global cropping intensity at 30 m resolution [47]. The 50% NDVI amplitude threshold was used as a key node in identifying crop growth cycles in this algorithm. However, multiple crop rotations and intercropping systems in various regions will affect the accuracy and reliability of the results. Thus, there is a need to develop a new phenology-based algorithm to identify cropping intensity in complex planting areas.
At present, mapping the cropping intensity over any sizeable region is still a challenge, especially for complex planting areas (e.g., the Henan Province in PRC (details on this region are given in the methods section below)). Firstly, the frequency and dates of highquality images vary greatly in time and space due to the influence of the atmosphere, clouds, cloud shadows, etc. Under normal circumstances, irregular image time series are usually difficult to support the development of classification models [48,49]. Secondly, crop management measures in different regions lead to high intra-class variabilities of crop spectrum and phenology in the Henan Province. For example, in Qixian and Tongxu, Henan Province, corn is sown in some areas in mid-April, and a considerable part of corn is sown in early June. In many of these areas, crop rotation management is implemented, and the dates of crop planting and harvesting vary greatly in the same year and in different years. Therefore, only using a specific time window as the study interval cannot accurately capture the complete dynamic changes of crop growth [30]. Thirdly, winter crops usually have two peaks in one growing season. Therefore, only extracting the peak numbers of VIs' time series to quantify cropping intensity can easily lead to overestimation of the results.
Considering the limitations of using remote sensing to identify cropping intensity in complex planting areas, we developed a phenological-based algorithm on the Google Earth Engine (GEE) cloud platform, and applied the algorithm to generate a cropping intensity product with 10 m spatial resolution in the Henan Province in 2020. The algorithm accurately captures the subtle changes of planting systems in cropland, which can be used to produce cropping intensity products with high-resolution, and provides a reference for the monitoring of highly heterogeneous farmland vegetation.

Study Area
The Henan Province (31 • 23 ~36 • 22 N, 110 • 21 ~116 • 39 E) is located in the middleeastern part of PRC, in the middle and lower reaches of the Yellow River, with an area of 6444 km 2 ( Figure 1). The overall terrain is high in the west and low in the east; surrounded by mountains in the south is the North China Plain in the middle east and Nanyang Basin in the southwest. The north of the Henan Province is located in the warm temperate zone, and a small part of the south is located in the subtropical zone, belonging to the continental monsoon climate from the north subtropical to the warm temperate zone. The annual average temperature is 12-16 • C, and the annual average precipitation is 500-900 mm. It crosses the four major water systems of the Haihe River, the Yellow River, the Huaihe River, and the Yangtze River. Abundant light, heat, and water have laid a good foundation for the development of agriculture in the Henan Province. In 2019, the production of grain in the Henan province was 66.953 million tons, accounting for 10.09% of the national grain output (https://data.cnki.net/yearbook/Single/N2021010028 accessed on 1 May 2021). As an important agricultural product production base of PRC, the Henan Province's main crops are wheat, corn, rice and rape, etc. The cropping system is mainly single cropping and double cropping. Its land use and grain yield are not only related to Henan's own development, but also of great significance to the national food supply.

Sentinel-2A/B Imagery
Sentinel-2 is an Earth observation mission of the European Union Copernicus Program. It consists of two satellites, which were launched in June 2015 (Sentinel-2A) and March 2017 (Sentinel-2B). The revisit period of the dual-satellite synchronization is 5 days, and the Multispectral Instrument (MSI) with a spatial resolution of 10 m, 20 m, and 60 m has 13 spectral bands from visible light to short-wave infrared [50]. The GEE cloud platform contains two types of data: Top-of-Atmosphere Reflectance (TOA) and Surface Reflectance (SR). Since TOA data are sensitive to changes in atmospheric composition and have certain limitations, we used all available SR data in the GEE platform from 2019 to 2021, which correspond to "COPERNICUS/S2_SR" in GEE, because the cropping intensity in 2020 included winter crops from 2019 to 2020 and 2020 to 2021.
The presence of clouds, cloud shadows, and snow in the image will affect the spectral bands, and then affect the calculation results of relevant indexes [51][52][53][54]. The quality assessment band (QA60) in the Sentinel-2 metadata was used to mask the image quality and remove clouds and cloud shadows [38].
We collected a total of 8883 images in the Henan Province during the study period (1 July 2019-1 July 2021). After preprocessing the image according to the above steps, we counted the total observations and good-quality observations of each pixel in the Henan Province, respectively (Figure 2a,b). Most pixels had >120 total observations and >40 goodquality observations in the Henan Province. Moreover, we counted available images for the study (Figure 2c) and observations per pixel in the Henan Province every month (Figure 2d). During the study period, there were >300 images per month, and the total observation and good-quality observations of each pixel were more than 7 times and 2 times per month, respectively.  . The revisit period of the dual-satellite synchronization is 5 days, and the Multispectral Instrument (MSI) with a spatial resolution of 10 m, 20 m, and 60 m has 13 spectral bands from visible light to short-wave infrared [50]. The GEE cloud platform contains two types of data: Top-of-Atmosphere Reflectance (TOA) and Surface Reflectance (SR). Since TOA data are sensitive to changes in atmospheric composition and have certain limitations, we used all available SR data in the GEE platform from 2019 to 2021, which correspond to "COPERNICUS/S2_SR" in GEE, because the cropping intensity in 2020 included winter crops from 2019 to 2020 and 2020 to 2021.
The presence of clouds, cloud shadows, and snow in the image will affect the spectral bands, and then affect the calculation results of relevant indexes [51][52][53][54]. The quality assessment band (QA60) in the Sentinel-2 metadata was used to mask the image quality and remove clouds and cloud shadows [38].
We collected a total of 8883 images in the Henan Province during the study period (1 July 2019-1 July 2021). After preprocessing the image according to the above steps, we counted the total observations and good-quality observations of each pixel in the Henan Province, respectively (Figure 2a,b). Most pixels had >120 total observations and >40 good-quality observations in the Henan Province. Moreover, we counted available images for the study (Figure 2c) and observations per pixel in the Henan Province every month (Figure 2d). During the study period, there were >300 images per month, and the total observation and good-quality observations of each pixel were more than 7 times and 2 times per month, respectively.

Ground Reference Data
We obtained the reference datasets from three aspects to verify the classification results of this paper. Firstly, based on the Google Earth images with 1 m spatial resolution, the stratified random sampling method was used to obtain single, double, and triple cropping validation samples sets. Secondly, during the period from March to June 2020, we carried out seven field works in Zhengzhou City, Xinxiang City, Kaifeng City, Anyang City, Hebi City, Shangqiu City, Zhoukou City, Xuchang City, Luohe City, and Xinyang City in the Henan Province and collected geo-referenced field photos of different cropping intensities. These field photos include main planting systems such as double cropping of

Ground Reference Data
We obtained the reference datasets from three aspects to verify the classification results of this paper. Firstly, based on the Google Earth images with 1 m spatial resolution, the stratified random sampling method was used to obtain single, double, and triple cropping validation samples sets. Secondly, during the period from March to June 2020, we carried out seven field works in Zhengzhou City, Xinxiang City, Kaifeng City, Anyang City, Hebi City, Shangqiu City, Zhoukou City, Xuchang City, Luohe City, and Xinyang City in the Henan Province and collected geo-referenced field photos of different cropping intensities.
These field photos include main planting systems such as double cropping of wheat-corn rotation, single cropping of rice, and triple cropping of vegetable. In addition, during field works, we used drones to obtain multispectral images, which were used as auxiliary data in visual interpretation. The spatial resolution of these multi-spectral images can reach 0.1 m, which clearly shows the different planting systems of the land. Finally, we obtained a total of 1219 samples for single cropping, 2131 samples for double cropping, and 22 samples for triple cropping (Figure 3a).
Furthermore, polygonal buffers were generated from the obtained sample points according to the actual shape and size of the farmland (Figure 3b-d). Through random stratified sampling, we collected a total of 1219 polygons (10,971 pixels) for the single cropping system, 2131 polygons (19,179 pixels) for the double cropping system, and 22 polygons (198 pixels) for the triple cropping system. These polygonal buffers were superimposed on Google Earth images, combined with geo-referenced field photos and drone multi-spectral images, and finally obtained single, double, and triple cropping validation grids through visual interpretation. wheat-corn rotation, single cropping of rice, and triple cropping of vegetable. In addition, during field works, we used drones to obtain multispectral images, which were used as auxiliary data in visual interpretation. The spatial resolution of these multi-spectral images can reach 0.1 m, which clearly shows the different planting systems of the land. Finally, we obtained a total of 1219 samples for single cropping, 2131 samples for double cropping, and 22 samples for triple cropping (Figure 3a). Furthermore, polygonal buffers were generated from the obtained sample points according to the actual shape and size of the farmland (Figure 3b-d). Through random stratified sampling, we collected a total of 1219 polygons (10,971 pixels) for the single cropping system, 2131 polygons (19,179 pixels) for the double cropping system, and 22 polygons (198 pixels) for the triple cropping system. These polygonal buffers were superimposed on Google Earth images, combined with geo-referenced field photos and drone multispectral images, and finally obtained single, double, and triple cropping validation grids through visual interpretation.

Other Data
In order to obtain the cropping intensity map accurately, it is necessary to remove the non-cropland pixels in the study area. We used the cropland pixels in the 10 m resolution global land cover product (FROM-GLC10) to mask the constructed good-quality time series datasets, which were published by the Department of Earth System Sciences, Tsinghua University [55]. The map can be freely downloadable from (http://data.ess.tsinghua.edu.cn accessed on 1 May 2021). At the same time, we used statistical data from 17 cities in the Henan Province to compare and analyze the classification results of this paper. The statistical data come from the website of the National Bureau of Statistics (http://www.stats.gov.cn accessed on 1 May 2021). Figure 4 shows the process of extracting cropping intensity in this study. Firstly, we obtained all available Sentinel-2A/B SR data from July 2019 to July 2021 from the GEE

Other Data
In order to obtain the cropping intensity map accurately, it is necessary to remove the non-cropland pixels in the study area. We used the cropland pixels in the 10 m resolution global land cover product (FROM-GLC10) to mask the constructed good-quality time series datasets, which were published by the Department of Earth System Sciences, Tsinghua University [55]. The map can be freely downloadable from (http://data.ess. tsinghua.edu.cn accessed on 1 May 2021). At the same time, we used statistical data from 17 cities in the Henan Province to compare and analyze the classification results of this paper. The statistical data come from the website of the National Bureau of Statistics (http://www.stats.gov.cn accessed on 1 May 2021). Figure 4 shows the process of extracting cropping intensity in this study. Firstly, we obtained all available Sentinel-2A/B SR data from July 2019 to July 2021 from the GEE platform, and used the quality assessment band (QA60) to remove the bad-quality observations in the Sentinel-2A/B SR data. Secondly, we used a land cover product (FROM-GLC10) with a spatial resolution of 10 m to mask the land in the Henan Province and obtained data on cropland. Thirdly, we calculated the NDVI and LSWI indices of the Henan Province pixel by pixel, and reconstructed the NDVI and LSWI time series, including the 10-day composite, linear interpolation, and Savitzky-Golay filter smoothing (only for NDVI time series). Fourth, combining the LSWI time series and the growth period length (GPL), we eliminated the "false peaks" in the NDVI time series. Finally, according to the number of effective peaks obtained, we determined the cropping intensity of the Henan Province pixel by pixel, and obtained a cropping intensity map of the Henan Province with a spatial resolution of 10 m. platform, and used the quality assessment band (QA60) to remove the bad-quality observations in the Sentinel-2A/B SR data. Secondly, we used a land cover product (FROM-GLC10) with a spatial resolution of 10 m to mask the land in the Henan Province and obtained data on cropland. Thirdly, we calculated the NDVI and LSWI indices of the Henan Province pixel by pixel, and reconstructed the NDVI and LSWI time series, including the 10-day composite, linear interpolation, and Savitzky-Golay filter smoothing (only for NDVI time series). Fourth, combining the LSWI time series and the growth period length (GPL), we eliminated the "false peaks" in the NDVI time series. Finally, according to the number of effective peaks obtained, we determined the cropping intensity of the Henan Province pixel by pixel, and obtained a cropping intensity map of the Henan Province with a spatial resolution of 10 m.

Index Calculation and Time Series Reconstruction
The NDVI value of crops gradually increased from the green-up stage, reaching the maximum value at the peak growth stage, and then decreased from the mature stage to harvest ( Figure 5). Therefore, the NDVI index can be used to indicate surface vegetation coverage and growth status. The change trend in the LSWI value during the crop growth period was roughly the same as that of the NDVI value ( Figure 5). Soil moisture is much

Index Calculation and Time Series Reconstruction
The NDVI value of crops gradually increased from the green-up stage, reaching the maximum value at the peak growth stage, and then decreased from the mature stage to harvest ( Figure 5). Therefore, the NDVI index can be used to indicate surface vegetation coverage and growth status. The change trend in the LSWI value during the crop growth period was roughly the same as that of the NDVI value ( Figure 5). Soil moisture is much lower than crops, so LSWI values are usually low before crops are planted and after harvest. This can be used to judge the condition of the cropland before sowing and after harvest [56,57]. The calculation formulas of NDVI [58] and LSWI [59] are as follows: where the ρ NIR , ρ RED , and ρ SWIR represent the near infrared band, the red band, and the shortwave infrared band, respectively. lower than crops, so LSWI values are usually low before crops are planted and after harvest. This can be used to judge the condition of the cropland before sowing and after harvest [56,57]. The calculation formulas of NDVI [58] and LSWI [59] are as follows: where the ρ NIR , ρ RED , and ρ SWIR represent the near infrared band, the red band, and the shortwave infrared band, respectively. Since the angle of sunlight, the angle of satellite observation, clouds, etc. change greatly over time [60], there is a lot of noise in the original NDVI and LSWI time series data. Thus, the data need to be further processed. We used the following three steps to reconstruct the time series data: (1) Image compositing can reduce the influence of cloud and temporal uneven observation [61,62]. In order to obtain time series of equal lengths and intervals, we calculated the maximum value of the NDVI and the average value of the LSWI every 10 days to generate a composite image for 10 days. The maximum value of the NDVI is selected because the maximum value of the NDVI can reflect the true growth conditions of vegetation during this period, and the average value of the LSWI is selected because the LSWI is very sensitive to changes in soil and vegetation moisture, and the average value of the LSWI is closer to the growth changes in vegetation during that period. After 10 days of compositing processing, we obtained a total of 66 data points (the black dots in Figure 5). (2) Affected by factors such as severe weather, cloud or snow, etc., good-quality observations may not be available in some areas. When there were no good-quality observations in 10 days, based on good-quality observations around 10 days, we used the linear interpolate [63] to fill the gaps in the time series of the NDVI and LSWI. Since the angle of sunlight, the angle of satellite observation, clouds, etc. change greatly over time [60], there is a lot of noise in the original NDVI and LSWI time series data. Thus, the data need to be further processed. We used the following three steps to reconstruct the time series data: (1) Image compositing can reduce the influence of cloud and temporal uneven observation [61,62]. In order to obtain time series of equal lengths and intervals, we calculated the maximum value of the NDVI and the average value of the LSWI every 10 days to generate a composite image for 10 days. The maximum value of the NDVI is selected because the maximum value of the NDVI can reflect the true growth conditions of vegetation during this period, and the average value of the LSWI is selected because the LSWI is very sensitive to changes in soil and vegetation moisture, and the average value of the LSWI is closer to the growth changes in vegetation during that period. After 10 days of compositing processing, we obtained a total of 66 data points (the black dots in Figure 5). (2) Affected by factors such as severe weather, cloud or snow, etc., goodquality observations may not be available in some areas. When there were no good-quality observations in 10 days, based on good-quality observations around 10 days, we used the linear interpolate [63] to fill the gaps in the time series of the NDVI and LSWI. After linear interpolation processing, the available data increased to 77 data points. (3) Even if the NDVI time series data are composited for 10 days, the false high value, the continuous low NDVI value during the synthesis period, and other noise residuals still exist [60], which may cause jagged fluctuations in the NDVI time series. Savitzky-Golay filtering [64] was used to smooth the NDVI dataset (the solid line in Figure 5). Since the growth period of most crops exceeds 90 days, we used a moving window of size 9 and a 2 order of filter [15]. The LSWI is more sensitive to surface moisture, so we did not smooth the LSWI time series.
Based on geo-referenced field photos and Google Earth images, we selected four typical sites representing single, double, and triple cropping intensity. We established time series of NDVI and LSWI values at these sites to obtain phenological indicators of different cropping intensities ( Figure 5). Single-cropping crops have a growth cycle within a year, and their NDVI time series usually has only one peak in the year (Figure 5a). Doublecropping crops reach their peak growth stage twice in a year, and their NDVI time series usually has two (Figure 5c) or three peaks (Figure 5b) during the year. For double-cropping rice (Figure 5c) or other non-winter crops double-cropping rotation systems, the NDVI time series usually shows two normal peaks. However, for winter wheat-corn or winter crops and other crops in the double-cropping rotation mode, the NDVI time series has a lower false peak before pre-winter (1 October 2020-1 February 2021) (Figure 5b). Figure 5 shows that from sowing to harvest, the NDVI curve of the crop first rises to the peak, and then gradually drops to the bottom. This process occurs once in a year for the single cropping area, twice in a year for the double cropping area, and three times in a year for the triple cropping area. Therefore, a growth cycle usually corresponds to a peak (except for winter crops). Then, we use peak detection methods to identify potential peaks (PP) and potential trough (PT) from the crop profile of the NDVI time series in a year. Before and after crops planting, the cropland is bare soil. The LSWI of bare soil is much lower than that of vegetation. Therefore, the LSWI value during the bare soil period can be used to identify the complete growth cycle of the crop [57,65] and eliminate the extra peak value of winter crops.

Algorithm for Crop Intensity Mapping
However, it was obviously unreasonable to determine the cropping intensity of the pixel only based on the number of peaks obtained in the above steps. It is not suitable for the intercropping crop [66]. Crops in intercropping systems usually exhibit two or more crop rotations, and the latter crop is planted when the former crop has not yet been harvested. Therefore, it is impossible to distinguish the two growth cycles of crops using the LSWI to detect bare soil. The double cropping in the intercropping system was incorrectly identified as single cropping. Therefore, we tried to resolve the problem that crops in the intercropping systems are not easy to identify, specifically, the adjacent peaks and the trough between the two peaks in the NDVI time series.
Moreover, the independent growth period of major crops in PRC is more than 90 days [67]. In other words, effective continuous peek recognition must meet 90 days during the crop growth season. The final effective peak (EP) was obtained through the above restriction conditions, and the number of EPs was used to determine the cropping intensity of the pixel. Specific steps are as follows:

1.
Identified PP and PT through the moving window. In the NDVI time series, when the NDVI value of a certain point was higher (lower) than its two neighboring values, it was judged as the peak (trough) value.

2.
Marked the PT between two consecutive PPs. In order to avoid identifying the autumn growth of winter wheat and other crops as a separate growth cycle, we used the LSWI value of the crop at PT as one of the limiting conditions for extracting PP. After the crops were harvested, bare soil appeared on the surface pixels, and a growth cycle of the crops ended. Thus, bare soil was an important indicator for judging the complete growth cycle. According to previous studies, due to the low water content of soil in northern PRC, when the LSWI value of a pixel was less than 0 at PTs [16,65], it was recognized as bare soil Equation (3). If the bare soil signal is detected between two consecutive peaks, the two consecutive peaks are identified as two growth cycles.
where LSWI PT is the LSWI value at the potential trough of NDVI.

3.
Eliminate invalid PPs. Since it is usually difficult to detect LSW in the bare soil period for intercropping crops, we applied a supplementary algorithm to identify its growth cycle. When the NDVI value was between 0.2 and 0.5, the pixels were usually considered to be a mixture of vegetation and bare soil. When the NDVI value was greater than 0.5, the pixels were usually considered as complete vegetation [36,68]. Therefore, for neighboring PPs, if both PPs are greater than 0.5, and the PT between the two PPs is less than 0.5, then both PPs are recorded. Otherwise, only one is recorded; see Equation (4).
where NDVI PPA and NDVI PPB are the neighboring potential peak values of the NDVI, and NDVI PTC is the potential trough value of the NDVI between the NDVI PPA and NDVI PPB , respectively.

4.
Used the growth period length (GPL) to filter the effective PPs further. For crops in intercropping systems, the crop growth cycle should be greater than 90 days. Therefore, when the time span of the wave of the NDVI time series is greater than 90 days, it is recorded as an effective growth cycle, and its peak value is identified by Equation (5). GPL > 90 (5) where GPL is the growth period length of crops.

5.
Based on the number of effective peaks (EP) obtained above, the cropping intensity of each pixel was determined. If there was only 1 EP in a pixel in 2020, it would be marked as a single cropping. If there are 2 or 3 EPs, it would be marked as double cropping or triple cropping. Accordingly, a 10 m spatial resolution cropping intensity map in the Henan Province in 2020 is generated.

Accuracy Assessment and Intercomparison
The four evaluation indicators of mapping accuracy (PA), user accuracy (UA), overall accuracy (OA), and the Kappa coefficient [69] were used to verify the accuracy and application potential of the classification method in this study. We use the polygon buffers obtained in Section 2.2.2 and the classification results of this study to construct a confusion matrix, and calculate the PA, UA, OA, and Kappa coefficients, respectively. At the same time, in order to verify the accuracy of the extracted results in this paper indirectly, we compared the planting area extracted from remote sensing data with the crop planting area of the national prefecture-level city and county statistics data.

Annual Map of Crop Intensity in 2020
As described in Section 2.3.2, we determined the effective peak number of all cropland pixels in 2020 and drew a cropping intensity map with 10 m resolution in the Henan Province in 2020 ( Figure 6). The planting systems in the Henan Province are dominated by single cropping (52,236.9 km 2 ) and double cropping (74,334.1 km 2 ), and the distribution of triple cropping is very small (1927.1 km 2 ). A single cropland is mostly concentrated in the southern, western, and northwestern edges of the study area, accounting for 40.65% of the cropland in the Henan Province. Double cropland is mainly distributed in the central and eastern part of the study area, accounting for 57.85% of the cropland in the Henan Province. The distribution of triple cropland is scattered and irregular, mainly due to human influence (e.g., people's preferences and needs), accounting for only 1.49% of the cropland in the Henan Province. 40.65% of the cropland in the Henan Province. Double cropland is mainly distributed in the central and eastern part of the study area, accounting for 57.85% of the cropland in the Henan Province. The distribution of triple cropland is scattered and irregular, mainly due to human influence (e.g., people's preferences and needs), accounting for only 1.49% of the cropland in the Henan Province.

Accuracy Assessment
We used the validation samples set described in Section 2.2.2 to evaluate the accuracy of the cropping intensity map with 10 m spatial resolution in the Henan Province in 2020. Using the confusion matrix between the validation samples set and resultant map, we obtained that the OA and Kappa coefficient of the intensity map were 90.95% and 0.81, respectively ( Table 1). The UA of the single cropping system is greater than the PA, which indicates that the single cropping system is more likely to be omitted, resulting in the underestimation of the planting area of the single cropping system. At the same time, the PA of the double cropping system and triple cropping system is greater than that of UA, which indicates that the double cropping system and triple cropping system are more likely to be misclassified, leading to overestimation of the planting area of the double-and triple-cropping system.

Accuracy Assessment
We used the validation samples set described in Section 2.2.2 to evaluate the accuracy of the cropping intensity map with 10 m spatial resolution in the Henan Province in 2020. Using the confusion matrix between the validation samples set and resultant map, we obtained that the OA and Kappa coefficient of the intensity map were 90.95% and 0.81, respectively ( Table 1). The UA of the single cropping system is greater than the PA, which indicates that the single cropping system is more likely to be omitted, resulting in the underestimation of the planting area of the single cropping system. At the same time, the PA of the double cropping system and triple cropping system is greater than that of UA, which indicates that the double cropping system and triple cropping system are more likely to be misclassified, leading to overestimation of the planting area of the double-and triple-cropping system. Furthermore, we selected three areas with mixed planting systems. Each area is a square with a side length of 1 km. They were visually interpreted through ground reference data and generated polygons with different types of crop intensity. These polygons were used to evaluate the classification accuracy of our cropping intensity map quantitatively. Figure 7 intuitively describes the results of the accuracy assessment. Shrubs, roads, and houses are the main areas of misclassification. This is mainly because the side length of our pixel is 10 m, which is usually less than the width of the rural road. Overall, the consistency between our classification results and ground reference data is very obvious. Furthermore, we selected three areas with mixed planting systems. Each area is a square with a side length of 1 km. They were visually interpreted through ground reference data and generated polygons with different types of crop intensity. These polygons were used to evaluate the classification accuracy of our cropping intensity map quantitatively. Figure 7 intuitively describes the results of the accuracy assessment. Shrubs, roads, and houses are the main areas of misclassification. This is mainly because the side length of our pixel is 10 m, which is usually less than the width of the rural road. Overall, the consistency between our classification results and ground reference data is very obvious. We also compared the crop statistics of cities (including 17 cities in total) in the Henan Province with our results (Figure 8). There is a significant linear relationship between the sown area derived from our results and the sown area derived from the statistical data (slope = 0.7692, R 2 = 0.9717, p < 0.001, RMSE = 1715.9 km 2 ). The sown area was smaller in We also compared the crop statistics of cities (including 17 cities in total) in the Henan Province with our results (Figure 8). There is a significant linear relationship between the sown area derived from our results and the sown area derived from the statistical data (slope = 0.7692, R 2 = 0.9717, p < 0.001, RMSE = 1715.9 km 2 ). The sown area was smaller in the cropping intensity maps than in the national statistical data for almost all cities, mainly because the cropland area was smaller in FROM-GLC10 than in the national statistical yearbook. the cropping intensity maps than in the national statistical data for almost all cities, mainly because the cropland area was smaller in FROM-GLC10 than in the national statistical yearbook.

Spatial Distribution of Cropping Intensity in Henan Province
The cropping intensity and crop types of cropland are not only affected by natural processes, but also by anthropogenic drivers [70]. The hydrothermal conditions in the south of the Henan Province are more abundant than those in the north. Therefore, in theory, the cropping intensity in the south should be higher than that in the north. However, our results indicate that most areas in the south (such as Xinyang City) are single cropping systems. The Xinyang area has sufficient rainfall, dense ponds, and high temperatures, which provide a good growth environment for rice growth. Through our field investigation, we found that the main factor affecting Xinyang's planting system is the loss of population. As more and more young people migrate from rural to urban areas, the rural population of Xinyang has dropped sharply. Therefore, most cropland has changed from a double rice cropping system or triple rice cropping system to a single rice cropping system. Different from the single rice cropping system in Xinyang, the cropping system in the mountainous areas around Nanyang is mainly single peanuts. This is because these areas have high terrain, poor soil quality, small cropland, and insufficient irrigation conditions. Therefore, farmers plant more drought-tolerant peanuts. In addition, the special growth characteristics of peanuts prevent it from being planted twice in the same area within a year, which is also the reason for the implementation of the single peanut cropping system in these areas. In Kaifeng City, the main cropping systems are wheat-peanut, garlic-corn, wheat-watermelon rotation. This is mainly because most of the cropland in Kaifeng City is sandy soil, with low soil water content and poor water retention capacity, which is suitable for planting drought-tolerant crops.

Algorithm Improvement
The peak point detection method has been widely used to map cropping intensity [31,32,49]. This method is simple and feasible, but it is more sensitive to noise in the time series. For example, the vegetation index time series of winter crops usually have two peaks in one growth cycle (Figure 5b). Therefore, it is necessary to use a supplementary algorithm to eliminate false peaks [5,71]. Gray et al. (2014) used the calendar-year thermal growth season to eliminate the extra peak of winter crops in autumn [32]. However, this method does not consider additional peaks caused by weather factors in other seasons. For example, summer crops may also have multiple peaks between April and October.

Spatial Distribution of Cropping Intensity in Henan Province
The cropping intensity and crop types of cropland are not only affected by natural processes, but also by anthropogenic drivers [70]. The hydrothermal conditions in the south of the Henan Province are more abundant than those in the north. Therefore, in theory, the cropping intensity in the south should be higher than that in the north. However, our results indicate that most areas in the south (such as Xinyang City) are single cropping systems. The Xinyang area has sufficient rainfall, dense ponds, and high temperatures, which provide a good growth environment for rice growth. Through our field investigation, we found that the main factor affecting Xinyang's planting system is the loss of population. As more and more young people migrate from rural to urban areas, the rural population of Xinyang has dropped sharply. Therefore, most cropland has changed from a double rice cropping system or triple rice cropping system to a single rice cropping system. Different from the single rice cropping system in Xinyang, the cropping system in the mountainous areas around Nanyang is mainly single peanuts. This is because these areas have high terrain, poor soil quality, small cropland, and insufficient irrigation conditions. Therefore, farmers plant more drought-tolerant peanuts. In addition, the special growth characteristics of peanuts prevent it from being planted twice in the same area within a year, which is also the reason for the implementation of the single peanut cropping system in these areas. In Kaifeng City, the main cropping systems are wheat-peanut, garlic-corn, wheat-watermelon rotation. This is mainly because most of the cropland in Kaifeng City is sandy soil, with low soil water content and poor water retention capacity, which is suitable for planting drought-tolerant crops.

Algorithm Improvement
The peak point detection method has been widely used to map cropping intensity [31,32,49]. This method is simple and feasible, but it is more sensitive to noise in the time series. For example, the vegetation index time series of winter crops usually have two peaks in one growth cycle (Figure 5b). Therefore, it is necessary to use a supplementary algorithm to eliminate false peaks [5,71]. Gray et al. (2014) used the calendar-year thermal growth season to eliminate the extra peak of winter crops in autumn [32]. However, this method does not consider additional peaks caused by weather factors in other seasons. For example, summer crops may also have multiple peaks between April and October. Recently, Liu et al. (2020b) also used the peak point detection method to identify the peak of the crop growth curve, to extract the SOS and EOS on the crop growth curve, and finally map the cropping intensity of all seven typical grain production areas in China [46]. However, his algorithm has certain problems in identifying continuous peaks. Through our field investigation, we found that even the same crops in the same area have very different planting times. For example, in Kaifeng City, Henan Province, some garlic was planted in August of last year and harvested at the end of May of the next year, but a considerable part of garlic was planted in mid-September of last year and harvested in June of the next year. In Xinxiang City, Henan Province, some peanuts were sown in May and harvested in September, and some peanuts were sown at the end of March and harvested in August. Therefore, it is possible to identify continuous peaks only by limiting the NDVI trough from June to September. The length of the growth cycle is not limited by the sowing and harvesting time of crops. In this study, we not only used NDVI and LSWI time series to determine the potential growth cycle of the crop, but also further filtered the potential growth cycle according to the length of the crop growth cycle, which solves the problem of different sowing and harvest times in the same planting system in the Henan Province.
Moreover, the spatial resolution of 30 m is difficult to capture some small farmland and mixed agricultural landscapes (Figure 9). Higher spatial resolution of cropping intensity information not only helps decision makers to implement more accurate land management, but also helps to understand the dynamic changes of cropland under human influence [70]. In this study, we used Sentinel-2 images with a spatial resolution of 10 m to map cropping intensity in the Henan Province. Compared with previous studies using AVHRR (1000 m), MODIS (250 m to 1000 m), or Landsat (30 m) data, our results provide more spatial details (10 m) and successfully achieved fine mapping of large areas of complex planting areas. Higher resolution cropping intensity data can also improve the accuracy of many surface and crop models based on farmland planting systems [72][73][74].

Uncertainty
Although the method proposed in this study has successfully mapped fine (10 m) cropping intensity maps in complex planting areas, our analysis also pointed out several issues worthy of attention. Firstly, the accessibility of reliable cropland boundary data has always been a key constraint for cropland-based applications. In this study, we used the 10 m resolution global land cover product FROM-GLC10 [55] as the base cropland/noncropland map. Although our research shows that the FROM-GLC10 product accurately represents the actual situation of cropland in the Henan Province, there are still certain classification errors, especially the rural roads and surrounding farmland in urban and rural residential areas (Figure 7). These errors are liable to be propagated to the final cropping intensity outputs. Meanwhile, the cropland data show the land cover in 2017, which may lead to a certain gap between our results (2020) and the actual cropland. With the release of the latest land use data, the cropland data synchronized with the study period can be used for more accurate remote sensing mapping in the future.
Secondly, in the process of constructing the time series, we used linear interpolation based on neighboring pixels to fill the blanks. If there is a lack of good-quality observation data when the peak (trough) is reached, then the real peak (trough) cannot be obtained by interpolation. In addition, due to cloud contamination, atmospheric variability, and bidirectional effects, noises appear in the NDVI time series. For these noises, we used the Maximum Value Composite (MVC) technique [75] and Savitzky-Golay filtering [76]. Figure 5 shows that the Savitzky-Golay filtering captures curve details and risk of underfitting. However, it is necessary to evaluate various smoothing methods (e.g., least-squares fits to asymmetric Gaussian functions [64] and double logistic functions [77], Whittaker smoother [78], wavelet transforms [79] and the best index slope extraction [80], etc.), which are expected to find more suitable methods to improve the accuracy of the results. Moreover, the microwave images from Sentinel-1 can easily pass through clouds, and can identify bare soil according to the degree of depolarization attained by the illuminated signal at the incident plane. In future research, remote sensing mapping can be combined with Sentinel-1A/B data not affected by weather.
Thirdly, the identification of the peak is also associated with assumptions that should be further verified. In some areas of Xinxiang City, the cropland near both sides of the road is usually a mixture of crops and evergreen trees. Since the LSWI of evergreen tree pixels always exceeds 0, the most mixed pixels of crops and evergreen trees are identified as single cropping. To exclude the extra peak of crops, our algorithm requires the length of each growth cycle to exceed 90 days. This threshold assumption, however, may be violated for special crops. For example, the growth period of lettuce and spinach is usually 20-40 days, so the growth cycle of these crops will be wrongly eliminated.

Uncertainty
Although the method proposed in this study has successfully mapped fine (10 m) cropping intensity maps in complex planting areas, our analysis also pointed out several issues worthy of attention. Firstly, the accessibility of reliable cropland boundary data has always been a key constraint for cropland-based applications. In this study, we used the 10 m resolution global land cover product FROM-GLC10 [55] as the base cropland/noncropland map. Although our research shows that the FROM-GLC10 product accurately

Conclusions
In this study, we developed a new phenological-based cropping intensity extraction method, and proposed a cropping intensity mapping system based on Sentinel-2A/B data. The system was used to map cropping intensity with 10 m resolution in the Henan Province in 2020. The overall accuracy is 90.95%, and the Kappa coefficient is 0.81. Compared with the cropping intensity map obtained from MODIS images, Landsat images, and images' fusion of multi-source data in previous studies, the cropping intensity map in this study has more detailed spatial distribution and higher accuracy, and shows greater potential in estimating cropping intensity in complex planting areas. This dataset can provide quantitative information about the changes in the farming systems and improve our ability to predict grain yield.