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

Abstract

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 fields and multi-cropping at fine resolution are still lacking. Landsat observation is the primary source of remote sensing data that has continuously mapped regional rice paddy fields 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 fields 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 flooding 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 coefficient of the rice paddy planting areas were 93.66% and 0.85, respectively. The overall accuracy and kappa coefficient 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 fields and using the multi-cropping index in complex landscapes in southern China.

1. 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 global population rapidly increasing, there has been growing pressure on agricultural systems due to the rising demand for food, especially rice [2]. Within the context of climate change and rapid economic development and urbanization, however, the increasing area of rice paddy fields has been threatened. This has led to an intensified rice paddy cropping cycle in some areas with limited rice paddy lands [4,5], while it has caused a decrease in the rice paddy planting areas of some developing countries [6]. In addition, the continuous modification and intensive management of rice paddy agriculture has also substantially influenced the biogeochemical and hydrologic cycles, climate change, and biodiversity [7]. Thus, achieving an increase in rice paddy yield to ensure food security and minimize the negative environmental influence of rice paddy agriculture is still one of the current challenges [8,9]. Timely and accurate information about where, when, and how much paddy rice is planted is a key step in food security assessment, agricultural resource management, and ecosystem sustainability.

Remote sensing with satellite imagery is considered an effective and valuable tool for rice paddy mapping [10,11], typically superior to conventional ground-based surveys [6,12]. Various efforts to map global and regional crop areas have employed coarse resolution MODIS satellite images, owing to their high revisit frequencies (i.e., daily) [5,13]. However, some studies have noted that MODIS and AVHRR pixels generally contain different crop types in the areas with patchy and fragmented fields, particularly in countries such as China, and thus, small rice paddy field patches are ignored [13,14,15]. Comparatively, long-term Landsat satellite images with 30-m spatial resolution provide unmatched resources to monitor and map rice paddy fields in fragmented landscapes. However, the accuracies of recent methods for mapping rice paddies largely depend on the availability of Landsat time series images [13], which are primarily limited by the long revisit interval (16 days) and by weather conditions, particularly those that occur during the rice growing season, which is often contaminated by frequent cloud coverage [12,13]. Overall, it is difficult for MODIS satellite images alone or Landsat satellite images alone to meet the requirements for mapping complex and fragmented rice paddy areas [14]. Data fusion algorithms provide a method for bridging between medium-resolution Landsat images and high temporal resolution MODIS images [16,17,18]. A spatial and temporal adaptive reflectance fusion model (STARFM) and its updated version—the enhanced STARFM (ESTARFM) approach—have been widely used to synthesize high spatio-temporal resolution fused Landsat–MODIS imagery [19,20]. One important advantage of ESTARFM is the application of a conversion coefficient to improve the accuracy of the reconstruction in areas with heterogeneous landscapes [20,21].

In recent decades, a variety of rice paddy mapping approaches have been developed (e.g., the reflectance data and image statistic-based approach and the phenology-based approach) [13]. The phenology-based approach can monitor different crop growth properties and detect the unique signals of the rice paddies through their differences in their respective crop schedule or management [22]. Recently, phenology-based approaches have been increasingly applied to map rice paddy planting areas [6,10,23]. Rice paddy fields are a mixture of surface water and planted rice during the flooding and open-canopy periods (i.e., after rice transplanting) and exert unique biophysical characteristics, which could be captured by the relationship between the normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and land surface water index (LSWI) [10,24]. Based on the NDVI, EVI, and LSWI, a transplanting-based algorithm can effectively detect phenological information about the rice paddy and recognize the key phenological periods from time series imagery [10]. This straightforward algorithm has been widely used and become one of the most representative phenology-based approaches for rice paddy mapping [13,25]. However, the periods of flooding and transplanting are very short, and many regions are in the rainy season during this period, particularly in southern China; thus, transplanting-based algorithms require adequate availability of high-quality observations to achieve high accuracy [13,25]. Additionally, a recent fusion and phenology-based rice paddy mapping algorithm has been used in northeast Asia, where there is one single cropping season [14]. However, previous studies noted that the extended application of the mapping algorithm to southern China required additional case studies and algorithmic improvement [24].

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.

2. Study Area and Data Sources

