Mapping Rice Cropping Systems in Vietnam Using an NDVI-Based Time-Series Similarity Measurement Based on DTW Distance

1 State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China; guanxd@lreis.ac.cn (X.G.); liugh@lreis.ac.cn (G.L.); liuqs@lreis.ac.cn (Q.L.) 2 College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China 3 Key Laboratory of Ecosystem Network Observation and Modeling, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China 4 Department of Geography and Anthropology, Louisiana State University, Baton Rouge, LA 70803, USA; smeng@lsu.edu 5 Jiangsu Center for Collaborative Innovation in Geographic Information Resource Development and Application, Nanjing 210023, China * Correspondence: huangch@lreis.ac.cn; Tel.: +86-10-6488-8372; Fax: +86-10-6488-9630


Introduction
Rice is the most important crop in Vietnam and accounts for 40% of gross agricultural production in the country [1][2][3].As a result of its suitable climate, Vietnam has increased its rice production continually since 1995 (FAO, 2012), eventually becoming the second-largest rice exporting country in 2012 [4].Effective rice monitoring is essential for supporting policy making with accurate information concerning rice-growing areas, estimation of production, and seasonal cropping patterns to achieve increased economic development and food security.
Remote sensing represents a reliable and effective approach to assessing crop growth and production at various scales [5][6][7][8][9][10][11].To map the distribution of rice fields, moderate-scale analyses often use high to medium spatial resolution (6 to 30 m) optical data acquired from Landsat [12,13], Systeme Probatoire d'Observation de la Terre (SPOT) [14], and Indian Remote Sensing Satellite (IRS) [15].Regional-and global-scale analyses commonly apply coarser spatial resolution imagery with larger swath widths such as imagery from the Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer (MODIS), and SPOT VEGETATION (SPOT VGT) [16].Various classification methods, such as unsupervised classification and supervised classification, have been adopted to map rice systems [17,18].The phase analysis method using time-series VI indicates the greenness of vegetation growing in different phases.Notably, Xiao et al. [19,20], using LSWI, EVI, and NDVI information derived from MODIS data to map rice systems in China and Southeast Asia based on the fact that the land surface includes a mixture of surface water and green rice plants during the flooding and rice transplanting period, found that crops generally exhibit distinct growth characteristics during different seasons [21].Gumma et al. [22] used a spectral similarity value to evaluate the similarities between the NDVI time series and rice pixels grouped using the 500-m resolution MODIS data and mapped the Bangladesh rice cultivation distribution.Chen et al. [23] applied a linear mixture model with 250-m resolution MODIS NDVI data to classify the rice-cropping systems in the Mekong Delta of Vietnam.Kontgis et al. [24] recently mapped rice paddy extent and intensification in the Mekong River Delta with dense time stacks of Landsat data (high to medium spatial resolution data); in this project, their methods drew upon other methods developed by Xiao et al. using MODIS time-series data [19,20].
Several factors related to rice cropping and the local environment present challenges for rice mapping using low-temporal-resolution remote sensing data such as TM data and SPOT data.First, the suitable sunlight, heat, and water conditions in Southeast Asia allow for a flexible rice planting schedule.The variation in rice-growing periods produces mixed pixels in remote sensing data.Second, these satellites' long revisit cycles of 16 days to 26 days impose limitations on obtaining real-time images.Third, the cloudy and rainy weather conditions of Southeast Asia often produce heavy cloud coverage in images and significantly decrease the availability of high-quality images.According to the statistics in this study, within the Vietnamese territory, the probability of obtaining images with less than 20% cloud cover is less than 5% (according to cloud detection applied to TM images in Vietnam in 2010) and during the monsoon season of Vietnam (May to November, namely, days 120 to days 304), the mean cloud coverage was 47% (Figure 1).As a result, remote sensing classification methods based on a single temporal vegetation index are limited.However, crops of the same type have the same phenology and share certain NDVI curve characteristics.Identifying crop cycles from high-temporal-resolution images continues to represent a promising approach [25].To mine the temporal characteristics of rice growth from high-temporal-resolution data, we applied a Dynamic Time Warping (DTW) similarity measure to time-series NDVI data derived from MODIS imagery to extract rice cultivated areas.
The time-series similarity measure is widely used in many fields primarily to describe the characteristics of data variations over time.Time-series similarity research has become a core component of time-series data mining.Time-series similarity also has numerous prospects in such remote sensing applications as the extraction of vegetation areas, and such data have many applications.For example, Geerken and Zaitchik [25] classified pasture vegetation using time-series similarity based on the Fourier transform method with MODIS NDVI time series.Gumma et al. [22] mapped seasonal rice cropland using 500-m-resolution MODIS data based on a spectral matching technique.DTW (dynamic time warping) is commonly used to compare the similarity of two given time series under certain conditions.The DTW distance, which is the measure of the DTW-based similarity of two time series, was originally applied to the matching of text data in the field of speech processing; visual pattern recognition, curve similarity measurement, and time-series classification and clustering in the field of speech recognition; medical analysis; and moving object identification.Research has shown that this algorithm, which is based on nonlinear bending, can obtain high matching accuracy [26,27].The DTW distance can also be integrated into many classification techniques such as nearest neighbor classifiers, support vector machines, and neural networks to replace the Euclidean distance and thus improve classification accuracy [28][29][30][31][32][33][34].Petitjean et al. tested the K-means algorithm with DTW and DTW barycenter averaging (DBA) on satellite image time-series analysis and used a K-means-based algorithm and DTW distance to cluster pixels from a SPOT image time series [35,36].Izakian et al. evaluated Fuzzy C-means (FCM) clustering and Fuzzy C-medoids (FCMdd) clustering and improved both with a proposed hybrid technique [37].Jeong et al. proposed a new kernel function based on a weighted DTW distance for support vector machines for classifying non-aligned time-series data [38].
Remote Sens. 2016, 8, 19 3/25 classification and clustering in the field of speech recognition; medical analysis; and moving object identification.Research has shown that this algorithm, which is based on nonlinear bending, can obtain high matching accuracy [26,27].The DTW distance can also be integrated into many classification techniques such as nearest neighbor classifiers, support vector machines, and neural networks to replace the Euclidean distance and thus improve classification accuracy [28][29][30][31][32][33][34].Petitjean et al. tested the K-means algorithm with DTW and DTW barycenter averaging (DBA) on satellite image time-series analysis and used a K-means-based algorithm and DTW distance to cluster pixels from a SPOT image time series [35,36].Izakian et al. evaluated Fuzzy C-means (FCM) clustering and Fuzzy C-medoids (FCMdd) clustering and improved both with a proposed hybrid technique [37].Jeong et al. proposed a new kernel function based on a weighted DTW distance for support vector machines for classifying non-aligned time-series data [38].The main objective of this paper is to map rice cropping systems in Vietnam for the year 2010 according to the similarity between the NDVI time series of every pixel and the standard rice growth NDVI time series.We use the DTW distance, which finds the optimal alignment (or coupling) between two sequences of numerical values and captures flexible similarities by aligning the coordinates inside both sequences, to address irregular rice growth NDVI time-series data (caused by the flexible rice planting schedule in Vietnam).We verified the performance of the DTW-distance-based similarity method applied to mapping rice-growing systems in Vietnam, and the results demonstrated that this method can be applied to the determination of rice cropping systems over wide areas in which rice planting schedules are flexible.
The study is structured as follows.Section 2 introduces the study area and utilized dataset and discusses why we chose the MOD09A1 data set and NDVI time series.Section 3 introduces three filtering algorithms in the TIMESAT software and shows how we compare them.Section 4 discusses how we constructed the standard rice growth NDVI time-series and the DTW theory and how we applied the DTW similarity measure to the NDVI time-series data to extract rice cultivation areas.Section 5 presents the results.We map the distribution of single rice cropping, double rice cropping, and triple rice cropping in Vietnam using a threshold in the DTW distance map.Then we use the statistical data combined with the sampled data from field surveys and Google Earth The main objective of this paper is to map rice cropping systems in Vietnam for the year 2010 according to the similarity between the NDVI time series of every pixel and the standard rice growth NDVI time series.We use the DTW distance, which finds the optimal alignment (or coupling) between two sequences of numerical values and captures flexible similarities by aligning the coordinates inside both sequences, to address irregular rice growth NDVI time-series data (caused by the flexible rice planting schedule in Vietnam).We verified the performance of the DTW-distance-based similarity method applied to mapping rice-growing systems in Vietnam, and the results demonstrated that this method can be applied to the determination of rice cropping systems over wide areas in which rice planting schedules are flexible.
The study is structured as follows.Section 2 introduces the study area and utilized dataset and discusses why we chose the MOD09A1 data set and NDVI time series.Section 3 introduces three filtering algorithms in the TIMESAT software and shows how we compare them.Section 4 discusses how we constructed the standard rice growth NDVI time-series and the DTW theory and how we applied the DTW similarity measure to the NDVI time-series data to extract rice cultivation areas.Section 5 presents the results.We map the distribution of single rice cropping, double rice cropping, and triple rice cropping in Vietnam using a threshold in the DTW distance map.Then we use the statistical data combined with the sampled data from field surveys and Google Earth sampling data to analyze the accuracy of the rice distribution map.Section 6 includes some conclusions and discusses future work.

