GIS/EM4 - An Advanced Shape-Of-Country Classifier:

An Advanced Shape-Of-Country Classifier:

Extraction of surface features from DEMs

GIS/EM4 No. 205

Lee Herrington
Gerald Pellegrini

Abstract

A number of environmental analyses require the classification of digital elevation models (DEMs) into shape-of-country features. The 10 terrain feature classes used here are ridge, valley, peak, pit, saddle, saddle slope, hill slope, concave slope, convex slope, and flat. These feature classifications are used in hydrological, ecological and other studies. Our specific use is for the development of knowledge based soil property estimators, in particular trafficability analyses, and for the refinement of limited soils information. The technique described here, a combination of surface smoothing in the spatial frequency domain and the determination of the feature class of each cell through the fitting of a mathematical surface to a 3x3 roving window passed over a grid results in 100% of the grids being classified and a classification accuracy of approximately 90%. The Fast Fourier Transform (FFT) is used to produce an image of the DEM in the spatial frequency domain. A low-pass filter with a square wave roll-off in the spatial frequency domain is used to remove high spatial frequency information. This is followed by transformation back into the space domain. The resulting smoothed surface is analyzed by passing a 3x3 roving window over the grid. Eigenvectors, containing all the information needed to represent the feature class of the window, are generated using the Jacobi method from a polynomial representation of the surface in the window. Obtaining the best estimators of the feature class requires that surface cells be classified according to the most significant feature occurring anywhere within the central cell of the window. For example, ridges are more significant than slopes. What we have added to this type of analysis is the determination of the significant features that may occur anywhere in the central cell. Heuristics are then used to refine the classification. The results of applying the classifier to several mathematical surfaces and to real DEMs from Central New York are presented.

Keywords

Shape of Country, Terrain Analysis, Spatial frequency domain, Filters, Digital elevation models, DEM, Smoothing.

Introduction

Terrain Analysis combines the qualitative and quantitative descriptions of the patterns and structure of hypsographic surfaces (Mitchel, 1973; Way, 1973). These surfaces are often represented by a digital elevation model (DEM). The usual way to model the surface is to create a regular square tessellation or grid wherein each cell in the grid has the value of the elevation of the underlying point on the landscape. One result of terrain analysis is a map or image of the DEM's "Shape-of-Country". This term applies to a second grid, derived from the DEM, in which each cell contains a value representing the nature of the underlying surface. Thus instead of elevation values the Shape-of-Country grid contains nominal values representing the nature of the surface. Is the surface horizontal? A concave slope? A saddle, a peak, a pit, etc?

The nature of the surface, its shape, is important to many environmental and ecological processes and problems. For example, the floristic composition of ecosystems can be related to these measures. In forest ecosystems, foresters use site index curves to estimate the productivity of the forest. These species dependent curves are dependent on soil, climate, and topography (Spurr and Barnes, 1980). The hypsography can be broken down into measures of slope, aspect, and location on the slope (top, mid slope, bottom). These elements do not act independently: they all influence climate and soil, for example. Hydrologic modeling provides another example of the importance of shape-of-country since the nature and shape of the surface determine the flow of water over the landscape (Wang and Yin, 1997).

The problem which created the need for this work was the need to be able to estimate physical soil properties in areas where no direct measurements or even estimates existed. The specific area of interest was the Adirondack Mountains of upstate New York. The lands in the Adirondack Preserve are "forever wild" and thus there has been little need for soils information. One of our colleagues, Dr. Phillip Craul, worked with a Defense Department Contractor to develop a technique to estimate soil tractability based on knowledge of shape-of-country and satellite derived cover information. The contractor built an expert system that had a high-level of classification accuracy using Dr. Craul's knowledge of the soils relationship to the shape-of-country. The system they used to derive shape-of-country was, unfortunately, classified so we had to reinvent it. Thus this project was born.

Problem statement

The problem is quite simply stated: develop a technique to accurately convert a DEM into an grid where each and every cell has a nominal value representing the shape-of-country classifications shown in Table 1.

Class Descritpion Class Description

Peak A surface location that is higher in elevation than all surrounding cells Concave Hillside Hillside that is concave upward

Pit A surface location that is lower than all surrounding cells Convex hillside Hillside that is convex upward

Ridge Contiguous lines of cells that occur when the cells are higher than the cells on either side Saddle hillside Hillside that has positive curvature in one direction and negative curature in the other, orthogonal direction

