r.tri added option to use parallelized r.mapcalc.tiled during the cal…

…culation

grass7/raster/r.tri/r.tri.html

@@ -1,17 +1,55 @@
 <h2>DESCRIPTION</h2>
 
-<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.
+<p><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.</p>
 
 <h2>NOTES</h2>
-<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.
 
-Like in many GRASS GIS algorithms, cell padding is not performed automatically 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. To minimize this effect the flag -p can be set to grow the DEM by the chosen radius prior to the TRI calculation.
+<p><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 that can potentially
+    better highlight finer-scale terrain features.
+</p>
+
+<p>
+    Similar to many other GRASS GIS algorithms, cell padding is not performed automatically which will leave null values
+    at the boundaries of the output raster relative to the size of the input raster. To minimize this effect the DEM
+    needs to be buffered/grown prior to using r.tri.
+</p>
+
+<p>
+    To reduce computational times for large raster datasets, setting <em>n_jobs</em> &gt 1 will use a parallelized and
+    tiled calculations that is spread across multiple processing cores. This option requires the
+    GRASS addon <a href="r.mapcalc.tiled.html">r.mapcalc.tiled</a> to be installed.
+</p>
 
 <h2>EXAMPLE</h2>
-r.tri input=srtm size=1 output=tri
+
+<div class="code">
+ <pre>
+     d.rast map=elev_lid792_1m@PERMANENT
+     g.region raster=elev_lid792_1m@PERMANENT -a
+     r.tri input=elev_lid792_1m@PERMANENT size=9 output=tri
+ </pre>
+</div>
+
+<center>
+    <img src="tri.png" alt="Terrain Ruggedness Index">
+</center>
 
 <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.
+
+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>SEE ALSO</h2>
+
+<em>
+<a href="r.mapcalc.html">r.mapcalc</a>
+<a href="r.mapcalc.tiled.html">r.mapcalc.tiled</a>
+</em>
 
 <h2>AUTHOR</h2>
 Steven Pawley

grass7/raster/r.tri/r.tri.py

@@ -37,37 +37,47 @@
 #% guisection: Required
 #%end
 
+#%option
+#% key: n_jobs
+#% type: integer
+#% description: Number of processes cores for computation
+#% required: no
+#% answer: 1
+#% end
+
 #%flag
 #% key: c
 #% description: Use circular neighborhood
 #%end
 
-from grass.pygrass.modules.shortcuts import general as g
-from grass.pygrass.modules.shortcuts import raster as r
-from grass.pygrass.gis.region import Region
-from grass.pygrass.raster import RasterRow
 import sys
 import atexit
 import random
 import string
+import multiprocessing as mp
+import math
 import grass.script as gs
+import grass.script.core as gcore
+from grass.pygrass.modules.shortcuts import raster as grast
+from grass.pygrass.gis.region import Region
 
 TMP_RAST = []
 
+
 def cleanup():
     gs.message("Deleting intermediate files...")
 
     for f in TMP_RAST:
-        gs.run_command(
-            "g.remove", type="raster", name=f, flags="f", quiet=True)
+        gs.run_command("g.remove", type="raster", name=f, flags="f", quiet=True)
 
 
 def temp_map(name):
-    tmp = name + ''.join(
-        [random.choice(string.ascii_letters + string.digits) for n in range(4)])
+    tmp = name + "".join(
+        [random.choice(string.ascii_letters + string.digits) for n in range(4)]
+    )
     TMP_RAST.append(tmp)
 
-    return (tmp)
+    return tmp
 
 
 def focal_expr(radius, window_square=False):
@@ -94,57 +104,90 @@ def focal_expr(radius, window_square=False):
     # generate a list of spatial neighbourhood offsets for the chosen radius
     # ignoring the centre cell
     if window_square:
