[GRASS-SVN] r71359 - grass-addons/grass7/raster/r.terrain.texture
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Aug 9 10:40:37 PDT 2017
Author: spawley
Date: 2017-08-09 10:40:37 -0700 (Wed, 09 Aug 2017)
New Revision: 71359
Modified:
grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.py
Log:
enhancement to r.terrain.texture using resampling filter to calculate pit and peak density
Modified: grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.py
===================================================================
--- grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.py 2017-08-09 17:15:22 UTC (rev 71358)
+++ grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.py 2017-08-09 17:40:37 UTC (rev 71359)
@@ -6,12 +6,12 @@
# AUTHOR(S): Steven Pawley
#
# PURPOSE: Calculates the density of pits and peaks in a DEM
-# Based on the methodology of Iwahashi & Pike (2007)
-# Automated classifications of topography from DEMs by an unsupervised
-# nested-means algorithm and a three-part geometric signature.
-# Geomorphology. 86, 409-440
+# Based on the methodology of Iwahashi & Pike (2007)
+# Automated classifications of topography from DEMs by an unsupervised
+# nested-means algorithm and a three-part geometric signature.
+# Geomorphology. 86, 409-440
#
-# COPYRIGHT: (C) 2015 Steven Pawley, and by the GRASS Development Team
+# COPYRIGHT: (C) 2017 Steven Pawley and by the GRASS Development Team
#
##############################################################################
@@ -21,7 +21,7 @@
#%option G_OPT_R_INPUT
#% description: Input elevation raster:
-#% key: elevation
+#% key: input
#% required : yes
#%end
@@ -34,7 +34,7 @@
#%end
#%option
-#% key: window
+#% key: size
#% type: double
#% description: Size of counting window:
#% answer: 21
@@ -43,59 +43,80 @@
#%option G_OPT_R_OUTPUT
#% description: Output Texture Image:
-#% key: pitdensity
+#% key: output
#% required : yes
#%end
import sys
import random
import string
-import grass.script as grass
+from subprocess import PIPE
+import atexit
+import grass.script as gs
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.script.utils import parse_key_val
+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 main():
- elevation = options['elevation']
- thres = options['thres']
- window = options['window']
- pitdensity = options['pitdensity']
+ elevation = options['input']
+ thres = float(options['thres'])
+ window = int(options['size'])
+ pitdensity = options['output']
- # Internal grid calculations - to be removed at end of process
- median = 'tmp_' + ''.join([random.choice(string.ascii_letters + string.digits) for n in range(8)])
- pitpeaks = 'tmp_' + ''.join([random.choice(string.ascii_letters + string.digits) for n in range(8)])
- pitcount = 'tmp_' + ''.join([random.choice(string.ascii_letters + string.digits) for n in range(8)])
- pitdensity_tmp = 'tmp_' + ''.join([random.choice(string.ascii_letters + string.digits) for n in range(8)])
+ # current region settings
+ current_reg = parse_key_val(
+ g.region(flags='pg', stdout_=PIPE).outputs.stdout)
+ del current_reg['projection']
+ del current_reg['zone']
+ del current_reg['cells']
# Smooth the DEM
- grass.message("Smooth DEM with a 3x3 median filter:")
- grass.run_command("r.neighbors", input = elevation, method = "median", size = 3, output = median)
-
+ gs.message("Smoothing input DEM with a 3x3 median filter...")
+ filtered_dem = 'tmp_filtred_dem' + ''.join(
+ [random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(filtered_dem)
+
+ gs.run_command("r.neighbors", input = elevation, method = "median",
+ size = 3, output = filtered_dem, flags='c', quiet=True)
+
# Extract the pits and peaks based on the threshold
- grass.message("Extract pits and peaks:")
- grass.mapcalc('{x} = if ( abs({dem}-{median})>{thres}, 1, null())'.format(x=pitpeaks, dem=elevation, thres=thres, median=median))
+ pitpeaks = 'tmp_pitpeaks' + ''.join(
+ [random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(pitpeaks)
- # Count the number of pits and peaks
- grass.message("Count the number of pits and peaks in a moving window:")
- grass.run_command("r.neighbors", input = pitpeaks, method = "count", size = window, output = pitcount, flags = 'c')
+ gs.message("Extracting pits and peaks with difference > thres...")
+ gs.mapcalc(
+ '{x} = if ( abs({dem}-{median})>{thres}, 1, 0)'.format(
+ x=pitpeaks, dem=elevation, thres=thres, median=filtered_dem),
+ quiet=True)
+
+ # calculate density of pits and peaks
+ gs.message("Using resampling filter to create terrain texture...")
+ window_radius = (window-1)/2
+ y_radius = float(current_reg['ewres'])*window_radius
+ x_radius = float(current_reg['nsres'])*window_radius
+
+ r.resamp_filter(input=pitpeaks, output=pitdensity, filter=['bartlett','gauss'],
+ radius=[x_radius,y_radius])
+
+ # return to original region
+ g.region(**current_reg)
- # Convert the pit and peak counts into a percentage
- grass.message("Convert the pit and peak counts into a percentage:")
- size = int(window)*int(window)
- grass.mapcalc('{x} = 100*{pitcount}/float({size})'.format(x=pitdensity_tmp, pitcount=pitcount, size = size))
+ # set colors
+ r.colors(map=pitdensity, color='haxby', quiet=True)
- # Mask the raster to remove the background
- grass.message("Mask the terrain surface texture with the input DEM:")
- grass.run_command("r.mask", raster = elevation, maskcats = "*", layer = "1")
- grass.mapcalc('{x} = {pitdensity_tmp}'.format(x=pitdensity, pitdensity_tmp=pitdensity_tmp))
- grass.run_command("r.mask", raster = elevation, maskcats = "*", layer = "1", flags = 'r')
-
- # Clean-up
- grass.message("Deleting intermediate files:")
- grass.run_command("g.remove", type="raster", name=median, flags="f", quiet=True)
- grass.run_command("g.remove", type="raster", name=pitpeaks, flags="f", quiet=True)
- grass.run_command("g.remove", type="raster", name=pitcount, flags="f", quiet=True)
- grass.run_command("g.remove", type="raster", name=pitdensity_tmp, flags="f", quiet=True)
-
return 0
if __name__ == "__main__":
- options, flags = grass.parser()
+ options, flags = gs.parser()
+ atexit.register(cleanup)
sys.exit(main())
More information about the grass-commit
mailing list