Granulometric Analysis on Remote Sensing Images: Application to Mapping Retrospective Changes in the Sahelian Ligneous Cover

Abstract

This paper illustrates how the use of mathematical morphology can be a powerful tool for the mapping of ligneous cover in semi-arid lands. Ligneous cover plays a fundamental role in Sahel semi-arid regions since this resource is vital to the resilience of rural societies and can be used as an indicator of socio-environmental conditions. Grey tone vertical images from Sahelian villages in 1975 and 2010/2011 have been selected to perform a diachronic analysis to test the method. Granulometric profiles have been calculated for each pixel and then an unsupervised classification has been performed to obtain k classes that account for ligneous patches of different sizes. This method is particularly successful when the most recent images are used, given that these have better contrast and sharpness. Nested classifications were required to accomplish the ligneous mapping of images from 1975. The accuracy assessment for the most recent images classifications shows satisfactory results. The classification of ligneous cover according to different sizes is important for a better understanding of the ligneous dynamics.

1. Introduction

This paper presents a methodology for ligneous cover mapping by means of mathematical morphology. Combined use of granulometric analysis on grey tone images and unsupervised classification is proposed to detect ligneous patches of different sizes in semi-arid lands.

High-resolution imagery over the Sahelian region has been used to test the method. In semi-arid regions, the ligneous populations are often perceived as permanent inter-annual species so that changes in their dynamics are deeply linked to environmental conditions [1]. In the Sahel the ligneous species play an important role in the resilience of rural societies: they increase soil fertility, provide food and fodder, prevent from soil erosion and contribute to the diversification of tasks within the community [2]. Thereby, the mapping and comprehension of spatial changes in the ligneous distribution is an important challenge. The analysis of the size and the distribution of the ligneous patches from high resolution images is an appropriate way of understanding its dynamics [3].

A variety of methods have been proposed to automatically detect the tree-crowns from passive remote sensing images [4,5]. These methods propose the detection and delineation of tree-crowns by classical image processing techniques such as binarization [6] or object-oriented classification [7,8]. However, mathematical morphology has rarely been suggested for this task [9,10], and it could be useful not only for the detection of ligneous cover but also to classify it depending on its size following a mathematical logic [9].

Granulometry is a well-known technique to provide the distribution of object sizes of an image [11,12]. In this paper, successive closings by reconstruction have been performed in order to extract black patches from grey tones images. This technique permits to account both for the spatial and spectral characteristics of the image structures [13].

Some studies have proposed the use of granulometric profiles (or morphological profiles) for image classification [14,15,16,17]. A method for macro-texture mapping of vegetation density from granulometrics has been proposed [1,18]. However, the majority of the studies using morphological profiles have been based on urban imagery [14,15,16,17,19]. A successful application of mathematical morphology to tree-crown detection and characterization will be demonstrated throughout this paper.

An unsupervised k-means classification [20] has been performed over the morphological profiles extracted from the grey tone images. Although some redundancies can be expected from the classification of morphological profiles [15] it will be shown that choosing the adequate parameters can lead to satisfactory results.

The main aim of this fieldwork has been to find an appropriate method to analyze the retrospective changes of the ligneous cover dynamics in the Sahel using aerial and satellite images. The method shows satisfactory results and highlights the benefits of mathematical morphology for ligneous cover mapping.

This article is structured in sections. In the first section, a brief overview of the data, the morphological operators and the unsupervised classification used is described. In the following section, the results of this method applied to two different Sahelian villages are presented. In the last section a brief discussion on the topic is proposed.

2. Materials and Methods

2.1. Data and Background

High-resolution vertical images provided by aerial missions and satellites are an essential source of information for understanding the dynamics of the ligneous cover. The images presented in this article are scanned aerial images from an IGN mission (40/600) during March 1975 (scale 1:62,500) and screen captures from Google Earth Pro (this option is available directly from the Google Earth Pro menu). To allow comparison between both images the screen capture from Google Earth Pro (sub-meter resolution) has been resampled to the spatial resolution of the 1975 image. The final spatial resolution of the images is circa 1.3 m. This research is part of the results of a mission to southwest Niger held in November 2009 [21], which was intended to give some clues about the spatialisation of forage resources. Ligneous populations are a fundamental piece of the social-ecological resilience in this region.

2.2. Morphological Analysis

The technique presented in this section generates granulometric profiles for each image pixel. The granulometric profiles take into account the spatial and the spectral characteristics of the image structures.

