Calculating Topographic Ruggedness Index in ArcGIS Desktop?

Does anyone know how to calculate Topographic Ruggedness Index in ArcGIS Desktop without access to command line ArcInfo Workstation?

"The topographic ruggedness index (TRI) is a measurement developed by Riley, et al. (1999) to express the amount of elevation difference between adjacent cells of a digital elevation grid. The process essentially calculates the difference in elevation values from a center cell and the eight cells immediately surrounding it. Then it squares each of the eight elevation difference values to make them all positive and averages the squares. The topographic ruggedness index is then derived by taking the square root of this average, and corresponds to average elevation change between any point on a grid and it’s surrounding area." -- from an aml arcscript by Jeffrey Evans

I would recommend to look outside ArcGIS) Very easy using the free gdal software: http://www.gdal.org/gdaldem.html

gdaldem TRI input_dem output_TRI_map

Or if you'd prefer it in saga gis: http://www.saga-gis.org/saga_modules_doc/ta_morphometry/ta_morphometry_16.html

• +1 I always appreciate seeing non-ArcGIS solutions to ArcGIS problems :-). This is a matter of principle, not antagonism towards ArcGIS in particular. One should avoid being locked in to a single software solution: not only is it professionally risky, it's intellectually stifling. – whuber Feb 11 '11 at 15:03
• I know I asked for an arcgis specific solution, but I'm accepting this one because of it's directness. GDAL Utilities are easy to acquire and install, universally acknowledged as best in class, and the command to generate this particular product is the definition of simplicity. – matt wilkie Feb 14 '11 at 17:45

Let's do a little (just a little) algebra.

Let x be the value in the central square; let x_i, i = 1, .., 8 index the values in the neighboring squares; and let r be the topographic ruggedness index. This recipe says r^2 equals the sum of (x_i - x)^2. Two things we can compute easily are (i) the sum of the values in the neighborhood, equal to s = Sum{ x_i } + x; and (ii) the sum of squares of the values, equal to t = Sum{ x_i^2 } + x^2. (These are focal statistics for the original grid and for its square.)

Expanding the squares gives

r^2 = Sum{ (x_i - x)^2 }

= Sum{ x_i^2 + x^2 - 2*x*x_i }

= Sum{ x_i^2 } + 8*x^2 - 2*x*Sum{x_i}

= [Sum{ x_i^2 } + x^2] + 7*x^2 - 2*x*[Sum{ x_i } + x - x]

= t + 7*x^2 - 2*x*[Sum{ x_i } + x] + 2*x^2

= t + 9*x^2 - 2*x*s.

For example, consider a neighborhood

1 2 3
4 5 6
7 8 9

Here, x = 5, s = 1+2+...+9 = 45, and t = 1+4+9+...+81 = 285. Then

(1-5)^2 + (2-5)^2 + ... + (9-5)^2 = 16 + 9 + 4 + 1 + 1 + 4 + 9 + 16 = 60 = r^2

and the algebraic equivalence says

60 = r^2 = 285 + 9*5^2 -2*5*45 = 285 + 225 - 450 = 60, which checks.

The workflow therefore is:

Given a DEM.

• Compute s = Focal sum (over 3 x 3 square neighborhoods) of [DEM].

• Compute DEM2 = [DEM]*[DEM].

• Compute t = Focal sum (over 3 x 3 square neighborhoods) of [DEM2].

• Compute r2 = [t] + 9*[DEM2] - 2*[DEM]*[s].

Return r = Sqrt([r2]).

This consists of 9 grid operations in toto, all of which are fast. They are readily carried out in the raster calculator (ArcGIS 9.3 and earlier), the command line (all versions), and Model Builder (all versions).

BTW, this is not an "average elevation change" (because elevation changes can be positive and negative): it is a root mean square elevation change. It is not equal to the "topographic position index" described at http://arcscripts.esri.com/details.asp?dbid=14156 , which (according to the documentation) equals x - (s - x)/8. In the example above, the TPI equals 5 - (45-5)/8 = 0 whereas the TRI, as we saw, is Sqrt(60).