Valley Contiguous lines of cells that occur when the cells are lower than the cells on either sided Flat hillside Hillside with no curvature

Saddle point Intersection of ridge and valley Flat Surface with no slope

Table 1. Classification categories for shape-of-country.

All cells in a DEM representing terrain can be put into one of these 10 classes. The general problem was reduced to two sub-problems: smoothing the DEM and classifying the image. Shape-of-country images should represent broad features of the terrain like hillsides, valleys, etc. and not the roughly 2x cell dimension variations in elevation that result from error in the DEM(s) and/or the naturally occurring shortwave variations in elevation between nearby cells. Thus a technique had to be found that would smooth the surface so that the resulting shape-of-country image would represent these broad features of the terrain. There are several techniques for classifying DEMs into shape-of-country. The more common approach is to use discontinuous techniques, based on cell to cell variations, to identify terrain features. There appear to be no methods of this type that can effectively generate the 10 classes we required. Therefore the technique had to be based on continuous techniques which consider the grid elevations to be samples that represent a continuous surface.

Background

Smoothing

Extracting topographic classes or tokens requires removing short wavelength variation in terrain (roughness or texture) (Mark, 1975) from the DEM. This is because the roughness will obscure the longer wavelength shape-of-country features (grain) when using a relatively small n x n roving window. In terms of sampling, roughness features with a minimum wavelength Fourier component of 2x the spatial distance between cells cannot be resolved. Grain features larger than the dimensions of overall grid cannot be resolved (Blackman and Tukey, 1958). This means that to resolve grain the short wavelength variations have to be removed from the DEM by filtering or smoothing before processing.

The most common way to smooth a surface to reduce texture is to pass a low-pass filter over the DEM. Low-pass filters are usually a 3x3 cell (or other odd numbered pair) roving window that is passed over the DEM so that a new grid is created in which each cell contains the average of the original values in the 9 cells of the window. A DEM can be thought of as the summation of sinusoids of various wavelengths or spatial frequencies (Fourier Analysis). Since the averaging process does not remove all the short wavelength (high spatial frequency) information the filtered surface does not have these roughness features totally eliminated.

If the Fourier transform is used to transform the spatial data (units:meters) into the spatial frequency domain (units:1/meters) then more specific filters can be made. In the spatial frequency domain the DEM becomes an image where grid values represent a sinusoid of a specific frequency (radial distance) of variation acting in a direction (angle) and with an amplitude that is the value of the cell. In this domain a very sharp cutoff filter can be made by multiplying the DEM by another grid that has a central disk of 1's whose radius is the cutoff frequency. This results in all the high spatial frequency elements of the image being removed. The inverse transform converts the data back into the spatial domain. The resulting DEM will be very smooth. Pellegrini (1995) provides a detailed description of this process.

Mathematical basis for classification

There are basically two methods used to attempt to classify DEMs into shape-of-country information: discontinuous and continuous. Discontinuous methodology ignores the continuous nature of elevation surfaces and bases shape classification on qualitative evaluations of geometric relationships. Flow routing algorithms are an example of this methodology. Although these methods often provide satisfactory results there are a number of problems with their use (Pellegrini,1995).

The continuous methodology is better suited to the recognition of surface structure and is the methodology used here. Continuous methods determine terrain features by fitting a mathematical function to a neighborhood of elevation samples. The neighborhood between sample points is assumed to be continuous so a windowed operator provides a functional approximation to the original continuous surface. Under ideal circumstances the general form of the surface model equation minimizes both prediction error and equation complexity. Slope, aspect, and curvature describe the mathematical relationships between points within a neighborhood in the DEM ( Haralick et al, 1983; Anton, 1988; Larson and Hostetler, 1990). Slope is the first derivative of the surface fitting equation, aspect is the compass direction of maximum slope, and curvature is the second derivative of the surface function. Curvature is evaluated in two orthogonal directions: profile curvature measures how rapidly the surface is changing in the direction of the aspect while plan curvature measures this orthogonal to the aspect. The calculus of surface derivatives provides a rigorous definition of surface structure. Once basic features such as ridgelines have been mathematically proven to exist, then terrain shapes can be reconstructed by comparing the mathematical properties of each pixel to known surface structures. Assembling individual cells into recognizable landforms can then proceed using higher order reasoning (Frank et al, 1986; Lazaroff et al, 1990; Graff and Usery, 1993). Zevenbergen and Thorne's (1987) partial quartic equation, Figure 1, was used to fit a surface to a 3 x 3 roving window.

