Oil Palm Tree Detection with High Resolution Multi-Spectral Satellite Imagery

1.1. The Need for Oil Palm Tree Detection

Oil palm is a perennial monocot plant with a long generation period of about 25 years. It is an important economic crop in tropical areas of Southeast Asia, Africa and South America. Oil palm is the most efficient oil seed crop in the world when compared to other oilseed crops such as soybean, sunflower, rapeseed or canola, peanut, cottonseed, coconut, olive, etc. One hectare of oil palm plantation can produce up to ten times more oil than other leading oilseeds [1]. Palm oil is used in a wide variety of food products and is a raw material in the production non-food products. Global demand for palm oil has grown exponentially during the last 50 years. Palm oil accounted for more than 30% of all edible oil production in 2012/2013.

Accurate inventories and monitoring of oil palm areas are critical for plantation management and plant area expansion. Remote sensing is one of the most reliable measurement tool for accurate monitoring over large areas. Tan et al. [2] demonstrated some techniques for studying the age of oil palm trees using UK-DMC 2 and ALOS PALSAR remote sensing data. Age is one of the important factors influencing the production of fruit bunches, which begins at the time a tree reaches maturity at about 3 years, with production reaching a peak between 6 and 10 years. Tan found that texture measurement and fraction of shadow are useful for studying the age of oil palm trees. Santoso et al. [3] investigated the potential of high resolution satellite imagery for detecting oil palms infected by basal stem rot disease, and for mapping the disease. Basal stem rot disease can cause of up to 50% loss of the palms in a plantation. The article reported that vegetation indices derived from remote sensing imagery can discriminate between healthy and infected oil palms with a mapping accuracy of 84%. Tan et al. [4] estimated oil palm productivity on a large scale using existing remote sensing based models. These models use Light Use Efficiency (LUE) logic. This study concluded that the existing remote sensing based models need to be modified to incorporate meteorological inputs, environmental constraints, and maximum LUE as a constant value, in order to accurately estimate oil palm productivity. Thenkabail et al. [5] developed a comprehensive methodological framework for oil palm carbon and biomass estimation using a multi-date wet and dry season high spatial resolution IKONOS dataset. The study site was West Africa, where oil palm trees are found primarily in shifting cultivation agro-forest farms. Morel et al. [6,7] tried to estimate and evaluate the potential for monitoring aboveground biomass in forest and oil palm plantation in Sabah, Malaysian Borneo using Landsat ETM+ and ALOS PALSAR data.

1.2. Related Methods

The purpose of oil palm tree detection is to locate or count oil palm trees on map, aerial image or satellite imagery. This information is the key factor for oil palm plantation management and monitoring in each plantation area. Knowing the amount of oil palm tree in each plantation area can facilitate the yield prediction. The simplest way to detect oil palm is to manually mark the oil palm trees on images or field surveying using GPS to collect the positions of oil palm trees and display their position on the image. However, if there a lot of oil palm plantation areas, which can be large and contain more than 1000 oil palm tress, manual detection and field surveying can be a time consuming process and expensive. This case is where remote sensing technique.

Broadly speaking, the oil palm tree detection can be considered as the tree detection problem which is an on-going research problem in remote sensing. The goal of tree detection is to extract tree positions or tree crowns from given data, i.e., photograph or laser point cloud. Existing methods for tree detection can be grouped in terms of data used. This work concentrates on using high-resolution image instead of LiDAR. A drawback of using LiDAR for detecting trees is the requirement of high point density, which is costly and requires large data storage. Moreover, using high resolution images, oil palm tree can be easily recognized and identified on images.

Aerial photograph and satellite imagery have been the most widely used data for the application of tree delineation or detection. Wulder et al. [8–11] performed a series of studies on the use of a local maximum filter to detect tree location in high spatial resolution imagery, i.e., airborne and satellite-borne images. The local maximum filter is applied to the intensity image to find local peaks which are hypothesized to be the likely positions of trees. Similarly, Gebreslasie et al. [12] used the concept of local peak detection for locating individual trees in plantation forest in KwaZulu-Natal, South Africa. Their method is based on the use of Gaussian filter to smooth satellite image in order to eliminate noise. Furthermore, extremely bright and dark areas are removed by using natural break classification. In order to determine the window size for local peak detection, they employed the semi-variogram. Szantoi et al. [13] used Sobel edge detection and several (color) band combinations to detect downed trees (palms and deciduous) from very high spatial resolution (0.3 m) aerial imagery. This resolution is very close to the one the authors use in this study.

