Digital Elevation Models in Geomorphology

Abstract

This chapter presents place of geomorphometry in contemporary geomorphology. The focus is on discussing digital elevation models (DEMs) that are the primary data source for the analysis. One has described the genesis and definition, main types, data sources and available free global DEMs. Then we focus on landform parameters, starting with primary morphometric parameters, then morphometric indices and at last examples of morphometric tools available in geographic information system (GIS) packages. The last section briefly discusses the landform classification systems which have arisen in recent years.

1. Introduction

Geomorphology, the study of the Earth’s physical land-surface features, such as landforms and landscapes and on-going creation and transformation of the Earth’s surface, is one of the most important research disciplines in Earth science. The term geomorphology was first used to describe the morphology of the Earth’s surface in the end of nineteenth century [1]. Geomorphological studies have focused on the description and classification of landforms (geometric shape, topologic attributes, and internal structure), on the dynamical processes characterizing their evolution and existence and on their relationship to and association with other forms and processes [2].

Geomorphology dates back to sixth to fifth century BC, when Xenophanes of Colophon (580–480 BC) speculated that, as seashells are found on the top of mountains, the surface of the Earth must have risen and fallen, or Herodotus (484–420 BC) thought that the lower part of Egypt was a former marine bay referring to the year-by-year accumulation of river-borne silt in the Nile delta region [3]. Geomorphology as an independent scientific discipline developed in the late nineteenth century [4]. The first modern theory of landscape evolution was the ‘geographical cycle’ expounded by Davis [5–7]. Another important theory, or rather a variation on Davis’ scheme was offered by Penck [8, 9] who stated that according to the Davis model, uplift and planation take place alternately.

Geomorphology in the twentieth century experienced rapid evolution and growth, and many overlapping phases of development can be identified. Contemporary geomorphology combines and consists of many individual fields of science, that is, geology, hydrology, meteorology, cartography, geographic information system (GIS), engineering, biology, archaeology, etc. This complexity and the long tradition associated with its development helped geomorphologist to use many research techniques and measurement tools according to the research requirements. It is often said about the four main research directions, which are historical, dynamic, structural, and climatic, in geomorphology. However, as rightly observed by Migoń [10], each area has its own history, dynamic nature of the processes, and is present in specific geological structures and with the participation of specific climatic conditions.

Herein, we will focus on research approach to landform as the main study subject in geomorphology. In general, over the last 50 years four primary approaches were distinguished in geomorphology.

1. Morphography: It describes physical appearance of landforms and is qualitative approach to the form. It is linked with the direct observation of forms in situ which allows specifying the appearance of the form and its morphographic classification (plain, hill, valley, ridge, etc.). These terms do not indicate the way of creation of forms rather only determine their external expressions.

2. Morphogenesis: It focuses on explaining the origin of the forms and determine the mechanisms of their contemporary development. Geomorphologists use different methods to determine the nature of the process in the past and the present form.

3. Morphochronology: It aims to specify age of the forms and the age-relationships between adjacent landforms. Geomorphologists examine both absolute as well as relative age between the forms.

4. Morphometry: It deals with establishing geometric features of the landforms on the basis of measurements. This chapter is just dedicated to morphometry that at the beginning was an element of quantitative geomorphology and later became an independent discipline in the Earth sciences—geomorphometry.

The quantitative phase in geomorphology was developed in 1940–1970 and reflected a broader trend within many of the Earth sciences disciplines towards enhanced use of sophisticated technologies (often derived from military purposes) to measure, describe and analyse the Earth’s surface features in number of categories. Horton’s publications [11, 12] on stream networks and drainage basin processes are classically identified as the precursor to this quantitative movement. Strahler [13] has limited a method for the quantitative analysis of the forms modelled by flowing water and gravitational movements over a longer period of time. The data were obtained mainly from measurements on detailed topographic maps and then applying statistical and morphometric analysis for various calculations and indicators of the fluvial relief within the drainage basin [14–18], karst [19], or glacial relief [20, 21] were made.

Further development of the geomorphometry resulted in working out of new theoretical and methodological basis. Lustig [22] distinguished two approaches for quantitative analysis and Earth’s surface characteristics which are (1) the characteristics of individual forms based largely on field surveys and (2) analysis of the surface area as a whole based on analysis of the map. Then Evans [23] proposed the division of geomorphometry into two forms: general and specific. In general, geomorphometry applies to and describes the continuous land surface. It provides a basis for the quantitative comparison of qualitatively different landscapes and it can adapt to the methods of surface analysis used outside geomorphology. Whereas specific geomorphometry applies to discrete landforms, describes selected relief or landform types as well as their geometry and laws of formation and development. Although such divisions may seem somewhat artificial, since the geometry of individual forms consists of the geometry of the surface as a whole, but they are useful to define more closely the principles of various types of morphometric analysis.

At the end of the 1980s and the beginning of the 1990s, with the personal computer revolution, algorithms have been implemented in many raster-based GIS packages (ArcInfo, MicroStation, MicroDEM, etc.) and it was possible to process digital elevation models (DEMs) over fairly large areas. Then began the new era of geomorphometry with new opportunities to visualize and compute land-surface parameters. Now-a-days, geomorphometry is an important component of terrain analysis and surface modelling including measurements of morphometry of continental ice surfaces, characterizing glacial troughs, mapping sea-floor terrain types, guiding missiles, assessing soil erosion, analysing wildfire propagation, and mapping ecoregions [24, 25]; it is a broad field that is important not only in various aspects of Earth sciences but also in engineering, biology, and medicine [26].

2. Digital elevation models (DEMs)

2.1. History and definition

