Geoinformatics & Geostatistics: An OverviewISSN: 2327-4581

Reach Us +1 850 900 2634

Research Article, Geoinfor Geostat An Overview Vol: 5 Issue: 2

Slope-Independent Landscape Roughness Attribute Provided by Measurement of Local Contour Line Density

Veronica Ochoa-Tejeda1* and Jean-François Parrot2

1SERTIT (Service Régional de Traitement d’Images et de Télédétection), Université de Strasbourg, Boulevard Sébastien Brant, 67412 Illkirch, France

2Instituto de Geografia, Universidad Nacional Autonoma de Mexico, Circuito Exterior, Ciudad Universitaria, Coyoacan, 04510, Mexico

Corresponding author : Veronica Ochoa-Tejeda
SERTIT, Université de Strasbourg, Boulevard Sébastien Brant, 67412 Illkirch, France
+ 33 (0) 786011897
E-mail: [email protected], [email protected], [email protected]

Received: February 02, 2017 Accepted: February 16, 2017 Published: February 21, 2017

Citation: Ochoa-Tejeda V, Parrot J (2017) Slope-Independent Landscape Roughness Attribute Provided by Measurement of Local Contour Line Density. Geoinfor Geostat: An Overview 5:3. doi: 10.4172/2327-4581.1000159


Various digital treatments have been proposed to describe the digital elevation model (DEM) surface roughness, depending on the geomorphologic purpose. Some of these treatments use the drainage network pattern directly or combine the former binary information with DEM data. Others take into account local characteristics of the DEM’s surface, such as the deviation of the perpendicular issued from each pixel or the computation of the surface curvature. The aim of the present study is to propose new geomorphic algorithms based on contour line density to calculate the DEM surface roughness. The contour line density is calculated within a moving window of n × n pixels, and the contour lines are extracted according to a given range of hypsometric intervals. The value of this total length inside the moving window reflects the local irregularities encountered in the studied zone. Meanwhile the slope-independent roughness parameter (SIRP) measures the total length of a limited number of contour lines inside a moving window, extracted according to the total hypsometric difference between the upper and lower altitudes encountered in the window. While the first calculation emphasizes the general features of the studied DEM in terms of its local surface roughness, the SIRP provides a roughness measurement independent of the slope. In the Sierra Norte de Puebla region (Central Mexico) where numerous landslides have occurred recently, slope instability is closely related with the smooth zones extracted by use of the SIRP attribute.

Keywords: DEM’s surface roughness; Contour line density; Slope independent roughness parameter; Geomorphic features; Sierra Norte de Puebla


Raster digital elevation models (DEMs) are nowadays widely used, since advances in computing technology and storage have increased strongly in recent years. The horizontal and vertical resolutions allow accurate calculation of parameters extracted from a DEM surface. Digital terrain analysis [1] allows definition of primary attributes calculated directly from the DEM surface (i.e. slope, aspect, plan and profile curvature) and secondary attributes whose calculation requires external information and that are more devoted to estimating the role of topography in the distribution of soil water or in the susceptibility of landscapes to erosion, for instance. Most of the primary attributes are computed directly with a finite-difference scheme or by fitting an interpolation function Z=f (x,y) to the DEM in order to calculate the derivatives of this function [2-5]. These attributes provide numerous geomorphologic indicators concerning the region under study [6]. They are used to describe the morphometry, catchment position, stream channels and other factors and to compute topographic attributes [7-12]. Several recent pedologic and geomorphologic applications are based on features automatically extracted from DEMs. The publications concerning this topic were surveyed by Pike [13]. Various softwares described in Gallant and Wilson [14] has been developed for the parametric study of gridded DEM surfaces.

Geomorphometric features such as the surface roughness of a DEM provide information about regional geomorphology, and many roughness parameters have been proposed. For instance, the fractal dimension can be considered as an indicator for terrain surface roughness in DEMs [15-18]. Other authors used parameters of curvature to analyse the surface locally [19-22]. Curvature attributes are useful for measuring flow convergence and divergence [23] and for defining geomorphic units [8], and they can be calculated by using a bi-quadratic surface to classify the types of surface [20,22,24].

Moreover, Teng and Tay [25] suggested defining a simple and easily computed surface roughness index RI based on the differences of heights between neighbouring pixels within a window. Dinesh [26] proposed a procedure to compute a scale-independent surface roughness parameter from DEMs and extended this method to compute the surface roughness of individual mountain objects extracted from DEMs by means of mathematical morphology [27]. With a similar purpose, Hebeler and Purves [28] assumed that topographic indices calculated by the software GLOBE [29] such as Altitude Roughness (standard deviation of altitude in a 3 × 3 neighbourhood) and the Slope Roughness (standard deviation of slope in a 3×3 neighbourhood) are sufficient to study the influence of elevation uncertainty on derivation of these indices. Gallay et al. [30] analysed the elevation error of high-resolution DEMs by using the roughness index defined by Fisher and Tate [31]. In this case, the area ratio roughness would correspond to the ratio of real surface area to the area of its orthogonal projection [32,33].

