Phenology-Based Rice Paddy Mapping Using Multi-Source Satellite Imagery and a Fusion Algorithm Applied to the Poyang Lake Plain, Southern China

: Accurate information about the spatiotemporal patterns of rice paddies is essential for the assessment of food security, management of agricultural resources, and sustainability of ecosystems. However, accurate spatial datasets of rice paddy ﬁelds and multi-cropping at ﬁne resolution are still lacking. Landsat observation is the primary source of remote sensing data that has continuously mapped regional rice paddy ﬁelds at a 30-m spatial resolution since the 1980s. However, Landsat data used for rice paddy studies reveals some challenges, especially data quality issues (e.g., cloud cover). Here, we present an algorithm that integrates time-series Landsat and MODIS (Moderate-resolution Imaging Spectroradiometer) images with a phenology-based approach (ILMP) to map rice paddy planting ﬁelds and multi-cropping patterns. First, a fusion of MODIS and Landsat data was used to reduce the cloud contamination, which added more information to the Landsat time series data. Second, the unique biophysical features of rice paddies during the ﬂooding and open-canopy periods (which can be captured by the dynamics of the vegetation indices) were used to identify rice paddy regions as well as those of multi-cropping. This algorithm was tested for 2015 in Nanchang County, which is located on the Poyang Lake plain in southern China. We evaluated the resultant map of the rice paddy and multi-cropping systems using ground-truth data and Google Earth images. The overall accuracy and kappa coe ﬃ cient of the rice paddy planting areas were 93.66% and 0.85, respectively. The overall accuracy and kappa coe ﬃ cient of the multi-cropping regions were 92.95% and 0.89, respectively. In addition, our algorithm was more capable of capturing detailed information about areas with fragmented cropland than that of the National Land Cover Dataset (NLCD) from 2015. These results demonstrated the great potential of our algorithm for mapping rice paddy ﬁelds and using the multi-cropping index in complex landscapes in southern China.


Introduction
Rice paddy fields account now for more than 12% of the global cultivated land area [1] and yield food for approximately half of the global population, particularly in monsoon Asia [2,3]. With the Southern China is one of the most famous rice production bases over the world. However, accurate mapping of the rice paddy fields at fine resolutions in southern China is still a challenge [13,26]. First, the cropland-filled are in southern China is quite small, and the landscape is heterogeneous. Meanwhile, several crops grow during the growing season of rice, which forms a very complex cropping structure in areas with fragmented and patchy fields, and the rice paddy areas are easily confused with other types of vegetation cover [23,27]. Second, the multiple cropping systems superimposed with prevailing cloud cover make it difficult to characterize a crop's phenological metrics [28]. Furthermore, the rice-wetland coexistence areas also present a challenge [22]. Both rice paddies and wetlands have similar spectral signals within growing seasons, which tends to result in misclassification of natural wetlands as rice paddies [10,29]. To overcome the limitations mentioned above, we present a feasible strategy of integrating a spatiotemporal fusion approach and phenology-based approach to map the rice paddy planting fields and multi-cropping patterns in southern China. The objectives of this study were to: (1) develop a robust phenology-assisted rice paddy planting area and multi-cropping system through the fused Landsat-MODIS time series images, and (2) assess the fusion-and phenology-based approaches for mapping rice paddies and cropping intensity in the fragmented areas of diverse land-cover types in southern China.

Study Area
The study area is located in Nanchang County, which is a part of the Poyang Lake plain in southern China (Figure 1). This plain is one of the famous (and important) grain production bases in China and is dominated by traditional double cropping systems. The study area neighbors Poyang Lake and is largely flat. It occupies an area of approximately 1800 km 2 and belongs to the subtropical climate type. Major annual crops in the study area include single-and double-cropping rice, rapeseed oil, oilseed rape, and lotus ponds. Meanwhile, smallholders dominate rice paddy agriculture in Nanchang County, resulting in a mixed agricultural landscape with small fields. According to the Jiangxi Statistical Yearbooks, the area of sown rice paddy accounts for above 90% of the total area of sown crops [30].
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 19 Southern China is one of the most famous rice production bases over the world. However, accurate mapping of the rice paddy fields at fine resolutions in southern China is still a challenge [13,26]. First, the cropland-filled are in southern China is quite small, and the landscape is heterogeneous. Meanwhile, several crops grow during the growing season of rice, which forms a very complex cropping structure in areas with fragmented and patchy fields, and the rice paddy areas are easily confused with other types of vegetation cover [23,27]. Second, the multiple cropping systems superimposed with prevailing cloud cover make it difficult to characterize a crop's phenological metrics [28]. Furthermore, the rice-wetland coexistence areas also present a challenge [22]. Both rice paddies and wetlands have similar spectral signals within growing seasons, which tends to result in misclassification of natural wetlands as rice paddies [10,29]. To overcome the limitations mentioned above, we present a feasible strategy of integrating a spatiotemporal fusion approach and phenology-based approach to map the rice paddy planting fields and multi-cropping patterns in southern China. The objectives of this study were to: (1) develop a robust phenologyassisted rice paddy planting area and multi-cropping system through the fused Landsat-MODIS time series images, and (2) assess the fusion-and phenology-based approaches for mapping rice paddies and cropping intensity in the fragmented areas of diverse land-cover types in southern China.

