[GRASS-SVN] r40887 - grass-addons/raster/r.stream.order

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Feb 9 03:46:50 EST 2010


Author: jarekj71
Date: 2010-02-09 03:46:49 -0500 (Tue, 09 Feb 2010)
New Revision: 40887

Modified:
   grass-addons/raster/r.stream.order/description.html
   grass-addons/raster/r.stream.order/global.h
   grass-addons/raster/r.stream.order/io.c
   grass-addons/raster/r.stream.order/main.c
   grass-addons/raster/r.stream.order/order.c
   grass-addons/raster/r.stream.order/orders.png
Log:
new ordering system: topological distance; some minor bugfixes

Modified: grass-addons/raster/r.stream.order/description.html
===================================================================
--- grass-addons/raster/r.stream.order/description.html	2010-02-09 08:30:16 UTC (rev 40886)
+++ grass-addons/raster/r.stream.order/description.html	2010-02-09 08:46:49 UTC (rev 40887)
@@ -14,8 +14,13 @@
 Also <em>stream</em> network and <em>accumulation</em> maps must have the same resolution. It is checked by default. If resolutions differ the module informs about it and stops. Region boundary
 and maps boundary may be differ but it may lead to unexpected results.</DD>
 <p>
+<DT><b>table</b></DT>
+<DD>Table where stream network topology can be stored. Because r.stream.order is prepared to work both with r.watershed and r.stream.extract, table by default is not attached to vector, but if stream network is produced by r.stream.extract it can be simply added to file using <a href="v.db.coonect.html">v.db.connect</a>. See DESCRIPTION for details</DD>
+
+<p>
 <DT><b>accum</b></DT>
 <DD>Flow accumulation (optional, not recommended): name of flow accumulation file produced by r.watershed or used in r.stream.extract. This map is an option only if Horton's or Hack's ordering is performed. Normally both Horton and Hack ordering is calculated on cumulative stream lrngth wchich is calculated internaly. Flow accumulation can be used if user want to calculate main stream as most accumulated stream. Flow accumulation map shall be of DCELL type, as is by default produced by r.watershed or converted do DCELL with r.mapcalc.</DD>
+
 <h2>OUTPUTS</h2>
 
 <P>At least one output map is required: </p>
@@ -29,8 +34,11 @@
 <DD>Name of Horton's stream order output map (require accum file): see notes for detail.</DD>
 
 <DT><b>hack</b></DT>
-<DD>Name of Hack's main streams output map (require accum file): see notes for detail.</DD>
+<DD>Name of Hack's main streams output map : see notes for detail.</DD>
 
+<DT><b>top</b></DT>
+<DD>Name of topological dimensions streams output map: see notes for detail.</DD>
+
 </DL>
 
 <h2>DESCRIPTION</h2>
@@ -72,10 +80,30 @@
 <BR>
 <B>Advantages and disadvantages of main stream ordering:</B>
 The biggest advantage of that method is the possibility to compare and analyze topology upstream, according to main streams. Because all tributaries of main channel have order of 2, streams can be quickly and easily filtered and its proprieties and relation to main stream determined. The main disadvantage of that method is the problem with the comparison of subcatchment topology of the same order. Subcatchments of the same order may be both highly branched and widespread in the catchment area and a small subcatchment with only one stream. 
+<H4>Topological dimension streams order</H4>
+This method of ordering calculates topological distance of every stream from catchment outlet. The topopological distance is defined as the number of segments which separates the current segment from the outlet basin 
+<BR>
 
+
+<H4>Stream network topology table description</H4>
+	<li><b>cat</b> integer: category;
+	<li><b>stream</b>integer: stream number, usually equal to cat;
+	<li><b>next_stream</b> integer: stream to which contribute current stream (downstream);
+	<li><b>prev_streams</b>; two or more contributing streams (upstream);
+	<li><b>strahler</b> integer: Strahler's stream order:
+	<li><b>horton</b> integer: Hortons's stream order:
+	<li><b>shreve</b> integer: Shreve's stream magnitude;
+	<li><b>hack</b> integer: Hack's main streams order;
+	<li><b>topo</b> integer: Topological dimension streams order;
+	<li><b>length</b> double precision: stream length;
+	<li><b>cum_length</b> double precision: length of stream from source;
+	<li><b>out_dist</b> double precision: distance of current stream init from outlet;
+	<li><b>stright</b> double precision: length of stream as stright line;
+	<li><b>fractal</b> double precision: fractal dimention: stream length/stright stream length
+
 <h2>NOTES</H2>
 <P>
