[GRASS-SVN] r64747 - grass-addons/grass7/raster/r.catchment
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Feb 25 11:53:35 PST 2015
Author: martinl
Date: 2015-02-25 11:53:35 -0800 (Wed, 25 Feb 2015)
New Revision: 64747
Added:
grass-addons/grass7/raster/r.catchment/r.catchment.py
Removed:
grass-addons/grass7/raster/r.catchment/r.catchment
Log:
r.catchment: fix compilation issue ( adds extension .py)
Deleted: grass-addons/grass7/raster/r.catchment/r.catchment
===================================================================
--- grass-addons/grass7/raster/r.catchment/r.catchment 2015-02-25 19:51:21 UTC (rev 64746)
+++ grass-addons/grass7/raster/r.catchment/r.catchment 2015-02-25 19:53:35 UTC (rev 64747)
@@ -1,320 +0,0 @@
-#!/usr/bin/env python
-#
-############################################################################
-#
-# MODULE: r.catchment
-# AUTHOR(S): Isaac Ullah, Arizona State University
-# PURPOSE: Creates a raster buffer of specified area around vector points
-# using cost distances. Module requires r.walk.
-# ACKNOWLEDGEMENTS: National Science Foundation Grant #BCS0410269
-# COPYRIGHT: (C) 2015 by Isaac Ullah, Arizona State University
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
-#
-#############################################################################
-
-
-#%Module
-#% description: Creates a raster buffer of specified area around vector points using cost distances using r.walk. NOTE: please run g.region first to make sure region boundaries and resolution match input elevation map.
-#%END
-
-
-#%option
-#% key: elevation
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Input elevation map (DEM)
-#% required : yes
-#%END
-#%option
-#% key: in_cost
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Input cost map (This will override the input elevation map, if none specified, one will be created from input elevation map with r.walk)
-#% required : no
-#%END
-#%option
-#% key: start_points
-#% type: string
-#% gisprompt: new,vector,vector
-#% description: Name of input vector site points map
-#% required : yes
-#%END
-#%option
-#% key: friction
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Optional map of friction costs. If no map selected, default friction=0 making output reflect time costs only
-#% answer:
-#% required : no
-#%END
-#%option
-#% key: a
-#% type: double
-#% description: Coefficients for walking energy formula parameters a,b,c,d
-#% answer: 0.72
-#% required : no
-#%END
-#%option
-#% key: b
-#% type: double
-#% description:
-#% answer: 6.0
-#% required : no
-#%END
-#%option
-#% key: c
-#% type: double
-#% description:
-#% answer: 1.9998
-#% required : no
-#%END
-#%option
-#% key: d
-#% type: double
-#% description:
-#% answer: -1.9998
-#% required : no
-#%END
-#%option
-#% key: lambda
-#% type: double
-#% description: Lambda value for cost distance calculation (for combining friction costs with walking costs)
-#% answer: 1
-#% required : no
-#%END
-#%option
-#% key: slope_factor
-#% type: double
-#% description: Slope factor determines travel energy cost per height step
-#% answer: -0.2125
-#% required : no
-#%END
-#%option
-#% key: buffer
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Output buffer map
-#% required : yes
-#%END
-#%option
-#% key: sigma
-#% type: double
-#% description: Slope threshold for mask
-#% answer:
-#% required : no
-#%END
-#%option
-#% key: area
-#% type: integer
-#% description: Area of buffer (Integer value to nearest 100 square map units)
-#% answer: 5000000
-#% required : yes
-#%END
-#%option
-#% key: map_val
-#% type: integer
-#% description: Integer value for output catchment area (all other areas will be Null)
-#% answer: 1
-#% required : yes
-#%END
-#%flag
-#% key: k
-#% description: -k Use knight's move for calculating cost surface (slower but more accurate)
-#%END
-#%flag
-#% key: c
-#% description: -c Keep cost surface used to calculate buffers
-#%END
-#%flag
-#% key: l
-#% description: -l Show a list of all cost surface values and the area of the catchment that they delimit
-#%END
-
-
-import sys
-import os
-import subprocess
-import tempfile
-import random
-grass_install_tree = os.getenv('GISBASE')
-sys.path.append(grass_install_tree + os.sep + 'etc' + os.sep + 'python')
-import grass.script as grass
-# first define a useful custom method
-
-# m is a grass/bash command that will generate some list of keyed info to stdout where the keys are numeric values, n is the character that separates the key from the data, o is a defined blank dictionary to write results to
-def out2dictnum(m, n, o):
- p1 = subprocess.Popen('%s' % m, stdout=subprocess.PIPE, shell='bash')
- p2 = p1.stdout.readlines()
- for y in p2:
- y0,y1 = y.split('%s' % n)
- y0num = float(y0)
- o[y0num] = y1.strip('\n')
-
-#main block of code starts here
-def main():
- pid = os.getpid()
- #setting up variables for use later on
- elevation = options["elevation"]
- start_points = options["start_points"]
- lmbda = options["lambda"]
- slope_factor = options["slope_factor"]
- a = options["a"]
- b = options["b"]
- c = options["c"]
- d = options["d"]
- sigma = options["sigma"]
- area = float(options["area"])
- buffer = options["buffer"]
- mapval = options["map_val"]
- w_coefs = a + ',' + b + ',' + c + ',' + d
- if "MASK" in grass.list_grouped('rast')[grass.gisenv()['MAPSET']] and bool(options["sigma"]) is True:
- grass.message('There is already a MASK in place, and you have also selected to mask slope values above %s.\n The high slope areas (slope mask) will be temporarily added to current MASKED areas for the calcualtion of the catchment geometry.\n The original MASK will be restored when the module finishes' % sigma)
- ismask = 2
- elif "MASK" in grass.list_grouped('rast')[grass.gisenv()['MAPSET']]:
- grass.message('There is a MASK in place. The areas MASKed out will be ignored while calculating catchment geometry.')
- ismask = 1
- else:
- ismask = 0
-
- grass.message("Wanted buffer area=%s\n" % int(area))
-
-####################################################
- if bool(options["in_cost"]) is True:
- grass.message('\n\nUsing input cost surface\n')
- cost = options["in_cost"]
- else:
- grass.message('\n\nstep 1 of 4: Calculating cost surface\n')
- cost = 'temporary.cost.%s' % pid
- if bool(options["friction"]) is True:
- grass.message('Calculating costs using input friction map\n')
- friction = options["friction"]
- else:
- grass.message('Calculating for time costs only')
- friction = "temporary.friction.%s" % pid
- grass.mapcalc("${out} = if(isnull(${rast1}), null(), 0)", overwrite = grass.overwrite(), quiet = True, out = friction, rast1 = elevation)
- if flags["k"] is True:
- grass.message('Using Knight\'s move\n')
- #NOTE! because "lambda" is an internal python variable, it is impossible to enter the value for key "lambda" in r.walk. It ends up with a python error.
- grass.run_command('r.walk', quiet = True, overwrite = grass.overwrite(), flags = 'k', elevation = elevation, friction = friction, output = cost, start_points = start_points, walk_coeff = w_coefs, slope_factor = slope_factor)
- else:
- grass.run_command('r.walk', quiet = True, overwrite = grass.overwrite(), elevation = elevation, friction = friction, output = cost, start_points = start_points, percent_memory = '100', walk_coeff = w_coefs, slope_factor = slope_factor)
- if bool(options["friction"]) is False:
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = friction)
-#################################################
- if bool(options["sigma"]) is True:
- grass.message('\n\nCreating optional slope mask\n')
- slope = "temporary.slope.%s" % pid
- grass.run_command('r.slope.aspect', quiet = True, overwrite = grass.overwrite(), elevation = elevation, slope = slope)
- if ismask == 2:
- grass.mapcalc("MASK=if(${rast1} <= ${sigma}, 1, if(${tempmask}, 1, null()))", overwrite = grass.overwrite(), quiet = True, sigma = sigma, rast1 = slope, tempmask = tempmask)
- else:
- grass.mapcalc("MASK=if(${rast1} <= ${sigma}, 1, null())", overwrite = grass.overwrite(), quiet = True, sigma = sigma, rast1 = slope)
- else:
- grass.message('No slope mask created')
-##################################################
- if flags["l"] is True:
- grass.message('\n\nCalculating list of possible catchment configurations\n')
- grass.message("cost value | catchment area")
- areadict = {}
- out2dictnum('r.stats -Aani input=' + cost + ' fs=, nv=* nsteps=255', ',', areadict)
- testarea = 0
- #start the loop, and list the values
- for key in sorted(areadict):
- testarea = testarea + int(float(areadict[key]))
- grass.message("%s | %s" % (int(key), testarea))
- if flags["c"] is True:
- if bool(options["in_cost"]) is False:
- grass.run_command('g.rename', overwrite = grass.overwrite(), quiet = True, rast = 'temporary.cost.%s,%s_cost_surface' % (pid, buffer))
- grass.message('Cleaning up...(keeping cost map)')
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
- else:
- grass.message('Cleaning up...1')
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
- else:
- if bool(options["in_cost"]) is False:
- grass.message('Cleaning up...2')
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s,temporary.cost.%s' % (pid, pid))
- else:
- grass.message('Cleaning up...3')
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
- if bool(options["sigma"]) is True:
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = slope)
- if ismask == 2:
- grass.message('Reinstating original MASK...')
- grass.run_command('g.rename', overwrite = grass.overwrite(), quiet = "True", rast = tempmask +',MASK')
- elif ismask == 0 and bool(options["sigma"]) is True:
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'MASK')
- elif ismask == 1:
- grass.message('Keeping original MASK')
- grass.message(' DONE!')
- return
- else:
- grass.message('\n\nCalculating buffer\n')
- areadict = {}
- out2dictnum('r.stats -Aani input=' + cost + ' fs=, nv=* nsteps=255', ',', areadict)
- tot_area = 0
- for key in sorted(areadict):
- tot_area = tot_area + int(float(areadict[key]))
- maxcost = key
- grass.message("Maximum cost distance value %s covers an area of %s square map units\n\nCommencing to find a catchment configuration.....\n\n" % (int(maxcost), tot_area))
- testarea = 0
- lastarea = 0
- lastkey = 0
- #start the loop, and home in on the target range
- for key in sorted(areadict):
- testarea = testarea + int(float(areadict[key]))
- if testarea >= area:
- break
- lastkey = key
- lastarea = testarea
- if (testarea - area) <= (area - lastarea):
- cutoff = key
- displayarea = testarea
- else:
- cutoff = lastkey
- displayarea = lastarea
- grass.message("Catchment configuration found! Cost cutoff %s produces a catchment of %s square map units." % (int(cutoff), displayarea))
- ####################################################
- grass.message('\n\nCreating output map\n')
- temp = tempfile.NamedTemporaryFile()
- temp.write('0 thru %s = %s\n' % (int(cutoff), mapval))
- temp.flush()
- grass.run_command('r.reclass', overwrite = grass.overwrite(), input = cost, output = 'cost.reclass.%s' % pid, rules = temp.name)
- temp.close()
- grass.mapcalc("${out}=if(isnull(${cost}), null(), ${cost})", overwrite = grass.overwrite(), quiet = True, cost = "cost.reclass.%s" % pid, out = buffer)
- grass.message("\nThe output catchment map will be named %s" % buffer)
- grass.run_command('r.colors', quiet = True, map = buffer, color = 'ryb')
- if flags["c"] is True:
- if bool(options["in_cost"]) is False:
- grass.run_command('g.rename', overwrite = grass.overwrite(), quiet = True, rast = 'temporary.cost,%s_cost_surface' % (buffer))
- grass.message('Cleaning up...(keeping cost map)')
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
- else:
- grass.message('Cleaning up...1')
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
- else:
- if bool(options["in_cost"]) is False:
- grass.message('Cleaning up...2')
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s,temporary.cost.%s' % (pid, pid))
- else:
- grass.message('Cleaning up...3')
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
- if bool(options["sigma"]) is True:
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = slope)
- if ismask == 2:
- grass.message('Reinstating original MASK...')
- grass.run_command('g.rename', overwrite = grass.overwrite(), quiet = "True", rast = tempmask +',MASK')
- elif ismask == 0 and bool(options["sigma"]) is True:
- grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name ='MASK')
- elif ismask == 1:
- grass.message('Keeping original MASK')
- grass.message(' DONE!')
- return
-
-# here is where the code in "main" actually gets executed. This way of programming is neccessary for the way g.parser needs to run in GRASS 7.
-if __name__ == "__main__":
- options, flags = grass.parser()
- main()
- exit(0)
Copied: grass-addons/grass7/raster/r.catchment/r.catchment.py (from rev 64744, grass-addons/grass7/raster/r.catchment/r.catchment)
===================================================================
--- grass-addons/grass7/raster/r.catchment/r.catchment.py (rev 0)
+++ grass-addons/grass7/raster/r.catchment/r.catchment.py 2015-02-25 19:53:35 UTC (rev 64747)
@@ -0,0 +1,320 @@
+#!/usr/bin/env python
+#
+############################################################################
+#
+# MODULE: r.catchment
+# AUTHOR(S): Isaac Ullah, Arizona State University
+# PURPOSE: Creates a raster buffer of specified area around vector points
+# using cost distances. Module requires r.walk.
+# ACKNOWLEDGEMENTS: National Science Foundation Grant #BCS0410269
+# COPYRIGHT: (C) 2015 by Isaac Ullah, Arizona State University
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+
+#%Module
+#% description: Creates a raster buffer of specified area around vector points using cost distances using r.walk. NOTE: please run g.region first to make sure region boundaries and resolution match input elevation map.
+#%END
+
+
+#%option
+#% key: elevation
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Input elevation map (DEM)
+#% required : yes
+#%END
+#%option
+#% key: in_cost
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Input cost map (This will override the input elevation map, if none specified, one will be created from input elevation map with r.walk)
+#% required : no
+#%END
+#%option
+#% key: start_points
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Name of input vector site points map
+#% required : yes
+#%END
+#%option
+#% key: friction
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Optional map of friction costs. If no map selected, default friction=0 making output reflect time costs only
+#% answer:
+#% required : no
+#%END
+#%option
+#% key: a
+#% type: double
+#% description: Coefficients for walking energy formula parameters a,b,c,d
+#% answer: 0.72
+#% required : no
+#%END
+#%option
+#% key: b
+#% type: double
+#% description:
+#% answer: 6.0
+#% required : no
+#%END
+#%option
+#% key: c
+#% type: double
+#% description:
+#% answer: 1.9998
+#% required : no
+#%END
+#%option
+#% key: d
+#% type: double
+#% description:
+#% answer: -1.9998
+#% required : no
+#%END
+#%option
+#% key: lambda
+#% type: double
+#% description: Lambda value for cost distance calculation (for combining friction costs with walking costs)
+#% answer: 1
+#% required : no
+#%END
+#%option
+#% key: slope_factor
+#% type: double
+#% description: Slope factor determines travel energy cost per height step
+#% answer: -0.2125
+#% required : no
+#%END
+#%option
+#% key: buffer
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Output buffer map
+#% required : yes
+#%END
+#%option
+#% key: sigma
+#% type: double
+#% description: Slope threshold for mask
+#% answer:
+#% required : no
+#%END
+#%option
+#% key: area
+#% type: integer
+#% description: Area of buffer (Integer value to nearest 100 square map units)
+#% answer: 5000000
+#% required : yes
+#%END
+#%option
+#% key: map_val
+#% type: integer
+#% description: Integer value for output catchment area (all other areas will be Null)
+#% answer: 1
+#% required : yes
+#%END
+#%flag
+#% key: k
+#% description: -k Use knight's move for calculating cost surface (slower but more accurate)
+#%END
+#%flag
+#% key: c
+#% description: -c Keep cost surface used to calculate buffers
+#%END
+#%flag
+#% key: l
+#% description: -l Show a list of all cost surface values and the area of the catchment that they delimit
+#%END
+
+
+import sys
+import os
+import subprocess
+import tempfile
+import random
+grass_install_tree = os.getenv('GISBASE')
+sys.path.append(grass_install_tree + os.sep + 'etc' + os.sep + 'python')
+import grass.script as grass
+# first define a useful custom method
+
+# m is a grass/bash command that will generate some list of keyed info to stdout where the keys are numeric values, n is the character that separates the key from the data, o is a defined blank dictionary to write results to
+def out2dictnum(m, n, o):
+ p1 = subprocess.Popen('%s' % m, stdout=subprocess.PIPE, shell='bash')
+ p2 = p1.stdout.readlines()
+ for y in p2:
+ y0,y1 = y.split('%s' % n)
+ y0num = float(y0)
+ o[y0num] = y1.strip('\n')
+
+#main block of code starts here
+def main():
+ pid = os.getpid()
+ #setting up variables for use later on
+ elevation = options["elevation"]
+ start_points = options["start_points"]
+ lmbda = options["lambda"]
+ slope_factor = options["slope_factor"]
+ a = options["a"]
+ b = options["b"]
+ c = options["c"]
+ d = options["d"]
+ sigma = options["sigma"]
+ area = float(options["area"])
+ buffer = options["buffer"]
+ mapval = options["map_val"]
+ w_coefs = a + ',' + b + ',' + c + ',' + d
+ if "MASK" in grass.list_grouped('rast')[grass.gisenv()['MAPSET']] and bool(options["sigma"]) is True:
+ grass.message('There is already a MASK in place, and you have also selected to mask slope values above %s.\n The high slope areas (slope mask) will be temporarily added to current MASKED areas for the calcualtion of the catchment geometry.\n The original MASK will be restored when the module finishes' % sigma)
+ ismask = 2
+ elif "MASK" in grass.list_grouped('rast')[grass.gisenv()['MAPSET']]:
+ grass.message('There is a MASK in place. The areas MASKed out will be ignored while calculating catchment geometry.')
+ ismask = 1
+ else:
+ ismask = 0
+
+ grass.message("Wanted buffer area=%s\n" % int(area))
+
+####################################################
+ if bool(options["in_cost"]) is True:
+ grass.message('\n\nUsing input cost surface\n')
+ cost = options["in_cost"]
+ else:
+ grass.message('\n\nstep 1 of 4: Calculating cost surface\n')
+ cost = 'temporary.cost.%s' % pid
+ if bool(options["friction"]) is True:
+ grass.message('Calculating costs using input friction map\n')
+ friction = options["friction"]
+ else:
+ grass.message('Calculating for time costs only')
+ friction = "temporary.friction.%s" % pid
+ grass.mapcalc("${out} = if(isnull(${rast1}), null(), 0)", overwrite = grass.overwrite(), quiet = True, out = friction, rast1 = elevation)
+ if flags["k"] is True:
+ grass.message('Using Knight\'s move\n')
+ #NOTE! because "lambda" is an internal python variable, it is impossible to enter the value for key "lambda" in r.walk. It ends up with a python error.
+ grass.run_command('r.walk', quiet = True, overwrite = grass.overwrite(), flags = 'k', elevation = elevation, friction = friction, output = cost, start_points = start_points, walk_coeff = w_coefs, slope_factor = slope_factor)
+ else:
+ grass.run_command('r.walk', quiet = True, overwrite = grass.overwrite(), elevation = elevation, friction = friction, output = cost, start_points = start_points, percent_memory = '100', walk_coeff = w_coefs, slope_factor = slope_factor)
+ if bool(options["friction"]) is False:
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = friction)
+#################################################
+ if bool(options["sigma"]) is True:
+ grass.message('\n\nCreating optional slope mask\n')
+ slope = "temporary.slope.%s" % pid
+ grass.run_command('r.slope.aspect', quiet = True, overwrite = grass.overwrite(), elevation = elevation, slope = slope)
+ if ismask == 2:
+ grass.mapcalc("MASK=if(${rast1} <= ${sigma}, 1, if(${tempmask}, 1, null()))", overwrite = grass.overwrite(), quiet = True, sigma = sigma, rast1 = slope, tempmask = tempmask)
+ else:
+ grass.mapcalc("MASK=if(${rast1} <= ${sigma}, 1, null())", overwrite = grass.overwrite(), quiet = True, sigma = sigma, rast1 = slope)
+ else:
+ grass.message('No slope mask created')
+##################################################
+ if flags["l"] is True:
+ grass.message('\n\nCalculating list of possible catchment configurations\n')
+ grass.message("cost value | catchment area")
+ areadict = {}
+ out2dictnum('r.stats -Aani input=' + cost + ' fs=, nv=* nsteps=255', ',', areadict)
+ testarea = 0
+ #start the loop, and list the values
+ for key in sorted(areadict):
+ testarea = testarea + int(float(areadict[key]))
+ grass.message("%s | %s" % (int(key), testarea))
+ if flags["c"] is True:
+ if bool(options["in_cost"]) is False:
+ grass.run_command('g.rename', overwrite = grass.overwrite(), quiet = True, rast = 'temporary.cost.%s,%s_cost_surface' % (pid, buffer))
+ grass.message('Cleaning up...(keeping cost map)')
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
+ else:
+ grass.message('Cleaning up...1')
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
+ else:
+ if bool(options["in_cost"]) is False:
+ grass.message('Cleaning up...2')
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s,temporary.cost.%s' % (pid, pid))
+ else:
+ grass.message('Cleaning up...3')
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
+ if bool(options["sigma"]) is True:
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = slope)
+ if ismask == 2:
+ grass.message('Reinstating original MASK...')
+ grass.run_command('g.rename', overwrite = grass.overwrite(), quiet = "True", rast = tempmask +',MASK')
+ elif ismask == 0 and bool(options["sigma"]) is True:
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'MASK')
+ elif ismask == 1:
+ grass.message('Keeping original MASK')
+ grass.message(' DONE!')
+ return
+ else:
+ grass.message('\n\nCalculating buffer\n')
+ areadict = {}
+ out2dictnum('r.stats -Aani input=' + cost + ' fs=, nv=* nsteps=255', ',', areadict)
+ tot_area = 0
+ for key in sorted(areadict):
+ tot_area = tot_area + int(float(areadict[key]))
+ maxcost = key
+ grass.message("Maximum cost distance value %s covers an area of %s square map units\n\nCommencing to find a catchment configuration.....\n\n" % (int(maxcost), tot_area))
+ testarea = 0
+ lastarea = 0
+ lastkey = 0
+ #start the loop, and home in on the target range
+ for key in sorted(areadict):
+ testarea = testarea + int(float(areadict[key]))
+ if testarea >= area:
+ break
+ lastkey = key
+ lastarea = testarea
+ if (testarea - area) <= (area - lastarea):
+ cutoff = key
+ displayarea = testarea
+ else:
+ cutoff = lastkey
+ displayarea = lastarea
+ grass.message("Catchment configuration found! Cost cutoff %s produces a catchment of %s square map units." % (int(cutoff), displayarea))
+ ####################################################
+ grass.message('\n\nCreating output map\n')
+ temp = tempfile.NamedTemporaryFile()
+ temp.write('0 thru %s = %s\n' % (int(cutoff), mapval))
+ temp.flush()
+ grass.run_command('r.reclass', overwrite = grass.overwrite(), input = cost, output = 'cost.reclass.%s' % pid, rules = temp.name)
+ temp.close()
+ grass.mapcalc("${out}=if(isnull(${cost}), null(), ${cost})", overwrite = grass.overwrite(), quiet = True, cost = "cost.reclass.%s" % pid, out = buffer)
+ grass.message("\nThe output catchment map will be named %s" % buffer)
+ grass.run_command('r.colors', quiet = True, map = buffer, color = 'ryb')
+ if flags["c"] is True:
+ if bool(options["in_cost"]) is False:
+ grass.run_command('g.rename', overwrite = grass.overwrite(), quiet = True, rast = 'temporary.cost,%s_cost_surface' % (buffer))
+ grass.message('Cleaning up...(keeping cost map)')
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
+ else:
+ grass.message('Cleaning up...1')
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
+ else:
+ if bool(options["in_cost"]) is False:
+ grass.message('Cleaning up...2')
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s,temporary.cost.%s' % (pid, pid))
+ else:
+ grass.message('Cleaning up...3')
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = 'cost.reclass.%s' % pid)
+ if bool(options["sigma"]) is True:
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name = slope)
+ if ismask == 2:
+ grass.message('Reinstating original MASK...')
+ grass.run_command('g.rename', overwrite = grass.overwrite(), quiet = "True", rast = tempmask +',MASK')
+ elif ismask == 0 and bool(options["sigma"]) is True:
+ grass.run_command('g.remove', quiet = True, flags = 'f', type = 'raster', name ='MASK')
+ elif ismask == 1:
+ grass.message('Keeping original MASK')
+ grass.message(' DONE!')
+ return
+
+# here is where the code in "main" actually gets executed. This way of programming is neccessary for the way g.parser needs to run in GRASS 7.
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
+ exit(0)
More information about the grass-commit
mailing list