Study Area
The study area is located in Nanchang County, which is a part of the Poyang Lake plain in southern China (Figure 1). This plain is one of the famous (and important) grain production bases in China and is dominated by traditional double cropping systems. The study area neighbors Poyang Lake and is largely flat. It occupies an area of approximately 1800 km 2 and belongs to the subtropical climate type. Major annual crops in the study area include single-and double-cropping rice, rapeseed oil, vegetable gardens, and lotus ponds. Meanwhile, smallholders dominate rice paddy agriculture in Nanchang County, resulting in a mixed agricultural landscape with small fields. According to the Jiangxi Statistical Yearbooks, the area of sown rice paddy accounts for above 90% of the total area of sown crops [30].   Landsat 7 (ETM+) and Landsat 8 (OLI) images that cover our chosen study area (path/row of Landsat imagery: 121/40) for 2015 were obtained from the USGS Earth Explorer data portal (http://earthexplorer.usgs.gov/). In total, 46 Landsat images were collected, which include 23 Landsat ETM+ images and 23 Landsat OLI images. These images were produced with terrain and atmospheric corrections, and the land surface reflectances were derived. Cloud cover and cloud shadow detection are critical concerns for optical remote sensing images processing. Therefore, we applied the Fmask software to generate cloud and cloud shadow layers to each processed image [31]. The SLC-off strips within the Landsat ETM images were classified as "no-observation" using this software. Once the clouds, cloud shadows, and SLC-off strips were excluded, we end up with high quality Landsat images suitable for our purposes. We also obtained eight-day 500-m MODIS reflectance images (MOD09A1) for 2015 from the USGS Earth Resources Observation and Science (EROS) data center (https://lpdaac.usgs.gov/).

Validation Data
Field surveys within the study area were performed on the period of May-August from 2015 to 2018. The survey questions included information on the land-cover type, crop type, and crop calendar of each site in 2015. The sampling distance between each site was approximately 2 km, and the width and length of each field site were larger than 100 m. Five classes of land-cover type were used: rice paddy (single-and double-cropping rice), wetland, lotus, dry land, and others. Water bodies, built-up or barren land, and forest were classified as "others." In the wetland locations, we did our surveying on the road running through the large wetland area. For the sites of croplands and other land-cover types, we travelled more than 60 m from the border in each direction in order to take accurate georeferenced photographs. The total number of field survey sites was 200 ( Figure 1).
Previous studies have indicated that combining field photographs with Google Earth images is a reliable technique for the validation of land cover classification [32]. Thus, we combined the field survey data with high-resolution images from Google Earth to generate validation points. A total of 2040 pixels were generated for this validation. Even though some field survey data were collected from 2016 to 2018, the survey in 2015 provided a reference for the validation points generated from the high-resolution images taken from Google Earth. The images we used in this study were mostly from April to July 2015. We divided all validation points into three categories: single-cropping rice, double-cropping rice, and "others" (e.g., wetland, forest, water pools, etc.).

Other Land Cover Datasets Used for Comparison
The National Land Cover Dataset (NLCD) for 2015 was obtained from the Resource and Environment Data Center, Chinese Academy of Sciences [33]. The NLCD datasets at a scale of 1:100,000 were derived from Landsat TM/ETM+/OLI images using the human-computer interactive interpretation method. The NLCD divides the land cover into six primary categories with 25 subcategories. The cropland includes "rice paddy land" and "dry cropland" categories, and the rice paddy land is primarily distributed in southern China.

Overview
An overview of the methodology for the current study is shown in Figure 2. The fusion-and phenology-based approaches for mapping rice paddies and multi-cropping lands included three core modules: (1) the remote sensing image preprocessing module (Section 3.2), (2) the pixel-and phenology-based rice paddy and multi-cropping mapping module (Section 3.3), and (3) the accuracy

Calculation of Spectral Indices
The reflectance data from the Landsat-ETM, Landsat OLI, and MODIS (MOD09A1) data were used to calculate the normalized difference vegetation index (NDVI) [34], enhanced vegetation index (EVI) [35], and land surface water index (LSWI) [10], using the land surface reflectance values of the blue ( ), red ( ), NIR ( ), and SWIR ( , 1570-1650 nm) bands. These indices have proven to be helpful in studies of vegetation canopies, water content, biomass, and phenology.

Calculation of Spectral Indices
The reflectance data from the Landsat-ETM, Landsat OLI, and MODIS (MOD09A1) data were used to calculate the normalized difference vegetation index (NDVI) [34], enhanced vegetation index (EVI) [35], and land surface water index (LSWI) [10], using the land surface reflectance values of the blue (ρ blue ), red (ρ red ), NIR (ρ NIR ), and SWIR (ρ SWIR , 1570-1650 nm) bands. These indices have proven to be helpful in studies of vegetation canopies, water content, biomass, and phenology.
Remote Sens. 2020, 12, 1022 6 of 18 Even after steps are performed to improve the data quality (e.g., maximum value composition method (MVC)) of the MODIS (MOD09A1) data [36]), atmospheric contamination in the images may still exist. We therefore employed a Savitzky-Golay filtering procedure to further reduce contamination by clouds, and to reconstruct the vegetation indices' (MOD09A1) time series [37].

