Optimization of Parallel K-means for Java Paddy Mapping Using Time-series Satelite Imagery

Spatiotemporal analysis of MODIS Vegetation Index Imagery widely used for vegetation seasonal mapping both on forest and agricultural site. In order to provide a long-terms of vegetation characteristic maps, a wide time-series images analysis is needed which require high-performance computer and also consumes a lot of energy resources. Meanwhile, for agriculture monitoring purpose in Indonesia, that analysis has to be employed gradually and endlessly to provide the latest condition of paddy field vegetation information. This research is aimed to develop a method to produce the optimized solution in classifying vegetation of paddy fields that diverse both spatial and temporal characteristics. The time-series EVI data from MODIS have been filtered using wavelet transform to reduce noise that caused by cloud. Sequential K-means and Parallel K-means unsupervised classification method were used in both CPU and GPU to find the efficient and the robust result. The developed method has been tested and implemented using the sample case of paddy fields in Java Island. The best system which can accommodate of the extend-ability, affordability, redundancy, energy-saving, maintainability indicators are ARM-based processor (Raspberry Pi), with the highest speed up of 8 and the efficiency of 60%.


Introduction
In Java Island Indonesia, Agricultural produces such as paddy, maize, soybean, groundnut and mung bean are cultivated in paddy field which divided into two or three cropping seasonal [1][2].As the fourth-highest rice consumption country after Myanmar, Vietnam, and Bangladesh, Indonesia consume 1.62 kg rice per capita per week in 2014 [3].Total rice consumption has been rising faster than production, as the growth rate of national rice area and yield has faltered.Based on U.S Department of Agriculture-Foreign Agricultural Service/USDA-FAS [4], Indonesia ranks 3 rd in the world in regards to total rice production but has also been the world 7 th largest rice importer over the past 5 years, on average requiring over 1.1 million tons of imports per year.Consequently, food security and the pursuit of national rice self sufficiency have become main concerns of the government in Indonesia.In the other hand, the recent modeling of climates reports that there is a potential that may threaten food security in the near future [5].
Crop phenology monitoring at a regional scale can provide useful information for agricultural management to enhance crop yield via irrigation regulation or adjustments in crop cultivation systems [6][7].Ozdogan et al. [8] have reviewed the method and the result of remote sensing in irrigated agriculture by the variant of local climate.Observation satellite that provides continuous data and open access worldwide: Landsat, ASTER, MODIS, MERIS, ACHRR.Some recent studies have been proved that time-series analysis of optical type satellite describe the phenology of terrestrial vegetation which represented the condition and the land used change based on trends and seasonal change [9][10][11][12][13][14][15][16].In order to reduce noise which caused by clouds and other aerosol problem, discrete wavelet transform (DWT) is used to denoising and extracting the agriculture vegetation information from the terrestrial surface [9], [17][18].Decision support system for weed control based on a precision agricultural approach that has been done Spatiotemporal analysis on time-series imagery in long-term data needs a lot of computational resources and a sophisticated computer to provide dynamic agriculture cropping maps.Meanwhile, monitoring the current condition can be done by that analysis gradually and continuously.Trends in computer architecture exponentially improve and multi -core machines are said to provide high-performance at low-energy cost via multi-threading [20].In the context of parallel and distributed computing, energy efficient voltage scheduling for multi-core processors is an important issue [21].The improvement of energy efficiency not only depending on the problem size but also the degree of multi-threading [22].Low-power consumption computer such as Raspberry Pi (3.5 Watt) can perform as a desktop pc and used for image processing and weed fractal dimension processing [23].
This study investigated the computational method both parallel and sequential of kmeans unsupervised classification for long-term spatiotemporal data using both CPU and GPU as a processing unit.Both CPU and GPU have a potential for an efficient outcome, but both processing units have limitation.To optimize those problems, the performance of sequential and parallel k-means using Multi-CPU with raspberry pi 2, Multi-core with Intel Core-i7 3th gen and Nvidia GTX 1060 6Gb GPU were compared.Energy consumption and processing time are considered with minimization goal and maximum processing time constraint is the period of satellite images.

Research Method 2.1. Study area
The study area on this research involves paddy field in Java Island which defines by Geospatial Information Agency (BIG) of Indonesia.The method was validated and tested for paddy field in Banten, West Java, Central Java, and East Java Province in Java Island.Ground check location selected by class distribution and the percentage area of paddy field inside the MODIS pixel.Paddy field map as shown in Figure 1 was used as the mask for pattern classification analysis.

Satellite Imagery
The MODIS product used in this study is the Vegetation Indices (VI) Composite 16-day Global 250 m SIN grid V006 or MOD13Q1 product [24], which provided vegetation coverage in Java Island.The Enhanced Vegetation Index (EVI) as in equation ( 1) is embedded in the MOD13Q1 product.The EVI developed to optimize the vegetation signal with improved sensitivity in high biomass regions and improved vegetation monitoring through a decoupling of the canopy background signal and a reduction in atmosphere influences [11].

Spatiotemporal Analysis
The collection and preprocessing geo-spatial data from MODIS imagery using MODIStsp library which provide by Busetton and Renghetti [25] was conducted.MODIS image from LP DAAC [24] collected from 2010 until 2016 for every 16-days.Image tiles for Java Island located in h28v09 and h29v09.Both of those images were merged to provide single raster image of Java Island.Research procedure involve RAW data acquisition, pre-processing, denoising, pattern classification, and spatial analysis as shown in Figure 2.