-Module can work only if direction map, stream map and region map has same settings. It is also required that stream map and direction map come from the same source. For lots of reason this limitation probably cannot be omitted.   this means if stream map comes from r.stream.extract also direction map from r.stream.extract must be used. If stream network was generated with MFD method also MFD direction map must be used. Nowadays f direction map comes from r.stream.extract  must be patched by direction map from r.watershed. (with r.patch). 
+Module can work only if direction map, stream map and region map has same settings. It is also required that stream map and direction map come from the same source. For lots of reason this limitation probably cannot be omitted. This means if stream map comes from r.stream.extract also direction map from r.stream.extract must be used. If stream network was generated with MFD method also MFD direction map must be used. Nowadays f direction map comes from r.stream.extract  must be patched by direction map from r.watershed. (with r.patch). 
 
 <h2>SEE ALSO</h2>
 
@@ -88,19 +116,21 @@
 </em>
 
 
-
 <h2>REFERENCES</h2>
+Claps, P., Fiorentino, M., Oliveto, G., (1994), <i>Informational entropy of fractal river networks</i>,
+Journal of Hydrology, 187(1-2), 145-156 .</p>
 Hack, J., (1957), <i>Studies of longitudinal stream profiles in Virginia and Maryland</i>, 
 <b>U.S. Geological Survey Professional Paper</b>, 294-B<p>
 Horton, R. E. (1945), <i>Erosional development of streams and their drainage basins: hydro-physical approach to quantitative morphology</i>, 
-<b>Geological Society of America Bulletin</b> 56 (3): 275-370<BR>
-<a h
-<p>
+<b>Geological Society of America Bulletin</b> 56 (3): 275-370<p>
 Shreve, R., <i>Statistical Law of Stream Numbers</i>, <b>J. Geol.</b>, 74, (1966), 17-37.<p>
 Strahler, A. N. (1952), <i>Hypsometric (area-altitude) analysis of erosional topology</i>, 
 <b>Geological Society of America Bulletin</b> 63 (11): 1117–1142<p>
 Strahler, A. N. (1957), <i>Quantitative analysis of watershed geomorphology</i>, 
 <b>Transactions of the American Geophysical Union</b> 8 (6): 913–920.<p>
+
+
+
 <h2>AUTHOR</h2>
 Jarek  Jasiewicz, Markus Metz
 

Modified: grass-addons/raster/r.stream.order/global.h
===================================================================
--- grass-addons/raster/r.stream.order/global.h	2010-02-09 08:30:16 UTC (rev 40886)
+++ grass-addons/raster/r.stream.order/global.h	2010-02-09 08:46:49 UTC (rev 40887)
@@ -49,6 +49,8 @@
     double length;
     double stright;
     double fractal;
+    double distance; /* distance to outlet */
+    int topo_dim; /* topological dimension */
     /* float length; */
 };
 
@@ -85,7 +87,7 @@
 
 GLOBAL struct Cell_head window;
 GLOBAL char *in_dirs, *in_streams, *in_vector, *in_table, *in_accum;	/* input dirrection and accumulation raster names */
-GLOBAL char *out_strahler, *out_shreeve, *out_hack, *out_horton;	/* output strahler, shreeve and class raster names */
+GLOBAL char *out_strahler, *out_shreeve, *out_hack, *out_horton, *out_topo;	/* output strahler, shreeve and class raster names */
 GLOBAL int out_zero; /* zero-value background */
 
 GLOBAL CELL **dirs, **streams;	/* matrix with input data */

Modified: grass-addons/raster/r.stream.order/io.c
===================================================================
--- grass-addons/raster/r.stream.order/io.c	2010-02-09 08:30:16 UTC (rev 40886)
+++ grass-addons/raster/r.stream.order/io.c	2010-02-09 08:46:49 UTC (rev 40887)
@@ -124,6 +124,8 @@
 	cur_mapset = G_find_cell2(in_streams, "");
 	G_read_range(in_streams,cur_mapset,&stream_range);
 	G_get_range_min_max(&stream_range, &c_min, &c_max);
+		if(c_max<1)
+	G_fatal_error("No streams found"	);
 	return c_max;
 
 }