Fusion of Landsat OLI and MODIS Images
Two sets of vegetation indices were generated from Landsat ETM+/OLI images; however, some images had serious cloud contamination, and SLC-off strips were not available. This made the temporal interval long in some months, especially during the growing season ( Figure 3). This phenomenon can result in incomplete crop phenology and can limit the accuracy of rice paddy mapping. To increase the available data for the growing season, this study applied the ESTARFM approach to disaggregate the 500-m vegetation index images (MOD09A1) to a 30-m pixel size to coincide with the vegetation indices from the Landsat OLI images [16,20,38,39]. Recently, several studies used the ESTARFM approach to generate synthetic imagery to track the phenological changes and found that fused Landsat-MODIS data could deduce the land surface phenology at sub-field scales more accurately than using the MODIS data alone [16,20,39,40]. method (MVC)) of the MODIS (MOD09A1) data [36]), atmospheric contamination in the images may still exist. We therefore employed a Savitzky-Golay filtering procedure to further reduce contamination by clouds, and to reconstruct the vegetation indices' (MOD09A1) time series [37].

Fusion of Landsat OLI and MODIS images
Two sets of vegetation indices were generated from Landsat ETM+/OLI images; however, some images had serious cloud contamination, and SLC-off strips were not available. This made the temporal interval long in some months, especially during the growing season ( Figure 3). This phenomenon can result in incomplete crop phenology and can limit the accuracy of rice paddy mapping. To increase the available data for the growing season, this study applied the ESTARFM approach to disaggregate the 500-m vegetation index images (MOD09A1) to a 30-m pixel size to coincide with the vegetation indices from the Landsat OLI images [16,20,38,39]. Recently, several studies used the ESTARFM approach to generate synthetic imagery to track the phenological changes and found that fused Landsat-MODIS data could deduce the land surface phenology at sub-field scales more accurately than using the MODIS data alone [16,20,39,40].
For the Landsat ETM+ images that had striping lines due to SLC failure, we only applied the fusion method to process the Landsat OLI images. The ESTARFM input format protocol dictated that all MODIS images were re-projected into the same projection as the Landsat images and resampled to a 30-m resolution by nearest-neighbor interpolation. We assumed that the clouds would not coincide completely in the different images. Therefore, we used more pairs of the observed Landsat/MODIS vegetation indices with available pixels for the same day to predict the vegetation indices on the same scale as the Landsat images for the other MODIS observation dates. We then compared the predicted vegetation indices at the Landsat-scale with the Landsat OLI indices for the same date and only replaced the pixels that emerged from clouds and cloud shadows with fusion data ( Figure 2). In this way, we tried to reduce the influences of the clouds and increased the available Landsat-scale pixels during the rice growing season.  For the Landsat ETM+ images that had striping lines due to SLC failure, we only applied the fusion method to process the Landsat OLI images. The ESTARFM input format protocol dictated that all MODIS images were re-projected into the same projection as the Landsat images and resampled to a 30-m resolution by nearest-neighbor interpolation. We assumed that the clouds would not coincide completely in the different images. Therefore, we used more pairs of the observed Landsat/MODIS vegetation indices with available pixels for the same day to predict the vegetation indices on the same scale as the Landsat images for the other MODIS observation dates. We then compared the predicted vegetation indices at the Landsat-scale with the Landsat OLI indices for the same date and only replaced the pixels that emerged from clouds and cloud shadows with fusion data (Figure 2). In this way, we tried to reduce the influences of the clouds and increased the available Landsat-scale pixels during the rice growing season.

