[GRASS-SVN] r72654 - grass-addons/grass7/raster/r.valley.bottom
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Apr 29 15:17:12 PDT 2018
Author: spawley
Date: 2018-04-29 15:17:12 -0700 (Sun, 29 Apr 2018)
New Revision: 72654
Modified:
grass-addons/grass7/raster/r.valley.bottom/r.valley.bottom.html
grass-addons/grass7/raster/r.valley.bottom/r.valley.bottom.py
Log:
partial rewrite of r.valley.bottom to fix output and eliminate shrinking
Modified: grass-addons/grass7/raster/r.valley.bottom/r.valley.bottom.html
===================================================================
--- grass-addons/grass7/raster/r.valley.bottom/r.valley.bottom.html 2018-04-29 14:19:55 UTC (rev 72653)
+++ grass-addons/grass7/raster/r.valley.bottom/r.valley.bottom.html 2018-04-29 22:17:12 UTC (rev 72654)
@@ -1,41 +1,53 @@
<h2>DESCRIPTION</h2>
-<p>
-<em>r.valley.bottom</em> calculates a Multi-resolution Valley Bottom
-Flatness (MrVBF) index. The user must specify the input <b>elevation
-raster</b> map.
-</p>
+<em>r.valley.bottom</em> calculates the Multi-resolution Valley Bottom Flatness (MRVBF) index (Gallant and Dowling, 2003). The MRVBF index assesses the flatness and lowness of terrain over multiple scales and DEM resolutions in order to identify valley bottoms, which represent areas that are flat across multiple scales, and remain low relative to the surrounding relief at coarser scales. The algorithm uses a sigmoid/logistic transform to rescale terrain slope angles and elevation percentile into a 0 to 1 range, and then combines these results across multiple levels of DEM smoothing and coarser grid resolutions. Although the resulting index represents a continuous value, values < 0.5 do not generally represent valley bottoms, values from 0.5 to 1.5 represent the steepest resolvable valley bottoms, and flatter/larger valley bottoms are represented by values > 1.5.
-<p>
-Default moving window is cells in 8 directions. With the <b>-s</b> flag
-a square moving window is used in calculations.
-</p>
+<h2>NOTES</h2>
+The user must specify the input <b>elevation</b> raster map as a required input. The output is given by the <b>mrvbf</b> argument. Optionally, the complementary Multiresolution Index of Ridge Top Flatness can be calculated by specifying the <b>mrrtf</b> argument. In addition, there are several parameters than can be used to change the behaviour of the argument, although note that in this case the results and their interpretation will differ from what was envisaged in the original paper. However, in practice, this is often required especially for high-resolution DEMs. The arguments are:
+
+<p></p>
+<li><em>t_slope</em> represents the initial threshold (t) for slope angle (in percentage). This specifies the slope angle that corresponds to a (logit) rescaled flatness value of 0.5. This means that slope angles lower than t_slope will be considered as flat areas, and slope angles higher than t_slope will be represented as non-flat areas. t_slope should be set based on the resolution of the input elevation dataset, and the algorithm was designed using with a 25 m DEM having a t_slope value of 16. Otherwise the t_slope value should by halved for every resolution step (a step consisting of a 3 x coarsening of resolution) above a 25 m resolution. For example, a 75 m DEM (3 x 25 m, 1 step) should have a t_slope value of 8, and a 250 m DEM (~2 resolution steps) should have a t_slope value of 4.</li>
+
+<li><em>p_slope</em> represents the shape parameter (p) for the sigmoid transformation. It defines the slope of the sigmoid function, i.e. how quickly changes in slope angle scale to being flat vs. non-flat areas. High p_slope values will cause a slow, smooth transition from flat areas to steep areas. Low p_slope values will result in much more rapid transitions that highlight more local vs. regional relief.</li>
+
+<li><em>t_pctl_v</em> represents the threshold for transformation of elevation percentile to evaluate lowness. This represents the threshold value for elevation percentile by which values less than this value will represent low areas. Elevation percentile represents the ratio of pixels of lower elevation relative to the total number of pixels in a moving-window neighborhood. Similarly <em>t_pctl_r</em> is the equivalent 'upness' threshold for the MRRTF index.</li>
+
+<li><em>p_pctl</em> represents the shape parameter (p) for the transformation of the elevation percentile. It defines the slope of the sigmoid function and governs how quickly transitions occur from low areas to upland areas.</li>
+
+<li><em>t_vf</em> and <em>t_rf</em> represent the thresholds for identifying valley bottoms (or ridge tops). Larger values indicate increasing valley bottom characteristics, with values less than 0.5 considered not to be in valley bottoms.</li>
+
+<p>The calculation of elevation percentile by default is performed using a circular window. With the <b>-s</b> flag a square moving window is used in calculations.</p>
+
+<p>In practice, the user does not usually need to alter the threshold-related parameters other than t_slope. However, changing the shape parameters can be useful for to emphasize more local vs. more regional variations in relief. The number of generalization steps can also be adjusted by the <em>levels</em> argument. Usually a setting between 4 to 6 provides the best results, and a minimum of three steps is required.</p>
+
<h2>EXAMPLE</h2>
+<p>Here we are going to use the GRASS GIS sample North Carolina data set as a basis to calculate the MRVBF index.</p>
+
<div class="code">
<pre>
- # align region to DEM and habitat vector
- g.region -a raster=DEM align=DEM
-
+ # align region to DEM
+ g.region -a raster=el_D782_6m
+
# run <em>r.valley.bottom</em>
- r.r.valley.bottom elevation=DEM
+ r.valley.bottom elevation=el_D782_6m mrvbf=mrvbf_el_D782_6m t_slope=40 levels=6 p_slope=3 p_pctl=2
+
+ # set colors
+ r.colors map=mrvbf_el_D782_6m color=bcyr -i
+
+ # display over a shaded relief map
+ r.relief input=el_D782_6m output=hs_D782_6m altitude=45 azimuth=315 zscale=4 scale=1
+ r.shade shade=hs_D782_6m color=mrvbf_el_D782_6m output=mrvbf_shade
+ d.rast map=mrvbf_shade
</pre>
</div>
-<h2>NOTES</h2>
+<center>
+<img src="mrvbf.png" alt="Multiresolution Index of Valley Bottom Flatness">
+</center>
-experimental
-<h2>TODO</h2>
-
-<p>
-Design the script independent of input elevation raster resolution. At
-the moment the input elevation raster map must have a resolution of 25m
-x 25m. Additionally the search window of the elevation percentile may be
-adapted.
-</p>
-
<h2>SEE ALSO</h2>
<em>
@@ -51,7 +63,7 @@
<h2>AUTHOR</h2>
-Helmut Kudrnovsky
+Helmut Kudrnovsky & Steven Pawley
<p>
<i>Last changed: $Date$</i>
Modified: grass-addons/grass7/raster/r.valley.bottom/r.valley.bottom.py
===================================================================
--- grass-addons/grass7/raster/r.valley.bottom/r.valley.bottom.py 2018-04-29 14:19:55 UTC (rev 72653)
+++ grass-addons/grass7/raster/r.valley.bottom/r.valley.bottom.py 2018-04-29 22:17:12 UTC (rev 72654)
@@ -5,14 +5,15 @@
MODULE: r.valley.bottom
AUTHOR(S): Helmut Kudrnovsky <alectoria AT gmx at>
+ Steven Pawley <dr.stevenpawley at gmail.com>
PURPOSE: Calculates Multi-resolution Valley Bottom Flatness (MrVBF) index
- Index is inspired by
+ Index is inspired by
John C. Gallant and Trevor I. Dowling 2003.
A multiresolution index of valley bottom flatness for mapping depositional areas.
WATER RESOURCES RESEARCH, VOL. 39, NO. 12, 1347, doi:10.1029/2002WR001426, 2003
-COPYRIGHT: (C) 2014 by the GRASS Development Team
+COPYRIGHT: (C) 2018 by the GRASS Development Team
This program is free software under the GNU General Public
License (>=v2). Read the file COPYING that comes with GRASS
@@ -27,657 +28,594 @@
#%option G_OPT_R_ELEV
#% key: elevation
-#% description: Name of elevation raster map
+#% description: Name of elevation raster map
#% required: yes
#%end
-#%flag
-#% key: s
-#% description: Use square moving window instead of moving window in 8 directions
+#%option G_OPT_R_OUTPUT
+#% key: mrvbf
+#% description: Name of output MRVBF raster map
+#% required: yes
#%end
+#%option G_OPT_R_OUTPUT
+#% key: mrrtf
+#% description: Name of output MRRTF raster map
+#% required: no
+#%end
- ##################################################################################
- # TODO: Design the script independent of input elevation raster resolution. At
- # the moment the input elevation raster map must have a resolution of 25m
- # x 25m. Additionally the search window of the elevation percentile may be
- # adapted.
+#%option
+#% key: t_slope
+#% description: Initial Threshold for Slope
+#% required: no
+#% answer: 16
+#% end
-import sys
-import os
-import grass.script as grass
-import math
+#%option
+#% key: t_pctl_v
+#% description: Threshold (t) for transformation of Elevation Percentile (Lowness)
+#% required: no
+#% answer: 0.4
+#% end
-if not os.environ.has_key("GISBASE"):
- grass.message( "You must be in GRASS GIS to run this program." )
- sys.exit(1)
+#%option
+#% key: t_pctl_r
+#% description: Threshold (t) for transformation of Elevation Percentile (Upness)
+#% required: no
+#% answer: 0.3
+#% end
-def main():
- r_elevation = options['elevation'].split('@')[0] # elevation raster
- r_slope_degree = r_elevation+'_slope_degree' # slope (percent) step 1
- moving_window_square = flags['s']
-
- ##################################################################################
- # Region settings
- grass.run_command('g.region', flags = 'a',
- raster = r_elevation,
- align = r_elevation)
- current_region = grass.region()
- Yres_base = current_region['nsres']
- Xres_base = current_region['ewres']
- grass.message( "Finest scale calculations" )
- grass.message( "Original resolution: %s x %s" % (Xres_base, Yres_base))
- grass.run_command('g.region', flags = 'p',
- save = "base_region_MrVBF",
- overwrite = True)
+#%option
+#% key: t_vf
+#% description: Threshold (t) for transformation of Valley Bottom Flatness
+#% required: no
+#% answer: 0.3
+#% end
- # Region settings step 3: resolution 3 x base resolution
+#%option
+#% key: t_rf
+#% description: Threshold (t) for transformation of Ridge Top Flatness
+#% required: no
+#% answer: 0.35
+#% end
- Yres_step3 = Yres_base * 3
- Xres_step3 = Xres_base * 3
+#%option
+#% key: p_slope
+#% description: Shape Parameter (p) for Slope
+#% required: no
+#% answer: 4
+#% end
- # Region settings step 4: resolution 4 x Xres_step3
+#%option
+#% key: p_pctl
+#% description: Shape Parameter (p) for Elevation Percentile
+#% required: no
+#% answer: 3
+#% end
- Yres_step4 = Yres_step3 * 3
- Xres_step4 = Xres_step3 * 3
-
- ##################################################################################
- # Step 1 (S1)
- # Finest-Scale Step
- # Calculation of slope step 1
- # 100 * tan (slope in degree)
- grass.message( "----" )
- grass.message( "Step 1: Calculation of slope" )
- grass.message( "..." )
- grass.run_command('r.resamp.rst', input = r_elevation,
- slope = r_slope_degree,
- ew_res = Xres_base,
- ns_res = Yres_base,
- quiet = True)
+#%option
+#% key: levels
+#% description: Number of generalization levels
+#% required: no
+#% answer: 4
+#% end
- grass.message( "..." )
-
- grass.mapcalc("$outmap = 100 * tan($slope) ",
- outmap = "r_slope_step1",
- slope = r_slope_degree)
+#%flag
+#% key: s
+#% description: Use square moving window instead of circular moving window
+#%end
- grass.message( "..." )
-
- grass.run_command("g.remove", flags="f", type="raster", name= r_slope_degree,
- quiet = True)
- grass.message( "Step 1: Calculation of slope done." )
- grass.message( "----" )
- # Calculation of flatness 1 (0 to 1)
- # Flatness F1 = 1 / (1 + pow ((slope / 16), 4)
- grass.message( "Step 1: Calculation of flatness F1" )
- grass.message( "..." )
- grass.mapcalc("F1 = 1.0 / (1.0 + pow( ( $slope / 16.0 ) ,4 ) )",
- slope = "r_slope_step1")
+from grass.pygrass.gis.region import Region
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+import sys
+import os
+import grass.script as grass
+import math
+import atexit
+import random
+import string
- grass.message( "Step 1: Calculation of flatness F1 done." )
- grass.message( "----" )
+if not os.environ.has_key("GISBASE"):
+ grass.message("You must be in GRASS GIS to run this program.")
+ sys.exit(1)
- # Calculation of elevation percentile PCTL1 step 1
- # elevation percentile: ratio of the number of points of lower
- # elevation to the total number of points in the surrounding region
- # elevation of percentile step 1: radius of 3 cells
- grass.message( "Step 1: Calculation of elevation percentile PCTL1" )
- grass.message( "..." )
- grass.mapcalc("elevation_step1 = $r_elevation",
- r_elevation = r_elevation)
+current_region = Region()
+TMP_RAST = []
- if moving_window_square :
+def cleanup():
+ grass.message("Deleting intermediate files...")
- offsets1 = [d
- for j in xrange(1,3+1)
- for i in [j,-j]
- for d in [(i,0),(i,1),(i,2),(i,3),(i,-1),(i,-2),(i,-3)]]
+ for f in TMP_RAST:
+ grass.run_command(
+ "g.remove", type="raster", name=f, flags="f", quiet=True)
+
+ Region.write(current_region)
- grass.message( "..." )
-
- terms = ["(elevation_step1[%d,%d] < elevation_step1)" % d
- for d in offsets1]
-
- expr1 = "PCTL1 = (100.0 / 48.0) * (%s)" % " + ".join(terms)
-
- else :
-
- offsets1 = [d
- for j in xrange(1,3+1)
- for i in [j,-j]
- for d in [(i,0),(0,i),(i,i),(i,-i)]]
- grass.message( "..." )
-
- terms = ["(elevation_step1[%d,%d] < elevation_step1)" % d
- for d in offsets1]
-
- expr1 = "PCTL1 = (100.0 / 24.0) * (%s)" % " + ".join(terms)
-
- grass.mapcalc( expr1 )
-
- grass.message( "Step 1: Calculation of elevation percentile PCTL1 done" )
- grass.message( "----" )
+def cell_padding(input, output, radius=3):
+ """Mitigates edge effect by growing an input raster map by radius cells
+
+ Args
+ ----
+ input, output : str
+ Names of GRASS raster map for input, and padded output
+ radius : int
+ Radius in which to expand region and grow raster
- # Elevation percentile is transformed to a local lowness value using
- # equation (1) with t = 0.4 and p = 3, and then combined with
- # flatness F1 to produce the preliminary valley flatness index
- # (PVF1) for the first step:
- grass.message( "Step 1: Calculation of preliminary valley flatness index PVF1" )
- grass.message( "..." )
-
- grass.mapcalc("PVF1 = F1 * (1.0 / (1.0 + pow ( ( PCTL1 / 0.4 ), 3 ) ) )")
+ Returns
+ -------
+ input_grown : str
+ GRASS raster map which has been expanded by radius cells"""
+
+ region = Region()
- grass.message( "Step 1: Calculation of preliminary valley flatness index PVF1 done." )
- grass.message( "----" )
+ g.region(n=region.north + (region.nsres * radius),
+ s=region.south - (region.nsres * radius),
+ w=region.west - (region.ewres * radius),
+ e=region.east + (region.ewres * radius))
- # Calculation of the valley flatness step 1 VF1
- # larger values of VF1 indicate increasing valley bottom character
- # with values less than 0.5 considered not to be in valley bottoms
- grass.message( "Step 1: Calculation of valley flatness VF1" )
- grass.message( "..." )
-
- grass.mapcalc("VF1 = 1 - (1.0 / (1.0 + pow ( ( PVF1 / 0.3 ), 4 ) ) )")
-
- grass.message( "Step 1: Calculation of valley flatness VF1 done." )
- grass.message( "----" )
+ r.grow(input=input,
+ output=output,
+ radius=radius+1,
+ quiet=True)
- ##################################################################################
- # Step 2 (S2)
- # base scale resolution
- grass.message( "Step 2" )
- grass.message( "base resolution: %s x %s" % (Xres_base, Yres_base))
- grass.message( "----" )
-
- # Calculation of flatness 1 (0 to 1)
- # The second step commences the same way with the
- # original DEM at its base resolution, using a slope threshold ts,2 half of ts,1:
- # Flatness F2 = 1 / (1 + pow ((slope / 8), 4)
- grass.message( "Step 2: Calculation of flatness F2" )
- grass.message( "..." )
- grass.mapcalc("F2 = 1.0 / (1.0 + pow( ( $slope / 8.0 ) ,4 ) )",
- slope = "r_slope_step1")
- grass.message( "Step 2: Calculation of flatness F2 done" )
- grass.message( "----" )
+ return (region)
- # Calculation of elevation percentile PCTL2 step 2
- # elevation of percentile step 2: radius of 6 cells
- grass.message( "Step 2: Calculation of elevation percentile PCTL2" )
- grass.message( "..." )
- if moving_window_square :
+def focal_expr(radius, window_square=False):
+ """Returns array offsets relative to centre pixel (0,0) for a matrix of
+ size radius
- offsets2 = [d
- for j in xrange(1,6+1)
- for i in [j,-j]
- for d in [(i,0),(i,1),(i,2),(i,3),(i,4),(i,5),(i,6),(i,-1),(i,-2),(i,-3),(i,-4),(i,-5),(i,-6)]]
+ Args
+ ----
+ radius : int
+ Radius of the focal function
+ window_square : bool. Optional (default is False)
+ Boolean to use a circular or square window
- grass.message( "..." )
-
- terms2 = ["(elevation_step1[%d,%d] < elevation_step1)" % d
- for d in offsets2]
-
- expr2 = "PCTL2 = (100.0 / 168.0) * (%s)" % " + ".join(terms2)
-
- else :
-
- offsets2 = [d
- for j in xrange(1,6+1)
- for i in [j,-j]
- for d in [(i,0),(0,i),(i,i),(i,-i)]]
+ Returns
+ -------
+ offsets : list
+ List of pixel positions (row, col) relative to the center pixel
+ ( 1, -1) ( 1, 0) ( 1, 1)
+ ( 0, -1) ( 0, 0) ( 0, 1)
+ (-1, -1) (-1, 0) (-1, 1)"""
- grass.message( "..." )
-
- terms2 = ["(elevation_step1[%d,%d] < elevation_step1)" % d
- for d in offsets2]
-
- expr2 = "PCTL2 = (100.0 / 48.0) * (%s)" % " + ".join(terms2)
-
- grass.mapcalc( expr2 )
- grass.run_command("g.remove", flags="f", type="raster", name= "elevation_step1",
- quiet = True)
-
- grass.message( "Step 2: Calculation of elevation percentile PCTL2 done" )
- grass.message( "----" )
+ offsets = []
- # PVF2 for step 2:
- grass.message( "Step 2: Calculation of preliminary valley flatness index PVF2" )
- grass.message( "..." )
-
- grass.mapcalc("PVF2 = F2 * (1.0 / (1.0 + pow ( ( PCTL2 / 0.4 ), 3 ) ) )")
+ # 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):
+ offsets.append((i, j))
+
+ else:
- grass.message( "Step 2: Calculation of preliminary valley flatness index PVF2 done." )
- grass.message( "----" )
+ for i in range(-radius, radius+1):
+ for j in range(-radius, radius+1):
+ row = i + radius
+ col = j + radius
- # Calculation of the valley flatness step 2 VF2
- grass.message( "Step 2: Calculation of valley flatness VF2" )
- grass.message( "..." )
-
- grass.mapcalc("VF2 = 1 - (1.0 / (1.0 + pow ( ( PVF2 / 0.3 ), 4 ) ) )")
-
- grass.message( "Step 2: Calculation of valley flatness VF2 done." )
- grass.message( "----" )
-
- # Calculation of weight W2
- grass.message( "Step 2: Calculation of weight W2" )
- grass.message( "..." )
+ if pow(row - radius, 2) + pow(col - radius, 2) <= \
+ pow(radius, 2) and (i, j) != (0,0):
+ offsets.append((j, i))
- grass.mapcalc("W2 = 1 - (1.0 / (1.0 + pow ( ( VF2 / 0.4 ), 6.68 ) ) )")
-
- grass.message( "Step 2: Calculation of weight W2 done." )
- grass.message( "----" )
-
- # Calculation of MRVBF2
- grass.message( "Step 2: Calculation of MRVBF2" )
- grass.message( "..." )
+ return offsets
- grass.mapcalc("MRVBF2 = (W2 * (1 + VF2)) + ((1 - W2) * VF1)")
-
- grass.message( "Step 2: Calculation of MRVBF2 done." )
- grass.message( "----" )
- # Calculation of combined flatness index CF2
- grass.message( "Step 2: Calculation of combined flatness index CF2" )
- grass.message( "..." )
+def get_percentile(L, input, radius=3, window_square=False):
+ """Calculates the percentile whichj is the ratio of the number of points of
+ lower elevation to the total number of points in the surrounding region
- grass.mapcalc("CF2 = F1 * F2")
-
- grass.message( "Step 2: Calculation of combined flatness index CF2 done." )
- grass.message( "----" )
+ Args
+ ----
+ L : int
+ Processing step (level)
+ input : str
+ GRASS raster map (elevation) to perform calculation on
+ radius : int
+ Neighborhood radius (in pixels)
+ window_square : bool. Optional (default is False)
+ Boolean to use square or circular neighborhood
+
+ Returns
+ -------
+ PCTL : str
+ Name of GRASS cell-padded raster map with elevation percentile for
+ processing step L"""
- # step 2 finished
- grass.message( "Step 2 calculations done." )
- grass.message( "----" )
+ PCTL = "PCTL{L}".format(L=L+1)
+ TMP_RAST.append(PCTL)
- ##################################################################################
- # remaining steps
- grass.message( "two remaining steps" )
- grass.message( "----" )
-
- ##################################################################################
- # Step 3
- grass.message( "Step 3" )
- # DEM smoothing 11 x 11 windows with Gaussian smoothing kernel (sigma) 3 x 3 -> change resolution (3 x base resolution) -> calculations
- grass.message( "----" )
- # DEM smoothing
- grass.message( "Step 3: DEM smoothing 11 x 11 windows with Gaussian smoothing kernel (sigma) 3" )
- grass.message( "..." )
- grass.run_command('r.neighbors', input = r_elevation,
- output = "DEM_smoothed_step3",
- size = 11,
- gauss=3)
-
- grass.message( "Step 3 DEM smoothing 11 x 11 windows with Gaussian smoothing kernel (sigma) 3 done." )
- grass.message( "----" )
+ # cell padding of input raster
+ input_grown = '_'.join(['tmp', input.split('@')[0], 'cell_padded']) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ region = cell_padding(input=input, output=input_grown, radius=radius)
+
+ # get offsets for given neighborhood radius
+ offsets = focal_expr(radius=radius, window_square=window_square)
+
+ # generate grass mapcalc terms and execute
+ n_pixels = float(len(offsets))
+ terms = ["({input}[%d,%d] <= {input})".format(input=input_grown) % d for d in offsets]
+ expr = "{x} = ({s}) / {n}".format(x=PCTL, s=" + ".join(terms), n=n_pixels)
+ r.mapcalc(expr)
- # Calculate slope step 3
- grass.message( "Step 3: Calculation of slope" )
- grass.message( "..." )
- grass.run_command('r.resamp.rst', input = "DEM_smoothed_step3",
- slope = "DEM_smoothed_step3_slope_degree",
- ew_res = Xres_base,
- ns_res = Yres_base,
- quiet = True)
+ # remove padded map
+ g.remove(type='raster', name=input_grown, flags='f', quiet=True)
- grass.message( "..." )
-
- grass.mapcalc("$outmap = 100 * tan($slope) ",
- outmap = "r_slope_step3",
- slope = "DEM_smoothed_step3_slope_degree")
+ # return to original region setting
+ Region.write(region)
- grass.message( "..." )
-
- grass.run_command('g.remove', flags='f', type='raster', name= "DEM_smoothed_step3_slope_degree",
- quiet = True)
- grass.message( "Step 3: Calculation of slope done." )
- grass.message( "----" )
+ return(PCTL)
- # change resolution
- grass.message( "Step 3: change to resolution of 3 x base resolution" )
- grass.run_command('g.region', res = Xres_step3,
- flags = 'a',
- raster = r_elevation)
- current_region_step3 = grass.region()
- Y_step3 = current_region_step3['nsres']
- X_step3 = current_region_step3['ewres']
- grass.message( "Resolution step 3: %s x %s" % (X_step3, Y_step3))
- grass.message( "----" )
+def get_slope(L, elevation):
+ """Calculation of r.slope.aspect using a cell-padded elevation"""
- # coarsening DEM to resolution step 3 and calculate PCTL3 step 3
+ # pad input dem by 2 cells
+ elevation_padded = '_'.join(['tmp', elevation.split('@')[0], 'cell_padded']) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ region = cell_padding(input=elevation, output=elevation_padded, radius=2)
- grass.message( "Step 3: Calculation of elevation percentile PCTL3" )
- grass.message( "..." )
- grass.mapcalc("DEM_smoothed_step3_coarsed = DEM_smoothed_step3")
+ # calculate slope
+ slope = 'tmp_slope_step{L}'.format(L=L+1) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(slope)
+ r.slope_aspect(
+ elevation=elevation_padded, slope=slope, format='percent', quiet=True)
+
+ # remove padded elevation
+ g.remove(type='raster', name=elevation_padded, flags='f', quiet=True)
- if moving_window_square :
+ # return region to original (non-padded) settings
+ Region.write(region)
- offsets3 = [d
- for j in xrange(1,6+1)
- for i in [j,-j]
- for d in [(i,0),(i,1),(i,2),(i,3),(i,4),(i,5),(i,6),(i,-1),(i,-2),(i,-3),(i,-4),(i,-5),(i,-6)]]
+ return(slope)
- grass.message( "..." )
-
- terms3 = ["(DEM_smoothed_step3_coarsed[%d,%d] < DEM_smoothed_step3_coarsed)" % d
- for d in offsets3]
-
- expr3 = "PCTL3 = (100.0 / 168.0) * (%s)" % " + ".join(terms3)
-
- else :
-
- offsets3 = [d
- for j in xrange(1,6+1)
- for i in [j,-j]
- for d in [(i,0),(0,i),(i,i),(i,-i)]]
- grass.message( "..." )
-
- terms3 = ["(DEM_smoothed_step3_coarsed[%d,%d] < DEM_smoothed_step3_coarsed)" % d
- for d in offsets3]
-
- expr3 = "PCTL3 = (100.0 / 48.0) * (%s)" % " + ".join(terms3)
+def get_flatness(L, slope, t, p):
+ """Calculates the flatness index
+ Flatness F1 = 1 / (1 + pow ((slope / t), p)
+
+ Equation 2 (Gallant and Dowling, 2003)"""
- grass.mapcalc( expr3 )
+ F = "tmp_F{L}".format(L=L+1) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(F)
+ grass.mapcalc("$g = 1.0 / (1.0 + pow(($x / $t), $p))",
+ g=F, x=slope, t=t, p=p)
+
+ return F
- grass.message( "Step 3: Calculation of elevation percentile PCTL3 done." )
- grass.message( "----" )
- # switch back to base resolution
- grass.message( "Step 3: switch back to base resolution." )
+def get_prelim_flatness(L, F, PCTL, t, p):
+ """Transform elevation percentile to a local lowness value using
+ equation (1) and combined with flatness F to produce the preliminary
+ valley flatness index (PVF) for the first step.
+
+ Equation 3 (Gallant and Dowling, 2003)"""
+
+ PVF = "tmp_PVF{L}" + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(PVF)
- grass.run_command('g.region', raster = r_elevation,
- align = r_elevation,
- flags = 'a')
- current_region_step3_switched_back = grass.region()
- Y_step3_switched_back = current_region_step3_switched_back['nsres']
- X_step3_switched_back = current_region_step3_switched_back['ewres']
- grass.message( "Resolution step 3: %s x %s" % (X_step3_switched_back, Y_step3_switched_back))
- grass.message( "----" )
+ grass.mapcalc("$g = $a * (1.0 / (1.0 + pow(($x / $t), $p)))",
+ g=PVF, a=F, x=PCTL, t=t, p=p)
+
+ return PVF
- # resample/refine PCTL3 to base resolution
- grass.message( "Step 3: refine PCTL3 to base resolution" )
- grass.message( "..." )
- grass.run_command('r.resamp.interp', input = "PCTL3",
- output = "PCTL3_refined_base_resolution",
- method = "bilinear")
-
- grass.message( "Step 3: refine PCTL3 to base resolution done." )
- grass.message( "----" )
+def get_prelim_flatness_rf(L, F, PCTL, t, p):
+ """Transform elevation percentile to a local upness value using
+ equation (1) and combined with flatness to produce the preliminary
+ valley flatness index (PVF) for the first step
+
+ Equation 3 (Gallant and Dowling, 2003)"""
+
+ PVF = "tmp_PVF_MRRTF{L}".format(L=L+1) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(PVF)
- # calculate flatness F3 at the base resolution
- grass.message( "Step 3: calculate F3 at base resolution" )
- grass.message( "..." )
+ grass.mapcalc("$g = $a * (1.0 / (1.0 + pow(((1-$x) / $t), $p)))",
+ g=PVF, a=F, x=PCTL, t=t, p=p)
+
+ return PVF
- grass.mapcalc("F3 = 1.0 / (1.0 + pow( ( $slope / 8.0 ) ,4 ) )",
- slope = "r_slope_step3")
- grass.message( "Step 3: calculate F3 at base resolution done." )
- grass.message( "----" )
+def get_valley_flatness(L, PVF, t, p):
+ """Calculation of the valley flatness step VF
+ Larger values of VF1 indicate increasing valley bottom character
+ with values less than 0.5 considered not to be in valley bottoms
+
+ Equation 4 (Gallant and Dowling, 2003)"""
- # calculate accumulated combined flatness CF3 at the base resolution
- grass.message( "Step 3: calculate combined flatness CF3 at base resolution" )
- grass.message( "..." )
+ VF = "tmp_VF{L}".format(L=L+1) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(VF)
+ grass.mapcalc("$g = 1 - (1.0 / (1.0 + pow(($x / $t), $p)))",
+ g=VF, x=PVF, t=t, p=p)
+
+ return VF
- grass.mapcalc("CF3 = CF2 * F3")
-
- grass.message( "Step 3: calculate combined flatness CF3 at base resolution done." )
- grass.message( "----" )
- # calculate preliminary valley flatness index PVF3 at the base resolution
- grass.message( "Step 3: calculate preliminary valley flatness index PVF3 at base resolution" )
- grass.message( "..." )
+def get_mrvbf(L, VF_Lminus1, VF_L, t):
+ """Calculation of the MRVBF index
+ Requires that L>1"""
- grass.mapcalc("PVF3 = CF3 * (1.0 / (1.0 + pow ( ( PCTL3_refined_base_resolution / 0.4 ), 3 ) ) )")
-
- grass.message( "Step 3: calculate preliminary valley flatness index PVF3 at base resolution done." )
- grass.message( "----" )
+ W = "tmp_W{L}".format(L=L+1) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ MRVBF = "MRVBF{L}".format(L=L+1) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
- # calculate valley flatness index VF3
- grass.message( "Step 3: calculate valley flatness index VF3 at base resolution" )
- grass.message( "..." )
+ # Calculation of weight W2 (Equation 9)
+ TMP_RAST.append(W)
+ p = (math.log10(((L+1) - 0.5) / 0.1)) / math.log10(1.5)
+ grass.mapcalc("$g = 1 - (1.0 / (1.0 + pow(($x / $t), $p)))",
+ g=W, x=VF_L, t=t, p=p)
+
+ # Calculation of MRVBF2 (Equation 8)
+ TMP_RAST.append(MRVBF)
+ grass.mapcalc("$MBF = ($W * ($L + $VF)) + ((1 - $W) * $VF1)",
+ MBF=MRVBF, L=L, W=W, VF=VF_L, VF1=VF_Lminus1)
+
+ return MRVBF
- grass.mapcalc("VF3 = 1 - (1.0 / (1.0 + pow ( ( PVF3 / 0.3 ), 4 ) ) )")
-
- grass.message( "Step 3: calculate valley flatness index VF3 at base resolution done." )
- grass.message( "----" )
- # Calculation of weight W3
- grass.message( "Step 3: Calculation of weight W3" )
- grass.message( "..." )
+def get_combined_flatness(L, F1, F2):
+ """Calculates the combined flatness index
- p3 = (math.log10((3 - 0.5) / 0.1)) / math.log10(1.5)
-
- grass.mapcalc("W3 = 1 - (1.0 / (1.0 + pow ( ( VF3 / 0.4 ), $p ) ) )", p = p3)
-
- grass.message( "Step 3: Calculation of weight W3 done." )
- grass.message( "----" )
+ Equation 13 (Gallant and Dowling, 2003)"""
- # Calculation of MRVBF3
- grass.message( "Step 3: Calculation of MRVBF3" )
- grass.message( "..." )
+ CF = "tmp_CF{L}".format(L=L+1) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(CF)
+ grass.mapcalc("$CF = $F1 * $F2", CF=CF, F1=F1, F2=F2)
- grass.mapcalc("MRVBF3 = (W3 * (3 - 1 + VF3)) + ((1 - W3) * MRVBF2)")
-
- grass.message( "Step 3: Calculation of MRVBF3 done." )
- grass.message( "----" )
-
- ##################################################################################
- # Step 4
- grass.message( "Step 4" )
+ return CF
- # change resolution
- grass.message( "Step 4: change to resolution of step 3" )
- grass.run_command('g.region', res = Xres_step3,
- flags = 'a',
- raster = r_elevation)
- current_region_step3 = grass.region()
- Y_step3 = current_region_step3['nsres']
- X_step3 = current_region_step3['ewres']
- grass.message( "Resolution 1 for step 4: %s x %s" % (X_step3, Y_step3))
- grass.message( "----" )
+def get_smoothed_dem(L, DEM):
+ """Smooth the DEM using an 11 cell averaging filter with gauss weighting of 3 radius"""
+
+ smoothed = "tmp_DEM_smoothed_step{L}".format(L=L+1) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(smoothed)
- # DEM smoothing 11 x 11 windows with Gaussian smoothing kernel (sigma) 3 x 3 -> change resolution (3 x base resolution) -> calculations
- grass.message( "----" )
- # DEM smoothing
- grass.message( "Step 4: DEM smoothing 11 x 11 windows with Gaussian smoothing kernel (sigma) 3" )
- grass.message( "..." )
- grass.run_command('r.neighbors', input = "DEM_smoothed_step3_coarsed",
- output = "DEM_smoothed_step4",
- size = 11,
- gauss=3)
-
- grass.message( "Step 4 DEM smoothing 11 x 11 windows with Gaussian smoothing kernel (sigma) 3 done." )
- grass.message( "----" )
+ #g = 4.3565 * math.exp(-(3 / 3.0))
+ r.neighbors(input=DEM, output=smoothed, size=11, gauss=3)
+
+ return smoothed
- # Calculate slope step 4
- grass.message( "Step 4: Calculation of slope" )
- grass.message( "..." )
- grass.run_command('r.resamp.rst', input = "DEM_smoothed_step4",
- slope = "DEM_smoothed_step4_slope_degree",
- ew_res = X_step3,
- ns_res = Y_step3,
- quiet = True)
- grass.message( "..." )
-
- grass.mapcalc("$outmap = 100 * tan($slope) ",
- outmap = "r_slope_step4",
- slope = "DEM_smoothed_step4_slope_degree")
+def refine(input, region, method='bilinear'):
+ """change resolution back to base resolution and resample a raster"""
+
+ input_padded = '_'.join(['tmp', input, '_padded']) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(input_padded)
+ cell_padding(input=input, output=input_padded, radius=2)
+
+ Region.write(region)
+ input = '_'.join(['tmp', input, 'refined_base_resolution']) + \
+ ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
+ TMP_RAST.append(input)
- grass.message( "..." )
-
- grass.run_command('g.remove', flags='f', type='raster', name= "DEM_smoothed_step4_slope_degree",
- quiet = True)
- grass.message( "Step 4: Calculation of slope done." )
- grass.message( "----" )
+ if method == 'bilinear':
+ r.resamp_interp(input=input_padded,
+ output=input,
+ method="bilinear")
- # change resolution
- grass.message( "Step 4: change to resolution to 3 x step 3 resolution" )
- grass.run_command('g.region', res = Xres_step4,
- flags = 'a',
- raster = r_elevation)
+ if method == 'average':
+ r.resamp_stats(input=input_padded,
+ output=input,
+ method='average',
+ flags='w')
+
+ return input
- current_region_step4 = grass.region()
- Y_step4 = current_region_step4['nsres']
- X_step4 = current_region_step4['ewres']
- grass.message( "Resolution step 4: %s x %s" % (X_step4, Y_step4))
- grass.message( "----" )
- # coarsening DEM to resolution step 4 and calculate PCTL step 4
+def main():
+ r_elevation = options['elevation']
+ mrvbf = options['mrvbf'].split('@')[0]
+ mrrtf = options['mrrtf'].split('@')[0]
+ t_slope = float(options['t_slope'])
+ t_pctl_v = float(options['t_pctl_v'])
+ t_pctl_r = float(options['t_pctl_r'])
+ t_vf = float(options['t_rf'])
+ t_rf = float(options['t_rf'])
+ p_slope = float(options['p_slope'])
+ p_pctl = float(options['p_pctl'])
+ moving_window_square = flags['s']
+ levels = int(options['levels'])
- grass.message( "Step 4: Calculation of elevation percentile PCTL4" )
- grass.message( "..." )
- grass.mapcalc("DEM_smoothed_step4_coarsed = DEM_smoothed_step4")
-
- if moving_window_square :
-
- offsets4 = [d
- for j in xrange(1,6+1)
- for i in [j,-j]
- for d in [(i,0),(i,1),(i,2),(i,3),(i,4),(i,5),(i,6),(i,-1),(i,-2),(i,-3),(i,-4),(i,-5),(i,-6)]]
-
- grass.message( "..." )
-
- terms4 = ["(DEM_smoothed_step4_coarsed[%d,%d] < DEM_smoothed_step4_coarsed)" % d
- for d in offsets4]
-
- expr4 = "PCTL4 = (100.0 / 168.0) * (%s)" % " + ".join(terms4)
+ ###########################################################################
+ # Some checks
+ if levels < 3:
+ grass.fatal('Number of generalization steps (levels) cannot be < 3')
+ if t_slope <=0 or t_pctl_v <=0 or t_pctl_r <=0 or t_vf <=0 or t_rf <=0 or \
+ p_slope <=0 or p_pctl <=0:
+ grass.fatal('Parameter values cannot be <= 0')
+ if levels > 10:
+ grass.warning('A large number (>10) processing steps are selected, recommended is between 3 to 8')
- else :
+ ###########################################################################
+ # Intermediate outputs
+ Xres_step, Yres_step, DEM = [], [], []
+ slope, F, PCTL, PVF, PVF_RF = [0]*levels, [0]*levels, [0]*levels, [0]*levels, [0]*levels
+ VF, VF_RF, MRVBF, MRRTF = [0]*levels, [0]*levels, [0]*levels, [0]*levels
- offsets4 = [d
- for j in xrange(1,6+1)
- for i in [j,-j]
- for d in [(i,0),(0,i),(i,i),(i,-i)]]
+ ###########################################################################
+ # Step 1 (L=0)
+ # Base scale resolution
+ L = 0
+ Xres_step.append(current_region.ewres)
+ Yres_step.append(current_region.nsres)
+ DEM.append(r_elevation)
+ radi = 3
- grass.message( "..." )
-
- terms4 = ["(DEM_smoothed_step4_coarsed[%d,%d] < DEM_smoothed_step4_coarsed)" % d
- for d in offsets4]
-
- expr4 = "PCTL4 = (100.0 / 48.0) * (%s)" % " + ".join(terms4)
-
- grass.mapcalc( expr4 )
+ g.message(os.linesep)
+ g.message("Step {L}".format(L=L+1))
+ g.message("------")
- grass.message( "Step 4: Calculation of elevation percentile PCTL4 done." )
- grass.message( "----" )
+ # Calculation of slope (S1) and calculation of flatness (F1) (Equation 2)
+ grass.message("Calculation of slope and transformation to flatness F{L}...".format(L=L+1))
+ slope[L] = get_slope(L, DEM[L])
+ F[L] = get_flatness(L, slope[L], t_slope, p_slope)
- # switch back to base resolution
- grass.message( "Step 4: switch back to base resolution." )
+ # Calculation of elevation percentile PCTL for step 1
+ grass.message("Calculation of elevation percentile PCTL{L}...".format(L=L+1))
+ PCTL[L] = get_percentile(L, DEM[L], radi, moving_window_square)
- grass.run_command('g.region', raster = r_elevation,
- align = r_elevation,
- flags = 'a')
- current_region_step4_switched_back = grass.region()
- Y_step4_switched_back = current_region_step4_switched_back['nsres']
- X_step4_switched_back = current_region_step4_switched_back['ewres']
- grass.message( "Resolution step 4: %s x %s" % (X_step4_switched_back, Y_step4_switched_back))
- grass.message( "----" )
+ # Transform elevation percentile to local lowness for step 1 (Equation 3)
+ grass.message("Calculation of preliminary valley flatness index PVF{L}...".format(L=L+1))
+ PVF[L] = get_prelim_flatness(L, F[L], PCTL[L], t_pctl_v, p_pctl)
+ if mrrtf != '':
+ grass.message("Calculation of preliminary ridge top flatness index PRF{L}...".format(L=L+1))
+ PVF_RF[L] = get_prelim_flatness_rf(L, F[L], PCTL[L], t_pctl_r, p_pctl)
- # resample/refine PCTL4 to base resolution
- grass.message( "Step 4: refine PCTL4 to base resolution" )
- grass.message( "..." )
+ # Calculation of the valley flatness step 1 VF1 (Equation 4)
+ grass.message("Calculation of valley flatness VF{L}...".format(L=L+1))
+ VF[L] = get_valley_flatness(L, PVF[L], t_vf, p_slope)
+ if mrrtf != '':
+ grass.message("Calculation of ridge top flatness RF{L}...".format(L=L+1))
+ VF_RF[L] = get_valley_flatness(L, PVF_RF[L], t_rf, p_slope)
- grass.run_command('r.resamp.interp', input = "PCTL4",
- output = "PCTL4_refined_base_resolution",
- method = "bilinear")
+ ##################################################################################
+ # Step 2 (L=1)
+ # Base scale resolution
+ L = 1
+ Xres_step.append(current_region.ewres)
+ Yres_step.append(current_region.nsres)
+ DEM.append(r_elevation)
+ t_slope /= 2.0
+ radi = 6
+
+ grass.message(os.linesep)
+ grass.message("Step {L}".format(L=L+1))
+ grass.message("------")
- grass.message( "Step 4: refine PCTL4 to base resolution done." )
- grass.message( "----" )
+ # Calculation of flatness for step 2 (Equation 5)
+ # The second step commences the same way with the original DEM at its base resolution,
+ # using a slope threshold ts,2 half of ts,1:
+ grass.message("Calculation of flatness F{L}...".format(L=L+1))
+ F[L] = get_flatness(L, slope[L-1], t_slope, p_slope)
- # calculate flatness F4 at the base resolution
- grass.message( "Step 4: calculate F4 at base resolution" )
- grass.message( "..." )
+ # Calculation of elevation percentile PCTL for step 2 (radius of 6 cells)
+ grass.message("Calculation of elevation percentile PCTL{L}...".format(L=L+1))
+ PCTL[L] = get_percentile(L, r_elevation, radi, moving_window_square)
- grass.run_command('r.resamp.interp', input = "r_slope_step4",
- output = "r_slope_step4_refined_base_resolution",
- method = "bilinear")
+ # PVF for step 2 (Equation 6)
+ grass.message("Calculation of preliminary valley flatness index PVF{L}...".format(L=L+1))
+ PVF[L] = get_prelim_flatness(L, F[L], PCTL[L], t_pctl_v, p_pctl)
+ if mrrtf != '':
+ grass.message("Calculation of preliminary ridge top flatness index PRF{L}...".format(L=L+1))
+ PVF_RF[L] = get_prelim_flatness_rf(L, F[L], PCTL[L], t_pctl_r, p_pctl)
+
+ # Calculation of the valley flatness VF for step 2 (Equation 7)
+ grass.message("Calculation of valley flatness VF{L}...".format(L=L+1))
+ VF[L] = get_valley_flatness(L, PVF[L], t_vf, p_slope)
+ if mrrtf != '':
+ grass.message("Calculation of ridge top flatness RF{L}...".format(L=L+1))
+ VF_RF[L] = get_valley_flatness(L, PVF_RF[L], t_rf, p_slope)
+
+ # Calculation of MRVBF for step 2
+ grass.message("Calculation of MRVBF{L}...".format(L=L+1))
+ MRVBF[L] = get_mrvbf(L, VF_Lminus1=VF[L-1], VF_L=VF[L], t=t_pctl_v)
+ if mrrtf != '':
+ grass.message("Calculation of MRRTF{L}...".format(L=L+1))
+ MRRTF[L] = get_mrvbf(L, VF_Lminus1=VF_RF[L-1], VF_L=VF_RF[L], t=t_pctl_r)
+
+ # Update flatness for step 2 with combined flatness from F1 and F2 (Equation 10)
+ grass.message("Calculation of combined flatness index CF{L}...".format(L=L+1))
+ F[L] = get_combined_flatness(L, F[L-1], F[L])
- grass.mapcalc("F4 = 1.0 / (1.0 + pow( ( $slope / 4.0 ) ,4 ) )",
- slope = "r_slope_step4_refined_base_resolution")
+ ##################################################################################
+ # Remaining steps
+ # DEM_1_1 refers to scale (smoothing) and resolution (cell size)
+ # so that DEM_L1_L-1 refers to smoothing of current step,
+ # but resolution of previous step
- grass.message( "Step 4: calculate F4 at base resolution done." )
- grass.message( "----" )
+ for L in range(2, levels):
- # calculate accumulated combined flatness CF3 at the base resolution
- grass.message( "Step 4: calculate combined flatness CF4 at base resolution" )
- grass.message( "..." )
+ t_slope /= 2.0
+ Xres_step.append(Xres_step[L-1] * 3)
+ Yres_step.append(Yres_step[L-1] * 3)
+ radi = 6
- grass.mapcalc("CF4 = CF3 * F4")
-
- grass.message( "Step 4: calculate combined flatness CF4 at base resolution done." )
- grass.message( "----" )
+ grass.message(os.linesep)
+ grass.message("Step {L}".format(L=L+1))
+ grass.message("------")
+
+ # Coarsen resolution to resolution of prevous step (step L-1) and smooth DEM
+ if L >= 3:
+ grass.run_command('g.region', ewres = Xres_step[L-1], nsres = Yres_step[L-1])
+ grass.message('Coarsening resolution to ew_res={e} and ns_res={n}...'.format(
+ e=Xres_step[L-1], n=Yres_step[L-1]))
+
+ grass.message("DEM smoothing 11 x 11 windows with Gaussian smoothing kernel (sigma) 3...")
+ DEM.append(get_smoothed_dem(L, DEM[L-1]))
- # calculate preliminary valley flatness index PVF3 at the base resolution
- grass.message( "Step 4: calculate preliminary valley flatness index PVF4 at base resolution" )
- grass.message( "..." )
+ # Calculate slope
+ grass.message("Calculation of slope...")
+ slope[L] = get_slope(L, DEM[L])
- grass.mapcalc("PVF4 = CF4 * (1.0 / (1.0 + pow ( ( PCTL4_refined_base_resolution / 0.4 ), 3 ) ) )")
-
- grass.message( "Step 4: calculate preliminary valley flatness index PVF4 at base resolution done." )
- grass.message( "----" )
+ # Refine slope to base resolution
+ if L >= 3:
+ grass.message('Resampling slope back to base resolution...')
+ slope[L] = refine(slope[L], current_region, method='bilinear')
- # calculate valley flatness index VF4
- grass.message( "Step 4: calculate valley flatness index VF4 at base resolution" )
- grass.message( "..." )
+ # Coarsen resolution to current step L and calculate PCTL
+ grass.run_command('g.region', ewres=Xres_step[L], nsres=Yres_step[L])
+ DEM[L] = refine(DEM[L], Region(), method = 'average')
+ grass.message("Calculation of elevation percentile PCTL{L}...".format(L=L+1))
+ PCTL[L] = get_percentile(L, DEM[L], radi, moving_window_square)
+
+ # Refine PCTL to base resolution
+ grass.message("Resampling PCTL{L} to base resolution...".format(L=L+1))
+ PCTL[L] = refine(PCTL[L], current_region, method='bilinear')
- grass.mapcalc("VF4 = 1 - (1.0 / (1.0 + pow ( ( PVF4 / 0.3 ), 4 ) ) )")
-
- grass.message( "Step 4: calculate valley flatness index VF4 at base resolution done." )
- grass.message( "----" )
+ # Calculate flatness F at the base resolution
+ grass.message("Calculate F{L} at base resolution...".format(L=L+1))
+ F[L] = get_flatness(L, slope[L], t_slope, p_slope)
- # Calculation of weight W4
- grass.message( "Step 4: Calculation of weight W4" )
- grass.message( "..." )
-
- p4 = (math.log10((4 - 0.5) / 0.1)) / math.log10(1.5)
-
- grass.mapcalc("W4 = 1 - (1.0 / (1.0 + pow ( ( VF4 / 0.4 ), $p ) ) )", p = p4)
-
- grass.message( "Step 4: Calculation of weight W4 done." )
- grass.message( "----" )
+ # Update flatness with combined flatness CF from the previous step
+ grass.message("Calculate combined flatness CF{L} at base resolution...".format(L=L+1))
+ F[L] = get_combined_flatness(L, F1=F[L-1], F2=F[L])
- # Calculation of MRVBF4
- grass.message( "Step 4: Calculation of MRVBF4" )
- grass.message( "..." )
+ # Calculate preliminary valley flatness index PVF at the base resolution
+ grass.message("Calculate preliminary valley flatness index PVF{L} at base resolution...".format(L=L+1))
+ PVF[L] = get_prelim_flatness(L, F[L], PCTL[L], t_pctl_v, p_pctl)
+ if mrrtf != '':
+ grass.message("Calculate preliminary ridge top flatness index PRF{L} at base resolution...".format(L=L+1))
+ PVF_RF[L] = get_prelim_flatness_rf(L, F[L], PCTL[L], t_pctl_r, p_pctl)
+
+ # Calculate valley flatness index VF
+ grass.message("Calculate valley flatness index VF{L} at base resolution...".format(L=L+1))
+ VF[L] = get_valley_flatness(L, PVF[L], t_vf, p_slope)
+ if mrrtf != '':
+ grass.message("Calculate ridge top flatness index RF{L} at base resolution...".format(L=L+1))
+ VF_RF[L] = get_valley_flatness(L, PVF_RF[L], t_rf, p_slope)
+
+ # Calculation of MRVBF
+ grass.message("Calculation of MRVBF{L}...".format(L=L+1))
+ MRVBF[L] = get_mrvbf(L, VF_Lminus1=MRVBF[L-1], VF_L=VF[L], t=t_pctl_v)
+ if mrrtf != '':
+ grass.message("Calculation of MRRTF{L}...".format(L=L+1))
+ MRRTF[L] = get_mrvbf(L, VF_Lminus1=MRRTF[L-1], VF_L=VF_RF[L], t=t_pctl_r)
- grass.mapcalc("MRVBF4 = (W4 * (3 - 1 + VF4)) + ((1 - W4) * MRVBF3)")
-
- grass.message( "Step 4: Calculation of MRVBF4 done." )
- grass.message( "----" )
+ # Output final MRVBF
+ grass.mapcalc("$x = $y", x = mrvbf, y=MRVBF[L])
- # clean up some temporay files and maps
- grass.message( "Some clean up ..." )
- grass.run_command('g.remove', flags = 'f', type = 'region', name = "base_region_MrVBF",
- quiet = True)
- grass.run_command('g.remove', flags = 'f', type = 'raster', name = ["F1,F2,F3,F4"],
- quiet = True)
- grass.run_command('g.remove', flags = 'f', type = 'raster', name = ["CF2,CF3,CF4"],
- quiet = True)
- grass.run_command('g.remove', flags = 'f', type = 'raster', name = ["PCTL1,PCTL2,PCTL3,PCTL3_refined_base_resolution,PCTL4,PCTL4_refined_base_resolution"],
- quiet = True)
- grass.run_command('g.remove', flags = 'f', type = 'raster', name = ["PVF1,PVF2,PVF3,PVF4"],
- quiet = True)
- grass.run_command('g.remove', flags = 'f', type = 'raster', name = ["VF1,VF2,VF3,VF4"],
- quiet = True)
- grass.run_command('g.remove', flags = 'f', type = 'raster', name = ["W2,W3,W4"],
- quiet = True)
- grass.run_command('g.remove', flags = 'f', type = 'raster', name = ["r_slope_step1,r_slope_step3,r_slope_step4,r_slope_step4_refined_base_resolution"],
- quiet = True)
- grass.run_command('g.remove', flags = 'f', type = 'raster', name = ["DEM_smoothed_step3,DEM_smoothed_step3_coarsed,DEM_smoothed_step4,DEM_smoothed_step4_coarsed"],
- quiet = True)
-
- grass.message( "Clean up done." )
- grass.message( "----" )
-
- # v.habitat.dem done!
- grass.message( "r.valley.bottom done!" )
+ if mrrtf != '':
+ grass.mapcalc("$x = $y", x = mrrtf, y=MRRTF[L])
-
if __name__ == "__main__":
options, flags = grass.parser()
+ atexit.register(cleanup)
sys.exit(main())
More information about the grass-commit
mailing list