[GRASS-SVN] r49993 - grass-addons/grass6/raster/r.basin
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Dec 30 14:51:18 EST 2011
Author: madi
Date: 2011-12-30 11:51:17 -0800 (Fri, 30 Dec 2011)
New Revision: 49993
Modified:
grass-addons/grass6/raster/r.basin/description.html
grass-addons/grass6/raster/r.basin/r.basin.py
Log:
general improvements, controls added
Modified: grass-addons/grass6/raster/r.basin/description.html
===================================================================
--- grass-addons/grass6/raster/r.basin/description.html 2011-12-30 17:49:00 UTC (rev 49992)
+++ grass-addons/grass6/raster/r.basin/description.html 2011-12-30 19:51:17 UTC (rev 49993)
@@ -3,6 +3,11 @@
starting from the digital elevation model and the coordinates of
the basin's closing section.
+<h3>Important:</h3>
+<p>The closing section's coordinates must belong to the river network generated by
+<a href="r.stream.extract.html">r.stream.extract</a>. It is suggested to run it
+before.</p>
+
<h3>Morphometric parameters of basin</h3>
<p>The main parameters are:
<ul>
Modified: grass-addons/grass6/raster/r.basin/r.basin.py
===================================================================
--- grass-addons/grass6/raster/r.basin/r.basin.py 2011-12-30 17:49:00 UTC (rev 49992)
+++ grass-addons/grass6/raster/r.basin/r.basin.py 2011-12-30 19:51:17 UTC (rev 49993)
@@ -5,10 +5,12 @@
# MODULE: r.basin.py
# AUTHOR(S): Margherita Di Leo, Massimo Di Stefano
# PURPOSE: Morphometric characterization of river basins
-# COPYRIGHT: (C) 2010 by the GRASS Development Team
+# COPYRIGHT: (C) 2010 by Margherita Di Leo & Massimo Di Stefano
+# dileomargherita at gmail.com
#
# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
+# License (>=v3.0) and comes with ABSOLUTELY NO WARRANTY.
+# See the file COPYING that comes with GRASS
# for details.
#
#############################################################################
@@ -55,7 +57,7 @@
#%end
#%flag
#% key: a
-#% description: Use default Threshold (1km^2)
+#% description: Use default threshold (1km^2)
#%END
#%flag
#% key: c
@@ -93,7 +95,7 @@
r_shreve = prefix+'_shreve'
r_horton = prefix+'_horton'
r_hack = prefix+'_hack'
- r_distance = prefix+'_distance'
+ r_distance = prefix+'_dist2out'
r_hillslope_distance = prefix+'_hillslope_distance'
r_height_average = prefix+'_height_average'
r_aspect_mod = prefix+'_aspect_mod'
@@ -106,16 +108,26 @@
r_average_hillslope = prefix+'_average_hillslope'
r_mainchannel_dim = prefix+'_mainchannel_dim'
r_outlet = prefix+'_r_outlet'
- v_outlet = prefix+'_v_outlet'
+ v_outlet = prefix+'_outlet'
v_basin = prefix+'_basin'
v_centroid1 = prefix+'_centroid1'
v_mainchannel = prefix+'_mainchannel'
v_mainchannel_dim = prefix+'_mainchannel_dim'
v_network = prefix+'_network'
- v_stream_e = prefix+'_stream_e'
v_ord_1 = prefix+'_ord_1'
grass.run_command('g.remove' , rast = 'MASK')
+
+ # Save current region
+ 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 = 'a',
+ overwrite = True)
+
# Managing flag
if autothreshold :
info_region = grass.read_command('g.region', flags = 'p', rast = '%s' % (r_elevation))
@@ -126,300 +138,446 @@
else :
th = options['threshold']
- # Watershed SFD
- grass.run_command('r.watershed', elevation = r_elevation , threshold = th , accumulation = r_accumulation , drainage = r_drainage , stream = r_stream+'_nothin' , convergence = 5, flags = 'a')
-
# Stream extraction
- grass.run_command('r.stream.extract', elevation = r_elevation, accumulation = r_accumulation, threshold = th, d8cut = 'infinity', mexp = 0, stream_rast = r_stream_e, stream_vect = v_stream_e, direction = r_drainage_e, overwrite = True)
+ grass.run_command('r.stream.extract', elevation = r_elevation,
+ accumulation = r_accumulation,
+ threshold = th,
+ d8cut = 'infinity',
+ mexp = 0,
+ stream_rast = r_stream_e,
+ direction = r_drainage_e,
+ overwrite = True)
+
+ try:
- # Delineation of basin
- grass.run_command('r.stream.basins', dir = r_drainage , basins = r_basin , coors = '%s,%s' % (east , north) )
- grass.message( "Delineation of basin done" )
+ # Delineation of basin
+ grass.run_command('r.stream.basins', dir = r_drainage,
+ basins = r_basin,
+ coors = '%s,%s' % (east , north),
+ overwrite = True)
+ grass.message( "Delineation of basin done" )
+
+ # Mask and cropping
+ elevation_name = r_elevation = r_elevation.split('@')[0]
+ grass.mapcalc("$r_mask = $r_basin / $r_basin",
+ r_mask = r_mask,
+ r_basin = r_basin)
+ grass.mapcalc("$r_accumulation = $r_accumulation / $r_mask",
+ r_accumulation = r_accumulation,
+ r_mask = r_mask)
+ grass.mapcalc("$r_drainage = $r_drainage / $r_mask",
+ r_drainage = r_drainage,
+ r_mask = r_mask)
+ grass.mapcalc("$r_elevation_crop = $r_elevation * $r_mask",
+ r_mask = r_mask,
+ r_elevation = r_elevation,
+ r_elevation_crop = 'r_elevation_crop')
+ grass.mapcalc("$r_drainage_e = $r_drainage_e * $r_mask",
+ r_mask = r_mask,
+ r_drainage_e = r_drainage_e)
+ grass.mapcalc("$r_stream_e = $r_stream_e * $r_mask",
+ r_mask = r_mask,
+ r_stream_e = r_stream_e)
+ grass.run_command('r.thin', input = r_stream_e,
+ output = r_stream_e+'_thin',
+ overwrite = True)
+ grass.run_command('r.to.vect', input = r_stream_e+'_thin',
+ output = v_network,
+ feature = 'line',
+ overwrite = True)
- # Backup and mask
- elevation_name = r_elevation = r_elevation.split('@')[0]
- grass.run_command('g.copy', rast = r_elevation+','+elevation_name+'_crop')
- grass.run_command('r.mapcalculator', amap = r_basin , formula = '%s/%s' % (r_basin, r_basin) , outfile = r_mask, overwrite = True)
- grass.run_command('r.mapcalculator', amap = r_accumulation , formula = '%s*%s' % (r_accumulation, r_mask) , outfile = r_accumulation, overwrite = True)
- grass.run_command('r.mapcalculator', amap = r_drainage , formula = '%s*%s' % (r_drainage, r_mask) , outfile = r_drainage, overwrite = True)
- grass.run_command('r.mapcalculator', amap = r_stream+'_nothin' , formula = '%s*%s' % (r_stream+'_nothin', r_mask) , outfile = r_stream+'_nothin', overwrite = True)
- grass.run_command('r.mapcalculator', amap = r_elevation , formula = '%s*%s' % (r_elevation, r_mask) , outfile = r_elevation+'_crop', overwrite = True)
- grass.run_command('r.mapcalculator', amap = r_drainage_e , formula = '%s*%s' % (r_drainage_e, r_mask) , outfile = r_drainage_e, overwrite = True)
- grass.run_command('r.mapcalculator', amap = r_stream_e , formula = '%s*%s' % (r_stream_e, r_mask) , outfile = r_stream_e, overwrite = True)
- grass.run_command('r.thin', input = r_stream_e , output = r_stream_e+'_thin' )
- grass.run_command('r.to.vect', input = r_stream_e+'_thin' , output = v_network , feature = 'line' )
-
- # Creation of slope and aspect maps
- grass.run_command('r.slope.aspect', elevation = r_elevation+'_crop' , slope = r_slope , aspect = r_aspect )
+ # Creation of slope and aspect maps
+ grass.run_command('r.slope.aspect', elevation = 'r_elevation_crop',
+ slope = r_slope,
+ aspect = r_aspect,
+ overwrite = True)
- # Basin mask (vector)
- grass.run_command('r.to.vect', input = r_basin , output = v_basin , feature = 'area' )
- grass.run_command('v.to.db', map = v_basin , type = 'line,boundary' , layer = 1, qlayer = 1, option = 'perimeter' , units = 'kilometers', columns = 'cat', qcolumn = 'cat')
+ # Basin mask (vector)
+ grass.run_command('r.to.vect', input = r_basin,
+ output = v_basin,
+ feature = 'area',
+ overwrite = True)
+ grass.run_command('v.to.db', map = v_basin,
+ type = 'line,boundary',
+ layer = 1,
+ qlayer = 1,
+ option = 'perimeter',
+ units = 'kilometers',
+ columns = 'cat',
+ qcolumn = 'cat',
+ overwrite = True)
- # Creation of order maps: strahler, horton, hack, shreeve
- grass.message( "Creating %s" % r_hack )
- grass.run_command('r.stream.order', stream = r_stream_e , dir = r_drainage_e , strahler = r_strahler , shreve = r_shreve , horton = r_horton , hack = r_hack )
+ # Creation of order maps: strahler, horton, hack, shreeve
+ grass.message( "Creating %s" % r_hack )
+ grass.run_command('r.stream.order', stream = r_stream_e,
+ dir = r_drainage_e,
+ strahler = r_strahler,
+ shreve = r_shreve,
+ horton = r_horton,
+ hack = r_hack,
+ overwrite = True)
- # Distance to outlet
- grass.write_command('v.in.ascii', output = v_outlet, stdin = "%s|%s|9999" % (east, north) )
- grass.run_command('v.to.rast', input = v_outlet, output = r_outlet , use = 'cat' , type = 'point,line,area' , layer = 1 , value = 1 , rows = 4096)
- grass.run_command('r.stream.distance', stream = r_outlet , dir = r_drainage , flags = 'o' , distance = r_distance )
+ # Distance to outlet
+ grass.write_command('v.in.ascii', output = v_outlet,
+ stdin = "%s|%s|9999" % (east, north),
+ overwrite = True)
+ grass.run_command('v.to.rast', input = v_outlet,
+ output = r_outlet,
+ use = 'cat',
+ type = 'point',
+ layer = 1,
+ value = 1,
+ rows = 4096,
+ overwrite = True)
+ grass.run_command('r.stream.distance', stream = r_outlet,
+ dir = r_drainage,
+ flags = 'o',
+ distance = r_distance,
+ overwrite = True)
- # Width Function
- grass.run_command('r.wf.py', map = r_distance, image = prefix)
+ # Width Function
+ grass.run_command('r.wf.py', map = r_distance,
+ image = prefix)
- # Ipsographic curve
- grass.run_command('r.ipso.py', map = r_elevation+'_crop', image = prefix, flags = 'ab')
+ # Ipsographic curve
+ grass.run_command('r.ipso.py', map = 'r_elevation_crop',
+ image = prefix,
+ flags = 'ab')
- # Creation of map of hillslope distance to river network
- grass.run_command('r.stream.distance', stream = r_stream+'_nothin' , dir = r_drainage , dem = r_elevation+'_crop' , distance = r_hillslope_distance )
+ # Creation of map of hillslope distance to river network
+ grass.run_command('r.stream.distance', stream = r_stream_e,
+ dir = r_drainage,
+ dem = 'r_elevation_crop' ,
+ distance = r_hillslope_distance,
+ overwrite = True)
- # Mean elevation
- grass.run_command('r.average' , base = r_basin , cover = r_elevation+'_crop' , output = r_height_average)
- mean_elev = float(grass.read_command('r.info' , flags = 'r' , map = r_height_average).split('\n')[0].split('=')[1])
+ # Mean elevation
+ grass.run_command('r.average' , base = r_basin,
+ cover = 'r_elevation_crop',
+ output = r_height_average,
+ overwrite = True)
+ mean_elev = float(grass.read_command('r.info', flags = 'r',
+ map = r_height_average).split('\n')[0].split('=')[1])
- # In Grass the aspect categories represent the number degrees of east and they increase counterclockwise:
- # 90deg is North, 180 is West, 270 is South 360 is East. The aspect value 0 is used to indicate undefined aspect in flat areas with slope=0.
- # We calculate the number of degree from north, increasing counterclockwise.
- grass.run_command('r.mapcalculator', amap = r_aspect , formula = 'if(%s>90,450-%s,90-%s)' % (r_aspect, r_aspect, r_aspect) , outfile = r_aspect_mod)
+ # In Grass, aspect categories represent the number degrees of east and they increase
+ # counterclockwise: 90deg is North, 180 is West, 270 is South 360 is East.
+ # The aspect value 0 is used to indicate undefined aspect in flat areas with slope=0.
+ # We calculate the number of degree from north, increasing counterclockwise.
+ grass.mapcalc("$r_aspect_mod = if($r_aspect > 90, 450 - $r_aspect, 90 - $r_aspect)",
+ r_aspect = r_aspect,
+ r_aspect_mod = r_aspect_mod)
- # Centroid and mean slope
- baricenter_slope_baricenter = grass.read_command('r.volume', data = r_slope ,clump = r_basin , centroids = v_centroid1)
- baricenter_slope_baricenter = baricenter_slope_baricenter.split()
- mean_slope = baricenter_slope_baricenter[28]
- baricenter_aspect_baricenter = grass.read_command('r.volume', data = r_aspect ,clump = r_basin )
- baricenter_aspect_baricenter = baricenter_aspect_baricenter.split()
- basin_average_aspect = baricenter_aspect_baricenter[28]
+ # Centroid and mean slope
+ baricenter_slope_baricenter = grass.read_command('r.volume', data = r_slope,
+ clump = r_basin,
+ centroids = v_centroid1,
+ overwrite = True)
+ baricenter_slope_baricenter = baricenter_slope_baricenter.split()
+ mean_slope = baricenter_slope_baricenter[28]
+ baricenter_aspect_baricenter = grass.read_command('r.volume', data = r_aspect,
+ clump = r_basin )
+ baricenter_aspect_baricenter = baricenter_aspect_baricenter.split()
+ basin_average_aspect = baricenter_aspect_baricenter[28]
- # Rectangle containing basin
- basin_east = baricenter_slope_baricenter[31]
- basin_north = baricenter_slope_baricenter[32]
- info_region_basin = grass.read_command('g.region', vect = options['prefix']+'_'+mapname[0]+'_basin' ,flags = 'm')
- dict_region_basin = dict(x.split('=', 1) for x in info_region_basin.split('\n') if '=' in x)
- basin_resolution = float(dict_region_basin['nsres'])
- x_massimo = float(dict_region_basin['n']) + (basin_resolution * 10)
- x_minimo = float(dict_region_basin['w']) - (basin_resolution * 10)
- y_massimo = float(dict_region_basin['e']) + (basin_resolution * 10)
- y_minimo = float(dict_region_basin['s']) - (basin_resolution * 10)
- nw = dict_region_basin['w'] , dict_region_basin['n']
- se = dict_region_basin['e'] , dict_region_basin['s']
+ # Rectangle containing basin
+ basin_east = baricenter_slope_baricenter[31]
+ basin_north = baricenter_slope_baricenter[32]
+ info_region_basin = grass.read_command('g.region',
+ vect = options['prefix']+'_'+mapname[0]+'_basin',
+ flags = 'm')
+ dict_region_basin = dict(x.split('=', 1) for x in info_region_basin.split('\n') if '=' in x)
+ basin_resolution = float(dict_region_basin['nsres'])
+ x_massimo = float(dict_region_basin['n']) + (basin_resolution * 10)
+ x_minimo = float(dict_region_basin['w']) - (basin_resolution * 10)
+ y_massimo = float(dict_region_basin['e']) + (basin_resolution * 10)
+ y_minimo = float(dict_region_basin['s']) - (basin_resolution * 10)
+ nw = dict_region_basin['w'], dict_region_basin['n']
+ se = dict_region_basin['e'], dict_region_basin['s']
- # Area and perimeter of basin
- report_basin = grass.read_command('v.report', map = v_basin , layer = 1 , option = 'area' , units = 'kilometers')
- report_basin = report_basin.replace('\n',' ')
- report_basin = report_basin.split()
- area_basin = float(report_basin[1].split('|')[3])
- report_basin = grass.read_command('v.report', map = v_basin , layer = 1 , option = 'length' , units = 'meters')
- perimeter_basin = float(report_basin.split('\n')[1].split('|')[0])
+ # Area and perimeter of basin
+ report_basin = grass.read_command('v.report', map = v_basin,
+ layer = 1,
+ option = 'area',
+ units = 'kilometers')
+ report_basin = report_basin.replace('\n',' ')
+ report_basin = report_basin.split()
+ area_basin = float(report_basin[1].split('|')[3])
+ report_basin = grass.read_command('v.report', map = v_basin,
+ layer = 1,
+ option = 'length',
+ units = 'meters')
+ perimeter_basin = float(report_basin.split('\n')[1].split('|')[0])
- # Directing vector
- delta_x = abs(float(basin_east) - east)
- delta_y = abs(float(basin_north) - north)
- L_orienting_vect = math.sqrt((delta_x**2)+(delta_y**2)) / 1000
+ # Directing vector
+ delta_x = abs(float(basin_east) - east)
+ delta_y = abs(float(basin_north) - north)
+ L_orienting_vect = math.sqrt((delta_x**2)+(delta_y**2)) / 1000
- # Prevalent orientation
- prevalent_orientation = math.atan(delta_y/delta_x)
+ # Prevalent orientation
+ prevalent_orientation = math.atan(delta_y/delta_x)
- # Compactness coefficient
- C_comp = perimeter_basin / ( 2 * math.sqrt( area_basin / math.pi))
+ # Compactness coefficient
+ C_comp = perimeter_basin / ( 2 * math.sqrt( area_basin / math.pi))
- # Circularity ratio
- R_c = ( 4 * math.pi * area_basin ) / (perimeter_basin **2)
+ # Circularity ratio
+ R_c = ( 4 * math.pi * area_basin ) / (perimeter_basin **2)
- # Mainchannel
- grass.run_command('r.mapcalculator', amap = r_hack , formula = 'if(%s==1,1,null())' % (r_hack) , outfile = r_mainchannel)
- grass.run_command('r.thin', input = r_mainchannel , output = r_mainchannel+'_thin' )
- grass.run_command('r.to.vect', input = r_mainchannel+'_thin' , output = v_mainchannel , feature = 'line' , flags = 'v')
- param_mainchannel = grass.read_command('v.what', map = v_mainchannel , east_north = '%s,%s' % (east,north) )
- dict_mainchannel = dict(x.split(':', 1) for x in param_mainchannel.split('\n')if ':' in x)
- mainchannel = float(dict_mainchannel['Length']) / 1000
+ # Mainchannel
+ grass.mapcalc("$r_mainchannel = if($r_hack==1,1,null())",
+ r_hack = r_hack,
+ r_mainchannel = r_mainchannel)
+ grass.run_command('r.thin', input = r_mainchannel,
+ output = r_mainchannel+'_thin',
+ overwrite = True)
+ grass.run_command('r.to.vect', input = r_mainchannel+'_thin',
+ output = v_mainchannel,
+ feature = 'line',
+ flags = 'v',
+ overwrite = True)
+ param_mainchannel = grass.read_command('v.what', map = v_mainchannel,
+ east_north = '%s,%s' % (east,north) )
+ dict_mainchannel = dict(x.split(':', 1) for x in param_mainchannel.split('\n')if ':' in x)
+ mainchannel = float(dict_mainchannel['Length']) / 1000
- # Topological Diameter
- grass.run_command('r.mapcalculator', amap = r_shreve , formula = '-(%s - %s)' % (r_mainchannel,r_shreve) , outfile = r_mainchannel_dim)
- grass.run_command('r.thin', input = r_mainchannel_dim , output = r_mainchannel_dim+'_thin' )
- grass.run_command('r.to.vect', input = r_mainchannel_dim+'_thin' , output = v_mainchannel_dim , feature = 'line' , flags = 'v')
- try:
- D_topo = float(grass.read_command('v.info', map = v_mainchannel_dim, layer = 1, flags = 't').split('\n')[2].split('=')[1])
- except:
- D_topo = 1
- grass.message( "Topological Diameter = WARNING" )
+ # Topological Diameter
+ #grass.run_command('r.mapcalculator', amap = r_shreve , formula = '-(%s - %s)' % (r_mainchannel, r_shreve) , outfile = r_mainchannel_dim)
+ grass.mapcalc("$r_mainchannel_dim = -($r_mainchannel - $r_shreve) + 1",
+ r_mainchannel_dim = r_mainchannel_dim,
+ r_shreve = r_shreve,
+ r_mainchannel = r_mainchannel)
+ grass.run_command('r.thin', input = r_mainchannel_dim,
+ output = r_mainchannel_dim + '_thin',
+ overwrite = True)
+ grass.run_command('r.to.vect', input = r_mainchannel_dim + '_thin',
+ output = v_mainchannel_dim,
+ feature = 'line',
+ flags = 'v',
+ overwrite = True)
+ try:
+ D_topo1 = grass.read_command('v.info', map = v_mainchannel_dim,
+ layer = 1,
+ flags = 't')
+ D_topo = float(D_topo1.split('\n')[2].split('=')[1])
+ except:
+ D_topo = 1
+ grass.message( "Topological Diameter = WARNING" )
- # Mean slope of mainchannel
- grass.run_command('v.to.points', flags='n', input = v_mainchannel_dim, output = v_mainchannel_dim+'_point', 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)
- pendenze = []
- for i in range(len(vertex)):
- x, y = float(vertex[i].split('|')[0]) , float(vertex[i].split('|')[1])
- vertice = grass.read_command('r.what', verbose = True , input = r_elevation+'_crop', east_north='%s,%s' % (x,y)).replace('\n','').replace('||','|').split('|')
- nodi[i,0],nodi[i,1], nodi[i,2] = float(vertice[0]), float(vertice[1]), float(vertice[2])
- for i in range(0,len(vertex)-1,2):
- dist = math.sqrt(math.fabs((nodi[i,0] - nodi[i+1,0]))**2 + math.fabs((nodi[i,1] - nodi[i+1,1]))**2)
- deltaz = math.fabs(nodi[i,2]-nodi[i+1,2])
- # Control to prevent float division by zero (dist=0)
- try:
- pendenza = deltaz / dist
- pendenze.append(pendenza)
- mainchannel_slope = sum(pendenze) / len(pendenze) * 100
- except :
- grass.message( "Mean slope of mainchannel = WARNING" )
+ # Mean slope of mainchannel
+ grass.run_command('v.to.points', flags='n',
+ input = v_mainchannel_dim,
+ output = v_mainchannel_dim+'_point',
+ type = 'line',
+ overwrite = True)
+ vertex = grass.read_command('v.out.ascii', verbose = True,
+ input = v_mainchannel_dim+'_point').strip().split('\n')
+ nodi = zeros((len(vertex),4),float)
+ pendenze = []
+ for i in range(len(vertex)):
+ x, y = float(vertex[i].split('|')[0]) , float(vertex[i].split('|')[1])
+ vertice1 = grass.read_command('r.what', verbose = True,
+ input = 'r_elevation_crop',
+ east_north='%s,%s' % (x,y))
+ vertice = vertice1.replace('\n','').replace('||','|').split('|')
+ nodi[i,0],nodi[i,1], nodi[i,2] = float(vertice[0]), float(vertice[1]), float(vertice[2])
+ for i in range(0,len(vertex)-1,2):
+ dist = math.sqrt(math.fabs((nodi[i,0] - nodi[i+1,0]))**2 + math.fabs((nodi[i,1] - nodi[i+1,1]))**2)
+ deltaz = math.fabs(nodi[i,2] - nodi[i+1,2])
+ # Control to prevent float division by zero (dist=0)
+ try:
+ pendenza = deltaz / dist
+ pendenze.append(pendenza)
+ mainchannel_slope = sum(pendenze) / len(pendenze) * 100
+ except :
+ pass
+ #grass.message( "Mean slope of mainchannel = WARNING" )
- # Elongation Ratio
- R_al = (2 * math.sqrt( area_basin / math.pi) ) / mainchannel
+ # Elongation Ratio
+ R_al = (2 * math.sqrt( area_basin / math.pi) ) / mainchannel
- # Shape factor
- S_f = area_basin / mainchannel
+ # Shape factor
+ S_f = area_basin / mainchannel
- # Characteristic altitudes
- height_basin_average = grass.read_command('r.what' , input = r_height_average , cache = 500 , east_north = '%s,%s' % (east , north ))
- height_basin_average = height_basin_average.replace('\n','')
- height_basin_average = float(height_basin_average.split('|')[-1])
- minmax_height_basin = grass.read_command('r.info' , flags = 'r' , map = r_elevation+'_crop')
- minmax_height_basin = minmax_height_basin.strip().split('\n')
- min_height_basin , max_height_basin = float(minmax_height_basin[0].split('=')[-1]) , float(minmax_height_basin[1].split('=')[-1])
- H1 = max_height_basin
- H2 = min_height_basin
- HM = H1 - H2
+ # Characteristic altitudes
+ height_basin_average = grass.read_command('r.what', input = r_height_average ,
+ cache = 500 ,
+ east_north = '%s,%s' % (east , north ))
+ height_basin_average = height_basin_average.replace('\n','')
+ height_basin_average = float(height_basin_average.split('|')[-1])
+ minmax_height_basin = grass.read_command('r.info', flags = 'r',
+ map = 'r_elevation_crop')
+ minmax_height_basin = minmax_height_basin.strip().split('\n')
+ min_height_basin, max_height_basin = float(minmax_height_basin[0].split('=')[-1]), float(minmax_height_basin[1].split('=')[-1])
+ H1 = max_height_basin
+ H2 = min_height_basin
+ HM = H1 - H2
- # Concentration time (Giandotti, 1934)
- t_c = ((4 * math.sqrt(area_basin)) + (1.5 * mainchannel)) / (0.8 * math.sqrt(HM))
+ # Concentration time (Giandotti, 1934)
+ t_c = ((4 * math.sqrt(area_basin)) + (1.5 * mainchannel)) / (0.8 * math.sqrt(HM))
- # Mean hillslope length
- grass.run_command('r.average', cover = r_stream_e, base = r_mask, output = r_average_hillslope)
- mean_hillslope_length = float(grass.read_command('r.info' , flags = 'r' , map = r_average_hillslope).split('\n')[0].split('=')[1])
+ # Mean hillslope length
+ grass.run_command('r.average', cover = r_stream_e,
+ base = r_mask,
+ 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])
- # Magnitudo
- grass.run_command('r.mapcalculator', amap = r_strahler , formula = 'if(%s==1,1,null())' % (r_strahler) , outfile = r_ord_1)
- grass.run_command('r.thin', input = r_ord_1 , output = r_ord_1+'_thin', iterations = 200 )
- grass.run_command('r.to.vect', input = r_ord_1+'_thin' , output = v_ord_1 , feature = 'line', flags = 'v' )
- magnitudo = float(grass.read_command('v.info', map = v_ord_1, layer = 1, flags = 't').split('\n')[2].split('=')[1])
+ # Magnitudo
+ grass.run_command('r.mapcalculator', amap = r_strahler , formula = 'if(%s==1,1,null())' % (r_strahler) , outfile = r_ord_1)
+ 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)
+ grass.run_command('r.to.vect', input = r_ord_1+'_thin',
+ output = v_ord_1,
+ feature = 'line',
+ flags = 'v',
+ overwrite = True)
+ magnitudo = float(grass.read_command('v.info', map = v_ord_1,
+ layer = 1,
+ flags = 't').split('\n')[2].split('=')[1])
- # First order stream frequency
- FSF = magnitudo/area_basin
+ # First order stream frequency
+ FSF = magnitudo / area_basin
- # Statistics
- stream_stats = grass.read_command('r.stream.stats', stream = r_strahler, dir = r_drainage_e, dem = r_elevation+'_crop' )
- grass.message( "output r.stream.stats %s" % stream_stats )
- stream_stats_summary = stream_stats.split('\n')[4].split('|')
- stream_stats_mom = stream_stats.split('\n')[8].split('|')
- Max_order , Num_streams , Len_streams , Stream_freq = stream_stats_summary[0] , stream_stats_summary[1] , stream_stats_summary[2] , stream_stats_summary[5]
- Bif_ratio , Len_ratio , Area_ratio , Slope_ratio = stream_stats_mom[0] , stream_stats_mom[1] , stream_stats_mom[2] , stream_stats_mom[3]
- drainage_density = float(Len_streams) / float(area_basin)
+ # Statistics
+ stream_stats = grass.read_command('r.stream.stats', stream = r_strahler,
+ dir = r_drainage_e,
+ dem = 'r_elevation_crop' )
+ grass.message( "output r.stream.stats %s" % stream_stats )
+ stream_stats_summary = stream_stats.split('\n')[4].split('|')
+ stream_stats_mom = stream_stats.split('\n')[8].split('|')
+ Max_order , Num_streams , Len_streams , Stream_freq = stream_stats_summary[0] , stream_stats_summary[1] , stream_stats_summary[2] , stream_stats_summary[5]
+ Bif_ratio , Len_ratio , Area_ratio , Slope_ratio = stream_stats_mom[0] , stream_stats_mom[1] , stream_stats_mom[2] , stream_stats_mom[3]
+ drainage_density = float(Len_streams) / float(area_basin)
- # Cleaning up
- grass.run_command('g.remove', rast = r_height_average)
- grass.run_command('g.remove', rast = r_aspect_mod)
- grass.run_command('g.remove', rast = r_mainchannel)
- grass.run_command('g.remove', rast = r_stream_e)
- grass.run_command('g.remove', rast = r_drainage_e)
- grass.run_command('g.remove', rast = r_mask)
- grass.run_command('g.remove', rast = r_ord_1)
- grass.run_command('g.remove', rast = r_average_hillslope)
- grass.run_command('g.remove', rast = r_mainchannel_dim)
- grass.run_command('g.remove', rast = r_outlet)
- grass.run_command('g.remove', rast = r_basin)
- grass.run_command('g.remove', rast = prefix+'_mainchannel_thin')
- grass.run_command('g.remove', rast = prefix+'_mainchannel_dim_thin')
- grass.run_command('g.remove', rast = prefix+'_ord_1_thin')
- grass.run_command('g.remove', rast = prefix+'_stream_nothin')
- grass.run_command('g.remove', rast = prefix+'_stream_e_thin')
- grass.run_command('g.remove', vect = v_mainchannel_dim+'_point')
- grass.run_command('g.remove', vect = v_centroid1)
- grass.run_command('g.remove', vect = v_mainchannel_dim)
- grass.run_command('g.remove', vect = v_ord_1)
+ # Cleaning up
+ grass.run_command('g.remove', rast = 'r_elevation_crop')
+ grass.run_command('g.remove', rast = r_height_average)
+ grass.run_command('g.remove', rast = r_aspect_mod)
+ grass.run_command('g.remove', rast = r_mainchannel)
+ grass.run_command('g.remove', rast = r_stream_e)
+ grass.run_command('g.remove', rast = r_drainage_e)
+ grass.run_command('g.remove', rast = r_mask)
+ grass.run_command('g.remove', rast = r_ord_1)
+ grass.run_command('g.remove', rast = r_average_hillslope)
+ grass.run_command('g.remove', rast = r_mainchannel_dim)
+ grass.run_command('g.remove', rast = r_outlet)
+ grass.run_command('g.remove', rast = r_basin)
+ grass.run_command('g.remove', rast = prefix+'_mainchannel_thin')
+ grass.run_command('g.remove', rast = prefix+'_mainchannel_dim_thin')
+ grass.run_command('g.remove', rast = prefix+'_ord_1_thin')
+ grass.run_command('g.remove', rast = prefix+'_stream_e_thin')
+ grass.run_command('g.remove', vect = v_mainchannel_dim+'_point')
+ grass.run_command('g.remove', vect = v_centroid1)
+ grass.run_command('g.remove', vect = v_mainchannel_dim)
+ grass.run_command('g.remove', vect = v_ord_1)
- if nomap :
- grass.run_command('g.remove', vect = v_outlet)
- grass.run_command('g.remove', vect = v_basin)
- grass.run_command('g.remove', vect = v_mainchannel)
- grass.run_command('g.remove', vect = v_stream_e)
- grass.run_command('g.remove', rast = r_accumulation)
- grass.run_command('g.remove', rast = r_drainage)
- grass.run_command('g.remove', rast = r_aspect)
- grass.run_command('g.remove', rast = r_strahler)
- grass.run_command('g.remove', rast = r_shreve)
- grass.run_command('g.remove', rast = r_horton)
- grass.run_command('g.remove', rast = r_hack)
- grass.run_command('g.remove', rast = r_distance)
- grass.run_command('g.remove', rast = r_hillslope_distance)
- grass.run_command('g.remove', rast = r_slope)
+ if nomap :
+ grass.run_command('g.remove', vect = v_outlet)
+ grass.run_command('g.remove', vect = v_basin)
+ grass.run_command('g.remove', vect = v_mainchannel)
+ grass.run_command('g.remove', rast = r_accumulation)
+ grass.run_command('g.remove', rast = r_drainage)
+ grass.run_command('g.remove', rast = r_aspect)
+ grass.run_command('g.remove', rast = r_strahler)
+ grass.run_command('g.remove', rast = r_shreve)
+ grass.run_command('g.remove', rast = r_horton)
+ grass.run_command('g.remove', rast = r_hack)
+ grass.run_command('g.remove', rast = r_distance)
+ grass.run_command('g.remove', rast = r_hillslope_distance)
+ grass.run_command('g.remove', rast = r_slope)
- ####################################################
+ ####################################################
- parametri_bacino = {}
- parametri_bacino["mean_slope"] = float(mean_slope)
- parametri_bacino["mean_elev"] = float(mean_elev)
- parametri_bacino["basin_average_aspect"] = float(basin_average_aspect)
- parametri_bacino["basin_east"] = float(basin_east)
- parametri_bacino["basin_north"] = float(basin_north)
- parametri_bacino["basin_resolution"] = float(basin_resolution)
- parametri_bacino["nw"] = nw
- parametri_bacino["se"] = se
- parametri_bacino["area_basin"] = float(area_basin)
- parametri_bacino["perimeter_basin"] = float(perimeter_basin)
- parametri_bacino["L_orienting_vect"] = float(L_orienting_vect)
- parametri_bacino["prevalent_orientation"] = float(prevalent_orientation)
- parametri_bacino["C_comp"] = float(C_comp)
- parametri_bacino["R_c"] = float(R_c)
- parametri_bacino["mainchannel"] = float(mainchannel)
- parametri_bacino["D_topo"] = float(D_topo)
- parametri_bacino["mainchannel_slope" ] = float(mainchannel_slope)
- parametri_bacino["R_al"] = float(R_al)
- parametri_bacino["S_f"] = float(S_f)
- parametri_bacino["H1"] = float(H1)
- parametri_bacino["H2"] = float(H2)
- parametri_bacino["HM"] = float(HM)
- parametri_bacino["t_c"] = float(t_c)
- parametri_bacino["mean_hillslope_length"] = float(mean_hillslope_length)
- parametri_bacino["magnitudo"] = float(magnitudo)
- parametri_bacino["Max_order"] = float(Max_order)
- parametri_bacino["Num_streams"] = float(Num_streams)
- parametri_bacino["Len_streams"] = float(Len_streams)
- parametri_bacino["Stream_freq"] = float(Stream_freq)
- parametri_bacino["Bif_ratio"] = float(Bif_ratio)
- parametri_bacino["Len_ratio"] = float(Len_ratio)
- parametri_bacino["Area_ratio"] = float(Area_ratio)
- parametri_bacino["Slope_ratio"] = float(Slope_ratio)
- parametri_bacino["drainage_density"] = float(drainage_density)
- parametri_bacino["FSF"] = float(FSF)
- grass.message( "\n" )
- grass.message( "##################################" )
- grass.message( "Morphometric parameters of basin :" )
- grass.message( "##################################\n" )
- grass.message( "Easting Centroid of basin : %s " % basin_east )
- grass.message( "Northing Centroid of Basin : %s " % basin_north )
- grass.message( "Rectangle containing basin N-W : %s , %s " % nw )
- grass.message( "Rectangle containing basin S-E : %s , %s " % se )
- grass.message( "Area of basin [km^2] : %s " % area_basin )
- grass.message( "Perimeter of basin [km] : %s " % perimeter_basin )
- grass.message( "Max Elevation [m s.l.m.] : %s " % H1 )
- grass.message( "Min Elevation [m s.l.m.]: %s " % H2 )
- grass.message( "Elevation Difference [m]: %s " % HM )
- grass.message( "Mean Elevation [m s.l.m.]: %s " % mean_elev )
- grass.message( "Mean Slope : %s " % mean_slope )
- grass.message( "Length of Directing Vector [km] : %s " % L_orienting_vect )
- grass.message( "Prevalent Orientation [degree from north, counterclockwise] : %s " % prevalent_orientation )
- grass.message( "Compactness Coefficient : %s " % C_comp )
- grass.message( "Circularity Ratio : %s " % R_c )
- grass.message( "Topological Diameter : %s " % D_topo )
- grass.message( "Elongation Ratio : %s " % R_al )
- grass.message( "Shape Factor : %s " % S_f )
- grass.message( "Concentration Time (Giandotti, 1934) [hr] : %s " % t_c )
- grass.message( "Length of Mainchannel [km] : %s " % mainchannel )
- grass.message( "Mean slope of mainchannel [percent] : %f " % mainchannel_slope )
- grass.message( "Mean hillslope length [m] : %s " % mean_hillslope_length )
- grass.message( "Magnitudo : %s " % magnitudo )
- grass.message( "Max order (Strahler) : %s " % Max_order )
- grass.message( "Number of streams : %s " % Num_streams )
- grass.message( "Total Stream Length [km] : %s " % Len_streams )
- grass.message( "First order stream frequency : %s " % FSF )
- grass.message( "Drainage Density [km/km^2] : %s " % drainage_density )
- grass.message( "Bifurcation Ratio (Horton) : %s " % Bif_ratio )
- grass.message( "Length Ratio (Horton) : %s " % Len_ratio )
- grass.message( "Area ratio (Horton) : %s " % Area_ratio )
- grass.message( "Slope ratio (Horton): %s " % Slope_ratio )
- grass.message( "\n" )
- grass.message( "Done!")
+ parametri_bacino = {}
+ parametri_bacino["mean_slope"] = float(mean_slope)
+ parametri_bacino["mean_elev"] = float(mean_elev)
+ parametri_bacino["basin_average_aspect"] = float(basin_average_aspect)
+ parametri_bacino["basin_east"] = float(basin_east)
+ parametri_bacino["basin_north"] = float(basin_north)
+ parametri_bacino["basin_resolution"] = float(basin_resolution)
+ parametri_bacino["nw"] = nw
+ parametri_bacino["se"] = se
+ parametri_bacino["area_basin"] = float(area_basin)
+ parametri_bacino["perimeter_basin"] = float(perimeter_basin)
+ parametri_bacino["L_orienting_vect"] = float(L_orienting_vect)
+ parametri_bacino["prevalent_orientation"] = float(prevalent_orientation)
+ parametri_bacino["C_comp"] = float(C_comp)
+ parametri_bacino["R_c"] = float(R_c)
+ parametri_bacino["mainchannel"] = float(mainchannel)
+ parametri_bacino["D_topo"] = float(D_topo)
+ parametri_bacino["mainchannel_slope" ] = float(mainchannel_slope)
+ parametri_bacino["R_al"] = float(R_al)
+ parametri_bacino["S_f"] = float(S_f)
+ parametri_bacino["H1"] = float(H1)
+ parametri_bacino["H2"] = float(H2)
+ parametri_bacino["HM"] = float(HM)
+ parametri_bacino["t_c"] = float(t_c)
+ parametri_bacino["mean_hillslope_length"] = float(mean_hillslope_length)
+ parametri_bacino["magnitudo"] = float(magnitudo)
+ parametri_bacino["Max_order"] = float(Max_order)
+ parametri_bacino["Num_streams"] = float(Num_streams)
+ parametri_bacino["Len_streams"] = float(Len_streams)
+ parametri_bacino["Stream_freq"] = float(Stream_freq)
+ parametri_bacino["Bif_ratio"] = float(Bif_ratio)
+ parametri_bacino["Len_ratio"] = float(Len_ratio)
+ parametri_bacino["Area_ratio"] = float(Area_ratio)
+ parametri_bacino["Slope_ratio"] = float(Slope_ratio)
+ parametri_bacino["drainage_density"] = float(drainage_density)
+ parametri_bacino["FSF"] = float(FSF)
+ grass.message( "\n" )
+ grass.message( "##################################" )
+ grass.message( "Morphometric parameters of basin :" )
+ grass.message( "##################################\n" )
+ grass.message( "Easting Centroid of basin : %s " % basin_east )
+ grass.message( "Northing Centroid of Basin : %s " % basin_north )
+ grass.message( "Rectangle containing basin N-W : %s , %s " % nw )
+ grass.message( "Rectangle containing basin S-E : %s , %s " % se )
+ grass.message( "Area of basin [km^2] : %s " % area_basin )
+ grass.message( "Perimeter of basin [km] : %s " % perimeter_basin )
+ grass.message( "Max Elevation [m s.l.m.] : %s " % H1 )
+ grass.message( "Min Elevation [m s.l.m.]: %s " % H2 )
+ grass.message( "Elevation Difference [m]: %s " % HM )
+ grass.message( "Mean Elevation [m s.l.m.]: %s " % mean_elev )
+ grass.message( "Mean Slope : %s " % mean_slope )
+ grass.message( "Length of Directing Vector [km] : %s " % L_orienting_vect )
+ grass.message( "Prevalent Orientation [degree from north, counterclockwise] : %s " % prevalent_orientation )
+ grass.message( "Compactness Coefficient : %s " % C_comp )
+ grass.message( "Circularity Ratio : %s " % R_c )
+ grass.message( "Topological Diameter : %s " % D_topo )
+ grass.message( "Elongation Ratio : %s " % R_al )
+ grass.message( "Shape Factor : %s " % S_f )
+ grass.message( "Concentration Time (Giandotti, 1934) [hr] : %s " % t_c )
+ grass.message( "Length of Mainchannel [km] : %s " % mainchannel )
+ grass.message( "Mean slope of mainchannel [percent] : %f " % mainchannel_slope )
+ grass.message( "Mean hillslope length [m] : %s " % mean_hillslope_length )
+ grass.message( "Magnitudo : %s " % magnitudo )
+ grass.message( "Max order (Strahler) : %s " % Max_order )
+ grass.message( "Number of streams : %s " % Num_streams )
+ grass.message( "Total Stream Length [km] : %s " % Len_streams )
+ grass.message( "First order stream frequency : %s " % FSF )
+ grass.message( "Drainage Density [km/km^2] : %s " % drainage_density )
+ grass.message( "Bifurcation Ratio (Horton) : %s " % Bif_ratio )
+ grass.message( "Length Ratio (Horton) : %s " % Len_ratio )
+ grass.message( "Area ratio (Horton) : %s " % Area_ratio )
+ grass.message( "Slope ratio (Horton): %s " % Slope_ratio )
+ grass.message( "\n" )
+ grass.message( "Done!")
+
+ except:
+ grass.message( "\n" )
+ grass.message( "##################################" )
+ grass.message( "\n" )
+ grass.message( "ERROR" )
+ grass.message( "Outlet coordinates must belong to the river network!" )
+ grass.message( "Please run r.stream.extract and choose coordinates matching with the extracted stream map" )
+
+ # Set region to original
+ grass.read_command('g.region', flags = 'p', region = 'original')
if __name__ == "__main__":
options, flags = grass.parser()
More information about the grass-commit
mailing list