@@ -131,8 +133,8 @@
 int write_maps(void)
 {
     int r, c;
-    CELL *strahler_buf, *shreeve_buf, *hack_buf, *horton_buf;
-    int out_str_fd, out_shr_fd, out_hck_fd, out_hrt_fd;	/* output file descriptors: outstr_fd - strahler... etc */
+    CELL *strahler_buf, *shreeve_buf, *hack_buf, *horton_buf, *topo_buf;
+    int out_str_fd, out_shr_fd, out_hck_fd, out_hrt_fd, out_topo_fd;	/* output file descriptors: outstr_fd - strahler... etc */
     struct History history;	/* holds meta-data (title, comments,..) */
 
     /* routine check if legal file names and we able to opuen files */
@@ -160,6 +162,12 @@
 	    G_fatal_error(_("Unable to create raster map <%s>"), out_horton);
 	horton_buf = G_allocate_cell_buf();
     }
+    
+     if (out_topo) {
+	if ((out_topo_fd = G_open_raster_new(out_topo, CELL_TYPE)) < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), out_topo);
+	topo_buf = G_allocate_cell_buf();
+    }
 
     G_message(_("Writing maps..."));
 
@@ -174,6 +182,8 @@
 	    G_set_c_null_value(hack_buf, ncols);
 	if (out_horton)
 	    G_set_c_null_value(horton_buf, ncols);
+	if (out_topo)
+	    G_set_c_null_value(topo_buf, ncols);    
 
 	for (c = 0; c < ncols; ++c) {
 	    if (!out_zero) {
@@ -186,6 +196,8 @@
 			hack_buf[c] = s_streams[streams[r][c]].hack;
 		    if (out_horton)
 			horton_buf[c] = s_streams[streams[r][c]].horton;
+				if (out_topo)
+			topo_buf[c] = s_streams[streams[r][c]].topo_dim;
 		}		/* end if streams */
 	    }
 	    else if (out_zero) {
@@ -197,6 +209,9 @@
 		    hack_buf[c] = s_streams[streams[r][c]].hack;
 		if (out_horton)
 		    horton_buf[c] = s_streams[streams[r][c]].horton;
+		if (out_topo)
+		    topo_buf[c] = s_streams[streams[r][c]].topo_dim;
+		    
 	    }
 	}			/* end for cols */
 
@@ -208,6 +223,9 @@
 	    G_put_c_raster_row(out_hck_fd, hack_buf);
 	if (out_horton)
 	    G_put_c_raster_row(out_hrt_fd, horton_buf);
+	if (out_topo)
+	    G_put_c_raster_row(out_topo_fd, topo_buf);    
+	    
     }				/* end for r */
 
     G_percent(r, nrows, 2);
@@ -242,6 +260,14 @@
 	G_write_history(out_horton, &history);
 	G_message(_("%s Done!"),out_horton);
     }
+    
+  if (out_topo) {
+	G_close_cell(out_topo_fd);
+	G_free(topo_buf);
+	G_short_history(out_topo, "raster", &history);
+	G_write_history(out_topo, &history);
+	G_message(_("%s Done!"),out_topo);
+    }  
 
     return 0;
 
@@ -265,7 +291,8 @@
 	dbString table_name, db_sql, val_string;
 	char buf[1000];
 	char ins_prev_streams[50]; /* insert */
-	
+	char* cat_col_name="cat";	
+
 	/* table definition */
 	char* tab_cat_col_name="cat integer";
 	char* tab_stream="stream integer";
@@ -279,11 +306,13 @@
 	char* tab_cumlength="cum_length double precision";
 	char* tab_stright="stright double precision";
 	char* tab_fractal="fractal double precision";
+	char* tab_distance="out_dist double precision";
+	char* tab_topo_dim="topo_dim integer";
 	
 	G_message("Adding table...");
 	
 	/* trib num */
-	for (i=0;i<stream_num;++i) {
+	for (i=0;i<stream_num+1;++i) {
 			if (s_streams[i].trib_num>max_trib) 
 		max_trib=s_streams[i].trib_num;
 	}
@@ -301,9 +330,9 @@
 		else
 	sprintf(out_table, "%s_new",in_vector);
 	*/
-	
+
 	out_table=in_table;
-	
+
 	/* init */
 	
 	db_init_string(&db_sql);
@@ -340,7 +369,7 @@
 	}
 		
 	