Crop Calendar
The croplands in the study area were mainly planted with rice, and oilseed rape cultivation uses a different cropping calendar from other crops ( Figure 4). Rice seeds are generally germinated in small seedbed nurseries; rice seedlings are then transplanted to the paddy fields in one month.
Flooding is one of the key practices of rice paddy agriculture and often occurs a few weeks before rice transplantation [41]. The rice paddies are often planted with single-cropping rice and double-cropping rice. Double-cropping rice systems consist of early rice and late rice; the early rice is sown in late March, transplanted in late April, and harvested in late July. This is followed by the transplantation of the late rice (sowing in small seedbed nurseries in late June), which is harvested in late October. Compared with the double-cropping rice, the single-cropping rice plants were transplanted once per year in late May and harvested in early October ( Figure 4). Oilseed rape plants are sown in the autumn of the first year and then harvested in the spring of the next year. In addition to the crops mentioned above, there were many lotus ponds in the region, which began to germinate in early April and withered in late September.
a different cropping calendar from other crops ( Figure 4). Rice seeds are generally germinated in small seedbed nurseries; rice seedlings are then transplanted to the paddy fields in one month. Flooding is one of the key practices of rice paddy agriculture and often occurs a few weeks before rice transplantation [41]. The rice paddies are often planted with single-cropping rice and doublecropping rice. Double-cropping rice systems consist of early rice and late rice; the early rice is sown in late March, transplanted in late April, and harvested in late July. This is followed by the transplantation of the late rice (sowing in small seedbed nurseries in late June), which is harvested in late October. Compared with the double-cropping rice, the single-cropping rice plants were transplanted once per year in late May and harvested in early October (Figure 4). Rapeseed oil plants are sown in the autumn of the first year and then harvested in the spring of the next year. In addition to the crops mentioned above, there were many lotus ponds in the region, which began to germinate in early April and withered in late September.
Our study area is located in a subtropical region, where the forests are dominated by a mixture of evergreen and deciduous forests which (along with wetlands) have a longer growing season than rice. The plants in the natural wetlands begin to germinate in mid-March. When wetlands are in the leafing stage and have closed canopies in early or middle April, the early and single-cropping rice paddies are still not transplanted. Meanwhile, natural wetlands are often inundated in the summer when the rainy season occurs. The asynchrony of rice paddy and wetland growth phases makes it possible to distinguish between them in satellite images that cover multiple times. . Phenology of land-cover types found in our study area. Rice: 1-sowing, 2-seedling, 3flooding/transplanting, 4-tillering to stem elongation, 5-panicle initiation to heading, 6filled/milky, 7-maturity and harvest. Oilseed rape: 1-sowing, 2-seedling, 3-fifth true leaf, 4transplanting, 5-flower bud, 6-bolting, 7-flowering, 8-mature. Deciduous forest: 1-dormancy; 2-budburst, 3-leaf unfolding, 4-total leaf expansion, 5-leaf coloration, and 6-defoliation. Natural wetland: 1-dormancy, 2-germination, 3-trophophase, 4-flowering, 5-withering period. Lotus root: 1-dormancy, 2-germination, 3-trophophase, 4-flowering, 5-fruiting, 6withering period.

Maps of Rice Paddy Fields and Multi-cropping Areas
There are several land-use types that are not croplands, which include built-up and bare land, evergreen vegetation, natural deciduous vegetation, water bodies, and natural wetlands. To minimize the potential impacts of these regions on our rice paddy maps, we initially designed a procedure to mask these factors (see Figure 2). The built-up and bare land have a LSWI < 0 throughout the year [42]. However, the built-up area is often interspersed with vegetation (e.g., forest or grassland), which results in LSWI values above 0 in the mixed pixels. We used the "trial and error" method and found that the conditions with LSWI < 0.1 could extract the built-up area well. We calculated the number of observations of LSWI < 0.1 in one year and divided them by the total number of high quality observations in the pixels; i.e., the yearly frequency of observations with LSWI < 0.1. The pixels with frequencies of >50% were classified as built-up and/or bare land. Evergreen plants and trees have green leaves throughout the year and have a LSWI > 0 in all high quality records [43]. Phenology of land-cover types found in our study area. Rice: 1-sowing, 2-seedling, 3-flooding/transplanting, 4-tillering to stem elongation, 5-panicle initiation to heading, 6-filled/milky, 7-maturity and harvest. Oilseed rape: 1-sowing, 2-seedling, 3-fifth true leaf, 4-transplanting, 5-flower bud, 6-bolting, 7-flowering, 8-mature. Deciduous forest: 1-dormancy; 2-budburst, 3-leaf unfolding, 4-total leaf expansion, 5-leaf coloration, and 6-defoliation. Natural wetland: 1-dormancy, 2-germination, 3-trophophase, 4-flowering, 5-withering period. Lotus root: 1-dormancy, 2-germination, 3-trophophase, 4-flowering, 5-fruiting, 6-withering period.
Our study area is located in a subtropical region, where the forests are dominated by a mixture of evergreen and deciduous forests which (along with wetlands) have a longer growing season than rice. The plants in the natural wetlands begin to germinate in mid-March. When wetlands are in the leafing stage and have closed canopies in early or middle April, the early and single-cropping rice paddies are still not transplanted. Meanwhile, natural wetlands are often inundated in the summer when the rainy season occurs. The asynchrony of rice paddy and wetland growth phases makes it possible to distinguish between them in satellite images that cover multiple times.

