[GRASS-SVN] r39697 - grass-addons/raster/r.stream.order
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Nov 9 09:39:19 EST 2009
Author: jarekj71
Date: 2009-11-09 09:39:19 -0500 (Mon, 09 Nov 2009)
New Revision: 39697
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
Log:
New wersion of r.stream.order
Modified: grass-addons/raster/r.stream.order/description.html
===================================================================
--- grass-addons/raster/r.stream.order/description.html 2009-11-09 14:33:09 UTC (rev 39696)
+++ grass-addons/raster/r.stream.order/description.html 2009-11-09 14:39:19 UTC (rev 39697)
@@ -15,7 +15,7 @@
and maps boundary may be differ but it may lead to unexpected results.</DD>
<p>
<DT><b>accum</b></DT>
-<DD>Flow accumulation (optional): name of flow accumulation file produced by r.watershed or r.stream.extract. This map is requied only if Horton's or Hack's ordering is performed. Flow accumulation map shall be of DCELL type, as is by default produced by r.watershed or r.stream.extract.</DD>
+<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>
@@ -55,7 +55,7 @@
Strahler's stream order has a good mathematical background. All catchments with streams in this context are directed graphs, oriented from the root towards the leaves. Equivalent definition of the Strahler number of a tree is that it is the height of the largest complete binary tree that can be homeomorphically embedded into the given tree; the Strahler number of a node in a tree is equivalent to the height of the largest complete binary tree that can be embedded below that node. The disadvantage of that methods is the lack of distinguishing a main channel which may interfere with the analytical process in highly elongated catchments
<H4>Horton's stream order</H4>
-Horton's stream order applies to the stream as a whole but not to segments or links since the order on any channel remains unchanged from source till it "dies" in the higher order stream or in the outlet of the catchment. The main segment of the catchment gets the order of the whole catchment, while its tributaries get the order of their own subcatchments. The main difficulties of the Horton's order are criteria to be considered to distinguish between "true" first order segments and extension of higher order segments. Thqat is the reason why Horton's ordering has rather historical sense and is substituted by the more unequivocal Strahler's ordering system. There are no natural algorithms to order stream network according to Horton' paradigm. The algorithm used in r.stream.order requires to first calculate Strahler's stream order (downstream) and next recalculate to Horton ordering (upstream). To make a decision about proper ordering it uses first Strahler ordering, and next, if both branches have the same orders it uses flow accumulation to choose the actual link. The algorithm starts with the outlet, where the outlet link is assigned the corresponding Strahler order. Next it goes upstream and determines links according to Strahler ordering. If the orders of tributaries differ, the algorithm proceeds with the channel of highest order, if all orders are the same, it chooses that one with higher flow accumulation rate (higher catchment area). When it reaches the initial channel it goes back to the last undetermined branch, assign its Strahler order as Horton order and goes upstream to the next initial links. In that way stream orders remain unchanged from the point where Horton's order have been determined to the source.
+Horton's stream order applies to the stream as a whole but not to segments or links since the order on any channel remains unchanged from source till it "dies" in the higher order stream or in the outlet of the catchment. The main segment of the catchment gets the order of the whole catchment, while its tributaries get the order of their own subcatchments. The main difficulties of the Horton's order are criteria to be considered to distinguish between "true" first order segments and extension of higher order segments. Thqat is the reason why Horton's ordering has rather historical sense and is substituted by the more unequivocal Strahler's ordering system. There are no natural algorithms to order stream network according to Horton' paradigm. The algorithm used in r.stream.order requires to first calculate Strahler's stream order (downstream) and next recalculate to Horton ordering (upstream). To make a decision about proper ordering it uses first Strahler ordering, and next, if both branches have the same orders it uses flow accumulation to choose the actual link. The algorithm starts with the outlet, where the outlet link is assigned the corresponding Strahler order. Next it goes upstream and determines links according to Strahler ordering. If the orders of tributaries differ, the algorithm proceeds with the channel of highest order, if all orders are the same, it chooses that one with higher flow length rate or higher catchment area if accumulation is used. When it reaches the initial channel it goes back to the last undetermined branch, assign its Strahler order as Horton order and goes upstream to the next initial links. In that way stream orders remain unchanged from the point where Horton's order have been determined to the source.
<BR>
<B>Advantages and disadvantages of Horton's ordering:</B>
@@ -68,7 +68,7 @@
The algorithm is very similar to Strahler's algorithm, it proceeds downstream, and at every node the stream magnitude is the sum of its tributaries.
<P>
<H4>Hack's main streams order</H4>
-This method of ordering calculates main streams of main catchment and every subcatchments. Main stream of every catchment is set to 1, and consequently all its tributaries receive order 2. Their tributaries receive order 3 etc. The order of every stream remains constant up to its initial link. The route of every main stream is determined according to the maximum flow accumulation value of particular streams. So the main stream of every subcatchment is the stream with high accumulation. In most cases the main stream is the longest watercourse of the catchment, but in some cases, when a catchment consists of both rounded and elongated subcatchments these rules may not be maintained. The algorithm assigns 1 to every outlets stream and goes upstream according to maximum flow accumulation of every branch. When it reaches an initial stream it step back to the first unassigned confluence. It assigns order 2 to unordered tributaries and again goes upstream to the next initial stream. The process runs until all branches of all outlets are ordered.
+This method of ordering calculates main streams of main catchment and every subcatchments. Main stream of every catchment is set to 1, and consequently all its tributaries receive order 2. Their tributaries receive order 3 etc. The order of every stream remains constant up to its initial link. The route of every main stream is determined according to the maximum flow length value of particular streams. So the main stream of every subcatchment is the longest stream or strean with highest accumulation rate if accumulation map is used. In most cases the main stream is the longest watercourse of the catchment, but in some cases, when a catchment consists of both rounded and elongated subcatchments these rules may not be maintained. The algorithm assigns 1 to every outlets stream and goes upstream according to maximum flow accumulation of every branch. When it reaches an initial stream it step back to the first unassigned confluence. It assigns order 2 to unordered tributaries and again goes upstream to the next initial stream. The process runs until all branches of all outlets are ordered.
<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.
Modified: grass-addons/raster/r.stream.order/global.h
===================================================================
--- grass-addons/raster/r.stream.order/global.h 2009-11-09 14:33:09 UTC (rev 39696)
+++ grass-addons/raster/r.stream.order/global.h 2009-11-09 14:39:19 UTC (rev 39697)
@@ -44,8 +44,14 @@
/* int *trib; */
double accum;
+ /* float length; */
};
+#define INITS struct init
+INITS {
+ int r;
+ int c;
+};
/* functions.c */
@@ -53,6 +59,7 @@
int trib_nums(int r, int c);
int init_streams(int stream_num);
int find_nodes(int stream_num);
+int do_cum_length (void);
int strahler(void);
int horton(void);
int shreeve(void);
@@ -82,5 +89,6 @@
GLOBAL int nrows, ncols;
GLOBAL STREAM *s_streams; /* stream structure all parameters we have here */
+GLOBAL INITS *s_inits;
int *springs, *outlets;
int springs_num, outlets_num;
Modified: grass-addons/raster/r.stream.order/io.c
===================================================================
--- grass-addons/raster/r.stream.order/io.c 2009-11-09 14:33:09 UTC (rev 39696)
+++ grass-addons/raster/r.stream.order/io.c 2009-11-09 14:39:19 UTC (rev 39697)
@@ -116,6 +116,7 @@
int stream_number(void)
{
+
char *cur_mapset;
CELL c_min, c_max;
struct Range stream_range;
@@ -124,6 +125,18 @@
G_read_range(in_streams,cur_mapset,&stream_range);
G_get_range_min_max(&stream_range, &c_min, &c_max);
return c_max;
+/*
+
+int r, c;
+int max=0;
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+ if (streams[r][c]>max)
+ max=streams[r][c];
+ }
+ }
+ return max;
+*/
}
int write_maps(void)
Modified: grass-addons/raster/r.stream.order/main.c
===================================================================
--- grass-addons/raster/r.stream.order/main.c 2009-11-09 14:33:09 UTC (rev 39696)
+++ grass-addons/raster/r.stream.order/main.c 2009-11-09 14:39:19 UTC (rev 39697)
@@ -64,7 +64,7 @@
in_acc_opt->answer = NULL;
in_acc_opt->gisprompt = "old,cell,raster";
in_acc_opt->description =
- "Name of accumulation file (r.watershed or r.stream.extract)";
+ "(Not recommended) Name of accumulation file (r.watershed or r.stream.extract)";
/* output option - at least one is reqired */
@@ -93,7 +93,7 @@
out_hrt_opt->answer = NULL;
out_hrt_opt->gisprompt = "new,cell,raster";
out_hrt_opt->description =
- "Name of Horton's stream order output map (require accum file)";
+ "Name of Horton's stream order output map";
out_hrt_opt->guisection = _("Output options");
out_hck_opt = G_define_option();
@@ -103,7 +103,7 @@
out_hck_opt->answer = NULL;
out_hck_opt->gisprompt = "new,cell,raster";
out_hck_opt->description =
- "Name of Hack's main streams output map (require accum file)";
+ "Name of Hack's main streams output map";
out_hck_opt->guisection = _("Output options");
/* Define flags */
@@ -119,8 +119,7 @@
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"));
- if (out_hrt_opt->answer && !in_acc_opt->answer)
- G_fatal_error(_("If you select horton or hack, accum is neccessary"));
+
/* stores input options to variables */
in_dirs = in_dir_opt->answer;
@@ -155,10 +154,13 @@
ncols = G_window_cols();
create_base_maps();
- stream_num = stream_number();
- stack_max = stream_num; /* stack's size depends on number of streams */
+ stream_num = stream_number();
+
+ 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_accum)
+ do_cum_length();
if (out_strahler || out_horton)
strahler();
if (out_shreeve)
Modified: grass-addons/raster/r.stream.order/order.c
===================================================================
--- grass-addons/raster/r.stream.order/order.c 2009-11-09 14:33:09 UTC (rev 39696)
+++ grass-addons/raster/r.stream.order/order.c 2009-11-09 14:39:19 UTC (rev 39697)
@@ -24,14 +24,14 @@
G_fatal_error(_("Error finding nodes. Stream and direction maps probably do not match..."));
if (trib > 3)
G_warning(_("Stream network may be too dense..."));
- return trib;
+
+ return trib;
} /* end trib_num */
int init_streams(int stream_num)
{
int i;
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;
@@ -39,7 +39,8 @@
s_streams[i].shreeve = -1;
s_streams[i].hack = -1;
s_streams[i].horton = -1;
- s_streams[i].accum = 0.0;
+ s_streams[i].accum = 0.;
+ /* s_streams[i].length = 0.0; */
s_streams[i].trib[0] = 0;
s_streams[i].trib[1] = 0;
s_streams[i].trib[2] = 0;
@@ -65,6 +66,7 @@
outlets = (int *)G_malloc((stream_num) * sizeof(int));
springs = (int *)G_malloc((stream_num) * sizeof(int));
+ s_inits = (INITS *)G_malloc((stream_num) *sizeof(INITS));
for (r = 0; r < nrows; ++r) {
for (c = 0; c < ncols; ++c) {
@@ -102,23 +104,27 @@
G_fatal_error(_("Error finding nodes. Stream and direction maps probably do not match..."));
s_streams[cur_stream].trib_num = 0;
+ s_inits[springs_num].r=r;
+ s_inits[springs_num].c=c;
springs[springs_num++] = cur_stream; /* collecting springs */
+
}
if (trib_num > 1) { /* adding tributuaries */
s_streams[cur_stream].trib_num = trib_num;
+
for (i = 1; i < 9; ++i) {
-
if (trib > 4)
G_fatal_error(_("Error finding nodes. Stream and direction maps probably do not match..."));
if (r + nextr[i] < 0 || r + nextr[i] > (nrows - 1) ||
c + nextc[i] < 0 || c + nextc[i] > (ncols - 1))
continue;
+
j = (i + 4) > 8 ? i - 4 : i + 4;
if (streams[r + nextr[i]][c + nextc[i]] > 0 &&
dirs[r + nextr[i]][c + nextc[i]] == j) {
- if (in_accum) { /* only if horton and hack */
+ if (in_accum) { /* only if accum map is selected */
s_streams[streams[r + nextr[i]]
[c + nextc[i]]].accum =
fabs(accum[r + nextr[i]][c + nextc[i]]);
@@ -135,6 +141,76 @@
return 0;
}
+int do_cum_length (void) {
+
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ int i, s, d; /* s - streams index, d - direction */
+ int done = 1;
+ int r, c;
+ int next_r, next_c;
+ int cur_stream, next_stream;
+ float cur_northing, cur_easting;
+ float next_northing, next_easting;
+ double cur_length=0.;
+
+ G_begin_distance_calculations();
+
+ for (s=0; s<springs_num; ++s) { /* main loop on springs */
+ r = s_inits[s].r;
+ c = s_inits[s].c;
+ cur_stream=streams[r][c];
+ cur_length=0;
+ done=1;
+
+ while(streams && done) {
+
+ cur_northing = window.north - (r + .5) * window.ns_res;
+ cur_easting = window.west + (c + .5) * window.ew_res;
+ d=dirs[r][c];
+ if (d<1 || /* end of route: sink, border or something else */
+ r + nextr[d] < 0 || r + nextr[d] > (nrows - 1) ||
+ c + nextc[d] < 0 || c + nextc[d] > (ncols - 1)) {
+
+ next_northing = window.north - (r + .5) * window.ns_res;
+ next_easting = window.west + (c + .5) * window.ew_res;
+ cur_length=G_distance(next_easting, next_northing, cur_easting, cur_northing);
+ s_streams[cur_stream].accum += cur_length;
+ break;
+ }
+
+ next_r=r+nextr[dirs[r][c]];
+ next_c=c+nextc[dirs[r][c]];
+ next_northing = window.north - (next_r + .5) * window.ns_res;
+ next_easting = window.west + (next_c + .5) * window.ew_res;
+ cur_length = G_distance(next_easting, next_northing, cur_easting, cur_northing);
+ s_streams[cur_stream].accum += cur_length;
+ r=next_r;
+ c=next_c;
+
+ if (streams[next_r][next_c] != cur_stream) {
+ next_stream=streams[next_r][next_c];
+
+ if (s_streams[next_stream].accum<s_streams[cur_stream].accum)
+ s_streams[next_stream].accum=s_streams[cur_stream].accum;
+
+ for (i=0; i<s_streams[next_stream].trib_num;++i){
+ if (s_streams[s_streams[next_stream].trib[i]].accum==0)
+ done=0;
+ }
+
+ cur_stream=next_stream;
+ }/* end if */
+ } /* end while */
+ } /* end for */
+
+/* for (i=1;i<stream_num+1;++i) {
+ G_message(_("%d: %f"),i,s_streams[i].length);
+} */
+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.
More information about the grass-commit
mailing list