As Pike [27] noted, a numerical description of the ground surface is helpful in addressing many geomorphological problems. In the definition of the subject of this chapter appears the word ‘model’, that is, a representation, generally in miniature, to show the construction or appearance of something. Meyer [28] neatly expressed it—reality scaled down and converted to a form which we can comprehend. In this case it is used to represent the original situation—approximation of topography. The term and concept of ‘terrain model’ was first described by Miller and Laflamme [29] but it did not come into general use until the 1960s and even later because of the technological limitations (computers, processors, etc). They defined digital terrain model (DTM) as just a statistical representation of continuous Earth’s surface which consists of many points of known co-ordinates x, y, z in the arbitrary co-ordinate system. In the following years, a number of terms related to the presentation of the Earth’s surface by numerical methods, that is, Burrough [30] described DEM as a regular gridded matrix representation of the continuous variation of relief over space. Among many terms, the most common and accepted in geomorphometric and GIS terminologies is digital elevation model (DEM) or digital terrain model (DTM). Some authors make distinctions between them, that is, DEM is ‘an ordered array of numbers that represent the spatial distribution of elevations above some arbitrary datum in a landscape’. DTM is ‘an ordered array of numbers that represent the spatial distribution of terrain attributes’; therefore DEM is a subset of DTM. Moore et al. [31, 32], Weibel and Heller [33], and Li et al. [34] defined DTM as a digital (numerical) representation of the terrain. Herein both terms will be treated as synonyms but, in fact, the concept of DEM is broader and more universal. In summary, we can say that digital elevation model is the set of digital data describing elevation values of Earth’s ground surface (or any other surface) which contains additional information about the character of this surface (i.e. structural lines, break lines, water bodies, etc.) and interpolation algorithm, which is the best for approximation (modelling) of the real topography. A DEM is a complete representation of a land surface which means that heights are available at each point in the area of interest [35].

2.2. Types

Due to geometry and data organization there are several basic types of models: raster (GRID), vector triangulated irregular network (TIN), point and hybrid [34]. Herein we will skip point and hybrid models because they are rarely used and absent in the most popular GIS software packages (i.e. ArcGIS, QGIS, MapInfo, Surfer).

Raster model (GRID) is the most widely used digital height data structure during the last years because of its simplicity [36]. The simplicity consists, in fact, that elevations are stored using a regular square grid that is consistent in each part of the study area. All square grids form regular matrix of heights for which plain co-ordinates (x, y) can be easily calculated due to the regular spacing of the grid points (Figure 1). This kind of DEM has many advantages, such as simple elevation matrices that record topological relations between data points implicitly and ease of computer implementation [31, 37]. Moreover, DEM is considerably easier to design land-surface parameters and objects using grids because simpler algorithms can be used; grids have a uniform spatial structure and almost all properties of gridded DEMs are defined by a single characteristic which is cell size; a grid model is more suited to the computer models used in image processing and for printing [35]. For this reason, some DEM software packages accept only grid data.

Figure 1.

Example of the GRID model (A), and TIN model (B).

Raster DEMs also have disadvantages, such as grids present under-sampled topography in areas where the topography is complex, and they over-sample smooth topography; re-projection of a grid is slow and sometimes leads to a loss of accuracy (because the initial grid loses its regular structure in a new projection and so it has to be re-calculated); the different distances between grid centres in cardinal and diagonal directions have a negative impact on the precision of many hydrological modelling [35]; the computed upslope flow paths will tend to have zigzag pattern across the landscape and increase the difficulty of calculating specific catchment areas accurately [38, 31]; square grids cannot handle abrupt changes in elevation easily and they will often skip important details of the land surface in flat areas [39].

The second most popular and widely used is a vector model named triangulated irregular network (TIN). TIN is proposed by Hormann [40] which devised a TIN idea, linking selected points on divides, drainage lines and breaks in slope to interrelate height, slope gradient and aspect. TIN is based on triangular elements (facets that vary in shape and size) with vertices at the sample height points [31]. These facets consist of planes joining the three adjacent points in the network and are usually constructed using Delaunay triangulation [33]. TIN can be interpolated directly from surveyed points or discrete features that are extracted manually from maps or by computer from a grid or contour DEM [27]. The triangle may be regarded as the most basic and universal unit in all geometrical patterns (because of their great flexibility in terms of shape and size), since a regular grid of square or rectangular cells or any polygon with any shape can be decomposed into a series of triangles [34].

The advantages of such a model are: relevance to gravitational movements, and especially to hydrological applications [41]; TINs are widely used in perspective representations of surfaces, especially in dynamic fly-through displays where the foreground is represented in full detail but the background can be simplified as larger triangles; for areas with high relief or rougher surfaces, irregular DEMs can use smaller spacing between points, and larger spacing where relief is lower or where the surface is smoother and in this way they can more accurately describe geological faults and other sharp elevation changes using the same number of points as grids [35]; the TIN structure gives best reflecting processes of erosion and deposition mimics paths of steepest gradient [41]. Of course we can use the algorithms to convert the TIN into grid and vice versa, but usually at the loss on the quality and accuracy of the model.

2.3. Data sources

DEM source means the data acquired by different techniques and from different sources. In general, such data can be acquired from topo maps and from field surveys.

Manual or semi-automated digitizing of contour lines and every single elevation point on topographic maps on the computer screen or digitizer can be done relatively cheaper without the need for resurveying. Derived DEMs represent the underlying terrain without surface vegetation and buildings. Despite many critical voices regarding the use of this method [42, 43], it is a quick and effective method to get quite good DEM and till the end of 1990s this was the most common source of elevation data for DEMs.

The second group of data sources are surveys, such as ground survey techniques, digital photogrammetry, and remote sensing. Field surveys are carried out using total stations, theodolite, levelling instruments, global positioning system (GPS), and differential global positioning system (DGPS). These kinds of surveys have high accuracy (sometimes even less than 1 cm), flexibility (the measurement density can be varied, depending on the terrain), and very little processing is required after the measurements have been taken, but especially suited for measuring small areas. Digital photogrammetry relies on the stereoscopic interpretation of aerial photographs or satellite imagery using manual or automatic stereoplotters [39, 33]. Remote sensing surveys consist of Airborne laser scanning (ALS) or satellite platforms and include laser-ranging altimetry (LiDAR—light detection and ranging), synthetic aperture radar interferometry (InSAR or IfSAR). The final results will often include tree-top canopies and buildings. This gives higher elevation values, rough surfaces, and high slope values [43]. The advantages of the LiDAR method are production time that is typically shorter than that for photogrammetrically generated DEMs [44] and great spatial and vertical accuracy (<0.5 m). The main disadvantage of LiDAR data is point cloud which produces a very dense and detailed land-surface model that could be difficult to handle during the production process and sometimes the accuracy of the readings vary according to the characteristics of the terrain (i.e. very steep slopes) [45]. However, these days LiDAR is definitely the best method of the DEM production. Several countries have already produced national LiDAR DEM/DSM (e.g. Belgium, the Netherlands at resolutions of 2–5 m, Poland at 1 m).