Maps of Rice Paddy Fields and Multi-Cropping Areas
There are several land-use types that are not croplands, which include built-up and bare land, evergreen vegetation, natural deciduous vegetation, water bodies, and natural wetlands. To minimize the potential impacts of these regions on our rice paddy maps, we initially designed a procedure to mask these factors (see Figure 2). The built-up and bare land have a LSWI < 0 throughout the year [42]. However, the built-up area is often interspersed with vegetation (e.g., forest or grassland), which results in LSWI values above 0 in the mixed pixels. We used the "trial and error" method and found that the conditions with LSWI < 0.1 could extract the built-up area well. We calculated the number of observations of LSWI < 0.1 in one year and divided them by the total number of high quality observations in the pixels; i.e., the yearly frequency of observations with LSWI < 0.1. The pixels with frequencies of >50% were classified as built-up and/or bare land. Evergreen plants and trees have green leaves throughout the year and have a LSWI > 0 in all high quality records [43]. Our study area is located in a subtropical region, and the forests are dominated by a mixture of evergreen and deciduous forests. Based on random sampling, we ascertained the characteristics of the yearly curves of the forest vegetation indices (Figures 5a and 6a), and then designed a rule to exclude the forest land-cover type as follows. First, we counted the number of observations in which a pixel had LSWI > 0.1 from the set of all high quality observations, and then divided it by the total number of high Remote Sens. 2020, 12, 1022 8 of 18 quality records. The pixels with a ratio of > 95% were subsequently identified as forest. The water pool usually displayed low NDVI and EVI values and high LSWI values (Figures 5b and 6b). The water pools were extracted by capturing the pixels with NDVI < 0.1 and LSWI > NDVI (or EVI) for each Landsat image [10]. The water pools were labeled as persistent if these pixels have a ratio of > 80% in all high quality observations [24].
Of all the dominant land-cover types, rice paddies, lotus root ponds, and natural wetlands are sensitive to flooding/inundation events in our study area. Several studies suggested that LSWI ≥ NDVI or LSWI ≥ EVI are representative of flooding events [10,44]. These phenomena are also displayed in Figures 5b-f and 6b-f. Pixels that were identified as flooded throughout a year might be pure water or mixtures of water and vegetation in which the water dominates the pixel (Figures 5b and 6b). Our study area is located in a subtropical region, and the forests are dominated by a mixture of evergreen and deciduous forests. Based on random sampling, we ascertained the characteristics of the yearly curves of the forest vegetation indices (Figures 5a and 6a), and then designed a rule to exclude the forestland-cover type as follows. First, we counted the number of observations in which a pixel had LSWI > 0.1 from the set of all high quality observations, and then divided it by the total number of high quality records. The pixels with a ratio of > 95% were subsequently identified as forest. The water pool usually displayed low NDVI and EVI values and high LSWI values (Figures  5b and 6b). The water pools were extracted by capturing the pixels with NDVI < 0.1 and LSWI > NDVI (or EVI) for each Landsat image [10]. The water pools were labeled as persistent if these pixels have a ratio of > 80% in all high quality observations [24]. Of all the dominant land-cover types, rice paddies, lotus root ponds, and natural wetlands are sensitive to flooding/inundation events in our study area. Several studies suggested that LSWI ≥ NDVI or LSWI ≥ EVI are representative of flooding events [10,44]. These phenomena are also displayed in Figure 5b-f and Figure 6b-f. Pixels that were identified as flooded throughout a year might be pure water or mixtures of water and vegetation in which the water dominates the pixel (Figures 5b and 6b).   Figure 6. As in Figure 5, but for the seasonal dynamics of NDVI, EVI, and LSWI of the major landcover types from fused Landsat-MODIS data.
After excluding natural wetland and lotus ponds from the flooding or inundated pixels, most of the remaining research area is constituted by rice paddies. Using the remaining pixels combined with the phenology of different trice paddy types, we designed a set of rules to differentiate between single-and double-cropping rice paddies. Single-cropping rice was identified by pixels with LSWI > NDVI (or EVI) from mid-May to mid-June and NDVI > 0.8 in August (Figures 5e and 6e). In Figure 6. As in Figure 5, but for the seasonal dynamics of NDVI, EVI, and LSWI of the major land-cover types from fused Landsat-MODIS data.
Pixels that were identified as seasonally flooded included natural wetlands, lotus ponds, and rice paddies. Based on the remaining research area, and after excluding the built-up and barren land, forest, and pure water from the research region, we extracted the seasonally flooded pixels by identifying the pixels with LSWI > NDVI (or EVI) for each Landsat image if these pixels were obtained in all high quality observations from March to October. The natural wetlands and lotus ponds displaced the water body phenomena (Figures 5c-d and 6c-d) after September, in which a single grain of rice appears filled/milky, and later rice shows panicle initiation phenology (Figure 4) without flooding. We counted the number of records in which a pixel had LSWI > EVI out of all the high quality records occurring after September. The pixels that showed a ratio (frequency) value of >50% were identified as natural wetlands and lotus fields.
After excluding natural wetland and lotus ponds from the flooding or inundated pixels, most of the remaining research area is constituted by rice paddies. Using the remaining pixels combined with the phenology of different trice paddy types, we designed a set of rules to differentiate between singleand double-cropping rice paddies. Single-cropping rice was identified by pixels with LSWI > NDVI (or EVI) from mid-May to mid-June and NDVI > 0.8 in August (Figures 5e and 6e). In comparison, double-cropping rice was identified by pixels with LSWI > NDVI (or EVI) from mid-April to mid-May and NDVI > 0.8 in September (Figures 5f and 6f).