The classical morphological operators, erosion and dilation, are applied to an image using a set of known shape and size called the structuring element (SE). The size of the SE is chosen according to the size of the structures to be extracted from the image. The closing is a morphological operator that intends to recover the image structures sizes after dilation. This is accomplished by eroding the dilated image. The closing by reconstruction is a special morphological operator that preserves the shapes of input images [22]. This is achieved by dilating the image with a SE of size λ followed by erosion with a SE of size 1. Then erosions of size 1 are repeated on the resulting maximum of the initial image and the eroded image until idempotence is reached. On high-resolution grey tone images, the ligneous cover in semi-arid lands is characterized by a peculiar punctiform texture. A closing by reconstruction with a SE being a disk of size λ should erase the dark spots smaller than λ during the dilation not recovering them during the subsequent reconstruction. Thus, successive closings by reconstruction can be used to retrieve image structures of different size. The granulometric density (Equation (1)) provides an efficient description of the size of the image structures [1].

The granulometric analysis is performed by successive closings by reconstruction with SE of increasing size expressed by Equation (1):

$$g{c}_{\mathsf{\lambda }}=\frac{[{\mathsf{\varphi }}_{\left(\mathsf{\lambda }\right)}(\mathrm{I})-{\mathsf{\varphi }}_{(\mathsf{\lambda }-1)}(\mathrm{I})]}{\mathrm{I}}\forall \mathsf{\lambda }\ge 1$$

(1)

where ${\mathsf{\varphi }}_{(\mathsf{\lambda })}(\mathrm{I})$ is a closing by reconstruction of an image I with a structuring element equal to a disk of radius λ (being λ an integer number). The operations in this equation are performed pixel-wise. The result is multiplied by 100 to work in percentages. The granulometric profile of a pixel p can be written as:

G(p) = [gc₁(p), gc₂(p), …, gc_n(p)]

(2)

The size of the G vector is determined by the number of granulometric levels (n), this parameter is chosen depending on the size of the image structures.

The next step of the method proposed consists in classifying the pixels by the k-means method. The resulting image is a k-classes image, where each class corresponds to a set of pixels having a similar granulometric profile (Figure 1).

2.3. k-Means Classification

The k-means is a clustering algorithm aiming at partitioning a data set into a set of meaningful subsets. The k-means algorithm fixes a set of k centroids that minimizes the squared distance from each pixel to its nearest center.

Firstly, the k initial centroids are seeded (k is the final number of classes or clusters chosen) in the multidimensional space (in this study the number of dimensions is equal to the number of levels of granulometry used, or the G vector size). The program uses an improved seeding more efficient than the conventional random seeding (Matlab documentation), the first centroid is chosen randomly, then the next centroids are iteratively chosen with a probability proportional to their distance from the already chosen centroids [23]. Secondly, all the distances from each pixel to each centroid are computed and each pixel is assigned to the closest centroid. Lastly, the average of pixels of each cluster is computed to obtain a new centroid location. Then the second and third steps are repeated until there are no more migrations of pixels to a different centroid, assuring a nearly equal classification from the same input at each test. A nested k-means classification is applied when the overlap of classes is obvious after the first classification.

Although several methods to choose an appropriate number of k have been proposed [24], it seems however difficult to do this without human intervention given that the number of classes depends on the image structures and the human image interpretation cannot be easily substituted. As it will be seen image characteristics such as the contrast or sharpness will have an impact on the choice of k.

3. Results

3.1. Guileyni

Guileyni is a village (380 inhabitants in 2001) located 60 km to the east of Niamey. It is located in the lowest part of a fossil valley. The largest trees are the Faedherbia albida: found around the village, these trees are important for the population as they improve the yield. Smaller trees and shrubs are scattered all over the landscape (Figure 2). The area surrounding the village is occupied by crops and fallows.

The most recent image of the Guileyni area used in this test (January 2011) is a screen capture from a very high-resolution image provided by Digital Globe. It has been resampled to match the spatial resolution of the 1975 image (Figure 3). The image used to perform the diachronic analysis is an aerial photograph taken in a mission of the Institut Géographique National (IGN) in March 1975. The image has been scanned at 1200 dots per inch (dpi) and the spatial resolution of the resulting image is circa 1.3 m. Since the two images are from the dry season, there are no significant phenological differences.

