[GRASS-SVN] r40122 - grass-addons/raster/r.stream.extract
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Dec 23 09:35:07 EST 2009
Author: mmetz
Date: 2009-12-23 09:35:06 -0500 (Wed, 23 Dec 2009)
New Revision: 40122
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/load.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:
output full flow direction
Modified: grass-addons/raster/r.stream.extract/close.c
===================================================================
--- grass-addons/raster/r.stream.extract/close.c 2009-12-23 12:56:05 UTC (rev 40121)
+++ grass-addons/raster/r.stream.extract/close.c 2009-12-23 14:35:06 UTC (rev 40122)
@@ -236,7 +236,6 @@
CELL *cell_buf1, *cell_buf2;
unsigned int thisindex;
struct History history;
- struct ddir draindir, *founddir;
/* cheating... */
stream_fd = dir_fd = -1;
@@ -266,17 +265,13 @@
if (stream[thisindex] > 0) {
if (stream_rast)
cell_buf1[c] = stream[thisindex];
- if (dir_rast) {
- draindir.pos = thisindex;
- if ((founddir =
- rbtree_find(draintree, &draindir)) != NULL) {
- cell_buf2[c] = founddir->dir;
- }
- else {
- cell_buf2[c] = 0;
- }
+ }
+ if (dir_rast) {
+ if (!G_is_c_null_value(&ele[thisindex])) {
+ cell_buf2[c] = asp[thisindex];
}
}
+
}
if (stream_rast)
G_put_raster_row(stream_fd, cell_buf1, CELL_TYPE);
Modified: grass-addons/raster/r.stream.extract/description.html
===================================================================
--- grass-addons/raster/r.stream.extract/description.html 2009-12-23 12:56:05 UTC (rev 40121)
+++ grass-addons/raster/r.stream.extract/description.html 2009-12-23 14:35:06 UTC (rev 40122)
@@ -90,13 +90,11 @@
location of confluences.
<p>
<dt><em>direction</em>
-<dd>Output raster map with flow direction for extracted streams. Flow
-direction is of D8 type with a range of 1 to 8. Multiplying values with
-45 gives degrees CCW from East. Flow direction was adjusted during
-thinning, taking shortcuts and skipping cells that were eliminated by
-the thinning procedure. Non-stream cells are set to NULL. A full,
-corrected flow direction map can be created by patching the
-<em>direction</em> output map with the flow direction map of r.watershed.
+<dd>Output raster map with flow direction for all non-NULL cells in
+input elevation. Flow direction is of D8 type with a range of 1 to 8.
+Multiplying values with 45 gives degrees CCW from East. Flow direction
+was adjusted during thinning, taking shortcuts and skipping cells that
+were eliminated by the thinning procedure.
</dl>
<h2>NOTES</h2>
Modified: grass-addons/raster/r.stream.extract/do_astar.c
===================================================================
--- grass-addons/raster/r.stream.extract/do_astar.c 2009-12-23 12:56:05 UTC (rev 40121)
+++ grass-addons/raster/r.stream.extract/do_astar.c 2009-12-23 14:35:06 UTC (rev 40122)
@@ -21,12 +21,11 @@
int do_astar(void)
{
int r, c, r_nbr, c_nbr, ct_dir;
- struct ast_point astp;
+ unsigned int astp;
int count, is_in_list, is_worked;
int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
CELL ele_val, ele_up, ele_nbr[8];
- char asp_val;
unsigned int thisindex, nindex;
/* sides
* |7|1|4|
@@ -38,6 +37,7 @@
double dx, dy, dist_to_nbr[8], ew_res, ns_res;
double slope[8];
struct Cell_head window;
+ int skip_me;
count = 0;
@@ -80,7 +80,7 @@
astar_pts[first_cum] = astp;
first_cum--;
- thisindex = astp.idx;
+ thisindex = astp;
r = thisindex / ncols;
c = thisindex - r * ncols;
@@ -91,43 +91,51 @@
r_nbr = r + nextdr[ct_dir];
c_nbr = c + nextdc[ct_dir];
slope[ct_dir] = ele_nbr[ct_dir] = 0;
+ skip_me = 0;
/* check that neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 && c_nbr < ncols) {
is_in_list = FLAG_GET(in_list, r_nbr, c_nbr);
is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ nindex = INDEX(r_nbr, c_nbr);
if (!is_worked) {
- nindex = INDEX(r_nbr, c_nbr);
ele_nbr[ct_dir] = ele[nindex];
slope[ct_dir] =
get_slope2(ele_val, ele_nbr[ct_dir],
dist_to_nbr[ct_dir]);
+ }
+ if (!is_in_list) {
/* avoid diagonal flow direction bias */
- if (ct_dir > 3) {
- if (slope[nbr_ew[ct_dir]]) {
+ if (ct_dir > 3 && slope[ct_dir] > 0) {
+ if (slope[nbr_ew[ct_dir]] > 0) {
/* slope to ew nbr > slope to center */
if (slope[ct_dir] <
get_slope2(ele_nbr[nbr_ew[ct_dir]],
ele_nbr[ct_dir], ew_res))
- is_in_list = 1;
+ skip_me = 1;
}
- if (!is_in_list && slope[nbr_ns[ct_dir]]) {
+ if (!skip_me && slope[nbr_ns[ct_dir]] > 0) {
/* slope to ns nbr > slope to center */
if (slope[ct_dir] <
get_slope2(ele_nbr[nbr_ns[ct_dir]],
ele_nbr[ct_dir], ns_res))
- is_in_list = 1;
+ skip_me = 1;
}
}
}
- if (is_in_list == 0) {
- nindex = INDEX(r_nbr, c_nbr);
+ if (is_in_list == 0 && skip_me == 0) {
ele_up = ele[nindex];
- asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
- heap_add(r_nbr, c_nbr, ele_up, asp_val);
+ asp[nindex] = drain[r_nbr - r + 1][c_nbr - c + 1];
+ heap_add(r_nbr, c_nbr, ele_up, asp[nindex]);
FLAG_SET(in_list, r_nbr, c_nbr);
}
+ else if (is_in_list && is_worked == 0) {
+ /* neighbour is edge in list, not yet worked */
+ if (asp[nindex] < 0) {
+ asp[nindex] = drain[r_nbr - r + 1][c_nbr - c + 1];
+ }
+ }
}
} /* end neighbours */
FLAG_SET(worked, r, c);
@@ -161,7 +169,7 @@
{
unsigned int child, child_added, parent;
CELL elep;
- struct ast_point childp;
+ unsigned int childp;
child = start;
child_added = astar_added[child];
@@ -170,7 +178,7 @@
while (child > 1) {
parent = get_parent(child);
- elep = ele[astar_pts[parent].idx];
+ elep = ele[astar_pts[parent]];
/* child < parent */
if (heap_cmp(elec, child_added, elep, astar_added[parent]) == 1) {
@@ -210,8 +218,7 @@
G_fatal_error(_("heapsize too large"));
astar_added[heap_size] = nxt_avail_pt;
- astar_pts[heap_size].idx = INDEX(r, c);
- astar_pts[heap_size].asp = asp;
+ astar_pts[heap_size] = INDEX(r, c);
nxt_avail_pt++;
@@ -240,13 +247,13 @@
parent = 1;
while ((child = get_child(parent)) <= heap_size) {
- elec = ele[astar_pts[child].idx];
+ elec = ele[astar_pts[child]];
if (child < heap_size) {
childr = child + 1;
i = child + 3;
while (childr <= heap_size && childr < i) {
- eler = ele[astar_pts[childr].idx];
+ eler = ele[astar_pts[childr]];
if (heap_cmp
(eler, astar_added[childr], elec,
@@ -269,7 +276,7 @@
astar_added[parent] = astar_added[heap_size];
astar_pts[parent] = astar_pts[heap_size];
- elec = ele[astar_pts[parent].idx];
+ elec = ele[astar_pts[parent]];
/* sift up last swapped point, only necessary if hole moved to heap end */
sift_up(parent, elec);
}
Modified: grass-addons/raster/r.stream.extract/load.c
===================================================================
--- grass-addons/raster/r.stream.extract/load.c 2009-12-23 12:56:05 UTC (rev 40121)
+++ grass-addons/raster/r.stream.extract/load.c 2009-12-23 14:35:06 UTC (rev 40122)
@@ -26,7 +26,7 @@
int load_maps(int ele_fd, int acc_fd, int weight_fd)
{
int r, c, thisindex;
- char asp_value;
+ char asp_value, *aspp;
void *ele_buf, *ptr, *acc_buf = NULL, *acc_ptr = NULL, *weight_buf =
NULL, *weight_ptr = NULL;
CELL *loadp, ele_value;
@@ -37,7 +37,6 @@
int is_worked;
size_t ele_size, acc_size = 0, weight_size = 0;
int ele_map_type, acc_map_type = 0, weight_map_type = 0;
- CELL *streamp;
DCELL *accp, *weightp;
if (acc_fd < 0 && weight_fd < 0)
@@ -85,7 +84,6 @@
in_list = flag_create(nrows, ncols);
loadp = ele;
- streamp = stream;
accp = acc;
weightp = accweight;
@@ -122,8 +120,6 @@
FLAG_UNSET(worked, r, c);
FLAG_UNSET(in_list, r, c);
- *streamp = 0;
-
/* check for masked and NULL cells */
if (G_is_null_value(ptr, ele_map_type)) {
FLAG_SET(worked, r, c);
@@ -177,7 +173,6 @@
loadp++;
accp++;
- streamp++;
ptr = G_incr_void_ptr(ptr, ele_size);
if (acc_fd >= 0)
acc_ptr = G_incr_void_ptr(acc_ptr, acc_size);
@@ -203,8 +198,8 @@
}
astar_pts =
- (struct ast_point *)G_malloc((n_points + 1) *
- sizeof(struct ast_point));
+ (unsigned int *)G_malloc((n_points + 1) *
+ sizeof(unsigned int));
/* astar_heap will track astar_pts in ternary min-heap */
/* astar_heap is one-based */
@@ -216,11 +211,13 @@
/* load edge cells to A* heap */
G_message(_("set edge points"));
loadp = ele;
+ aspp = asp;
for (r = 0; r < nrows; r++) {
G_percent(r, nrows, 2);
for (c = 0; c < ncols; c++) {
+ *aspp = 0;
is_worked = FLAG_GET(worked, r, c);
if (is_worked)
@@ -249,6 +246,7 @@
thisindex = INDEX(r, c);
ele_value = ele[thisindex];
heap_add(r, c, ele_value, asp_value);
+ asp[thisindex] = asp_value;
FLAG_SET(in_list, r, c);
continue;
}
@@ -266,12 +264,13 @@
thisindex = INDEX(r, c);
ele_value = ele[thisindex];
heap_add(r, c, ele_value, asp_value);
+ asp[thisindex] = asp_value;
FLAG_SET(in_list, r, c);
break;
}
}
-
+ aspp++;
}
}
G_percent(nrows, nrows, 2); /* finish it */
Modified: grass-addons/raster/r.stream.extract/local_proto.h
===================================================================
--- grass-addons/raster/r.stream.extract/local_proto.h 2009-12-23 12:56:05 UTC (rev 40121)
+++ grass-addons/raster/r.stream.extract/local_proto.h 2009-12-23 14:35:06 UTC (rev 40122)
@@ -14,12 +14,6 @@
int dir;
};
-struct ast_point
-{
- int idx;
- char asp;
-};
-
struct point
{
int r, c;
@@ -37,7 +31,7 @@
} *stream_node;
int nrows, ncols;
-struct ast_point *astar_pts;
+unsigned int *astar_pts;
unsigned int n_search_points, n_points, nxt_avail_pt;
unsigned int heap_size, *astar_added;
unsigned int n_stream_nodes, n_alloc_nodes;
@@ -45,8 +39,8 @@
unsigned int n_outlets, n_alloc_outlets;
DCELL *acc, *accweight;
CELL *ele;
+char *asp;
CELL *stream;
-int *strahler, *horton; /* strahler and horton order */
FLAG *worked, *in_list;
extern char drain[3][3];
unsigned int first_cum;
Modified: grass-addons/raster/r.stream.extract/main.c
===================================================================
--- grass-addons/raster/r.stream.extract/main.c 2009-12-23 12:56:05 UTC (rev 40121)
+++ grass-addons/raster/r.stream.extract/main.c 2009-12-23 14:35:06 UTC (rev 40122)
@@ -124,7 +124,7 @@
output.dir_rast = G_define_standard_option(G_OPT_R_OUTPUT);
output.dir_rast->key = "direction";
output.dir_rast->description =
- _("Output raster map with flow direction for streams");
+ _("Output raster map with flow direction");
output.dir_rast->required = NO;
output.dir_rast->guisection = _("Output options");
@@ -233,9 +233,8 @@
/* allocate memory */
ele = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
- /* TODO: allocate acc and stream memory only after A* Search */
+ asp = (char *) G_malloc(nrows * ncols * sizeof(char));
acc = (DCELL *) G_malloc(nrows * ncols * sizeof(DCELL));
- stream = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
if (input.weight->answer)
accweight = (DCELL *) G_malloc(nrows * ncols * sizeof(DCELL));
else
@@ -258,15 +257,12 @@
if (do_accum(d8cut) < 0)
G_fatal_error(_("could not calculate flow accumulation"));
}
- else {
- /* load accumulation */
- }
+ stream = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
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)
G_free(accweight);
@@ -286,7 +282,9 @@
output.dir_rast->answer) < 0)
G_fatal_error(_("could not write output maps"));
+ G_free(ele);
G_free(stream);
+ G_free(asp);
exit(EXIT_SUCCESS);
}
Modified: grass-addons/raster/r.stream.extract/streams.c
===================================================================
--- grass-addons/raster/r.stream.extract/streams.c 2009-12-23 12:56:05 UTC (rev 40121)
+++ grass-addons/raster/r.stream.extract/streams.c 2009-12-23 14:35:06 UTC (rev 40122)
@@ -260,17 +260,11 @@
for (killer = 1; killer <= n_points; killer++) {
G_percent(killer, n_points, 1);
- thisindex = astar_pts[killer].idx;
+ thisindex = astar_pts[killer];
r = thisindex / ncols;
c = thisindex - r * ncols;
- aspect = astar_pts[killer].asp;
+ aspect = asp[thisindex];
- /* do not distribute flow along edges */
- if (aspect < 0) {
- G_debug(3, "edge");
- continue;
- }
-
if (aspect) {
dr = r + asp_r[abs((int)aspect)];
dc = c + asp_c[abs((int)aspect)];
@@ -287,7 +281,6 @@
/***************************************/
/* get weights for flow distribution */
-
/***************************************/
max_weight = 0;
@@ -359,7 +352,6 @@
/************************************/
/* distribute and accumulate flow */
-
/************************************/
/* MFD, A * path not included, add to mfd_cells */
@@ -433,6 +425,7 @@
{
int r, c, dr, dc;
CELL is_swale, ele_val, ele_nbr;
+ CELL *streamp;
DCELL value, valued;
struct Cell_head window;
int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
@@ -497,15 +490,24 @@
/* reset worked flag */
flag_clear_all(worked);
+ /* initialize streams */
+ streamp = stream;
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ *streamp = 0;
+ streamp++;
+ }
+ }
+
workedon = 0;
for (killer = 1; killer <= n_points; killer++) {
G_percent(killer, n_points, 1);
- thisindex = astar_pts[killer].idx;
+ thisindex = astar_pts[killer];
r = thisindex / ncols;
c = thisindex - r * ncols;
- aspect = astar_pts[killer].asp;
+ aspect = asp[thisindex];
/* do not distribute flow along edges */
if (aspect < 0) {
@@ -537,8 +539,8 @@
dc = c;
}
- r_max = dr;
- c_max = dc;
+ r_nbr = r_max = dr;
+ c_nbr = c_max = dc;
value = acc[thisindex];
@@ -556,7 +558,7 @@
is_null = 0;
edge = 0;
flat = 1;
- /* this loop is needed to get the sum of weights */
+ /* find main drainage direction */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r_nbr, c_nbr for neighbours */
r_nbr = r + nextdr[ct_dir];
@@ -630,6 +632,10 @@
outlets[n_outlets].r = r;
outlets[n_outlets].c = c;
n_outlets++;
+ if (asp[thisindex] > 0) {
+ aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+ asp[thisindex] = aspect;
+ }
}
continue;
}
@@ -657,6 +663,11 @@
max_side = np_side;
}
+ /* update aspect */
+ if ((r_max != dr || c_max != dc) && (r_max != r || c_max != c)) {
+ asp[thisindex] = drain[r - r_max + 1][c - c_max + 1];
+ }
+
is_swale = stream[thisindex];
/**********************/
@@ -717,24 +728,26 @@
/* continue stream */
/*********************/
- /* can't continue stream, add outlet point */
- if (is_swale > 0 && r_max == r && c_max == c) {
- G_debug(2, "can't continue stream at %d %d", r, c);
+ if (is_swale > 0) {
+ /* can't continue stream, add outlet point */
+ if (r_max == r && c_max == c) {
+ G_debug(2, "can't continue stream at %d %d", r, c);
- if (n_outlets >= n_alloc_outlets) {
- n_alloc_outlets += stream_node_step;
- outlets =
- (struct point *)G_malloc(n_alloc_outlets *
- sizeof(struct point));
+ if (n_outlets >= n_alloc_outlets) {
+ n_alloc_outlets += stream_node_step;
+ outlets =
+ (struct point *)G_malloc(n_alloc_outlets *
+ sizeof(struct point));
+ }
+ outlets[n_outlets].r = r;
+ outlets[n_outlets].c = c;
+ n_outlets++;
}
- outlets[n_outlets].r = r;
- outlets[n_outlets].c = c;
- n_outlets++;
+ else {
+ continue_stream(is_swale, r, c, r_max, c_max, thisindex,
+ &stream_no, min_length);
+ }
}
- else if (is_swale > 0) {
- continue_stream(is_swale, r, c, r_max, c_max, thisindex,
- &stream_no, min_length);
- }
FLAG_SET(worked, r, c);
}
@@ -744,6 +757,7 @@
flag_destroy(worked);
G_free(dist_to_nbr);
+ G_free(astar_pts);
G_debug(1, "%d outlets", n_outlets);
G_debug(1, "%d nodes", n_stream_nodes);
Modified: grass-addons/raster/r.stream.extract/thin.c
===================================================================
--- grass-addons/raster/r.stream.extract/thin.c 2009-12-23 12:56:05 UTC (rev 40121)
+++ grass-addons/raster/r.stream.extract/thin.c 2009-12-23 14:35:06 UTC (rev 40122)
@@ -58,6 +58,7 @@
last_r = r_nbr;
last_c = c_nbr;
draindir.pos = INDEX(last_r, last_c);
+ asp[draindir.pos] = founddir->dir;
thinned = 1;
}
More information about the grass-commit
mailing list