Monitoring Spatio-Temporal Distribution of Rice Planting Area in the Yangtze River Delta Region Using MODIS Images

Abstract

A large-area map of the spatial distribution of rice is important for grain yield estimations, water management and an understanding of the biogeochemical cycling of carbon and nitrogen. In this paper, we developed the Normalized Weighted Difference Water Index (NWDWI) for identifying the unique characteristics of rice during the flooding and transplanting period. With the aid of the ASTER Global Digital Elevation Model and the phenological data observed at agrometeorological stations, the spatial distributions of single cropping rice and double cropping early and late rice in the Yangtze River Delta region were generated using the NWDWI and time-series Enhanced Vegetation Index data derived from MODIS/Terra data during the 2000–2010 period. The accuracy of the MODIS-derived rice planting area was validated against agricultural census data at the county level. The spatial accuracy was also tested based on a land use map and Landsat ETM+ data. The decision coefficients for county-level early and late rice were 0.560 and 0.619, respectively. The MODIS-derived area of late rice exhibited higher consistency with the census data during the 2000–2010 period. The algorithm could detect and monitor rice fields with different cropping patterns at the same site and is useful for generating spatial datasets of rice on a regional scale.

1. Introduction

Paddy rice fields provide essential food for more than half of the population of the entire world [1]. Rice is widely cultivated in Asian countries, especially China. Recent FAO (Food and Agricultural Organization) estimates indicate that to satisfy the projected demand of the year 2050, global agricultural production must increase 60 percent above the level of 2005–2007 [2]. Large-area assessments of potential food production regions and their impact on biogeochemical cycling require the acquisition of the best possible information on the distribution of paddy rice fields [3].

The official statistical data on rice sowing areas have been generated based on ground sample surveys and extrapolated to the provincial and national scales. Large-scale census data cannot provide accurate spatial distributions of paddy rice, and a time lag is present in the datasets. Huke developed Asian rice datasets using agricultural statistical datasets collected at the sub-country level [4]. Leff et al. generated a global rice map at a spatial resolution of five arcminutes as part of a global cropland product using satellite-derived land cover data and agricultural census data [5]. Frolking et al. generated 0.5°-resolution maps of the distribution of rice agriculture in mainland China using a combination of county-scale agricultural census data and land cover maps derived from Landsat images collected during the 1995–1996 period [3]. A thematic land use map of China at a scale of 1:100,000 was generated via the visual interpretation of Landsat TM (Thematic Mapper) data [6]. A classification system of 25 land use categories, including paddy rice, was used in this work. The land use map was converted to 1-km gridded data. However, more updated datasets of annual rice distribution with finer resolution are needed at the regional scale.

Approximately half of the cropland in China is multi-cropped each year, and this land has a significant influence on the biogeochemical cycling of carbon and nitrogen. To date, many studies have been conducted to map paddy rice using fine-resolution satellite images, such as Landsat MSS, TM, ETM+ and NOAA/AVHRR images, by applying image classification procedures, but few of these studies have provided detailed information regarding the locations of multi-cropping [7,8,9,10]. Furthermore, because of the fine resolution of these images, it is difficult to obtain more comprehensive images covering an entire region simultaneously over a large area. Rice distribution maps have also been produced via multi-temporal analysis of NOAA/AVHRR and SPOT4/VEGETATION data with a resolution of ~1 km, which is rather coarse for rice mapping [11,12,13].

