[GRASS-SVN] r30780 - grass-addons/vector

[svn_grass at osgeo.org]

Author: annalisapg
Date: 2008-03-28 11:05:22 -0400 (Fri, 28 Mar 2008)
New Revision: 30780

Added:
   grass-addons/vector/v.mainchannel
Log:
Main channel individuation on a vector stream network

Added: grass-addons/vector/v.mainchannel
===================================================================
--- grass-addons/vector/v.mainchannel	                        (rev 0)
+++ grass-addons/vector/v.mainchannel	2008-03-28 15:05:22 UTC (rev 30780)
@@ -0,0 +1,113 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE:	v.mainchannel
+# AUTHOR(S):	Annalisa Minelli, Ivan Marchesini
+# PURPOSE:	Find the main channel of a vector stream network (useful to compare with a network ordered by Horton-Strahler)
+# COPYRIGHT:	(C) 2005 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
+#		for details.
+#
+###########################################################################
+#%Module
+#%  description: Find the main channel of a vector stream network (useful to compare with a network ordered by Horton-Strahler)
+#%  keywords: strahler, main_channel, vector
+#%End
+#%option
+#% key: dem
+#% type: string
+#% key_desc: dem
+#% gisprompt: old,cell,raster
+#% description: Name of DEM raster map
+#% required : yes
+#%END
+#%option
+#% key: xcoor
+#% type: double
+#% description: x coord of outlet (read by dem in the pixel where is the vector file cross section)
+#% required : yes
+#%end
+#%option
+#% key: ycoor
+#% type: double
+#% description: y coord of outlet (read by dem in the pixel where is the vector file cross section)
+#% required : yes
+#%end
+#%option
+#% key: cost
+#% type: integer
+#% description: Distance between isolines (separated by <comma>)
+#% required : yes
+#%end
+#%option
+#% key: input
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Vector map of Strahler-ordered streams
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Name of the main_channel vector map
+#% required : yes
+#%END
+
+if  [ -z "$GISBASE" ] ; then
+    echo "You must be in GRASS GIS to run this program." >&2
+    exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+    exec g.parser "$0" "$@"
+fi
+
+dem=$GIS_OPT_DEM
+xcoor=$GIS_OPT_XCOOR
+ycoor=$GIS_OPT_YCOOR
+cost=$GIS_OPT_COST
+input=$GIS_OPT_INPUT
+output=$GIS_OPT_OUTPUT
+
+### setup enviro vars ###
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+
+g.remove vect=end_basin,out_1del,out_1add,out_basin,out_isodistance,out_2,out_2del,out_2points,out_2addp,out_3D,out_path
+
+echo "$xcoor|$ycoor|1" | v.in.ascii out=end_basin
+thr=`g.region -p | grep nsres | cut -f2 -d':' | tr -d ' '`
+thresh=`echo "$thr*1.5" | bc -l | cut -f1 -d'.'`
+v.category input=$input output=out_1del type=line option=del cat=1 layer=1 step=1
+v.category input=out_1del output=out_1add type=line option=add cat=1 layer=1 step=1
+v.net input=out_1add points=end_basin output=out_basin operation=connect alayer=1 nlayer=2 thresh="$thresh"
+v.net.iso input=out_basin output=out_isodistance type=line alayer=1 nlayer=2 ccats=1-99999999 costs=$cost
+v.category input=out_isodistance opt=report -g > temp1
+cat temp1 | grep line | cut -f5 -d ' ' > temp2
+v.extract input=out_isodistance output=out_2 type=line layer=1 new=-1 file=temp2
+v.category input=out_2 output=out_2del type=point,line,boundary,centroid,area option=del cat=1 layer=1 step=1
+v.net input=out_2del output=out_2points operation=nodes alayer=1 nlayer=2
+v.category input=out_2points output=out_2addp type=point option=add cat=1 layer=2 step=1
+v.drape input=out_2addp type=point,line rast=$dem method=nearest output=out_3D
+v.db.addtable map=out_3D table=tab_coord layer=2 'columns=cat integer, x double precision, y double precision, z double precision'
+v.to.db map=out_3D type=point layer=2 qlayer=1 option=coor units=meters 'column=x, y, z'
+v.univar -g map=out_3D type=point column=z layer=2 percentile=90 > temp
+max=`cat temp | grep max | cut -f2 -d'='`
+db_driver=`db.connect -p | grep driver | cut -f2 -d':'`
+db_database=`db.connect -p | grep database | cut -f2 -d':'`
+x_max=`db.select -c table=tab_coord database=$db_database driver=$db_driver 'sql=select x from tab_coord where z='$max''`
+y_max=`db.select -c table=tab_coord database=$db_databese driver=$db_driver 'sql=select y from tab_coord where z='$max''`
+echo "1 $x_max $y_max $xcoor $ycoor" |  v.net.path input=out_isodistance output=$output type=line alayer=1 nlayer=2
+
+#plot the result of elaboration in relation to the strahler ordered stream network
+d.erase
+d.vect $input display=shape,cat
+d.vect map=$output color=255:0:0 lcolor=0:0:0 fcolor=170:170:170 display=shape type=line icon=basic/x size=5 width=2 layer=1 lsize=8
+d.barscale tcolor=0:0:0 bcolor=255:255:255  at=2,2
+
+g.remove vect=end_basin,out_1del,out_1add,out_basin,out_2,out_2del,out_2points,out_2addp,out_3D
+#you can modify the "cost" distances in order to have main channel including the max_ordered channel


Property changes on: grass-addons/vector/v.mainchannel
___________________________________________________________________
Name: svn:executable
   + *