Figure 1 . Zevenbergen and Thorn's (1987) polynomial used to fit a continuous surface to a 3x3 roving window

This polynomial provides a good fit to the elevation points in the 3x3 window as long as there are no discontinuities (cliffs or overhangs) (Pellegrini, 1995). The equations for determining slope, aspect, and the two curvature measures can be found in Pellegrini (1995) or Zevenbergen and Thorn (1987). The second directional derivative for the central cell of the window is given by Pellegrini (1995) as:

(Notation 1)

From the perspective of surface analysis this derivative provides all the information necessary to correctly classify the shape of the window. It is assumed that the significant changes in surface shape are most easily identified when evaluating the second derivative at locations where the surface's curvature is maximum. Solving the second derivative for its extreme values requires determining the angle that maximizes the second derivative. This task is complicated by the function's cross product term which prohibits the angle theta from being isolated. Notation 1 can be rewritten as:

(Notation 2) Notation 1

where x^tis the transpose of x. This quadratic form of the derivative provides a way for the individual components of the second derivative to be isolated. The symmetric 2x2 matrix H is a Hessian matrix which contains all the function's curvature information (Strang, 1976). The roots of these equations are the directions of the maximum and minimum surface curvature. Eigenvalue analysis is used to generate the first estimate of shape classifications (see Table 1) (Strang, 1976; Anton,1991).

Thus the second derivative of the neighborhood function can be used to generate a theoretically robust description of surface shapes. The curvature analysis was used to develop computer code capable of classifying terrain in a single pass. However, the usual resolutions of DEMs are unlikely to capture the smooth elevation contrasts that are the goal of "shape of country" macrostructure landform analysis. In addition, the resolution of the DEMs is unlikely to always result in the highest value classes, ridges and valleys, for example, being located in the center of the center cell in the 3x3 roving windows used in this kind of analysis.

Methods

The problem of grid features not being aligned with major features of the terrain is dealt with mathematically. If a cell that contains a ridge line that does not pass through the central point the procedure described above will classify the cell as a side hill even though it should be classed as a ridge since ridges contain more surface information than side hills. The solution to this problem was to use the continuous fit property of the neighborhood operator to determine if there is a change in slope within the cell. If surface slope remains constant within the cell then there is no change in surface structure within the cell. If there is a slope change within the cell then the type of inflection point can only be identified by evaluating the signs of the eigenvalues calculated at the exact location of the inflection point (Haralick et al,1983). Determining if there is a change in slope within the cell's area requires obtaining at least two additional measures of slope at points within the cell. The logical locations for testing sign change are between the central point and the cell's edges in the directions of the maximum and minimum curvature. If there is a change in sign between either or both of the pairwise comparisons then the location of the zero crossing(s) in the direction of the first derivative must be determined. Once the zero crossing(s) have been found then the Hessian matrix H is again used to compute two new eigenvalues that measure surface curvature. The cell is then classified according to the flow chart in Figure 2. The tables used in the processing are showing in Figure 3.

Figure 2. Flow chart used in determing the final shape class of 3x3 roving window

Figure 3. The three sets of hueristic tables refered to in Figure 2.

Findings

The process described was applied to several ideal mathematical surfaces, the same surfaces with random noise added, and 2 real DEMs. Figure 4 A is an IDRISI orthographic view of a mathematical saddle. This surface, if correctly classified, should contain four unique surface features: a single ridge line, a single valley line, a saddle point, and hillsides. As can be seen in Figure 4B the procedure correctly identified all of these features. No smoothing was applied since the surface was very smooth as generated. Figure 5A is a real, 30m DEM. Figure 5B shows the classification of this image. Although the DEM looks smooth in the rendering it is clear that it is not smooth from the classification shown. Microsurface features dominate although ridge and valley lines can be seen. However, the classification is not very satisfactory. Figure 6A shows the DEM after being filtered and Figure 6B shows the classification of the smoothed surface. This classification is much better. However, close observation of the original DEM and the classification shows that the ridge and valley lines are not always exactly where they should be. The cause of this is that the filtering has removed the short wavelength information from the image and that has shifted the locations of the some of the primary feature lines. Close study of the classification relative to the smoothed DEM shows it to be correct, however.

Figure 4 Classification of a mathematical saddle . A is an orthographic IDRISI view and B is its classification.

