[GRASS-SVN] r30146 - grass-addons/v.strahler
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Feb 14 04:03:03 EST 2008
Author: florian
Date: 2008-02-14 04:03:02 -0500 (Thu, 14 Feb 2008)
New Revision: 30146
Added:
grass-addons/v.strahler/r.strahler.sh
Modified:
grass-addons/v.strahler/description.html
grass-addons/v.strahler/main.c
grass-addons/v.strahler/strahler.c
grass-addons/v.strahler/strahler.h
grass-addons/v.strahler/write.c
Log:
Added contributions made by Annalisa Minelli and Ivan Marchesini
Modified: grass-addons/v.strahler/description.html
===================================================================
--- grass-addons/v.strahler/description.html 2008-02-14 08:42:16 UTC (rev 30145)
+++ grass-addons/v.strahler/description.html 2008-02-14 09:03:02 UTC (rev 30146)
@@ -1,31 +1,87 @@
<h2>DESCRIPTION</h2>
-<em>v.strahler</em> calculates the Strahler Order for all lines of a given dendritic network. The input vector map must be free of cycles. More than one tree in the input data is allowed. No given flow direction is needed. To find the outlet of each tree, a DEM must be given.
+<em>v.strahler</em> calculates the Strahler Order for all lines of a given dendritic network. The input vector map must be free of cycles. For the elaboration a new imported network or a network extracted from DEM by using <em><a HREF="r.watershed.html">r.watershed</a></em> can be used (in this case the topology has to be cleaned manually). More than one tree in the input data is allowed. No given flow direction is needed. To find the outlet of each tree, a DEM must be given.
<h2>NOTES</h2>
<h3>Problems</h3>
-This program is in alpha status. It has the following shortcomings: <br />
-- Robust calculation of Strahler Order is not guaranteed <br />
-- Strahler Order is not attached as a layer, it is merely written to an ASCII file <br />
-- The input data has to be topologically clean. Use sloppy=0.0 <br />
+This program is in beta status. It has the following shortcomings: <br />
+- The input data has to be topologically clean. Use <em>sloppy</em>=0.0 <br />
- Source code comments are not doxydized. <br />
-...
+
<h3>How the algorithm works</h3>
In a first step, <em>v.strahler</em> identifies all separate networks in the input dataset. That means, all connected lines are assigned a common <em>Basin ID</em>.<br />
Consequently, the lowest leaf of each tree is identified as outlet. An auxiliary DEM is needed at this point. <br />
-Strahler ordering begins at each leaf of the tree with order N=1 (excluding the outlet). At a confluence, the order N(r) of the resulting stream is equal to the highest order N(max) of the joining streams or is raised by one if there are two or more joining streams of order N(max).
+Strahler ordering begins at each leaf of the tree with order N=1 (excluding the outlet). At a confluence, the order N(r) of the resulting stream is equal to the highest order N(max) of the joining streams or is raised by one if there are two or more joining streams of order N(max). <br />
+The algorithm returns an ASCII text file with columns: <em>Category</em> (from input map), <em>Line</em> (topology), <em>Basin</em>, <em>Order</em>; and the output map has the Strahler Order value instead of "category" for each line, and no connection with the database.
<h2>EXAMPLE</h2>
-
-TODO
-
+The input map (vector on DEM): <br />
+<br />
+<img src="images/example_input.jpg"> <br />
+<br />
+<br />
+An example of the sintax in GRASS shell: <br />
+<br />
+<div class="code"><pre>
+GRASS 6.3.0RC1 (ED50_Z33):~ > v.strahler input=esp_1 at mapset output=esp_1 dem=dem_20 at mapset
+txout=/home/mapset/esp_1 sloppy=0 layer=1 --overwrite 2>pippo
+</pre></div>
+<br />
+<br />
+An extract from the resultant text file: <br />
+<br />
+<div class="code"><pre>
+== Result of Strahler Order ==
+ Category: Line: Basin: Order:
+ 22 1 1 4
+ 73 2 1 1
+ 25 3 1 4
+ 27 4 1 4
+ 39 5 1 1
+ 48 6 1 1
+ 56 7 1 4
+ 55 8 1 2
+ 88 9 1 4
+ 59 10 1 1
+ 60 11 1 1
+ 61 12 1 2
+ 83 13 1 2
+ 82 14 1 1
+ 91 15 1 1
+ 87 16 1 4
+ 95 17 1 2
+ 96 18 1 1
+ 106 19 1 1
+ 102 20 1 1
+ 104 21 1 4
+ 111 22 1 1
+ 112 23 1 3
+ 121 24 1 3
+ 126 25 1 1
+ 128 26 1 3
+ 137 27 1 3
+</pre></div>
+<br />
+<br />
+The output map (vector on DEM) form a snapshot of GRASS display: <br />
+<br />
+<img src="images/example_output.jpg"> <br />
+<br />
+<br />
+The output map (vector only) from Qgis (different colors for different orders): <br />
+<br />
+<img src="images/Qgis.jpg"> <br />
+<br />
+<br />
<h2>SEE ALSO</h2>
+<em><a href="r.watershed.html">r.watershed</a></em><BR>
-
<h2>AUTHOR</h2>
-Florian Kindl, Univ. Innsbruck
+Florian Kindl, Univ. Innsbruck.<br />
+<br />
+Modified by: Ivan Marchesini and Annalisa Minelli, Univ. Perugia.<br />
-<p><i>Last changed: $Date: 2006/08/05 15:29:16 $</i>
+<p><i>Last changed: $Date: 2007/11/25 15:29:16 $</i>
Modified: grass-addons/v.strahler/main.c
===================================================================
--- grass-addons/v.strahler/main.c 2008-02-14 08:42:16 UTC (rev 30145)
+++ grass-addons/v.strahler/main.c 2008-02-14 09:03:02 UTC (rev 30146)
@@ -1,8 +1,10 @@
-/****************************************************************
+/****************************************************************************
*
* MODULE: v.strahler
*
* AUTHOR(S): Florian Kindl, Norbert Pfeifer, The GRASS development team
+ * Modified by: Ivan Marchesini <marchesini unipg.it> and
+ Annalisa Minelli <annapatri tiscali.it>
*
* PURPOSE: Assign Strahler order to dendritic network
*
@@ -13,15 +15,12 @@
* Read the file COPYING that comes with GRASS
* for details.
*
- **************************************************************/
+ ****************************************************************************/
/*! \file main.c
\brief Assign Strahler order to dendritic network
\author Florian Kindl
-\todo Add clean way to write output to *Out
- \li Outpout is now text-only
- \li Add fields in table for Strahler Order and BasinID
\todo Deal with poor topology
\li implemement sloppy mode or
\li force clean topology
@@ -47,6 +46,7 @@
#define EBUG 1
#define REATEDBLAH 1
+int ret, type; /* return and type related to Vect_read_line and Vect_write_line */
int fdrast; /* file descriptor for raster file is int */
int ntrees; /* number of trees in dataset, returned by StrahForestToTrees */
double sloppy;
@@ -54,7 +54,7 @@
int main(int argc, char **argv)
{
int line, node;
- char buf[1024];
+ char buf[1024];
int nnodes, nlines, tlines;
DBBUF *dbbuf;
NODEV *nodev;
@@ -66,20 +66,25 @@
FILE *txout;
struct Cell_head window; /* for cell sampling */
- extern int ntrees; /* number of trees calculated by StrahForestToTrees */
+ extern int ntrees; /* number of trees calculated by StrahForestToTrees */
extern int fdrast;
- extern double sloppy;
+ extern double sloppy;
/* Attribute table (from v.net.path/path.c) */
+ struct line_cats *Cats; /* introduced and initialized the structures line_cats and line_pnts */
+ static struct line_pnts *Points;
+ Points = Vect_new_line_struct ();
+ Cats = Vect_new_cats_struct ();
dbString sql;
dbDriver *driver;
dbColumn *column;
struct field_info *Fi;
- int field;
-
+ int field;
+
/* Initialize the GIS calls */
+
G_gisinit (argv[0]) ;
input = G_define_standard_option(G_OPT_V_INPUT);
@@ -101,9 +106,7 @@
txout_opt->description = _("Path to ASCII file where results will be written");
-/*
-option sloppy for bad topology
-*/
+ /* option sloppy for bad topology */
sloppy_opt = G_define_option();
sloppy_opt->key = "sloppy";
@@ -126,19 +129,22 @@
G_fatal_error (_("Could not find input input <%s>"), input->answer);
}
- /* open datasets at topology level 2 */
+ /* open datasets at topology level 2 */
+
Vect_set_open_level(2);
Vect_open_old (&In, input->answer, inputset);
G_debug(1, "Input vector opened");
- /* Open new vector, make 3D if input is 3D */
+ /* Open new vector, make 3D if input is 3D */
+
if (1 > Vect_open_new( &Out, output->answer, Vect_is_3d( &In ))) {
Vect_close( &In );
G_fatal_error("Failed opening output vector file %s", output->answer);
}
/* Open DEM */
+
demset = G_find_cell2( dem_opt->answer, NULL);
if ( demset == NULL ) {
G_fatal_error ("DEM file %s not found", dem_opt->answer);
@@ -148,6 +154,7 @@
}
/* Open text file for results if given */
+
if ( txout_opt->answer ) {
txout = fopen( txout_opt->answer, "w" );
if ( txout == NULL )
@@ -159,19 +166,18 @@
/* write history */
- Vect_copy_head_data( &In, &Out );
+ /*Vect_copy_head_data( &In, &Out );
Vect_hist_copy( &In, &Out );
- Vect_hist_command( &Out );
+ Vect_hist_command( &Out );*/
#ifdef COPYBLAH
/* this works, we want to continue here and ADD two columns to the &Out */
/* Copy input to output (from v.clean/main.c) */
- G_debug( 1, "Copying vector lines to Output" );
+ /*G_debug( 1, "Copying vector lines to Output" );
/* This works for both level 1 and 2 */
- Vect_copy_map_lines ( &In, &Out );
- Vect_copy_tables ( &In, &Out, 0 );
- /* 0: copy all fields, else field number */
+ /*Vect_copy_map_lines ( &In, &Out ); */
+ /* 0: copy all fields, else field number */
#endif
@@ -260,7 +266,7 @@
/* initialize properly */
for ( line=1; line <= nlines; line++) {
- dbbuf[line].line = dbbuf[line].bsnid = dbbuf[line].sorder = 0;
+ dbbuf[line].category = dbbuf[line].line = dbbuf[line].bsnid = dbbuf[line].sorder = 0;
}
for ( node=1; node <= nnodes; node++) {
nodev[node].node = nodev[node].degree = nodev[node].visited = 0;
@@ -313,8 +319,27 @@
db_close_database_shutdown_driver( driver );
#endif
+
+ /* writing StrahOrder instead of category in the output */
+
+ for ( line=1; line<=nlines; line++) {
+ type=Vect_read_line(&In, Points, Cats, line);
+ ret = Vect_cat_del (Cats, field);
+ G_debug(4, "deleted categories, ret=%d", ret);
+ Vect_cat_set (Cats, field, dbbuf[line].sorder);
+ Vect_write_line ( &Out, type, Points, Cats );
+ G_debug(4, "category written for line %d", line);
+ }
- /* Write text file for results if given */
+ /* assign category values */
+
+ for ( line=1; line<=nlines; line++) {
+ dbbuf[line].category = Vect_get_line_cat ( &In, line, field );
+ G_debug(4, "added category %d for line %d", dbbuf[line].category, dbbuf[line].line);
+ }
+
+ /* Write text file for results if given */
+
if ( txout_opt->answer )
StrahWriteToFile( dbbuf, nlines, txout );
Added: grass-addons/v.strahler/r.strahler.sh
===================================================================
--- grass-addons/v.strahler/r.strahler.sh (rev 0)
+++ grass-addons/v.strahler/r.strahler.sh 2008-02-14 09:03:02 UTC (rev 30146)
@@ -0,0 +1,174 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.strahler.sh
+# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
+# PURPOSE: Create a vector map of strahler ordered streems of a single
+# basin, starting from a DEM.
+# The script extracts from a DEM the network (by using
+# r.watershed), converts the network to vector lines,
+# cleans the topology, lets you decide an outlet for your basin,
+# and executes v.strahler for the upstream basin.
+#
+# COPYRIGHT: (C) 2008 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.
+#
+# TODO: solve stability problems for low area threshold
+#############################################################################
+#%Module
+#% description: Create a vector map of strahler ordered streems of a single basin starting from a DEM
+#% keywords: strahler, streems, 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: thr
+#% type: integer
+#% description: minimum size of exterior watershed basin
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Name of streems ordered map created
+#% required : yes
+#%end
+#%option
+#% key: textoutput
+#% type: string
+#% gisprompt: new_file,file,output
+#% key_desc: name
+#% description: Name for output text file
+#% required: yes
+#%end
+#%option
+#% key: bkgrmap
+#% type: string
+#% key_desc: bkgrmap
+#% gisprompt: old,cell,raster
+#% description: Name of background map to plot
+#% required : no
+#%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
+thr=$GIS_OPT_THR
+output=$GIS_OPT_OUTPUT
+textoutput=$GIS_OPT_TEXTOUTPUT
+bkgrmap=$GIS_OPT_BKGRMAP
+
+#check presence of raster MASK, put it aside
+MASKFOUND=0
+eval `g.findfile element=cell file=MASK`
+if [ "$file" ] ; then
+ g.message "Raster MASK found, temporarily disabled"
+ g.rename MASK,"${TMPNAME}_origmask" --quiet
+ MASKFOUND=1
+fi
+
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET
+
+#uncomment this line if you need to remove some maps remaining from previous runs
+#g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
+
+#get resolution of DEM
+res=`g.region -p | grep nsres | cut -f2 -d':' | tr -d ' '`
+
+#find raster streems and plot them
+r.watershed elevation=$dem threshold=$thr stream=rnetwork drainage=drainage --overwrite
+r.null map=rnetwork setnull=0
+r.thin input=rnetwork output=rnetwork_th iterations=200 --overwrite
+r.mapcalc "rnetwork_thin=rnetwork_th/rnetwork_th"
+d.erase
+d.rast map=$bkgrmap
+d.rast -o map=rnetwork_thin
+
+#ask the user for zooming and choosing the outlet cell
+#-------
+echo "
+
+
+Now, please zoom to the area where you want to put the outlet cross section
+
+
+"
+
+d.zoom
+
+echo "
+
+
+Now, please click on the cell representing the cross section
+
+
+"
+
+#prevent from choosing a wrong cell outside from the streems
+cat=2
+while [ "$cat" != "1" ]
+do
+result=`d.what.rast -t -1 rnetwork_thin`
+coor=`echo $result | cut -f1 -d ' '`
+x=`echo $coor | cut -f1 -d':'`
+y=`echo $coor | cut -f2 -d':'`
+cat=`echo $result | cut -f3 -d ' ' | tr -d :`
+done
+
+#come back to the previous zoom
+d.zoom -r
+
+#find the basin and set the mask
+r.water.outlet drainage=drainage basin=basin easting=$x northing=$y
+d.rast basin
+r.null map=basin setnull=0
+r.mapcalc "MASK=basin"
+d.rast MASK
+
+#clean raster map and convert it to vector
+g.region rast=MASK
+r.to.vect input=rnetwork_thin output=vnetwork feature=line --overwrite
+d.vect map=vnetwork
+v.type input=vnetwork output=vnetwork_b type=line,boundary --overwrite
+v.centroids input=vnetwork_b output=vnetwork_b_add option=add layer=1 cat=1 step=1 --overwrite
+newres=`echo "$res/4" | bc -l`
+g.region res=$newres
+v.to.rast input=vnetwork_b_add output=rnetwork_b use=val type=area layer=1 value=1 rows=4096 --overwrite
+v.to.rast input=vnetwork output=rnetwork_l use=val type=line layer=1 value=1 rows=4096 --overwrite
+r.patch input=rnetwork_b,rnetwork_l output=rnetwork_patch --overwrite
+r.thin input=rnetwork_patch output=rnetwork_thin2 iterations=200 --overwrite
+r.to.vect input=rnetwork_thin2 output=vnetwork2 feature=line --overwrite
+
+#remove dangles and create polylines
+soglia=`echo "sqrt(2*($res^2))+1" | bc -l`
+v.clean input=vnetwork2 output=vnetwork_dangle type=line tool=rmdangle thresh=$soglia --overwrite
+v.build.polylines input=vnetwork_dangle output=vnetwork_poly cats=first --overwrite
+d.vect vnetwork_poly
+
+#enlarge region to ensure v.strahler get raster values at streems ends and use v.strahler
+g.region -a res=$res n=n+$soglia s=s-$soglia e=e+$soglia w=w-$soglia
+v.strahler input=vnetwork_poly output=$output dem=$dem sloppy=0 layer=1 txout=$textoutput
+
+#remove temporary files
+g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
+
Property changes on: grass-addons/v.strahler/r.strahler.sh
___________________________________________________________________
Name: svn:executable
+ *
Modified: grass-addons/v.strahler/strahler.c
===================================================================
--- grass-addons/v.strahler/strahler.c 2008-02-14 08:42:16 UTC (rev 30145)
+++ grass-addons/v.strahler/strahler.c 2008-02-14 09:03:02 UTC (rev 30146)
@@ -15,7 +15,7 @@
nnodes = Vect_get_num_nodes( In );
- z = z_init = 27000000.0; /* is it safe to initialize lowest z with the height of Olympus Mons above Martian Datum in millimeters? */
+ z = z_init = 27000000.0; /* is it safe to initialize lowest z with the height of Olympus Mons above Martian Datum in millimeters? */
outlets = (OUTLETS *) G_malloc ( (ntrees + 1) * ((int)sizeof (OUTLETS)) );
@@ -41,6 +41,7 @@
if ( degr == 1 ) {
unode = node;
+ nodev[node].visited = 1; /*set the "visited" parameter to 1... can't be less*/
/*aline = abs( Vect_get_node_line( In, unode, 0 ) );*/
aline = abs( StrahGetNodeLine( In, unode, 0 ) );
/*Vect_get_line_nodes( In, aline, &unode, &dnode);*/
@@ -88,10 +89,10 @@
int StrahOrder( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev )
{
int nlines, line;
- int fnode, tnode, unode, dnode; /* from-, to-, up-, down- node */
+ int fnode, tnode, unode, dnode; /* from-, to-, up-, down- node */
int degr, d, aline, aorder; /* find adjacent lines */
- int corder, norder, dorder, dline; /* assign order */
- int cline, rfinish; /* follow one line until this stop-condition is 1 */
+ int corder, norder, dorder, dline; /* assign order */
+ int cline, rfinish; /* follow one line until this stop-condition is 1 */
struct dbbuf;
struct nodev;
@@ -101,9 +102,9 @@
G_debug( 1, "reached StrahOrder");
for (line = 1; line <= nlines; line++) {
- if (dbbuf[line].sorder == 1) { /* get lines of order 1 */
+ if (dbbuf[line].sorder == 1) { /* get lines of order 1 */
- cline = line; /* and start run downwards */
+ cline = line; /* and start run downwards */
rfinish = 0;
while (rfinish == 0) {
G_debug( 3, "reached line %d", cline);
@@ -160,8 +161,8 @@
G_debug( 4, "deg for dnode %d is %d", dnode, degr );
- corder = dbbuf[cline].sorder; /* current order - result won't be less than that */
- norder = dorder = 0; /* no order - how many lines have none?, highest order of others */
+ corder = dbbuf[cline].sorder; /* current order - result won't be less than that */
+ norder = dorder = 0; /* no order - how many lines have none?, highest order of others */
for (d = 0; d < degr; d++) {
/*aline = abs( Vect_get_node_line ( In, dnode, d ) );*/
aline = abs( StrahGetNodeLine( In, dnode, d ) );
@@ -176,27 +177,27 @@
}
}
}
-
- if (norder > 1 || norder == 0) { /* if (norder > 1: node indeterminate OR norder == 0 : subtree finished) */
+
+ if (norder > 1 || norder == 0) { /* if (norder > 1: node indeterminate OR norder == 0 : subtree finished) */
G_debug( 3, "Finished run at line %d because norder=%d", cline, norder);
G_debug( 4, "visiting unode %d for line %d\n", unode, cline);
- nodev[unode].visited += 1; /* visit unode, for the sake of completeness */
- rfinish = 1; /* finish this run */
+ nodev[unode].visited += 1; /* visit unode, for the sake of completeness */
+ rfinish = 1; /* finish this run */
- } else { /* or */
+ } else { /* or */
if (dorder > corder) {
- continue; /* assign highest order of alines */
+ ; /* assign highest order of alines */
} else if (dorder < corder) {
dorder=corder; /* or keep order of cline */
} else {
- dorder++; /* or raise order by one */
+ dorder++; /* or raise order by one */
}
- dbbuf[dline].sorder = dorder; /* assign order */
+ dbbuf[dline].sorder = dorder; /* assign order */
G_debug( 4, "visiting dnode %d and unode %d for line %d", dnode, unode, cline);
- nodev[dnode].visited += 1; /* visit dnode */
- nodev[unode].visited += 1; /* visit unode, for the sake of completeness */
+ nodev[dnode].visited += 1; /* visit dnode */
+ nodev[unode].visited += 1; /* visit unode, for the sake of completeness */
cline = dline; /* and continue with this line */
Modified: grass-addons/v.strahler/strahler.h
===================================================================
--- grass-addons/v.strahler/strahler.h 2008-02-14 08:42:16 UTC (rev 30145)
+++ grass-addons/v.strahler/strahler.h 2008-02-14 09:03:02 UTC (rev 30146)
@@ -20,6 +20,7 @@
*/
typedef struct { /* collect lines in Strahler order */
+ int category;
int line; /* line number */
int bsnid; /* basin ID */
int sorder; /* Strahler order */
Modified: grass-addons/v.strahler/write.c
===================================================================
--- grass-addons/v.strahler/write.c 2008-02-14 08:42:16 UTC (rev 30145)
+++ grass-addons/v.strahler/write.c 2008-02-14 09:03:02 UTC (rev 30146)
@@ -6,9 +6,9 @@
G_debug( 1, "Reached WriteToFile - writing %d lines", nlines);
fprintf( txout, "== Result of Strahler Order ==\n" );
- fprintf( txout, " Line: Basin: Order:\n" );
+ fprintf( txout, " Category: Line: Basin: Order:\n" ); /*added cat column*/
for ( line=1; line<=nlines; line++) {
- fprintf( txout, "%8d%8d%8d\n",dbbuf[line].line, dbbuf[line].bsnid, dbbuf[line].sorder);
+ fprintf( txout, "%9d%8d%8d%8d\n", dbbuf[line].category, dbbuf[line].line, dbbuf[line].bsnid, dbbuf[line].sorder );
}
return 1;
}
More information about the grass-commit
mailing list