Paddy field classification with MODIS-terra multi-temporal image transformation using phenological approach in Java Island

Received Jul 3, 2018 Revised Sept 6, 2018 Accepted Oct 11, 2018 This paper presents the paddy field classification model using the approach based on periodic plant life cycle events and how these elevations in climate as well as habitat factors, such as elevation. The data used are MODIS-Terra two tiles of H28v09 and H29v09 of 2016, consist of 46 series of 8-daily data, with 500 meter resolution in Java region. The paddy field classification method based on the phenological model is done by Maximum Likelihood on the transformed annual multi-temporal image of the reflectance data, index data, and the combination of reflectance and index data. The results of the study showed that, with the reference of the Paddy Field Map from the Ministry of Agriculture (MoA), the overall accuracies of the paddy field classification results using the combination of reflectance and index data provide the highest (85.4%) among the reflectance data (83.5%) and index data (81.7%). The accuracy levels were varied; these depend on the slope and the types of paddy fields. Paddy fields on the slopes of 0-2% could be well identified by MODIS-Terra data, whereas it was difficult to identify the paddy fields on the slope >2%. Rain-fed lowland paddy field type has a lower user accuracy than irrigated paddy fields. This study also performed correlation (r2) between the analysis results and the statistical data based on district and provincial boundaries were >0.85 and >0.99 respectively. These correlations were much higher than the previous study results, which reached 0.49-0.65 (hilly-flat areas of county-level), and 0.80-0.88 (hilly-flat areas of provincial level) for China, and reached 0.44 for Indonesia.

are satellite images that have become public domain, have the spatial resolution of 250 metres, 500 metres, and 1 km with four recording times per day at the same location on earth [3]. However, the geographic condition of Indonesia in the tropics reduces the availability of MODIS image quality because the territory of Indonesia is covered by MODIS image which is often covered by clouds and haze.
Many researchers have developed cloud extraction techniques from various satellite imagery [5], developing Mosaic Pixel Based for Landsat data [6], and Mosaic Tile Based of Landsat-8 to obtain minimal cloud cover [7]. The University of Maryland has performed time-series transformation of MODIS and Landsat data in Congo Basin [8], Landsat data from 1985 to 2012 in Eastern Europe's [9]. Image transformation is done by histogram-based metrics approach, such as average, and by sequential metric approach [8]- [10].
The collaborations among LAPAN, Ministry of Environment and Forestry (MEF), and World Resources Institute (WRI) have tried to transform time-series images with MODIS-Terra and Landsat-8 OLI data. The data used are MODIS from 2000 to 2017, and Landsat from 2015 to 2017. The experiment was conducted with histogram-based metrics approach, such as average, and time-sequential metrics approach, such as regression. Data time-series level used is single image transformation, annual image transformation (metric level-1), and inter-annual image transformation (metric level-2). The results show that image transformation can be used to identify changes in forest coverage, without having to analyse individual change overlays [11].
The development of multi-temporal image is still conventionally done, is by analyzing annual data individually based on reflectance, and then compared with the annual data in different time to obtain the phenomenon of land use change. Trend or land use change analysis using these conventional methods takes longer and requires specific application skills.
The use of the image transformation approach to analysing paddy field mapping has begun to be used evolutively. In the first generation, paddy field mapping was done with the use of category one algorithm, such as data reflectance and image statistic-based approaches. The next development emerged as the second generation using vegetation index and enhanced image statistic-based approaches. In the third generation development, the vegetation index or RADAR back-scatter-based temporal analysis is used [12]. Recent developments began using the phenological of paddy through remote sensing recognition of key growth approach phases [13]- [16].
A phenological model is an approach based on periodic plant life cycle events, and how these are influenced by seasonal and inter-annual variations in climate, as well as habitat factors, such as elevation. Variables used for paddy field classification were developed from data reflectance, Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) [17], [16], [14]. Normalised Difference Water Index (NDWI) [18]- [20], and Land Surface Water Index (LSWI) [17], [19]. The use of vegetation indices (e.g. NDVI, EVI, LSWI, NDWI) make the accuracy of land cover/use increased compared with the original reflectance [21], [22].
Phases of data used were developed from single image in early growing season before transplanting (Fang in 1998), multi-images from seedling and ripening stages [23], multi images in growing season [24], multi images in transplanting and tillering stages [25], all available images [26], and multi images in early rice growing season [27].
Xiangming Xiao (2005) had studied the paddy fields in South and Southeast Asia using multitemporal MODIS images. He had mapped for 13 South and South East Asian countries with a MODIS 500-meter spatial resolution over 8-day data in 2002. Phenological models were used in the study. Paddy rice fields were characterised by an initial period of flooding and transplanting, during which a mixture of surface water and rice seedlings exists. He applied a paddy field mapping algorithm that uses a time-series of MODIS-derived vegetation indices to identify the initial period of flooding and transplanting in paddy fields, based on the increased surface moisture. The resultant MODIS-derived paddy field map was compared to national agricultural statistical data at national and sub national levels. The results show a similarity with the location of the paddy field as a whole, but there are variations in some of the locations on the topic. Although the results still need to be done further research, the method and use of MODIS data provide potential [16], [15], [14].
From the various developments of the above research, paddy field mapping or classification using phenological approach resulting from image transformation with the combination of reflectance and MODIS-Terra annual multi-temporal image index has not been done yet. The success of finding fast, precise and accurate procedures to assist in monitoring a land area and rice production will greatly assist in the planning and implementation of food resilience and independence programs [28], [2].
The Java Island was selected as the study area. This is the most populous island with 146,675,400 inhabitants, or 56.7% of the Indonesia total population of 258,704,900 in 2016. This island is a buffer area of rice production with a production capacity of 38,970,026 tons (51.7%) of all rice production in Indonesia  [29]. This island is the most dynamic island among large islands in Indonesia, due to population density and rapid development of the region [1]. For that, the existence of a good monitoring tool becomes more necessary.
This main objective of this study is to test the paddy field classification by using MODIS-Terra multi-temporal image with a combination of reflectance (Red, NIR, SWIR-1) and index (NDVI, NOAI), based on phenological approach. This research is different from what has been done by Xiangming Xiao (2005), he only uses parameter index (NDVI, EVI, and NDWI). The similarity with Xiao Xiangming's research is the use of category four that is phenological-based models through remote sensing recognition of key growth phases with the annual period. The phenology of rice plant in Java Island, characterized by a) at the beginning of planting always flooded with water, b) existence of up and down trend in vegetation index and open index area, and c) change of dynamic land cover every year, water phase, vegetation and fallow land is shown by the variance of reflectance and index values.