Validation and Comparison
The accuracy of the ESTARFM-fusion images was evaluated by comparing them with the corresponding real Landsat images in the 3000 random samples, where the observation had good quality. The correlation coefficient (r), mean absolute error (MAE), and root mean square error (RMSE) were estimated as evaluation indicators to assess the accuracy of the fused NDVI/EVI/LSWI values, according to the following functions: In which N indicates the number of pixels; P and O denote the value of each pixel in the fused images and the corresponding real image, respectively; and P and O represent the average values of P and O, respectively.
After mapping the rice paddy areas and multi-cropping regions, we used the confusion matrix accuracy assessment to validate the product. The validation points, including the 200 field survey sites and the 2040 sites extracted from the high-resolution images in Google Earth, were divided into three categories (single-cropping rice fields, double-cropping rice fields, and non-rice paddy fields). Then, these validation points were applied to evaluate the accuracy of the resultant rice paddy maps and tillage characteristics maps. We also calculated the user accuracy and producer accuracy of the map categories, overall accuracy and kappa coefficient [24,45].
Finally, we compared the map of rice paddies derived from the fused Landsat-MODIS data in our study with the rice paddy land in the NLCD data (hereafter referred to as ILMP-Rice and NLCD-Rice) for 2015 in Nanchang County. The purpose of the comparison was to verify whether multi-source satellite imagery and fusion algorithm could effectively map rice paddy planting areas compared to the existing data of rice paddy areas in southern China.

Landsat-MODIS Fusion Results
We used the ESTARFM approach to fuse the Landsat and MODIS images. Figure 7 shows the comparison of the fused Landsat-MODIS data and the actual Landsat data for the three indices (NDVI, EVI, and LSWI). The scatter plots display high correlations and small biases for the three indices between the fused data and actual data. The r values were above 0.93 for the NDVI and EVI values and 0.84 for the LSWI value; the MAE values for the three indices were within 0.07; and the RMSE values for the three indices were not higher than 0.10. Additionally, in these scatter plots, we also found that the data for the three indices (NDVI, EVI, and LSWI) were close to the 1:1 line. Figure 8 shows that the influences of clouds and cloud shadows were almost eliminated by utilizing the fused data. These results indicate that the ESTARFM approach can synthesize fine resolution images with high accuracy. Therefore, the fusion results can add more available data to improve the ability to discern between the targeted land-cover types.

Landsat-MODIS Fusion Results
We used the ESTARFM approach to fuse the Landsat and MODIS images. Figure 7 shows the comparison of the fused Landsat-MODIS data and the actual Landsat data for the three indices (NDVI, EVI, and LSWI). The scatter plots display high correlations and small biases for the three indices between the fused data and actual data. The r values were above 0.93 for the NDVI and EVI values and 0.84 for the LSWI value; the MAE values for the three indices were within 0.07; and the RMSE values for the three indices were not higher than 0.10. Additionally, in these scatter plots, we also found that the data for the three indices (NDVI, EVI, and LSWI) were close to the 1:1 line. Figure  8 shows that the influences of clouds and cloud shadows were almost eliminated by utilizing the fused data. These results indicate that the ESTARFM approach can synthesize fine resolution images with high accuracy. Therefore, the fusion results can add more available data to improve the ability to discern between the targeted land-cover types.

Landsat-MODIS Fusion Results
We used the ESTARFM approach to fuse the Landsat and MODIS images. Figure 7 shows the comparison of the fused Landsat-MODIS data and the actual Landsat data for the three indices (NDVI, EVI, and LSWI). The scatter plots display high correlations and small biases for the three indices between the fused data and actual data. The r values were above 0.93 for the NDVI and EVI values and 0.84 for the LSWI value; the MAE values for the three indices were within 0.07; and the RMSE values for the three indices were not higher than 0.10. Additionally, in these scatter plots, we also found that the data for the three indices (NDVI, EVI, and LSWI) were close to the 1:1 line. Figure  8 shows that the influences of clouds and cloud shadows were almost eliminated by utilizing the fused data. These results indicate that the ESTARFM approach can synthesize fine resolution images with high accuracy. Therefore, the fusion results can add more available data to improve the ability to discern between the targeted land-cover types.