Brief Introduction to the Study Area
The study area of this research includes the entire country of Vietnam.Vietnam is located in the eastern region of the Indo-China peninsula, extending from 8 ˝10 1 to 23 ˝24 1 N and from 102 ˝09 1 to 109 ˝30 1 E and covering an area of 329,556 km 2 .Vietnam is long and narrow, with high elevations in the west and low elevations in the east.Vietnam experiences a tropical monsoon climate characterized by high temperatures and rainy seasons.The annual average temperature is approximately 24 ˝C, and the average annual rainfall ranges from 1500 mm to 2000 mm.
Vietnam is a traditional agricultural country; the agricultural population accounts for approximately 75 percent of the total population.Cultivated land and forest land constitute 60% of the total area.Food crops include rice, corn, potatoes, sweet potatoes, and cassava, and economic crops include coffee, rubber, cashew nuts, peanuts, silk, and tea.
According to the topography, soil, and climate conditions, agricultural activity can be divided into eight agricultural ecological sections: North Central Coast, North East, Mekong River Delta, Red River Delta, South East, South Central Coast, North West, and Central Highlands (Figure 2).In Vietnam, much of the rice cultivation is concentrated in two river deltas: the Mekong (over half of the country's rice area) and Red deltas.The planted area in 2000 was approximately 7.6 million ha, and the cropping intensity (ratio of sown area to land area for a given crop) is the highest in the world [39].The high cropping intensity is largely due to the planting of the triple rice crops that are common in much of the Mekong Delta area [40].According to Son et al. [41], five cropping seasons per year are observed in Vietnam: the rainy season (July-August to December-January), the winter-spring season (November-December to February-March), the spring-summer season (March-April to May-June), the summer-autumn season (April-May to July-August), and the autumn-winter season (July-September to October-December).The rain-fed rice cropping system using long-duration varieties (160-180 days) is invariably applied in areas that are subject to soil and water constraints, whereas irrigated rice cropping systems using short-duration varieties (90-100 days) are applied in areas where irrigation conditions are favorable.

MODIS NDVI Time-Series Data
MODIS imagery was acquired from NASA's (National Aeronautics and Space Administration) website [42]; specifically, the surface reflectance MOD09A1 eight-day composite products named "MODIS/Terra surface reflectance 8-day Global 500 m SIN GRID V005" for 2010 were used.We obtained the first two bands of the MOD09A1 data for NDVI calculation as well as the corresponding Quality Assessment (QA) data.Five tiles (i.e., h27v06, h27v07, h28v06, h28v07, and h28v08) covered the entire territory of Vietnam.The MOD09A1 data products were systematically

MODIS NDVI Time-Series Data
MODIS imagery was acquired from NASA's (National Aeronautics and Space Administration) website [42]; specifically, the surface reflectance MOD09A1 eight-day composite products named "MODIS/Terra surface reflectance 8-day Global 500 m SIN GRID V005" for 2010 were used.We obtained the first two bands of the MOD09A1 data for NDVI calculation as well as the corresponding Quality Assessment (QA) data.Five tiles (i.e., h27v06, h27v07, h28v06, h28v07, and h28v08) covered the entire territory of Vietnam.The MOD09A1 data products were systematically corrected for the effects of gaseous and aerosol scattering.Each pixel contains the optimal quality L2G observation data for eight days with information about the highest observation range and minimum observation angle (without clouds or cloud shadows) [43].
The MOD13 datasets provided the level 2G (gridded) surface reflectance and temporally composited these values to generate the 16-day, 250/500-m, or 1-km VI products.However, certain key periods in the rice-growing process (e.g., transplanting) last for only a few days, and thus, the 16-day composite VI data might be insufficient for catching critical growing phases.Therefore, the NDVI time-series similarity measure was applied to the MOD09A1 data.
As an effective measurement of vegetation with high temporal resolution from MODIS, NDVI data facilitate meaningful comparisons of seasonal and inter-annual changes in vegetation growth and activity.The strength of NDVI lies in its rationing concept, which reduces many forms of noise (illumination differences, cloud shadows, atmospheric attenuation, and certain topographic variations) present in multiple bands.In addition, a comparison of sensitivity to the dynamics of vegetative cover shows that NDVI data are more sensitive because they enable more accurate change detection [44] as a result of NDVI data presenting higher variability and abrupt changes in the growing season [45].Reflectance values in band 1 and band 2 in the MOD09A1 product are used to calculate NDVI, as shown in Equation (1): where ρ nir is the reflectance value of the NIR red band of the MOD09A1 data and ρ red is the reflectance value of the red band of the MOD09A1 data.

Digital Elevation Model (DEM) Data
DEM data were used to mask the territory where rice planting was unlikely to occur.We chose the Shuttle Radar Topography Mission (SRTM) 90-m DEM version 4 data products.The SRTM was provided by The CGIAR Consortium for Spatial Information (CGIAR-CSI) [46].Eight images covering Vietnam (lines 57 to 58 and columns 8 to 11) were downloaded.

Ancillary Data
Ground-truth data were collected for class identification and result verification, and we performed field investigations at 229 locations from 17 March to 3 April 2013 (Figure 2).Among the 229 locations, 98 locations are long-term cultivated rice fields.We obtained information on the irrigation system used for rice cropping by speaking with local residents at 53 rice growth locations.In addition, we used Google Earth to select large rice fields in the high-spatial-resolution image and obtained a total of 1468 locations, for which imagery was obtained during October 2009-March 2010.Among the 1468 Google Earth selected points, 268 pure rice field points were used to build a standard rice growth NDVI time-series base, combined with the 98 field survey rice field points.Another 1200 points were reserved for validation purposes.
In addition, statistical data from the General Statistics Office of VIETNAM [47] provided data on the planted areas of spring rice, autumn rice, and winter rice by province for 2010.

Data Preprocessing
Even after being corrected, the effects of the BRDF (Bidirectional Reflectance Distribution Function), clouds, ozone, dust, and other aerosols generally decreased the near infrared reflectance and affected the VI time series in the form of noise.In particular, due to the cloudy and rainy climate characteristic sf Vietnam, cloud cover seriously affects the quality of the NDVI time series, especially during the monsoon season.We used the internal cloud algorithm flag to calculate the possibility of cloud appearance in Vietnam [48].As shown in Figure 1, among the 46 MODIS images, the cloud coverage of 33 MODIS images exceeded 20%, and that of 16 images exceeded 40%.During the monsoon season of Vietnam (May to November, namely, days 120 to days 304), the mean cloud coverage was 47%.Thus, we first de-noised the NDVI time series to reduce the influence of the atmosphere.
Various noise reduction methods have been adapted to temporal NDVI data reconstruction [49][50][51][52][53]. Hird and McDermid [54] advised caution in the application of noise reduction techniques and also argued that the strength and nature of the noise present in a dataset should be considered when selecting a method for NDVI time-series reconstruction.We manually created some NDVI time series with cloud gaps from a clear sky rice growth NDVI time series and using three noise reduction methods (the adaptive Savizky-Golay filter, the Asymmetric Gaussian function, and the Double logistic function) in the TIMESAT software [55] to reconstruct the gap NDVI time series with consideration of the QA data, we found that the Logistic algorithm's overall fitting degree was the highest but that the S-G algorithm provided the best fit at the turning points caused by changes in rice phenology.However, it is difficult to produce a quantitative evaluation of the best algorithm for fitting every curve in every phase.Thus, we assessed the characteristics of each algorithm, as shown in Figure 3.The curves fitted by the Logistic and Gaussian algorithms are smooth, and the Logistic algorithm is the "smoothest" curve-fitting method.However, the curves fit by the S-G algorithm preserve more of the original NDVI values at the turning points.Finally, we reconstructed the NDVI time-series data using S-G filtering combined with the QA data.Adaptive Savizky-Golay filtering was used to replace each data value yi, i = 1,…, N with a linear combination of nearby values in a window, as shown in Equation ( 2): (2) Finally, we reconstructed the NDVI time-series data using S-G filtering combined with the QA data.Adaptive Savizky-Golay filtering was used to replace each data value y i , i = 1, . . ., N with a linear combination of nearby values in a window, as shown in Equation (2): where c j is the weight of the ith NDVI value of the filter, which is replaced by the weighted QA value of the NDVI.The data value y i is replaced by the average of the values in the window.The window width n denotes the filter size and is the half-width of the smoothing window.The index j is the running index of the original ordinate data table.
During the fitting process, universal settings were applied to each algorithm, namely, the adaptation strength was equal to 4.0, the number of envelope iterations was equal to 2, and the same QA data layers were used.In addition, in the case of Savizky-Golay filtering, the size of the moving window was set to 4 [56].

Methods
After the NDVI time series were reconstructed via S-G filtering, we applied DTW-based NDVI time-series similarity measurements for detecting rice from time-series MODIS data.The method consists of two major procedures: (i) building a standard NDVI time-series base for rice growth cycles through field sampling data; and (ii) extracting rice fields based on the DTW distance of the standard rice-growth NDVI time-series data.The DTW distance was applied to sequences of NDVI values of every pixel.Each pixel had two attributes, namely, its coordinates (x, y) and an NDVI sequence defined as < NDVI 1 , NDVI 2 , . . ., NDVI n >.Then, a sequence was identified as <NDVI 1 (x, y), NDVI 2 (x, y), . . ., NDVI n (x, y)>.Because the MOD09A1 annual period includes 46 images, n equals 46.We then calculated the DTW distance of the sequence with the standard rice growth NDVI time series of each region, and a DTW distance was calculated for each pixel with the standard rice growth NDVI time series, thereby determining the similarity of each pixel's NDVI time series with the standard rice growth NDVI time series.Using this similarity measure, we distinguished the rice from other crops.

Building Standard NDVI Time-Series Base
After the NDVI data were subjected to noise filtering, we combined the 98 long-term-cultivated pure rice field points obtained from field survey data with the 268 pure rice field points obtained from Google-Earth-selected rice field points (a total of 366 points) and composed a standard rice growth NDVI time-series base as a reference NDVI time series.Among the 366 sampled rice field points, 92 were located in the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.We averaged the NDVI time series with similar shapes and obtained a typical NDVI curve for each region.In addition, we used the same standard rice growth NDVI time series for the South Central Coast, Central Highlands, and South East regions because of the similar terrain, light, and temperature conditions characterizing these regions.Because the field survey data included 53 points containing information about the irrigation system of the rice growth points, we indicated the single, double, and triple rice by the shape of the NDVI time series of these 53 rice growth points.The Google-Earth-selected rice field points do not include information about the irrigation system in Vietnam, and we judged the rice cropping situation based on the shapes of curves, similar to the 53 NDVI time-series points from the field-collected information.

Dynamic Time Warping Based on Time-Series Similarity Measurements
DTW uses a dynamic programming technique to find the minimal distance between two time series, where sequences are warped by stretching or shrinking the time dimension [57].The DTW algorithm is described below.
Suppose that we have two time series, i.e., C = {c 1 ,c 2 , . . .,c m } and Q = {q 1 ,q 2 , . . .,q n }, with respective lengths of m and n.We construct an m ˆn matrix A mˆn and define the distance between each element as a ij = d(C i ,Q j ) = c ´Ci ´Qj ¯2.In the matrix A mˆn , a winding path is set by a group of adjacent matrix elements, and notes for W = {w 1 , w 2 ,..., w k } and the kth element in W are defined as w k = (a ij ) k ; this path must satisfy the following conditions: This element must satisfy the condition 0 ď i ´i', 0 ď j ´j' ď 1, and thus, DTW (C,Q) = min The DTW algorithm can be summarized by applying ideal dynamic programming to find the best (i.e., least bending) cost path, as shown in Equation ( 3):  In this paper, C refers to NDVI time series in the standard profile (constructed NDVI time-series base), and Q refers to the NDVI time series of the pixel.ci and qj denote the absolute value of the NDVI value in the NDVI time series.In Vietnam, the crop phenological stages are not fixed, and the planting dates depend on rainfall conditions.The planting date for irrigated rice exhibits greater variance, and geographical factors extend the variance of the length of each phenological stage.Thus, the absolute value between the NDVI time series of rice growth pixels may differ from the reference NDVI time series.As Figure 5 illustrates, suppose that NDVI1 is the reference NDVI time series and that NDVI2 is the NDVI time series of a pixel.Because the rice planting data for the pixels In this paper, C refers to NDVI time series in the standard profile (constructed NDVI time-series base), and Q refers to the NDVI time series of the pixel.c i and q j denote the absolute value of the NDVI value in the NDVI time series.In Vietnam, the crop phenological stages are not fixed, and the planting dates depend on rainfall conditions.The planting date for irrigated rice exhibits greater variance, and geographical factors extend the variance of the length of each phenological stage.Thus, the absolute value between the NDVI time series of rice growth pixels may differ from the reference NDVI time series.As Figure 5 illustrates, suppose that NDVI1 is the reference NDVI time series and that NDVI2 is the NDVI time series of a pixel.Because the rice planting data for the pixels correspond to dates after those of the rice planting data of the reference NDVI time series, the phases are not aligned with the reference data; however, the DTW algorithm aligns the phases between these two NDVI time series by stretching or shrinking the time dimension.As a result, the DTW distance similarity measure allows for the identification of similar rice types despite dislocations of the NDVI time series because the DTW distance is small, and when applying a proper threshold to the DTW distance, they are the same type.
In this paper, C refers to NDVI time series in the standard profile (constructed NDVI time-series base), and Q refers to the NDVI time series of the pixel.ci and qj denote the absolute value of the NDVI value in the NDVI time series.In Vietnam, the crop phenological stages are not fixed, and the planting dates depend on rainfall conditions.The planting date for irrigated rice exhibits greater variance, and geographical factors extend the variance of the length of each phenological stage.Thus, the absolute value between the NDVI time series of rice growth pixels may differ from the reference NDVI time series.As Figure 5 illustrates, suppose that NDVI1 is the reference NDVI time series and that NDVI2 is the NDVI time series of a pixel.Because the rice planting data for the pixels correspond to dates after those of the rice planting data of the reference NDVI time series, the phases are not aligned with the reference data; however, the DTW algorithm aligns the phases between these two NDVI time series by stretching or shrinking the time dimension.As a result, the DTW distance similarity measure allows for the identification of similar rice types despite dislocations of the NDVI time series because the DTW distance is small, and when applying a proper threshold to the DTW distance, they are the same type.

Standard NDVI Time-Series Base and DTW Distance
For the purpose of getting standard rice growth NDVI time-series, we acquired 366 rice field points and composed a standard rice growth NDVI time-series base as the reference NDVI time series.Among the 366 sampled rice field points, a total of 92 sampled rice field points were located in the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table 1.After obtaining the standard NDVI time series, we calculated the DTW distance between the NDVI time-series of every pixel to the standard NDVI cycle at different regions, as shown in Figure 6.It is acknowledged in the DTW algorithm that the NDVI cycle shapes exhibit greater similarity if the DTW distances are shorter.To keep this paper reasonably concise, we demonstrate six DTW distances of NDVI time series for the Mekong River Delta.The six types of rice growth NDVI time series for the Mekong River Delta in Table 1 correspond to the six maps of DTW distances shown in Figure 6.The calculation is performed in MATLAB using a loop statement.
Figure 6a shows that the DTW distances for the irrigated single rice NDVI time series are large in the wild region in the southern region of the Mekong River Delta, which denotes the dissimilarity between the NDVI time series of the pixels with the standard NDVI time series of the irrigated single rice.At the river mouth and southern region, the similarity between pixels with irrigated single rice is relatively high because the DTW distance is small.Figure 6b,c show the DTW distance with the irrigated double rice cropping I.Note that the similarity is high in certain regions in the middle north region; however, in Figure 6f (DTW distance with triple rice cropping III), the DTW distances are smaller, and thus the large regions in this area are growing triple rice.

Irrigated Double Rice
Cropping II After obtaining the standard NDVI time series, we calculated the DTW distance between the NDVI time-series of every pixel to the standard NDVI cycle at different regions, as shown in Figure 6.It is acknowledged in the DTW algorithm that the NDVI cycle shapes exhibit greater similarity if the DTW distances are shorter.To keep this paper reasonably concise, we demonstrate six DTW distances of NDVI time series for the Mekong River Delta.The six types of rice growth NDVI time series for the Mekong River Delta in Table 1 correspond to the six maps of DTW distances shown in Figure 6.The calculation is performed in MATLAB using a loop statement.
Figure 6a shows that the DTW distances for the irrigated single rice NDVI time series are large in the wild region in the southern region of the Mekong River Delta, which denotes the dissimilarity between the NDVI time series of the pixels with the standard NDVI time series of the irrigated single rice.At the river mouth and southern region, the similarity between pixels with irrigated single rice is relatively high because the DTW distance is small.Figure 6b,c show the DTW distance with the irrigated double rice cropping I.Note that the similarity is high in certain regions in the middle north region; however, in Figure 6f (DTW distance with triple rice cropping III), the DTW distances are smaller, and thus the large regions in this area are growing triple rice.

DTW Threshold and Rice Distribution Map
The DTW distance can reflect the degree of similarity as well as dissimilarity of the standard rice growth NDVI time-series with the pixels' NDVI time-series.It is acknowledged in the DTW algorithm that the NDVI cycle shapes exhibit greater similarity if the DTW distances are shorter.We perform an interactive examination of the NDVI time series and ground-truth data to determine the cut-off DTW distance thresholds to identify single-crop rice and multiple-crop rice.The threshold values were set as shown in Table 2.It is assumed that the pixels' NDVI time-series' DTW distance greater than the threshold shows in Table 2 are not likely to be a rice planted field.

DTW Threshold and Rice Distribution Map
The DTW distance can reflect the degree of similarity as well as dissimilarity of the standard rice growth NDVI time-series with the pixels' NDVI time-series.It is acknowledged in the DTW algorithm that the NDVI cycle shapes exhibit greater similarity if the DTW distances are shorter.We perform an interactive examination of the NDVI time series and ground-truth data to determine the cut-off DTW distance thresholds to identify single-crop rice and multiple-crop rice.The threshold values were set as shown in Table 2.It is assumed that the pixels' NDVI time-series' DTW distance greater than the threshold shows in Table 2 are not likely to be a rice planted field.Each pixel's DTW distance to the standard rice-growing NDVI curve was calculated, and pixels with a DTW distance for a crop type that exceeded these thresholds were considered unlikely to be the dominant crop type across the extent of the pixel.Overlaid pixels with smaller DTW distances within the threshold value were caused by the presence of mixed pixels (for example, if the DTW distance of triple rice cropping III in Mekong River Delta for a pixel is 2.7 and if the DTW distance of triple rice cropping I in Mekong River Delta is 2.8 for the same pixel, where both values are below the threshold), and we determined the crop type based on the lowest DTW distance, namely triple rice cropping III in the Mekong River Delta.Finally, the spatial distribution of rice in Vietnam is shown in Figure 7.
The application of the DTW-based similarity measure to the NDVI time series for the Vietnam rice cropping system for 2010 produces a rice crop type map for Vietnam as shown in Figure 7.The Red River Delta and Mekong River Delta are the most intensive rice cropping regions in Vietnam, and the respective cropping intensities are approximately 105.5% and 101.6%.Because the statistical records of the rice cropping areas by province in 2010 do not discriminate plantation times during the year, we unified the statistical data with the MODIS-derived data by summing up the extracted single rice areas; if there were double rice cropping areas, we counted the area twice and added it to the single rice area, so as to obtain the triple rice cropping area.For example, in the Red River Delta, the total crop area of double rice is 7704.5 thousand ha, single rice cropping in the Red River Delta represents an area of 165 thousand ha, and the total sown area is 7704.5 ˆ2 + 165 = 15,574 thousand ha.The Mekong River Delta exhibits mixed types of single, double, and triple cropped rice, primarily being composed of double and triple cropped rice.The South East Coast contains single, double, and triple cropped rice.Other regions in Vietnam, including the North East, North West, North East Coast, Central Highland, and Southeast regions, are mountainous and contain mostly single cropped rice in the costal lowlands.The areas of single rice, double rice, and triple rice's area derived by MODIS data and the statistical data are shown in Table 3.

Accuracy Assessment
Figure 8 shows a comparison between the statistical data and the MODIS-derived rice areas.The statistical records of the total rice cropped areas in 2010 by province are presented.This comparison shows that at the national scale, as demonstrated in Figure 8a, the R 2 value is 0.809, indicating that the two data points fit well.However, at the province scale, the results are not ideal.We separate the results into three components using the statistical data.The first component is rice-planted areas smaller than 500 km 2 in the statistical data, representing 21 provinces.The second component includes rice-planted areas from 500 km 2 to 1000 km 2 in the statistical data, representing 18 provinces, and the third component is the statistical rice-planted areas from 1000 km 2 to 7000 km 2 , representing 22 provinces.The comparison between the statistical data and the MODIS-derived rice areas for 21 provinces of the first component, where rice planting might be scattered or small, is shown in Figure 8b; the R 2 between the two datasets in this component is 0.086, which shows that they are significantly unrelated.The comparison between the statistical data and the MODIS-derived rice areas for 18 provinces of the second component is shown in Figure 8c; the R 2 of the comparison for this component is 0.226, which is greater than that in the first comparison; however, the result remains unsatisfactory.For the third component for 22 provinces where rice planting is relatively intensive, the comparison between the statistical data and the MODIS-derived rice areas is shown in Figure 8d, and the R 2 between the two datasets was 0.661.

Accuracy Assessment
Figure 8 shows a comparison between the statistical data and the MODIS-derived rice areas.The statistical records of the total rice cropped areas in 2010 by province are presented.This comparison shows that at the national scale, as demonstrated in Figure 8a, the R 2 value is 0.809, indicating that the two data points fit well.However, at the province scale, the results are not ideal.We separate the results into three components using the statistical data.The first component is rice-planted areas smaller than 500 km 2 in the statistical data, representing 21 provinces.The second component includes rice-planted areas from 500 km 2 to 1000 km 2 in the statistical data, representing 18 provinces, and the third component is the statistical rice-planted areas from 1000 km 2 to 7000 km 2 , representing 22 provinces.The comparison between the statistical data and the MODIS-derived rice areas for 21 provinces of the first component, where rice planting might be scattered or small, is shown in Figure 8b; the R 2 between the two datasets in this component is 0.086, which shows that they are significantly unrelated.The comparison between the statistical data and the MODIS-derived rice areas for 18 provinces of the second component is shown in Figure 8c; the R 2 of the comparison for this component is 0.226, which is greater than that in the first comparison; however, the result remains unsatisfactory.For the third component for 22 provinces where rice planting is relatively intensive, the comparison between the statistical data and the MODIS-derived rice areas is shown in Figure 8d, and the R 2 between the two datasets was 0.661.Due to the notably low extraction accuracy of the rice fields in the provinces, we examined whether the MOD09A1 data have sufficient spatial resolution to discriminate rice fields in these provinces.Figure 9 shows a fractal rice field in the North East in which rice was planted near the river.In Figure 9a, the rice cropping is scattered due to the mountainous terrain.In Figure 9c, the DTW distances of the three rice fields with standard rice field NDVI are 0.42, 0.53, and 0.57, exceeding the threshold in Table 2 so that it would not be recognized as a rice field.Due to the notably low extraction accuracy of the rice fields in the provinces, we examined whether the MOD09A1 data have sufficient spatial resolution to discriminate rice fields in these provinces.Figure 9 shows a fractal rice field in the North East in which rice was planted near the river.In Figure 9a, the rice cropping is scattered due to the mountainous terrain.In Figure 9c, the DTW distances of the three rice fields with standard rice field NDVI are 0.42, 0.53, and 0.57, exceeding the threshold in Table 2 so that it would not be recognized as a rice field.The parcel area of paddy fields is related to classification accuracy, and to evaluate the heterogeneity of land cover we performed multi-resolution segmentation (the scales of the multi-resolution segmentation are the same) inside MODIS pixels on the TM image, and in each region, one TM image in the DEM masked the map to exclude water and permanent woodland.The relationship between the land cover heterogeneity and DTW-extracted rice field accuracy is shown in Figure 10; the figure shows the decreasing trends of extraction accuracy when the number of patches per MODIS pixel is increased, which denotes a greater heterogeneity of land cover.The extraction accuracy is computed as (Ae-As)/As, where Ae is the extracted rice area and as is the statistical rice area.The number of patches per MODIS pixel is computed as the average number of multi-resolution segmented patches of TM images in MODIS pixels.This result shows that the rice extraction accuracy of a small rice-cropping region is decreased due to the coarse spatial resolution of MOD09A1 data.

Conclusions
We have developed a new systematic method for detecting the area of a paddy rice system from time-series MODIS data.The method consists of three procedures: (i) filtering the time-series NDVI data using the S-G filtering; (ii) building a standard rice-growing NDVI time series base through Google earth sample points and ground-truth data; and (iii) calculating the DTW distance of the The parcel area of paddy fields is related to classification accuracy, and to evaluate the heterogeneity of land cover we performed multi-resolution segmentation (the scales of the multi-resolution segmentation are the same) inside MODIS pixels on the TM image, and in each region, one TM image in the DEM masked the map to exclude water and permanent woodland.The relationship between the land cover heterogeneity and DTW-extracted rice field accuracy is shown in Figure 10; the figure shows the decreasing trends of extraction accuracy when the number of patches per MODIS pixel is increased, which denotes a greater heterogeneity of land cover.The extraction accuracy is computed as (Ae-As)/As, where Ae is the extracted rice area and as is the statistical rice area.The number of patches per MODIS pixel is computed as the average number of multi-resolution segmented patches of TM images in MODIS pixels.This result shows that the rice extraction accuracy of a small rice-cropping region is decreased due to the coarse spatial resolution of MOD09A1 data.

Conclusions
We have developed a new systematic method for detecting the area of a paddy rice system from time-series MODIS data.The method consists of three procedures: (i) filtering the time-series NDVI data using the S-G filtering; (ii) building a standard rice-growing NDVI time series base through Google earth sample points and ground-truth data; and (iii) calculating the DTW distance of the The extraction accuracy is computed as (Ae-As)/As, where Ae is the extracted rice area and as is the statistical rice area.The number of patches per MODIS pixel is computed as the average number of multi-resolution segmented patches of TM images in MODIS pixels.This result shows that the rice extraction accuracy of a small rice-cropping region is decreased due to the coarse spatial resolution of MOD09A1 data.

Conclusions
We have developed a new systematic method for detecting the area of a paddy rice system from time-series MODIS data.The method consists of three procedures: (i) filtering the time-series NDVI data using the S-G filtering; (ii) building a standard rice-growing NDVI time series base through Google earth sample points and ground-truth data; and (iii) calculating the DTW distance of the standard rice-growth NDVI time series of each pixel and using a threshold to extract the rice field area.Our method is able to discriminate between single-rice-crop and multi-rice-crop systems in Vietnam.
In the field of pattern recognition, many algorithms have been developed to compare the similarities between curves.The DTW is a classic curve-matching technique and is defined between sequences of points.We applied the classic DTW to the recognition of NDVI curve similarity for rice field recognition.Additionally, we set a simple threshold to decide whether a pixel belongs to a specific type of rice growth NDVI time series.The selection of the threshold is performed via interactive examination with the ground-truth data and Google Earth sampling data.The accuracy assessment based on statistical data at the province level and ground-truth data obtained from field surveys and Google Earth sampling shows that our method is able to discriminate between single-rice-crop and multi-rice-crop systems in Vietnam at a national scale and that the classification accuracy is satisfactory; however, at the province scale, the accuracy is decreased.Through analysis of wrongly classified rice field points, we proved that the MOD09A1 spatial resolutions were too coarse for rice field extraction in the areas in which small rice fields are predominant and in which rice planting is scattered and fragmented.
This study confirms that the DTW approach can be applied to the determination of rice cropping systems over wide connecting areas in which rice planting scheduling is flexible and rice growth periods in different phases are uncertain.We believe that this work opens the door to a number of research directions.First, one can combine the DTW with many other mature classification approaches such as nearest neighbor and neural networks to classify satellite time-series data.Second, the DTW approach provides a more rational way to calculate the similarity or dissimilarity between curves, thereby providing a method of quantifying variations in crop cultivation.However, more scientific approaches to determining thresholds should be studied.Another problem is that the DTW distance of the standard NDVI curve shapes, which exhibit smoother and lesser dynamic changes in rice growth phenology, tend to have smaller DTW distances.As a result, the discrimination between the three types of rice is vague.Therefore, normalization between the standard NDVI curves should be studied in greater detail.
At the moment, various limitations remain, including the mixed-pixel effect.As an objective and repeatable methodology, this method is sufficiently robust for addressing complex forms of rice cropping systems in Southeast Asia.

Figure 1 .
Figure 1.Cloud coverage in images.Among the 46 MODIS images, the cloud coverage of 33 MODIS images exceeds 20%, and that of 16 images exceeds 40%.During the rainy season of Vietnam (May to November, namely, day 120 to day 304), the mean cloud coverage was 47%.

Figure 1 .
Figure 1.Cloud coverage in images.Among the 46 MODIS images, the cloud coverage of 33 MODIS images exceeds 20%, and that of 16 images exceeds 40%.During the rainy season of Vietnam (May to November, namely, day 120 to day 304), the mean cloud coverage was 47%.

Figure 3 .
Figure 3.Comparison between S-G algorithm filtering, Logistic algorithm filtering, and Gaussian algorithm filtering, where the black line denotes the original NDVI curves, which are sampled from a clear sky rice growth point; the red line denotes the S-G-algorithm-fitted NDVI curve with our manually generated gap NDVI curve; the Logistic fitted curve is denoted by the blue line; and the Gaussian-algorithm-fitted simulated gap NDVI is denoted by the green line.

Figure 3 .
Figure 3.Comparison between S-G algorithm filtering, Logistic algorithm filtering, and Gaussian algorithm filtering, where the black line denotes the original NDVI curves, which are sampled from a clear sky rice growth point; the red line denotes the S-G-algorithm-fitted NDVI curve with our manually generated gap NDVI curve; the Logistic fitted curve is denoted by the blue line; and the Gaussian-algorithm-fitted simulated gap NDVI is denoted by the green line.

Figure 5 .
Figure 5. Illustration of DTW-corrected distances.The NDVI2 cycle's planting date is later than that of the NDVI1 cycle.

Figure 5 .
Figure 5. Illustration of DTW-corrected distances.The NDVI2 cycle's planting date is later than that of the NDVI1 cycle.
standard NDVI time series, we calculated the DTW distance between the NDVI time-series of every pixel to the standard NDVI cycle at different regions, as shown in Figure standard NDVI time series, we calculated the DTW distance between the NDVI time-series of every pixel to the standard NDVI cycle at different regions, as shown in Figure standard NDVI time series, we calculated the DTW distance between the NDVI time-series of every pixel to the standard NDVI cycle at different regions, as shown in Figure standard NDVI time series, we calculated the DTW distance between the NDVI time-series of every pixel to the standard NDVI cycle at different regions, as shown in Figure standard NDVI time series, we calculated the DTW distance between the NDVI time-series of every pixel to the standard NDVI cycle at different regions, as shown in Figure standard NDVI time series, we calculated the DTW distance between the NDVI time-series of every pixel to the standard NDVI cycle at different regions, as shown in Figure

Figure 6 .
Figure 6.The results for the DTW distance: (a) The DTW distance of single rice cropping in the Mekong River Delta; (b) the DTW distance of double rice cropping I in the Mekong River Delta; (c) the DTW distance of double rice cropping II in the Mekong River Delta; (d) the DTW distance of triple rice cropping I in the Mekong River Delta; (e) the DTW distance of triple rice cropping II in the Mekong River Delta; (f) the DTW distance of triple rice cropping III in the Mekong River Delta.

Figure 7 .
Figure 7. Result of rice cultivation area extraction using DTW-based similarity measure.

Figure 7 .
Figure 7. Result of rice cultivation area extraction using DTW-based similarity measure.

Figure 10 .
Figure 10.Decreasing trends of extraction accuracy when the patches per MODIS pixel are increased, which denotes a greater heterogeneity of the land cover.

Figure 9 .Figure 9 .
Figure 9. Fractal rice field in the North East region, in which rice was planted nearby the river: (a) the Landsat TM image for 5 April 2009, in the North East region.The grid in black is the corresponding extent of a pixel in MODIS pixel.;(b) The Google Earth image obtained on 18 December 2013; (c) the NDVI curve of three points in the rice area; (d) the segmentation result based on multiresolution segmentation and the TM image.

Figure 10 .
Figure 10.Decreasing trends of extraction accuracy when the patches per MODIS pixel are increased, which denotes a greater heterogeneity of the land cover.

Figure 10 .
Figure 10.Decreasing trends of extraction accuracy when the patches per MODIS pixel are increased, which denotes a greater heterogeneity of the land cover.

Table 1 .
Standard rice growing NDVI time series and calculated DTW distance.Among the 366 sampled rice field points, a total of 92 sampled rice field points were located in the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table1.

Table 1 .
Standard rice growing NDVI time series and calculated DTW distance.

Table 1 .
Among the 366 sampled rice field points, a total of 92 sampled rice field points were located in the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table1.Standard rice growing NDVI time series and calculated DTW distance.

Table 1 .
Among the 366 sampled rice field points, a total of 92 sampled rice field points were located in the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table1.Standard rice growing NDVI time series and calculated DTW distance.

Table 1 .
Among the 366 sampled rice field points, a total of 92 sampled rice field points were located in the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table1.Standard rice growing NDVI time series and calculated DTW distance.

Table 1 .
series.Among the 366 sampled rice field points, a total of 92 sampled rice field points were located in the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table1.Standard rice growing NDVI time series and calculated DTW distance.

Table 1 .
Irrigated Double Rice Cropping series.Among the 366 sampled rice field points, a total of 92 sampled rice field points were located in the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table 1.Standard rice growing NDVI time series and calculated DTW distance.

Table 1 .
Among the 366 sampled rice field points, a total of 92 sampled rice field points were located in the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table1.Standard rice growing NDVI time series and calculated DTW distance.

Table 1 .
North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table 1.Standard rice growing NDVI time series and calculated DTW distance. the

Table 1 .
Irrigated Double Rice Cropping the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table 1.Standard rice growing NDVI time series and calculated DTW distance.

Table 1 .
Irrigated Triple Rice Cropping the North Central Coast, 86 in the Mekong River Delta, 67 in the South Central Coast, 41 in the Red River delta, 33 in the Central Highlands, 27 in the South East, 21 in the North East, and 20 in the North West.NDVI time-series with similar shapes were combined and averaged (average NDVI at each layer to get a new NDVI time-series) to obtain the typical NDVI curves for each region.The standard rice growth NDVI time-series of each region is shown in Table 1.Standard rice growing NDVI time series and calculated DTW distance.

Table 2 .
Thresholds of DTW distance applied to multi-cultivated rice in each region.

Table 3 .
Statistics of rice cropping areas in each province derived from MODIS.