Data
Primary data used in this research are the 8-day's reflectance of MODIS-Terra of 2016, with pixel resolution 500 meter of data from NASA's Land Processes Distributed Active Archive Centre (LP DAAC). This data consists of 46 series data of paths H28v09 and H29v09 of Java Island, covers band Red, NIR, and SWIR-1. While the two types of secondary data from the Ministry of Agriculture (MoA) were used, in the form of Paddy Field Map in 2012 as shown in Figure 1 and paddy field area base on district and province statistical data report in 2015. The paddy field map obtained from the delineation of high-resolution satellite imagery, but the annual statistical report from field estimation.

Image classification
The method used in paddy field mapping for this study is an analysis of annual multi-temporal imagery with a phenological approach. The analysis process includes two stages: a) the steps to extract the 8-day's MODIS data into multi-temporal feature information image, performed using image transformation, and b) the step to classify the transformed multi-temporal feature image with the Maximum Likelihood Classification (MLC) approach [30]. The pre-processing of MODIS image is executed before a multitemporal transformation, which includes cloud masking, time-series filtering, and interpolation of blank data due to the cloud. Pre-processing is done to minimize the image of the cloud cover [31], [32]. Diagrammatically, the illustration of the paddy field classification model with a phenological model using MODIS-Terra multi-temporal image transformation is illustrated in Figure 2.

Image transformation
There are three types of reflectance images with the MODIS-Terra multi-temporal image transformation processed in this study, those are three reflectance images of SWIR-1, NIR, and Red; and two index images, namely Normalized Difference Vegetation Index (NDVI), and Normalized Open Area Index (NOAI), with the following formulas. The selection of index formula that is the minimum of reflectance of NIR and the minimum reflectance SWIR-1 done to highlight paddy field when flooded at the time of planting. The ups and downs of NDVI indicate the phenology that during the rice growing period, the NDVI value will increase during flooding to the vegetative stage, while at the maturation stage of the paddy, the value of NDVI will decrease.
where SWIR_1, NIR and red are reflectance in shortwave infrared 1, near infrared and red, respectively.
The three reflectance images and two index images are then transformed by ten algorithms. The result of a multi-temporal image transformation in the form of the new image is called metric image. An extraction algorithm for obtaining metric image can be calculated based on the statistic value, regardless of recording time or with respect to recording time sequence. Image metric is an image (feature) that contains information in accordance with the needs of the application. There are ten types of algorithm formula to obtain image feature metric, that is: 1) Average all clear pixels 2) Deviation standard all clear pixels 3) Minimum of all clear pixels

6)
A deviation standard of 70% clear pixels