Pre-processing and denoising
Preparation of MODIS image data begins with image mosaicking on tile h28v09 and h29v09.Clipping the image of the merging is done with the help of Java vector map obtained from Geospatial Information Agency (BIG) of Indonesia.The imagery at each time frame is arranged in order to produce raster data with z axis matrix as time domain.The matrix manipulation of raster data in the form of three dimensional matrices begins with the separation of wetland data with the help of wetland rice maps from BIG as Mask.Matrix with X and Y dimensions as location information, while Z as data in time domain is modified into two dimensional matrix.The modified matrix is a data table with information on the dynamics of the vegetation index in each column for each location.
Wavelet transforms are used to reduce the interference caused by clouds and other weather disturbance.Wavelets are employed at each data location to correct the dynamics of the vegetation value in the time domain.Mother wavelet used is a coiflet with expected results in the form of approximation function that describes the data at low frequencies.
where  is mother wavelet function, ) (x f j is approximation and ) (x g j is detail.

Parallel k-means
K-means clustering is used to categorize a data set based on its average value (ISODATA).Categorization of larger data sets takes a lot of time and intensive calculations.Modification of serial processing into parallel processing by employing MPI library can speed up processing time by optimizing the application of multithreading [26].The process starts from the initial centroid determination.Then the serialized determination of the data sketch against the centroid of the K groups and the calculation of the new centroid.The process will be done unt il the number of cluster changes per amount of data smaller than the threshold.Pseudo-code from K-Means serial shows in Algorithm 1. )) , ce( min(Distan 1 (4) From equation 4, k-means algorithm finds k data points on the instance space such that the mean square error (that is, the total distance of all instances to the nearest cluster center) is minimized [27].The sequential version of the k-means algorithm is depicted in Algorithm 1. Parallel modifications to k-means are performed on data sharing according to the number of processors involved [26].Calculations are performed on each processor, then collected for results.The proposed method for processing of k-means using MPI is shown in the Figure 3.

Optimization
The selection of the best system for image data classification computing is performed based on the performance of each parallel method.Parameters to be considered for optimization calculations include: processing times, speedup, efficiency, durability, extendability, affordability, redundancy, power consumption, maintainability, and simplicity.Increased speedup ( Sp) on the utilization of parallel processing is approximated by the Amdahl's Law equation.In the equation 5, speedup can be computed from the serial fraction.So that parallel efficiency ( ) can be obtained from the equation 6.In this study, the Amdahl's law calculation is performed on RPi and x86 systems because the number of processors on both systems can be determined based on the number of available devices.While on CUDA/GPU, the number of processors is determined by the type of device and performed with single card configuration not dual GPU or SLI.

Parallel benchmarking
Performance test results from parallel computing using raspberry pi (RPi) and x86 computer desktop are shown in Figure 4.The speedup improves by adding more processors/cores both on RPi and x86 system.These results show that the ability of the RPi is still far below x86, but the increase in speed of the RPi is higher.The speedup reaches 8 times with 14 processors working in parallel.While the speedup on x86 reached 7 times with 16 processors working in parallel.The result of the Amdahl's law approach used to assess the performance of parallel processing shows that RPi has higher efficiency than x86 as shown in Figure 5.The result of the performance projection based on the Amdahl's equation as shown in Figure 6 shows the parallel utilization opportunity in RPi higher than in x86.Then, the comparison is done with the three types of devices used.In addition, the memory capabilities of the GPU are limited, so it needs to be done by comparing the ability of each GPU in processing small and large capacity data.The results of this test are presented in Table 1.

Conclusion
The method for parallel processing of k-means clustering using MPI has been developed, tested and implemented.The utilization of parallel processing in k -means clustering for Java paddy mapping gives some advantages such as faster computational time, more robust system by redundant configuration, and higher efficiency by using an ARM based system.All of those architectures are can be use to perform fast computation for paddy mapping.The performance of parallel k-means can be improved by adding more processors or processing cores.From this research, the best system architecture which can satisfy almost all of the performance indicator is ARM-based processor (Raspberry Pi), which achieved the highest speed up and the efficiency.It is important to highlight that both CPU-based and GPU-based (CUDA) systems can be employed for parallel processing to classify the paddy fields patterns.

Figure 1 .
Figure 1.Paddy field maps of Java Island define by Geospatial Information Agency of Indonesia

Figure 2 .
Figure 2. Research framework in spatiotemporal analysis to provide annual cropping maps

Algorithm 1 K-Means clustering 1 :
Select initial centroid from K class 2: w hile d/n data < threshold do 3: Find nearest cluster for K class to the centroid 4: Recompute new cluster center 5: Count cluster change as d 6:

Figure 3 .
Figure 3.The proposed method for parallel processing of k -means clustering using MPI Optimization of Parallel K-means for Java Paddy Mapping Using... (Alvin Fatik hunnada)

Figure 4 .Figure 5 .
Figure 4. Performance of parallel k-means clustering on (a) Intel core i7 and (b) Raspberry Pi 2 using MPI

Figure 6 .
Figure 6.Performance projection based on the Amdahl's equation both in Intel core i7 and Raspberry Pi 2 using MPI

Table 1 .
Performance of parallel k-means clustering using CUDA accelerator