2.4. Free global DEMs

Currently, many elevation data with better accuracy and parameters are freely available in the world including LiDAR data. Individual countries or institutions offer different kinds of models. In Table 1, there are models (and their properties) available for the whole world, completely free of charge.

Model	Spatial resolution (grid size)	Accuracy vertical/horizontal	Co-ordinate System	Data format	Data extent	Source	Institution, publication
GTOPO30	30 × 30″ (~1000 m)	10–300/150 m	WGS-84	DTED, USGS DEM, DCW	90°N to 90°S	https://lta.cr.usgs.gov/GTOPO30	USGS (USA) and NASA, UNEP/GRID, USAID, INEGI (Mexico), GSI (Japan), MWLR (New Zealand), SCAR 1996
GLOBE DEM	30 × 30″ (~1000 m)	250/160 m	WGS-84	DTED	90°N to 90°S	https://www.ngdc.noaa.gov/mgg /topo/gltiles.html	NOAA, 1999
SRTM v3	1 × 1″ (~30 m); 3″ × 3″ (~90 m)	10/13 m	WGS-84/EGM96	DTED, BIL, GeoTIFF, hgt	60°N to 54°S	https://earthexplorer.usgs.gov/	NASA, NGA, 2013
ETOPO1	1 × 1″ (~30 m)	?	WGS-84	netCDF, GRD98, xyz, GeoTIFF	90°N to 90°S	https://www.ngdc.noaa.gov/mgg/global/global.html	NOAA, 2009
ASTGTM v02; AST14DEM v03	1 × 1″ (~30 m); 1 × 1″ (~30 m)	20/30 m	WGS-84/EGM96; WGS-84/UTM	GeoTIFF	83°N to 83°S	https://search.earthdata.nasa.gov/search?q=ASTGTM+V002	NASA, METI, 2009, 2011
GMTED2010	7.5′ × 7.5″ (~250 m); 15 × 15′ (~500 m); 30″ × 30′ (~1000 m)	26–30 m/?; 29–32 m/?; 25–42 m/?	WGS-84	ESRI ArcGrid,; GeoTIFF	84°N to 56°S	https://topotools.cr.usgs.gov/gmted_viewer/viewer.htm	USGS, NGA (USA), 2011
AW3D30	1 × 1″ (~30 m)	5/5 m	ITRF97/GRS80	GeoTIFF	82°N to 82°S	http://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm	JAXA, 2016

Table 1.

Basic information about the free global DEMs.

GTOPO30 (Global 30 arc-second Elevation): This is a global DEM with a horizontal grid spacing of 30 arc seconds (approx. 1 km) resulting in a product having dimension of 21,600 rows and 43,200 columns. GTOPO30 was derived from eight several raster and vector sources of topographic information which are digital terrain elevation data (DTED), Digital Chart of the World, USGS 1-degree DEM’s, army map service (AMS) 1:1,000,000, International Map of the World 1:1,000,000, Peru Map 1:1,000,000, New Zealand DEM, and Antarctic Digital Database [46]. The horizontal accuracy is ±150 m linear error at 90% confidence and vertical accuracy is 10–300 m at the 90% confidence level [47–49].

GLOBE DEM (Global Land 1-km base elevation DEM): It is one of the first global DEM, which was initially spearheaded in 1990. This data set covers 180^°W to 180^°E longitude and 90^°N to 90^°S latitude. The horizontal resolution is 0.5 arc-minute in latitude and longitude, resulting in dimensions of 21,600 rows and 43,200 columns [50]. GLOBE version 1.0 has 11 broad sources of information: 6 gridded DEMs and 5 cartographic sources, which were adapted for use in GLOBE. These were digital terrain elevation data (DTED), Digital Chart of the World, Australian DEM, Antarctic Digital Database, Brazil Maps 1:1,000,000, DEM for Greenland, army map service (AMS) Maps 1:1,000,000, DEM for Japan, DEM for Italy, DEM for New Zealand, and Peru Map 1:1,000,000 [51]. Vertical accuracy expressed as ±30 m linear error at 90% confidence can also be described as a root mean square error (RMSE) of 18 m. Linear error distribution is 97 m for RMSE and it can be expressed as ±160 m linear error at 90% confidence.

Shuttle radar topography mission v3 (SRTM): It was flown to aboard the space shuttle Endeavour February 11–22, 2000. The National Aeronautics and Space Administration (NASA) and the National geospatial-intelligence agency (NGA) participated in an international project to acquire radar data that was used to create near-global set of land elevations. SRTM successfully collected radar data over 80% of the Earth’s land surface between 60°N and 56°S latitude with data points posted every 1 arc-second (approximately 30 m). Absolute height error of SRTM data sets is 5–10 m and absolute geolocation error is 7–13 m [52]. The level of processing and the resolution of the data vary by SRTM data set:

1. SRTM Non-Void Filled. This version was edited or finalized by the NGA to delineate and flatten water bodies, better define coastlines, remove spikes and wells, and fill small voids. Data were sampled at 1 arc-seconds (USA) and 3 arc-seconds (rest of world);

2. SRTM Void Filled. This elevation data are the result of additional processing to address areas of missing data or voids in the SRTM Non-Void Filled collection. The resolution for SRTM is 1 arc-seconds (USA) and 3 arc-seconds (rest of world);

3. SRTM 1 Arc-Second Global. This elevation data offer worldwide coverage of void filled data at 1 arc-second (30 m) resolution and provides open distribution of this high-resolution global data set. Some tiles may still contain voids. Please note that tiles above 50°N and below 50°S latitude are sampled at a resolution of 2 by 1 arc-second [53].

ETOPO1 (Global Relief Model at 1 arc-minute resolution): It was made in 2008 by the National Geophysical Data Centre (NGDC), an office of the National Oceanic and Atmospheric Administration (NOAA). This model was developed as an improvement to the ETOPO2v2 Global Relief Model. ETOPO1 is available in two versions: ‘Ice Surface’ (top of Antarctic and Greenland ice sheets) and ‘Bedrock’ (base of the ice sheets). These versions of ETOPO1 were generated from diverse (regional and global) digital data sets, which were shifted to common horizontal and vertical datum. Next steps were evaluated and edited as needed [54]. Shoreline, bathymetric, topographic, integrated bathymetric–topographic, and bedrock digital data sets (13 sources) were obtained from several U.S. government agencies, international agencies, and academic institutions.

