Mapping Ratoon Rice Planting Area in Central China Using Sentinel-2 Time Stacks and the Phenology-Based Algorithm

Mapping rice cropping systems is important for grain yield prediction and food security assessments. Both singleand double-season rice are the dominant rice systems in central China. However, because of increasing labor shortages and high costs, there has been a gradual decline in double-season rice. Ratoon rice (RR) has been proposed as an alternative system that balances the productivity, cost, and labor requirements of rice cultivation. RR has been expanding in central China, encouraged by the improved cultivars, machinery, and favorable policies. However, to our knowledge, the distribution of RR has not been mapped with remote sensing techniques. This study developed a phenology-based algorithm to map RR at a 10 m resolution in Hubei Province, Central China, using dense time stacks of Sentinel-2 images (cloud cover <80%) in 2018. The key in differentiating RR from the other rice cropping systems is through the timing of maturity. We proposed to use two contrast vegetation indices to identify RR fields. The newly-developed yellowness index (YI) calculated with the reflectance of blue, green, and red bands was used to detect the ripening phase, and the enhanced vegetation index (EVI) was used to detect the green-up of the second-season crop to eliminate the misclassification caused by stubbles left in the field. The RR map demonstrated that RR was mainly distributed in the low alluvial plains of central and southern Hubei Province. The total planting area of RR in 2018 was 2225.4 km2, accounting for 10.03% of the total area of paddy rice fields. The overall accuracy of RR, non-RR rice fields, and non-rice land cover types was 0.76. The adjusted overall accuracy for RR and non-RR was 0.91, indicating that the proposed YI and the phenology-based algorithm could accurately identify RR fields from the paddy rice fields.


Introduction
Rice (Oryza sativa) cropping systems are of great importance for global food security [1]. China is the largest rice-producing country, accounting for 28% of global rice production [2]. Middle-and double-season rice are the dominant rice systems in central China [3]. However, intensive cropping systems result in adverse impacts on the environment [4] as well as in higher production costs. During recent decades, double-season rice systems in central China have progressively shifted into single-season rice mainly due to high costs and an increase in labor shortages [5]. Hence, there is Remote Sens. 2020, 12, 3400 2 of 13 a growing interest in identifying an alternative rice system that can increase the productivity while reducing the environmental impact and the labor cost of existing cultivation practices.
Ratooning is a management option in rice systems that allows for rice regrowth leading to a second or ratoon crop harvest following that of the main crop [6]. Compared to a double-season rice system, the ratoon rice (RR) system does not require the additional labor for transplanting a second rice crop. Ratoon crop yield can be as high as 50% of that produced from the main crop [7]. The yield of RR is generally higher than single-season rice and lower than double-season rice [5]. Worldwide, ratooning may help adapt rice production to climate uncertainty [8]. RR is currently practiced in many countries including China, India, Japan, Philippines, Puerto Rico, Brazil, Thailand, and the U.S., but only in limited areas [9,10]. In China, RR has been traditionally practiced widely since the 1950s. In recent years, RR planted area increased rapidly in central China, taking advantage of the development of new rice cultivars with high ratooning ability, the advancements in crop and water management, and mechanical harvesting techniques [5]. Because of the increasing planted area of RR and the need to balance agricultural labor costs and food security, the sustainable and effective management of rice production in central China requires accurate mapping of the RR planted area. Given the fragmented crop fields in central China, remote sensing provides the most practical and cost-effective way to map RR in this region.
Many efforts have been made to map the extent or area of rice paddy fields using different remotely sensed data including optical and microwave remote sensing images [11][12][13][14][15] and various algorithms (e.g., supervised classifiers, unsupervised classifiers, and knowledge-based approaches) [16][17][18][19]. A thorough review of regional to global paddy rice mapping methods was provided by Dong and Xiao [20]. As rice spends a portion of its phenological cycle submerged in water, most remote sensing methods map it by extracting the flooded paddy area during the initial growth stages of rice. Xiao et al. [21] found that the flooding and transplanting signals can be extracted based on the relationship between the normalized difference vegetation index (NDVI) or the enhanced vegetation index (EVI) and the land surface water index (LSWI), in which the LSWI values are temporarily greater than NDVI or EVI values during these phases. This phenology-based method has been widely used to map paddy rice at regional and national scales [12,22,23]. Improvements have been made to adapt the original method to different study areas and increase the accuracy of paddy rice mapping [20,24].
Due to the complexity of the rice cropping system in the Southeast Asia, a few studies have focused on differentiating the single-, double-, and triple-season rice in this region with various methods. Rice paddy fields in the Mekong River Delta were mapped by adapting the phenology-based method and using Landsat time stacks by Kontgis et al. [24]. They differentiated between single-, double-, and triple-season rice paddies by counting the peaks in EVI trajectories. In Guan et al. [25]'s study, rice cropping systems in Vietnam were mapped based on the similarity between the NDVI time series of every pixel and the standard NDVI time series of different rice cropping systems. The single-and double-season rice in Hanoi, Vietnam were mapped using a Random Forest Algorithm and a time series of Sentinel-1 SAR imagery by Lasko et al. [26]. The rice intensity in the Senegal River Valley was analyzed using the PhenoRice algorithm and MODIS time series by Busetto et al. [15].
The current methods for mapping rice cropping systems mainly used the flooding signals, the peaks of EVI/NDVI, or the similarities between the reference times series of EVI/NDVI and those of the targets. However, these methods are not suitable for mapping RR fields. The challenge of differentiating the RR from the other rice cropping system is that the rice stubble left after the harvest of the main crop may minimize the flooding signal of the second crop in the remotely sensed images. In addition, the NDVI and EVI time series of RR may resemble those of double-season rice, both of which have two phenological cycles. RR would be easily classified as double-season rice with current methods developed for mapping rice cropping systems.
The phenology of RR is slightly different from the other rice systems. Its ripening stage starts later than that of the early-season rice and earlier than that of the middle-season rice [5]. We hypothesized that the spectral characteristics of rice at the ripening phase can be used to differentiate RR from Remote Sens. 2020, 12, 3400 3 of 13 the other rice systems. Therefore, the objective of this study was to develop a phenology-based algorithm to map the RR planted area in Hubei Province, Central China in 2018. The novelty of this phenology-based algorithm was to identify RR fields with two contrasting vegetation indices: the newly-developed yellowness index (YI) to detect the ripening phase and the EVI to detect the green-up of the second-season crop. Since the crop fields in central China are highly fragmented, we used Sentinel-2A and 2B time stacks at a 10-m spatial resolution to obtain a detailed and accurate map of an RR planted area.

