Mapping paddy rice in Jiangsu Province, China, based on phenological parameters and a decision tree model

Abstract

Timely and accurate mapping of rice planting areas is crucial under China’s current cropping structure. This study proposes a new paddy rice mapping method by combining phenological parameters and a decision tree model. Six phenological parameters were developed to identify paddy rice areas based on the analysis of the Moderate Resolution Imaging Spectroradiometer (MODIS) Enhanced Vegetation Index (EVI) time series and the Land Surface Water Index (LSWI) time series. The six phenological parameters considered the performance of different land cover types during specific phenological phases (EVI₁ and EVI₂), one-half of or the entire rice growing cycle (LSWI₁ and LSWI₂), and the shape of the LSWI time series (Kurtosis_LSWI and Skewness_LSWI). A hierarchical decision tree model was designed to classify paddy rice areas according to the potential separability of different land cover types in paired phenological parameter spaces. Results showed that the decision tree model was more sensitive to LSWI₁, LSWI₂, and Skewness_LSWI than the other phenological parameters. A paddy rice map of Jiangsu Province for 2015 was generated with an optimal threshold set of (0.4, 0.42, 9, 19, 1.5, –1.7, 0.0) with a total accuracy of 93.9%. The MODIS-derived paddy rice map generally agreed with the paddy land fraction map from the National Land Cover Dataset project, but there were regional discrepancies because of their different definitions of land use and the inability of MODIS to map paddy rice at a fragmental level. The MODIS-derived paddy rice map showed high correlation (R²=0.85) with county-level agricultural statistics. The results of this study indicate that the phenological parameter-based paddy rice mapping algorithm could be applied at larger spatial scales.

Keywords

phenological parameter / paddy rice / MODIS / EVI / LSWI / decision tree

Introduction

Rice is among the three major cereal crops in China, and paddy rice is widely planted in southern China in many types of crop rotations, such as rapeseed-rice, winter wheat-rice, single rice, double rice, and triple rice (^{Frolking et al., 2002}). Paddy rice consumes much more water during its growing cycle than other crops (^{Bouman, 2009}; ^{Dvorak, 2012}). Under the current changing climate, the spatial-temporal distribution of water resources will significantly change, and extreme climatic disasters may frequently occur (^{Matthews and Caldeira, 2007}; Robock et al., 2008). These factors may have a significant impact on paddy rice production (^{Pongratz et al., 2012}). Furthermore, because of economic demands, many farmers have changed rice cultivation to grow cash crops or reduced their rice cultivation frequency. Currently, China is working to adjust its cropping structure, and the approach to adjusting rice planting structure is highly dependent on the temporal intensity and geographical extent of rice cropping. Furthermore, an accurate distribution map of paddy rice is important to understand rice production and develop agricultural policies. Understanding rice distribution is also critical for food trade, rice growth monitoring, rice yield prediction, and global gaseous emission research (^{Labus et al., 2002}; ^{Zou et al., 2007}; ^{Wardlow and Egbert, 2008}; ^{Sakamoto et al., 2009}; ^{Qiu et al., 2015}).

Phenology-based approaches using remote-sensing times-series data have been increasingly used for mapping paddy rice (^{Xiao et al., 2005}, ²⁰⁰⁶; ^{Peng et al., 2011}; ^{Qin et al., 2015}; ^{Dong and Xiao, 2016}). The algorithms for identifying paddy rice fields are mainly based on 1) vegetation indices that are sensitive to the dynamics of the canopy and surface water content, such as the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI) and Land Surface Water Index (LSWI) (^{Zhang et al. 2015}); and 2) distinct physical features that the paddy rice LSWI during flooding and transplanting is higher than its NDVI or EVI (^{Dong and Xiao, 2016}). However, the period of image composition for most vegetation index products is longer than eight days (e.g., Moderate Resolution Imaging Spectroradiometer [MODIS] EVI/NDVI, SPOT NDVI, and AVHRR NDVI), while the periods for flooding and transplanting are quite short. Moreover, during this period, most regions in southern China are in the rainy season, and the LSWI for paddy rice is always reduced because of rainfall and clouds. This greatly restricts the applicability of this method (^{Dong and Xiao, 2016}).

In southern China, many crops grow during the rice growing season, resulting in a much more complex cropping structure. Crops growing during the corresponding period have similar phenological dynamics and vegetation cover change patterns. Hence, it is difficult to directly distinguish rice from other crops using only one phenological feature. Combining different phenological parameters developed from key phenological periods can better reflect a crop’s biophysical characteristics. Many previous studies have showed that phenological parameters will greatly improve the accuracy of land cover classification (^{Loveland and Belward, 1997}; ^{Bradley and Mustard, 2008}). Recently, phenological parameters have been used to estimate net primary production (^{Kimball et al., 2004}), map crop types (^{Wardlow and Egbert, 2008}; Pringle et al., 2012; ^{Liu et al., 2014b}), map crop proportions (^{Pan et al., 2012}), and determine cropping intensities (^{Liu et al., 2012}). Phenological parameters integrate useful information and suppress non-critical signals. Using phenological parameters may also help discriminate between contemporaneous crops. This approach also reduces the volume of data, which enhances its applicability to different regions (^{Potgieter et al., 2010}, 2013).