ASTGTM (ASTER global digital elevation model v002 and AST14DEM—ASTER digital elevation model v003): These DEMs were developed jointly by the US NASA and Japan’s Ministry of Economy, Trade and Industry (METI). ASTER is capable of collecting in-track stereo using nadir and aft-looking near infrared cameras. Since 2001, these stereo pairs have been used to produce single-scene (60 × 60 km) DEMs having vertical RMSE accuracies generally between 10 and 25 m [55]. The ASTER GDEM covers land surfaces between 83°N and 83°S and is comprised of 22,702 tiles. This model is distributed as georeferenced tagged image file format (GeoTIFF) files. The data have resolution at 1 arc-second (approximately 30 m) grid and referenced to the 1984 World Geodetic System (WGS84)/1996 Earth Gravitational Model (EGM96) geoid. Although the ASTER GDEM v. 002 is better model than ASTER GDEM v. 001, users have to know that the data still may contains anomalies and artefacts. One should know that these mistakes can introduce large elevation errors on local scales [56].

GMTED2010 (Global Multi-resolution Terrain Elevation Data 2010): This model was made by collaborating USGS and the NGA, which replaced GTOPO30 model. The new elevation data set has been generated at three separate horizontal resolutions of 30 (about 1 km), 15 (about 500 m), and 7.5 arc-seconds (about 250 m). The global aggregated vertical accuracy of GMTED2010 can be summarized in root mean square error (RMSE). At 30 arc-seconds, the RMSE range is 25–42 m; at 15 arc-seconds, range is 29–32 m, and at 7.5 arc-seconds, range is in between 26 and 30 m. This new product suite provides global coverage of all land areas from latitude 84°N to 56°S for most products, and coverage from 84°N to 90°S for several products [57]. An additional advantage of the new multi-resolution model over GTOPO30 is that at each resolution seven new raster elevation data sets are available. The new models have been produced using the various aggregation methods: minimum elevation, maximum elevation, mean elevation, median elevation, standard deviation of elevation, systematic subsample, and breakline emphasis. GMTED2010 is based on data derived from different raster-based elevation sources, such as SRTM DTED, non-SRTM DTED, Canadian digital elevation data (CDED), Satellite Pour l’Observation de la Terre (SPOT 5), Reference3D, National elevation dataset (NED), GEODATA, an Antarctica satellite radar and laser altimeter DEM, and Greenland satellite radar altimeter DEM [57].

AW3D30 (ALOS Global Digital Surface Model): The global digital surface model (DSM) at 1 arc-second (approximately 30 m) resolution that was released by Japan aerospace exploration agency (JAXA). This model has been compiled with images acquired by the advanced land observing satellite (ALOS). The elevation data are published based on the DSM data set (5-m mesh version) of the ‘World 3D Topographic Data’ [58]. The source data were a huge amount of stereo-pairs images derived from satellite mission in the years 2006–2011. Next, they were processed semi-automatically to provide digital surface model (DSM). The height accuracy of the data set is approximately <5 m from the evaluation with ground control points (GCPs) or reference DSMs derived from the LiDAR [59].

In addition to the above global DEMs, there are many elevation data with much better accuracy, but for smaller areas, for example, for Europe EU-DEM with grid size 25 m, for Spain MDT05/MDT05-LIDAR (5 m), for Netherlands AHN3 DTM (0.5 m), for all Slovenia (LiDAR data), or LiDAR data for Haiti (1 m).

3. Morphometric landform properties

The main object of the studies in geomorphology is the relief (of the Earth’s surface). The term relief of the Earth’s surface was used to describe the vertical dimension or amplitude of topography [60] or as the elevation difference over a pre-determined area or collective elevations and their inequalities of a land surface [61]. One can also say that the relief is the complexity of the Earth’s surface shapes, and these shapes form the landforms. If we would present a continuous Earth’s surface as a matrix of discrete points with a defined distance interval from one another, the shape dimensions would express in the changes of height, distance, and direction between these points. These three vectors and the relationships between them are the basis for most morphometric landform properties.

3.1. Morphometric parameters

Elementary data used in geomorphometry are DEM. Strength of DEM is a suggestive (plastic and three-dimensional) visual message and the ability of quantifying the topography. Quantification of the Earth’s surface is expressed by topographic attributes which can be determined from the derivatives of the topographic surface. Generally, one can say, these derivatives measure rate at which elevation changes in response to changes in location. Deng et al. [62] noted that local terrain shape (as the continuous variation of elevation values over the terrain surface from point to point) has an enormous impact on local terrain attributes, but this role is influenced by data and computational factors.

Theoretical assumptions concerning the morphometric properties of the surface come from the 1950s [63–65], and even earlier. Stone and Dugundji [66] and Hobson [67] stated that measures of landforms can be considered as a kind of the roughness of the surface. In general, roughness refers to the irregularity of a topographic surface and cannot be completely defined by any single measure but must be represented by a roughness vector or set of parameters. With roughness, concept of wavelength and amplitude ideas was related. The significant wavelengths of topography are termed grain or texture, while amplitudes associated with these wavelengths correspond to the concept of relief [60]. Texture and grain are terms that have been used to indicate in some way the scale of horizontal variations in the topography. Texture is used to refer the shortest significant wavelength in the topography and grain is used for the longest significant wavelength. Texture is related to the smallest landform elements, while the grain is related to the size of area over which one measures other parameters. Wood and Snell [68] defined grain as the size of area over which the other factors are to be measured. It is dependent on the spacing of major ridges and valleys and thus indicates texture of topography. Texture refers to the shortest significant topographic wavelength.

The systems proposed by Evans [23] and Krcho [69] for field variables are local-based, relating to a vanishingly small area around each point. They divided topographic properties into primary geomorphometric parameters (height, slope, aspect, profile, and planar curvature) and statistical measures derived in square-matrices out of the gridded DEMs (such as relief, standard deviation, and skewness of heights).