2.1. 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² and belongs to the subtropical climate type. Major annual crops in the study area include single- and double-cropping rice, oilseed rape, 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].

2.2. Data Sources

2.2.1. Remote Sensing Images

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/).

2.2.2. 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.).

2.2.3. 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.

3. Methodology

3.1. 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 assessment (Section 3.4). To study the methodology, we selected one county-level administrative region (path/row of Landsat imagery: 121/40) on the Poyang Lake plain in China as the study area.

3.2. Preprocessing of Remote Sensing Image

3.2.1. 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 (${\rho }_{blue}$), red (${\rho }_{red}$), NIR (${\rho }_{NIR}$), and SWIR (${\rho }_{SWIR}$, 1570–1650 nm) bands. These indices have proven to be helpful in studies of vegetation canopies, water content, biomass, and phenology.

$$NDVI=\frac{{\rho }_{NIR}-{\rho }_{red}}{{\rho }_{NIR }+{\rho }_{red}},$$

(1)

$$EVI=2.5\times \frac{{\rho }_{NIR}-{\rho }_{red}}{{\rho }_{NIR }+6\times {\rho }_{red}-7.5\times {\rho }_{blue}+1},$$

(2)

$$LSWI=\frac{{\rho }_{NIR}-{\rho }_{SWIR}}{{\rho }_{NIR+}{\rho }_{SWIR}},$$

(3)

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].

3.2.2. 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.

3.3. Pixel- and Phenology-based Paddy Rice Mapping Algorithms

3.3.1. 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.

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.

3.3.2. 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 (Figure 5a and Figure 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 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 (Figure 5b and Figure 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 (Figure 5b and Figure 6b).

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 (Figure 5c–d and Figure 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 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 (Figure 5e and Figure 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 (Figure 5f and Figure 6f).

3.4. 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:

$$r=\frac{{\sum }_{i=1}(P_i-\overline{P}) \left(O_i-\overline{O}\right)}{\sqrt{{\sum }_{i=1}{\left(P_i-\overline{P}\right)}^2{\left(O_i-\overline{O}\right)}^2}},$$

(4)

$$MAE=\frac{1}{N}{\displaystyle \sum }_{i=1}^n│P_i-O_i│,$$

(5)

$$RSME=\sqrt{\frac{1}{N}{{\displaystyle \sum }}_{i=1}^n{\left(P_i-O_i\right)}^2}$$

(6)

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.

4. Results

4.1. 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.

4.2. 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². 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², 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.

4.3. 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 (

Table 1. Confusion matrix between ILMP-Rice in Nanchang County and the validation points from the high-resolution images via Google Earth and in-situ data.

		Ground Truth Pixels		Classified Pixels	User Accuracy
		Rice Paddies	Non-Rice
	Rice paddies	1500	133	1633	91.86%
	Non-rice	9	598	607	98.52%
Ground truth pixels		1509	731	2240	OA = 93.66%
Producer accuracy		99.4%	81.81%		Kappa = 0.85

and

Table 2. Confusion matrix between ILMP-Rice in Nanchang County and validation points from high-resolution images via Google Earth and in-situ data.

		Ground Truth Pixels			Classified Pixels	User Accuracy
		Single-Cropping Rice	Double-Cropping Rice	Non-Rice
Classified results	Single-cropping rice	652	14	56	722	90.3%
	Double-cropping rice	2	832	77	911	91.33%
	Non-rice	4	5	598	607	98.52%
Ground truth pixels		658	851	731	2240	OA = 92.95%
Producer accuracy		99.09%	97.77%	81.81%		Kappa = 0.89

). 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.

4.4. 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 (Table 1 and

Table 3. Confusion matrix between NLCD-Rice in Nanchang County and validation points from high-resolution images via Google Earth and in-situ data.

		Ground Truth Pixels		Classified Pixels	User Accuracy
		Rice Paddies	No-Rice
	Rice paddies	1526	107	1633	93.45%
	Non-rice	155	452	607	74.46%
Ground truth pixels		1681	559		OA = 88.30%
Producer accuracy		90.78%	80.86%		Kappa = 0.70

). 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.

5. 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 (Figure 7 and Figure 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. (2018, 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 fusion- and 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.

6. 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.