Study Area
This study was conducted in Hubei Province, Central China, covering 185,900 km 2 ( Figure 1). The Jianghan Plain takes up most of central and southern Hubei, while western Hubei is more mountainous. Hubei has a humid subtropical climate, with precipitation following a monsoon pattern. In 2018, the mean temperature was 17.1 • C and the mean precipitation was 1119 mm [27].
Hubei is one of the most important rice producing provinces in China, with a total annual rice production of 17 Mt, representing approximately 8% of the total national rice production [28]. A gradual transition from middle-and double-season rice to RR has occurred in Hubei [29]. With the increasing shortage of labor and the improved RR cultivars and management strategies, additional planting areas of middle-and double-season rice are expected to be converted to RR.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 14 RR from the other rice systems. Therefore, the objective of this study was to develop a phenologybased algorithm to map the RR planted area in Hubei Province, Central China in 2018. The novelty of this phenology-based algorithm was to identify RR fields with two contrasting vegetation indices: the newly-developed yellowness index (YI) to detect the ripening phase and the EVI to detect the green-up of the second-season crop. Since the crop fields in central China are highly fragmented, we used Sentinel-2A and 2B time stacks at a 10-m spatial resolution to obtain a detailed and accurate map of an RR planted area.

Study Area
This study was conducted in Hubei Province, Central China, covering 185,900 km 2 ( Figure 1). The Jianghan Plain takes up most of central and southern Hubei, while western Hubei is more mountainous. Hubei has a humid subtropical climate, with precipitation following a monsoon pattern. In 2018, the mean temperature was 17.1 °C and the mean precipitation was 1119 mm [27].
Hubei is one of the most important rice producing provinces in China, with a total annual rice production of 17 Mt, representing approximately 8% of the total national rice production [28]. A gradual transition from middle-and double-season rice to RR has occurred in Hubei [29]. With the increasing shortage of labor and the improved RR cultivars and management strategies, additional planting areas of middle-and double-season rice are expected to be converted to RR.

Sentinel-2 Data and Preprocessing
We used all available Level 1 scenes of the 2018 growing season (March to November) with cloud cover of less than 80% in order to capture the flooding signals of the complex rice cropping systems