However, for most of these treatments, the expression of the surface roughness is too closely related to the slope. Therefore, we tried to develop a surface roughness measurement that is independent of the slope values. The calculation of this new attribute is based on the measurement of the total length of a given number of contour lines extracted inside a moving window of nn pixels. As the number of extracted contour lines remains unchanged, the range of the hypsometric intervals varies according to the slope value, i.e. the total length of the contour lines depends only on the design they draw according to the roughness of the moving window surface.

To illustrate the geomorphologic information provided by this new parameter, the treatment has been applied to the region of La Soledad, Sierra Norte de Puebla, Central Mexico, where numerous mass movements occurred in 1999 and 2005. A close relationship occurs between the position of the landslide starting points and the smooth zones supplied by the new roughness attribute called the Slope-Independent Roughness Parameter (SIRP).

Ochoa-Tejeda [34] studied the distribution of all the landslides that occurred in this region in 1999 and measured the relative influence of various factors related to mass movements: geological units, geomorphologic parameters and primary attributes such as slope, aspect, convexity, energy of the relief, drainage density and vegetation density [35] as well as the local 3D fractal dimension [18]. Another approach related the morphology of the landslide traces with the nature of the material remobilized during the event [36].

In 1999, Tropical Depression N° 11 coming from the Gulf of Mexico triggered more than 3000 mass movements related to the torrential and sudden rainfalls of 4-5 October [37]. Data from the Mexican National Meteorological Service indicate that the amount of rain during the first days of October considerably exceeded the monthly and annual average values. For instance, in Zacapoaxtla, the accumulated precipitation during the 4 days of 3-6 October was 844 mm, i.e. 60% of the annual average of 1421.2 mm.

In 2005, from 1 to 4 October, hurricane Stan produced intense rains in the northern region of Puebla with a peak of 273 mm during 24 h. In addition, Tropical Depression N° 40 interacted with the hurricane from 4 to 6 October causing rains of up to 305 mm in 24 h. This phenomenon was responsible for severe soil saturation and water flows [38]. The effect of rain depends basically on the intensity, duration and distribution of the storm. Manzini and Rabufetti indicate that when the superficial formations are constituted by clastic and colluvial soils, the precipitation threshold that produces shallow landslides is related to the inclination of the slope, the infiltration and the loss of apparent cohesion.

The results obtained in this test region show that the SIRP clearly emphasizes the different roughness zones according to the size of the moving window.

Material and Methods

Derivation of the roughness attribute

To calculate the contour-line density inside a moving window of size n×n, it is first necessary to take out, one by one, each contour line extracted from the DEM according to the hypsometric interval retained by the user. At each step, all the pixels that are lower than or equal to a given elevation are reported in a temporary image. Then a treatment in this matrix ensures the extraction of the contour lines from the corresponding zone. A plain procedure consists of considering that a pixel defines a contour line (Figure 1A) when at least one of the cardinal neighbouring pixels does not belong to its corresponding extracted zone (Figure 1A). When necessary, the remaining corners are eliminated in order to obtain an 8-vicinity curve, i.e. a curve formed by segments of 1 to n pixels linked only by their corners.

Figure 1: Example of contour line extraction.

The length of this contour line is calculated as follows. One can consider, within the pixel surface, the value of one-half of the straight line that links the pixel centre and the centre of the neighbouring pixels. In this case (Figure 2), and according to the configuration encountered in a test 3×3 pixel window, the two lines values dz1 and dz2 encountered within the pixel surface are respectively equal to ½ and to of the cell size (Ps). The value dz1 is applied to the cardinal neighbours and the value dz2 is used when the pixels are diagonal neighbours.

Figure 2: Length calculation within the pixel surface.

Taking into account the cell size Ps, the algorithm calculates the lengths lz1 and lz2 encountered in the 3 × 3 pixel window. The computation is done in this way:

image (1)

Where lz1 = dz1× Ps and lz2 = dz2×Ps

When strong escarpments are present, it is not rare for a pixel to register more than one contour line, since for a low hypsometric interval value successive contour lines may pass through the same point. The total length Lz corresponds to the sum of all the local values Lz obtained. According to the raw order, this calculation is done for each pixel of each extracted contour line in order to avoid any coalescence between contour lines such as they appear in Figure 1. The value obtained is reported in a float matrix and eventually summed to a pre-existing value.

In the example of the three contour lines reported in Figure 1, Lz is 29.796. For instance, if the value of the pixel side Ps is 25 m, the total length for the three contour lines encountered in the moving window is 744.9 m.

The size of the moving window side n corresponds to an owe value (i.e. 11 × 11 pixels in this example). The surface S where the three curves are emplaced (Figure 1) is equal to n-1×n-1, (here 10 × 10 pixels), because when a pixel is on the window border there is no information about the pixels neighbouring emplaced throughout this window. In this case, the surface S is equal to:

image (2)

For the pixel side Ps formerly proposed, S corresponds to 62 500 m2 in such a way that the Density of the Contour Lines in the moving window is equal to:

image (3)

For example, Cd is 1.146 in the case of the contour lines reported in the Figure 1. At the end of the treatment, normalization (0, 100) can be operated in order to obtain the Cd percentage reported in an image of the Contour Line Density (CLD).

In order to compare the results obtained using the calculation of the contour line length inside a moving window on a training zone with the slope values, the mean slope angle Msa has been calculated in a moving window that has the same size:

image (4)

where n is the number of the pixels encountered in the moving window and A the angular value of each pixel in this window.

Slope-independent approach

As the CLD is related mainly to the terrain slope, another approach has been developed. The calculation of the “Slope-Independent Roughness Parameter” or SIRP consists of applying the former algorithm inside a n × n moving window taking into account the local difference observed between the upper and lower altitudes inside the moving window. Although this calculation is computationally intensive, its results are more consistent with the “roughness” notion. The upper and lower hypsometric values (Hmax and Hmin) observed inside the moving window allow calculation of the local hypsometric interval H used to extract the contour lines dividing the difference between Hmax and Hmin by ξ+1, where ξ corresponds to the number of required contour lines.

image (5)

The value of I allows the same number of contour lines to be drawn in each successive moving window.

Following the procedure described above (Figure 1) all the pixels that are lower than or equal to the succeeding altitude levels are reported in a temporary n×n matrix; the pixels describing the perimeter of the consecutive slices thus obtained correspond to a contour line. The length of the line crossing each pixel belonging to the contour line is calculated according to the formula (1) and the total length Lz of this contour line corresponds to the sum of all the local values lz obtained. As previously, the surface where the three curves are emplaced is equal to n-1 × n-1, and the SIRP is calculated according to the equation (3).

Case Study

Validation and control of the results used a DEM of one region from the Sierra Norte de Puebla (Central Mexico) that has a surface area of ~135 km2. This region was affected by torrential rainfalls in October 1999, producing numerous landslides and flooding [39-41]. In the whole region, about 3000 mass movements occurred, consisting of rock and soil slides and slips, debris flows and avalanches. Capra et al. [41] identified two types of mass movements in the Teziutlan volcanic area: a) small detrital fans especially in the ignimbrite formations; and b) soil slide/debris flow in the volcanic sequence. In the region of Zapotitlán de Méndez [42] rotational landslides affected folded shales. The study region (19° 53’-20° 00’ N and 97° 24’-97° 25’ W) lies between the Sierra Madre Oriental geological province and the northern border of the oriental part of the Trans- Mexican Volcanic Belt (Figure 3). The main geological units [43] are the Paleozoic formations mainly formed by mylonitized schist, discordant Jurassic formations (sandstones and marls in the lower sequences and limestone in the Upper Jurassic), Cretaceous limestone of Neocomian age and Tertiary and Quaternary volcanic flows (Pliocene andesites and Quaternary ignimbrites). With the exception of these volcanic formations, numerous necks, dykes and sills of rhyolitic composition intrude throughout the stratigraphic sequence.

Figure 3: Localization of the study region on the geological sketch map of Mexico.

The DEM results from an interpolation [44,45] taking into account the contour lines digitized from the INEGI topographic map of the region (1/50 000). The resulting DEM has 2374 rows and 2270 columns (Figure 4). The pixel size is 5 m and the elevation is in decimetres. Even if this DEM presents some artefacts inherent to this type of interpolation that over-estimates the number of pixels describing the contour lines, a light smoothing blurs these defaults because the hypsometric scale is in centimetres and the result is good enough to calculate in an accurate way the indices described above.

Figure 4: Shadowed DEM of the study region.

Results and Discussion

General contour line density

The result does not present the artifacts observed in the case of different primary attributes such as, for example see [1] the “deviation from mean elevation” (formula 5, 6 and 7) that is supposed to emphasize the local surface roughness or the “elevation range” (formula 8) that produces visible overlapping circular regions (Figure 5). These attributes are respectively calculated as follows where the set of cells included in the circle about the point being considered is denoted by C, and the number of cells in the circle is nC.

Figure 5: Results of applying A. the “range elevation” attribute inside a circle (radius 50 pixels) and B. the “contour line density” index calculated by using the same radius value.

image (6)

Where the Standard deviation SD and the mean ͞Z (Z corresponds to the altitude inside the circle) are calculated as:

image (7)

image (8)

image (9)

But even if the resulting image underlines some geological features, the CLD remains too closely related to the slope values. In fact, on the scatter plot CLD versus Slope (Figure 6), the normal trend has a high correlation coefficient of 0.7319.

Figure 6: Scatter plot CLD versus Slope

Slope-independent roughness parameter

The calculation that takes into account the local hypsometry inside the moving window offers a completely different result. As the number of contour lines extracted inside the moving window depends on the hypsometric difference between the upper and lower altitude, and remains always equal to a given number (3 for example), the contour line complexity is related only to the roughness of the DEM surface observed inside the moving window. For instance, a flat zone will have the same total length Lz and the same contour line density regardless of the slope value.

Figure 7 illustrates the result obtained using 11×11 moving window size. The scatter plot “Slope Independent Roughness Parameter [SIRP]” versus Slope (Figure 8) clearly shows that the measurement is not related to slope values. The corresponding coefficient of correlation is in this case –0.043.

Figure 7: Slope-independent roughness parameter” calculated in a 11 × 11 moving window.

Figure 8: Scatter plot SIRP versus Slope.


The SIRP is a new morphometric parameter directly issued from the DEM surface. As shown above, this primary attribute is not related to the corresponding local slope and can be used to provide information about the different types of landslides observed in the study region.

In Figure 9, the grey tone (from light to dark) shows the slope value (11°-20°, 21°-30°, 31°-40°, >41°) inside the segment 14-25 of the SIRP histogram (Figure 10). Lower values of the SIRP histogram are essentially associated with the La Soledad Lake.

Figure 9: Relationships between smooth zones and landslides. Filled squares correspond to the starting point of the 1999 landslides(142) and empty squares to the 2005 landslides (97). Inside the smoothed zones, dark grey tones correspond to slopes greater than 41° and the lightest grey tones to slopes of 11- 20°.

Figure 10: Distribution of the landslide starting points according to SIRP values. A. Global: 1= rough zones, 2=smooth zones. B. Distribution taking into account the slope range inside the smooth zones: 1=slope 11-20°, 2=slope 21-30°, 3=slope 31-40°, 4=slope>41°.

A close relationship exists between the low values of the SIRP and the landslides distribution. More than 80% of the 1999 landslides and 75% of the 2005 landslides are in smooth zones (Figure 10). Hence, roughness might somehow represent an obstacle to the triggering of shallow landslides. The landslides related to very rough zones represent 5.75% in 1999 and 5.20% in 2005. They mainly affected the rhyolites in 1999 and the Upper Jurassic cliffs in 2005, producing rock avalanches and debris flows.

Within the smooth zones, the distribution of landslide starting points in relation to slope values (Figure 10) was similar between 1999 and 2005 except for the lowest slope vales, where the percentage was greater for 2005. This slight difference may be related to the nature of the tropical depressions that affected the mountainous zone of the state of Puebla, mainly in its northern region, in October 1999 and October 2005, as well as the different types of material involved.

According to these observations, it seems that soil saturation in October 2005 was responsible for the development of mass movements affecting essentially superficial formations. This observation is corroborated by the landslide distribution according to the geological units where the increased percentage in 2005 concerns the Quaternary ignimbrites and the superficial formations related to the Middle Jurassic and the metamorphic formations.


The aim was to define a roughness parameter derived from the gridded surface of the DEM in order to obtain morphometric information concerning geomorphologic studies. The Slope- Independent Roughness Parameter is a new primary attribute that is not closely related to the slope and thus represents a powerful indicator of the DEM surface roughness.

Its application to the distribution of mass movements that occurred in the Sierra Norte de Puebla (Mexico) during October 1999 and October 2005 shows that the superficial and shallow landslides are closely related to the smooth zones. They are essentially debris flows that are generated when hillside colluvium or landslide material becomes rapidly saturated, as well as mudflows, creeps and eventually lateral spreads. Falls involving rocks and boulders detached from steep slopes or cliffs, and rock avalanches occurred in zones of high roughness because the mobilization of the material is mainly related to discontinuities such as fractures, joints and bedding planes. They produce coarse deposits in talus and talus breccias and sometimes colluvium easily remobilized by torrential rainfalls.

The SIRP is a useful parameter that can be used in tandem with other pertinent parameters such as slope angle, aspect or concavity in order to map the susceptibility of hillslopes to mass movements.


Track Your Manuscript

Share This Page