In stead of detecting only the position of trees, some approaches also delineating tree crowns. In [14], Culvenor presented a method for automatically delineating tree crowns in high spatial resolution image. The method first detects the tree crown centroid, (i.e., the image coordinates of the likely crown peak) by detecting local radiometric maxima. The local maximum pixels are then used to construct boundaries between all likely tree crowns. The digital number/brightness of the maxima pixels is used to compute the threshold for aggregating contiguous pixels within the same boundary defined by local minima pixels into a unique tree crown. The seed points for the aggregation are the maxima pixels. This work was tested on imagery of native Eucalypt forests in Australia.

Jiang and Lin [15] introduced a method for the delineation of individual deciduous trees in leaf-off aerial Ultra-High Spatial Resolution Remotely Sensed (UHSRRS) imagery. Individual tree locations are first detected by using mathematical morphology. Given the potential tree locations, tree boundaries are then identified by the marker-controlled watershed segmentation algorithm. The segmented tree regions are then classified as deciduous or coniferous by using a distribution free support vector machine, i.e., least-square support vector machine. The drawback of this method is that it required ultra-high resolution image which can be impractical.

Ardila et al. [16] presented an active contour based method for the tree detection in urban areas using very high resolution image. The process starts with the computation of a normalized difference vegetation index image. The local peaks on the image are then detected. The likely boundaries of the tree are defined as ellipses and extracted by fitting a two-dimensional Gaussian equation around the detected local peak. The obtained ellipses are then refined by using active contour model to get more realistic tree boundaries. Since the active contour is used to refine the shapes of tree crown boundaries, this approach can be impractical when ultra-high resolution image is not used.

Shafri et al. [17] proposed a method for detecting oil palm trees in an image of an oil palm plantation site. The input data was airborne hyperspectral imagery which provides rich information for removing non-oil palms from the image. To detect oil palm trees, this method employed several processing algorithms. Texture analysis, i.e., grey-level co-occurrence matrix was used to make the oil palm area more homogeneous. Edge enhancement, particularly the Sobel edge filter, was then applied to delineate edges. Oil palm trees were then extracted by using a threshold technique. The shapes of the extracted oil palm trees were refinned by morphology reconstruction and the number of oil palm trees was counted by blob analysis. The drawback of this approach is that only hyperspectral image is required.

The study reported in this paper deals with the problem of oil palm tree detection from an image of an oil palm plantation area on high resolution multi-spectral imagery. That is, a novel algorithm is proposed. An advantage of the proposed method is that only high resolution multispectral satellite image is required in contrast to [15–17] which require very high resolution or hyperspectral images. The concept of the proposed method is based on the local peak detection on a vegetation index image. Since there are various existing vegetation indices, we propose to use the vegetation index having the most power in discrimination between oil palm and non-oil palm objects. With the feature selection process, the proposed method is versatile because some indices cannot be extracted from some imagery or when the radiometric information is transformed. For example, the radiometric information of pan-sharpened image is transformed; radiance transformation cannot be computed. As a result, NDVI cannot show the true phenology of trees.

The remaining of this paper is organized as following; in Section 2, the proposed method and also the dataset used in the experiment are introduced. The experimental result and performance evaluation are reported in Section 3. The detection accuracy and discussion will be in Section 4. The conclusion of this work will be presented in Section 5.

2.1. Dataset

In our experiments, to evaluate the performance of our proposed palm detection method, QuickBird imagery of Phang-nga province in the southern region of Thailand centered at 8.4397°N, 98.5181°E was utilized. This area has been traditionally cultivated intensively with oil palm and rubber tree. The pan sharpened imagery has ground sampling distance about 60 cm per pixel and 4 bands including blue, green, red and near-infrared. With this spatial resolution, it can be recognized and identified an individual tree crown whose radius is more than 2 or 3 m. We manually selected the palm plantation areas and marked their boundaries. Moreover, pan-sharpened WorldView-2 imagery of Phuket province centered at 7.8900°N, 98.3983°E was also used to demonstrate the performance of the proposed method. Both of the imagery used in this work were already ortho-rectified without radiance transformation.

2.2. Workflow of the Oil Palm Tree Detection Method