Sentinel-2 Data and Preprocessing
We used all available Level 1 scenes of the 2018 growing season (March to November) with cloud cover of less than 80% in order to capture the flooding signals of the complex rice cropping systems in the study area. The preprocessing involved the default atmospheric and topographic correction using Remote Sens. 2020, 12, 3400 4 of 13 Sen2Cor in the Sentinel-2 toolbox in SNAP (vers.6.0, ESA, 2016). This study used bands 2, 3, 4, 8, and 11. The central wavelength, bandwidths, and spatial resolutions of these bands are shown in Table 1. Band 11 was resampled to 10-m resolution, using the cubic convolution method. Cloud-cover and cloud-shadow pixels were identified as those pixels with blue reflectance higher than 0.25 [18]. In addition to the reflectance data, we calculated the EVI and the LSWI with the following equations and spectral bands: where B4, B8, and B11 are the surface reflectance of Band 4, Band 8, and Band 11 in the Sentinel-2A and 2B sensors.
To remove the small fluctuations in the EVI times series, we used the improved Savitzky-Golay filter developed by Chen et al. [30]. Since the LSWI is sensitive to the change in vegetation water and soil moisture [31], the LSWI time series was not smoothed in this study.

Reference Data
Reference data of RR, non-RR paddy rice, and non-rice land cover types were obtained from field surveys and high-resolution satellite images. We visited 60 villages and collected ground truths data with a portable GPS (Garmin MTK 6735, Garmin Ltd., Switzerland) and an unmanned aerial vehicle (Mavic Pro, SZ DJI Technology Co., Ltd., Shenzhen, China) during the growing season of 2018 and 2019. Additionally, we visually identified non-RR paddy rice samples and non-rice land cover samples from two villages in GF-2 multispectral images. We collected 984 reference data, including 355 RR samples, 412 non-RR paddy rice samples, and 217 non-rice land cover types. Figure 1 shows the spatial distribution of the 60 villages in total where reference data were collected in the study area.

FROM-GLC10
FROM-GLC10 is the 10-m resolution global land cover map based on 2017 Sentinel-2 data [32]. The dataset includes ten land cover types, namely cropland, forest, grassland, shrubland, wetland, water, tundra, impervious surface, barren, and snow/ice. The overall accuracy of the global land cover map is 72.76%. We downloaded FROM-GLC10 of the study area from http://data.ess.tsinghua.edu.cn and the cropland in Hubei province derived was used to screen the paddy rice field pixels.