-        
-        for i in range(-radius, radius+1):
-            for j in range(-radius, radius+1):
-                if (i,j) != (0,0):
+
+        for i in range(-radius, radius + 1):
+            for j in range(-radius, radius + 1):
+                if (i, j) != (0, 0):
                     offsets.append((i, j))
-    
+
     else:
 
-        for i in range(-radius, radius+1):
-            for j in range(-radius, radius+1):
+        for i in range(-radius, radius + 1):
+            for j in range(-radius, radius + 1):
                 row = i + radius
                 col = j + radius
 
-                if pow(row - radius, 2) + pow(col - radius, 2) <= \
-                    pow(radius, 2) and (i, j) != (0,0):
+                if pow(row - radius, 2) + pow(col - radius, 2) <= pow(radius, 2) and (
+                    i,
+                    j,
+                ) != (0, 0):
                     offsets.append((j, i))
 
     return offsets
 
 
 def main():
-    dem = options['input']
-    tri = options['output']
-    size = int(options['size'])
-    circular = flags['c']
+    dem = options["input"]
+    tri = options["output"]
+    size = int(options["size"])
+    n_jobs = int(options["n_jobs"])
+    circular = flags["c"]
 
-    radius = int((size-1)/2)
+    radius = int((size - 1) / 2)
+
+    if n_jobs != 1:
+        if not gcore.find_program("r.mapcalc.tiled", "--help"):
+            msg = "Parallelized computation requires the extension 'r.mapcalc.tiled' to be installed."
+            msg = " ".join([msg, "Install it using 'g.extension r.mapcalc.tiled'"])
+            gs.fatal(msg)
+
+        if n_jobs == 0:
+            gs.fatal(
+                "Number of processing cores for parallel computation must not equal 0"
+            )
+
+        if n_jobs < 0:
+            system_cores = mp.cpu_count()
+            n_jobs = system_cores + n_jobs + 1
 
     # calculate TRI based on map calc statements
     gs.message("Calculating the Topographic Ruggedness Index...")
 
     # generate a list of spatial neighbourhood offsets for the chosen radius
     # ignoring the center cell
-    offsets = focal_expr(radius = radius, window_square=not circular)
+    offsets = focal_expr(radius=radius, window_square=not circular)
     terms = []
     for d in offsets:
-        valid = ','.join(map(str, d))        
-        invalid = ','.join([str(-d[0]), str(-d[1])])
+        valid = ",".join(map(str, d))
+        invalid = ",".join([str(-d[0]), str(-d[1])])
         terms.append(
             "if(isnull({dem}[{d}]), if(isnull({dem}[{e}]), (median({dem})-{dem})^2, ({dem}[{e}]-{dem}^2)), ({dem}[{d}]-{dem})^2)".format(
-            dem=dem, d=valid, e=invalid))
+                dem=dem, d=valid, e=invalid
+            )
+        )
 
     # define the calculation expression
-    ncells = len(offsets)+1
-    expr = "$tri = sqrt((%s" % " + ".join(terms) + ") / $ncells)"
+    ncells = len(offsets) + 1
+    terms = " + ".join(terms)
+
+    # expr = "$tri = sqrt((%s" % " + ".join(terms) + ") / $ncells)"
+
+    expr = "{tri} = sqrt(({terms}) / {ncells})".format(
+        tri=tri, ncells=ncells, terms=terms
+    )
 
     # perform the r.mapcalc calculation with the moving window
-    gs.mapcalc(expr, tri=tri, dem=dem, ncells=ncells)
+    if n_jobs == 1:
+        gs.mapcalc(expr, tri=tri, dem=dem, ncells=ncells)
+    else:
+        reg = Region()
+        width = math.ceil(reg.cols / n_jobs)
+        height = math.ceil(reg.rows / n_jobs)
+        grast.mapcalc_tiled(expr, width=width, height=height, processes=n_jobs, overlap=radius)
 
     return 0
 
+
 if __name__ == "__main__":
     options, flags = gs.parser()
     atexit.register(cleanup)