8)
Maximum of 70% clear pixels where x ̅ =average; x_i=the value of the ith sample; 〖sx〗_i=the value of the ith sample after being sorted from small to large; n=number of samples; δ=deviation standard; However, in this study only feature metric algorithm correlated with paddy field growth phenology in Java are selected, that is 1) average of 70% clear pixels, 2) standard deviation of 70% clear pixels, 3) minimum of 70% clear pixels, 4) maximum of 70% clear pixels, 5) amplitude up, and 6) amplitude down. The process of image transformation and its algorithmic representation is shown in Figure 3 and  The reflectance image of band Red, NIR, and SWIR-1 are classified by MLC to obtain the distribution of paddy (rice field) and non-paddy (non-rice field) field. While the index images of NDVI and NOAI are also classified with MLC to obtain the distribution of paddy field and non-paddy field. The combination of reflectance images of band Red, NIR, SWIR-1 and the index images of NDVI and NOAI are also classified with MLC to obtain the distribution of paddy field and non-paddy field. The results of the three classifications with the input of metric reflectance image, metric index image, and also the combination of reflectance image and index image were calculated by its categories and compared its accuracy with the reference of Paddy Field Map of Java Island in 2012 from the MoA.

RESULTS AND DISCUSSIONS 3.1. Analysis of the image transformation result
There are several examples of the annual metric reflectance image and the annual metric index image of MODIS-Terra multi-temporal image transformation results in 2016. From the RGB color composite of the transformed reflectance image analysis in Figure 5 it is known that the blue in Figure 5(a) shows the dominance of the water content for a year in the multi-temporal image. While the blue in Figure 5(b) shows the dominance of paddy fields when inundated in the multi-temporal image. While from the analysis the transformed index image in Figure 5(c), it is known that the green shows the increase of NOAI index, while the red value indicates the decrease of NOAI index value, and the yellow indicates the rising and the falling value at different locations. From the Image in Figure 5(d) shows that the red indicates a high variation in NDVI index value, which means that land cover changes in one-year intervals are likely caused by the changes in land cover during paddy cultivation. The result of classification of the reflectance image using bands Red, NIR, and SWIR-1 by MLC to obtain the distribution of paddy field and non-paddy field is shown in The Image in Figure 6(a). The result of classification with the index image input of NDVI and NOAI is presented in The Image in Figure 6(b). While the classification result with the combination of image reflectance band Red, NIR, SWIR-1 and NDVI and NOAI index images is presented in The Image in Figure 6(c).

