[GRASS-SVN] r39513 - grass-addons/raster/r.stream.extract
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Oct 13 14:19:38 EDT 2009
Author: mmetz
Date: 2009-10-13 14:19:38 -0400 (Tue, 13 Oct 2009)
New Revision: 39513
Added:
grass-addons/raster/r.stream.extract/del_streams.c
Modified:
grass-addons/raster/r.stream.extract/close.c
grass-addons/raster/r.stream.extract/description.html
grass-addons/raster/r.stream.extract/do_astar.c
grass-addons/raster/r.stream.extract/local_proto.h
grass-addons/raster/r.stream.extract/main.c
grass-addons/raster/r.stream.extract/streams.c
grass-addons/raster/r.stream.extract/thin.c
Log:
optimized A* Search, added option to delete short stream segments
Modified: grass-addons/raster/r.stream.extract/close.c
===================================================================
--- grass-addons/raster/r.stream.extract/close.c 2009-10-13 17:15:48 UTC (rev 39512)
+++ grass-addons/raster/r.stream.extract/close.c 2009-10-13 18:19:38 UTC (rev 39513)
@@ -56,6 +56,9 @@
thisindex = INDEX(r, c);
stream_id = stream[thisindex];
+ if (!stream_id)
+ continue;
+
Vect_reset_line(Points);
Vect_reset_cats(Cats);
@@ -194,7 +197,7 @@
sprintf(buf, "insert into %s values ( %d, \'%s\', %d )",
Fi->table, i,
- (stream_node[i].n_trib > 0 ? "start" : "intermediate"),
+ (stream_node[i].n_trib > 0 ? "intermediate" : "start"),
(stream_node[i].n_trib > 0));
db_set_string(&dbsql, buf);
Added: grass-addons/raster/r.stream.extract/del_streams.c
===================================================================
--- grass-addons/raster/r.stream.extract/del_streams.c (rev 0)
+++ grass-addons/raster/r.stream.extract/del_streams.c 2009-10-13 18:19:38 UTC (rev 39513)
@@ -0,0 +1,225 @@
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+/* delete short stream segments according to threshold */
+int del_streams(int min_length)
+{
+ int i;
+ int n_deleted = 0;
+ unsigned int thisindex, next_stream_pos = -1;
+ int curr_stream, stream_id, other_trib, tmp_trib;
+ int slength;
+
+ G_message(_("delete stream segments shorter than %d cells..."), min_length);
+
+ /* go through all nodes */
+ for (i = 1; i <= n_stream_nodes; i++) {
+ G_percent(i, n_stream_nodes, 2);
+
+ /* not a stream head */
+ if (stream_node[i].n_trib > 0)
+ continue;
+
+ /* already deleted */
+ thisindex = INDEX(stream_node[i].r, stream_node[i].c);
+ if (stream[thisindex] == 0)
+ continue;
+
+ /* get length counted as n cells */
+ if ((slength = seg_length(i, &next_stream_pos)) >= min_length)
+ continue;
+
+ stream_id = i;
+
+ /* check n sibling tributaries */
+ if ((curr_stream = stream[next_stream_pos]) != stream_id) {
+ /* only one sibling tributary */
+ if (stream_node[curr_stream].n_trib == 2) {
+ if (stream_node[curr_stream].trib[0] != stream_id)
+ other_trib = stream_node[curr_stream].trib[0];
+ else
+ other_trib = stream_node[curr_stream].trib[1];
+
+ /* other trib is also stream head */
+ if (stream_node[other_trib].n_trib == 0) {
+ /* use shorter one */
+ if (seg_length(other_trib, NULL) < slength) {
+ tmp_trib = stream_id;
+ stream_id = other_trib;
+ other_trib = tmp_trib;
+ }
+ }
+ del_stream_seg(stream_id);
+ n_deleted++;
+
+ /* update downstream IDs */
+ update_stream_id(curr_stream, other_trib);
+ }
+ /* more than one other tributary */
+ else {
+ del_stream_seg(stream_id);
+ n_deleted++;
+ }
+ }
+ /* stream head is also outlet */
+ else {
+ del_stream_seg(stream_id);
+ n_deleted++;
+ }
+ }
+
+ G_debug(1, "%d stream segments deleted", n_deleted);
+
+ return n_deleted;
+}
+
+/* get stream segment length */
+int seg_length(int stream_id, unsigned int *new_stream_pos)
+{
+ int r, c, r_nbr, c_nbr;
+ int slength = 1;
+ struct ddir draindir, *founddir;
+ int curr_stream;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ r = stream_node[stream_id].r;
+ c = stream_node[stream_id].c;
+
+ /* get next downstream point */
+ draindir.pos = INDEX(r, c);
+ if (new_stream_pos)
+ *new_stream_pos = draindir.pos;
+ while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
+ r_nbr = r + asp_r[(int)founddir->dir];
+ c_nbr = c + asp_c[(int)founddir->dir];
+
+ /* user-defined depression */
+ if (r_nbr == r && c_nbr == c)
+ break;
+ /* outside region */
+ if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
+ break;
+ /* next stream */
+ if ((curr_stream = stream[INDEX(r_nbr, c_nbr)]) != stream_id) {
+ if (new_stream_pos)
+ *new_stream_pos = INDEX(r_nbr, c_nbr);
+ break;
+ }
+ slength++;
+ r = r_nbr;
+ c = c_nbr;
+ if (new_stream_pos)
+ *new_stream_pos = draindir.pos;
+ draindir.pos = INDEX(r, c);
+ }
+
+ return slength;
+}
+
+/* delete stream segment */
+int del_stream_seg(int stream_id)
+{
+ int i, r, c, r_nbr, c_nbr;
+ unsigned int this_index;
+ struct ddir draindir, *founddir;
+ int curr_stream;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ r = stream_node[stream_id].r;
+ c = stream_node[stream_id].c;
+ this_index = INDEX(r, c);
+ stream[this_index] = 0;
+ curr_stream = stream_id;
+
+ /* get next downstream point */
+ draindir.pos = this_index;
+ while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
+ r_nbr = r + asp_r[(int)founddir->dir];
+ c_nbr = c + asp_c[(int)founddir->dir];
+
+ /* user-defined depression */
+ if (r_nbr == r && c_nbr == c)
+ break;
+ /* outside region */
+ if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
+ break;
+ /* next stream */
+ if ((curr_stream = stream[INDEX(r_nbr, c_nbr)]) != stream_id)
+ break;
+ r = r_nbr;
+ c = c_nbr;
+ this_index = INDEX(r, c);
+ stream[this_index] = 0;
+ rbtree_remove(draintree, &draindir);
+ draindir.pos = this_index;
+ }
+
+ /* update tributaries */
+ if (curr_stream != stream_id) {
+ for (i = 0; i < stream_node[curr_stream].n_trib; i++) {
+ if (stream_node[curr_stream].trib[i] == stream_id) {
+ stream_node[curr_stream].trib[i] = stream_node[curr_stream].trib[stream_node[curr_stream].n_trib - 1];
+ stream_node[curr_stream].n_trib--;
+ break;
+ }
+ }
+ }
+
+ return 1;
+}
+
+/* update downstream id */
+int update_stream_id(int stream_id, int new_stream_id)
+{
+ int i, r, c, r_nbr, c_nbr;
+ unsigned int this_index;
+ struct ddir draindir, *founddir;
+ int curr_stream;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ r = stream_node[stream_id].r;
+ c = stream_node[stream_id].c;
+ this_index = INDEX(r, c);
+ stream[this_index] = new_stream_id;
+ curr_stream = stream_id;
+
+ /* get next downstream point */
+ draindir.pos = this_index;
+ while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
+ r_nbr = r + asp_r[(int)founddir->dir];
+ c_nbr = c + asp_c[(int)founddir->dir];
+
+ /* user-defined depression */
+ if (r_nbr == r && c_nbr == c)
+ break;
+ /* outside region */
+ if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
+ break;
+ /* next stream */
+ if ((curr_stream = stream[INDEX(r_nbr, c_nbr)]) != stream_id)
+ break;
+ r = r_nbr;
+ c = c_nbr;
+ this_index = INDEX(r, c);
+ stream[this_index] = new_stream_id;
+ draindir.pos = this_index;
+ }
+
+ /* update tributaries */
+ if (curr_stream != stream_id) {
+ for (i = 0; i < stream_node[curr_stream].n_trib; i++) {
+ if (stream_node[curr_stream].trib[i] == stream_id) {
+ stream_node[curr_stream].trib[i] = new_stream_id;
+ break;
+ }
+ }
+ }
+
+ return curr_stream;
+}
Modified: grass-addons/raster/r.stream.extract/description.html
===================================================================
--- grass-addons/raster/r.stream.extract/description.html 2009-10-13 17:15:48 UTC (rev 39512)
+++ grass-addons/raster/r.stream.extract/description.html 2009-10-13 18:19:38 UTC (rev 39513)
@@ -65,6 +65,11 @@
Larger <em>mexp</em> values decrease the number of streams in flat areas
and increase the number of streams in steep areas. If <em>weight</em> is
given, the weight is applied first.
+<p>
+<dt><em>stream_length</em>
+<dd>Minimum stream length in number of cells for first-order (head/spring)
+stream segments. All first-order stream segments shorter than
+<em>stream_length</em> will be deleted.
<p>
<dt><em>stream_rast</em>
@@ -148,4 +153,4 @@
<h2>AUTHOR</h2>
Markus Metz
-<p><i>Last changed: $Date: 2008-05-16 21:09:06 +0200 (Fri, 16 May 2008) $</i>
+<p><i>Last changed: $Date$</i>
Property changes on: grass-addons/raster/r.stream.extract/description.html
___________________________________________________________________
Added: svn:keywords
+ Date
Modified: grass-addons/raster/r.stream.extract/do_astar.c
===================================================================
--- grass-addons/raster/r.stream.extract/do_astar.c 2009-10-13 17:15:48 UTC (rev 39512)
+++ grass-addons/raster/r.stream.extract/do_astar.c 2009-10-13 18:19:38 UTC (rev 39513)
@@ -50,10 +50,10 @@
astar_pts[first_cum] = astp;
first_cum--;
- r = astp.r;
- c = astp.c;
+ thisindex = astp.idx;
+ r = thisindex / ncols;
+ c = thisindex - r * ncols;
- thisindex = INDEX(r, c);
ele_val = ele[thisindex];
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -103,7 +103,7 @@
int sift_up(unsigned int start, CELL elec)
{
- unsigned int child, child_added, parent, nindex;
+ unsigned int child, child_added, parent;
CELL elep;
struct ast_point childp;
@@ -114,10 +114,8 @@
while (child > 1) {
parent = get_parent(child);
- nindex = INDEX(astar_pts[parent].r, astar_pts[parent].c);
+ elep = ele[astar_pts[parent].idx];
- elep = ele[nindex];
-
/* child < parent */
if (heap_cmp(elec, child_added, elep, astar_added[parent]) == 1) {
/* push parent point down */
@@ -156,8 +154,7 @@
G_fatal_error(_("heapsize too large"));
astar_added[heap_size] = nxt_avail_pt;
- astar_pts[heap_size].r = r;
- astar_pts[heap_size].c = c;
+ astar_pts[heap_size].idx = INDEX(r, c);
astar_pts[heap_size].asp = asp;
nxt_avail_pt++;
@@ -188,13 +185,13 @@
parent = 1;
while ((child = get_child(parent)) <= heap_size) {
- elec = ele[INDEX(astar_pts[child].r, astar_pts[child].c)];
+ elec = ele[astar_pts[child].idx];
if (child < heap_size) {
childr = child + 1;
- i = child + 3; /* change the number, GET_CHILD() and GET_PARENT() to play with different d-ary heaps */
+ i = child + 3;
while (childr <= heap_size && childr < i) {
- eler = ele[INDEX(astar_pts[childr].r, astar_pts[childr].c)];
+ eler = ele[astar_pts[childr].idx];
if (heap_cmp
(eler, astar_added[childr], elec,
@@ -217,7 +214,7 @@
astar_added[parent] = astar_added[heap_size];
astar_pts[parent] = astar_pts[heap_size];
- elec = ele[INDEX(astar_pts[parent].r, astar_pts[parent].c)];
+ elec = ele[astar_pts[parent].idx];
/* sift up last swapped point, only necessary if hole moved to heap end */
sift_up(parent, elec);
}
Modified: grass-addons/raster/r.stream.extract/local_proto.h
===================================================================
--- grass-addons/raster/r.stream.extract/local_proto.h 2009-10-13 17:15:48 UTC (rev 39512)
+++ grass-addons/raster/r.stream.extract/local_proto.h 2009-10-13 18:19:38 UTC (rev 39513)
@@ -16,7 +16,7 @@
struct ast_point
{
- int r, c;
+ int idx;
char asp;
};
@@ -62,12 +62,18 @@
int do_astar(void);
unsigned int heap_add(int, int, CELL, char);
+/* streams.c */
+int do_accum(double);
+int extract_streams(double, double, int, int);
+
/* thin.c */
int thin_streams(void);
-/* streams.c */
-int do_accum(double);
-int extract_streams(double, double, int);
+/* del_streams.c */
+int del_streams(int);
+int seg_length(int, unsigned int *);
+int del_stream_seg(int);
+int update_stream_id(int, int);
/* close.c */
int close_maps(char *, char *, char *);
Modified: grass-addons/raster/r.stream.extract/main.c
===================================================================
--- grass-addons/raster/r.stream.extract/main.c 2009-10-13 17:15:48 UTC (rev 39512)
+++ grass-addons/raster/r.stream.extract/main.c 2009-10-13 18:19:38 UTC (rev 39513)
@@ -30,6 +30,7 @@
struct Option *ele, *acc, *weight;
struct Option *threshold, *d8cut;
struct Option *mont_exp;
+ struct Option *min_stream_length;
} input;
struct
{
@@ -40,6 +41,7 @@
struct GModule *module;
int ele_fd, acc_fd, weight_fd;
double threshold, d8cut, mont_exp;
+ int min_stream_length = 0;
char *mapset;
G_gisinit(argv[0]);
@@ -95,6 +97,16 @@
input.mont_exp->description =
_("Montgomery: accumulation is multiplied with pow(slope,mexp) and then compared with threshold.");
+ input.min_stream_length = G_define_option();
+ input.min_stream_length->key = "stream_length";
+ input.min_stream_length->type = TYPE_INTEGER;
+ input.min_stream_length->required = NO;
+ input.min_stream_length->answer = "0";
+ input.min_stream_length->label =
+ _("Delete stream segments shorter than stream_length cells.");
+ input.min_stream_length->description =
+ _("Applies only to first-order stream segments (springs/stream heads).");
+
output.stream_rast = G_define_standard_option(G_OPT_R_OUTPUT);
output.stream_rast->key = "stream_rast";
output.stream_rast->description =
@@ -121,7 +133,6 @@
/***********************/
/* check options */
-
/***********************/
/* input maps exist ? */
@@ -168,6 +179,16 @@
else
mont_exp = 0;
+ /* Minimum stream segment length */
+ if (input.min_stream_length->answer) {
+ min_stream_length = atoi(input.min_stream_length->answer);
+ if (min_stream_length < 0)
+ G_fatal_error(_("Minimum stream length must be positive or zero but is %d"),
+ min_stream_length);
+ }
+ else
+ min_stream_length = 0;
+
/* Check for some output map */
if ((output.stream_rast->answer == NULL)
&& (output.stream_vect->answer == NULL)) {
@@ -176,7 +197,6 @@
/*********************/
/* preparation */
-
/*********************/
/* open input maps */
@@ -213,6 +233,7 @@
/* allocate memory */
ele = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
+ /* TODO: allocate acc and stream memory only after A* Search */
acc = (DCELL *) G_malloc(nrows * ncols * sizeof(DCELL));
stream = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
if (input.weight->answer)
@@ -226,7 +247,6 @@
/********************/
/* processing */
-
/********************/
/* sort elevation and get initial stream direction */
@@ -237,16 +257,15 @@
if (acc_fd < 0) {
if (do_accum(d8cut) < 0)
G_fatal_error(_("could not calculate flow accumulation"));
- if (extract_streams
- (threshold, mont_exp, (input.weight->answer != NULL)) < 0)
- G_fatal_error(_("could not extract streams"));
}
else {
- if (extract_streams
- (threshold, mont_exp, (input.weight->answer != NULL)) < 0)
- G_fatal_error(_("could not extract streams"));
+ /* load accumulation */
}
+ if (extract_streams
+ (threshold, mont_exp, (input.weight->answer != NULL), min_stream_length) < 0)
+ G_fatal_error(_("could not extract streams"));
+
G_free(ele);
G_free(acc);
if (input.weight->answer)
@@ -256,6 +275,12 @@
if (thin_streams() < 0)
G_fatal_error(_("could not extract streams"));
+ /* delete short streams */
+ if (min_stream_length) {
+ if (del_streams(min_stream_length) < 0)
+ G_fatal_error(_("could not extract streams"));
+ }
+
/* write output maps */
if (close_maps(output.stream_rast->answer, output.stream_vect->answer,
output.dir_rast->answer) < 0)
Modified: grass-addons/raster/r.stream.extract/streams.c
===================================================================
--- grass-addons/raster/r.stream.extract/streams.c 2009-10-13 17:15:48 UTC (rev 39512)
+++ grass-addons/raster/r.stream.extract/streams.c 2009-10-13 18:19:38 UTC (rev 39513)
@@ -38,8 +38,8 @@
return result;
}
-int continue_stream(CELL is_swale, int r, int c, int r_max, int c_max,
- unsigned int thisindex, int *stream_no)
+int continue_stream(CELL stream_id, int r, int c, int r_max, int c_max,
+ unsigned int thisindex, int *stream_no, int min_length)
{
char aspect;
int curr_stream;
@@ -68,134 +68,142 @@
if (curr_stream <= 0) {
/* no confluence, just continue */
G_debug(2, "no confluence, just continue stream");
- stream[INDEX(r_max, c_max)] = is_swale;
+ stream[INDEX(r_max, c_max)] = stream_id;
+ return 0;
}
- else {
- G_debug(2, "confluence");
- /* new confluence */
- if (stream_node[curr_stream].r != r_max ||
- stream_node[curr_stream].c != c_max) {
- G_debug(2, "new confluence");
- /* set new stream id */
- curr_stream = stream[INDEX(r_max, c_max)] = ++(*stream_no);
- /* add stream node */
- if (*stream_no >= n_alloc_nodes - 1) {
- n_alloc_nodes += stream_node_step;
- stream_node =
- (struct snode *)G_realloc(stream_node,
- n_alloc_nodes *
- sizeof(struct snode));
- }
- stream_node[*stream_no].r = r_max;
- stream_node[*stream_no].c = c_max;
- stream_node[*stream_no].id = *stream_no;
- stream_node[*stream_no].n_trib = 0;
- stream_node[*stream_no].n_trib_total = 0;
- stream_node[*stream_no].n_alloc = 0;
- stream_node[*stream_no].trib = NULL;
- stream_node[*stream_no].acc = NULL;
- n_stream_nodes++;
- /* debug */
- if (n_stream_nodes != *stream_no)
- G_warning(_("stream_no %d and n_stream_nodes %d out of sync"),
- *stream_no, n_stream_nodes);
+ G_debug(2, "confluence");
+ /* delete short stream segments */
+ /* if (min_length && stream_node[stream_id].n_trib == 0) {
+ if (seg_length(stream_id, NULL) < min_length) {
+ del_stream_seg(stream_id);
+ return 0;
+ }
+ }
+ */
+
+ /* new confluence */
+ if (stream_node[curr_stream].r != r_max ||
+ stream_node[curr_stream].c != c_max) {
+ G_debug(2, "new confluence");
+ /* set new stream id */
+ curr_stream = stream[INDEX(r_max, c_max)] = ++(*stream_no);
+ /* add stream node */
+ if (*stream_no >= n_alloc_nodes - 1) {
+ n_alloc_nodes += stream_node_step;
+ stream_node =
+ (struct snode *)G_realloc(stream_node,
+ n_alloc_nodes *
+ sizeof(struct snode));
+ }
+ stream_node[*stream_no].r = r_max;
+ stream_node[*stream_no].c = c_max;
+ stream_node[*stream_no].id = *stream_no;
+ stream_node[*stream_no].n_trib = 0;
+ stream_node[*stream_no].n_trib_total = 0;
+ stream_node[*stream_no].n_alloc = 0;
+ stream_node[*stream_no].trib = NULL;
+ stream_node[*stream_no].acc = NULL;
+ n_stream_nodes++;
- /* add all tributaries */
- G_debug(2, "add all tributaries");
- for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r_nbr, c_nbr for neighbours */
- r_nbr = r_max + nextdr[ct_dir];
- c_nbr = c_max + nextdc[ct_dir];
- /* check that neighbour is within region */
- if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
- c_nbr < ncols) {
- draindir.pos = INDEX(r_nbr, c_nbr);
- if ((founddir =
- rbtree_find(draintree, &draindir)) != NULL) {
- if (r_nbr + asp_r[(int)founddir->dir] == r_max &&
- c_nbr + asp_c[(int)founddir->dir] == c_max) {
+ /* debug */
+ if (n_stream_nodes != *stream_no)
+ G_warning(_("stream_no %d and n_stream_nodes %d out of sync"),
+ *stream_no, n_stream_nodes);
- /* add tributary to stream node */
- if (stream_node[curr_stream].n_trib >=
- stream_node[curr_stream].n_alloc) {
- size_t new_size;
+ /* add all tributaries */
+ G_debug(2, "add all tributaries");
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r_nbr, c_nbr for neighbours */
+ r_nbr = r_max + nextdr[ct_dir];
+ c_nbr = c_max + nextdc[ct_dir];
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
+ c_nbr < ncols) {
+ draindir.pos = INDEX(r_nbr, c_nbr);
+ if ((founddir =
+ rbtree_find(draintree, &draindir)) != NULL) {
+ if (r_nbr + asp_r[(int)founddir->dir] == r_max &&
+ c_nbr + asp_c[(int)founddir->dir] == c_max) {
- stream_node[curr_stream].n_alloc += 2;
- new_size =
- stream_node[curr_stream].n_alloc *
- sizeof(int);
- stream_node[curr_stream].trib =
- (int *)G_realloc(stream_node[curr_stream].
- trib, new_size);
- new_size =
- stream_node[curr_stream].n_alloc *
- sizeof(double);
- stream_node[curr_stream].acc =
- (double *)
- G_realloc(stream_node[curr_stream].acc,
- new_size);
- }
+ /* add tributary to stream node */
+ if (stream_node[curr_stream].n_trib >=
+ stream_node[curr_stream].n_alloc) {
+ size_t new_size;
- stream_node[curr_stream].
- trib[stream_node[curr_stream].n_trib] =
- stream[draindir.pos];
- stream_node[curr_stream].
- acc[stream_node[curr_stream].n_trib++] =
- acc[draindir.pos];
+ stream_node[curr_stream].n_alloc += 2;
+ new_size =
+ stream_node[curr_stream].n_alloc *
+ sizeof(int);
+ stream_node[curr_stream].trib =
+ (int *)G_realloc(stream_node[curr_stream].
+ trib, new_size);
+ new_size =
+ stream_node[curr_stream].n_alloc *
+ sizeof(double);
+ stream_node[curr_stream].acc =
+ (double *)
+ G_realloc(stream_node[curr_stream].acc,
+ new_size);
}
+
+ stream_node[curr_stream].
+ trib[stream_node[curr_stream].n_trib] =
+ stream[draindir.pos];
+ stream_node[curr_stream].
+ acc[stream_node[curr_stream].n_trib++] =
+ acc[draindir.pos];
}
}
}
+ }
- /* update stream IDs downstream */
- G_debug(2, "update stream IDs downstream");
- r_nbr = r_max;
- c_nbr = c_max;
+ /* update stream IDs downstream */
+ G_debug(2, "update stream IDs downstream");
+ r_nbr = r_max;
+ c_nbr = c_max;
+ draindir.pos = INDEX(r_nbr, c_nbr);
+
+ while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
+ if (asp_r[(int)founddir->dir] == 0 &&
+ asp_c[(int)founddir->dir] == 0)
+ G_fatal_error("no valid stream direction");
+ r_nbr = r_nbr + asp_r[(int)founddir->dir];
+ c_nbr = c_nbr + asp_c[(int)founddir->dir];
draindir.pos = INDEX(r_nbr, c_nbr);
-
- while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
- if (asp_r[(int)founddir->dir] == 0 &&
- asp_c[(int)founddir->dir] == 0)
- G_fatal_error("no valid stream direction");
- r_nbr = r_nbr + asp_r[(int)founddir->dir];
- c_nbr = c_nbr + asp_c[(int)founddir->dir];
- draindir.pos = INDEX(r_nbr, c_nbr);
- if (stream[INDEX(r_nbr, c_nbr)] <= 0)
- G_fatal_error("stream id not set");
- else
- stream[INDEX(r_nbr, c_nbr)] = curr_stream;
- }
+ if (stream[INDEX(r_nbr, c_nbr)] <= 0)
+ G_fatal_error("stream id not set");
+ else
+ stream[INDEX(r_nbr, c_nbr)] = curr_stream;
}
- else {
- /* stream node already existing here */
- G_debug(2, "existing confluence");
- /* add new tributary to stream node */
- if (stream_node[curr_stream].n_trib >=
- stream_node[curr_stream].n_alloc) {
- size_t new_size;
+ }
+ else {
+ /* stream node already existing here */
+ G_debug(2, "existing confluence");
+ /* add new tributary to stream node */
+ if (stream_node[curr_stream].n_trib >=
+ stream_node[curr_stream].n_alloc) {
+ size_t new_size;
- stream_node[curr_stream].n_alloc += 2;
- new_size = stream_node[curr_stream].n_alloc * sizeof(int);
- stream_node[curr_stream].trib =
- (int *)G_realloc(stream_node[curr_stream].trib, new_size);
- new_size = stream_node[curr_stream].n_alloc * sizeof(double);
- stream_node[curr_stream].acc =
- (double *)G_realloc(stream_node[curr_stream].acc,
- new_size);
- }
-
- stream_node[curr_stream].trib[stream_node[curr_stream].n_trib++] =
- stream[thisindex];
+ stream_node[curr_stream].n_alloc += 2;
+ new_size = stream_node[curr_stream].n_alloc * sizeof(int);
+ stream_node[curr_stream].trib =
+ (int *)G_realloc(stream_node[curr_stream].trib, new_size);
+ new_size = stream_node[curr_stream].n_alloc * sizeof(double);
+ stream_node[curr_stream].acc =
+ (double *)G_realloc(stream_node[curr_stream].acc,
+ new_size);
}
- /* end new confluence */
- G_debug(2, "%d tribs", stream_node[curr_stream].n_trib);
- if (stream_node[curr_stream].n_trib == 1)
- G_warning("stream node %d only 1 trib: %d", curr_stream,
- stream_node[curr_stream].trib[0]);
+ stream_node[curr_stream].trib[stream_node[curr_stream].n_trib++] =
+ stream[thisindex];
}
+ G_debug(2, "%d tribs", stream_node[curr_stream].n_trib);
+ if (stream_node[curr_stream].n_trib == 1)
+ G_warning("stream node %d only 1 trib: %d", curr_stream,
+ stream_node[curr_stream].trib[0]);
+
return 1;
}
@@ -252,12 +260,11 @@
for (killer = 1; killer <= n_points; killer++) {
G_percent(killer, n_points, 1);
- r = astar_pts[killer].r;
- c = astar_pts[killer].c;
+ thisindex = astar_pts[killer].idx;
+ r = thisindex / ncols;
+ c = thisindex - r * ncols;
aspect = astar_pts[killer].asp;
- thisindex = INDEX(r, c);
-
/* do not distribute flow along edges */
if (aspect < 0) {
G_debug(3, "edge");
@@ -422,7 +429,7 @@
/*
* extracts streams for threshold, accumulation is provided
*/
-int extract_streams(double threshold, double mont_exp, int use_weight)
+int extract_streams(double threshold, double mont_exp, int use_weight, int min_length)
{
int r, c, dr, dc;
CELL is_swale, ele_val, ele_nbr;
@@ -495,12 +502,11 @@
for (killer = 1; killer <= n_points; killer++) {
G_percent(killer, n_points, 1);
- r = astar_pts[killer].r;
- c = astar_pts[killer].c;
+ thisindex = astar_pts[killer].idx;
+ r = thisindex / ncols;
+ c = thisindex - r * ncols;
aspect = astar_pts[killer].asp;
- thisindex = INDEX(r, c);
-
/* do not distribute flow along edges */
if (aspect < 0) {
G_debug(3, "edge");
@@ -536,11 +542,10 @@
value = acc[thisindex];
- /************************************/
+ /**********************************/
/* find main drainage direction */
+ /**********************************/
- /************************************/
-
max_acc = -1;
max_side = np_side = -1;
mfd_cells = 0;
@@ -609,7 +614,7 @@
is_swale = stream[thisindex];
- /* do not distribute flow along edges, this causes artifacts */
+ /* do not continue streams along edges, these are artifacts */
if (edge) {
G_debug(3, "edge");
if (is_swale) {
@@ -656,7 +661,6 @@
/**********************/
/* start new stream */
-
/**********************/
/* weight map has precedence over Montgomery */
@@ -670,7 +674,6 @@
G_warning
("Can't use Montgomery's method, no stream direction found");
else {
-
ele_nbr = ele[INDEX(r_max, c_max)];
slope = (double)(ele_val - ele_nbr) / ele_scale;
@@ -679,7 +682,6 @@
slope /= diag;
value *= pow(fabs(slope), mont_exp);
-
}
}
@@ -713,7 +715,6 @@
/*********************/
/* continue stream */
-
/*********************/
/* can't continue stream, add outlet point */
@@ -732,7 +733,7 @@
}
else if (is_swale > 0) {
continue_stream(is_swale, r, c, r_max, c_max, thisindex,
- &stream_no);
+ &stream_no, min_length);
}
FLAG_SET(worked, r, c);
@@ -742,7 +743,6 @@
workedon, n_points);
flag_destroy(worked);
-
G_free(dist_to_nbr);
G_debug(1, "%d outlets", n_outlets);
Modified: grass-addons/raster/r.stream.extract/thin.c
===================================================================
--- grass-addons/raster/r.stream.extract/thin.c 2009-10-13 17:15:48 UTC (rev 39512)
+++ grass-addons/raster/r.stream.extract/thin.c 2009-10-13 18:19:38 UTC (rev 39513)
@@ -100,6 +100,9 @@
thisindex = INDEX(r, c);
stream_id = stream[thisindex];
+ if (stream_id == 0)
+ continue;
+
/* add root node to stack */
G_debug(2, "add root node");
top = 0;
More information about the grass-commit
mailing list