In this paper, a method for detecting individual oil palm trees using high-resolution multi-spectral images, particularly in plantation areas is proposed. The images used in this methods are already ortho-rectified and can be either RGB or RGB with NIR. After the boundary of the plantation area has been manually marked, a local peak detection approach is applied. In this approach, the tree crown apex is assumed to have locally maximal digital number, e.g., reflectance or vegetation index, greater than that of surrounding neighbours within the window.

The local peak detection is applied on a vegetation index image derived from the multi-spectral image. There are several predefined indices. The best index is chosen, where best is defined as the index that maximizes inter-class dissimilarity oil palm trees and non-oil palm samples.

As explained previously, the proposed method is based the use of local peak detection which can be implemented by using 2D window function. The determination of the window size is hence very important. We propose to use a fixed-window size for any one image, which is estimated by using 2D semi-variogram technique. This approach is used because the pattern of oil palm plantation area is normally well structured i.e., the distances between an one oil palm tree and its neighbors are close to constant.

After the appropriate vegetation index feature and window size are determined, the method proceeds to the enhancement of the index image. In order to maximize the distinction between oil palm and non-oil palm trees, a rank transformation is applied to the vegetation index image. The local peaks are then detected with a non-maximal suppression algorithm using the window size obtained from the 2D semi-variogram technique. The work flow of the proposed method is summarized in the Figure 1.

2.3. Calculation of Vegetation Indices

To detect objects of interest on an image, the first step is commonly to reduce the amount of information by using image features. Since the proposed method is oriented toward green segmentation, i.e., palm trees, we believe that vegetation indices are good features for this purpose. In literature, there exist various spectral indices defined for inferring surface properties, for example, vegetation health, burned area and fire severity. In other words, these indices were usually defined for specific tasks [18,19]. To choose an appropriate feature for palm tree detection, some vegetation indices presented in the literature are reviewed. The discussion on how to select an appropriate index is deferred to Section 3.2.2. Note that the purpose of using vegetation indices in this work is not to find the true phenology of tree but just to differentiate oil palm tree and background. Therefore, the radiance transformation of the images is not required in this work because the true phenology of the tree is not needed in this work.

In [20], Guijarro et al. utilized various indices in the segmentation of agricultural images. Such indices include Excess Green (ExG), Excess Red (ExR), Excess Blue (ExB) and Excess Green minus excess Red (ExGR). Given an image in the RGB color space, the computation of those indices starts from the normalization of spectral red (R), green (G) and blue(B) components at a certain pixel:

As a consequence, the normalized spectral r, g and b are in the range [0, 1]. The indices are then computed as following:

In [21], Woebbecke et al. used green and red bands to compute Normalized Difference Index (NDI) for distinguishing plants from the soil:

The aforementioned indices are very useful for segmenting vegetation on an image when only red, green and blue spectral bands are available. However most remote sensing satellites provide images in the RGB-NIR space. Given an image with higher spectral resolution, the NIR spectrum can be taken into account in defining vegetation indices. The simplest vegetation index using NIR spectrum is the simple ratio [22]:

It is the most commonly used vegetation index. Its values vary between −1 and 1. The NDVI can be modified to obtained Transformed Vegetation Index (TVI):

In general, the radiance transformation should be performed before computing vegetation indices in order to obtain true vegetation indices. However, in this work, the vegetation indices are not used for monitoring the true phenology of the oil palm tree; they are actually used for differentiating oil palm tree and backgrounds. The true vegetation indices are not hence required.

2.4. Feature Selection

As it was discussed previously, different vegetation indices have been defined for specific tasks. Their relative performance in palm tree detection problem can vary. Another reason for performing feature selection is that some index cannot be always computed. For example, when NIR band is not available, NDVI, which is widely used in remote sensing for separating soil and vegetation, cannot be computed. The feature selection process will select a suitable index from all of the indices that can be computed from the image. The input data for the proposed method is hence not limited to only multi-spectral or satellite imagery. That is, RGB image or aerial photograph can be used. Moreover, the feature selection can make the proposed more robust when the radiometric contents are distorted or transformed, i.e., the case of pan-sharpened images. In this case, the computed indices are not the true ones; the commonly used index, e.g., NDVI is a suitable one for classification.