Mark [60] stated, that probably the most important single class of processes which has shaped the Earth’s surface could be divided into geometrical properties that involve the relationships among dimensional properties, such as elevation, lengths, areas, volumes, and topological properties which relate numbers of objects in the drainage network (e.g. the bifurcation ratio). Then Moore and Thornes [70] developed LEAP-land erosion analysis programs to examine the spatial distribution of slope length, slope steepness, and the plan curvature for assessing topographic erosion potential.

Pike [71] noted that terrain can be abstracted using geometric signature—a set of measurements that describes topographic form sufficiently well to distinguish geomorphically disparate landscapes. Then he developed the concept of a geometric signature, a multi-variate description of topography using a suite of measures, and later expanded the concept with a listing of 49 variables that could be grouped into 22 attributes [72, 73]. He considered roughness and height, the two most important attributes, with two measures of texture at seventh and eleventh position. Fifteen different variables contribute to roughness. Pike et al. [74] also referred to grain concept. Relief at topographic grain is an estimate of local relief optimized by varying unit-cell size. In homogeneous terrain, local relief within nested circles increases with circle size and then levels off at a diameter termed grain, a measure of characteristic local ridgeline-to-channel spacing. Topographic grain is the characteristic horizontal spacing of major ridges and valleys and is an important descriptor of meso-scale topographic texture. The grain concept arose from the need for a variable and non-arbitrary unit-cell size and also enables to calculate local relief parameter.

Moore and Grayson [41] and Moore et al. [31, 32] described terrain attributes adapted from Speight [75, 76], and then often mentioned and cited [3, 36, 77–80]. They divided terrain attributes into primary and secondary (compound) attributes. Primary attributes, hydrologically related, are directly calculated from elevation data and include variables, such as altitude, upslope height, aspect, slope, upslope slope, dispersal slope, catchment slope, upslope area, dispersal area, catchment area, specific catchment area, flow path length, upslope length, dispersal length, catchment length, profile curvature, and plan curvature. It is interesting that Bork and Rohdenburg [81] have developed a digital relief model for estimating and displaying the distribution of these morphographic parameters.

Schmidt and Dikau [78] subdivided primary geomorphometric parameters into three types: simple, complex, and combined. Simple primary geomorphometric parameters (usually derived through a filter operation within a 3 × 3 moving window) are height, slope, aspect, profile curvature, contour curvature, drainage direction, and real area of pixel. Complex geomorphometric parameters are derived through the analysis of the whole matrix of a DEM. They contain structural information about the surrounding morphometry and consist of contributing area, mean slope of contributing area, average, and variance of primary parameters in contributing area, length of flowpath to outlet, average and variance of primary parameters in flowpath to outlet, length of flowpath to stream, average and variance of primary parameters in flowpath to stream, x- and y-co-ordinates of corresponding stream point, height of corresponding stream point, height distance to corresponding stream point, length of minimum flowpath to watershed, relative slope position after minimum, slope length, relative slope position after maximum, and minimum and maximum slope length.

Interesting approach to the morphometric parameters was presented by Shary et al. [82], who stated that morphometric variables describe not the land surface itself, but rather the system, land surface + vector field, where vector fields of common interest are gravitational field and solar irradiation. He divided morphometric variables and concepts into field-specific (may refer to this system description) and field-invariant (invariant with respect to any vector field, that is, describing the land surface geometrical form). On the other hand, morphometric variables may divide into local, regional, or global (when height data of all the Earth is needed for their determination).

Basso [79] divided topographic attributes into: (i) local (calculated from a small neighbouring area surrounding the DEM cell, usually 3 × 3), (ii) regional (calculated using considerably larger geometric area than the local attributes, less sensitive to the DEM resolution), (iii) catchment oriented (related to the whole catchment area, and are the measurement of certain catchment characteristics), and (iv) process oriented (describe or characterize the spatial variability of a simple representation of specific processes that occur on the landscape).

Goodwin and Tarboton [84] described five categories of the morphometric parameters of drainage basin which are size properties (provide measures of scale that can be used to compare the magnitudes of two or more drainage basins), surface properties (quantities depicted by fields comprising a value at each point within a domain, drainage basin), shape properties (i.e. length, width, perimeter, and more complex function), relief properties (total basin relief, relief ratio [85], and hypsometric curve), and texture properties (amount of landscape dissection by a channel network).

Olaya [86] presented division of the morphometric parameters into local and regional. Local parameters consist of geometric (slope, aspect, curvatures, visibility, visual exposure, and visibility index) and statistical (i.e. average value, standard deviation, skewness coefficient, kurtosis coefficient, range of values, etc.). Regional parameters are connected with hydrological properties (catchment area, height, slope, proximity, etc. [86]).

The latest approach to the classification system of the fundamental geomorphometric variables is presented by Evans and Minár [87]. They proposed field and object variables. Field variables include specific to gravity field (local point-based, local area-based, and regional), specific to other fields and field-invariant variables (local point- and regional-based). Object variables differ between areal, linear, and point features.

Thus the optimal number of variables depends on spatial scale, resolution of the source data (DEM), and the requirements of the research problem; and many measures describe the same attribute of surface form and thus are redundant; for this reason new geomorphometric parameters are very rare.

3.2. Morphometric indices

Morphometric parameters, which were discussed in the previous section, showed that there are many different classifications of these parameters. In this section, we want to look at some popular morphometric indices called secondary or composite attributes [3, 32, 36, 62, 80, 89], combined or compound geomorphometric parameters [78, 83], statistical parameters [23, 86] or process oriented [79]. These indices are combinations of the primary morphometric attributes and describe or characterize the spatial variability of specific processes occurring on the landscape. Sometimes these morphometric indices can be derived empirically but it is preferable to develop them through the application and simplification of the underlying physics of the processes. With the index approach we simplify some physical sophistication to allow improved estimates of spatial patterns in the landscape [31].

In Table 2, some selected geomorphometric indices commonly used and their definitions were listed.

