[GRASS-SVN] r62672 - grass-addons/grass7/raster/r.basin
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Nov 8 09:55:08 PST 2014
Author: neteler
Date: 2014-11-08 09:55:08 -0800 (Sat, 08 Nov 2014)
New Revision: 62672
Modified:
grass-addons/grass7/raster/r.basin/r.basin.html
grass-addons/grass7/raster/r.basin/r.basin.py
Log:
r.basin Addon: various bugfixes + example added
Modified: grass-addons/grass7/raster/r.basin/r.basin.html
===================================================================
--- grass-addons/grass7/raster/r.basin/r.basin.html 2014-11-08 16:54:08 UTC (rev 62671)
+++ grass-addons/grass7/raster/r.basin/r.basin.html 2014-11-08 17:55:08 UTC (rev 62672)
@@ -115,33 +115,47 @@
<h2>EXAMPLE</h2>
+North Carolina sample dataset example:
+
<div class="code"><pre>
-r.basin map=elevation prefix=example coordinates=636639,218816 dir=/home/yourusername/Desktop threshold=1000
+g.region rast=elevation -p
+r.basin map=elevation prefix=my_basin coord=637500.0,221750.0 dir=/tmp/bla threshold=1000
+
+# visualize some results
+d.mon wx0
+d.rast my_basin_elevation_hack
+
+d.rast my_basin_elevation_dist2out
+
+d.his i=aspect h=my_basin_elevation_dist2out
</pre></div>
<h3>Dependencies</h3>
<ul>
<li>Matplotlib</li>
+<li>r.hypso</li>
<li>r.stream.basins</li>
<li>r.stream.extract</li>
<li>r.stream.stats</li>
<li>r.stream.distance</li>
<li>r.stream.order</li>
<li>r.width.funct</li>
-<li>r.hypso</li>
</ul>
<h2>SEE ALSO</h2>
<em>
-<a href="r.stream.basins.html">r.stream.basins</a>,
-<a href="r.stream.distance.html">r.stream.distance</a>,
+<a href="r.hypso.html">r.hypso</a>,
+<a href="http://grass.osgeo.org/grass70/manuals/addons/r.stream.channel.html">r.stream.channel</a> (Addon),
+<a href="http://grass.osgeo.org/grass70/manuals/addons/r.stream.distance.html">r.stream.distance</a> (Addon),
<a href="r.stream.extract.html">r.stream.extract</a>,
-<a href="r.stream.order.html">r.stream.order</a>,
-<a href="r.stream.stats.html">r.stream.stats</a>,
-<a href="r.width.funct.html">r.width.funct</a>,
-<a href="r.hypso.html">r.hypso</a>,
-<a href="r.watershed.html">r.watershed</a>
+<a href="http://grass.osgeo.org/grass70/manuals/addons/r.stream.order.html">r.stream.order</a> (Addon),
+<a href="http://grass.osgeo.org/grass70/manuals/addons/r.stream.segment.html">r.stream.segment</a> (Addon),
+<a href="http://grass.osgeo.org/grass70/manuals/addons/r.stream.slope.html">r.stream.slope</a> (Addon),
+<a href="http://grass.osgeo.org/grass70/manuals/addons/r.stream.snap.html">r.stream.snap</a> (Addon),
+<a href="http://grass.osgeo.org/grass70/manuals/addons/r.stream.stats.html">r.stream.stats</a> (Addon),
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.width.funct.html">r.width.funct</a>
</em>
<h2>CITE AS</h2>
Modified: grass-addons/grass7/raster/r.basin/r.basin.py
===================================================================
--- grass-addons/grass7/raster/r.basin/r.basin.py 2014-11-08 16:54:08 UTC (rev 62671)
+++ grass-addons/grass7/raster/r.basin/r.basin.py 2014-11-08 17:55:08 UTC (rev 62672)
@@ -78,13 +78,26 @@
grass.message( "You must be in GRASS GIS to run this program." )
sys.exit(1)
+# check requirements
+def check_progs():
+ for prog in ('r.hypso', 'r.stream.basins', 'r.stream.distance', 'r.stream.extract',
+ 'r.stream.order','r.stream.snap','r.stream.stats', 'r.width.funct'):
+ if not grass.find_program(prog, '--help'):
+ grass.fatal(_("'%s' required. Please install '%s' first using 'g.extension %s'") % (prog, prog, prog))
+
def main():
+ # check dependencies
+ check_progs()
+
r_elevation = options['map'].split('@')[0]
mapname = options['map'].replace("@"," ")
mapname = mapname.split()
mapname[0] = mapname[0].replace(".","_")
coordinates = options['coordinates']
directory = options['dir']
+ # Check if directory exists
+ if not os.path.isdir(directory):
+ os.makedirs(directory)
autothreshold = flags['a']
nomap = flags['c']
prefix = options['prefix']+'_'+mapname[0]
@@ -122,15 +135,14 @@
# Save current region
- grass.read_command('g.region', flags = 'p', save = 'original', overwrite = True)
+ grass.read_command('g.region', flags = 'p', save = 'original')
# Watershed SFD
grass.run_command('r.watershed', elevation = r_elevation,
accumulation = r_accumulation,
drainage = r_drainage,
convergence = 5,
- flags = 'am',
- overwrite = True)
+ flags = 'am')
# Managing flag
if autothreshold :
@@ -147,26 +159,18 @@
d8cut = 1000000000,
mexp = 0,
stream_rast = r_stream_e,
- direction = r_drainage_e,
- overwrite = True)
+ direction = r_drainage_e)
- try:
-
- # Delineation of basin
- ## check if r.stream.basins addon is installed
- if not grass.find_program('r.stream.basins', '--help'):
- grass.fatal(_("The 'r.stream.basins' module was not found, install it first:") +
- "\n" +
- "g.extension r.stream.basins")
-
+ try:
+ # Delineation of basin
# Create outlet
grass.write_command('v.in.ascii', output = v_outlet,
input = "-",
sep = ",",
- stdin = "%s,9999" % (coordinates),
- overwrite = True)
+ stdin = "%s,9999" % (coordinates))
# Snap outlet to stream network
+ # TODO: does snap depend on the raster resolution? hardcoded 30 below
grass.run_command('r.stream.snap', input = v_outlet,
output = v_outlet_snap,
stream_rast = r_stream_e,
@@ -177,14 +181,12 @@
use = 'cat',
type = 'point',
layer = 1,
- value = 1,
- overwrite = True)
+ value = 1)
grass.run_command('r.stream.basins', direction = r_drainage_e,
basins = r_basin,
- points = v_outlet_snap,
- overwrite = True)
+ points = v_outlet_snap)
grass.message( "Delineation of basin done" )
@@ -214,8 +216,7 @@
grass.mapcalc("$r_elevation_crop = $r_elevation * $r_mask",
r_mask = r_mask,
r_elevation = r_elevation,
- r_elevation_crop = 'r_elevation_crop',
- overwrite = True)
+ r_elevation_crop = 'r_elevation_crop')
grass.mapcalc("tmp = $r_drainage_e * $r_mask",
r_mask = r_mask,
@@ -235,27 +236,23 @@
grass.run_command('g.rename', rast = ('tmp',r_stream_e))
grass.run_command('r.thin', input = r_stream_e,
- output = r_stream_e+'_thin',
- overwrite = True)
+ output = r_stream_e+'_thin')
grass.run_command('r.to.vect', input = r_stream_e+'_thin',
output = v_network,
- type = 'line',
- overwrite = True)
+ type = 'line')
# Creation of slope and aspect maps
grass.run_command('r.slope.aspect', elevation = 'r_elevation_crop',
slope = r_slope,
- aspect = r_aspect,
- overwrite = True)
+ aspect = r_aspect)
# Basin mask (vector)
# Raster to vector
grass.run_command('r.to.vect', input = r_basin,
output = v_basin,
type = 'area',
- flags = 'sv',
- overwrite = True)
+ flags = 'sv')
# Add two columns to the table: area and perimeter
grass.run_command('v.db.addcolumn', map = v_basin,
@@ -271,8 +268,7 @@
qlayer = 1,
option = 'perimeter',
units = 'kilometers',
- columns = 'perimeter',
- overwrite = True)
+ columns = 'perimeter')
# Read perimeter
tmp = grass.read_command('v.to.db', map = v_basin,
@@ -292,8 +288,7 @@
qlayer = 1,
option = 'area',
units = 'kilometers',
- columns = 'area',
- overwrite = True)
+ columns = 'area')
# Read area
tmp = grass.read_command('v.to.db', map = v_basin,
@@ -314,26 +309,19 @@
strahler = r_strahler,
shreve = r_shreve,
horton = r_horton,
- hack = r_hack,
- overwrite = True)
+ hack = r_hack)
# Distance to outlet
grass.run_command('r.stream.distance', stream_rast = r_outlet,
direction = r_drainage_e,
flags = 'o',
- distance = r_distance,
- overwrite = True)
+ distance = r_distance)
# hypsographic curve
grass.message( "------------------------------" )
- #### check if r.hypso addon is installed
- if not grass.find_program('r.hypso', '--help'):
- grass.fatal(_("The 'r.hypso' module was not found, install it first:") +
- "\n" +
- "g.extension r.hypso")
grass.run_command('r.hypso', map = 'r_elevation_crop',
image = os.path.join(directory,prefix), flags = 'ab')
@@ -343,11 +331,6 @@
grass.message( "------------------------------" )
- #### check if r.width.funct addon is installed
- if not grass.find_program('r.width.funct', '--help'):
- grass.fatal(_("The 'r.width.funct' module was not found, install it first:") +
- "\n" +
- "g.extension r.width.funct")
grass.run_command('r.width.funct', map = r_distance,
image = os.path.join(directory,prefix))
@@ -358,16 +341,14 @@
grass.run_command("r.stream.distance", stream_rast = r_stream_e,
direction = r_drainage,
elevation = 'r_elevation_crop',
- distance = r_hillslope_distance,
- overwrite = True)
+ distance = r_hillslope_distance)
# Mean elevation
grass.run_command("r.stats.zonal", base = r_basin,
cover = "r_elevation_crop",
method = "average",
- output = r_height_average,
- overwrite = True)
+ output = r_height_average)
grass.message("r.stats.zonal done")
@@ -387,8 +368,7 @@
# Centroid and mean slope
baricenter_slope_baricenter = grass.read_command("r.volume", input = r_slope,
- clump = r_basin,
- overwrite = True)
+ clump = r_basin)
grass.message("r.volume done")
@@ -441,13 +421,11 @@
r_mainchannel = r_mainchannel)
grass.run_command("r.thin", input = r_mainchannel,
- output = r_mainchannel+'_thin',
- overwrite = True)
+ output = r_mainchannel+'_thin')
grass.run_command('r.to.vect', input = r_mainchannel+'_thin',
output = v_mainchannel,
type = 'line',
- verbose = True,
- overwrite = True)
+ verbose = True)
# Get coordinates of the outlet (belonging to stream network)
@@ -473,7 +451,7 @@
param_mainchannel = grass.read_command('v.what', map = v_mainchannel,
coordinates = '%s,%s' % (east_o,north_o),
- distance = 5 )
+ distance = 5)
tmp = param_mainchannel.split('\n')[7]
mainchannel = float(tmp.split()[1]) / 1000 # km
@@ -483,14 +461,12 @@
r_shreve = r_shreve,
r_mainchannel = r_mainchannel)
grass.run_command('r.thin', input = r_mainchannel_dim,
- output = r_mainchannel_dim + '_thin',
- overwrite = True)
+ output = r_mainchannel_dim + '_thin')
grass.run_command('r.to.vect', input = r_mainchannel_dim + '_thin',
output = v_mainchannel_dim,
type = 'line',
flags = 'v',
- verbose = True,
- overwrite = True)
+ verbose = True)
try:
D_topo1 = grass.read_command('v.info', map = v_mainchannel_dim,
layer = 1,
@@ -505,8 +481,7 @@
grass.run_command('v.to.points',
input = v_mainchannel_dim,
output = v_mainchannel_dim+'_point',
- type = 'line',
- overwrite = True)
+ type = 'line')
vertex = grass.read_command('v.out.ascii', verbose = True,
input = v_mainchannel_dim+'_point').strip().split('\n')
nodi = zeros((len(vertex),4),float)
@@ -559,24 +534,21 @@
grass.run_command("r.stats.zonal", cover = r_stream_e,
base = r_mask,
method = "average",
- output = r_average_hillslope,
- overwrite = True)
- mean_hillslope_length = float(grass.read_command('r.info', flags = 'r',
- map = r_average_hillslope).split('\n')[0].split('=')[1])
+ output = r_average_hillslope)
+ mean_hillslope_length = float(grass.read_command('r.info', flags = 'r',
+ map = r_average_hillslope).split('\n')[0].split('=')[1])
- # Magnitudo
+ # Magnitude
grass.mapcalc("$r_ord_1 = if($r_strahler==1,1,null())",
r_ord_1 = r_ord_1,
r_strahler = r_strahler)
grass.run_command('r.thin', input = r_ord_1,
output = r_ord_1+'_thin',
- iterations = 200,
- overwrite = True)
+ iterations = 200)
grass.run_command('r.to.vect', input = r_ord_1+'_thin',
output = v_ord_1,
type = 'line',
- flags = 'v',
- overwrite = True)
+ flags = 'v')
magnitudo = float(grass.read_command('v.info', map = v_ord_1,
layer = 1,
flags = 't').split('\n')[2].split('=')[1])
@@ -789,7 +761,7 @@
+ [Slope_ratio])
-# Import table "summary", attaches it to "outlet_snap", then drops it
+ # Import table "summary", attaches it to "outlet_snap", then drops it
grass.run_command("db.in.ogr", dsn = csvfileT,
output = "summary")
@@ -849,7 +821,7 @@
grass.message( "------------------------------" )
grass.message( "\n" )
grass.message( "An ERROR occurred running r.basin" )
- grass.message( "Please try with another pairs of outlet coordinates" )
+ grass.message( "Please check for error messages above or try with another pairs of outlet coordinates" )
# Set region to original
More information about the grass-commit
mailing list