The 2011 image presents significantly higher contrast (0.6357, calculated simply as the normalized difference between the 1% brightest pixels and the 1% darkest pixels) and sharpness (14.0493, calculated as the sum of all gradient norms divided by the number of pixels) than the 1975 image (0.2663 and 2.0253 respectively). A higher contrast and sharpness lead to easier discrimination of classes. The contrast enhancement and histogram matching of the 1975 image did not work well as it made it more complicated to discriminate the classes. In addition, both images present illumination problems: it is known that shadows can artificially increase the size of trees and behave as outliers in a k-means classification. Therefore, the images have been corrected at particularly low values to avoid outliers. The minimum pixel value was empirically set to 50 for all the images presented in this work. When taking into account shadows the exact surface of trees cannot be estimated but the impact on the estimation of the number of trees and the variety of trees sizes is much lower.

To perform the 2011 image classification, 12 levels of granulometry have been used. This choice has been determined by the image resolution and the size of the elements of interest. These are the bigger trees around the village and the smaller trees and shrubs. The largest crown of some grouped trees can reach a diameter superior to 30 m (more than 20 pixels) therefore 12 levels of granulometry is a good choice to detect them without adding useless information.

To perform a k-means classification in 3 classes, 12 granulometric levels have been employed. The result shown in Figure 4 presents suitable ligneous cover detection and a net discrimination between smaller and bigger trees. The third class corresponding to the soil or background has been made transparent to facilitate the visualization. This method allows distinguishing trees from the darker soils (fallows), which could be hardly achieved applying a thresholding technique. A little pond to the north of the village has not been classed as blue since its spectral response slightly differs from that of the trees.

The graphic in Figure 5 displays the centroid locations of the 3 classes depending on the size of the SE. Good discrimination between classes has been observed, being the overlapping of the classes minimal. The blue class shows a peak value for a SE of size 5, corresponding to the medium size of trees detected by this class. The green class presents a maximum for a SE of size 1 and 2, accounting for the smallest trees. This maximum is inferior to that of the blue class because the spectral response of the smallest trees can be easily merged with that of the background, thus increasing the mean value of these pixels (and reducing the difference between the original image and the reconstructed image). The red class (corresponding to the soil background in the image) logically shows a response close to zero for any size of the SE.

A classification in 7 classes using 20 granulometric levels would display an obvious class overlapping as shown in Figure 6. Blue and green classes or cyan and red classes account for objects of the same size but different spectral response. The red class displays a maximum for a SE of size 5. The objects detected by this class are brighter than the largest trees spotted by the cyan class, corresponding mainly to dark soils (fallows). As observed in the Figure 6 little information would be added by using granulometric levels superior to 12.

The 1975 image classification has been performed using 4 classes and 12 levels of granulometry. As this image has a significantly lower contrast and sharpness than the 2011 image, a higher k was needed for a better discrimination of classes. The number of k was chosen empirically: a k inferior to 4 would have led to an important intra-class confusion while a k superior to 4 would have led to a noteworthy inter-classes overlapping. A k equal to 4 was found to be a good compromise for the discrimination of classes, even though some problems remained.

While the number of granulometric levels (SE size) is easily selected according to the size of the structure of the images, the problem of systematically choosing an appropriate k remains. The correct k relies on the structure of the images and on the human interpretation of them. The elbow method [25] to determine k has been used during this study but was inconclusive.

Figure 7 displays the 1975 classification in 4 classes. The low contrast and sharpness leads to maxima values for larger SE sizes. Overlapping exists between the magenta and the green classes and between the magenta and the blue ones (Figure 8). The blue and magenta classes depict the largest trees while the green one accounts essentially for smaller trees and dark spots linked to soil conditions, which have brighter grey values than the trees.

To improve the classification a k-means classification in 3 classes was performed over the green class. This nested classification allows to distinguish the smaller trees of the green class from the first classification. As shown in Figure 9, the nested classification for the green class displays a maxima for SE of size 5, depicting a higher centroid value even though the corresponding SE size is smaller than the one for the other two classes. The brown and black classes from the nested classification has been eliminated as they are related to dark spots not linked to the ligneous cover. The final classification is displayed in Figure 10. The total number of smallest trees and shrubs is highly underestimated in this classification as compared to 2011 classification, being this directly linked to the image characteristics (low contrast and sharpness). The classification is mostly satisfactory leading to 3 properly differentiated classes.

The comparison between the two final classifications (1975 and 2011) shows some changes during the last decades. A higher number of final classes could suggest a higher ligneous diversity in 1975. This result agrees with the observed decrease in ligneous diversity in this region during the last decades. A lower diversity does not imply less ligneous cover: in fact the analysis shows a higher presence of large and small trees in 2011. This could be explained by a sustainable management of the local population and a better integration of ligneous cover and agriculture. The improved species selection by the peasants and the 80’s droughts could explain a declining diversity.