Figure 5. Classification of an unsmoothed DEM from Central New York. A is an orthographic IDRISI view and B is its classification.

Figure 6. The smoothed DEM of Figure 5 (A) and its classification (B).

Conclusion

It appears that this shape-of-country classifier does a very good job of classifying real DEMs into useful information. However, there is room for further research.

Recommendations for future research

There are several lines of research that can be explored. This classifier does not perform perfectly since the smoothing needed to remove roughness also causes shifts in the exact location of some of the primary features. The spatial frequency domain filter also tends to change elevations and create pits which can cause problems in hydrological modeling. Wechsler(2000) has shown that the range of spatial dependence as determined by Krieging varies over some DEMs. This means that the filtering needed to remove roughness should have a cutoff frequency that depends on the nature of the surface in the vicinity of the analysis. Thus it would appear that there is room for further research in the area of effective DEM smoothing. Until better filters are found there is not much reason for trying to improve on the classification procedure reported here.

References used

Anton H. 1991. Elementary Linear Algebra. New York: J Wiley.

Anton H. 1988. Calculus with Analytic Geometry. New York: J Wiley.

Blackman RB, Tukey JW. 1958. The Measurement of Power Spectra. New York: Dover.

Frank AU, Palmer B, Robinson VB. 1986. Formal Methods for the Accurate Definition of Some Fundamental Terms in Physical Geography. International Geographic Union, Commission of Geographical Data Sending and Processing. p 583-599.

Graff LH, Usery EL. 1993. Automated Classification of Generic Terrain Features from Digital Elevation Models. Photogram Eng Rem Sens 59 (9):1409-1417.

Haralick RM, Watson LT, Laffey TJ. 1983. The Topographic Primal Sketch. Int J Rob Res 2 (1): 50-72.

Larson RE, Hostetler RP. 1990. Calculus with Analytic Geometry. Lexington: DC Heath.

Lazaroff MB, DeFeo NJ, Stein MC, Craul PJ. 1990. Automatic Inference of Soil Characteristics from Topography, Surficial Geology and Remotely Sensed Imagery. ACSM-ASPRS Annual Convention (Technical Papers) IV. Baltimore: p 260-269.

Mark DM. 1975. Geomorphometric Parameters, A Review and Evaluation. In: Geografiska Annaler. 57A 3-4: 179-188.

Mitchel C. 1973. Terrain Evaluation: An Introductory Handbook to the History, Principles and Methods of Practical Terrain Assessment. London: Longman.

Pellegrini GJ. 1995. Terrain Shape Classification of Digital Elevation Models using Eigenvectors and Fourier Transforms.[dissertation]. Syracuse (NY): State University of New York, College of Environmental Science and Forestry. 207 p. Available from: University Microfilms, Ann Arbor MI; AAI9543009.

Spurr SH, Barnes BV. 1973. Forest Ecology. New York: Ronald Pr.

Strahler A. 1964. Quantitative Geomorphology of Drainage Basins and Channel Networks. In: Chow, editor. Handbook of Applied Hydrology. New York: McGraw-Hill.

Strang G. 1976. Linear Algebra and its Applications. New York: Academic Pr.

Wang X, Yin Z.1997. An Evaluation of Using ARC/INFO to Extract Basin Physiographic Parameters from DEMs. In: Proc Environmental Systems Research Institute User Conference. San Diego: ESRI.

Way DS. 1973. Terrain Analysis. Stroudsburg (PA): Douden, Hutchinson Ross.

Wechsler SP. 2000. Effect of DEM Uncertainity on Topographic Parameters, DEM Scale and Terrain Evaluation. [dissertation]. Syracuse (NY):State University of New York, College of Environmental Science and Forestry. 380p.

Zevenbergen LW, Thorne CR. 1987. Quantitative Analysis of Land Surface Topography. Earth Surface Processes Landforms 12:47-56.

Authors

Lee Herrington,Ph.D., Professor of Resources Information Management
State University of New York, College of Environmental Science and Forestry
1 Forestry Drive, Syracuse, New York 13210
Email: lpherrin@syr.edu Tel: +1-315-470-6674, Fax: +1-315-470-6535

Gerald Pellegrini, Ph.D., Engineering Specialist, Advanced Development
BAE Systems, 16250 Technology Dr., San Diego CA, USA, 92104
Email: Gerald.Pellegrini@marconi-is.com, Tel: +1-858-592-5556, Fax: +1858-592-5320