Digital Elevation Model (DEM)
DEM data (GDEMDEM) with a spatial resolution of 30 m was downloaded from Geospatial Data Cloud (http://www.gscloud.cn/search). The 30-m DEM data was resampled to 10 m using the nearest neighbor resampling method. Slope was calculated based on the 10-m DEM data in ArcGIS10.2 (Environmental Systems Research Institute, Inc., Red Lands, CA, USA).

Rice Planting Area from Rural Statistical Yearbook
Hubei Rural Statistical Yearbook publishes annual reports on the planting areas of the main crops in each city. The statistical yearbook that reports the rice planting area in 2018 was downloaded from https://navi.cnki.net/knavi/Yearbook.html. The rice planting area of 17 prefectural-level cities was used to validate the rice planting area derived from this study at the city level.

Methods
The RR planning area was mapped based on the paddy rice field map, which was generated by detecting the flooding signal of cropland during the growing season. RR was identified based on the differences in the timing of ripening phase from that of the middle-and double-season rice system in the study area. We proposed a YI to capture the ripening signal. The EVI values after the ripening period were used to remove the confusing pixels that had stubbles from the early-season rice left in the field. The accuracy of the RR map was evaluated with the reference data. An overview of the workflow is given in Figure 2. DEM data (GDEMDEM) with a spatial resolution of 30 m was downloaded from Geospatial Data Cloud (http://www.gscloud.cn/search). The 30-m DEM data was resampled to 10 m using the nearest neighbor resampling method. Slope was calculated based on the 10-m DEM data in ArcGIS10.2 (Environmental Systems Research Institute, Inc., Red Lands, CA, USA).

Rice Planting Area from Rural Statistical Yearbook
Hubei Rural Statistical Yearbook publishes annual reports on the planting areas of the main crops in each city. The statistical yearbook that reports the rice planting area in 2018 was downloaded from https://navi.cnki.net/knavi/Yearbook.html. The rice planting area of 17 prefectural-level cities was used to validate the rice planting area derived from this study at the city level.

Methods
The RR planning area was mapped based on the paddy rice field map, which was generated by detecting the flooding signal of cropland during the growing season. RR was identified based on the differences in the timing of ripening phase from that of the middle-and double-season rice system in the study area. We proposed a YI to capture the ripening signal. The EVI values after the ripening period were used to remove the confusing pixels that had stubbles from the early-season rice left in the field. The accuracy of the RR map was evaluated with the reference data. An overview of the workflow is given in Figure 2.

Mapping Paddy Rice Field
Rice has three distinct growing phases, including the sowing-transplanting phase, the growing phase, and the after-harvest phase [33]. Each phase has a distinct spectral signature that can be explored for mapping paddy rice with remote sensing techniques. Especially, during the sowingtransplanting phase, paddy rice fields are flooded and canopies are open. In this study, the phenology-based method proposed by Xiao et al. [21] was adopted to map paddy rice fields and the

Mapping Paddy Rice Field
Rice has three distinct growing phases, including the sowing-transplanting phase, the growing phase, and the after-harvest phase [33]. Each phase has a distinct spectral signature that can be explored for mapping paddy rice with remote sensing techniques. Especially, during the sowing-transplanting phase, paddy rice fields are flooded and canopies are open. In this study, the phenology-based method proposed by Xiao et al. [21] was adopted to map paddy rice fields and the difference between the EVI and the LSWI (EVI-LSWI) was calculated for each scene to detect the flooding signal of paddy rice fields. The phenology-based method detects flooding and transplanting signals of paddy rice using the criteria LSWI > EVI or LSWI > NDVI. The study area has early-, middle-, and late-season rice systems. Rice is transplanted no earlier than mid-March and no later than mid-August. Thus, a pixel was identified as a potential paddy rice field as long as the value of (EVI-LSWI) was smaller than 0 at any point from March 15th to August 15th. There are many small ponds in the agricultural fields in the central and southern Hubei Province. These ponds are mainly used to grow aquatic plants or for fish farming. The ponds that grow aquatic plants exhibit a flooding signal before aquatic plants thrive, and their EVI variations are thus similar to those of paddy rice. Comparing the time series of the EVI and the LSWI of ponds that grow aquatic plants and the paddy rice fields, we found that the EVI over the ponds increased sharply in June and reached higher EVI values than over paddy rice fields. The average June EVI values of the pond samples were calculated and used as the threshold to differentiate the paddy rice fields from ponds.
Moreover, as the paddy rice is mainly planted in the plain area due to flooding requirements during the transplanting stage, land with a slope higher than 7 • was excluded from this study. The threshold to define the sloping land was derived from the slope of the 95% quantile of the ground paddy rice samples.
Therefore, a pixel was classified as the paddy rice field only if these four conditions were met: where CROP == 1 means the pixel is identified as cropland in FROM-GLC10, SLOPE is the slope calculated based on the 10-m DEM data, EVI i and LSWI i are the EVI and LSWI values on any day during the study period, and EVI June is the average of the EVI values in June.

Mapping RR
As rice matures, its canopy turns 'golden' and the reflectance of the green and red bands increases. To capture the maturity signal of RR and differentiate RR from the middle-and double-season rice systems, we proposed the YI calculated as the contrast between blue light and the combination of red and green light: where B2, B3, and B4 are the surface reflectance of Band 2, Band 3, and Band 4 in the Sentinel-2A and 2B sensors.
Since the timing of maturity may vary across the study area, usually between July 20th and August 10th, we firstly calculated the mean YI (YI mean ) during July 20th to August 10th for each pixel. To enhance the contrast between the mature RR and the non-RR, the Otsu's method was employed to classify the new dataset into two types [34]. The pixels with the higher YI mean were considered to be potential RR pixels, and they were assigned with the maximal YI (YI max ) during July 20th to August 10th. The pixels with the lower YI mean were assigned with the minimal YI (YI min ) during July 20th-August 10th. Thus, a new dataset with the enhanced YI values was composed.
Based on this new dataset, pixels with YI values higher than the threshold (Th MI ) were identified as potential RR. The Th MI was determined based on the YI max of the 50% samples randomly selected from the RR reference data: where YI max is the average YI max of RR samples and YI std is the standard deviation of the YI in the RR samples. However, in the study area, stubbles of the early-season rice were usually left in the field, and their spectra confused with those of the mature RR. Therefore, pixels identified as potential RR fields also needed to meet the criterion that the maximum EVI (EVI max ) after one month of the ripening phase of Remote Sens. 2020, 12, 3400 7 of 13 RR, namely from August 11th to September 11th, had to be higher than a threshold (Th EVI ). The Th EVI was determined based on the EVI max of 50 % of samples randomly selected from the RR reference data: where EVI max is the average of EVI max , and EVI std is the standard deviation of the EVI max of RR samples during August 11th and September 11th. Moreover, to differentiate the RR from the double-season rice, the (EVI-LSWI) should be higher than 0 during the period from July 20th to August 10th, because late-season rice was usually transplanted in the field and thus had a flooding signal while the RR was maturing.

Validation
The rice planting area extracted from the Hubei Rural Statistical Yearbook was used to validate the rice planting area at the city level. Unfortunately, the Hubei Rural Statistical Yearbook does not provide the statistics of RR planting area. A confusion matrix, the overall accuracy, and the Kappa coefficient based on the ground reference data were calculated to evaluate the accuracy of the RR map.

Phenological Characteristics of RR
The time series of the smoothed EVI of middle-season rice, double-season rice, and RR showed that the ripening phase of RR provided a good opportunity to differentiate RR from the other rice systems. RR and middle-season rice can be easily differentiated in this window as middle-season rice was 'green' while RR was 'yellow'. Both RR and double-season rice had two peaks in the time series of the EVI, but the early-season rice reached the peak earlier than RR and also matured earlier than RR (Figure 3a). While RR was maturing, early-season rice was harvested and late-season rice was transplanted in the field. Thus, in the ripening phase of RR, the (EVI-LSWI) was smaller than 0 over double season rice (Figure 3b). However, in the study area, stubbles of the early-season rice were usually left in the field, and their spectra confused with those of the mature RR. Therefore, pixels identified as potential RR fields also needed to meet the criterion that the maximum EVI (EVImax) after one month of the ripening phase of RR, namely from August 11th to September 11th, had to be higher than a threshold (ThEVI). The ThEVI was determined based on the EVImax of 50 % of samples randomly selected from the RR reference data: where is the average of EVImax, and EVIstd is the standard deviation of the EVImax of RR samples during August 11th and September 11th.
Moreover, to differentiate the RR from the double-season rice, the (EVI-LSWI) should be higher than 0 during the period from July 20th to August 10th, because late-season rice was usually transplanted in the field and thus had a flooding signal while the RR was maturing.

Validation
The rice planting area extracted from the Hubei Rural Statistical Yearbook was used to validate the rice planting area at the city level. Unfortunately, the Hubei Rural Statistical Yearbook does not provide the statistics of RR planting area. A confusion matrix, the overall accuracy, and the Kappa coefficient based on the ground reference data were calculated to evaluate the accuracy of the RR map.

Phenological Characteristics of RR
The time series of the smoothed EVI of middle-season rice, double-season rice, and RR showed that the ripening phase of RR provided a good opportunity to differentiate RR from the other rice systems. RR and middle-season rice can be easily differentiated in this window as middle-season rice was 'green' while RR was 'yellow'. Both RR and double-season rice had two peaks in the time series of the EVI, but the early-season rice reached the peak earlier than RR and also matured earlier than RR (Figure 3a). While RR was maturing, early-season rice was harvested and late-season rice was transplanted in the field. Thus, in the ripening phase of RR, the (EVI-LSWI) was smaller than 0 over double season rice (Figure 3b).

Identification of RR
The key in differentiating RR from other rice cropping systems is its different timing of maturity. To capture the maturity signal, we proposed a YI and composited the enhanced YI image to identify the potential RR. Figure 4 shows an example of a true-color Sentinel-2A image of a cropland dominated by RR fields, its YI mean image, the enhanced YI image, and the segmentation of the potential RR Remote Sens. 2020, 12, 3400 8 of 13 fields. From the visual inspection, the YI mean showed the degree of yellowness of the fields, while the enhanced YI image increased the contrast between the mature RR and the other rice systems or crops. However, potential RR fields still had miss-classified pixels mainly due to the stubbles left in the field after the harvest of the early-season rice. The increase of the EVI after the maturity was a requirement to determine the RR fields.

Identification of RR
The key in differentiating RR from other rice cropping systems is its different timing of maturity. To capture the maturity signal, we proposed a YI and composited the enhanced YI image to identify the potential RR. Figure 4 shows an example of a true-color Sentinel-2A image of a cropland dominated by RR fields, its YI mean image, the enhanced YI image, and the segmentation of the potential RR fields. From the visual inspection, the YI mean showed the degree of yellowness of the fields, while the enhanced YI image increased the contrast between the mature RR and the other rice systems or crops. However, potential RR fields still had miss-classified pixels mainly due to the stubbles left in the field after the harvest of the early-season rice. The increase of the EVI after the maturity was a requirement to determine the RR fields.

RR Map and Accuracy Assessment
The paddy rice field map was validated with both the statistical data at the city level and the reference data collected in the field. Figure 5 illustrates that the predicted rice planting area of 17 prefectural-level cities was close to the rice planting area extracted from the Hubei Rural Statistical Yearbook. The total area of the paddy rice fields predicted in this study was 22,196.20 km 2 , slightly lower than the area from the statistical data (22,372.21 km 2 ), mainly because sloping crop pixels were excluded in this study. The overall accuracy of paddy rice fields based on the field data was 0.76. Table 2 demonstrates the confusion matrix of RR, non-RR rice systems, and non-rice land cover types based on the reference data collected in the field. The RR category had a producer's accuracy of 0.95 and user's accuracy of 0.76. The accuracy was mainly influenced by the misclassification of non-RR rice systems as non-rice land cover types. If the accuracy was adjusted to the identification of RR from the correctly classified paddy rice fields, the overall accuracy increased to 0.91.
The RR and paddy rice field map illustrated that rice was the major summer crop in Hubei Province (Figure 6), and RR was mainly distributed in the low alluvial plain in central and southern Hubei Province. The total planting area of RR in 2018 was 2225.4 km 2 , accounting for 10.03% of the total area of paddy rice fields.

RR Map and Accuracy Assessment
The paddy rice field map was validated with both the statistical data at the city level and the reference data collected in the field. Figure 5 illustrates that the predicted rice planting area of 17 prefectural-level cities was close to the rice planting area extracted from the Hubei Rural Statistical Yearbook. The total area of the paddy rice fields predicted in this study was 22,196.20 km 2 , slightly lower than the area from the statistical data (22,372.21 km 2 ), mainly because sloping crop pixels were excluded in this study. The overall accuracy of paddy rice fields based on the field data was 0.76. Table 2 demonstrates the confusion matrix of RR, non-RR rice systems, and non-rice land cover types based on the reference data collected in the field. The RR category had a producer's accuracy of 0.95 and user's accuracy of 0.76. The accuracy was mainly influenced by the misclassification of non-RR rice systems as non-rice land cover types. If the accuracy was adjusted to the identification of RR from the correctly classified paddy rice fields, the overall accuracy increased to 0.91.
The RR and paddy rice field map illustrated that rice was the major summer crop in Hubei Province (Figure 6), and RR was mainly distributed in the low alluvial plain in central and southern Hubei Province. The total planting area of RR in 2018 was 2225.4 km 2 , accounting for 10.03% of the total area of paddy rice fields.

Advantages of the Dense Time Stacks of Sentinel-2 Images
This study indicated the feasibility and reliability of identifying RR from the other rice systems in the plain areas using all available Sentinel-2 images (clouds cover <80%). The high spatial and temporal resolution of Sentinel-2 enabled the capture of important phenological information to identify RR in fragmented croplands. As the key in differentiating RR from the non-RR rice systems was the timing of its maturity, the data availability during the period July 20th-August 10th had a great impact on the accuracy of the RR map. Luckily, in 2018 at least three images were available within this window in most of the paddy rice planting areas. An alternative to improve the data availability is to incorporate Landsat 8 images. Efforts have been made to successfully map crops or crop intensity by integrating Sentinel-2 and Landsat images [31,35]. However, due to the different spectral bands and spatial resolution of these sensors, the methods to integrate Sentinel-2 and Landsat images and the applications of the integrated data to mapping rice systems still need further investigations.

Identification of RR Using the YI
To our knowledge, this is the first study to map RR. As RR is a rice cropping system, we compared our study to previous studies regarding the methods and the accuracy of classification of rice cropping systems. Previous studies identified different rice cropping systems either using the number of EVI/NDVI peaks [24,36,37] or the similarities in NDVI/EVI time series between the reference and the target [13,25]. Due to the frequent cloud cover in remotely sensed images in the rice planting areas, the classification accuracy derived from these methods were dependent on the smoothing of NDVI/EVI time series [36,37]. Instead of using one index (NDVI/EVI) to identify rice cropping systems, this study highlights the advantages of using two contrasting indices. The proposed YI was used to detect the ripening phase for identifying potential RR fields, and the EVI was used to detect the green-up of the second season of RR for eliminating the false signals from the stubbles left in the field. The applications of two indices for the detection of different phenological phases may provide inspirations in decreasing the misclassifications in rice cropping systems, which are caused by missing observations during the significant phenological phase (e.g., the heading phase).
The proposed YI was calculated with reflectance of the blue, green, and red bands. The values in the blue band have been known to be sensitive to the unstable weather conditions [38]. In this study, having performed the atmospheric correction of the data, we did not see a substantial influence of the variability in the blue band on the identification of the RR, possibly because the blue band was used only during a very short period (July 20th-August 10th). As the middle and lower reaches of Yangzi River has multiple rice systems, the proposed YI could also be combined with a greenness index, such as EVI, to map middle-, early-, and late-season rice planting area.
The overall accuracy for identifying RR, non-RR rice systems, and non-rice land cover types in this study was 0.76, comparable to previous studies of rice cropping systems. Rice cropping systems in Mekong Delta, Vietnam were mapped using Moderate Resolution Imaging Spectroradiometer (MODIS) data by Son et al. [37], and they achieved the overall accuracy of 80.6% for 2006 and 85.5% for 2012. Kontgis et al. [24]'s study achieved the overall accuracies of 71% and 69% for the circa 2000 and circa 2010 maps of rice cropping systems in the Mekon River Delta, Vietnam, using dense time stacks of Landsat data. In Guan et al. [25]'s study, rice cropping systems were mapped in Vietnam using an NDVI-based time series similarity measurement, with the R 2 between their results and the statistical data of 0.809.

Sources of Errors in the RR Map in Hubei Province
There were two types of misclassifications of RR: (1) RR misclassified as the non-RR rice systems and (2) RR misclassified as the non-rice land cover types. The overall accuracy of the RR map given that the paddy rice fields were correctly identified was 0.91, much higher than the overall accuracy including both RR and non-RR rice fields. This implied that the failure to identify RR was mainly due to the misclassification of paddy rice fields in the first place. As pointed out by previous studies that mapped rice with a phenology-based method, the lack of data during the transplanting period causes the flooding signal in the time series of remotely sensed images to be missed [18,20,24,25]. Mixed pixels of paddy rice fields and non-rice fields (roads, trees, and vegetable fields) may also cause the flooding signal to be missed. Additionally, the irrigation or flooding events in non-rice land cover types also affect the classification results. Although our algorithm includes a rule to reduce the confusion from ponds of aquatic plants, some of the ponds were still classified as the paddy rice fields because of their similar temporal properties of EVI.
If the paddy rice fields were correctly classified, the misclassification of RR was mainly due to the variable transplanting timing across the province. Differences in weather conditions between north and south and between the mountain and the plain led to the variations in seedling transplantation timing across Hubei province. Late transplanting of the early-season rice may result in an overlap of the maturity timing with that of RR. The stubbles from the early-season rice were generally left in the fields, and if the late-season rice were not planted due to the labor shortages, the weeds thrived in late August and resulted in an EVI peak. Therefore, fields with the early-and double-season rice were easily misclassified as RR when they presented similar phenologies in the remotely sensed data.

Conclusions
This study developed a phenology-based algorithm to map RR at a 10-m resolution in Hubei Province, Central China, using dense time stacks of Sentinel-2 images (cloud cover<80%) in 2018. The novelty of this phenology-based algorithm was to identify RR with two contrasting indices, which included the newly-developed YI to detect the ripening phase and the EVI to detect the green-up of the second-season crop. The total planting area of RR in 2018 was 2225.4 km 2 , accounting for 10.03% of the total area of paddy rice fields. The overall accuracy including the RR, non-RR rice fields, and non-rice land cover types was 0.76. The adjusted overall accuracy for RR and non-RR was 0.91. The algorithm proposed by this study provided accurate and up-to-date maps of RR fields, which may facilitate the prediction of grain yield and the evaluation of the environmental impact of rice production. Furthermore, monitoring the shift from non-RR fields to RR fields could be of interest for agricultural policy-making and the socio-economic impacts assessment.