[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