Parameter	Formula	Description	Source
Drainage density	DD= L/A; where: L is sum of the channel lengths and A is basin area	The sum of the channel lengths divided by basin area. It is important indicator of the linear scale of landform elements in a drainage basin, indicates the closeness of spacing of channels, thus providing a quantitative measure of the average length of stream channel for the whole basin.	Horton [12]
Form factor	Rf= A(lb)^2; where: A is area of basin, Lb is the basin length	The ratio of the basin area to the square of the basin length. Indicates the flow intensity of a basin of a defined area. The form factor value should be always less than 0.7854 (the value corresponding to a perfectly circular basin).	Horton [12]
Terrain ruggedness index (TRI)	TRI = Y[Σ(xij−x00)^2]^1/2; where: xij is elevation of neighbour cell to cell (0,0)	Measure of topographic heterogeneity (amount of elevation difference between adjacent cells of a DEM).	Riley et al. [90]
Stream power index	Ω = pgq.tanβ); where: pg is the unit weight of water, q is the discharge per unit width and β is the representative slope angle	This is the time rate of energy expenditure and has been used extensively in studies of erosion, sediment transport and geomorphology as a measure of the erosive power of flowing water.	Moore et al. [31]
Elevation-relief ratio	E = (Hmean − Hmin)/(Hmax − Hmin); where: H is elevation	Expresses relative proportion of upland to lowland within a sample region. Usually ranges from 0.15 to 0.85. Low values occur in terrains characterized by isolated relief features standing above extensive level surfaces, and high values describe broad and level surfaces, broken by occasional depressions.	Pike and Wilson [91] after Wood and Snell [68]
Topographic openness algorithm	Positive openness; φL = (0φL+45φL+…+315φL)/8; negative openness; ⵖL =(0ⵖL+45ⵖL+…+315ⵖL)/8; where: L is specified distance	Describes the degree of enclosure of a location on an irregular Earth’s surface.; Positive values expressing openness over surface topography is high for convex forms, whereas negative values expressing and describing attribute below the surface topography and are high concave forms.	Yokoyama et al. [92]
Terrain shape index	_; TSI = Z/R; where: Z is mean elevation of the sample plot boundary, and R is plot radius measured in the units used for elevation	This index is equivalent to the mean slope gradient of the plot boundary as viewed from the plot centre, with units of meters change in elevation per meter of plot radius.; Typical TSI values for mountain landforms ranged from −0.24 to −0.12, on convex upper slopes near ridge tops, and from +0.09 to +0.17 on concave lower slopes.	McNab [93]
Landform index	LI = H^°/100; where: H^° is vertical gradient to the horizon	The landform index is the average vertical gradient to the topographic horizon, divided by 100 to convert percent to a decimal value.; The index is dimensionless and the effects of height and distance to the landform are compensating factors.	McNab [94]
Compound topographic index (CTI)	CTI = ln (Af/tanβ); where: Af is the specific catchment area draining through the point, and β is the representative slope angle	Ratio between slope and catchment area; quantification of catenary topographic convergence represented by slope angle and catchment. For the same contributing area CTI values are higher for pixels with lower slopes —this means that CTI primarily reflects accumulation processes.	Moore et al. [32]
Heat load index	HLI = [1 − cos(θ−45)]/2; where: θ is aspect in degrees east of north	Quantitative measure of aspect and steepness of slope. The equations applied to 0–60° north latitude, slopes from 0 to 90° and all aspects.	McCune and Keon [95]
Anisotropy index	ANI = Ȓmin/Ȓmax; where: Ȓmin is smallest estimated range parameter, and Ȓmax is highest estimated range parameter in various directions	Ratio between the minimum and maximum range parameter of spatial dependence, fitted for various directions	Bishop and Minasny [96]
Shape complexity index	​SCI = P / (2rπ ) , r = ​√ __ ​ A __ π ​ ​​; where: P is the perimeter of the polygon, A its area, r is the radius of the circle with the same area	Index which is used to describe polygons on DEM slices. Indicates how compact (or oval) a feature is.	Hengl et al. [97]
Basin relief ratio	Rh = H/L where: H is total basin relief, and L is basin length	Ratio between total basin relief (difference in elevation of basin mouth and summit) and basin length, measured as the longest dimension of the drainage basin. Indicates overall slope of the watershed surface. It is a dimensionless number, readily correlated with other measures that do not depend on total drainage basin dimensions.	Schumm [98]
Relative relief	Rp = H/P where H is total basin relief, and P is basin perimeter	Ratio between total basin relief and drainage basin perimeter.	Melton [99, 100]
Drainage basin compactness	Bc = P/A; where: P is drainage basin perimeter, and A is drainage basin area	Ratio between perimeter and area of drainage basin. Higher values correspond to the basins of developing the long-term share erosion running in conditions of relative peace tectonic or are typical for catchment formed in low resistance rocks.	Engstrom [101]
Drainage basin shape ratio	Bs = Bl/Bw; where: Bl is max length of the drainage basin, and Bw is max width of the drainage basin	Ratio between maximum length and maximum width of drainage basin. Higher values correspond to more elongated basins and also indicate a relatively higher tectonic activity of the area.	Cannon [102], Ramírez-Herrera [103]
Surface area ratio	SAR = A/Ap; where: A is cell’s surface area, and Ap is planimetric area of that cell	Ratio between surface area and planimetric area; surface ratios will always be greater than or equal to 1.	Jennes [104]
Valley height-width ratio	Vf = 2Vfw/[(Eld − Esc) + (Erd − Esc)]; where: Vfw is width of the valley floor, Eld/Erd height of the left/right watershed and Esc height of the valley floor	Ratio of the width of the valley bottom and the average height of the slopes; allows comparison of erosional patterns between watersheds. Low index value Vf characteristics are deeply cut river valleys developing under high uplift area, while high values correspond to the valleys in the wide areas of thin tectonic activity.	Bull and McFadden [105]
Dissection index	DI = RR/AR; where: RR is relative relief, and AR is absolute relief	Ratio between relative and absolute relief which always varies between zero (complete absence of dissection) and one (extreme case, vertical cliff). Assess the degree of incision of a landscape.	Sharma [106]
Relief index (RI)	RI = CL/AP; where: CL is total length of contour lines, and AP is planar surface area	Ratio of the summary length of the contour lines and the planar surface area at which they occur; shows relief variability. RI is based on a combination of local relief (number of contour lines and elevational changes) and degree of surface cut (length and shape of the contour lines) with reference to the planar surface area.	Szypuła [107]

Table 2.

Selected geomorphometric indices derived from DEMs.

3.3. Morphometric tools available in GIS packages