The Moderate-Resolution Imaging Spectroradiometer (MODIS) aboard the Terra and Aqua satellites, with its advantages of a high revisit period, moderate spatial resolution, wide field of view (FOV) and free access, has been applied for paddy rice mapping. Decision tree algorithms and spectral matching techniques were used to map rice-growing areas using temporal MODIS data [14,15]. A MODIS time-series analysis of spectral indices was found to be more useful for monitoring the phenological variations of paddy rice over a long period [16,17]. A paddy rice field is typically prepared by flooding a few days before the rice seedlings are transplanted. The wet growing season is regarded as a unique and significant characteristic of rice compared with other crops [18]. Thus, the flooding period is recognized as the best phase for rice identification. Spectral indices and bands that are sensitive to water and green vegetation are needed for monitoring the flooding and transplanting period of rice crops. These spectral indices are always calculated using two or more spectral bands to enhance the contrast between target and background and to reduce the effects of the atmosphere and solar illumination geometry. The Normalized Difference Vegetation Index (NDVI), developed using the red and near-infrared bands, is correlated with the Leaf Area Index and chlorophyll content and has been widely used for crop yield estimations and the detection of changes in land use/cover [19,20]. The Enhanced Vegetation Index (EVI) was proposed because of the saturation of the NDVI in high-biomass regions to adjust for residual atmospheric contamination and background reflectance [21,22]. The infrared range is useful for estimating the water content of vegetation and in discriminating water from land. The Normalized Difference Water Index (NDWI) was developed using the reflectance in the near-infrared and green bands to enhance the detection of water features while eliminating soil and terrestrial vegetation features [23]. The modified NDWI (MNDWI) substitutes the near-infrared band with a middle-infrared band, such as Band 5 of Landsat TM, to efficiently enhance open-water signals and suppress or remove the signals from built-up land, as well as vegetation and soil [24]. However, because the reflectance of rice pixels during the flooding and transplanting stage is a mixture of water and vegetation, the sensitivity of the spectral index to flooding features should be further improved for rice mapping. The Land Surface Water Index (LSWI), which was formulated by combining the red and shortwave infrared channels of MODIS, has been used for the identification of rice pixels [25]. However, the threshold between the LSWI and the EVI was determined by considering local practices and rice cropping systems. Qiu et al. proposed a method for mapping rice planting areas by considering the vegetation phenology and surface water variations. The ratios of the changes in amplitude of the LSWI to the two-band Enhanced Vegetation Index 2 (EVI2) during the period from the tillering to the heading stage were used as one indicator to discriminate rice from non-rice fields [26]. Mosleh and Hassan developed a method for mapping “Boro” rice in Bangladesh using the MODIS-derived 16-day composite NDVI at a spatial resolution of 250 m [27]. The ISODATA clustering and the formulation of the mathematical model were the key procedures of this algorithm.

The objectives of the present study are to: (1) develop a Normalized Weighted Difference Water Index for identifying the flooding period of paddy rice fields; (2) map the early, single cropping and double cropping late rice distributions of the Yangtze River Delta region in the 2000–2010 period; and (3) validate the results using land use maps, Landsat ETM+ data and agricultural statistical data.

2. Study Area and Data

2.1. Study Area

The study area is the Yangtze River Delta region, which is one of the major rice-producing areas in China and spans three provinces (Figure 1). This region extends from 118°50ʹ5″E to 134°46ʹ26″E in longitude and from 38°43ʹ15″N to 53°33ʹ39″N in latitude, with a territory of 2.1 × 10⁵ km². The climate of the Yangtze River Delta is humid subtropical and is largely controlled by the East Asian monsoon [28]. Rice is the major food crop in the study area, with a high level of production and a wide distribution. The cropping system in Jiangsu Province and Shanghai City consists essentially of one crop of rice and another crop of winter wheat or oil rape, whereas single and double rice cropping systems are the two major planting patterns in Zhejiang Province.

Figure 1. Location and ASTER Global Digital Elevation Model (GDEM) of the study area.

Figure 1. Location and ASTER Global Digital Elevation Model (GDEM) of the study area.

2.2. Data Acquisition

2.2.1. Field Data

Field experiments can yield accurate data under controlled conditions. A field experiment was conducted at the Experimental Farm of Zhejiang University, Hangzhou, Zhejiang Province, from June to October in 2004. Two rice cultivars (i.e., Xieyou 9308 and Xiushui 110) were planted in 18 plots with three different nitrogen fertilization treatments: 0, 140 and 240 kg/ha. Each treatment was repeated three times. Rice seedlings were transplanted into the field on 8 July 2004, and the canopy reached full closure in August. The rice canopy reflectance of each plot was acquired using an Analytical Spectral Devices (ASD) Field Spec Pro Full Range (350–2500 nm) spectroradiometer on 20 July, 8 August, 28 August, 22 September, 5 October and 27 October. At each plot, 10 reflectance measurements were acquired with a nadir view of 25° from a height of 1.0 m above the rice. The spectrum of each plot was recorded as the average of the 10 measurements.

2.2.2. Satellite Data

The MODIS sensor records data in 36 spectral bands and products at spatial resolutions of 250 m, 500 m and 1000 m. In this study, MODIS/Terra eight-day composite surface reflectance products (MOD09A1) were chosen for the mapping of rice planting regions.

