# Correcting for local topography in calculating TRI from DEM in ArcGIS

I'm currently investigating the use of different equipment in calculating TRI from elevation data using ArcGIS. From the methods available (http://gis4geomorphology.com/roughness-topographic-position/) the ruggedness calculated is influenced by overall elevation, i.e. a slope can be smooth but have a higher TRI than a bumpy flat surface.

Is there a way to correct for this influence of local topography? I'd ideally like to create a relative ruggedness index where the raster layer is purely based on ruggedness rather than elevation?

I would like to calculate the ruggedness of a terrain, specifically a volcano with a solidified lava flow. In theory (photographs of the site) the lava flow is considered rougher than the volcano. However using the methods of calculating ruggedness the volcano is considered rough and the lava flow considered smooth. I would like to remove this local topography to calculate the roughness based on, essentially, the geology alone.

• The ruggedness is not influenced by "overall elevation" at all, as the equations clearly show (they are based entirely on elevation differences). However, this "ruggedness" is not a local measure relative to a fitted plane, so a nonzero slope will contribute positively to the ruggedness. This is easily adjusted--but then you will not be computing the TRI, you will be computing something else altogether. What we need to know, then, is what exactly do you propose this new "ruggedness" should be? What physical properties of the landscape is your personal ruggedness intended to quantify? – whuber Feb 3 '15 at 16:52
• @whuber I've added some additional information of the context of the assignment. You mentioned on a previous question: gis.stackexchange.com/questions/6056/… - that you had an idea of removing local slope to calculate relative ruggedness - I believe this is what value I am after. – babbitt Feb 3 '15 at 17:11
• We are working on a paper detailing a hierarchically decomposed metric that is standardized by the variance in the previous scale(s). The intent of the metric is to account for bias associated with rapidly changing surface variation (eg., mountainous to large valley). Currently we are mucking around with a nonparametric version but have stumbled across an ESRI bug, in FocalStatistics, that is causing notable issues. We should have the parametric version of the model available in the next release of our "Geomorphometry & Gradient Metrics ArcGIS Toolbox". – Jeffrey Evans Feb 3 '15 at 17:59
• @Jeffrey That must be a relatively new bug. Could you indicate what it is and which stats it affects? – whuber Feb 3 '15 at 18:07
• We found that for large window sizes (eg., 27x27), using the median statistic, the function drops cells past a certain number of cells in the window, almost like there is a cell limitation to the median statistic. It is now a registered ESRI bug. This is at version 10.2.1. – Jeffrey Evans Feb 3 '15 at 18:12

The question (as clarified in a comment) asks how to

remove local slope to calculate relative ruggedness.

There is a simple way to do this. It relies on computing the slope using the same local data as the ruggedness (which usually is a 3 by 3 square neighborhood). I recall verifying that ArcGIS computes slope (s) and aspect in exactly this manner: its algorithm is equivalent to fitting a plane to the nine neighborhood values using ordinary least squares and then reporting the slope and aspect of that plane.

We can therefore exploit results of least squares fitting to conclude that the following is the solution for a grid of cellsize c:

To obtain the local ruggedness, subtract 6 * c^2 * s^2 from the topographic ruggedness.

(The constant "6" is derived from 9--the amount of local data--minus 3, the number of parameters used to describe the local plane (a mean elevation, a slope, and an aspect). This formula can be extended to other neighborhoods simply by replacing "6" by three less than the number of cells within those neighborhoods, assuming a common neighborhood is used for the topographic ruggedness, the mean elevation, and the slope calculations.)

The "local ruggedness" in this sense is the sum of squares of deviations between the elevation data and the corresponding elevations of the locally fitted plane, whereas the topographic ruggedness is the sum of squares of deviations between the data and the average elevation within the neighborhood.

(Using sums in this manner produces values that depend on the cellsize and neighborhood dimensions. These arbitrary elements of the calculation make it somewhat less than scientific or meaningful. It would be better to compute mean squared deviations and to standardize the radius of the neighborhood that is used. But those issues should be discussed elsewhere; here let us just accept that the sums are desired and recognize that it is straightforward to convert such sums into mean values.)

### Example

Let the cellsize be 3 meters and suppose the elevations in a neighborhood, also in meters, are

``````55 58 57
54 56 51
52 53 59
``````

These average is 55 + 58 + ... + 59 = 55. The elevations relative to this average are

`````` 0  3  2
-1  1 -4
-3 -2  4
``````

Its "topographic ruggedness" is the sum of the squares of these residuals, equal to (0^2 + 3^2 + ... + 4^2) = 60.

The least squares plane has these fitted values:

``````55 56 57
54 55 56
53 54 55
``````

The differences between the data and this fit are

`````` 0  2  0
0  1 -5
-1 -1  4
``````

The local ruggedness is the sum of their squares, equal to 0^2 + 2^2 + ... + 4^2 = 48. It did not have to be separately computed, though, because the squared slope is 2/9. (The GIS will compute this. As a check of the example, look at the fit, which obviously tilts from upper right to lower left across the diagonal line of 55's. It drops by 2 meters across the diagonal of each cell: from 57 down to 55 down to 53 as you go from upper right to lower left. The Pythagorean Theorem tells us the square of the diagonal is 3^2 + 3^2. The squared slope consequently equals 2^2 / (3^2 + 3^2) = 2/9.)

Therefore, all we have to do is calculate

48 = 60 - 6 * 3^2 * 2/9

### Things to watch out for

It is important to compute the slope correctly. As used here, it is the usual "rise/run" ratio--not an angle. Either compute the slope in percent and divide by 100, or compute the tangent of the slope-as-angle. Make sure to use the same units of measure for the positional coordinates and the elevations.

Using a highly accurate projection, although desirable, is not necessary: the distortions induced by projections are locally linear and that will not change the planar fit because the plane is given by a linear function of the coordinates. The formula for the plane will be wrong and the slope will be wrong, but everything cancels out! Accordingly, it is probably best to work with the original data rather than reprojecting them, thereby eliminating resampling error. (Resampling almost always reduces ruggedness; the more often you resample, the more you bias the ruggedness estimates downwards.) You can then reproject the results in any way you like.

I should have mentioned that the original ruggedness calculation is described at https://gis.stackexchange.com/a/6059/664. It requires two focal sums and some arithmetic combinations ("map algebra").