• Thank you Bill. I appreciate seeing the specifics of how a tool or operation works. From this anyone with a suitable investment of time and intellectual energy can build a new apparatus to carry this job out using the tools they have at hand. It's information like this which will make GIS.se a useful service in the long haul. – matt wilkie Feb 14 '11 at 17:50
• +1 Great explanation. I guess this means a steep but smooth surface could have a higher TRI than a flat but bumpy surface. – Kirk Kuykendall Feb 14 '11 at 21:39
• @Kirk That is correct. There are ways to remove the effect of the local slope in order to obtain an index of "relative" ruggedness if you like. Although I haven't worked out the details, I believe that subtracting some universal multiple of (c*a)^2 from r2--where c is the cellsize and a is the slope (as rise/run, not as an angle or percent)--ought to do the trick. – whuber Feb 20 '13 at 17:09
• @whuber as always your answers contain amazing amount of knowledge!! Just one question please: does this mean that it is not possible to calculate the TRI of the cells which are located at the very edges of the raster? Due to them not being surrounded by neighboring cells all around? – marco Sep 14 '16 at 19:28
• @marco The TRI can be estimated even at the boundary cells. As indicated in the question, it should be expressed as an average, rather than a sum, by dividing the values I give here by 9. At boundary cells the value "9" in the formula and in the denominator needs to be replaced by the number of non-null values within their 3X3 neighborhoods: 6 for edge cells, 4 for corner cells. A grid of such values can be obtained from the focal sum of the indicator grid of the original values (it has 1's at all non-NoData cells and 0's elsewhere). Use that grid in place of the constant "9" in the formulas. – whuber Sep 15 '16 at 14:17

The Riley et al., (1999) TRI is square root of the summed squared deviations. This is very close to unscaled variance. If you want an implementation of Riley's TRI then please follow the methodology outlined by @whuber (the methodology provided by @user3338736 generalized the metric to the maximum in the window and does not represent the cell by cell variation).

I have a variation of TRI in our Geomorphometry & Gradient Metrics ArcGIS Toolbox that is the variance of a specified window. I find this more flexible and justifiable. There are also some other surface configuration metrics including rugosity and dissection.

• thanks Jeffrey. For some reason that page is blank except for the title in Firefox, thought it's ok in Chrome; might be one of my extensions. I'm happy to report at least that the scripts work unchanged in 10.2.2 (ones I tested anyway). – matt wilkie Jun 10 '14 at 18:23

-Edit: the information below is incorrect. Please see the post by whuber explaining the correct process.....

TRI (Riley 1999) and TPI (Jenness 2002) are similar, but different.

To calculate TRI and TPI using ArcGIS 10.x...

Step 1: Use the Focal Statistics tool to make 2 new raster datasets from a DEM.

Raster 1 "MAX") Neighborhood: Rectangle, Height: 3, Width: 3, Units: Cell, Statistics type: Maximum

Raster 2 "MIN") Neighborhood: Rectangle, Height: 3, Width: 3, Units: Cell, Statistics type: Minimum

Step 2: Use the Raster Calculator to perform the following functions on the 2 raster datasets you just created.

For TRI: SquareRoot(Abs((Square("%MAX%") - Square("%MIN%"))))

For TPI: ("%Input DEM%" - "%MIN%") / ("%MAX%" - "%MIN%")

Here is sample Python code exported from a model I built for TRI....

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# script.py
# Created on: 2014-03-06 08:56:13.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: script <Input_raster> <TRI_Raster>
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Script arguments
Input_raster = arcpy.GetParameterAsText(0)

TRI_Raster = arcpy.GetParameterAsText(1)
if TRI_Raster == '#' or not TRI_Raster:
TRI_Raster = "C:\\Users\\Documents\\ArcGIS\\Default.gdb\\rastercalc1" # provide a default value if unspecified

# Local variables:
MIN = Input_raster
MAX = Input_raster

# Process: 3x3Max
arcpy.gp.FocalStatistics_sa(Input_raster, MAX, "Rectangle 3 3 CELL", "MAXIMUM", "DATA")

# Process: 3x3Min
arcpy.gp.FocalStatistics_sa(Input_raster, MIN, "Rectangle 3 3 CELL", "MINIMUM", "DATA")

# Process: Raster Calculator
arcpy.gp.RasterCalculator_sa("SquareRoot(Abs((Square(\"%MAX%\") - Square(\"%MIN%\"))))", TRI_Raster)
• This is not the TRI described in the question. In fact, it cannot be thought of as measuring "ruggedness" at all, because it changes when you merely shift the vertical datum. For instance, your TRI of a 3x3 neighborhood with values (1,2,...,9) would be sqrt(9^2-1^2)=8.9, but adding 100 to the values (which just changes the datum without changing the shape of the surface at all) gives sqrt(109^2-101^2)=41. – whuber Jun 10 '14 at 17:32

This sounds very much like the Topographic Position Index, a process I used recently for one of my projects. There's an ArcScript on the ESRI support page, a Topography toolbox on the ESRI Resource Center page, and some more info on the process on the Jenness Enterprises page.

• TPI is a very different metric than roughness. Please, let's not go down the road of using them interchangeably. I believe that Topographic Position Index is traditional calculated as [dem - focalmean(dem)]. – Jeffrey Evans Jun 10 '14 at 16:03