The eight-day composite surface reflectance products were routinely processed for atmospheric and radiometric correction for the effects of aerosols and cirrus clouds, as well as to select the best observation and the lowest value in the blue band for each pixel over the eight-day period [29]. Three tiles (h27v05, h27v06 and h28v06) for the 2000–2010 period were acquired from the project website (https://lpdaac.usgs.gov/). The downloaded MODIS data were then mosaicked and reprojected to Albers equal-area conic projection using the MODIS Reprojection Tool (MRT).

A Landsat ETM+ image acquired on 13 May 2000 (path/row: 118/41) was downloaded from the International Scientific Data Service Platform (http://datamirror.csdb.cn/). The region spanned by the image covers the main rice-producing zones in Wenzhou City. According to the phenological data recorded at the local agricultural meteorological station, rice seedlings were generally transplanted into the fields in early May. Radiometric calibration was applied to the Landsat ETM+ image. The fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) model was selected for atmospheric correction. The Landsat ETM+ image was resized to a 90 km × 90 km subset and reprojected to the Albers equal-area conic projection.

In addition, ASTER Global Digital Elevation Model (GDEM) data covering the study area were freely obtained from the Earth Remote Sensing Data Analysis Center of Japan (http://gdem.ersdac.jspacesystems.or.jp/).

2.2.3. Ancillary Data

A digital administrative map of China was obtained from the National Fundamental Geographic Information System. A land use map of Wenzhou City in 2005 was obtained from the Land and Resources Bureau of Wenzhou City. The annual sowing areas of paddy rice for each county in the study area during the 2000–2010 period were provided by the bureau of statistics.

3. Methodology

3.1. Spectral Characteristics of Rice during the Flooding and Transplanting Period

Canopy reflectance data collected in the field at a spectral resolution of 1 nm were used to simulate the reflectance in the first seven bands of the MODIS sensor (ρ_MOD) based on its spectral response function (Figure 2). The reflectance in the near-infrared and shortwave infrared wavelength bands was very low, whereas the reflectance in the visible bands (Bands 1, 3 and 4) was greater than in other growth periods. During the transplanting period, the water in the rice field was found to absorb most of the incident radiant flux, especially in the shortwave infrared region. It was also observed that the reflectance in Band 6 was lower than that in Band 4 (green band). With an increase in the tiller number and leaf area index, the reflectance in the visible bands decreased on 8 August 2004; however, the reflectance in the near-infrared and shortwave infrared bands increased significantly. The reflectance in Band 6 became higher than that in Band 4. In the previous literature, many water indices have been developed using the visible and infrared bands [23,24,25]. The green spectral range is highly sensitive to the Chl-a concentration over a wide range of variation and, thus, is helpful for the remote sensing of vegetation [30]. The near-infrared and shortwave infrared regions are the best wavelength regions for discriminating land from water. Because 1 of the 20 detectors in Terra MODIS Band 5 is noisy, there are stripes in the image. Band 6 was selected as the band sensitive to water, and Band 4 was used as the band sensitive to the presence of green seedlings.

Xu proposed the MNDWI (Equation (1)) to enhance the features of open water in remotely-sensed imagery [24]. Water pixels will have positive values of this index; however, pixels corresponding to flooded rice fields, built-up land and vegetation will have negative values. To enhance the features of rice pixels in the transplanting stage, we introduced a weight in the green band. Thus, we developed a Normalized Weighted Difference Water Index (NWDWI) based on the MNDWI. The ρ_band4 and ρ_band6 values described in Equation (2) denote the reflectances in MODIS Bands 4 and 6, respectively. A threshold of zero was applied to the NWDWI to separate flooded rice pixels from vegetation pixels. As shown in Equation (3), the values of the Ratio Vegetation Index (RVI) during different rice growth periods were calculated using the ρ_MOD values of Bands 6 and 4. Figure 3 shows that the RVI was the lowest during the transplanting stage and reached its peak at the heading stage, after which the RVI slowly deceased.

$$\text{MNDWI=}\frac{{\text{ρ}}_{\text{green}}-{\text{ρ}}_{\text{MIR}}}{{\text{ρ}}_{\text{green}}\text{+}{\text{ρ}}_{\text{MIR}}}$$

(1)

$$\text{NWDWI=}\frac{{\text{ρ}}_{\text{band}6}-a\cdot {\text{ρ}}_{\text{band4}}}{{\text{ρ}}_{\text{band}6}\text{+}a\cdot {\text{ρ}}_{\text{band4}}}$$

(2)

$$\text{RVI= }\frac{{\text{ρ}}_{\text{band6}}}{{\text{ρ}}_{\text{band4}}}$$

(3)

Figure 2. Simulated MODIS reflectances throughout the entire growth cycle of rice.

Figure 2. Simulated MODIS reflectances throughout the entire growth cycle of rice.

Figure 3. Statistical analysis of the Ratio Vegetation Index (RVI) throughout the entire growth cycle of a rice paddy in the field.

Figure 3. Statistical analysis of the Ratio Vegetation Index (RVI) throughout the entire growth cycle of a rice paddy in the field.

The acquisition date of the ETM+ image was consistent with the flooding and transplanting date for rice in Wenzhou, Zhejiang Province. Figure 4 is a false-color composite image of the ETM+ data. The paddy rice binary map was obtained from the ETM+ image using the maximum likelihood method and then degraded to the same resolution as the MODIS data using the pixel aggregate method. Because the overpass time of Landsat 7 is close to that of Terra, the aggregated rice map could be used as ground-truth data for validation. The MOD09A1 data (day of year: 2000129) were used to validate the performance of the NWDWI. When a = 1, flooded rice field and vegetation pixels both had positive NWDWI values. When a = 1.5, 54.8% of the rice pixels had negative values, and 21.2% of the negative pixels in the MODIS-derived NWDWI image were labeled as rice pixels in the reference map. When a = 2, approximately 91.4% of the rice pixels had negative NWDWI values, and 24.7% of the negative pixels in the MODIS-derived NWDWI image were labeled as rice pixels in the reference map. The results of a simple density slice classification for the NWDWI (a = 2) demonstrated its ability to discriminate water pixels (Figure 5). In the ocean areas, the pixels had the lowest NWDWI values, whereas in the non-water areas, the pixels had positive values. In a comparison with the rice map derived from the ETM+ image, the omitted rice pixels were mainly distributed at the edges of the rice fields because of the mixed-pixel phenomenon and the uncertainty of edge pixels near large tracts of rice fields. When a = 2.5, approximately 98.5% of the rice pixels had negative NWDWI values. However, there were only 18.1% of pixels with negative values in the NWDWI image labeled as rice pixels in the reference map. Many of the forest, shrub and bare land pixels had negative value, as well as flooded rice pixels. Figure 3 also shows that the RVI was greater than 2 from the tillering stage to the harvest stage. Therefore, when a = 2, NWDWI ≤0 can be used to identify possible flooded rice pixels. Figure 5 also indicates that built-up pixels and pixels corresponding to natural water bodies also had low NWDWI values, which should allow them to be distinguished from rice pixels.

Figure 4. Color composite image of ETM+ data at the test site (R: Band 5; G: Band 4; B: Band 3).

Figure 4. Color composite image of ETM+ data at the test site (R: Band 5; G: Band 4; B: Band 3).

Figure 5. Spatial Normalized Weighted Difference Water Index (NWDWI) distribution of the MODIS image.

Figure 5. Spatial Normalized Weighted Difference Water Index (NWDWI) distribution of the MODIS image.

3.2. Reconstruction of the Spectral Index Profile

Although the 8-day composite surface reflectance data were routinely processed, some pixels were still affected by clouds. The pixels with Band 3 reflectances of greater than 10% were labeled as cloud pixels and removed as abnormal data [16]. Cloud masks for each MOD09A1 image were generated individually. Cloud-free data are important requirements for the operational monitoring of rice distributions using optical sensors [31]. To fill in the gaps in the EVI and NWDWI time series caused by clouds, the conditional temporal interpolation method was used in this study [32]. Compared with wavelet analysis and the Savitzky–Golay filter, the advantage of this method is its ability to retain the values of good pixels and repair the bad ones using valid pixels in the previous and subsequent images. If a pixel was contaminated in all three adjacent images, it was removed for further analysis.

3.3. Algorithm for Mapping Rice Planting Areas Using Time Series MODIS Data

The Yangtze River Delta region can be separated into two zones: one is the single rice planting area, including Jiangsu and Shanghai, and the other is Zhejiang Province [33]. In hilly regions, the elevation and slope are considered to be two important geographical factors for improving the accuracy and stability of classification in rice mapping [34]. In the study area, rice is grown in regions with elevations of less than 800 m and slopes of less than 10°; thus, the elevation and slope were used to exclude non-paddy rice regions.

According to the rice growth calendar for the period of 2000–2010 collected from agrometeorological stations, the transplanting stages for early rice, single cropping rice and double cropping late rice occurred in late April–early May, mid-June–early July and mid-to-late July, respectively. The MODIS data corresponding to the transplanting periods were used to identify flooded rice pixels using the NWDWI to reduce the interference of other wetland plants or crops with short-term precipitation.

The time series of spectral indices are essential for analyzing the annual variability of vegetation activityFigure 6 shows the seasonal EVI and NWDWI profiles of various types of land cover in 2005. Figure 6a–c shows the seasonal spectral index profiles of three major cropping systems in the Yangtze River Delta region. It is obvious that the NWDWI decreased significantly when a pixel was labeled as a cloud pixel. The time series of the EVI and NWDWI revealed the growth stages of the crops. When the rice seedlings were transplanted into the field, the NWDWI was less than zero, because the reflectance of the rice pixels was dominated by water. However, the seasonal NWDWIs of rain-fed crops maintained consistently positive values. Forests in the study area exhibited high EVI and NWDWI values throughout the entire growth period (Figure 6e). Natural water and built-up pixels could be distinguished by their long periods of consistently low NWDWI and EVI values. In the study area, 40 single cropping rice samples and 35 double cropping rice samples were selected to perform a decision tree classification algorithm for mapping the planting areas of early, single cropping and double cropping late rice. For early rice, a pixel with a negative NWDWI and EVI <0.26 was labeled as a potential rice pixel during the transplanting stage. According to the characteristics of the growth period of early rice, the maximum EVI throughout the entire growth period was greater than 0.35, and the EVI decreased below 0.35 in the eleventh 8-day period after that identified as the transplanting stage. Because of the longer growth period of single cropping rice, the EVI decreased to less than 0.35 in the fifteenth 8-day period after the rice seedlings were transplanted. Late rice was transplanted to the same field after the early rice was harvested. If a pixel had NWDWI <0.05 and EVI <0.35, it was recognized as a transplanted rice pixel. The maximum EVI throughout the entire growth period was higher than 0.35, and in the twelfth 8-day period after the transplanting stage, the EVI decreased below 0.35. Natural water bodies and built-up pixels could be excluded by removing pixels with an EVI that was less than 0.3 for fourteen consecutive 8-day periods between April and September. Figure 7 shows the decision tree for the mapping of single cropping, early and double cropping late rice planting areas. Finally, the spatial distribution maps of early, single cropping and double cropping late rice were generated.

Figure 6. Spectral index time series corresponding to various types of land cover at the test sites in 2005: (a) winter wheat and single cropping rice, (b) winter wheat and rain-fed crops, (c) early and late rice, (d) lakes, (e) forests and (f) built-up regions.

Figure 6. Spectral index time series corresponding to various types of land cover at the test sites in 2005: (a) winter wheat and single cropping rice, (b) winter wheat and rain-fed crops, (c) early and late rice, (d) lakes, (e) forests and (f) built-up regions.

Figure 7. Flow chart for the extraction of rice planting regions.

Figure 7. Flow chart for the extraction of rice planting regions.

3.4. Accuracy Assessment

The classification results for the Landsat ETM+ imagery at the studied site were used as the reference rice map for validation. An error confusion matrix was applied to evaluate the agreement of the MODIS-derived rice map with the reference rice map. The commission error, omission error, user’s accuracy and producer’s accuracy were calculated as follows:

$$\text{Commission error (\%) =}\frac{N_{\text{commit}}}{N_{\text{MODIS}}}\times 100\%$$

(4)

User’s accuracy (%) = 100 − commission error

(5)

$$\text{Omission error (\%) =}\frac{N_{\text{omit}}}{N_{\text{ETM}}}\times 100\%$$

(6)

Producer’s accuracy (%) = 100 − omission error

(7)

Here, N_commit and N_omit represent the numbers of committed and omitted rice pixels, respectively, in the MODIS-derived result, and N_MODIS and N_ETM represent the numbers of rice pixels in the MODIS-derived map and the aggregated reference rice map, respectively.

The error matrix was analyzed at the pixel level. Furthermore, because of the edge effects originating from the spatial aggregation of the ETM+ data and the geometric mismatch between the ETM+ and MODIS data, the error matrix was calculated using a 3 × 3 moving window [35].

4. Results and Discussion

4.1. Spatial and Temporal Distribution of Rice Planting Areas in the Yangtze River Delta Region

The spatial distributions of early, single cropping and double cropping late rice planting areas in the Yangtze River Delta region from 2000 to 2010 were generated using the presented algorithm, and the results are shown in Figure 8, Figure 9 and Figure 10. As shown in Figure 8, single cropping rice was mainly distributed in Jiangsu, Shanghai, Hangzhou-Jiaxing-Huzhou plain, Jinhua-Quzhou basin, Ningbo-Shaoxing plain and the coastal plain of southeastern Zhejiang Province. Double cropping rice was mainly distributed in the Jinhua-Quzhou basin, Ningbo-Shaoxing plain and coastal plain of southeastern Zhejiang Province (Figure 9 and Figure 10). The complexity of the terrain posed a considerable challenge in the extraction of scattered rice fields using MODIS 500-m data because of the mixed-pixel phenomenon. The rice fields were scattered throughout hilly regions, and most of them were distributed along rivers or in terraced planting regions.

The total planting area of early rice decreased from one year to the next from 2000 to 2003 and then remained stable afterward. Single cropping rice began to be planted instead of double cropping rice in some areas of the Yangtze River Delta region. Figure 10 shows that the early rice area decreased significantly in Ningbo-Shaoxing plain and Jinhua-Quzhou basin. The results reveal the change in the cropping systems used in the Yangtze River Delta region over the studied decade.

Figure 8. Spatio-temporal distribution of single cropping rice in the Yangtze River Delta region during the period of 2000–2010.

Figure 8. Spatio-temporal distribution of single cropping rice in the Yangtze River Delta region during the period of 2000–2010.

Figure 9. Spatio-temporal distribution of early rice in the Yangtze River Delta region during the period of 2000–2010.

Figure 9. Spatio-temporal distribution of early rice in the Yangtze River Delta region during the period of 2000–2010.

Figure 10. Spatio-temporal distribution of double cropping late rice in the Yangtze River Delta region during the period of 2000–2010.

Figure 10. Spatio-temporal distribution of double cropping late rice in the Yangtze River Delta region during the period of 2000–2010.

4.2. Comparison of the Estimated Paddy Rice Planting Areas with Agricultural Census Data

It is a time-consuming and labor-intensive task to implement a large-scale regional survey of rice planting region and to obtain an annual spatial map of the study region. Agricultural census data were used as reference data to test the accuracy and stability of our algorithm. Table 1 presents a comparison of the annual total paddy planting areas derived from MODIS and agricultural census data in the Yangtze River Delta region. The absolute errors of the extracted annual total rice areas were less than 15%, except for 2007 and 2010. The MODIS-derived single cropping rice areas were underestimated in 2007 and 2010. The relative error was highest in 2007. The rice planting areas in this year were severely underestimated, especially those in the south of Jiangsu Province, Shanghai City and the north of Zhejiang Province. In 2010, underestimation mainly occurred in northern Nantong City, Taizhou, western Wuxi, northeastern Suzhou and central-southern Zhejiang. In these regions, the cloud occurrence frequency was greater than 60% during the transplanting stage of single cropping rice. Cloud cover during the rainy season may obscure optical observations. The existence of clouds and cloud shadows can result in abnormal changes in the spectral index. Continuous cloud contamination during the transplanting period was the major cause of the underestimation. Although eight-day composite surface reflectance products were generated by selecting the date within the eight-day window with the clearest atmospheric conditions for each pixel, the effects of cloud contamination cannot be neglected. In this study, the conditional temporal interpolation method was applied to reconstruct invalid pixels contaminated by clouds, but if three consecutive eight-day composite data points were all invalid during the flooding and transplanting period of the rice crop, that pixel was eliminated from further analysis. Radar images are a potential alternative means of rice mapping in these regions, especially during the rainy season, because they are independent on the time of day and unimpaired by weather conditions.

Table 1. Comparison of rice planting areas derived using the MODIS algorithm and from the agricultural census data.

Year	Census Data (kha)	* RiceMOD (kha)	Relative Error (%)
2000	3928.50	3882.06	−11.8
2001	3476.30	3753.93	7.99
2002	3239.83	3472.60	7.18
2003	2842.23	3161.56	11.24
2004	3285.70	3576.88	8.86
2005	3409.26	3609.66	5.88
2006	3439.45	3638.51	5.79
2007	3338.49	2140.38	−35.89
2008	3359.37	2893.90	−13.86
2009	3335.40	3546.10	6.32
2010	3309.10	2636.50	−20.33

* Rice_MOD denotes the rice planting area derived from MODIS.

Table 1. Comparison of rice planting areas derived using the MODIS algorithm and from the agricultural census data.

Year	Census Data (kha)	* RiceMOD (kha)	Relative Error (%)
2000	3928.50	3882.06	−11.8
2001	3476.30	3753.93	7.99
2002	3239.83	3472.60	7.18
2003	2842.23	3161.56	11.24
2004	3285.70	3576.88	8.86
2005	3409.26	3609.66	5.88
2006	3439.45	3638.51	5.79
2007	3338.49	2140.38	−35.89
2008	3359.37	2893.90	−13.86
2009	3335.40	3546.10	6.32
2010	3309.10	2636.50	−20.33

* Rice_MOD denotes the rice planting area derived from MODIS.

Furthermore, county-level validation of the rice planting area extraction results was performed. Because of the different standards for agricultural census data collected at the county level, single cropping rice and double cropping late rice were combined and treated simply as late rice for the comparison. The comparison results for early and late rice are shown in Figure 11. The solid line in the plot is the 1:1 line. The points in the plot are clustered near the 1:1 line, indicating that the MODIS-derived area of early rice is well correlated with the agricultural census data at the county level. The decision coefficients (R²) for early rice and late rice are 0.560 and 0.619, respectively. The MODIS-derived area of late rice demonstrates a higher consistency with the census data during the 2000–2010 period, and the extracted early-rice area exhibits greater bias than that of late rice.

Figure 11. Correlation between areas of (a) early rice and (b) late rice derived using the MODIS algorithm and from the census data at the county level for the 2000–2010 period.

Figure 11. Correlation between areas of (a) early rice and (b) late rice derived using the MODIS algorithm and from the census data at the county level for the 2000–2010 period.

The topography in Zhejiang Province is very complicated, including plains, hills and mountains. Xu and Wang studied regionalization for rice yield estimation in Zhejiang Province by considering the local rice cropping systems, agroclimates, landforms, surface feature structures and rice yield levels. The county borders were treated as the region boundaries in the regionalization [36]. According to the regionalization map, the MODIS-derived early-rice area was close to that indicated by the census data in counties dominated by plains, but a large error was still observed in counties that grew less rice (Figure 12). The MODIS-derived late-rice area in counties dominated by plains was very close to that indicated by the census data. The Jinhua-Quzhou basin, located in central Zhejiang Province, is the major rice-producing region in the study area. The planting areas of early and late rice derived from the MODIS data were underestimated in the counties located in the Jinhua-Quzhou basin region (Figure 13). The results were unsatisfactory because of the influence of the terrain on the land surface reflectance. In counties located in mountainous and hilly regions, the MODIS-derived areas of early and late rice were underestimated to different extents (Figure 14). The rice planting areas derived using the MODIS algorithm were severely underestimated in counties located in mountainous and hilly regions, where the rice fields were typically fragmentary and smaller than a MODIS pixel. Because the spatial resolution of the MODIS data used in the study was 500 m, it was unfeasible to recognize a pixel with a low abundance of rice as a rice pixel. The rice fields were not successfully identified in regions with complicated topographies. However, the MODIS-derived results are still useful for developing large-scale, timely and relatively accurate spatial datasets of paddy rice fields, especially in plain regions, and for providing vital information for yield estimation, growth monitoring, water management and greenhouse gas emission estimation.

Figure 12. Comparison of MODIS-derived areas with census data in counties dominated by plains in the 2000–2010 period.

Figure 12. Comparison of MODIS-derived areas with census data in counties dominated by plains in the 2000–2010 period.

Figure 13. Comparison of MODIS-derived areas with census data in counties located in the Jinhua-Quzhou basin in the 2000–2010 period.

Figure 13. Comparison of MODIS-derived areas with census data in counties located in the Jinhua-Quzhou basin in the 2000–2010 period.

Figure 14. Comparison of MODIS-derived areas with census data in counties located in mountainous and hilly regions in the 2000–2010 period.

Figure 14. Comparison of MODIS-derived areas with census data in counties located in mountainous and hilly regions in the 2000–2010 period.

4.3. Spatial Comparison of Extracted Paddy Rice Planting Areas

In addition to the accuracy of the total estimated area, the spatial matching of the MODIS-derived results is also very important for practical applications. Table 2 summarizes the accuracy assessment of the MODIS-derived rice results. At the pixel level, the commission and omission errors were 26.30% and 22.67%, respectively. The user’s and producer’s accuracies were 73.70% and 77.33%, respectively. It is difficult to co-register Landsat ETM+ imagery with MODIS data because of the extremely large difference in spatial resolution between the two datasets. The pixels at the edges of discriminated rice pixels can give rise to considerable error when the accuracy validation is conducted for individual pixels.

Table 2. Accuracy assessment of MODIS-derived results at the studied sites.

Level	Commission Error (%)	User’s Accuracy (%)	Omission Error (%)	Producer’s Accuracy (%)
Pixel level	26.30	73.70	22.67	77.33
3 × 3 window	3.23	96.77	0.04	99.96

Table 2. Accuracy assessment of MODIS-derived results at the studied sites.

Level	Commission Error (%)	User’s Accuracy (%)	Omission Error (%)	Producer’s Accuracy (%)
Pixel level	26.30	73.70	22.67	77.33
3 × 3 window	3.23	96.77	0.04	99.96

Therefore, the commission and omission errors for moving windows of 3 × 3 pixels were also calculated. If a pixel were identified as a rice pixel in the MODIS-derived result, but the eight pixels surrounding it were all labeled as non-rice pixels in the aggregated reference rice map, we considered it to be a committed pixel. If a pixel were labeled as a rice pixel in the aggregated rice map, but no pixels in the moving 3 × 3 pixel window surrounding it were identified as rice in the MODIS-derived result, it was considered to be an omitted pixel. In this analysis, the user’s and producer’s accuracies were found to be 96.77% and 99.96%, respectively.

Figure 15. Irrigated paddies in the land use map (a) and the rice distribution map derived from MODIS (b) for Wenzhou City in 2005.

Figure 15. Irrigated paddies in the land use map (a) and the rice distribution map derived from MODIS (b) for Wenzhou City in 2005.

Although an irrigated paddy indicated in the land use map may be used not only for planting paddy rice, but also for planting aquatic plants, such as reeds and lotus roots, any MODIS-derived rice region should be located in an irrigated paddy. Therefore, a pixel was considered to be misclassified if it was a MODIS-derived rice pixel located in a non-irrigated paddy region. The irrigated paddy regions in the land use map of Wenzhou City for 2005 were extracted for the validation of the spatial matching of the results derived using the MODIS algorithm. In Figure 15a, it is seen that the majority of irrigated paddies were concentrated in the eastern coastal region, especially in Yueqing, Rui’an, Pingyang and Cangnan, and that the total area of irrigated paddies in Wenzhou City in 2005 was 109.7 kha. The two maps shown in Figure 15 were overlaid to examine the agreement between them. The number of pixels extracted using the MODIS algorithm was 4309, 88.95% of which were located in irrigated paddy regions. As indicated by the spatial matching analysis at the county level, the accuracy was lowest in counties with less than 5 kha (Table 3).

Table 3. Consistency analysis of MODIS-derived rice pixels with irrigated paddies in Wenzhou City, 2005.

City	Number of Rice Pixels Located in Irrigated Paddies	Number of RiceMOD *	Overlap Proportion (%)	Area of Irrigated Paddies (kha)	Area of RiceMOD (km2)
Yueqing	832	888	93.69	21.38	19.06
Pingyang	586	625	93.76	16.50	13.42
Wencheng	55	93	59.14	3.51	2.00
Yongjia	178	238	74.79	7.99	5.11
Taishun	5	53	9.43	0.66	1.14
Dongtou	3	52	5.77	0.06	1.12
Wenzhou	371	440	84.32	12.44	9.44
Rui’an	810	878	92.26	25.08	18.85
Cangnan	993	1042	95.30	22.10	22.37
Total	3833	4309	88.95	109.70	92.50

* Rice_MOD denotes the rice pixels derived from MODIS.

Table 3. Consistency analysis of MODIS-derived rice pixels with irrigated paddies in Wenzhou City, 2005.

City	Number of Rice Pixels Located in Irrigated Paddies	Number of RiceMOD *	Overlap Proportion (%)	Area of Irrigated Paddies (kha)	Area of RiceMOD (km2)
Yueqing	832	888	93.69	21.38	19.06
Pingyang	586	625	93.76	16.50	13.42
Wencheng	55	93	59.14	3.51	2.00
Yongjia	178	238	74.79	7.99	5.11
Taishun	5	53	9.43	0.66	1.14
Dongtou	3	52	5.77	0.06	1.12
Wenzhou	371	440	84.32	12.44	9.44
Rui’an	810	878	92.26	25.08	18.85
Cangnan	993	1042	95.30	22.10	22.37
Total	3833	4309	88.95	109.70	92.50

* Rice_MOD denotes the rice pixels derived from MODIS.

5. Conclusions

In this study, the NWDWI was proposed to enhance the signal of flooding regions in remotely-sensed images. The algorithm for the identification of flooded pixels was evaluated based on spectral data measured in the field, as well as ETM+ and MODIS data. Built-up regions and natural watersheds could be readily separated using the NWDWI and EVI. The spatial distribution maps of rice with different cropping patterns (i.e., early rice, single cropping rice and double cropping late rice) in the period of 2000–2010 were generated using a decision tree classification algorithm. The accuracy of the extracted annual total rice area was greater than 85%, except for 2007 and 2010, for which it was poorer, because of the large areas of cloud masking during the transplanting period. The identified rice areas were also validated at the county level. The MODIS-derived area of late rice demonstrated a higher consistency with the census data during the period of 2000–2010. The user’s and producer’s accuracies for moving windows of 3 × 3 pixels were both greater than 95%. The algorithm also revealed the interannual variations in single and double cropping rice in the Yangtze River Delta region.

However, there were several factors that may have affected the accuracy of the results, such as cloud contamination, spatial resolution and topography. The accuracy was not satisfactory in counties with complex terrain. The value of a in the formula for the NWDWI used in this study was determined by field experimental data and the ETM+ image and was appropriate for discriminating rice pixels at a spatial resolution of 500 m during the transplanting stage. The higher value of a in NWDWI could lead to misclassifying numbers of non-rice pixels as rice pixels. Otherwise, rice pixels could be omitted by using a lower value of a in NWDWI. The developed algorithm was found to be unsuitable for use in regions with continuous rainy weather.

Despite the uncertainties in this algorithm, MODIS data are a suitable choice for generating rice distribution maps at large scales, which are useful for long-term grain yield estimations and the detection of changes in land use/cover change. The application of Aqua/MODIS data in combination with Terra/MODIS data could improve the accuracy of our algorithm.