[GRASS-SVN] r71398 - grass-addons/grass7/raster/r.terrain.texture
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Aug 14 01:07:36 PDT 2017
Author: spawley
Date: 2017-08-14 01:07:36 -0700 (Mon, 14 Aug 2017)
New Revision: 71398
Added:
grass-addons/grass7/raster/r.terrain.texture/classification_example.png
grass-addons/grass7/raster/r.terrain.texture/convexity_example.png
grass-addons/grass7/raster/r.terrain.texture/flowchart.png
grass-addons/grass7/raster/r.terrain.texture/texture_example.png
Modified:
grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.html
grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.py
Log:
r.terrain.texture completion of surface form and classification
Added: grass-addons/grass7/raster/r.terrain.texture/classification_example.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.terrain.texture/classification_example.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/raster/r.terrain.texture/convexity_example.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.terrain.texture/convexity_example.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/raster/r.terrain.texture/flowchart.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.terrain.texture/flowchart.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.html
===================================================================
--- grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.html 2017-08-13 14:31:56 UTC (rev 71397)
+++ grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.html 2017-08-14 08:07:36 UTC (rev 71398)
@@ -1,42 +1,58 @@
<h2>DESCRIPTION</h2>
-<i>r.terrain.texture </i>calculates the Terrain Surface Texture of
-Iwahashi and Pike (2007). The calculation measures the frequency of peaks and
-pits in a DEM over a user specified window size. In practice, this measure
-represents terrain 'texture' or 'grain'.
+<em>r.terrain.texture</em> calculates the nested-means terrain classification of Iwahashi and Pike (2007). This classification procedure relies on three surface-form metrics consisting of: (a) terrain surface texture which is represented by the spatial frequency of peaks and pits; (b) surface-form convexity and concavity which are represented by the spatial frequency of convex/concave locations; and (c) slope-gradient which should be supplied by <em>r.slope.aspect</em> or <em>r.param.scale</em>. These metrics are combined using the mean of each variable as a dividing measure into a 8, 12 or 16 category classification of the topography.
<p>
-The calculation attempts to follow the description in Iwahashi and Pike
-(2007). Pits and peaks are firstly identified in the DEM by passing the
-DEM through a median filter with a 3x3 neighbourhood, and the pits and
-peaks are identified based on the difference between this smoothed DEM
-and the original terrain surface. By default, the algorithm uses a
-threshold of zero (i.e., any difference is identified as a pit or
-peak), however a user specified 'flatness' threshold can also be
-specified.
+The calculation follows the description in Iwahashi and Pike (2007). Terrain surface texture is calculated by smoothing the input <i>elevation</i> using a median filter over a neighborhood specified by the <i>filter_size</i> parameter (in pixels). Second, pits and peaks are extracted based on the difference between the smoothed DEM and the original terrain surface. By default the algorithm uses a threshold of 1 (> 1 m elevation difference) to identify pits and peaks. This threshold is set in the <i>flat_thres</i> parameter. The spatial frequency of pits and peaks is then calculated using a Gaussian resampling filter over a neighborhood size specified in the <i>counting_filter</i> parameter (default is 21 x 21 pixels, as per Iwahashi and Pike (2007).
+</p>
-The final terrain surface texture calculation involves <i>r.neighbors</i> to
-count the number of pits and peaks occurring within a user specified window
-size (default is 21 x 21, as per Iwahashi and Pike (2007)), and then this
-is converted into a percentage.
+<p>
+Surface-form convexity and concavity are calculated by first using a Laplacian filter approximating the second derivative of elevation to measure surface curvature. The Laplacian filter neighborhood size is set by the <i>filter_size</i> parameter (in pixels). This yields positive values in convex-upward areas, and negative values in concave areas, and zero on planar areas. Pixels are identified as either convex or concave when their values exceed the <i>curv_thres</i>. Similarly to terrain surface texture, the spatial frequency of these locations is then calculated to yield terrain surface convexity and concavity measures.
+</p>
+<p>
+Optionally, these surface-form metrics along with slope-gradient can be used to produce a nested-means classification of the topography. Each class is estimated based on whether the pixels values for each variable exceed the mean of that variable. The classification sequence follows:</p>
+
+<center>
+<img src="flowchart.png" alt="Terrain classification flowchart">
+</center>
+
+<p>A single iteration of this decision tree is completed for the 8-category classification. However for the 12 category classification, classes 1-4 remain but pixels that otherwise relate to classes 5-8 are passed to a second iteration of the decision tree, except that the mean of the gentler half of the area is used as the decision threshold, to produce 8 additional classes. Similarly for the 16 category classification, pixels that otherwise relate to classes 8-12 are passed onto a third iteration of the decision tree, except that the mean of the gentler quarter of the area is used as the decision threshold. This iterative subdivision of terrain properties acts to progressively partition the terrain into more gentle terrain features.</p>
+
<h2>NOTES</h2>
-<i>r.terrain.texture</i> performs best when using larger window sizes, either
-21 x 21 or 11 x 11 represent reasonable parameters.
+In the original algorithm description, SRTM data was smoothed using a fixed 3 x 3 size median filter and the spatial frequency of extracted features were measured over a 21 x 21 sized counting window. However, a larger smoothing filter size (~15 x 15) is often required to extract meaningful terrain features from higher resolution topographic data such as derived from LiDAR, and therefore both <i>filter_size</i> and <i>counting_size</i> parameters were exposed in the GRASS implementation. Further, if a large filter size is used then the counting window size should be increased accordingly.
-<h2>TODO</h2>
-The same approach was also used by Iwahashi and Pike (2007) to calculate
-terrain surface convexity, but using convexities and concavities rather
-than pits and peaks. The ability to calculate this will be implemented in
-the future. Also, these two morphometric variables can serve as a basis
-to classify the landscape into generic geomorphic categories.
-
<h2>EXAMPLE</h2>
+<p>Here we are going to use the GRASS GIS sample North Carolina data set as a basis to perform a terrain classification. First we set the computational region to the elev_state_500m dem, and then generate shaded relief (for visualization) and slope-gradient maps:</p>
+
<div class="code"><pre>
-r.terrain.texture.py elevation=DEM thres=0 pitdensity=Texture
+g.region raster=elev_state_500m at PERMANENT
+r.relief input=elev_state_500m at PERMANENT output=Hillshade_state_500m altitude=45 azimuth=315
+r.slope.aspect elevation=elev_state_500m at PERMANENT slope=slope_state_500m
</pre></div>
+<p>Then we produce the terrain classification:</p>
+
+<div class="code"><pre>
+r.terrain.texture elevation=elev_state_500m at PERMANENT slope=slope_state_500m at PERMANENT texture=texture_state_500m convexity=convexity_state_500m concavity=concavity_state_500m features=classification_state_500m
+</pre></div>
+
+<p>Terrain surface texture:</p>
+<center>
+<img src="texture_example.png" alt="Terrain surface texture result">
+</center>
+
+<p>Surface-form convexity:</p>
+<center>
+<img src="convexity_example.png" alt="Terrain convexity result">
+</center>
+
+<p>8-category terrain classification:</p>
+<center>
+<img src="classification_example.png" alt="8-category nested-means classification result">
+</center>
+
<h2>REFERENCES</h2>
Iwahashi, J., and Pike, R.J. 2007. Automated classifications of topography
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-13 14:31:56 UTC (rev 71397)
+++ grass-addons/grass7/raster/r.terrain.texture/r.terrain.texture.py 2017-08-14 08:07:36 UTC (rev 71398)
@@ -1,11 +1,11 @@
#! /usr/bin/env python
##############################################################################
#
-# MODULE: Pit and peak density
+# MODULE: Unsupervised nested-means algorithm for terrain classification
#
# AUTHOR(S): Steven Pawley
#
-# PURPOSE: Calculates the density of pits and peaks in a DEM
+# PURPOSE: Divides topography based on terrain texture and curvature.
# 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.
@@ -16,7 +16,7 @@
##############################################################################
#%module
-#% description: Pit and peak density
+#% description: Unsupervised nested-means algorithm for terrain classification
#% keyword: raster
#% keyword: terrain
#% keyword: classification
@@ -24,35 +24,88 @@
#%option G_OPT_R_INPUT
#% description: Input elevation raster:
-#% key: input
-#% required : yes
+#% key: elevation
+#% required: yes
#%end
+#%option G_OPT_R_INPUT
+#% description: Input slope raster:
+#% key: slope
+#% required: no
+#%end
+
#%option
-#% key: thres
+#% key: flat_thres
#% type: double
#% description: Height threshold for pit and peak detection:
#% answer: 1
-#% required: yes
+#% required: no
#%end
#%option
-#% key: size
+#% key: curv_thres
#% type: double
+#% description: Curvature threshold for convexity and concavity detection:
+#% answer: 0
+#% required: no
+#%end
+
+#%option
+#% key: filter_size
+#% type: integer
+#% description: Size of smoothing filter window:
+#% answer: 3
+#% guisection: Optional
+#%end
+
+#%option
+#% key: counting_size
+#% type: integer
#% description: Size of counting window:
#% answer: 21
#% guisection: Optional
#%end
+
+#%option
+#% key: classes
+#% type: integer
+#% description: Number of classes in nested terrain classification:
+#% options: 8,12,16
+#% answer: 8
+#% guisection: Optional
+#%end
+
#%option G_OPT_R_OUTPUT
-#% description: Output Texture Image:
-#% key: output
-#% required : yes
+#% description: Output terrain texture:
+#% key: texture
+#% required: yes
#%end
+#%option G_OPT_R_OUTPUT
+#% description: Output terrain convexity:
+#% key: convexity
+#% required: yes
+#%end
+
+#%option G_OPT_R_OUTPUT
+#% description: Output terrain concavity:
+#% key: concavity
+#% required: yes
+#%end
+
+#%option G_OPT_R_OUTPUT
+#% description: Output terrain classification:
+#% key: features
+#% required : no
+#%end
+
+import os
import sys
import random
import string
+import math
+import numpy as np
from subprocess import PIPE
import atexit
import grass.script as gs
@@ -64,59 +117,456 @@
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 parse_tiles(tiles):
+ # convert r.tileset output into list of dicts
+ tiles = [i.split(';') for i in tiles]
+ tiles = [[parse_key_val(i) for i in t] for t in tiles]
+
+ tiles_reg = []
+ for index, tile in enumerate(tiles):
+ tiles_reg.append({})
+ for param in tile:
+ tiles_reg[index].update(param)
+
+ return(tiles_reg)
+
+
+def laplacian_matrix(w):
+ s = """TITLE Laplacian filter
+ MATRIX {w}
+ """.format(w=w)
+
+ x = np.zeros((w,w))
+ x[:] = -1
+ x[w/2, w/2] = (np.square(w))-1
+ x_str=np.array2string(a=x)
+ x_str = x_str.replace(' [', '')
+ x_str = x_str.replace('[', '')
+ x_str = x_str.replace(']', '')
+ x_str = x_str.replace('.', '')
+ s += x_str
+ s += """
+ DIVISOR 1
+ TYPE P"""
+ return s
+
+
+def categories(nclasses):
+ if nclasses == 8:
+ s = """1|steep, high convexity, fine-textured
+ 2|steep, high convexity, coarse-textured
+ 3|steep, low convexity, fine-textured
+ 4|steep, low convexity, coarse-textured
+ 5|gentle, high convexity, fine-textured
+ 6|gentle, high convexity, coarse-textured
+ 7|gentle, low convexity, fine-textured
+ 8|gentle, low convexity, coarse-textured"""
+ elif nclasses == 12:
+ s = """1|steep, high convexity, fine-textured
+ 2|steep, high convexity, coarse-textured
+ 3|steep, low convexity, fine-textured
+ 4|steep, low convexity, coarse-textured
+ 5|moderately steep, high convexity, fine-textured
+ 6|moderately steep, high convexity, coarse-textured
+ 7|moderately steep, low convexity, fine-textured
+ 8|moderately steep, low convexity, coarse-textured
+ 9|gentle, high convexity, fine-textured
+ 10|gentle, high convexity, coarse-textured
+ 11|gentle, low convexity, fine-textured
+ 12|gentle, low convexity, coarse-textured"""
+ elif nclasses == 16:
+ s = """1|steep, high convexity, fine-textured
+ 2|steep, high convexity, coarse-textured
+ 3|steep, low convexity, fine-textured
+ 4|steep, low convexity, coarse-textured
+ 5|moderately steep, high convexity, fine-textured
+ 6|moderately steep, high convexity, coarse-textured
+ 7|moderately steep, low convexity, fine-textured
+ 8|moderately steep, low convexity, coarse-textured
+ 9|slightly steep, high convexity, fine-textured
+ 10|slightly steep, high convexity, coarse-textured
+ 11|slightly steep, low convexity, fine-textured
+ 12|slightly steep, low convexity, coarse-textured
+ 13|gentle, high convexity, fine-textured
+ 14|gentle, high convexity, coarse-textured
+ 15|gentle, low convexity, fine-textured
+ 16|gentle, low convexity, coarse-textured"""
+ return s
+
+def colors(nclasses):
+ if nclasses == 8:
+ s = """
+ 1 139:92:20
+ 2 255:0:0
+ 3 255:165:0
+ 4 255:255:0
+ 5 0:0:255
+ 6 30:144:255
+ 7 0:128:0
+ 8 0:255:3
+ nv 255:255:255
+ default 255:255:255
+ """
+ elif nclasses == 12:
+ s = """
+ 1 165:42:42
+ 2 255:0:0
+ 3 255:165:0
+ 4 255:255:0
+ 5 0:128:0
+ 6 11:249:11
+ 7 144:238:144
+ 8 227:255:227
+ 9 0:0:255
+ 10 30:144:255
+ 11 0:255:231
+ 12 173:216:230
+ nv 255:255:255
+ default 255:255:255
+ """
+ elif nclasses == 16:
+ s = """
+ 1 165:42:42
+ 2 255:0:0
+ 3 255:165:0
+ 4 255:255:0
+ 5 0:128:0
+ 6 11:249:11
+ 7 144:238:144
+ 8 227:255:227
+ 9 128:0:128
+ 10 244:7:212
+ 11 234:109:161
+ 12 231:190:225
+ 13 0:0:255
+ 14 30:144:255
+ 15 0:255:231
+ 16 173:216:230
+ nv 255:255:255
+ default 255:255:255
+ """
+ return s
+
+
def main():
- elevation = options['input']
- thres = float(options['thres'])
- window = int(options['size'])
- pitdensity = options['output']
+ elevation = options['elevation']
+ slope = options['slope']
+ flat_thres = float(options['flat_thres'])
+ curv_thres = float(options['curv_thres'])
+ filter_size = int(options['filter_size'])
+ counting_size = int(options['counting_size'])
+ nclasses = int(options['classes'])
+ texture = options['texture']
+ convexity = options['convexity']
+ concavity = options['concavity']
+ features = options['features']
- # current region settings
+ # remove mapset from output name in case of overwriting existing map
+ texture = texture.split('@')[0]
+ convexity = convexity.split('@')[0]
+ concavity = concavity.split('@')[0]
+ features = features.split('@')[0]
+
+ # error checking
+ if flat_thres < 0:
+ gs.fatal('Parameter thres cannot be negative')
+
+ if filter_size % 2 == 0 or counting_size % 2 == 0:
+ gs.fatal('Filter or counting windows require an odd-numbered window size')
+
+ if filter_size >= counting_size:
+ gs.fatal('Filter size needs to be smaller than the counting window size')
+
+ # store 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
- 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)
+ # check for existing mask
+ mask_test = gs.list_grouped(
+ type='rast', pattern='MASK')[gs.gisenv()['MAPSET']]
+ if mask_test:
+ # backup old mask
+ original_mask = temp_map('tmp_original_mask')
+ g.copy(raster=['MASK', original_mask])
+
+ # create mask of input elevation
+ elevation_mask = temp_map('tmp_elevation_mask')
+ tmp = gs.tempfile()
+ rc = open('%s' % (tmp), 'wt')
+ rc.write('* = 0')
+ rc.close()
+
+ r.reclass(input=elevation, output=elevation_mask, rules=tmp)
+ new_mask = temp_map('tmp_new_mask')
+ r.mapcalc(expression='{a}={b}'.format(a=new_mask, b=elevation_mask))
+
+ # Terrain Surface Texture -------------------------------------------------
+ # smooth the dem
+ gs.message("Calculating terrain surface texture...")
+ gs.message("1. Smoothing input DEM with a {n}x{n} median filter...".format(n=filter_size))
+ filtered_dem = temp_map('tmp_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
- pitpeaks = 'tmp_pitpeaks' + ''.join(
- [random.choice(string.ascii_letters + string.digits) for n in range(4)])
- TMP_RAST.append(pitpeaks)
+ size = filter_size, output = filtered_dem, flags='c',
+ quiet=True)
- 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),
+ # extract the pits and peaks based on the threshold
+ pitpeaks = temp_map('tmp_pitpeaks')
+ gs.message("2. Extracting pits and peaks with difference > thres...")
+ r.mapcalc(expression='{x} = if ( abs({dem}-{median})>{thres}, 1, 0)'.format(
+ x=pitpeaks, dem=elevation, thres=flat_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
+ gs.message("3. Using resampling filter to create terrain texture...")
+ window_radius = (counting_size-1)/2
y_radius = float(current_reg['ewres'])*window_radius
x_radius = float(current_reg['nsres'])*window_radius
+
+ resample = temp_map('tmp_density')
+ r.resamp_filter(input=pitpeaks, output=resample, filter=['bartlett','gauss'],
+ radius=[x_radius,y_radius], quiet=True)
+
+ # convert to percentage
+ if mask_test:
+ r.mask(raster=new_mask, overwrite=True, quiet=True)
+ else:
+ r.mask(raster=elevation, overwrite=True, quiet=True)
+
+ gs.message("4. Converting to percentage...")
+ r.mapcalc(expression='{x} = float({y} * 100)'.format(x=texture, y=resample),
+ quiet=True)
+ if mask_test:
+ r.mask(raster=original_mask, overwrite=True, quiet=True)
+ else:
+ r.mask(flags='r', quiet=True)
+
+ # set colors
+ r.colors(map=texture, color='haxby', quiet=True)
+
+ # Terrain convexity/concavity ---------------------------------------------
+ # surface curvature using lacplacian filter
+ gs.message("Calculating terrain convexity and concavity...")
+ gs.message("1. Calculating terrain curvature using laplacian filter...")
- r.resamp_filter(input=pitpeaks, output=pitdensity, filter=['bartlett','gauss'],
- radius=[x_radius,y_radius])
+ # grow the map to remove border effects
+ dem_grown = temp_map('tmp_elevation_grown')
+ g.region(n=float(current_reg['n']) + (float(current_reg['nsres']) * filter_size),
+ s=float(current_reg['s']) - (float(current_reg['nsres']) * filter_size),
+ w=float(current_reg['w']) - (float(current_reg['ewres']) * filter_size),
+ e=float(current_reg['e']) + (float(current_reg['ewres']) * filter_size))
+ r.grow(input=elevation, output=dem_grown, radius=filter_size, quiet=True)
- # return to original region
+ # write matrix filter to tempfile
+ tmp = gs.tempfile()
+ lmatrix = open('%s' % (tmp), 'wt')
+ lmatrix.write(laplacian_matrix(filter_size))
+ lmatrix.close()
+
+ laplacian = temp_map('tmp_laplacian')
+ r.mfilter(input=dem_grown, output=laplacian, filter=tmp, quiet=True)
+
+ # extract convex and concave pixels
+ gs.message("2. Extracting convexities and concavities...")
+ convexities = temp_map('tmp_convexities')
+ concavities = temp_map('tmp_concavities')
+
+ r.mapcalc(
+ expression='{x} = if({laplacian}>{thres}, 1, 0)'\
+ .format(x=convexities, laplacian=laplacian, thres=curv_thres),
+ quiet=True)
+ r.mapcalc(
+ expression='{x} = if({laplacian}<-{thres}, 1, 0)'\
+ .format(x=concavities, laplacian=laplacian, thres=curv_thres),
+ quiet=True)
+
+ # calculate density of convexities and concavities
+ gs.message("3. Using resampling filter to create surface convexity/concavity...")
+ resample_convex = temp_map('tmp_convex')
+ resample_concav = temp_map('tmp_concav')
+ r.resamp_filter(input=convexities, output=resample_convex,
+ filter=['bartlett','gauss'], radius=[x_radius,y_radius],
+ quiet=True)
+ r.resamp_filter(input=concavities, output=resample_concav,
+ filter=['bartlett','gauss'], radius=[x_radius,y_radius],
+ quiet=True)
+
+ # convert to percentages
g.region(**current_reg)
+ if mask_test:
+ r.mask(raster=new_mask, overwrite=True, quiet=True)
+ else:
+ r.mask(raster=elevation, overwrite=True, quiet=True)
+ gs.message("4. Converting to percentages...")
+ r.mapcalc(expression='{x} = float({y} * 100)'.format(x=convexity, y=resample_convex),
+ quiet=True)
+ r.mapcalc(expression='{x} = float({y} * 100)'.format(x=concavity, y=resample_concav),
+ quiet=True)
+ if mask_test:
+ r.mask(raster=original_mask, overwrite=True, quiet=True)
+ else:
+ r.mask(flags='r', quiet=True)
# set colors
- r.colors(map=pitdensity, color='haxby', quiet=True)
-
+ r.colors_stddev(map=convexity, quiet=True)
+ r.colors_stddev(map=concavity, quiet=True)
+
+ # Terrain classification Flowchart-----------------------------------------
+ if features != '':
+ gs.message("Performing terrain surface classification...")
+ # output result as final classification if nclasses == 8
+ # else pass temporary result to next threshold
+ if nclasses >= 8:
+ if nclasses == 8:
+ cla1_8=features
+ else:
+ cla1_8 = temp_map('tmp_classes1_8')
+
+ # gather univariate statistics for terrain maps
+ slope_mean = parse_key_val(r.univar(
+ map=slope, flags='g', stdout_=PIPE).outputs.stdout)['mean']
+ convex_mean = parse_key_val(r.univar(
+ map=convexity, flags='g', stdout_=PIPE).outputs.stdout)['mean']
+ texture_mean = parse_key_val(r.univar(
+ map=texture, flags='g', stdout_=PIPE).outputs.stdout)['mean']
+
+ # calculate first threshold
+ r.mapcalc(
+ expression='{x} = if({s}>{smean}, if({c}>{cmean}, if({t}<{tmean}, 1, 2), if({t}<{tmean}, 3,4)), if({c}>{cmean}, if({t}<{tmean}, 5, 6), if({t}<{tmean}, 7, 8)))'\
+ .format(x=cla1_8,
+ s=slope, smean=slope_mean,
+ t=texture, tmean=texture_mean,
+ c=convexity, cmean=convex_mean))
+
+ if nclasses >= 12:
+ cla9_12 = temp_map('tmp_classes9_12')
+
+ # output result as final classification if nclasses == 12
+ # else pass temporary result to next threshold
+ if nclasses == 12:
+ cla1_12 = features
+ else:
+ cla1_12 = temp_map('tmp_classes1_12')
+
+ # get projection and create tiles
+ proj = g.proj(flags='jf', stdout_=PIPE).outputs.stdout.strip(os.linesep)
+ tiles = r.tileset(
+ sourceproj=proj,
+ maxcols=int(math.ceil(float(current_reg['cols'])+1 / 2)),
+ maxrows=int(math.ceil(float(current_reg['rows'])+2)),
+ flags='g', quiet=True, stdout_=PIPE).outputs.stdout.split(os.linesep)[:-1]
+ tiles_reg = parse_tiles(tiles)
+
+ # get image means in each tile
+ slope_stats, convex_stats, texture_stats = [],[],[]
+ for tile in tiles_reg:
+ g.region(**tile)
+ slope_stats.append(parse_key_val(r.univar(
+ map=slope, flags='g', stdout_=PIPE).outputs.stdout)['mean'])
+ convex_stats.append(parse_key_val(r.univar(
+ map=convexity, flags='g', stdout_=PIPE).outputs.stdout)['mean'])
+ texture_stats.append(parse_key_val(r.univar(
+ map=texture, flags='g', stdout_=PIPE).outputs.stdout)['mean'])
+
+ # calculate second threshold classes
+ g.region(**current_reg)
+ r.mapcalc(
+ expression='{x} = if({s}>{smean}, if({c}>{cmean}, if({t}<{tmean}, 5, 6), if({t}<{tmean}, 7,8)), if({c}>{cmean}, if({t}<{tmean}, 9, 10), if({t}<{tmean}, 11, 12)))'\
+ .format(x=cla9_12,
+ s=slope, smean=min(slope_stats),
+ t=texture, tmean=min(texture_stats),
+ c=convexity, cmean=min(convex_stats)))
+
+ r.mapcalc(
+ expression='{x} = if({a}>5, {b}, {a})'.format(x=cla1_12, a=cla1_8, b=cla9_12))
+
+ if nclasses == 16:
+ cla9_16 = temp_map('tmp_classes9_16')
+
+ # get projection and create tiles
+ tiles = r.tileset(
+ sourceproj=proj,
+ maxcols=int(math.ceil(float(current_reg['cols'])+1 / 2)),
+ maxrows=int(math.ceil(float(current_reg['rows'])+1 / 2)),
+ flags='g', quiet=True, stdout_=PIPE).outputs.stdout.split(os.linesep)[:-1]
+ tiles_reg = parse_tiles(tiles)
+
+ # get image means in each tile
+ slope_mean, convex_mean, texture_mean = [],[],[]
+ for tile in tiles_reg:
+ g.region(**tile)
+ slope_mean.append(parse_key_val(r.univar(
+ map=slope, flags='g', stdout_=PIPE).outputs.stdout)['mean'])
+ convex_mean.append(parse_key_val(r.univar(
+ map=convexity, flags='g', stdout_=PIPE).outputs.stdout)['mean'])
+ texture_mean.append(parse_key_val(r.univar(
+ map=texture, flags='g', stdout_=PIPE).outputs.stdout)['mean'])
+
+ # third threshold
+ g.region(**current_reg)
+ r.mapcalc(
+ expression='{x} = if({s}>{smean}, if({c}>{cmean}, if({t}<{tmean}, 9, 10), if({t}<{tmean}, 11, 12)), if({c}>{cmean}, if({t}<{tmean}, 13, 14), if({t}<{tmean}, 15, 16)))'\
+ .format(x=cla9_16,
+ s=slope, smean=min(slope_mean),
+ t=texture, tmean=min(texture_mean),
+ c=convexity, cmean=min(convex_mean)))
+
+ r.mapcalc(
+ expression='{x} = if({a}>8, {b}, {a})'.format(x=features, a=cla1_12, b=cla9_16))
+
+
+ # Write metadata ----------------------------------------------------------
+ history = 'r.terrain.texture '
+ for key,val in options.iteritems():
+ history += key + '=' + str(val) + ' '
+
+ r.support(map=texture,
+ title=texture,
+ description='generated by r.terrain.texture',
+ history=history)
+ r.support(map=convexity,
+ title=convexity,
+ description='generated by r.terrain.texture',
+ history=history)
+ r.support(map=concavity,
+ title=concavity,
+ description='generated by r.terrain.texture',
+ history=history)
+
+ if features != '':
+ r.support(map=features,
+ title=features,
+ description='generated by r.terrain.texture',
+ history=history)
+
+ # write color and category rules to tempfiles
+ tmp_col = gs.tempfile()
+ col = open('%s' % (tmp_col), 'wt')
+ col.write(colors(nclasses))
+ col.close()
+
+ tmp_cat = gs.tempfile()
+ cat = open('%s' % (tmp_cat), 'wt')
+ cat.write(categories(nclasses))
+ cat.close()
+
+ r.category(map=features, rules=tmp_cat, separator='pipe')
+ r.colors(map=features, rules=tmp_col, quiet=True)
+
return 0
if __name__ == "__main__":
Added: grass-addons/grass7/raster/r.terrain.texture/texture_example.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.terrain.texture/texture_example.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
More information about the grass-commit
mailing list