The decision tree classifier is widely used in remote-sensing data analysis, feature selection, and classification (^{Friedl and Brodley, 1997}; ^{Pal and Mather, 2003}; ^{Tooke et al., 2009}). In addition to its rapid computation ability, it offers several other advantages when compared to other classifiers as it requires no statistical assumptions, can handle data measured at different scales, and can process non-linear relationships between features (^{Otukei and Blaschke, 2010}). Therefore, it has great flexibility in data analysis and is the favorite method for complex classification problems (^{Friedl and Brodley, 1997}). Many regional land use and land cover maps, and even several global land cover products, have been produced using the decision tree classifier (^{Loveland and Belward, 1997}; ^{Friedl et al., 2010}). However, the challenge in using the decision tree classification is the determination of thresholds.

This study attempted to map paddy rice areas in Jiangsu Province, China, based on multiple phenological parameters. MODIS EVI and LSWI time series data were used to retrieve six phenological parameters characterizing different land cover types. A hierarchical decision tree classification model was developed according to the analysis of the six phenological parameters. In addition, the sensitivity of the decision tree model to each parameter was also explored to obtain the optimal thresholds for mapping the paddy rice. The resultant paddy rice map could be used to support cropping system adjustment, cultivated land consolidation, rural land circulation, large-scale rice production, agricultural-related land use and land cover analysis, and water resources management.

Study area and data

　Study area

Jiangsu Province is on the eastern coast of China (Fig. 1). It is among the thirteen major grain-producing provinces in China. Its longitude ranges from 116°18′E to 121°57′E, and its latitude ranges from 30°45′N to 35°20′N. Its total land area is 102,600 km². Most regions in the study area are less than 50 m in altitude. According to the China Statistical Yearbook, the total cropland area in Jiangsu Province was 45,749 km² in 2015. The central part of the province has extensive plains, while low mountains and hills occur in the north and southwest. The region is situated in a transitional zone from a subtropical to a warm temperate climate. It is hot and rainy during summer and cold and dry during winter. The geographical and climatic conditions in Jiangsu Province are suitable for rice cultivation. A rotation of winter wheat-paddy rice occurs mostly in the northern study area, and a rotation of rapeseed-paddy rice occurs mostly in the central and southern areas. The main contemporary crops of rice are maize, soybean, and cotton.

　Data and preprocessing

　MODIS data

For this study, we obtained the 8-day composite MOD09A1 product with seven reflectance bands and a 500-m spatial resolution from June 2 to November 25, 2015. Tile h27v05 and h28v05 were mosaicked to cover the entire study area. MOD09A1 had already undergone atmospheric correction, and the downloaded data were reprojected to the Albers conical equal area projection using the MODIS Reprojection Tool (MRT). The EVI and the LSWI were calculated from the MODIS blue band (459–479 nm), red band (620–670 nm), near-infrared band (841–875 nm), and short-infrared band (1628–1652 nm) using Eqs. (1) and (2) as follows:

(1)\mathrm{EVI}=\mathrm{2.5}\times {\displaystyle \frac{{\rho }_{\mathrm{nir}}-{\rho }_{\mathrm{red}}}{{\rho }_{\mathrm{nir}}+\mathrm{6}\times {\rho }_{\mathrm{red}}-\mathrm{7.5}\times {\rho }_{\mathrm{blue}}+\mathrm{1}}},

(2)\mathrm{LSWI}={\displaystyle \frac{{\rho }_{\mathrm{nir}}-{\rho }_{\mathrm{swir}}}{{\rho }_{\mathrm{nir}}+{\rho }_{\mathrm{swir}}}}.

Unlike the NDVI, the EVI adds the blue band to the calculation to correct the effects of soil background and aerosol scattering, thus enhancing its sensitivity to the canopy biomass estimation of dense vegetation (^{Huete et al., 2002}). The LSWI is sensitive to vegetation water content and soil moisture. Although the MOD09A1 product chooses the optimal observation from each 8-day period to eliminate the effects of poor atmospheric conditions, other sources of data noise remain. Thus, the original EVI and the LSWI data were unfit to be directly used for land use/land cover analysis. The changing weight filter method developed by ^{Zhu et al. (2012)} was used to smooth the EVI and LSWI time series. The smoothed time series curves can more clearly reflect seasonal changes in different land cover types.

　Training and validation samples