The 2011 image classification was compared to ground truth sites for accuracy assessment. The ground truth sites were digitized manually by visual interpretation of the 2011 image using the original Google Earth image for verification (Figure 11). Only the small trees and shrubs that were distinguishable (after verification) in the 2011 resampled image were digitized. Three sites are representative of croplands, the bottom of the valley and the village were chosen for the accuracy assessment (Figure 12). Unfortunately, the same evaluation cannot be performed with the 1975 classification. In general, the classification shows lower level of errors of commission in tree detection while the level of errors of omission is higher (

Table 1. Error matrix analysis for three sites (on a pixel basis). Ground truth in columns, classification in rows. Guileyni.

1. CROPS	Background	Small trees	Big trees	ErrorC
Background	33,535	152	3	0.0046
Small trees	13	1802	18	0.0169
Big trees	0	1	181	0.0055
ErrorO	0.004	0.0783	0.1040
2. VALLEY	Background	Small trees	Big trees	ErrorC
Background	31,483	77	13	0.0029
Small trees	16	2105	0	0.0075
Big trees	0	10	2380	0.0042
ErrorO	0.0005	0.0397	0.0054
3. VILLAGE	Background	Small trees	Big trees	ErrorC
Background	30,385	3	57	0.0020
Small trees	1224	567	2	0.6838
Big trees	711	0	3329	0.1760
ErrorO	0.0599	0.0053	0.0174

ErrorO: error of omission. ErrorC: error of commission. The errors are expressed as proportions.

). This is particularly true for the croplands and the bottom of the valley sites; at these sites, the omission error for the class representing the small trees is particularly high (0.0783 and 0.0397 respectively). Some small isolated trees and shrubs were omitted by the classification, as a consequence of the image resolution used for the classification (1.3 m). In the small fraction of the image occupied by the village the commission errors are much higher than omission errors, this is due to shadows and dark houses that are classified as small or big trees.

As pointed out in [26], trees occupying a surface close to or inferior to the pixel surface (around 2 m²) can hardly be detected. This is an unavoidable drawback of the spatial resolution used. It is assumed that the same problem is present in the 1975 classification, thus underestimating the presence of isolated ligneous cover.

3.2. Kirib Kaina

Kirib Kaina (1068 habitants in 2001) is located 75 km to the northeast of Niamey. This village, located in the lower-right corner of the image in Figure 13, is close to a fossil valley and a degraded plateau that has been over-exploited during the last decades for fodder supply (Figure 14). In these photos, the tiger bush is characterized by banded ensembles of trees and shrubs following topoedaphic conditions. The images show that the tiger bush has nearly completely disappeared between 1975 and 2010 leading to bare soil. The changes in land use have triggered the degradation of the hydrological system [27]. Isolated trees, which were randomly scattered throughout the area in 1975, are generally found to be organized linearly following a system of gullies in 2010. The hydrological degradation seems to favor an increase of spatial heterogeneity and the fragmentation of landscape.

The same procedure as for the Guileyni village has been used. The most recent image of Kirib Kaina used in this test (January 2010) is a screen capture from a very high-resolution image provided by Digital Globe. This image was resampled to the resolution of a scanned aerial photograph from an IGN mission in March 1975. The resulting spatial resolution is circa 1.3 m.

As in the case of Guileyni, the most recent image presents higher contrast (0.6088) and sharpness (7.5778) than the 1975 image (0.3172 and 2.557 respectively). A high contrast and sharpness increases the discrimination capacity of the k-means method. Contrast enhancement has not been applied over the 1975 image as it would have increased the confusion between classes, particularly between the limits of the plateau (Figure 14) or the more humid soils, which are darker, and the tiger bush. Up to 25 levels of granulometry have been used for the classification of both images. The number of levels has been chosen according to the size of the elements characterizing the ligneous cover. Two opposed categories have been differentiated: the isolated trees and the tiger bush. As the tiger bush is at some places up to 60–70 m width (more than 45 pixels), a 25 level granulometry has been chosen.

The 2010 image classification was performed using the 25 levels resulting from the morphological analysis (Figure 15b). An automatic nested classification was tested, the final result is similar to the equivalent k-means classification. Firstly, a 4 classes k-means classification was used to automatically separate the ligneous cover from the background. Then the program detected the class with higher mean in the original image which corresponded to the background. Secondly, a 3 classes k-means was automatically performed on the pixels classified as ligneous cover to discriminate between the isolated trees and the tiger bush. The 3 final classes showed a good discrimination with no overlapping between classes (Figure 16). The blue class accounts for isolated shrubs and trees while the red class accounts for bigger trees with a peak for a SE of size 5. The green class is identified as the tiger bush and its granulometric profile is completely different from the other two classes. It shows several local maxima for a SE of size 13, 16, 21 and 23. The different structures existing in the image are well discriminated by the granulometric profiles.