Paddy field area and accuracy assessment
The three classified images with MLC approach were then calculated by their categories and compared their accuracies to the Paddy Field Map from the MoA of 2012. Comparison of classification images with reflectance inputs, index inputs, and the combinations of reflectance and index inputs are presented in Tables 1-3. From the table, it can be seen that the correctness of MLC classification with the inputs of reflectance images of bands Red, NIR, SWIR-1 shows that the paddy field is quite high (20.66%), and higher than the classification with inputs of the combination of reflectance images Red, NIR, SWIR-1 and index images of NDVI and NOAI (19.92%).
The accuracy of the classification analysis results with MLC on the three inputs is calculated and presented in Table 4. From the analysis results, it can be seen that the use of the combination of reflectance of band Red, NIR, SWIR-1, NDVI, and NOAI images gives higher accuracy of 85.43% (producer accuracy 93.40%), compared to the use of only reflectance images of band Red, NIR, SWIR-1, that is 83.50% (producer accuracy 86.43%). The value of producer accuracy indicates that the paddy field that is detected by MODIS-Terra multi-temporal image has high accuracy. While the user accuracy value indicates that there are still many unidentified paddy fields using MODIS-Terra data. As an example, there is still 39.09% of paddy field which could not be identified using input reflectance or 36.18% of paddy field cannot be identified using input index, and 39.78% of paddy field cannot be identified using the input of the combinations of reflectance and index.
The error of classification which should be classified as the paddy field, but is identified as nonpaddy fields, which are 3.24% for classification with reflectance image input and 1.41% in the classification with index image input. The type of paddy field affected classification accuracy, in addition to topography factor.
Superimpose classification results with reflectance inputs, index inputs, the combinations of reflectance and index, is shown in Figure 7. From that figure, where the classification of paddy fields using reflectance input (indicated by Red), index inputs (indicated by Blue), the combinations of reflectance and index (indicated by Green). The area of paddy field using reflectance input is 76.10%, the area using index input is 78.67%, and the area using the combination of reflectance and the index is 71.89%.
The causes of the low user accuracy and the number of paddy fields unidentified by the MODIS-Terra image are slope and the types of paddy fields, the higher the slope the lower the accuracy. Besides, the rain-fed fields have lower accuracy than irrigated paddy fields. The difference of the year of reference data that is Paddy Field Map from the MoA (2012) and image of MODIS-Terra classification (2016) is also influencing the results of comparative analysis  Figure 7. Combination of RGB of various paddy fields of classification results Note: Red is classification using reflectance Red, NIR, SWIR-1 (total area=76.10%); green is classification using both combinations of reflectance and index (total area=71.89%); blue is classification using index NDVI and NOAI (total area=78.67%) The higher the slope factor the lower user accuracy, as shown in Graph of Accuracy relation to the slope in Figure 8. From the graph, it can be seen that the higher the slope the lower the accuracy, the user accuracy and the producer accuracy. It can also be seen that user accuracy decreases linearly, at slope <17% user accuracy will be >80%; and also producer accuracy drops sharply between slope 1% to slope 2%, then down linearly. This indicates that on the slope of 0-2% paddy fields can be well identified by MODIS-Terra data, whereas on slope >2% is difficult to identify.
Types of paddy field will also influence the results of classification accuracy. The result of the calculation of the influence of these types of paddy fields is shown in Table 5. From the table, it is known that the variation of accuracy differences over the identified paddy field types is quite large, ranging from 23% to 78%. And a rain-fed lowland paddy is the most difficult type of paddy with user accuracy of 22.92% identified by this approach. From the result of comparative analysis between the results of area calculation based on the district boundary of three inputs of classification of paddy field and the Paddy Field Map from the MoA showed that all of the classification results performed accuracy more than 0.85. Those are the correlation (r 2 ) of the reflectance input is 0.87, the index input is 0.89, and the combination of reflectance and index input is 0.85. The correlation among those approaches is shown in Figure 9. The result of comparative analysis between the results of area calculation based on the province boundary of three inputs of classification of paddy field and the Paddy Field Map from the MoA, showed that all of the classification results performed correlation (r 2 ) more than 0.99. Those are the correlation (r 2 ) of the reflectance input is 0.998, the index input is 0.991, and the combination of reflectance and index input is 0.998. The correlation among those approaches is shown in Figure 10. The results of calculation per province of these phenological models are shown in Table 6. The result of this study provides a higher classification result than the previous study conducted by Xiao     By using the combination of reflectance and index metric input, the result of analysis using statistical data of 2016 only performed 64.1% of the identified paddy field area. However, by analysis, which, considering spatial distribution data as shown in Table 5 (% area of paddy field=21.33%) and Table 3 (% area of paddy field=33.08%), the result of accuracy base on statistical data is 21.33/33.08 or 64.4%.  Table 6 it is known that classified paddy fields vary from differences input in classification base on province regions. Compared with tabular data from the MoA in the same year of 2016, there are 2.3 million ha (71.8%) can be identified by reflectance metric input, 2.7 million ha (84.8%) can be identified by Index metric input, and just 2.1 million ha (64.1%) can be identified by the combination of reflectance and index metric input. It means that the analysis using the combination of reflectance and index, input without considering spatial distributions performed the low result. The analysis using the combination of reflectance and index, input with considers the spatial distributions, as mentioned in Table 4 performed the higher result (overall accuracy of 85.43%). In line with the producer accuracy in Table 4 the higher paddy fields area using the combination of reflectance and index metric performed the highest of producer accuracy of 93.4 %.

CONCLUSION
The development of the paddy field classification model using phenological model using annual multi-temporal satellite image transformation was carried out on the Java Island. The satellite data used was MODIS-Terra image of two tiles of H28v09 and H29v09 of 2016, consisting of 46 series data of 8-daily data, with the spatial resolution of 500 meters. The paddy field classification method is carried out by classifying Maximum Likelihood on the reflectance, index, and the combination of reflectance and index data of the transformed annual multi-temporal satellite image.
The paddy field classification results show that the use of annual multi-temporal image with a combination of reflectance and index data provides the highest classification accuracy among other approaches, with 85.4% of the overall accuracy, 60.91% of the user accuracy, and 93.4% of the producer accuracy. The classification using only the reflectance provides 83.5% of the overall accuracy, 60.22% of the user accuracy, and 86.43% of the producer accuracy. While the approach using index, input alone provides 81.7% of the overall accuracy, 63.82% of the user accuracy, and 80.52% of the producer accuracy. The low user accuracy indicates that many paddy fields cannot be identified using MODIS-Terra data with this approach.
The low user accuracy or the low areas of paddy fields that cannot be identified by the MODIS-Terra image with this approach is caused by two factors, namely slope, and type of paddy fields. The higher slopes have lower user accuracy, and in the 0-2% slopes can be well identified by MODIS-Terra data, whereas on the slope >2 is difficult to identify. Rain-fed lowland types have a lower user accuracy than irrigated paddy fields.
From the analysis with the reference of the Paddy Field Map from the MoA, it is known that the result of this study performed correlation (r 2 ) >0.85 based on the district boundary, and >0.99 based on provincial boundary. These correlations were much higher than the previous study, which reached 0.49-0.65 (hilly-flat area of county-level) and 0.80-0.88 (hilly-flat area of provincial level) for China, and reached 0.44 for Indonesia. It is concluded that the phenological model using image transformation could be used to predict the area of paddy field in a higher accuracy, for flat area as well as hilly area.