Reproducing terrain ruggedness index from the original paper example [closed]

Question

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:

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.

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

Any idea what am I missing?

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.