The 1975 image has lower contrast and sharpness and therefore the confusion between classes is more likely to happen. The number of k and the number of nested classifications must be higher because the discrimination of classes is more complex. There is a possible confusion between the tiger bush and the darker soils (fallows or natural savanna) and between the tiger bush and the slopes connecting the laterite plateau and the valley (dark, rocky soils). The tiger bush and the isolated trees are detected separately. A 6 k-means classification has been performed by using 25 granulometric levels. The granulometric profile shows an important overlap between classes (Figure 17): For example, the magenta, the orange and the green class are all linked to the tiger bush and at the same time they are linked to the plateau slopes. Successive nested classifications were applied to obtain the result shown in Figure 15a. There is still some confusion particularly between some dark soils and the tiger bush. The comparison between the two sets of data highlights the important ecological changes which have taken place during the last decades. The tiger bush has almost completely disappeared and the distribution of isolated trees has been modified.

To perform the accuracy assessment, two ground truth sites were digitized by visual interpretation of the 2010 image, as for Guileyni the original Google Earth image was used for verification (Figure 18). These two sites are located on the most significant landscape units existing in the image: the plateau and the drainage network associated to the plateau (Figure 19). Errors of omission are particularly high for small trees (

Table 2. Error matrix analysis for two sites (on a pixel basis). Ground truth in columns, classification in rows. Kirib Kaina.

PLATEAU	Big Trees	Tiger Bush	Small Trees	Background	ErrorC
Big trees	656	0	10	0	0.015
Tiger bush	0	8005	1	0	0.0001
Small trees	15	49	1069	0	0.0565
Background	0	29	108	61524	0.0022
ErrorO	0.0224	0.0096	0.1002	0.0000
DN	Big trees	Small trees	Background	ErrorC
Medium	2996	0	26	0.0086
Small trees	93	1837	28	0.0618
Background	81	107	66575	0.028
ErrorO	0.0549	0.0550	0.0008

DN: Drainage Network. ErrorO: error of omission. ErrorC: error of commission. The errors are expressed as proportions.

), which is directly linked to the spatial resolution of the image. Some confusion between classes is also significant, mainly the confusion between small trees and the tiger bush for the first site (tiger bush being classed as small trees) and between small trees and big trees for the second site.

4. Discussion and Conclusions

The results highlight that mathematical morphology is suitable for detecting and classifying ligneous cover in semi-arid lands. Mathematical morphology has been rarely used previously in ligneous cover detection. The method presented in this paper allows identifying ligneous cover while at the same time it classifies the vegetation patches according to their size following a mathematical logic linked to the richness and variety of the ligneous cover. This last characteristic is important to understand and quantify the organization and the dynamics of ligneous cover.

Results are particularly satisfactory when recent images that have both better contrast and sharpness are used. For these, the number of k classes of the final classification can be chosen intuitively from the number of types of structures perceived in the image. There is usually no need for a nested classification and the method is more efficient and more accurate. The images from aerial photographs of 1975 have lower contrast and sharpness for the same spatial resolution. In these images, the overlapping between different classes is frequently evident and nested classifications are required. The chosen k must then search a compromise between a high intra-class confusion (lower k) and a high inter-class confusion (higher k). If the appropriate number of granulometric levels can be easily inferred from the size of the structure of the images, the choice of the k number of classes is more complicated. This is particularly true when the global image quality is poor. Then the choice of a k relies on human interpretation and it can be rarely automated.

The assessment of the classification accuracy from comparison to a ground truth image shows a lower level of errors of commission while the level of errors of omission are more prominent. The errors of omission (principally small trees and shrubs that remain undetected) are linked to the image spatial resolution and are unavoidable. Even if there exists some confusion between classes the classifications of recent images are satisfactory. The method is not suitable for urban areas because of the impact of buildings and their shadows on the classification. The accuracy of the 1975 classifications cannot be measured, given that there are not higher resolution images and there is no field data. However, it can be assumed that errors of omission for small structures are higher for the 1975 classification.

The diachronic analysis between 1975 and recent years shows the important changes taking place during the last decades. Ligneous cover is the best witness of these changes as it is deeply linked to socio-environmental conditions. The mapping of ligneous cover in k classes accounting for different patches sizes is important for a better understanding of the dynamics observed.