[GRASS-SVN] r71538 - grass-addons/grass7/raster/r.tri

[svn_grass at osgeo.org]

Author: spawley
Date: 2017-10-07 12:46:51 -0700 (Sat, 07 Oct 2017)
New Revision: 71538

Modified:
   grass-addons/grass7/raster/r.tri/r.tri.html
   grass-addons/grass7/raster/r.tri/r.tri.py
Log:
r.tri enhancement - added option to perform cell padding to minimize edge effect

Modified: grass-addons/grass7/raster/r.tri/r.tri.html
===================================================================
--- grass-addons/grass7/raster/r.tri/r.tri.html	2017-10-07 19:31:47 UTC (rev 71537)
+++ grass-addons/grass7/raster/r.tri/r.tri.html	2017-10-07 19:46:51 UTC (rev 71538)
@@ -5,7 +5,7 @@
 <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 avoid this the input DEM can be grown by the chosen radius that is to be used for the TRI calculation.
+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.
 
 <h2>EXAMPLE</h2>
 r.tri input=srtm size=1 output=tri

Modified: grass-addons/grass7/raster/r.tri/r.tri.py
===================================================================
--- grass-addons/grass7/raster/r.tri/r.tri.py	2017-10-07 19:31:47 UTC (rev 71537)
+++ grass-addons/grass7/raster/r.tri/r.tri.py	2017-10-07 19:46:51 UTC (rev 71538)
@@ -37,19 +37,73 @@
 #% guisection: Required
 #%end
 
+#%flag
+#% key: p
+#% label: Cell padding
+#% description: Use cell padding to reduce edge effect
+#% guisection: Optional
+#%end
+
+
 import sys
+import atexit
+import random
+import string
+from subprocess import PIPE
 import grass.script as gs
+from grass.script.utils import parse_key_val
+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
 
+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)
+
+
+def temp_map(name):
+    tmp = name + ''.join(
+        [random.choice(string.ascii_letters + string.digits) for n in range(4)])
+    TMP_RAST.append(tmp)
+
+    return (tmp)
+
+
 def main():
     dem = options['input']
     tri = options['output']
     size = int(options['size'])
+    p = flags['p']
 
     radius = (size-1)/2
     neighcells = (size**2)-1
 
+    # store current region settings
+    current_reg = Region()
+
+    if p:
+        # check for existing mask
+        if RasterRow('MASK').exist():
+            gs.fatal('Cell padding option cannot be used with existing mask...please remove first')
+
+        # grow the map to remove border effects
+        gs.message('Growing DEM to avoid edge effects...')
+        dem_grown = temp_map('tmp_elevation_grown')
+        g.region(n=current_reg.north + (current_reg.nsres * (radius+2)),
+                 s=current_reg.south - (current_reg.nsres * (radius+2)),
+                 w=current_reg.west - (current_reg.ewres * (radius+2)),
+                 e=current_reg.east + (current_reg.ewres * (radius+2)))
+        r.grow(input=dem, output=dem_grown, radius=radius+2, quiet=True)
+        dem = dem_grown
+
     # calculate TRI based on map calc statements
-    gs.message("Calculating the Topographic Ruggedness Index:")
+    gs.message("Calculating the Topographic Ruggedness Index...")
 
     # generate a list of spatial neighbourhood indexs for the chosen radius
     # ignoring the center cell
@@ -66,10 +120,19 @@
     expr = "$tri = float(%s" % " + ".join(terms) + ") / $neighcells"
 
     # perform the r.mapcalc calculation with the moving window
-    gs.mapcalc(expr, tri=tri, dem=dem, neighcells=neighcells)
+    if p:
+        output_tmp = temp_map('tmp')
+        gs.mapcalc(expr, tri=output_tmp, dem=dem, neighcells=neighcells)
+        r.mask(raster=options['input'], quiet=True)
+        r.mapcalc('{x}={y}'.format(x=tri, y=output_tmp))
+        r.mask(flags='r', quiet=True)
+        Region.write(current_reg)
+    else:
+        gs.mapcalc(expr, tri=tri, dem=dem, neighcells=neighcells)
 
     return 0
 
 if __name__ == "__main__":
     options, flags = gs.parser()
+    atexit.register(cleanup)
     sys.exit(main())