Maps of Rice Paddy Planting Area and Cropping Intensity in 2015
We generated a map of rice paddy planting areas from the fused Landsat-MODIS data with a 30-m spatial resolution for Nanchang County in 2015 (Figure 9). The results showed that the rice paddy planting area was 540.59 km 2 . The rice paddy regions were mainly distributed on flat terrain with relative large patches in the northern part of Nanchang County, while a few were distributed in the southern part of the region where some topographic relief exists. Using these rice paddy planting areas, we extracted the cropping intensity. The areas of single-and double-cropping rice were 266.87 and 273.72 km 2 , respectively, indicating that the planting areas of single-cropping rice were close to those of double-cropping rice in 2015 and that the multiple cropping index was approximately 150%.
Spatially, the planting areas of the double-cropping rice were mainly located in the northern part of Nanchang County, while the single-cropping rice areas were primarily distributed in the middle of Nanchang County. paddy planting area was 540.59 km 2 . The rice paddy regions were mainly distributed on flat terrain with relative large patches in the northern part of Nanchang County, while a few were distributed in the southern part of the region where some topographic relief exists. Using these rice paddy planting areas, we extracted the cropping intensity. The areas of single-and double-cropping rice were 266.87 and 273.72 km 2 , respectively, indicating that the planting areas of single-cropping rice were close to those of double-cropping rice in 2015 and that the multiple cropping index was approximately 150%. Spatially, the planting areas of the double-cropping rice were mainly located in the northern part of Nanchang County, while the single-cropping rice areas were primarily distributed in the middle of Nanchang County.

Accuracy Assessment of Rice Paddy Planting Areas and Cropping Intensity in 2015
The ILMP-Rice map used here was compared with validation points from high-resolution images via Google Earth and the in-situ ground-truth data (Tables 1 and 2). The assessment results showed that the kappa coefficient and overall accuracy were 0.85% and 93.66%, respectively, and that the user and producer accuracy were 91.86% and 99.4% for rice paddy regions, respectively (Table 1). This result indicated a very high agreement between the ILMP-Rice map and ground-based data of rice paddy areas. We further found that the ILMP-Rice multi-cropping map also has high kappa coefficients and accuracies. The kappa coefficient and overall accuracy were 0.89% and 92.95%, respectively. The user and producer accuracies were 90.30% and 99.09%, respectively, for singlecropping rice and 91.33% and 97.77%, respectively, for double-cropping rice ( Table 2). These results

Accuracy Assessment of Rice Paddy Planting Areas and Cropping Intensity in 2015
The ILMP-Rice map used here was compared with validation points from high-resolution images via Google Earth and the in-situ ground-truth data (Tables 1 and 2). The assessment results showed that the kappa coefficient and overall accuracy were 0.85% and 93.66%, respectively, and that the user and producer accuracy were 91.86% and 99.4% for rice paddy regions, respectively (Table 1). This result indicated a very high agreement between the ILMP-Rice map and ground-based data of rice paddy areas. We further found that the ILMP-Rice multi-cropping map also has high kappa coefficients and accuracies. The kappa coefficient and overall accuracy were 0.89% and 92.95%, respectively. The user and producer accuracies were 90.30% and 99.09%, respectively, for single-cropping rice and 91.33% and 97.77%, respectively, for double-cropping rice ( Table 2). These results suggested that our mapping results were reasonable and could meet the requirements of user and producer accuracy in the complex landscapes of rice paddies.

Comparison between ILMP-Rice Map and NLCD-Rice Data in 2015
There are some notable differences between the ILMP-Rice map and the NLCD-Rice data. A comparison of the ILMP-Rice map and the NLCD-Rice map revealed that the NLCD-Rice map had a remarkable overestimation of the rice paddies compared to ILMP-Rice for 2015 in Nanchang County ( Figure 10). The analysis of a zoomed-in image of one region showed that the NLCD-Rice misclassified roads, village buildings, and ditches as rice paddies. As NLCD-Rice used images without or with less cloud contamination for interpretation, these images were often obtained in the autumn and winter seasons when the rice was harvested and the ditches were dry. The rice paddy fields, small village buildings, and ditches were difficult to separate from each other spectrally (Figure 10b-d). We utilized the same validation samples to analyze the NLCD-Rice data for 2015 and found that both the accuracy and kappa coefficient were lower than those of the ILMP-Rice map (Tables 1 and 3). Particularly, the relatively low producer accuracy and high misclassification of the NLCD-Rice data suggested that the abundant small village buildings, pools, roads, ditches, and wetlands were misclassified into rice paddy fields. The NLCD-Rice data in some areas were derived from the images during the peak vegetation growing season; however, there were two types of tillage (single-and double-cropping rice) in the study area. For example, the image used for the NLCD-Rice was obtained on July 12 of 2015 ( Figure 10b); at this time, the early rice in the double-cropping rice systems was harvested, but the later rice was not yet planted, so the area of the rice paddy field was bare, and the NLCD-Rice may have misclassified it as upland crops or other land use types.

Discussion
This study developed an improved rice paddy mapping approach by combining a spatiotemporal fusion algorithm with a phenology-based algorithm to map rice paddy fields in southern China. First, we employed the ESTARFM approach to fuse Landsat and MODIS images. The availability of higher-quality data is essential for mapping rice paddy fields and cropping intensity based on optical images at the regional scale. However, previous studies noted that Landsat ETM+ and OLI often overlook some rice paddy regions, mainly due to restrictions of data availability

