Python "Ruggedness" Analysis Code

[RachelRodriguez]

I need a code which essentially... "The terrain 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 terrain 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."

All I could figure out was....

ssdiff = ( ( sqr ( "in_elev"(0,0) - "in_elev"(-1,-1) ) ) + ( sqr ( "in_elev"(0,0) - "in_elev"(0,-1) ) ) + ( sqr ( "in_elev"(0,0) - "in_elev"(1,-1) ) ) + ( sqr ( "in_elev"(0,0) - "in_elev"(1,0) ) ) + ( sqr ( "in_elev"(0,0) - "in_elev"(1,1) ) ) + ( sqr ( "in_elev"(0,0) - "in_elev"(0,1) ) ) + ( sqr ( "in_elev"(0,0) - "in_elev"(-1,1) ) ) + ( sqr ( "in_elev"(0,0) - "in_elev"(-1,0) ) ) )

But how do I even try and translate it onto Python and into ArcMap. I have never tried to breach the ArcMap/Python boundary.

Please HELP!

Thanks, Rachel

[ChrisSnyder]

The sample you code you posted uses something called "DOCELL" notation (sometimes called neighborhood notation) to refer to neighboring cells using relative coordinates (for example, 1,-1). Unfortunately, ESRI still has not ported that functionality in ArcGIS. However, it is currently still supported in GRID via Workstation ArcInfo, so you could actually just run your code as is via Workstation. There are some pretty simple ways to do that via Python even (that is call grid.exe through Python).

This probably isn't a very good "beginner" Python project, but there are some ways of doing this in ArcGIS/Python:

1. Load the data into a Python dictionary or numpy array and write your own function (this would require some Python knowledge for sure.

2. A simpler idea is to make 8 separate grids that are each "shifted" in one of the eight directions. So for example, produce a grid shifted one pixel to the north, another shifted 1 pixel to the NE, another shifted to the E, etc. Use the "Shift" tool in ArcToolbox to do this. Name each shifted grid according to the direction of the shift - so you will have grids n, ne, e, se, s, sw, w, nw - also name your original non-shifted grid c (for center). Then use the Combine tool to overlay all the grids together. Hopefully your grids aren't too big; otherwise this might not work due to too many unique values in the combo grid... You may need to adjust  the max raster attribute limit in the AdvancedArcMapSettings.exe registry editor to accomplish this. Anyway, PRESTO! In the output combo grid there will then be fields for each of the direction grids (the values representing the elevations of each cardinal direction relative to "c"). Add a new field "TRI" and just calc it according to the TRI equation.

[curtvprice]

The sample you code you posted uses something called "DOCELL" notation (sometimes called neighborhood notation) to refer to neighboring cells using relative coordinates (for example, 1,-1). Unfortunately, ESRI still has not ported that functionality in ArcGIS.

Just to set the record straight, the "[,]" neighborhood notation is supported (and is much more efficient) outside of DOCELL blocks in ArcInfo GRID.

There are no plans I've heard to support neighborhood notation in ArcGIS; members of the raster team have also pointed me to raster to numpy array to accomplish this task.

(Are ya listening, your rasteroid guys and gals in Redlands? IMHO going to numpy is a bit kludgy to access what I feel is basic raster modeling functionality. Sure would be nice if you supported neighborhood notation as a method on the raster object, say: "myRaster.nbr(1,1)" or something like that.)

I just thought of another approach outside of shifting grids around - which is to calculate focal sums using custom kernel neighborhoods to access those neighborhood cells. This has the potential of performing better than using the shift tool, which most certainly would involve copying the whole grid to access neighborhood cell.

If you have access to ArcInfo GRID I'd look for an old person to help you do it there. 🙂

[JeffreyEvans]

The Riley et. al., (1999) index is merely unscaled variance. You can calculate an approximation using map algebra. Conveniently, this is scalable without modifying the neighborhood notation.

Variance in a 3x3 window
roughness = SQR( FOCALSTDV( dem, RECTANGLE, 3, 3 ) )

[MarkSappington]

Hi Rachel,

You could also try my script for a vector-based ruggedness measure with seems to decouple slope from ruggedness better than TRI. The script for ArcGIS 9.x and 10.x can be found here: http://resources.arcgis.com/gallery/file/geoprocessing/details?entryID=F65FF927-1422-2418-A02A-EE725...

Hope it helps.

Mark Sappington
GIS Specialist
Lake Mead National Recreation Area
mark_sappington@nps.gov

[TaylerHamilton]

I am trying to figure out the same problem as you, Rachel.

I looked into Jeffrey's post a little further because it seemed like a simple way to solve the TRI, but I have some questions about it:

If the TRI is: SQRT ((c-1)^2 + (c-2)^2 + (c-3)^2 + (c-4)^2 + (c-5)^2 + (c-6)^2 + (c-7)^2 + (c-8)^2)), where c is the focal centre within the 3x3 neighbourhood.

And standard deviation of a dataset is the square root of the variance, and variance is the expected value of the squared difference ( variance(x) = E((x-u)^2) ).

Then wouldn't TRI = (STDV[dem], rectangle, 3, 3)?

The difference from Jeffrey's is the lack of square rooting the equation, which I think is not necessary since the standard deviation already calculates the square root of the sum of the squared differences within the 3x3 grid.

If this is in fact the case, then the TRI can be calculated by using the Focal Statistics tool and solving with the standard deviation function within a 3x3 rectangle.

I am not sure if this I am correct, so any feedback is welcome.

[JeffreyEvans]

In grid algebra SQR is the square not the square root. The square of the standard deviation is the variance. Since the original TRI is only the sum of the squared deviations it is unscaled (not corrected by n-1) variance. So calculating a focal standard deviation and then taking the square of it returns an approximation of the TRI.

[TaylerHamilton]

My mistake, I meant square, not square root. The difference between the way I saw it and your equation was the lack of squaring after applying the standard deviation in mine.

I still don't understand how squaring it will scale the TRI. What does the approximation of the TRI mean?

Thank you for your help!!

[JeffreyEvans]

OK, so the original Riley et. al. (1999) paper does not take the square root, it is the sum of the squared deviations. Taking the square root of the TRI was an unpublished modification added later. If you take the sum of the squared deviations and divide by n-1 you get the sample variance. If you look at the code that the original poster pulled from the TRI AML you will notice that it does not take the square root. If you square of a focal standard deviation you get something quite similar to the original TRI (same resulting distributional shape with different values, thus an approximation). If you want the modification then just use the focal standard deviation.

The advantage of using the focal standard deviation (squared or not) is that, by changing the size of the focal window, the TRI becomes scalable (spatially). It is a bit erroneous to assume that a organism or process is going to respond to a scale represented only by a 3x3 window. This is something that should be tested across spatial scales.

Clear as mud, right?

[TaylerHamilton]

Awesome! Thanks for your help, that makes things more understandable.