-	sprintf(buf,"create table %s (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s)",
+	sprintf(buf,"create table %s (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s)",
 		out_table,
 		tab_cat_col_name,
 		tab_stream,
@@ -350,8 +379,10 @@
 		tab_horton,
 		tab_shreve,
 		tab_hack,
+		tab_topo_dim,
 		tab_length,
 		tab_cumlength,
+		tab_distance,
 		tab_stright,
 		tab_fractal);
 		
@@ -363,7 +394,7 @@
 		G_fatal_error("Cannot create table %s", db_get_string(&db_sql));
 	}
 	
-		if (db_create_index2(driver,out_table,tab_cat_col_name) !=DB_OK)
+		if (db_create_index2(driver,out_table,cat_col_name) !=DB_OK)
 	G_warning("cannot create index");
 	
 		if (db_grant_on_table(driver,out_table,DB_PRIV_SELECT, DB_GROUP |DB_PUBLIC) !=DB_OK)
@@ -371,7 +402,7 @@
 
 	db_begin_transaction(driver);
 
-	for (i=0;i<stream_num;++i) {
+	for (i=0;i<stream_num+1;++i) {
 		
 			if(s_streams[i].stream<0)
 		continue;
@@ -395,7 +426,7 @@
 	}
 	
 		sprintf(buf,"insert into %s  values	\
-			(%d, %d, %d, %s, %d, %d, %d, %d, %f, %f , %f, %f)",
+			(%d, %d, %d, %s, %d, %d, %d, %d, %d, %f, %f , %f, %f, %f)",
 		out_table,
 		i,
 		s_streams[i].stream,
@@ -405,8 +436,10 @@
 		s_streams[i].horton,
 		s_streams[i].shreeve,
 		s_streams[i].hack,
+		s_streams[i].topo_dim,
 		s_streams[i].length,
 		s_streams[i].accum,
+		s_streams[i].distance,
 		s_streams[i].stright,
 		s_streams[i].fractal);
 	

Modified: grass-addons/raster/r.stream.order/main.c
===================================================================
--- grass-addons/raster/r.stream.order/main.c	2010-02-09 08:30:16 UTC (rev 40886)
+++ grass-addons/raster/r.stream.order/main.c	2010-02-09 08:46:49 UTC (rev 40887)
@@ -25,7 +25,7 @@
 {
 
     struct GModule *module;	/* GRASS module for parsing arguments */
-    struct Option *in_dir_opt, *in_stm_opt, /* *in_vect_opt,*/ *in_table_opt, *in_acc_opt, *out_str_opt, *out_shr_opt, *out_hck_opt, *out_hrt_opt;	/* options */
+    struct Option *in_dir_opt, *in_stm_opt, /* *in_vect_opt,*/ *in_table_opt, *in_acc_opt, *out_str_opt, *out_shr_opt, *out_hck_opt, *out_hrt_opt, *out_topo_opt ;	/* options */
     struct Flag *out_back;	/* flags */
     int i;
 
@@ -69,6 +69,7 @@
     in_table_opt->key = "table";
     in_table_opt->type = TYPE_STRING;
     in_table_opt->required = NO;
+    in_table_opt->answer = NULL;
     in_table_opt->description =
 	"Name of new table to create";
     
@@ -120,6 +121,16 @@
     out_hck_opt->description =
 	"Name of Hack's main streams output map";
     out_hck_opt->guisection = _("Output options");
+    
+    out_topo_opt = G_define_option();
+    out_topo_opt->key = "topo";
+    out_topo_opt->type = TYPE_STRING;
+    out_topo_opt->required = NO;
+    out_topo_opt->answer = NULL;
+    out_topo_opt->gisprompt = "new,cell,raster";
+    out_topo_opt->description =
+	"Name of topological dimension output map";
+    out_topo_opt->guisection = _("Output options");
 
     /* Define flags */
     out_back = G_define_flag();
@@ -132,8 +143,8 @@
     G_get_window(&window);
 
     if (!out_str_opt->answer && !out_shr_opt->answer && !out_hck_opt->answer
-	&& !out_hrt_opt->answer)
-	G_fatal_error(_("You must select one or more output maps: strahler, horton, shreeve or hack"));
+	&& !out_hrt_opt->answer && !out_topo_opt->answer && !in_table_opt->answer)
+	G_fatal_error(_("You must select one or more output maps: strahler, horton, shreeve, hack, topo or insert the table name"));
     
 
     /* stores input options to variables */
@@ -148,6 +159,7 @@
     out_shreeve = out_shr_opt->answer;
     out_hack = out_hck_opt->answer;
     out_horton = out_hrt_opt->answer;
+    out_topo = out_topo_opt->answer;
     out_zero = (out_back->answer != 0);
 
 		if (out_strahler) {
@@ -176,7 +188,7 @@
 	stack_max = stream_num;	/* stack's size depends on number of streams */
     init_streams(stream_num);
     find_nodes(stream_num);
-    if (out_hack || out_horton || in_table)
+    if (out_hack || out_horton || in_table || out_topo)
    do_cum_length();
 
 	
@@ -186,7 +198,7 @@
 	shreeve();
     if (out_horton || in_table)
 	horton();
-    if (out_hack || in_table)
+    if (out_hack || out_topo || in_table)
 	hack();
   
   write_maps();

Modified: grass-addons/raster/r.stream.order/order.c
===================================================================
--- grass-addons/raster/r.stream.order/order.c	2010-02-09 08:30:16 UTC (rev 40886)
+++ grass-addons/raster/r.stream.order/order.c	2010-02-09 08:46:49 UTC (rev 40887)
@@ -31,7 +31,7 @@
 int init_streams(int stream_num)
 {
     int i;
-    s_streams = (STREAM *) G_malloc((stream_num + 1) * sizeof(STREAM));
+		s_streams = (STREAM *) G_malloc((stream_num + 1) * sizeof(STREAM));
         for (i = 0; i <= stream_num; ++i) {
 	s_streams[i].next_stream = -1;
 	s_streams[i].stream = -1;
@@ -44,6 +44,8 @@
 	s_streams[i].length = 0.;
 	s_streams[i].stright = 0.;
 	s_streams[i].fractal = 0.;
+	s_streams[i].distance=0.;
+	s_streams[i].topo_dim=0;
 	s_streams[i].trib[0] = 0;
 	s_streams[i].trib[1] = 0;
 	s_streams[i].trib[2] = 0;
@@ -234,6 +236,7 @@
 return 0;
 }
 
+
 /* 
    All algorithms used in analysis ar not recursive. For Strahler order and Shreve magnitude starts from initial channel and  proceed downstream. Algortitms try to assgin order for branch and if it is imposible start from next initial channel, till all branches are ordered.
    For Hortor and Hack ordering it proceed upstream and uses stack data structure to determine unordered branch. 
@@ -429,7 +432,7 @@
     return 0;
 }
 
-int hack(void)
+int hack(void) /* also calculate topological dimension */
 {
 
     int *stack;
@@ -437,17 +440,23 @@
     int cur_stream, cur_hack;
     double max_accum;
     int up_stream = 0;
+    
+    double cur_distance=0;
+    int cur_topo_dim=0;
 
-    G_message(_("Calculating Hack's main streams ..."));
+    G_message(_("Calculating Hack's main streams and topological dimension..."));
     stack = (int *)G_malloc(stack_max * sizeof(int));
 
     for (j = 0; j < outlets_num; ++j) {
+
 	cur_stream = s_streams[outlets[j]].stream;	/* outlet: init */
 	cur_hack = 1;
 	stack[0] = 0;
 	stack[1] = cur_stream;
 	top = 1;
-
+	
+	s_streams[cur_stream].topo_dim=top;
+	cur_distance=s_streams[cur_stream].distance=s_streams[cur_stream].length;
 	do {
 	    max_accum = 0;
 
@@ -473,21 +482,30 @@
 		}		/* end determining up_stream */
 
 		if (up_stream) {	/* at least one branch is not assigned */
+			
 		    if (s_streams[cur_stream].hack < 0) {
 			s_streams[cur_stream].hack = cur_hack;
-		    }
+				 }
 		    else {
 			cur_hack = s_streams[cur_stream].hack;
 			++cur_hack;
+			
 		    }
+
+				cur_distance=s_streams[cur_stream].distance;
+
 		    cur_stream = up_stream;
 		    stack[++top] = cur_stream;
-
+    
+				s_streams[cur_stream].distance=cur_distance+s_streams[cur_stream].length;
+				s_streams[cur_stream].topo_dim=top;
+			
+			
 		}
 		else {		/* all asigned, go downstream */
-
 		    cur_stream = stack[--top];
 
+
 		}		/* end up_stream */
 	    }			/* end spring/node */
 	} while (cur_stream);

Modified: grass-addons/raster/r.stream.order/orders.png
===================================================================
(Binary files differ)



More information about the grass-commit mailing list