Discussion
This study developed an improved rice paddy mapping approach by combining a spatiotemporal fusion algorithm with a phenology-based algorithm to map rice paddy fields in southern China. First, we employed the ESTARFM approach to fuse Landsat and MODIS images. The availability of higher-quality data is essential for mapping rice paddy fields and cropping intensity based on optical images at the regional scale. However, previous studies noted that Landsat ETM+ and OLI often overlook some rice paddy regions, mainly due to restrictions of data availability and quality [13,46]. This study tested the fusion method in southern China and confirmed that the fusion of MODIS and Landsat images was a feasible strategy for adding more available data (Figures 7 and 8) to improve the discernment capability of phenological information for rice paddy classification [26]. Second, pixel-and phenology-based algorithms were applied for mapping rice paddies. This algorithm has mostly been used in the cold temperate zones of northeast Asia, where there is one single-cropping season [3,14,24]. Our study extended the algorithm based on fused time series images into the subtropical region, where there are small cropland field sizes and multiple cropping systems employed. The accuracy assessment of the resultant map showed high accuracy and a high kappa coefficient (Table 1). Meanwhile, the comparison with the existing NLCD-Rice data revealed that the ILMP-Rice map displayed more details about the spatial distribution and lower misclassification in the fragmented areas with diverse land-cover types ( Figure 10 and Table 3). A similar algorithm has also been proven to perform better than the conventional phenology-based algorithm based on a single satellite sensor in areas with a single-cropping system [14]. These results suggest that the fusion-and phenology-based algorithms have great potential for mapping rice paddies in the complex landscapes of southern China.
Compared with previous algorithms used in areas with a single-cropping system [14], our approach could also be applied to improve the mapping of rice paddy cropping intensity in areas with complex cropping systems in southern China. Masking other land-cover types has proven essential for mapping rice paddies and cropping intensity [13], and improved mapping of rice paddy planting areas at high resolutions can improve mapping of cropping intensity [26]. Our study found that the multiple cropping indices were approximately 150%. Jiang et al. (2018Jiang et al. ( , 2019 noted that the planting area of double-cropping rice and multiple cropping indices considerably decreased in southern China, especially in the Poyang Lake plain and Hubei and Hunan plain [6,47]. They employed NLCD data to generate rice paddy masks. However, abundant small village buildings, pools, ditches, and wetlands were misclassified into rice paddies in the NLCD data ( Figure 10), which can lead to biased estimates of the sown areas, and hence the derived cropping intensity. In addition, they mapped the rice paddy planting area and cropping intensity using a phenology-based algorithm based on Landsat data alone. However, southern China has insufficient Landsat observations for annual analyses of rice paddy cropping intensity [48]. Our study has improved rice paddy mapping methods, as the accuracy assessment of the rice paddy cropping intensity map showed high accuracy and a high kappa coefficient at the county scale (Table 2). In recent years, the Google Earth Engine platform has provided unprecedented opportunities for geospatial analysis at the regional and global scales [49] and has been used to track changes in rice paddy planting area [32]. Thus, combining the fusionand phenology-based algorithm with the Google Earth Engine platform in future studies would assist a retrospective analysis of rice paddy cropping intensity change at the regional scale, such as in southern China.
The fusion algorithm could reasonably deduce the crop phenology metrics at the field scale by combining the VI data from multi-source satellite optical imagery. Recently, the free and open access Sentinel-1 and Sentinel-2 satellites have provided imagery with higher resolution than that of Landsat, which could be combined with the Landsat time series images for mapping rice paddies. The combined use of optical and radar imagery would further improve rice paddy mapping in the cloudy and rainy regions [50,51]. However, recent optical-SAR-synergistic rice mapping algorithms were generally applied to small study sites. Therefore, a retrospective analysis of rice paddy cropping intensity since the new millennium, at a large scale, with the fusion of multi-source optical time-series data, is still an effective strategy.

Conclusions
We presented a rice paddy mapping approach for the complex landscape of southern China by combining a spatiotemporal fusion algorithm with a phenology-based algorithm. This approach was tested for 2015 in Nanchang County, which is located on the Poyang Lake plain of southern China. The validation tests indicated a high accuracy of rice paddy fields and the multi-cropping map. The overall accuracy and kappa coefficient of the rice paddy planting area were 93.66% and 0.85, respectively; and the overall accuracy and kappa coefficient of the multi-cropping systems were 92.95% and 0.89, respectively. The comparison of our map with NLCD-Rice map yielded high levels of consistency and indicated that our products displayed more detailed information about the spatial distribution of rice paddies. With the improvement of cloud contamination detection and the fusion method, our improved mapping approach has great potential for rice paddy mapping and retrospective analysis of rice paddy cropping intensity in fragmented areas with diverse land-cover types found in subtropical regions.
Author Contributions: M.D. and L.L. conceived and designed the experiments; M.D. and L.L. wrote the paper; Q.G. and H.Z. did the satellite imagery pre-processing and data analysis; C.L. and L.Z. revised the paper and contributed to explanation of the results and the discussion. All authors have read and agreed to the published version of the manuscript.