The current availability of high speed computing platforms and high-resolution (less than 10 m spatial resolution) DEMs now provides the opportunity to perform quantitative analyses and calculating morphometric indices on new level. Many GIS packages offer tools to work with DEMs. On one hand we have the comWmercial software (ArcGIS, MapInfo, Surfer, Global Mapper, Terra Solid), and on other one, great free programs (SAGA, QGIS, MicroDEM), which offer a lot of interesting tools. In Table 3, we have presented useful tools for geomorphometric (or geomorphological) analysis. As examples we chose one commercial package (ArcGIS) and one free (SAGA). One should remember that rapid change of the GIS applications cause, next new versions a large number of greater tools.

SAGA 3.0.0	ArcGIS 10.5
Terrain analysis—Channels; Channel network; Channel network and drainage basins; Overland flow distance to channel network; Strahler order; Valley depth; Vertical distance to channel network; Watershed basins; Watershed basins (Extended)	3D analyst tools—Raster Surface; Aspect; Contour; Contour list; Contour with barriers; Curvature; Cut fill; Hillshade; Slope
Terrain analysis—compound analyses; Elevation; Analytical hillshading; Slope; Aspect; Plan curvature; Profile curvature; Convergence index; Closed depressions; Total catchment area; Topographic wetness index; LS-factor; Channel network; Drainage basins; Channel network base level; Channel network distance; Valley depth; Relative slope position	3D analyst tools—visibility; An overview of the visibility toolset; Construct sight lines; Intervisibility; Line of sight; Observer points; Skyline; Skyline barrier; Skyline graph; Sun shadow volume; Viewshed; Viewshed 2; Visibility; Spatial analyst tools—solar radiation; Area solar radiation; Points solar radiation; Solar radiation graphics; * Exemplary additional installed tools :
Terrain analysis—hydrology; Cell balance; Downslope area (interactive); Edge contamination; Flow accumulation (Flow Tracing); Flow accumulation (Mass-Flux Method); Flow accumulation (Recursive); Flow accumulation (Top-Down); Flow depth (interactive); Flow path length; Flow sinuosity (interactive); Flow width and specific catchment area; Isochrones constant speed (interactive); Isochrones variable speed (interactive); LS factor; LS-factor, field-based; Lake flood; Lake flood (interactive); Melton ruggedness number; SAGA wetness index; Slope length; Slope limited flow accumulation; Stream power index; TCI low; Topographic wetness index (TWI); Upslope area; Upslope area (interactive)	Hypsometry tools; Hypsometric curve; Area gradient tools; Max watershed; Surface area analysis; Surface area; Surface ratio; Flat area; Geomorphometric and gradient metrics; Classify aspect; Linear aspect; Mean slope; Slope impedance; Dissection; Hierarchical slope position; Landform curvature; Roughness; Slope position; Surface/area ratio; Surface relief ratio; Compound topographic index; Heat load index; Integrated moisture index; Site exposure index; Topographic radiation aspect index
Terrain analysis—lighting, visibility; Analytical hillshading; Potential incoming solar radiation; Sky view factor; Topographic correction; Topographic openness; Visibility (points); Visibility (single point) (interactive)	Topography tools; Beer’s aspect; McCune and Keon heat load index; Landform classification; PRISM data helper; Slope position classification; Solar illumination index; Topographic convergence/wetness index; Topographic position index
Terrain analysis—morphometry; Convergence index; Convergence index (Search radius); Curvature classification; Diurnal anisotropic heating; Downslope distance gradient; Effective air flow heights; Fuzzy landform element classification; Hypsometry; Land-surface temperature; Mass balance index; Morphometric features; Morphometric protection index; Multi-scale topographic position index (TPI); Multi-resolution index of valley bottom flatness (MRVBF); Real surface area; Relative heights and slope positions; Slope, aspect, curvature; Surface specific points; TPI-based landform classification; Terrain ruggedness index (TRI); Terrain surface classification (Iwahashi and Pike); Terrain surface convexity; Terrain surface texture; Topographic position index (TPI); Upslope and downslope curvature; Valley and ridge detection (Top hat approach); Vector ruggedness measure (VRM); Wind effect (Windward/Leeward Index); Wind exposition index	Local relief model; Terrain tools; Terrain ruggedness (VRM); Riparian topography tools; Calculate flooding height; Calculate inundation area; Height above river; Prepare HAR for flooding; Flow accumulation for both positive and negative values; Land facet analysis; Calculate density raster; Shannon’s diversity index; Identify termini polygons; Topographic position index tools; Mahalanobis distance tools

Table 3.

The list of exemplary morphometric tools available in SAGA and ArcGIS software (after [108–110]).

4. Landform classifications

Geomorphology studies the relief. If we want to understand the relief of the Earth’s surface (which is highly complex) we need to simplify and subdivide it into landforms. We have to focus on describing landforms, their spatial arrangement and the processes which led to their formation. Of course, many landforms can be delineated manually using photo-interpretation to assess their form, size, scale, adjacency, surface roughness, hydrological and contextual position but there is always problem with the boundary of the landform. And herein DEMs are helpful. Availability of global DEMs and high accuracy national LiDAR data made new possibilities of analysing these data to extract and classify geomorphic entities. The landform elements can be extracted automatically from DEMs by using land-surface parameters, such as slope, curvatures, catchment area, distance to streams, peaks and depression depth, etc. [111]. The goal of automated extraction of landforms and landform elements using semi-automated or fully-automated algorithms is to find their geometric signature, which Pike [72] defined as a set of measures that describe topographic form well enough to distinguish among geomorphically disparate landscapes.

There are many landform classification systems, and herein we take this issue very briefly and in general. The pioneer works in quantitative systems for landform classification were conducted by Hammond [112], Wood and Snell [68] and Anstey [113, 114] by using topographic contour maps. These studies were aimed primarily towards the systematization and logical interpretation of terrain data to assist in the determination of design criteria for materiel, and secondarily towards the development of a universal system for the quantification of landform data.

In the mid-1960s, Hammond [115] devised the three-level system of regional landform classification, based purely on geomorphometric parameters calculated in the chosen window size (Figure 2). For each window position the following parameters were calculated: (1) percentage of area where the ground is flat or gentle (less than 8% slope); (2) local relief (maximum minus minimum elevation); and (3) profile type (relative proportion of flat or gently sloping terrain that occurs in lowlands or uplands).

Figure 2.