In statistical pattern classification, a suitable feature is one that clearly separates or has high discriminating power between classes. In this case, oil palm trees are needed to be separated from background. To evaluate the discriminating power of a feature, let us consider this process as the binary classification problem. The aim of this problem is to classify two objects, i.e., palm tree and background. A basic approach for measuring the discrimination ability of a vegetation index is to use the dissimilarity between two histograms generated from the index values for the two classes to be distinguished (palm tree and background samples, in our study). Such dissimilarity can be measured by the dissimilarity/distance metric between two histograms.

In this paper, only the measures follow the metric properties are used in the evaluation because of the symmetry and ease of the interpretation. Measures that satisfy this constraint are listed in Table 1. The performance evaluation for each spectral index will be reported in Section 3 and the vegetation index having highest discriminating power will used in our experiment. Since there are various measures, the discriminating power, in this study, is therefore formulated as the summation of all histogram measures. The concept of using the summation of measures to evaluated the similarity between histograms is similar to the work of Stricker and Orengo [26] using the sum of the weighted distances of the histograms first three moments as the similarity measure. In [27], the using of multiple dissimilarity measures in classification problem, i.e., 2-class classification was analyzed. It was reported that a better classification performance can be obtained by the combination of dissimilarity measures.

2.5. Estimation of Appropriate Window Size

To estimate an appropriate window size for rank transformation and non-maximal suppression, it is assumed that the spatial pattern of palm tree is well structured, e.g., two consecutive palm trees are equally spaced. The semi-variogram technique is employed to extract the spatial pattern of the palm tree farm. As a result, an approximate distance between two neighbor palm trees can be obtained and used as the upper-bound for the window size of rank transform and non-maximal suppression.

The semi-variogram, which is referred to experimental variogram in this paper, is a technique for image texture characterization (Atkinson and Lewis [28]). The spatial arrangement of objects in image can be extracted including the spatial distance between objects, surface roughness or orientation of the pattern. As an example, Trias-Sanz [29] used a 2D semi-variogram to analyze the structures of plantation fields. In [30], Szantoi et al. used semi-variogram to estimate the window size for different texture features. The semi-variogram V of an n-band image I can be obtained by first computing:

The normalized semi-variogram V of the image I can then be computed by first finding the minimum and maximum of D while ignoring the values near the origin:

As a result, the semi-variogram V is then normalized to the range [0, 1]. The average distance between local peaks on the 2D semi-variogram V and its neighbors is assumed to be the average distance between palm trees and will be used as the window size for rank transform and non-maximal suppression.

2.6. Rank Transformation

The rank transformation was originally developed for the problem of stereo matching, by Zabih and Woodfill [31] and Hirschmuller and Scharstein [32]. Precisely, it was proposed to improve the robustness of template matching. The expected result from the rank filter is the removal of blurred object borders. Hence, this filter may provide benefits in our application because the discontinuity between the palm and background intensities will be increased.

The concept of the rank filter is based on the relative order of the intensity values in a local region. The rank filter replaces the intensity of a pixel p by the number of pixels in the neighborhood of the pixel p whose intensity is less than the intensity of the pixel p. It can be formally defined as:

One of the factors affecting the result of the rank transformation is the window size. We proposed to use the average distance between a local peak and its nearest neighbor, which is obtained from 2D semi-variogram, because when the window is located at a local peak, which is hypothesized to be the position of palm tree, the window will not cover other local peaks. Thus the intensity of the single peak where the window is located will be dominant. In other words, if the window size is too large, the window will cover many local peaks. As a consequence, the result from the rank transform will be a smooth image and local peaks can disappear.

Promising results from applying the rank filter to enhance our image features in palm tree detection will be illustrated in Section 3.

2.7. Local Peak Detection

The result from the feature enhancement methods is then used for the detection of palm trees. Since the possible region of palm trees are the regions with high digital number, a local peak detector is used. However, a local peak detector cannot be directly applied to the image because the shapes of some areas on the image surface can be flat-topped hill. That is, the peak cannot be well localized if we are finding local peak on flat-topped hill area. To diminish this problem, the Gaussian filter is first applied to the image. Moreover, the Gaussian filter, which is a low pass filter, also reduces the potential problem of multiple peaks a local area.

In this proposed work (Algorithm 1), a non-maximal suppression algorithm for detecting local peaks is employed. Particularly, the algorithm proposed by Neubeck and van Gool [33] is used. The algorithm seeks to find the local peak in (2n + 1) by (2n + 1) windows; the parameter n can be estimated such that 2n + 1 is close to the average distance between oil palm trees, as calculated by the semi-variogram technique explained in Section 2.5.

