Research Article
Extraction of Watersheds from Digital Elevation Models Using Mathematical Morphology
Science and Technology Research Institute for Defence (STRIDE), Ministry of Defence, KD Malaya, Pangkalan TLDM, 32100 Lumut, Perak, Malaysia
INTRODUCTION
A catchment is the topographic area from which all water runoffs finally reach one single given point, known as the pit. As shown in Fig. 1, watersheds refer to the topographic barriers that divide catchments from each other (Monkhouse, 1965; Band, 1986; Rodríguez-Iturbe and Rinaldo, 1997; Subramanya, 2006). Watersheds are important geomorphological features which play an important role in hydrological modeling. Many hydrogeological processes, such as soil erosion, mass movements, sediment transport and land cover changes, are strongly linked to this spatial reference unit.
Traditionally, watersheds were derived manually from topographic maps, which was a labor intensive activity. In recent times, extraction techniques have evolved from manual through computer assisted to automated methods; with Digital Elevation Models (DEMs) as the input data. DEMs are a popular source for hydrological modeling and watershed characterization because of their simple data structure and widespread availability and they lend themselves too many GIS processes and operations. In seeking the efficient extraction of watersheds from DEMs, various algorithms have been proposed. Extensive reviews of these algorithms are provided in Douglas (1986), Moore et al. (1991), Tribe (1992), Bertolo (2000) and Pike (2002). These algorithms can be classified into two categories; local morphology based and flow directional based.
Fig. 1: | The analogy of a watershed |
In the local morphology approach, points belonging to hydrological features are defined from local morphologies. For example, watershed pixels are considered as pixels that have a convex curvature coefficient that is higher than a given threshold (Johnston and Rosenfeld, 1975; Peucker and Douglas, 1975; Torikawa and Fukumara, 1978; Band, 1986; Soille and Ansoult, 1990; Vincent and Soille, 1991; Chockalingam and Sagar, 2003; Dinesh, 2007b). The disadvantage of this approach is there that is no absolute threshold to discriminate watershed pixels from other pixels, resulting errors in extracted watersheds.
In the flow directional approach (Mark, 1983, 1984; O`Callaghan and Mark, 1984; Marks et al., 1984; Jenson and Domingue, 1988; Fairfield and Laymarie, 1991; Freeman, 1991; Wichel et al., 1992; Benosky and Merry, 1995; Mackay and Band, 1998), a flow of water is simulated over the topographic surface. A commonly used watershed extraction algorithm is the flow direction based algorithm proposed in Jenson and Domingue (1988). Flow direction between neighboring pixels is described by a single direction according to the steepest downward slope in a size 3 square neighborhood. From the flow direction grid, a flow accumulation grid can be calculated and from this grid, the watersheds are be extracted. The disadvantage of this approach is that it is unable to operate effectively on flat areas (areas containing pixels of equal elevation) in DEMs, where flow routing is unable to occur. A number of drainage enforcement algorithms have been developed to determine flow direction in flat areas of DEMs (Jenson and Domingue, 1988; Hutchinson, 1989; Tribe, 1992; Soille and Gratin, 1994; Saunders and Maidment, 1995; Garbrecht and Martz, 1997; Martz and Garbrecht, 1998; Mackay and Band, 1998; Turcotte et al., 2001; Soille et al., 2003). However, drainage enforcement is a computationally expensive operation that can cause the formation of artifacts such as double streams and meanders (Soille et al., 2003).
In order to rectify these shortfalls, this study proposes a mathematical morphological based algorithm to extract watersheds from DEMs. Mathematical morphology is a branch of image processing that deals with the extraction of image components that are useful for representational and descriptional purposes. Mathematical morphology has a well developed mathematical structure that is based on set theoretic concepts. The effects of the basic morphological operations can be given simple and intuitive interpretations using geometric terms of shape, size and location. Mathematical morphology is well suited to the processing of elevation data because in morphology, any image is viewed as a topographic surface, the grey level of a pixel standing for its elevation (Soille and Ansoult, 1990). The fundamental morphological operators are discussed in Matheron (1975), Serra (1982) and Soille (2003). Morphological operators generally require two inputs; the input image A, which can be in binary or grayscale form and the kernel B, which is used to determine the precise effect of the operator. Each pixel in A is compared with B by moving B so that its center hits the pixel. Depending on the type of morphological operator employed, the pixel value is reset to the value or average value of one or more of its neighbors.
While the proposed mathematical morphological based watershed extraction algorithm is applicable to any DEM, in this study, the effectiveness of the proposed algorithm is tested by implementing it on the DEM in Fig. 2, which shows the area of Great Basin, Nevada, USA. The area is bounded by latitude 38° 15` to 42° N and longitude 118° 30` to 115° 30`W. The DEM was rectified and resampled to 925 m in both x and y directions. The DEM is a Global Digital Elevation Model (GTOPO30) and was downloaded from the USGS GTOPO30 website (http://edcwww.cr.usgs.gov/landdaac/gtopo30/gtopo30. html). GTOPO30 DEMs are available at a global scale, providing a digital representation of the Earth`s surface at a 30 arc-seconds sampling interval. The land data used to derive GTOPO30 DEMs are obtained from Digital Terrain Elevation Data (DTED), the 1-degree DEM for USA and the digital chart of the world (DCW). The accuracy of GTOPO30 DEMs varies by location according to the source data. The DTED and the 1-degree dataset have a vertical accuracy of ±30 m while the absolute accuracy of the DCW vector dataset is ±2000 m horizontal error and ±650 vertical error (Miliaresis and Argialas, 2002).
Fig. 2: | The GTOPO30 DEM of Great Basin. The elevation values of the terrain (minimum 1005 m and maximum 3651 m) are rescaled to the interval of 0 to 255 (the brightest pixel has the highest elevation). The scale is approximately 1:3,900,00 |
THE PROPOSED WATERSHED EXTRACTION ALGORITHM
DEM smoothening using morphological smoothening by reconstruction: The peaks of a terrain refer to the highest points of the mountains of the terrain while the pits of the terrain are the lowest points of the basins of the terrain. In Digital Elevation Models (DEMs), peaks are connected components that are completely surrounded by pixels of lower elevation while pits are connected components that are completely surrounded by pixels of higher elevation. The detection of peaks and pits is the first step in most techniques used to perform DEM characterization and to describe the general geomorphometry of a surface.
DEMs contain two types of peaks and pits; real peaks and pits and spurious peaks and pits. Spurious peaks and pits are errors cause by input data error, interpolation procedures and the limited horizontal and vertical resolution of DEMs. Spurious peaks and pits do not correspond to real landscape features and they cause errors in features extracted from DEMs. Hence, the removal of spurious peaks and pits from DEMs, known as DEM smoothening, is an important preprocessing step in DEM analysis. A number of smoothening operators have been employed to implement DEM smoothening, including mean smoothing (Mark, 1983, 1984; O`Callaghan and Mark, 1984), Gaussian blurring (Seemuller, 1989; Lee et al., 1992), Kalman filtering (Wang, 1998), the á Trous wavelet transform (John et al., 2002) and loess regression (Tate et al., 2005). However, the conventional DEM smoothening operators have several disadvantages:
• | When the kernel straddles an edge, it will interpolate new values for pixels on the edge and so will blur that edge (Fisher et al., 1994). In DEM processing, edges represent crenulations, which are used to extract hydrological features from DEMs. Convex crenulations are used to extract ridge networks while concave crenulations are used to extract drainage networks (Gilbert, 1909; Howard, 1994; Rodríguez-Iturbe and Rinaldo, 1997). The blurring of edges causes the unnecessary loss of information in the form of crenulations and hence, results in errors in extracted hydrological features. |
• | They operate with the assumption that DEMs contain equal number of spurious peaks and pits. However, since the noise (mean error in elevation) to signal (elevation) ratio in DEMs in higher in the lower elevation regions than in the higher elevation regions, DEMs actually contain more spurious pits than spurious peaks, Hence, conventional DEM smoothening methods cause the unnecessary loss of information in the form of real peaks. |
• | They alter all elevation values, including the pixels which are not marked as spurious peaks or pits. |
Dinesh (2006, 2007a) proposes a mathematical morphological based DEM smoothening algorithm, known as morphological smoothening by reconstruction, that is able to rectify these shortfalls.
Dilation sets the pixel values within the kernel to the maximum value of the pixel neighbourhood. The dilation operation is expressed as:
(1) |
Erosion sets the pixels values within the kernel to the minimum value of the kernel. Erosion is the dual operator of dilation:
(2) |
An opening (Eq. 3) is defined as an erosion followed by a dilation using the same kernel for both operations. A closing (Eq. 4) is defined as a dilation followed by an erosion using the same kernel for both operations.
(3) |
(4) |
Grayscale opening is used to darken small bright areas and to reduce sharp peaks in images, while grayscale closing is used to brighten small dark areas and to fill valleys in images. In DEM processing, opening is used to remove the spurious peaks and closing is used to remove the spurious pits. Morphological smoothening is implemented by performing opening followed by closing. The net result of implementing morphological smoothening on a DEM is that the spurious peaks and pits in the DEM are removed.
The biggest disadvantage of morphological smoothening is that the closing operation in the smoothening process can create big pits by filling the narrow valleys while the opening operation can create big peaks by enlarging valleys.
Morphological reconstruction allows for the isolation of certain features within an image based on the manipulation of a mask image X and a marker image Y. It is founded on the concept of geodesic transformations, where dilations or erosion of a marker image are performed until stability is achieved (represented by a mask image) (Vincent, 1993).
The geodesic dilation δG used in the reconstruction process is performed through iteration of elementary geodesic dilations δ (1) until stability is achieved.
(5) |
The elementary dilation process is performed using standard dilation of size one followed by an intersection.
(6) |
The operation in Eq. 6 is used for elementary dilation in binary reconstruction. In grayscale reconstruction, the intersection in the equation is replaced with a pointwise minimum (1993).
Morphological reconstruction can be used to maintain the spurious peak and pit removal effect of opening and closing, while avoiding their valley filling and enlargement effects (Vincent, 1993). Morphological smoothening by reconstruction (Dinesh, 2006, 2007a) is implemented using the following process:
Step 1: | Performing opening on the original image: This step removes the spurious peaks of the DEM. However, the valley enlargement effect of opening causes the formation of larger peaks. |
Step 2: | Performing morphological reconstruction on the opened image: Morphological reconstruction is performed using the opened image as the marker and the original image as the mask. This step is performed to avoid the valley enlargement caused by the opening operation. |
Step 3: | Performing closing on the opened reconstructed image: This step removes the spurious pits of the DEM. However, the valley narrowing effect of closing causes the formation of larger pits. |
Step 4: | Performing morphological reconstruction on the closed image: Morphological reconstruction is performed using the closed DEM as the marker and the opened reconstructed image as the mask. This step is performed to avoid the valley filling effect caused by the closing operation. The resultant DEM obtained from this process is known as a smoothened DEM. |
Morphological smoothening by reconstruction is implemented on the DEM of Great Basin using a disk kernel of size 3 to obtain the smoothened DEM (Fig. 3a). A total of 12,846 pixels (0.14%) are modified by the smoothening process (Fig. 3b).
The effectiveness of DEM smoothening using morphological smoothening by reconstruction is tested by performing drainage network extraction on the original and smoothened DEMs of Great Basin. In this study, drainage network extraction is performed using the drainage skeletonization algorithm proposed in Meisels et al. (1995). The algorithm extracts pixels that lie in high curvature contours starting from pixels of maximal elevation, elevation-level by elevation-level; the selection is based on a condition for a large enough number of higher elevation pixels in the immediate neighborhood of a pixel belonging to the elevation currently being processed. Complementary local conditions of connectivity are then used to connect all the pixels of the flow path.
As shown in Fig. 4a, the spurious peaks in the original DEM cause the formation of closed loops in the extracted drainage network. The spurious pits in the DEM cause a portion of the extracted the drainage network to be incomplete and disconnected. Drainage network extraction applied on the smoothened DEM (Fig. 4b) allows for the extraction of drainage networks that are loopless, complete and connected.
Pit extraction: Grayscale erosion can be used to remove bright areas in grayscale images. It causes small peaks in the image to disappear. However, it also causes valley widening which results in larger peaks. Morphological reconstruction can be used to maintain the peak removal effect of erosion while avoiding its the valley enlargement effect (Vincent, 1993). The peaks removed by erosion can be obtained by subtracting the reconstructed eroded image from the original image.
Fig. 3: | DEM smoothening using morphological smoothening by reconstruction (a) The smoothened DEM and (b) The mask of pixels modified by the smoothening process |
Fig. 4: | Implementation of drainage network extraction on the original and smoothened DEMs of Great Basin (a) The drainage network of the original DEM and (b) The drainage network of the smoothened DEM |
Fig. 5: | An example of the ultimate erosion operation (Source: Duchane and Lewis, 1996) |
In order to extract the peaks of a DEM, ultimate erosion is performed on the DEM. Ultimate erosion is implemented through the iterative erosion of the image until all objects vanish (images Xi) and the reconstruction of each eroded image using the eroded image, E(Xi), as the mask and the erosion of smaller size as the marker. The reconstructed images (images Yi) are subtracted from the corresponding eroded images to form the eroded sets (images Ui) (Duchane and Lewis, 1996). Figure 5 demonstrates the operation of ultimate erosion. The generated ultimate eroded set of the DEM forms the peaks of the DEM. The pits of the DEM are the peaks of the inverted DEM; pit extraction is implemented by performing ultimate erosion on the inverted DEM.
The pits of the smoothened DEM are extracted by implementing ultimate erosion on the inverted DEM. A total of 236 pits (Fig. 6) are extracted from the smoothened DEM. A total of 590 pixels (0.65%) are identified as pit pixels.
Watershed extraction: A skeleton by influence zone (SKIZ) is a skeletal structure that divides an image into regions, each of which contains just one of the distinct objects in the image. The boundaries are drawn such that all points within a particular boundary are closer to the binary object contained within that boundary than to any other (Fisher et al., 1994).
The catchments of a DEM are the influence zones of the pits of the DEM (Soille and Ansoult, 1990). Hence, the watershed enclosing a catchment can be viewed as the skeleton by influence zone of the pit of the drainage basin. Watershed extraction is performed by implementing the skeletonization by influence zone algorithm proposed in Fisher et al. (1994) on the dilated pits. First, skeletonization by morphological thinning is performed on the inverted image. Skeletonization by morphological thinning is defined as successive removal of outer layers of pixels from an object while retaining any pixels whose removal would alter the connectivity or shorten the legs of the skeleton. The process is converged or completed when no further pixels can be removed without altering the connectivity or shortening the skeletal legs (Jang and Chin, 1990). The barbs of the skeleton are then pruned off to obtain the watersheds.
A shown in Fig. 7a, a number of the extracted pits are too small for the purposes of segmentation. A number of extracted pits are too small for the purposes of segmentation. In order to smoothen these irregularities, the pits are dilated with a size 3 square kernel. Skeletonization by influence zone implemented on the dilated pits (dilated with a size 3 square kernel) allows for a proper portioning of the DEM (Fig. 7b).
Fig. 6: | A total of 236 pits are extracted from the smoothened DEM of Great Basin |
Fig. 7: | Skeletonization of the influence zone applied to the pits and dilated pits of the smoothened DEM of Great Basin (a) The skeleton by influence zone of the pits and (b) The skeleton by influence zone of the dilated pits. This forms the watersheds of the DEM |
Fig. 8: | Identification of the extracted catchments. The catchment count number is determined by the grey level; the brighter the grey level, the larger the count number |
Fig. 9: | Catchments extracted using the flow direction based algorithm proposed in Jenson and Domingue (1988) |
The catchments of the DEM are identified by performing the connected component labeling algorithm proposed in Pitas (1993) on the inverted watershed image. As shown in Fig. 8, the extracted watersheds partition the DEM into 167 catchments.
Evaluation of results: The results obtained using the flow directional based algorithm proposed in Jenson and Domingue (1988) algorithm is shown in Fig. 9. A total of 177 catchments are extracted from the DEM.
A good match was evident between the results obtained using proposed mathematical morphological based algorithm and the results obtained using Jenson`s flow direction algorithm. The disadvantage of Jenson`s algorithm is that flow direction algorithms are unable to operate effectively in flat areas of the DEM, resulting in the formation of spurious catchments. The proposed morphological based algorithm, which does not rely on flow direction algorithms, is able to operate effectively on flat areas in DEMs.
CONCLUSION
In this study, a mathematical morphological based algorithm to extract watersheds from DEMs was developed. In general, the proposed algorithm consists of 3 stages, namely DEM smoothening using morphological smoothening by reconstruction, pit extraction using ultimate erosion and skeletonization by influence zone of the dilated pits.
The advantage DEM smoothening using morphological smoothening by reconstruction is that it provides the user with the flexibility to choose the threshold level of spurious peak and pit removal by varying the number of opening and closing steps in each smoothening step and the size and type of kernel used. The reconstruction steps in the algorithm allow for the preservation of crenulation information in DEMs and ensure that only pixels marked as spurious peaks and pits are altered.
The disadvantage of many available pit extraction algorithms is that they are only able to extract shallow pits but fail to extract the deeper pits. The proposed pit extraction algorithm is able to solve this problem as it performs the ultimate erosion operator until convergence; allowing for all of the pits in the DEM to be extracted.
In general, the proposed watershed extraction algorithm, summarized in Fig. 10, is able to operate effectively on flat areas in DEMs and allows for the fast and effective segmentation of DEMs into their constituent catchments.
Apart from GIS applications, watersheds play an important role in the segmentation of raster images, because the grey level of an image can be regarded as the height of the terrain and many interesting features, such as cell walls in microscopic images, can be described using watersheds. The proposed watershed extraction algorithm is able to contribute to new insights into the use of watersheds in the field of image analysis and allows a user to resort to watershed extraction for solving complex segmentation problems. More experiments are currently being conducted to evaluate the feasibility of the proposed watershed extraction algorithm to perform the segmentation of various medical images.
Fig. 10: | The proposed watershed extraction algorithm |