[GRASS-SVN] r67161 - grass-addons/grass7/raster/r.tri
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Dec 16 00:33:09 PST 2015
Author: spawley
Date: 2015-12-16 00:33:09 -0800 (Wed, 16 Dec 2015)
New Revision: 67161
Modified:
grass-addons/grass7/raster/r.tri/r.tri.html
grass-addons/grass7/raster/r.tri/r.tri.py
Log:
Updated TRI calculation to use PyGRASS interface and extend calculation over any neighbourhood size
Modified: grass-addons/grass7/raster/r.tri/r.tri.html
===================================================================
--- grass-addons/grass7/raster/r.tri/r.tri.html 2015-12-15 22:09:34 UTC (rev 67160)
+++ grass-addons/grass7/raster/r.tri/r.tri.html 2015-12-16 08:33:09 UTC (rev 67161)
@@ -1,19 +1,19 @@
<h2>DESCRIPTION</h2>
-<i>r.tri </i>calculates the Terrain Ruggedness Index (TRI) of Riley et al. (1999). The index represents the sum change in elevation between a grid cell and its neighbours, over a user-specified moving window size. Unlike the original calculation in Riley et al., (1999) which used only a 3x3 neighbourhood, the TRI result is averaged over the number of neighbouring cells to allow the calculation to be extended over different window sizes. The calculation employed here produced exactly the same results as when using the terrain ruggedness focal function in the R raster package, which provides an alternative way to calculate the index using the R interface.
+<i>r.tri </i>calculates the Terrain Ruggedness Index (TRI) of Riley et al. (1999). The index represents the mean change in elevation between a grid cell and its neighbours, over a user-specified moving window size. The original calculation in Riley et al., (1999) used only a 3x3 neighbourhood and represented the sum of the absolute deviations between the center pixel and its immediate 8 neighbours. In r.tri, this calculation is modified so that the calculation can be extended over any scale by taking the mean of the absolute deviations. The results scale in a similar way to other packages that allow a multi-scale calculation (i.e., SAGA-GIS).
<h2>NOTES</h2>
-<i>r.tri</i> uses <i>r.mapcalc</i> to calculate the sum of differences between the centre pixel and its neighbours. <i>r.mapcalc</i> statements are limited in character length, so only a range of neighbourhood sizes are allowed based on the methodology employed. If other window sizes are required, the R package <i>raster</i> can be used via <i>rgrass7</i>.
+<i>r.tri</i> produces fairly similar results to the average deviation of elevation values, apart from the center pixel is used in place of the mean. In practice, this produces a slightly 'less smoothed' result, which in some cases can better highlight smaller-scale terrain features.
<h2>TODO</h2>
-Null values are not evaluated at present, so there will be an edge effect and the resulting TRI will not be calculated at the image edges so there will be missing pixels along the margins relative to the size of the input raster. There is potential to calculate TRI at the image edges by using the <i>is.null()</i> r.mapcalc function.
+Like in many GRASS GIS algorithms, null values are not evaluated at present and there will be an edge effect. The resulting TRI will not be calculated at the image edges so there will be missing pixels along the margins relative to the size of the input raster. There is potential to calculate TRI at the image edges by using a boolean operator to test for null values.
<h2>EXAMPLE</h2>
-r.tri dem=srtm wsize=3 output=tri
+r.tri dem=srtm radius=1 output=tri
<h2>REFERENCES</h2>
Riley, S. J., S. D. DeGloria and R. Elliot (1999). A terrain ruggedness index that quantifies topographic heterogeneity, Intermountain Journal of Sciences, vol. 5, No. 1-4, 1999.
<h2>AUTHOR</h2>
Steven Pawley
-<br><i>Last changed: Saturday 12 December 2015</i>
\ No newline at end of file
+<br><i>Last changed: Tuesday 15 December 2015</i>
\ No newline at end of file
Modified: grass-addons/grass7/raster/r.tri/r.tri.py
===================================================================
--- grass-addons/grass7/raster/r.tri/r.tri.py 2015-12-15 22:09:34 UTC (rev 67160)
+++ grass-addons/grass7/raster/r.tri/r.tri.py 2015-12-16 08:33:09 UTC (rev 67161)
@@ -2,21 +2,19 @@
#
##############################################################################
#
-# MODULE: Topographic Ruggedness Index
+# MODULE: Terrain Ruggedness Index
#
# AUTHOR(S): Steven Pawley
#
-# PURPOSE: Calculates the Topographic Ruggedness Index (TRI) of
+# PURPOSE: Calculates the Terrain Ruggedness Index (TRI) of
# Riley et al. (1999)
#
# COPYRIGHT: (C) 2015 Steven Pawley, and by the GRASS Development Team
#
-# DATE: May 5 2015
-#
##############################################################################
#%module
-#% description: Topographic Ruggedness Index
+#% description: Terrain Ruggedness Index
#%end
#%option G_OPT_R_INPUT
@@ -26,7 +24,7 @@
#%end
#%option G_OPT_R_OUTPUT
-#% description: Output Topographic Ruggedness Index (TRI)
+#% description: Output Terrain Ruggedness Index (TRI)
#% key: tri
#% required : yes
#%end
@@ -34,10 +32,9 @@
#%option
#% key: wsize
#% type: integer
-#% description: Size of neighbourhood (valid 3,5,7,9,11)
-#% answer: 3
+#% description: Radius of neighbourhood in cells
+#% answer: 1
#% guisection: Required
-#% options: 3,5,7,9,11
#%end
import sys
@@ -48,22 +45,28 @@
dem = options['dem']
tri = options['tri']
wsize = int(options['wsize'])
- neighcells = wsize^2 -1
+ neighcells = ((wsize*2+1)**2)-1
# calculate TRI based on map calc statements
- # NB could create script to calculate matrix but mapcalc is limited by statement length
grass.message("Calculating the Topographic Ruggedness Index:")
- if wsize == 3:
- grass.mapcalc('{x} = (abs({a}[0,0]-{a}[-1,-1])+abs({a}[0,0]-{a}[0,-1])+abs({a}[0,0]-{a}[1,-1])+abs({a}[0,0]-{a}[-1,0])+abs({a}[0,0]-{a}[0,0])+abs({a}[0,0]-{a}[1,0])+abs({a}[0,0]-{a}[-1,1])+abs({a}[0,0]-{a}[0,1])+abs({a}[0,0]-{a}[1,1]))/1'.format(x=tri, a=dem, b=neighcells))
- elif wsize == 5:
- grass.mapcalc('{x} = (abs({a}[0,0]-{a}[-2,-2])+abs({a}[0,0]-{a}[-1,-2])+abs({a}[0,0]-{a}[0,-2])+abs({a}[0,0]-{a}[1,-2])+abs({a}[0,0]-{a}[2,-2])+abs({a}[0,0]-{a}[-2,-1])+abs({a}[0,0]-{a}[-1,-1])+abs({a}[0,0]-{a}[0,-1])+abs({a}[0,0]-{a}[1,-1])+abs({a}[0,0]-{a}[2,-1])+abs({a}[0,0]-{a}[-2,0])+abs({a}[0,0]-{a}[-1,0])+abs({a}[0,0]-{a}[0,0])+abs({a}[0,0]-{a}[1,0])+abs({a}[0,0]-{a}[2,0])+abs({a}[0,0]-{a}[-2,1])+abs({a}[0,0]-{a}[-1,1])+abs({a}[0,0]-{a}[0,1])+abs({a}[0,0]-{a}[1,1])+abs({a}[0,0]-{a}[2,1])+abs({a}[0,0]-{a}[-2,2])+abs({a}[0,0]-{a}[-1,2])+abs({a}[0,0]-{a}[0,2])+abs({a}[0,0]-{a}[1,2])+abs({a}[0,0]-{a}[2,2]))/{b}'.format(x=tri, a=dem, b=neighcells))
- elif wsize == 7:
- grass.mapcalc('{x} = (abs({a}[0,0]-{a}[-3,-3])+abs({a}[0,0]-{a}[-2,-3])+abs({a}[0,0]-{a}[-1,-3])+abs({a}[0,0]-{a}[0,-3])+abs({a}[0,0]-{a}[1,-3])+abs({a}[0,0]-{a}[2,-3])+abs({a}[0,0]-{a}[3,-3])+abs({a}[0,0]-{a}[-3,-2])+abs({a}[0,0]-{a}[-2,-2])+abs({a}[0,0]-{a}[-1,-2])+abs({a}[0,0]-{a}[0,-2])+abs({a}[0,0]-{a}[1,-2])+abs({a}[0,0]-{a}[2,-2])+abs({a}[0,0]-{a}[3,-2])+abs({a}[0,0]-{a}[-3,-1])+abs({a}[0,0]-{a}[-2,-1])+abs({a}[0,0]-{a}[-1,-1])+abs({a}[0,0]-{a}[0,-1])+abs({a}[0,0]-{a}[1,-1])+abs({a}[0,0]-{a}[2,-1])+abs({a}[0,0]-{a}[3,-1])+abs({a}[0,0]-{a}[-3,0])+abs({a}[0,0]-{a}[-2,0])+abs({a}[0,0]-{a}[-1,0])+abs({a}[0,0]-{a}[0,0])+abs({a}[0,0]-{a}[1,0])+abs({a}[0,0]-{a}[2,0])+abs({a}[0,0]-{a}[3,0])+abs({a}[0,0]-{a}[-3,1])+abs({a}[0,0]-{a}[-2,1])+abs({a}[0,0]-{a}[-1,1])+abs({a}[0,0]-{a}[0,1])+abs({a}[0,0]-{a}[1,1])+abs({a}[0,0]-{a}[2,1])+abs({a}[0,0]-{a}[3,1])+abs({a}[0,0]-{a}[-3,2])+abs({a}[0,0]-{a}[-2,2])+abs({a}[0,0]-{a}[-1,2])+abs({a}[0,0]-{a}[0,2])+abs({a}[0,0]-{a}[1,2])+abs({a}[
0,0]-{a}[2,2])+abs({a}[0,0]-{a}[3,2])+abs({a}[0,0]-{a}[-3,3])+abs({a}[0,0]-{a}[-2,3])+abs({a}[0,0]-{a}[-1,3])+abs({a}[0,0]-{a}[0,3])+abs({a}[0,0]-{a}[1,3])+abs({a}[0,0]-{a}[2,3])+abs({a}[0,0]-{a}[3,3]))/{b}'.format(x=tri, a=dem, b=neighcells))
- elif wsize == 9:
- grass.mapcalc('{x} = (abs({a}[0,0]-{a}[-4,-4])+abs({a}[0,0]-{a}[-3,-4])+abs({a}[0,0]-{a}[-2,-4])+abs({a}[0,0]-{a}[-1,-4])+abs({a}[0,0]-{a}[0,-4])+abs({a}[0,0]-{a}[1,-4])+abs({a}[0,0]-{a}[2,-4])+abs({a}[0,0]-{a}[3,-4])+abs({a}[0,0]-{a}[4,-4])+abs({a}[0,0]-{a}[-4,-3])+abs({a}[0,0]-{a}[-3,-3])+abs({a}[0,0]-{a}[-2,-3])+abs({a}[0,0]-{a}[-1,-3])+abs({a}[0,0]-{a}[0,-3])+abs({a}[0,0]-{a}[1,-3])+abs({a}[0,0]-{a}[2,-3])+abs({a}[0,0]-{a}[3,-3])+abs({a}[0,0]-{a}[4,-3])+abs({a}[0,0]-{a}[-4,-2])+abs({a}[0,0]-{a}[-3,-2])+abs({a}[0,0]-{a}[-2,-2])+abs({a}[0,0]-{a}[-1,-2])+abs({a}[0,0]-{a}[0,-2])+abs({a}[0,0]-{a}[1,-2])+abs({a}[0,0]-{a}[2,-2])+abs({a}[0,0]-{a}[3,-2])+abs({a}[0,0]-{a}[4,-2])+abs({a}[0,0]-{a}[-4,-1])+abs({a}[0,0]-{a}[-3,-1])+abs({a}[0,0]-{a}[-2,-1])+abs({a}[0,0]-{a}[-1,-1])+abs({a}[0,0]-{a}[0,-1])+abs({a}[0,0]-{a}[1,-1])+abs({a}[0,0]-{a}[2,-1])+abs({a}[0,0]-{a}[3,-1])+abs({a}[0,0]-{a}[4,-1])+abs({a}[0,0]-{a}[-4,0])+abs({a}[0,0]-{a}[-3,0])+abs({a}[0,0]-{a}[-2,0])+abs({a}[0,0]-{a
}[-1,0])+abs({a}[0,0]-{a}[0,0])+abs({a}[0,0]-{a}[1,0])+abs({a}[0,0]-{a}[2,0])+abs({a}[0,0]-{a}[3,0])+abs({a}[0,0]-{a}[4,0])+abs({a}[0,0]-{a}[-4,1])+abs({a}[0,0]-{a}[-3,1])+abs({a}[0,0]-{a}[-2,1])+abs({a}[0,0]-{a}[-1,1])+abs({a}[0,0]-{a}[0,1])+abs({a}[0,0]-{a}[1,1])+abs({a}[0,0]-{a}[2,1])+abs({a}[0,0]-{a}[3,1])+abs({a}[0,0]-{a}[4,1])+abs({a}[0,0]-{a}[-4,2])+abs({a}[0,0]-{a}[-3,2])+abs({a}[0,0]-{a}[-2,2])+abs({a}[0,0]-{a}[-1,2])+abs({a}[0,0]-{a}[0,2])+abs({a}[0,0]-{a}[1,2])+abs({a}[0,0]-{a}[2,2])+abs({a}[0,0]-{a}[3,2])+abs({a}[0,0]-{a}[4,2])+abs({a}[0,0]-{a}[-4,3])+abs({a}[0,0]-{a}[-3,3])+abs({a}[0,0]-{a}[-2,3])+abs({a}[0,0]-{a}[-1,3])+abs({a}[0,0]-{a}[0,3])+abs({a}[0,0]-{a}[1,3])+abs({a}[0,0]-{a}[2,3])+abs({a}[0,0]-{a}[3,3])+abs({a}[0,0]-{a}[4,3])+abs({a}[0,0]-{a}[-4,4])+abs({a}[0,0]-{a}[-3,4])+abs({a}[0,0]-{a}[-2,4])+abs({a}[0,0]-{a}[-1,4])+abs({a}[0,0]-{a}[0,4])+abs({a}[0,0]-{a}[1,4])+abs({a}[0,0]-{a}[2,4])+abs({a}[0,0]-{a}[3,4])+abs({a}[0,0]-{a}[4,4]))/{b}'.format(x=tri, a=dem, b=
neighcells))
- elif wsize == 11:
- grass.mapcalc('{x} = (abs({a}[0,0]-{a}[-5,-5])+abs({a}[0,0]-{a}[-4,-5])+abs({a}[0,0]-{a}[-3,-5])+abs({a}[0,0]-{a}[-2,-5])+abs({a}[0,0]-{a}[-1,-5])+abs({a}[0,0]-{a}[0,-5])+abs({a}[0,0]-{a}[1,-5])+abs({a}[0,0]-{a}[2,-5])+abs({a}[0,0]-{a}[3,-5])+abs({a}[0,0]-{a}[4,-5])+abs({a}[0,0]-{a}[5,-5])+abs({a}[0,0]-{a}[-5,-4])+abs({a}[0,0]-{a}[-4,-4])+abs({a}[0,0]-{a}[-3,-4])+abs({a}[0,0]-{a}[-2,-4])+abs({a}[0,0]-{a}[-1,-4])+abs({a}[0,0]-{a}[0,-4])+abs({a}[0,0]-{a}[1,-4])+abs({a}[0,0]-{a}[2,-4])+abs({a}[0,0]-{a}[3,-4])+abs({a}[0,0]-{a}[4,-4])+abs({a}[0,0]-{a}[5,-4])+abs({a}[0,0]-{a}[-5,-3])+abs({a}[0,0]-{a}[-4,-3])+abs({a}[0,0]-{a}[-3,-3])+abs({a}[0,0]-{a}[-2,-3])+abs({a}[0,0]-{a}[-1,-3])+abs({a}[0,0]-{a}[0,-3])+abs({a}[0,0]-{a}[1,-3])+abs({a}[0,0]-{a}[2,-3])+abs({a}[0,0]-{a}[3,-3])+abs({a}[0,0]-{a}[4,-3])+abs({a}[0,0]-{a}[5,-3])+abs({a}[0,0]-{a}[-5,-2])+abs({a}[0,0]-{a}[-4,-2])+abs({a}[0,0]-{a}[-3,-2])+abs({a}[0,0]-{a}[-2,-2])+abs({a}[0,0]-{a}[-1,-2])+abs({a}[0,0]-{a}[0,-2])+abs({a}[0,0
]-{a}[1,-2])+abs({a}[0,0]-{a}[2,-2])+abs({a}[0,0]-{a}[3,-2])+abs({a}[0,0]-{a}[4,-2])+abs({a}[0,0]-{a}[5,-2])+abs({a}[0,0]-{a}[-5,-1])+abs({a}[0,0]-{a}[-4,-1])+abs({a}[0,0]-{a}[-3,-1])+abs({a}[0,0]-{a}[-2,-1])+abs({a}[0,0]-{a}[-1,-1])+abs({a}[0,0]-{a}[0,-1])+abs({a}[0,0]-{a}[1,-1])+abs({a}[0,0]-{a}[2,-1])+abs({a}[0,0]-{a}[3,-1])+abs({a}[0,0]-{a}[4,-1])+abs({a}[0,0]-{a}[5,-1])+abs({a}[0,0]-{a}[-5,0])+abs({a}[0,0]-{a}[-4,0])+abs({a}[0,0]-{a}[-3,0])+abs({a}[0,0]-{a}[-2,0])+abs({a}[0,0]-{a}[-1,0])+abs({a}[0,0]-{a}[0,0])+abs({a}[0,0]-{a}[1,0])+abs({a}[0,0]-{a}[2,0])+abs({a}[0,0]-{a}[3,0])+abs({a}[0,0]-{a}[4,0])+abs({a}[0,0]-{a}[5,0])+abs({a}[0,0]-{a}[-5,1])+abs({a}[0,0]-{a}[-4,1])+abs({a}[0,0]-{a}[-3,1])+abs({a}[0,0]-{a}[-2,1])+abs({a}[0,0]-{a}[-1,1])+abs({a}[0,0]-{a}[0,1])+abs({a}[0,0]-{a}[1,1])+abs({a}[0,0]-{a}[2,1])+abs({a}[0,0]-{a}[3,1])+abs({a}[0,0]-{a}[4,1])+abs({a}[0,0]-{a}[5,1])+abs({a}[0,0]-{a}[-5,2])+abs({a}[0,0]-{a}[-4,2])+abs({a}[0,0]-{a}[-3,2])+abs({a}[0,0]-{a}[-2,2])+abs({a}
[0,0]-{a}[-1,2])+abs({a}[0,0]-{a}[0,2])+abs({a}[0,0]-{a}[1,2])+abs({a}[0,0]-{a}[2,2])+abs({a}[0,0]-{a}[3,2])+abs({a}[0,0]-{a}[4,2])+abs({a}[0,0]-{a}[5,2])+abs({a}[0,0]-{a}[-5,3])+abs({a}[0,0]-{a}[-4,3])+abs({a}[0,0]-{a}[-3,3])+abs({a}[0,0]-{a}[-2,3])+abs({a}[0,0]-{a}[-1,3])+abs({a}[0,0]-{a}[0,3])+abs({a}[0,0]-{a}[1,3])+abs({a}[0,0]-{a}[2,3])+abs({a}[0,0]-{a}[3,3])+abs({a}[0,0]-{a}[4,3])+abs({a}[0,0]-{a}[5,3])+abs({a}[0,0]-{a}[-5,4])+abs({a}[0,0]-{a}[-4,4])+abs({a}[0,0]-{a}[-3,4])+abs({a}[0,0]-{a}[-2,4])+abs({a}[0,0]-{a}[-1,4])+abs({a}[0,0]-{a}[0,4])+abs({a}[0,0]-{a}[1,4])+abs({a}[0,0]-{a}[2,4])+abs({a}[0,0]-{a}[3,4])+abs({a}[0,0]-{a}[4,4])+abs({a}[0,0]-{a}[5,4])+abs({a}[0,0]-{a}[-5,5])+abs({a}[0,0]-{a}[-4,5])+abs({a}[0,0]-{a}[-3,5])+abs({a}[0,0]-{a}[-2,5])+abs({a}[0,0]-{a}[-1,5])+abs({a}[0,0]-{a}[0,5])+abs({a}[0,0]-{a}[1,5])+abs({a}[0,0]-{a}[2,5])+abs({a}[0,0]-{a}[3,5])+abs({a}[0,0]-{a}[4,5])+abs({a}[0,0]-{a}[5,5]))/{b}'.format(x=tri, a=dem, b=neighcells))
+ # generate a list of spatial neighbourhood indexs for the chosen radius
+ # ignoring the center cell
+ offsets = []
+ for j in range(-wsize, wsize+1):
+ for i in range(-wsize, wsize+1):
+ if (j,i) != (0,0):
+ offsets.append((j,i))
+
+ # define the calculation term
+ terms = ["abs($dem - $dem[%d,%d])" % d for d in offsets]
+
+ # define the calculation expression
+ expr = "$tri = (%s" % " + ".join(terms) + ") / $neighcells"
+
+ # perform the r.mapcalc calculation with the moving window
+ grass.mapcalc(expr, tri=tri, dem=dem, neighcells=neighcells)
+
return 0
if __name__ == "__main__":
More information about the grass-commit
mailing list