The crucial parameter for the non-maximal suppression is n. If n is too large, the searching window will contain many local peaks, only one of which will be selected by non-maximal suppression. Particularly, the non-maximal suppression is analogous to finding local peak in 2n by 2n window. Therefore, we suggest to used the half of the average distance between a local peak and its nearest neighbor obtained from the semi-variogram for the parameter n because the searching window is less likely to cover other local peaks which are hypothesized as the locations of other oil palm trees.

Algorithm 1:. Non-maximal suppression in 2 dimensions proposed in [33].

Algorithm 1:. Non-maximal suppression in 2 dimensions proposed in [33].

2.8. Performance Measures

To evaluate the performance of the proposed method, we compare the detection results from the proposed method with a reference solution based on manual interpretation. The comparison with the reference data can be quantitatively shown by common performance metrics. Before proceeding to the discussion about the performance metrics, let us first introduce some important terms which are. True Positive (TP), False Positive (FP) and False Negative (FN). A TP, in this context, is an oil palm tree that is correctly detected by the proposed method. A FN is an oil palm tree that is not detected. A FP represent a pixel that is identified as an oil palm tree but which is really something else.

For the object detection problem, common performance metrics for evaluating a detection strategy are the precision, recall and the F-measure [34–36]. These metrics are defined as following:

3.1. Feature Selection

As explained in Section 2.4, the approach for feature selection used in this paper is based on the dissimilarity metric comparing the histograms of palm tree and background, for the various vegetation indices which are considered. To generate the histograms of palm tree and non-palm-tree areas, palm trees and background from subset images are sampled and all of the indices presented in Section 2.4 for each pixel is computed. Then, the histograms of the computed indices are generated. Figures 2 and 3 show the histograms of each index using QuickBird imagery. Next, the metrics comparing the histogram of palm-tree and background pixels are computed. Note that all the metrics presented in Table 1 except the intersection metric give a large value when two histograms are dissimilar. The intersection metric will give a large number if two histograms are similar. For example, the intersection metric between two identical histograms will be 1. To facilitate comparison of the metrics, we subtracted the intersection metric from its maximum of 1 and use that transformed value in our tables.

Table 2 reports the dissimilarity metrics from QuickBird imagery. The analysis shows that the ExGR index produced the lowest total metric value. In this index, palm tree and background are very similar. The NDI index gave the highest metric and thus can be considered to be the best feature for use in the detection process for QuickBird imagery.

3.2. Oil Palm Detection Results

The appropriate window size by using semi-variogram approach are determined as described in Section 2.5. The lag u is defined as u ∈ {−u_m,...,u_m} × {−u_m,...,u_m} and the image I has 4 bands, i.e., blue, green, red and near-infrared. The parameter is set to be 32 which covers a sufficient number of oil palm trees because the approximate diameter of the oil palm canopy is about 3–4 m which corresponds to 8 pixels. This means to say, we want to analyze the spatial variation within the area of 64 by 64 square pixel which cover enough amount of oil palm tress. Figure 4 shows an example of the image used in the performance evaluation and its semi-variogram. The semi-variogram shows the approximate spatial variation of the palm trees. The average distance between a local peak and its nearest neighbor is used to determine the window size of rank transformation and the non-maximal suppression. For this image, the average distance is about 17 pixels. Remind that the detection of plantation area’s boundary is out of the scope of this paper, The boundary can be obtained using an existing automatic detection algorithm or GIS database.

In order to reduce noise, each band of the image is then smoothed with the Gaussian filter prior to computing the index. The index image is then transformed using rank transformation. An example result from the rank transformation on the index image is illustrated in Figure 5. The boundary of oil palm tree region is drawn with solid red line. As it can be observed, the rank transform enhances the discontinuity between hypothesized palm tree locations and background. After the rank transform is applied, the local peaks on the index image are then detected by using the non-maximal suppression. The oil palm trees detected from the image shown in Figure 4 are illustrated in Figure 6. It can be observed that the spatial variation of detected oil palm trees is close to the computed semi-variogram. The palm tree region is heterogeneous, i.e., consisting of non-palm trees and other features such as buildings which can be false positives in the detection results. The detection results from two of the data subsets using QuickBird imagery are illustrated in Figure 7. In Figure 8 shows the detection results from two of the data subsets using WorlView-2 imagery.