Elevation map (A), and simplified Hammond’s landform classification (B).

Peucker and Douglas [116] showed a few methods designed to detect pits, peaks, passes, ridges, ravines, and breaks, given an array of sampled, quantized terrain heights. Described methods used local analysis which means the results at each point do not depend on the results already achieved at other points.

As the hillslopes constitute a basic element of all landscapes and a fundamental component of geomorphologic systems, many subjective classifications of slope profiles that were intended to create conceptual classifications of hillslopes have been proposed. Ruhe [117] divided hillslopes into summits, shoulders, backslopes, footslopes, and toeslopes. Dalrymple et al. [118] and Conacher and Dalrymple [119] proposed nine unit classification of hillslopes: (i) interfluve, (ii) seepage slope, (iii) convex creep slope, (iv) fall face, (v) transportational midslope, (vi) colluvial footslope, (vii) alluvial toeslope, (viii) channel wall, and (ix) channel bed (after [120]).

Next, there were several approaches to standardize hillslope units using qualitative terms [121]. Most commonly, a hillslope was described by a series of basic units describing changes in slope, curvature, and processes along the hillslope profile [122].

As Young [123] noted the curvature can be classified into convex, concave, and rectilinear surfaces, Dikau [124] defined basic form elements of the landscape as the combination of three slope profile curvature characteristics and three plan curvature characteristics which lead to nine possible hillslope units (Figure 3). They deliver a disjunctive description of the hillslope surface into units of homogeneous curvature characteristics. Then Dikau [124] proposed digital geomorphological relief model (DGRM) which will generate form facets and elements (as basic relief units) for geomorphological mapping and simulations for the derivation of more complex relief units. Dikau et al. [125] developed and applied an automated method for classifying macro landform types from DEM that was based on analysis of variation in topographic measures within areas defined by moving windows.

Figure 3.

Classification of form elements by plan and profile curvature (after Dikau [124]).

The original classification of Pennock et al. [126] explicitly assumed that surface form, as described by curvature, could be directly related to surface processes and to relative landform position (divergent/convergent shoulder, backslope footslope and level). Thus, strong profile convexity was assumed to be indicative of upper, water-shedding slope positions; whereas strong profile concavity was associated with lower, water-receiving landform positions, and planar surfaces were associated with backslopes or flat areas. Of course, this pattern is not always adhered to and there are many instances where convex-concave patterns repeat over short distances along a longer hill slope [121].

Speight [121] described 10 morphological types of topographic landform positions: (i) crest, (ii) depression (open, closed), (iii) flat, (iv) slope, (v) simple slope, (vi) upper slope, (vii) mid-slope, (viii) lower slope, (ix) hillock, and (x) ridge (after [120]).

Very simple and interesting method for classifying relief on the base DEM is topographic position index (TPI). TPI is the difference between the elevation at a cell and the average elevation in a neighbourhood surrounding that cell (Figure 4). Positive values means that the cell is higher than its neighbours (indicate ridges, hills, etc.) while negative values means the cell is lower (indicate canyons, valleys, etc). TPI is a simplification of the landscape position index (LPI) described by Fels and Zobel [127] and was developed in detail by Weiss [128]. TPI values provide a simple and powerful means to classify the landscape into morphological classes [129]. TPI is naturally very scale-dependent, it means neighbourhood size (and shape) and DEM resolutions are critical to the final analysis, so the work with this index should be based on experiments with different threshold values.

Figure 4.

Elevation map (A), and TPI 3-category slope classification (B).

Drăguţ and Blaschke [130] presented an automated classification system of landform elements based on object-oriented image analysis. Firstly, one need to derive elevation, profile curvature, plan curvature, and slope gradient from DEM. Next, relatively homogenous objects are determined at several levels through segmentation of the image. These object primitives are classified as landform elements using a relative classification model, built both on the surface shape and on the altitudinal position of objects. The classification has nine classes which are peaks and toe slopes, steep and flat/gentle slopes, shoulders and negative contacts, head slopes, and side and nose slopes. Classes are defined using flexible fuzzy membership functions.

Iwahashi and Pike [88] developed an iterative procedure that classifies topography automatically into terrain types, grid cell by grid cell, on the basis of three morphometric variables, such as slope gradient, local convexity, and surface texture. They applied an unsupervised nested-means algorithm for 16 topographic types of the world. They noted the procedure is unsupervised and reflects frequency distributions of the input variables but not defined criteria. It causes that resulting classes are undefined and have to be calibrated empirically by subsequent analysis.

Within the context of defining landform units that maximize internal homogeneity and external differences, Minár and Evans [131] presented the concept of elementary forms (segments and units) defined by constant values of fundamental morphometric properties and limited by discontinuities of the properties. The basic system of form-defining properties represents altitude and its derivatives, constant values of which provide elementary forms with various types of homogeneity.

Drăguţ et al. [132] presented an algorithm to derive elementary forms from DEMs. Elementary forms were defined by constant values of fundamental morphometric properties and limited by discontinuities of these properties. A multi-resolution segmentation technique was customized to partition the layers of altitude derivatives into homogeneous divisions through a self-scalable procedure, which reveals the pattern encoded within the data. Layers were segmented successively, following the order of elevation derivatives (i.e. gradient, aspect, profile curvature, and plan curvature).

Jasiewicz and Stepinski [133] introduced a novel method for unsupervised classification and mapping of landforms from a DEM. This method involves the pattern recognition, not the differential geometry. The foundation of this idea is the concept of geomorphon (geomorphologic phonotypes), which is a simple ternary pattern that serves as an archetype of a particular terrain morphology. A finite number of 498 geomorphons constitute a comprehensive and exhaustive set of all possible morphological terrain types including standard elements of landscape as well as unfamiliar forms rarely found in natural terrestrial surfaces.

And one should remember that horizontal and vertical resolution of the elevation data used to present a terrain surface has a significant influence on the level of detail and the accuracy of portrayal of surface features and on the values of land-surface parameters that are computed from a DEM. One should first test out predictive efficiency for various DEM resolutions and neighbourhood sizes, and then objectively derive the most suitable resolution and search size [120].

Acknowledgments

The publication has been (partially) financed from the funds of the Leading National Research Centre (KNOW) received by the Centre for Polar Studies of the University of Silesia, Poland.