To identify different land cover types and analyze their EVI/LSWI time series features, 11 scenes of Landsat-8 Operational Land Imager (OLI) (Fig. 1) were downloaded from the Open Spatial Data Sharing System (http://ids.ceode.ac.cn/) for the study period. The detailed date and cloud cover information are listed in Table 1. Additionally, two multispectral images from the fourth satellite of the China-Brazil Earth Resource Satellite series (CBERS-04, http://www.cresda.com/CN/Satellite/4308.shtml) were used as fine-resolution reference data to validate the rice map. The multispectral image of CBERS-04 included four spectral bands with wavelengths ranging from 450 to 890 nm. The spatial resolution was 20 m, and the swath width was 120 km. All Landsat-8 OLI images and CBERS-04 images were geometrically registered and converted to the same projection as that of the MODIS data.

　Training and validation samples

Based on the field investigation, the main land cover types during the second half of the year were water bodies, artificial surfaces, natural vegetation, mid-season rice, late-season rice, summer maize, autumn maize, soybean, cotton, lotus, and cash crops. Generally, mid-season rice is seeded during May and harvested during late October while late-season rice is seeded during June and harvested during late November. Fifteen samples of each land cover type (240 samples in total, as shown in Fig. 1) were generated as training samples from the MODIS data using visual interpretation of the EVI/LSWI time series, the Landsat-8 OLI imagery, and the field investigation conducted during September 2015. These training samples were used to derive the average EVI and LSWI time series of each land cover type.

Another 1000 random MODIS pixels (not shown) were selected from the cloudless or thin cloud regions of the Landsat-8 OLI images. They were used as validation samples to analyze the sensitivity of the decision tree model and to obtain optimal thresholds (see section 3.3.2). Other than the training samples, most of the validation samples were mixed pixels. The true class (rice or non-rice) of each validation sample was visually interpreted using Landsat-8 OLI images and the associated EVI/LSWI time series. The land cover type was determined as the dominant land use type showing the maximum percentage.

Methodology

To map rice planting areas using phenological parameters, we first analyzed the EVI and LSWI time series features of the different land cover types. Then, six phenological parameters were developed to distinguish between different land cover types. Subsequently, a step-wise decision tree model was designed to exclude non-rice land cover types based on the six phenological parameters. The sensitivity of the decision tree model to each parameter was investigated using the 1000 validation samples. Ultimately, a set of thresholds with the highest overall accuracy was used to derive the final paddy rice planting map.

Analysis of the EVI and the LSWI time series

Different land cover types had different EVI dynamics. The EVI time series of each land cover type is shown in Fig. 2(a). Autumn maize and soybean were combined into one class because of their spectral similarity in the EVI and the LSWI time series (^{Liu et al., 2014b}). The EVI values of water bodies and artificial surfaces were clearly lower than those of the natural vegetation and crops throughout the entire season. The EVI time series curve of natural vegetation had a gradual descending tendency while that of cotton fluctuated at approximately 0.4. The EVI peak values of summer maize and lotus were greater than 0.4. Their EVI curves first ascended and then rapidly descended. From June to November, because of seasonal variations in crop growth, the EVI time series of mid-season rice, late-season rice, autumn maize, and cash crops were typical unimodal curves. The seeds are sown during mid-June for these crops; the EVI gradually increased after the seedlings’ emergence during mid-July, and the EVI gradually decreased after reaching the growth peak during the heading period of mid- or late-August. It was difficult to effectively distinguish between these four crops because the EVI curves are very similar.

Figure 2(b) shows the LSWI time series of different land cover types. However, they are difficult to discern in Fig. 2(b). Because water bodies, artificial surfaces, cotton, lotus, and natural vegetation can be separated based on changes in the EVI, only the LSWI time series curves of autumn maize, mid-season rice, late-season rice, and cash crops are drawn in Fig. 2(c). As seen in Fig. 2(c), these four crops had very similar LSWI dynamics prior to mid-August, but they appeared different after that time. The LSWI of autumn maize and cash crops sharply decreased, while the LSWI of rice remained high. The LSWI of paddy rice rapidly decreased during mid-October.

Based on the aforementioned analysis, it is obvious that the EVI can effectively distinguish vegetation from no vegetation, and it can separate natural vegetation, summer maize, cotton, and lotus from paddy rice. The LSWI can efficiently identify paddy rice from autumn maize and cash crops.

　Calculation of phenological parameters

Six phenological parameters were extracted for paddy rice identification based on the phenological features of different land cover types. Adjacent time intervals were considered because of the changes in phenological characteristics with latitude and the noise in the smoothed EVI/LSWI time series. The six phenological parameters were labelled as EVI₁, EVI₂, LSWI₁, LSWI₂, Kurtosis_LSWI, and Skewness_LSWI.

EVI₁ was defined as the minimum EVI during the sowing period from June 10 to June 26 (Fig. 2(a)). During this period, crops are at the beginning of the sowing phase, when cropland has relatively sparse vegetation cover, and the EVI of the cropland is very near that of bare land. Additionally, natural vegetation has high EVI and water bodies have low EVI, which can be effectively identified by EVI₁.

The maximum/peak EVI (EVI₂) was defined as the maximum EVI from August 13 to September 14 (Fig. 2(a)). During this period, rice is at the heading phase and has the relatively highest EVI value. The EVI of natural vegetation decreases slightly, whereas the EVI value of water bodies and artificial surfaces show nearly no change. Artificial surfaces, natural vegetation, cotton, lotus, and summer maize can be separately identified from paddy rice based on the variations between EVI₁ and EVI₂.

The moisture duration (LSWI₁) was defined as the number of periods when the LSWI values were continuously greater than 0.2 from June 26 to November 25. Water supply plays an important role during the rice growing season and causes the LSWI of paddy rice to have a higher value than that of other crops during the entire growing season. To characterize this trait, we developed this parameter to indicate the moisture duration. Generally, the moisture duration of paddy rice is longer than 10 periods, and the moisture duration of autumn maize and cash crops is usually shorter than 9 periods.

LSWI₂ was defined as the total LSWI from September 6 to November 25. Paddy rice is usually harvested one month (mid-season rice) or two months (late-season rice) later than autumn maize and cash crops. Therefore, during the second half of the growing season of paddy rice, it presented a higher LSWI than that of the maize and cash crops. As defined, LSWI₂ can effectively distinguish between rice and non-rice crops.

Kurtosis is a statistic that describes the degree of steepness of distributions in a population. Kurtosis_LSWI was computed using Eq. (3) as follows:

(3){\mathrm{Kurtosis}}_{\mathrm{LSWI}}={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}{{\displaystyle \sum _{i=\mathrm{1}}^n\left({\mathrm{LSWI}}_i-\overline{\mathrm{LSWI}}\right)}}^{\mathrm{4}}/{\sigma }^{\mathrm{4}}-\mathrm{3},

where n is the number of periods, {\mathrm{LSWI}}_i is the LSWI value at i^th period, \overline{\mathrm{LSWI}} is the mean LSWI, and \sigma is the standard deviation of the LSWI time series. A perfectly normal distribution has zero kurtosis, and the higher the kurtosis, the sharper the peak of the distribution (^{Jensen, 2004}). A negative kurtosis value implies that the overall data distribution is flat. As shown in Fig. 2, the LSWI time series of autumn maize had a sharper curve than that of mid-season rice, late-season rice, and cash crops.

Skewness depicts the symmetry of the probability distribution. Skewness_LSWI was calculated using Eq. (4) as follows:

(4){\mathrm{Skewness}}_{\mathrm{LSWI}}={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}{{\displaystyle \sum _{i=\mathrm{1}}^n\left({\mathrm{LSWI}}_i-\overline{\mathrm{LSWI}}\right)}}^{\mathrm{3}}/{\sigma }^{\mathrm{3}},

where n is the number of periods, {\mathrm{LSWI}}_i is the LSWI value at i^th period, \overline{\mathrm{LSWI}} is the mean LSWI, and \sigma is the standard deviation of the LSWI time series. A perfectly normal distribution has a skewness value of zero. Skewness values greater than zero have a right skewed data distribution. Conversely, a skewness value less than zero indicates a left skew. The LSWI time series of autumn maize and cash crops appeared more symmetric than that of mid-season rice and late-season rice (Fig. 2(c)).

　Decision tree classification

　Separation of different land cover types

To develop an effective separation strategy to classify rice pixels, we also checked the distributions of all training samples in the EVI₁-EVI₂, LSWI₁-LSWI₂, and Kurtosis_LSWI-Skewness_LSWI parameter space (Figs. 3(a)–3(c)). In the EVI₁-EVI₂ parameter space, water bodies and artificial surfaces occupied very different regions from paddy rice and had much lower EVI₁ and EVI₂ values. Summer maize, natural vegetation, cotton, and lotus were aggregated in the other corner from that of the paddy rice and had higher EVI₁ values. It was difficult to separate autumn maize and cash crops from paddy rice in the EVI₁-EVI₂ parameter space. In the LSWI₁-LSWI₂ parameter space (Fig. 3(b)), most samples of the other land cover types were separable from paddy rice, except natural vegetation and summer maize. As shown in Fig. 3(c), the land cover types that could not be distinguished from paddy rice were water bodies and natural vegetation. However, as shown in Figs. 3(a)–3(c), mid-season rice and late-season rice were not clearly separable. Though they appeared separable in the LSWI₁-LSWI₂ space, the distinction was not sufficiently strong. Therefore, we did not differentiate mid-season rice from late-season rice in this classification.

The gray regions in Fig. 3(a), Fig. 3(d), and Fig. 3(e) indicate the separation of paddy rice from other land cover types. Summer maize, natural vegetation, cotton, and lotus appear separable from paddy rice (as seen in Fig. 3(a)). Considering the complexity of the land surface, mixed pixels inevitably exist. To ensure all pixels containing paddy rice were retained at the beginning of the classification, we relaxed the criteria. If a condition was set as EVI₁<0.4 and EVI₂>0.4, all training samples that fell in the gray region of Fig. 3(a) are shown in Fig. 3(d).

Some significantly different land cover types could be excluded using the LSWI criteria (Fig. 3(d) with Fig. 3(b)). If a condition was set to 9<LSWI₁<19 and LSWI₂>1.5, the only interference type was natural vegetation. Finally, using Kurtosis_LSWI>‒1.7 and Skewness_LSWI<0.0 (Fig. 3(e)), paddy rice could be distinguished from natural vegetation.

　Sensitivity analysis of the decision tree model

We utilized a decision tree model to classify paddy rice in this study and carefully determined the thresholds for the parameters. Based on the aforementioned analysis, a hierarchical decision tree model was designed to generate paddy rice maps (Fig. 4). There were seven parameters (para₁–para₇) in our decision tree model that needed values to execute the classification process. A trial set of thresholds (para₁=0.4, para₂=0.4, para₃=9, para₄=19, para₅=1.5, para₆=1.7, para₇=0.0) were defined based on the training sample statistics from the aforementioned analysis.

There were abundant mixed pixels that could not be neglected during the mapping process. To obtain the optimal thresholds in mapping paddy rice, we explored the sensitivity of the decision tree model to each parameter to improve the classification accuracy. During the sensitivity analysis, the threshold for a specific parameter was changed while the thresholds for the other parameters remained the same as in the trial set in each test. Then, a new set of thresholds were used to classify the 1000 validation samples, and the overall accuracy (^{Story and Congalton, 1986}) was recorded. For example, the threshold for para₁ was changed from 0.35 to 0.45 at a step of 0.01 while para₂–para₇ were maintained as (0.4, 9, 19, 1.5, 1.7, 0.0), respectively, and the accuracy of each test was recorded. The sensitivity of decision tree classification depends on its accuracy using different parameters. If the classification accuracy of a decision tree greatly changes when the threshold of one parameter varies, the decision tree classification model could be considered to be sensitive to the parameter. In contrast, if the classification accuracy shows no or only a slight change when the threshold of one parameter varies, the model could be considered to be insensitive to the parameter.

Results and discussion

　Sensitivity the decision tree model

The total accuracy of each test is listed in Table 2. All tests achieved satisfactory results with total accuracies higher than 74%, while the threshold set with the lowest accuracy of 74.8% was (0.4, 0.4, 14, 19, 1.5, ‒1.7, 0.0). The threshold set with the highest accuracy of 93.9% was (0.4, 0.42, 9, 19, 1.5, ‒1.7, 0.0). The accuracies for the parameter tests first increased and then decreased with an increase in the threshold value among all the tests, except for para₄ and para₆. As the para₄ threshold value increased, the overall classification accuracy first increased and then remained unchanged. For para₆, the classification decreased with an increase in the threshold value. Therefore, the decision tree classification was sensitive to the values of para₄, para₅, and para₇, and it was insensitive to changes in para₁, para₂, para₃, and para₆. Finally, the threshold set (0.4, 0.42, 9, 19, 1.5, ‒1.7, 0.0) was used for the final classification of paddy rice.

　Distribution of paddy rice

The paddy rice map of Jiangsu Province in 2015 is shown in Fig. 5. Paddy rice in Jiangsu Province was mainly concentrated in the northeastern and central parts of the province with massive and contiguous areas of cropland. Paddy rice farming in Jiangsu was concentrated on flat plains and areas with abundant water resources. During recent years, paddy rice production in Jiangsu Province has shifted northward. The rice planting area in northern Jiangsu Province has slightly expanded, but the expanded area is mainly on saline land along the coast. In the southern hilly area of Jiangsu Province, paddy rice has been planted in flat areas around rivers and lakes in a relatively dispersed pattern.

The latest version of datasets from the NLCD project were produced based on remotely sensed images from 2015 (NLCD₂₀₁₅) (^{Liu et al., 2005}; ^{Liu et al., 2014a}), which we acquired as our reference map. NLCD₂₀₁₅ was mainly based on visual interpretation and digitization of Landsat-8 OLI images, with the land cover classes including six primary classes and 25 secondary classes. The cropland classes included paddy land and dry land. The fraction of paddy land was calculated by resampling the vector layer of the paddy land to a 1-km resolution. The NLCD classification system defines paddy land as cropland that has sufficient water supply and irrigation facilities for planting paddy rice, lotus, etc. The MODIS-derived results from this study were generally consistent with the NLCD₂₀₁₅ paddy field distribution. The total paddy land in Jiangsu Province was 40,198 km² from NLCD₂₀₁₅ and the MODIS-derived paddy rice planting area was 42,193 km² in 2015 (assuming each mapped rice pixel was a 100% fraction), accounting for 87.9% and 92.2% of the cropland area, respectively.

Comparison of paddy rice maps at a regional scale

Regions A–D are four typical inconsistent regions from the MODIS-derived paddy rice map and the NLCD₂₀₁₅ paddy rice map. Region A is in an intensive paddy rice cultivation area. Nearly all paddy rice was identified on the MODIS-derived paddy rice map, while the NLCD₂₀₁₅ showed that there were nearly no paddy fields (Fig. 6). Previous research has shown that this area has been planted with rice (^{Xiao et al., 2005}; ^{Qiu et al., 2015}). The CBERS-04 images were acquired on October 20, 2015, when the paddy rice was mature and the chlorophyll content of its leaves was gradually decreasing. The paddy rice appears yellow-green in the pseudo color synthesis of the CBERS-04 image while the emerging winter wheat appears bright green (Fig. 6(a)). We compared the MODIS-derived paddy rice distribution to the paddy fields in the NLCD₂₀₁₅ with rice areas derived from the Support Vector Machine (SVM) classification of the CBERS-04 image (Fig. 6(b)). The MODIS-derived paddy rice distribution (Fig. 6(c)) was generally consistent while many paddy fields were missed in the NLCD₂₀₁₅.

Region B is in the intertidal zone along the eastern coast of Jiangsu Province, an area with muddy tidal flats (Fig. 7). In the past, there was no rice planting in this area; however, currently, with progress in saline land transformation, salt-tolerant rice cultivation, and direct rice seeding technology in saline soil, more land in this area has been successfully planted with rice. Planting paddy rice in tidal flats can alleviate the challenges of an increasing population and limited cropland. Therefore, paddy rice planting areas and paddy rice yields have increased. However, NLCD₂₀₁₅ classified most of this region as dry land while the MODIS-derived results mapped this area as planted with paddy rice. On the Landsat-8 OLI image acquired on October 13, 2015, parcels with a green color indicated paddy rice, and the surrounding area in yellow suggested bare soil after the crop was harvested. Hence, the MODIS-derived results more accurately represented the area of planted paddy rice.

Based on the Landsat-8 OLI image acquired on October 13, 2015 (Fig. 7), it is obvious that the central area within Region C is water bodies. Paddy rice is rarely grown in this region because small land parcels in this area are inconvenient for agricultural management and rice is frequently drowned during the rainy season. This area is mainly used for aquaculture and aquatic crop cultivation, particularly lotus planting. According to the definition of paddy land and the intention of our paddy rice map, both the MODIS-derived maps and the NLCD₂₀₁₅ land use classification were correct.

Region D is in southern Jiangsu Province, in the middle and lower reaches of the Yangtze River. This region is among the most economically developed areas in China. Therefore, industrial construction and city development cover many areas of former cropland, resulting in fragmented arable land. Urban agriculture meets the needs of urban residents via a complex cropping structure. Therefore, grain production has gradually decreased in this region, and the mixed pixels make paddy rice identification and extraction challenging. Our result and the NLCD₂₀₁₅ results were inconsistent in this area. Because of the spatial resolution and the cloud cover of the Landsat-8 OLI image, rice could not be distinguished from the other crops. High spatial resolution or hyperspectral images are needed to discern different crops and validate the results for this region.

Comparison of paddy rice maps at the county level

The paddy rice map was also evaluated at the county level using agricultural statistics from the Survey Office of the National Bureau of Statistics in Jiangsu Province (Fig. 8). It showed a close correlation with an R²=0.85 between the MODIS-derived paddy rice area and the reported area, while the R² between the paddy rice area from the NLCD and the reported county-level rice planting area was 0.50; thus, the paddy rice area derived from this study was much closer to the actual paddy rice planting acreage. However, the slope of the derived regression equation between the MODIS-derived paddy rice area and the reported area was 0.41, not near 1.0; thus, MODIS-derived paddy rice map tends to overestimate the paddy rice area. This overestimation of paddy rice area was caused by the mixed pixels, as the size of the parcels in southern China are smaller than the MODIS pixel size (^{Ozdogan and Woodcock, 2006}).

Conclusions

This study proposed a new method to map paddy rice by combining phenological parameters with a decision tree model. Six phenological parameters were used to identify paddy rice based on an analysis of MODIS, EVI, and LSWI time series of different land cover types. The six parameters were minimum EVI during the sowing period (EVI₁), maximum EVI at the peak (EVI₂), moisture duration (LSWI₁), accumulated LSWI (LSWI₂), kurtosis of the LSWI time series (Kurtosis_LSWI), and skewness of the LSWI time series (Skewness_LSWI). These phenological parameters were extracted based on the performance of land cover types during specific phenological phases (EVI₁ and EVI₂), one- half or the entire rice growing cycle (LSWI₁ and LSWI₂), and the shape of the LSWI time series (Kurtosis_LSWI and Skewness_LSWI). According to the distributions of the training samples in the EVI₁-EVI₂, LSWI₁-LSWI₂, and Kurtosis_LSWI-Skewness_LSWI spaces, we designed a hierarchical decision tree model to identify the paddy rice areas. The results showed that the decision tree model was sensitive to the values of para₄, para₅, and para₇, but was relatively insensitive to changes in para₁, para₂, para₃, and para₆. A paddy rice map of Jiangsu Province from 2015 was generated with the optimal threshold set of (0.4, 0.42, 9, 19, 1.5, ‒1.7, 0.0) with a total accuracy of 93.9%. The MODIS-derived paddy rice map generally agreed with the paddy rice map derived from the NLCD project, but there were regional discrepancies because of their different definitions and the inability of MODIS to map paddy rice in a fragmented landscape. The MODIS-derived paddy rice map showed a high correlation (R² =0.85) with county-level agricultural statistics.

Two problems remain to be resolved in future studies. First, the coarse resolution of the MODIS data lowered the classification accuracy and resulted in overestimation of the rice planting area. We also found that the thresholds in the decision-tree classification model were important for classification accuracy because cropping systems differ from region to region, and even the same crop planted in different regions can show different phenological characteristics as it adjusts to the local climate. Although we could not provide universal optimal thresholds for other regions, the method described in this paper is quite adaptable to other regions, as the growth cycles of paddy rice are quite similar. In addition, the sensitivity analysis results can serve as an important guideline and help those who are interested in the method to obtain optimal results in their regions. When applying this method, one should pay more attention to setting the thresholds for para₄, para₅, and para₇. If there are sufficient ground samples, one can determine optimal thresholds using the method described in this study. When combined with available optimization algorithms and emerging deep learning technology, this method could be used to map rice planting areas at larger scales.

References

[1] Bouman B (2009). How much water does rice use? Rice Today, (1): 15

[2] Bradley B A, Mustard J F (2008). Comparison of phenology trends by land cover class: a case study in the Great Basin, USA. Glob Change Biol, 14(2): 334–346

[3] Dong J, Xiao X (2016). Evolution of regional to global paddy rice mapping methods: a review. ISPRS J Photogramm Remote Sens, 119: 214–227

[4] Dvorak W S (2012). Water use in plantations of eucalypts and pines: a discussion paper from a tree breeding perspective. Int Rev, 14(1): 110–119

[5] Friedl M A, Brodley C E (1997). Decision tree classification of land cover from remotely sensed data. Remote Sens Environ, 61(3): 399–409

[6] Friedl M A, Sulla-Menashe D, Tan B, Schneider A, Ramankutty N, Sibley A, Huang X (2010). MODIS Collection 5 global land cover: algorithm refinements and characterization of new datasets. Remote Sens Environ, 114(1): 168–182

[7] Frolking S, Qiu J, Boles S, Xiao X, Liu J, Zhuang Y, Li C, Qin X (2002). Combining remote sensing and ground census data to develop new maps of the distribution of rice agriculture in China. Global Biogeochem Cycles, 16(4): 1091

[8] Huete A, Didan K, Miura T, Rodriguez E P, Gao X, Ferreira L G (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens Environ, 83(1‒2): 195–213

[9] Jensen J R (2004). Introductory Digital Image Processing: A Remote Sensing Perspective (3nd ed). Englewood Cliffs: Prentice Hall

[10] Kimball J S, McDonald K C, Running S W, Frolking S E (2004). Satellite radar remote sensing of seasonal growing seasons for boreal and subalpine evergreen forests. Remote Sens Environ, 90(2): 243–258

[11] Labus M P, Nielsen G A, Lawrence R L, Engel R, Long D S (2002). Wheat yield estimates using multi-temporal NDVI satellite imagery. Int J Remote Sens, 23(20): 4169–4180

[12] Liu J, Kuang W, Zhang Z, Xu X, Qin Y, Ning J, Zhou W, Zhang S, Li R, Yan C, Wu S, Shi X, Jiang N, Yu D, Pan X, Chi W (2014a). Spatiotemporal characteristics, patterns, and causes of land-use changes in China since the late 1980s. J Geogr Sci, 24(2): 195–210

[13] Liu J, Liu M, Tian H, Zhuang D, Zhang Z, Zhang W, Tang X, Deng X (2005). Spatial and temporal patterns of China’s cropland during 1990–2000: an analysis based on Landsat TM data. Remote Sens Environ, 98(4): 442–456

[14] Liu J, Pan Y, Zhu X, Zhu W (2014b). Using phenological metrics and the multiple classifier fusion method to map land cover types. J Appl Remote Sens, 8(1): 083691

[15] Liu J, Zhu W, Cui X (2012). A Shape-matching Cropping Index (CI) mapping method to determine agricultural cropland intensities in China using MODIS time-series data. Photogramm Eng Remote Sensing, 78(8): 829–837

[16] Loveland T R, Belward A S (1997). The IGBP-DIS global 1 km land cover data set DISCover: first results. Int J Remote Sens, 18(15): 3289–3295

[17] Matthews H D, Caldeira K (2007). Transient climate-carbon simulations of planetary geoengineering. Proc Natl Acad Sci USA, 104(24): 9949–9954

[18] Otukei J R, Blaschke T (20 10). Land cover change assessment using decision trees, support vector machines and maximum likelihood classification algorithms. Int J Appl Earth Obs Geoinf, 12(Supplement 1): S27–S31

[19] Ozdogan M, Woodcock C E (2006). Resolution dependent errors in remote sensing of cultivated areas. Remote Sens Environ, 103(2): 203–217

[20] Pal M, Mather P M (2003). An assessment of the effectiveness of decision tree methods for land cover classification. Remote Sens Environ, 86(4): 554–565

[21] Pan Y, Li L, Zhang J, Liang S, Zhu X, Sulla-Menashe D (2012). Winter wheat area estimation from MODIS-EVI time series data using the crop proportion phenology index. Remote Sens Environ, 119(3): 232–242

[22] Peng D, Huete A R, Huang J, Wang F, Sun H (2011). Detection and estimation of mixed paddy rice cropping patterns with MODIS data. Int J Appl Earth Obs Geoinf, 13(1): 13–23

[23] Pongratz J, Lobell D, Cao L, Caldeira K (2012). Crop yields in a geoengineered climate. Nat Clim Chang, 2(2): 101–105

[24] Potgieter A B, Apan A, Hammer G, Dunn P (2010). Early-season crop area estimates for winter crops in NE Australia using MODIS. ISPRS J Photogramm Remote Sens, 65(4): 380–387

[25] Potgieter A B, Lawson K, Huete A R (2013). Determining crop acreage estimates for specific winter crops using shape attributes from sequential MODIS imagery. Int J Appl Earth Obs Geoinf, 23(8): 254–263

[26] Pringle M J, Denham R J, Devadas R (2012). Identification of cropping activity in central and southern Queensland, Australia, with the aid of MODIS MOD13Q1 imagery. Int J Appl Earth Obs Geoinf, 19(1): 276–285

[27] Qin Y, Xiao X, Dong J, Zhou Y, Zhu Z, Zhang G, Du G, Jin C, Kou W, Wang J, Li X (2015). Mapping paddy rice planting area in cold temperate climate region through analysis of time series Landsat 8 (OLI), Landsat 7 (ETM+) and MODIS imagery. ISPRS J Photogramm Remote Sens, 105: 220–233

[28] Qiu B, Li W, Tang Z, Chen C, Qi W (2015). Mapping paddy rice areas based on vegetation phenology and surface moisture conditions. Ecol Indic, 56: 79–86

[29] Robock A, Oman L, Stenchikov G L (2008). Regional climate responses to geoengineering with tropical and arctic SO₂ injections. J Geophys Res, 113: D16101

[30] Sakamoto T, Van Phung C, Kotera A, Nguyen K D, Yokozawa M (2009). Analysis of rapid expansion of inland aquaculture and triple rice-cropping areas in a coastal area of the Vietnamese Mekong Delta using MODIS time-series imagery. Landsc Urban Plan, 92(1): 34–46

[31] Story M, Congalton R G (1986). Accuracy assessment: a user’s perspective. Photogramm Eng Remote Sensing, 52(3): 397–399

[32] Tooke T R, Coops N C, Goodwin N R, Voogt J A (2009). Extracting urban vegetation characteristics using spectral mixture analysis and decision tree classifications. Remote Sens Environ, 113(2): 398–407

[33] Wardlow B D, Egbert S L (2008). Large-area crop mapping using time-series MODIS 250m NDVI data: an assessment for the U.S. Central Great Plains. Remote Sens Environ, 112(3): 1096–1116

[34] Xiao X, Boles S, Frolking S, Li C, Babu J Y, Salas W, Moore B III (2006). Mapping paddy rice agriculture in South and Southeast Asia using multi-temporal MODIS images. Remote Sens Environ, 100(1): 95–113

[35] Xiao X, Boles S, Liu J, Zhuang D, Frolking S, Li C, Salas W, Moore B III (2005). Mapping paddy rice agriculture in southern China using multi-temporal MODIS images. Remote Sens Environ, 95(4): 480–492

[36] Zhang G, Xiao X, Dong J, Kou W, Jin C, Qin Y, Zhou Y, Wang J, Menarguez M A, Biradar C (2015). Mapping paddy rice planting areas through time series analysis of MODIS land surface temperature and vegetation index data. ISPRS J Photogramm Remote Sens, 106: 157–171

[37] Zhu W, Pan Y, He H, Wang L, Mou M, Liu J (2012). A changing-weight filter method for reconstructing a high-quality NDVI time series to preserve the integrity of vegetation phenology. IEEE Trans Geosci Remote Sens, 50(4): 1085–1094

[38] Zou J, Huang Y, Zheng X, Wang Y (2007). Quantifying direct N₂O emissions in paddy fields during rice growing season in mainland China: dependence on water regime. Atmos Environ, 41(37): 8030–8042

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature