Terrain Roughness – 13 Ways

Topographic roughness (ruggedness) is defined differently depending on the calculation used. Topographic roughness may be based on standard deviation of slope, standard deviation of elevation, slope convexity, variability of plan convexity (contour curvature), or some other measure of topographic texture. Scale is important in any roughness analysis – both that of the DEM and that of the landscape features you intend to characterize. In many cases, starting with a smoothed DEM produces a more interpretable roughness raster (see Relative Topographic Position). Here are several ways to calculate topographic roughness.

1.) Relative Topographic Position
Relative topographic position (also called the Topographic Position Index) is a terrain ruggedness metric and a local elevation index (Jenness, 2002). Topographic position of each pixel is identified with respect to its local neighborhood, thus its relative position. The index is useful for identifying landscape patterns and boundaries that may correspond with rock type, dominant geomorphic process, soil characteristics, vegetation, or air drainage. Classify the final output raster into high, medium, low based on natural breaks or breaks that work for you. The index is applicable to bathymetric data, too (see refs).

Note: An elegant solution for TPI was recently developed by Arthur Crawford of ESRI. Here is a LINK to the ESRI Forum thread and an explanation of his method. It is an alternative to the steps listed immediately below.

1.) Create a minimum elevation raster using Spatial Analyst Tools > Neighborhood > Focal Statistics, input = DEM, output = “minDEM”, Rectangle, 10×10 cells, Min, check box for “Ignore NoData”.

2.) Create a maximum elevation raster using Spatial Analysis Tools > Neighborhood > Focal Statistics, input = DEM, output = “maxDEM”, Rectangle, 10×10 cells, Max, check box for “Ignore NoData”.

3.) Create a smoothed DEM using Spatial Analysis Tools > Neighborhood > Focal Statistics, input = DEM, output = “10×10”, Rectangle, 10×10 cells, Mean, check box for “Ignore NoData”.

* Smoothing the DEM often produces a more spatially coherent, interpretable roughness raster, but smoothing is not required. You may wish to vary the smoothing factor. Here I use 10-cell by 10-cell neighborhood, but 20×20 or 3×3 or 5×5 might prove better for you.

4.) In Raster Calculator, use buttons to enter the following equation. The output raster will contain values scaled between 0 and 1.

(“10×10” – “minDEM”) / (“maxDEM” – “minDEM”)

10×10 = name of smoothed elevation raster
minDEM = name of minimum elevation raster
maxDEM = name of maximum elevation raster

* Try displaying output in either a 5-class (red, pale red, white, pale blue, dark blue) or a 3-class (red, white, black) color scheme. Equal interval or natural breaks or custom breaks will work. Try something like 0.45, 0.55, and 1.00 and adjust from there.

Worked Example for Relative Topographic Position:
Iberian-Peninsula-Roughness-Example-Instructions_SRTM90m_Cooley_v2

roughness

Example roughness raster (low quality image).

 

 

 

 

 

 

 

 

 


2.) MAD Algorithm
S. Trevisani & M. Rocca (2015) Computers & Geosciences v. 81
MAD: Robust image texture analysis for applications in high-resolution geomorphometry
Source code on GitHub –> https://github.com/cageo/Trevisani-2015

S. Trevisani (2014)
MAD functions for surface/image texture characterization
Source code on GitHub –> https://github.com/strevisani/MADSurfaceTexture


3.) Standard Deviation of Elevation
Standard deviation of elevation is a measure of topographic roughness (Ascione et al., 2008).

1.) Create mean elevation raster, Spatial Analyst > Neighborhood > Focal Stats, input = DEM, output = “meanDEM”, Rectangle, 3×3 cell, Mean, check box for “Ignore NoData”.

2.) Create elevation range raster, Spatial Analyst > Neighborhood > Focal Stats, input = DEM, output = “rangeDEM”, Rectangle, 3×3 cell, Range, check box for “Ignore NoData”.

3.) In Spatial Analyst > Map Algebra > Raster Calculator enter the following equation. The values returned will scale between -1 and 1.

(“meanDEM” – “DEM”) / “rangeDEM”

meanDEM = Mean elevation raster
rangeDEM = raster containing range of elevation values
DEM = original elevation raster

 


 4.) Slope Variability
Best calculated using a slope raster and a relatively large roving window (>100m), slope variability (SV = Smax – Smin) is a measure of the “relief of slope” of a landscape (Ruszkiczay-Rudiger et al., 2009).

1.) Spatial Analyst Tools > Surface > Slope, using catchment DEM and percent for slope values.

2.) Spatial Analyst Tools > Neighborhood > Focal Statistics > Max, using slope raster as input and a relatively large roving window (>100x100m). Name output raster “Smax”.

4.) Spatial Analysis Tools > Neighborhood > Focal Statistics > Min, using slope raster and same size window as previous. Name output “Smin”.

5.) In Raster Calculator enter the following single-line equation. Brackets enclose name of raster.

“Smax” – “Smin”

SV = slope variability output raster that will be created
Smax = maximum slope value raster
Smin = minimum slope value raster

 


5.) Basin-scale Ruggedness

A simple index used to compare the relief of basins using streamlines and basin boundary polygons. Use segmented streamline shapefile as the input to drainage density calculation (see Drainage Density lesson). Use Spatial Analyst > Density > Line Density.

Rb = A / DD

Rb = Basin ruggedness index
A = Area of the basin within bounding polygon
DD = Drainage density

 


6.) 2D Area : 3D Area Ratio

The area of a sloping pixel is not the same size as one that lies flat. A sloped pixel occupies more area than a flat one. The planimetric-to-surface area ratio is an index of topographic roughness with has a wide range of uses, including habitat assessment to land use planning to geomorphology. The index is largely independent of scale.

In a study using a 90m SRTM DEM for the rugged states of Jammu and Kashmir, Rashid (2010) found the 3D and 2D areas differed by nearly 25% (296,513 km2 vs. 222,236 km2, respectively). 3D area is always equal to or greater than 2D area.

3D Area = 2D Area / cosine(Slope in degrees)

– First, find the 2D area of your watershed DEM:

Watershed Area = Number of pixels x Area of a single pixel

– Open Raster Calculator (ArcToolbox > Spatial Analyst > Map Algebra > Raster Calculator)…

– If your slope raster was created in percent rise (slope_per), you’ll need to convert it to degrees using Raster Calculator. Name the output raster slope_deg. Note that arctan is atan in Raster Calculator. The conversion is:

atan(“slope_per” / 100)

– Calculate the 3D area and name the output raster area_3D:

2D area / cos(“slope_deg”)

Example
2D Area of One Pixel: 28.40m x 28.40m pixel sloping 60 degrees = 806.56 m2
3D Area of Same Pixel: 806.56 m2 / cos(60) = 806.56/0.5 = 1613.12 m2

– The 3D area is the Maximum value of area_3D, found in summary statistics (Properties > Symbology). The area ratio between 3D and 2D will be something like 1.3, or a number in that ballpark, depending the topography of your particular watershed.

Berry (2007) walks us through an example starting with percent rise slope:
Given: Slope of 20%
Find Slope Ratio: 20 / 100 = 0.20
Find Slope Angle: arctan(0.20) = 11.30993247 Degrees
Find Surface Area Factor: cos(11.30993247) = 0.980580676
Surface Area: 1 ha / 0.980580676 = 1.019803903

Caution: While 3D surface area is in theory “true” area, the data from which it is derived (gridded elevation data) introduces error. DEMs do not resolve all the humps and hollows in a real landscape. DEM resolution matters (see referenced articles) and will affect your calculations. Expect calculations made for very small regions (<264 pixels) to yield erroneous results (Jenness, 2004).

 


7.) Melton Ruggedness Number
The Melton Ruggedness Number is the basin relief divided by the square root of basin area (Melton, 1965), or:

(Zmax-Zmin) / Sqrt(A)

Zmax = Maximum elevation of basin
Zmin = Minimum elevation of basin
A = Basin area in m2

Preliminary notes…need to check and revise Melton method:
– Zonal summary of the elevation by watershed. Summary will include min and max values.
– Find watershed area.
– Calculate tabular data: MRN = (ZMax – ZMin) / Sqrt(Area).
– Join the watershed zonal summary table to the watershed attribute table.
– Subtract your elevation grid from the resulting grid of “MaxZ”.
– If you want grids instead of table output, use the zonalmin and zonalmax functions.
– In Raster Calculator, create this statement:
(zonalmax(watersheds, dem) – zonalmin(watersheds, dem)) / sqrt(ZONALAREA(watersheds))

 


8.) TRI (Riley)
Terrain Ruggedness Index (TRI) is the difference between the value of a cell and the mean of an 8-cell neighborhood of surrounding cells. First create the two input neighborhood rasters from a DEM (use a 3×3 neighborhood for min and max), then run the equation in Raster Calculator. Syntax matters. Classify the resulting ruggedness index values using the categories of Riley et al. (1999). If not the life of Riley, at least his landscape. Your particular DEM may not yield the full range of values possible (fewer ruggedness categories may result).

** The Riley formula is currently under review. See this article and erratum: LINK **

TRI Categories
LEVEL = 0-80m
NEARLY LEVEL = 81-116m
SLIGHTLY RUGGED = 117-161m
INTERMEDIATELY RUGGED = 162-239m
MODERATELY RUGGED = 240-497m
HIGHLY RUGGED = 498-958m
EXTREMELY RUGGED = 959-4367m

ArcGIS Instructions
– ArcToolbox > Spatial Analyst Tools > Neighborhood > Focal Statistics:
Input = DEM, Neighborhood = Rectangle, Size = 3×3, Units = Cells, Statistics Type = Minimum, Output filename = 3x3min

– ArcToolbox > Spatial Analyst Tools > Neighborhood > Focal Statistics:
Input = DEM, Neighborhood = Rectangle, Size = 3×3, Units = Cells, Statistics Type = Maximum, Output filename = 3x3max

– Spatial Analyst > Map Algebra > Raster Calculator:
SquareRoot(Abs((Square(“3x3max”) – Square(“3x3min”))))

– Categorize the resulting values according to Riley’s seven TRI categories. Classify visually (Properties > Symbology tab) using 7 classes and Break Values at 80, 116, 161, 239, 497, 958, 4367 (or whatever the maximum value of the raster is). Classify more formally using the Reclassify tool and same class breaks (Spatial Analyst > Reclass > Reclassify tool).


9.) TRI (Nellemann’s Terrain Roughness Index)
Nelleman’s is a somewhat antiquated contour density (transect-and-contour map) approach with applications to arctic wildlife. See Nelleman et al. (2007), Nellemann and Fry (1995), Nelleman and Thomsen (1994) papers for methods. Nelleman et al. (2007), a paper on brown bears, classified TRI values for a study on Scandinavian brown bears into “Rugged”, “Flat”. They used a 1:100,000 scale DEM with a 10m contour interval, 4km transects within 4km x 4km grid cells: “Rugged” = TRI > 2.5, “Flat” = TRI < 2.5.

I am not going to provide a detailed explanation of or ArcGIS instructions for TRI here, but…

– You could place points at locations where contours intersect your transect line (Analysis toolbox > Overlay > Intersect tool, set output type to Point).

– You might also create an Aspect raster, turn your transect line into a set of evenly-spaced points, extract the aspect raster value associated with each (Spatial Analyst tools > Extraction > Extract Values to Points; append raster values to the output point attribute table), and summarize aspect fluctuations along the transect in a chart.

– You also might create a polygon shapefile sample grid (each grid cell is a Zone) and use Zonal Statistics to extract the total length of contours within each grid cell, then divide by area to get contour density.

TRI = (TNC x TNF) / (TNC + TNF)

TNC = Number of contour intercepts along transect
TNF = Number of aspect fluctuations along transect

 


10.) RIX Site Ruggedness Index
With application to wind power projects, the Site Ruggedness Index (RIX) was developed for use with Wind Atlas Analysis and Application Model (WAsP) (Mortensen et al., 2007). The RIX index is the fractional extent of the surrounding terrain steeper than a critical slope.

Parameters used by Bowen and Mortensen (1996, 2004):
Critical slope = 30% calculated along vectors radial to tower
Analysis radius = 3.5 km from wind generation tower

* Note: One way to simplify the method would be to substitute raster-generated slope values for those calculated along radial vectors from the tower position, then apply the same or slightly adjusted critical cut off. Raster slopes would always be equal to or slightly steeper than radially-calculated slopes.

 


11.) Standard Deviation of Slope
Use ArcToolbox > Spatial Analyst Tools > Neighborhood > Focal Statistics tool to create a 3×3 or 5×5 cell neighborhood with Standard Deviation as the operator. See Grohmann et al. (2011) for a comparison of several roughness techniques.

 


12.) 1st Derivative of Flow Direction
If you are looking for a qualitative way of revealing the “topographic grain” of a landscape, run the Flow Direction tool on a filled DEM, then run the Slope tool using FlowDir as the input. In theory, this would highlight the change in flow direction. The effect is to emphasize subtle features (like faults, gullies, ridgelines). I’ve found good success by classifying the display into 5 classes, equal interval, with this color scheme: White/Light blue/Black/Yellow/Red. The raster values are meaningless, but the layer may prove visually useful for reconnaissance.

 


13.) Hobson’s Surface Roughness Factor

Hobson (1972) introduced the Surface Roughness Factor (SRF). Xi, Yi, and Zi are vectors normal to the land surface. I have never used this. If you have, let me know how it works and what values are returned.

SRF = √ (∑ni=1 Xi)2 + (∑ni=1 Yi)2 + (∑ni=1 Zi)2 / n

Xi = sin(s) * cos(a)
Yi = sin(s) * sin(a)
Zi = cos(s)
n = number of cells in the analysis window (sample size)

 


Refs:
Ascione et al. (2008)
Berry (2007) Beyond Mapping III: Procedures and applications in GIS modeling, Berry & Associates
Bertuzzi et al. (1990) Soil Tillage and Research 17, p. 87-99
Bowen and Mortensen (2004), Riso National Lab report (Denmark)
Brasington et al. (2012) Water Resources Research 48
Brown and Hugenholtz (2013) Geosphere 9, p. 367-377
Brozovic et al. (1995)
Burbank & Anderson (2012, p. 307) Tectonic Geomorphology textbook
Chimi-Chiadjeu et al. (2013 in press) Computers and Geosciences
Cord et al. (2007) Icarus 191
Eitel et al. (2011) Catena 87
ESRI Forum Post: http://forums.esri.com/Thread.asp?c=93&f=995&t=227022
Felicisimo (1994a) Modelos digitales del terrano, Introduccion y aplicaciones en las cincias ambientales, Penalfa Ediciones, 222 pgs.
Frankel and Dolan (2007) Journal of Geophysical Research 112
Gillies et al. (2007) Boundary Layer Meteorology 122
Glenn et al. (2006) Geomorphology 73
Grohmann et al. (2011) IEEE Transactions on Geosciences and Remote Sensing 49 LINK
Guillobez (1998) Soil Tillage and Research 45, p. 419-432
Guisan et al. (1999) Plant Ecology 143
Haubrock et al. (2009) Catena 79
Henning and Radke (2008) Journal of Photogrammetry and Remote Sensing 63
Heritage and Hetherington (2005) in Walling and Horowitz (eds) IAHS-AISH v. 291
Helming et al. (1993) Soil Technology 6, p. 273-286
Heritage and Milan (2009) Geomorphology 113
Hobson (1972) in Chorley (editor), Spatial Analysis in Geomorphology
Huang and Bradford (1990) Soil Science Society of America 54
Jenness (2004) Wildlife Society Bulletin 32
Jeff Jenness/Jenness Enterprises (www.jenness.com)
Jester and Klik (2005) Catena 64
Jones et al. (2000) Environmental Monitoring and Assessment 64
Lundblad et al. (2006) Marine Geodesy 26
Mazzarini et al. (2008) Annals of Geophysics 51
Melton (1965) Journal of Geology 73
McKean and Roering (2004) Geomorphology 57
Morris et al. (2008) Journal of Geophyscial Research 113
Mortenson et al. (2007) <– Google search: “WAsP wind model”
Nellemann et al. (2007) Biological Conservation 138
Nellemann and Fry (1995), Arctic 48
Nellemann and Thomsen (1994), Arctic 47
Polleyea and Fairley (2011) Geosphere 8
Rashid (2010) Journal of Geomatics 4
Raupach (1992) Boundary Layer Meteorology 60
Rinehart et al. (2004)
Riley et al. (1999) Intermountain Journal of Sciences 5
Riley method explained for QGIS, LINK
Ruszkiczay-Rudiger et al. (2009)
Sankey et al. (2011) Geomorphology 135
Schmid et al. (2005) Remote Sensing and Spatial Information Sciences 36
Shepard et al. (2001) Journal of Geophysical Research 106 p. 32777-32795
Soulard et al. (2013) Earth Surface Processes and Landforms 38
Trevisani & Rocca (2015) Computers & Geosciences April 2015
Weiss (2001) Poster at 21st ESRI Int’l User Conference
Wilson et al. (2007), Marine Geodesy 30

Comments are closed.