I am not able to reproduce the terrain ruggedness index for the example used in the original Riley et. al paper.

The example is as follows:

enter image description here

When I try to reproduce the numbers in MS Excel I get different numbers.

The cell wise TRIs are calculated as follows maintaining the same grid order.

enter image description here

The TRIs for the whole terrain are also given below in the same order.

enter image description here

Any idea what am I missing?

  • 2
    You probably need to show us the formulas you used to calculate the TRIs
    – Ian Turton
    Jun 9, 2021 at 14:20

1 Answer 1


I hate the be the bearer of bad news but, Riley et al., (1999) is wrong. The ruggedness index itself is fine but, the example in the paper is incorrect. Obviously, the reviewers did not work through the math. I wrote the original TRI AML (Workstation Arc/Info) and caught this error shortly after the paper was published.

Functionally, TRI is representing a form of variance. Simply put, TRI is the square root of the summed squared deviations from the center cell whereas variance is the summed squared deviations from the mean. I have always though that variance is a more stable and interpretable measure of ruggedness. Also, variance is readily available as a focal function in most GIS software (via the squared standard deviation) whereas, TRI requires neighbor notation. Using standard focal functions allows one to easily derive/evaluate multiple spatial scales by simply changing the size of the kernel. If you need to exactly derive the published TRI, I have the original metric, along with option for a scaled version, available in the R package spatialEco.

Comparison of published results in Riley et al., (1999) in R

Fig 3-c published TRI = 100, correct TRI = 41.83. Here is the defined 3x3 matrix:

200  210  220
210  225  225
205  210  220

We can pull the center cell of the matrix, to use in calculating the deviations, by finding the center point of the unfolded matrix.

( x <- matrix(c(200,210,220,210,225,225,205,210,220), 3, 3) )
  ( x = as.vector(x) )
  i <- ceiling(length(x)/2) # position of center cell value
sqrt(sum((x[i] - x[-i])^2))

Here is the brute-force way that exactly replicates the DOCELL syntax presented in the paper.

x <-c(200,210,220,210,225,225,205,210,220)
sqrt( ((x[5] - x[1])^2) + ((x[5] - x[2])^2) + ((x[5] - x[3])^2) + 
      ((x[5] - x[4])^2) + ((x[5] - x[6])^2) + ((x[5] - x[7])^2) + 
      ((x[5] - x[8])^2) + ((x[5] - x[9])^2) ) 

Quick look at the population variance

1/length(x) * sum((x - mean(x))^2)
var(x) * (length(x) - 1) / length(x) # using R var function

Or, even better, mad(x) Median Absolute Deviation from Median (MAD), which would pick up the mode in the distribution and better represent peaks and valleys.

Not the answer you're looking for? Browse other questions tagged or ask your own question.