[GRASS-SVN] r44578 - grass-addons/raster/r.stream.extract

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Dec 10 03:18:01 EST 2010


Author: mmetz
Date: 2010-12-10 00:18:01 -0800 (Fri, 10 Dec 2010)
New Revision: 44578

Added:
   grass-addons/raster/r.stream.extract/flag.c
Removed:
   grass-addons/raster/r.stream.extract/flag_clr_all.c
   grass-addons/raster/r.stream.extract/flag_create.c
   grass-addons/raster/r.stream.extract/flag_destroy.c
   grass-addons/raster/r.stream.extract/flag_get.c
   grass-addons/raster/r.stream.extract/flag_set.c
   grass-addons/raster/r.stream.extract/flag_unset.c
Modified:
   grass-addons/raster/r.stream.extract/description.html
   grass-addons/raster/r.stream.extract/do_astar.c
   grass-addons/raster/r.stream.extract/flag.h
   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
Log:
add option for real depressions, clear up flag code

Modified: grass-addons/raster/r.stream.extract/description.html
===================================================================
--- grass-addons/raster/r.stream.extract/description.html	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/description.html	2010-12-10 08:18:01 UTC (rev 44578)
@@ -8,7 +8,6 @@
 
 <dl>
 <dt><em>elevation</em> 
-
 <dd>Input map, required: Elevation on which entire analysis is based.
 NULL (nodata) cells are ignored, zero and negative values are valid
 elevation data. Gaps in the elevation map that are located within the
@@ -29,6 +28,16 @@
 as the number of cells, it can also be the optionally adjusted or
 weighed total contributing area in square meters or any other unit. 
 <p>
+<dt><em>depression</em> 
+<dd>Input map, optional: All non-NULL and non-zero cells will be regarded
+as real depressions. Streams will not be routed out of depressions. If an
+area is marked as depression but the elevation model has no depression
+at this location, streams will not stop there. If a flow accumulation map
+and a map with real depressions are provided, the flow accumulation map
+must match the depression map such that flow is not distributed out of 
+the indicated depressions. It is recommended to use internally computed
+flow accumulation if a depression map is provided.
+<p>
 <dt><em>threshold</em>
 <dd>Required: <em>threshold</em> for stream initiation by overland flow:
 the minumum (optionally modifed) flow accumulation value that will initiate

Modified: grass-addons/raster/r.stream.extract/do_astar.c
===================================================================
--- grass-addons/raster/r.stream.extract/do_astar.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/do_astar.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -128,13 +128,16 @@
 		    ele_up = ele[nindex];
 		    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];
 		    }
+		    /* neighbour is inside real depression, not yet worked */
+		    if (asp[nindex] == 0 && ele_val <= ele[nindex]) {
+			asp[nindex] = drain[r_nbr - r + 1][c_nbr - c + 1];
+		    }
 		}
 	    }
 	}    /* end neighbours */

Copied: grass-addons/raster/r.stream.extract/flag.c (from rev 40659, grass-addons/raster/r.stream.extract/flag_create.c)
===================================================================
--- grass-addons/raster/r.stream.extract/flag.c	                        (rev 0)
+++ grass-addons/raster/r.stream.extract/flag.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -0,0 +1,49 @@
+#include <grass/gis.h>
+#include "flag.h"
+
+FLAG *flag_create(int nrows, int ncols)
+{
+    unsigned char *temp;
+    FLAG *new_flag;
+    register int i;
+
+    new_flag = (FLAG *) G_malloc(sizeof(FLAG));
+    new_flag->nrows = nrows;
+    new_flag->ncols = ncols;
+    new_flag->leng = (ncols + 7) / 8;
+    new_flag->array =
+	(unsigned char **)G_malloc(nrows * sizeof(unsigned char *));
+
+    temp =
+	(unsigned char *)G_malloc(nrows * new_flag->leng *
+				  sizeof(unsigned char));
+
+    for (i = 0; i < nrows; i++) {
+	new_flag->array[i] = temp;
+	temp += new_flag->leng;
+    }
+
+    return (new_flag);
+}
+
+int flag_destroy(FLAG * flags)
+{
+    G_free(flags->array[0]);
+    G_free(flags->array);
+    G_free(flags);
+
+    return 0;
+}
+
+int flag_clear_all(FLAG * flags)
+{
+    register int r, c;
+
+    for (r = 0; r < flags->nrows; r++) {
+	for (c = 0; c < flags->leng; c++) {
+	    flags->array[r][c] = 0;
+	}
+    }
+
+    return 0;
+}

Modified: grass-addons/raster/r.stream.extract/flag.h
===================================================================
--- grass-addons/raster/r.stream.extract/flag.h	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/flag.h	2010-12-10 08:18:01 UTC (rev 44578)
@@ -20,18 +20,18 @@
  ** FLAG *flags;
  **     sets all values in flags to zero.
  **
- ** flag_unset(flags, row, col)
+ ** FLAG_UNSET(flags, row, col)
  ** FLAG *flags;
  ** int row, col;
  **     sets the value of (row, col) in flags to zero.
  **
- ** flag_set(flags, row, col)
+ ** FLAG_SET(flags, row, col)
  ** FLAG *flags;
  ** int row, col;
  **     will set the value of (row, col) in flags to one.
  **
  ** int
- ** flag_get(flags, row, col)
+ ** FLAG_GET(flags, row, col)
  ** FLAG *flags;
  ** int row, col;
  **     returns the value in flags that is at (row, col).
@@ -55,23 +55,9 @@
 #define FLAG_GET(flags,row,col) \
 	(flags)->array[(row)][(col)>>3] & (1<<((col) & 7))
 
-/* flag_clr_all.c */
+/* flag.c */
 int flag_clear_all(FLAG *);
-
-/* flag_create.c */
 FLAG *flag_create(int, int);
-
-/* flag_destroy.c */
 int flag_destroy(FLAG *);
 
-/* flag_get.c */
-int flag_get(FLAG *, int, int);
-
-/* flag_set.c */
-int flag_set(FLAG *, int, int);
-
-/* flag_unset.c */
-int flag_unset(FLAG *, int, int);
-
-
 #endif /* __FLAG_H__ */

Deleted: grass-addons/raster/r.stream.extract/flag_clr_all.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_clr_all.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/flag_clr_all.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -1,14 +0,0 @@
-#include "flag.h"
-
-int flag_clear_all(FLAG * flags)
-{
-    register int r, c;
-
-    for (r = 0; r < flags->nrows; r++) {
-	for (c = 0; c < flags->leng; c++) {
-	    flags->array[r][c] = 0;
-	}
-    }
-
-    return 0;
-}

Deleted: grass-addons/raster/r.stream.extract/flag_create.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_create.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/flag_create.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -1,27 +0,0 @@
-#include <grass/gis.h>
-#include "flag.h"
-
-FLAG *flag_create(int nrows, int ncols)
-{
-    unsigned char *temp;
-    FLAG *new_flag;
-    register int i;
-
-    new_flag = (FLAG *) G_malloc(sizeof(FLAG));
-    new_flag->nrows = nrows;
-    new_flag->ncols = ncols;
-    new_flag->leng = (ncols + 7) / 8;
-    new_flag->array =
-	(unsigned char **)G_malloc(nrows * sizeof(unsigned char *));
-
-    temp =
-	(unsigned char *)G_malloc(nrows * new_flag->leng *
-				  sizeof(unsigned char));
-
-    for (i = 0; i < nrows; i++) {
-	new_flag->array[i] = temp;
-	temp += new_flag->leng;
-    }
-
-    return (new_flag);
-}

Deleted: grass-addons/raster/r.stream.extract/flag_destroy.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_destroy.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/flag_destroy.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -1,11 +0,0 @@
-#include <grass/gis.h>
-#include "flag.h"
-
-int flag_destroy(FLAG * flags)
-{
-    G_free(flags->array[0]);
-    G_free(flags->array);
-    G_free(flags);
-
-    return 0;
-}

Deleted: grass-addons/raster/r.stream.extract/flag_get.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_get.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/flag_get.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -1,6 +0,0 @@
-#include "flag.h"
-
-int flag_get(FLAG * flags, int row, int col)
-{
-    return (flags->array[row][col >> 3] & (1 << (col & 7)));
-}

Deleted: grass-addons/raster/r.stream.extract/flag_set.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_set.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/flag_set.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -1,8 +0,0 @@
-#include "flag.h"
-
-int flag_set(FLAG * flags, int row, int col)
-{
-    flags->array[row][col >> 3] |= (1 << (col & 7));
-
-    return 0;
-}

Deleted: grass-addons/raster/r.stream.extract/flag_unset.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_unset.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/flag_unset.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -1,8 +0,0 @@
-#include "flag.h"
-
-int flag_unset(FLAG * flags, int row, int col)
-{
-    flags->array[row][col >> 3] &= ~(1 << (col & 7));
-
-    return 0;
-}

Modified: grass-addons/raster/r.stream.extract/load.c
===================================================================
--- grass-addons/raster/r.stream.extract/load.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/load.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -23,12 +23,13 @@
  * gets start points for A* Search
  * start points are edges
  */
-int load_maps(int ele_fd, int acc_fd)
+int load_maps(int ele_fd, int acc_fd, int depr_fd)
 {
     int r, c, thisindex;
     char asp_value, *aspp;
     void *ele_buf, *ptr, *acc_buf = NULL, *acc_ptr = NULL;
     CELL *loadp, ele_value;
+    CELL *depr_buf;
     DCELL dvalue;
     int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
     int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
@@ -124,14 +125,19 @@
 		if (acc_fd < 0)
 		    *accp = 1;
 		else {
-		    if (acc_map_type == CELL_TYPE) {
+		    if (G_is_null_value(acc_ptr, acc_map_type))
+			G_fatal_error(_("Accumulation map does not match elevation map!"));
+
+		    switch (acc_map_type) {
+		    case CELL_TYPE:
 			*accp = *((CELL *) acc_ptr);
-		    }
-		    else if (acc_map_type == FCELL_TYPE) {
+			break;
+		    case FCELL_TYPE:
 			*accp = *((FCELL *) acc_ptr);
-		    }
-		    else if (acc_map_type == DCELL_TYPE) {
+			break;
+		    case DCELL_TYPE:
 			*accp = *((DCELL *) acc_ptr);
+			break;
 		    }
 		}
 
@@ -168,12 +174,24 @@
 
     nxt_avail_pt = heap_size = 0;
 
-    /* load edge cells to A* heap */
+    /* load edge cells and real depressions to A* heap */
+    if (depr_fd >= 0)
+	depr_buf = G_allocate_raster_buf(CELL_TYPE);
+    else
+	depr_buf = NULL;
+
     G_message(_("set edge points"));
     loadp = ele;
     for (r = 0; r < nrows; r++) {
+	G_percent(r, nrows, 2);
 
-	G_percent(r, nrows, 2);
+	if (depr_fd >= 0) {
+	    if (G_get_raster_row(depr_fd, depr_buf, r, CELL_TYPE) < 0) {
+		G_warning(_("could not read raster map at row <%d>"), r);
+		return -1;
+	    }
+	}
+
 	for (c = 0; c < ncols; c++) {
 
 	    is_worked = FLAG_GET(worked, r, c);
@@ -203,13 +221,13 @@
 
 		thisindex = INDEX(r, c);
 		ele_value = ele[thisindex];
+		asp[thisindex] = asp_value;
 		heap_add(r, c, ele_value, asp_value);
-		asp[thisindex] = asp_value;
-		FLAG_SET(in_list, r, c);
 		continue;
 	    }
 
 	    /* any neighbour NULL ? */
+	    asp_value = 0;
 	    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 		/* get r, c (r_nbr, c_nbr) for neighbours */
 		r_nbr = r + nextdr[ct_dir];
@@ -221,17 +239,33 @@
 		    asp_value = drain[r - r_nbr + 1][c - c_nbr + 1];
 		    thisindex = INDEX(r, c);
 		    ele_value = ele[thisindex];
+		    asp[thisindex] = asp_value;
 		    heap_add(r, c, ele_value, asp_value);
-		    asp[thisindex] = asp_value;
-		    FLAG_SET(in_list, r, c);
 
 		    break;
 		}
 	    }
+	    if (asp_value) /* some neighbour was NULL, point added to list */
+		continue;
+	    
+	    /* real depression ? */
+	    if (depr_fd >= 0) {
+		if (!G_is_c_null_value(&depr_buf[c]) && depr_buf[c] != 0) {
+		    thisindex = INDEX(r, c);
+		    ele_value = ele[thisindex];
+		    asp[thisindex] = 0;
+		    heap_add(r, c, ele_value, 0);
+		}
+	    }
 	}
     }
     G_percent(nrows, nrows, 2);	/* finish it */
 
+    if (depr_fd >= 0) {
+	G_close_cell(depr_fd);
+	G_free(depr_buf);
+    }
+
     G_debug(1, "%d edge cells", heap_size);
     G_debug(1, "%d non-NULL cells", n_points);
 

Modified: grass-addons/raster/r.stream.extract/local_proto.h
===================================================================
--- grass-addons/raster/r.stream.extract/local_proto.h	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/local_proto.h	2010-12-10 08:18:01 UTC (rev 44578)
@@ -30,27 +30,29 @@
     double *acc;
 } *stream_node;
 
-int nrows, ncols;
-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;
-struct point *outlets;
-unsigned int n_outlets, n_alloc_outlets;
-DCELL *acc;
-CELL *ele;
-char *asp;
-CELL *stream;
-FLAG *worked, *in_list;
+/* global variables */
+extern int nrows, ncols;
+extern unsigned int *astar_pts;
+extern unsigned int n_search_points, n_points, nxt_avail_pt;
+extern unsigned int heap_size, *astar_added;
+extern unsigned int n_stream_nodes, n_alloc_nodes;
+extern struct point *outlets;
+extern unsigned int n_outlets, n_alloc_outlets;
+extern DCELL *acc;
+extern CELL *ele;
+extern char *asp;
+extern CELL *stream;
+extern FLAG *worked, *in_list;
 extern char drain[3][3];
-unsigned int first_cum;
-char sides;
-int c_fac;
-int ele_scale;
-struct RB_TREE *draintree;
+extern unsigned int first_cum;
+extern char sides;
+extern int c_fac;
+extern int ele_scale;
+extern int have_depressions;
+extern struct RB_TREE *draintree;
 
 /* load.c */
-int load_maps(int, int);
+int load_maps(int, int, int);
 
 /* do_astar.c */
 int do_astar(void);

Modified: grass-addons/raster/r.stream.extract/main.c
===================================================================
--- grass-addons/raster/r.stream.extract/main.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/main.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -20,14 +20,33 @@
 #include <grass/glocale.h>
 #include "local_proto.h"
 
+/* global variables */
+int nrows, ncols;
+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;
+struct point *outlets;
+unsigned int n_outlets, n_alloc_outlets;
+DCELL *acc;
+CELL *ele;
+char *asp;
+CELL *stream;
+FLAG *worked, *in_list;
 char drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
+unsigned int first_cum;
+char sides;
+int c_fac;
+int ele_scale;
+int have_depressions;
+struct RB_TREE *draintree;
 
 
 int main(int argc, char *argv[])
 {
     struct
     {
-	struct Option *ele, *acc;
+	struct Option *ele, *acc, *depression;
 	struct Option *threshold, *d8cut;
 	struct Option *mont_exp;
 	struct Option *min_stream_length;
@@ -39,7 +58,7 @@
 	struct Option *dir_rast;
     } output;
     struct GModule *module;
-    int ele_fd, acc_fd;
+    int ele_fd, acc_fd, depr_fd;
     double threshold, d8cut, mont_exp;
     int min_stream_length = 0;
     char *mapset;
@@ -63,6 +82,13 @@
     input.acc->description =
 	_("Stream extraction will use provided accumulation instead of calculating it anew");
 
+    input.depression = G_define_standard_option(G_OPT_R_INPUT);
+    input.depression->key = "depression";
+    input.depression->label = _("Map with real depressions");
+    input.depression->required = NO;
+    input.depression->description =
+	_("Streams will not be routed out of real depressions");
+
     input.threshold = G_define_option();
     input.threshold->key = "threshold";
     input.threshold->label = _("Minimum flow accumulation for streams");
@@ -137,6 +163,14 @@
 	    G_fatal_error(_("Raster map <%s> not found"), input.acc->answer);
     }
 
+    if (input.depression->answer) {
+	if (!G_find_cell(input.depression->answer, ""))
+	    G_fatal_error(_("Raster map <%s> not found"), input.depression->answer);
+	have_depressions = 1;
+    }
+    else
+	have_depressions = 0;
+
     /* threshold makes sense */
     threshold = atof(input.threshold->answer);
     if (threshold <= 0)
@@ -203,6 +237,16 @@
     else
 	acc_fd = -1;
 
+    if (input.depression->answer) {
+	mapset = G_find_cell2(input.depression->answer, "");
+	depr_fd = G_open_cell_old(input.depression->answer, mapset);
+	if (depr_fd < 0)
+	    G_fatal_error(_("could not open input map %s"),
+			  input.depression->answer);
+    }
+    else
+	depr_fd = -1;
+
     /* set global variables */
     nrows = G_window_rows();
     ncols = G_window_cols();
@@ -215,7 +259,7 @@
     acc = (DCELL *) G_malloc(nrows * ncols * sizeof(DCELL));
 
     /* load maps */
-    if (load_maps(ele_fd, acc_fd) < 0)
+    if (load_maps(ele_fd, acc_fd, depr_fd) < 0)
 	G_fatal_error(_("could not load input map(s)"));
 
     /********************/

Modified: grass-addons/raster/r.stream.extract/streams.c
===================================================================
--- grass-addons/raster/r.stream.extract/streams.c	2010-12-10 03:59:25 UTC (rev 44577)
+++ grass-addons/raster/r.stream.extract/streams.c	2010-12-10 08:18:01 UTC (rev 44578)
@@ -217,7 +217,7 @@
     DCELL value, valued;
     int count;
     struct Cell_head window;
-    int mfd_cells, astar_not_set, is_null;
+    int mfd_cells, astar_not_set;
     double *dist_to_nbr, *weight, sum_weight, max_weight;
     double dx, dy;
     int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
@@ -257,6 +257,7 @@
     /* reset worked flag */
     flag_clear_all(worked);
 
+    /* distribute and accumulate */
     for (killer = 1; killer <= n_points; killer++) {
 	G_percent(killer, n_points, 1);
 
@@ -272,6 +273,12 @@
 	else {
 	    dr = r;
 	    dc = c;
+	    /* can only happen with real depressions */
+	    if (!have_depressions)
+		G_fatal_error(_("Bug in stream extraction"));
+	    FLAG_SET(worked, r, c);
+	    G_debug(1, "bottom of real depression");
+	    continue;
 	}
 
 	r_max = dr;
@@ -289,7 +296,6 @@
 	mfd_cells = 0;
 	astar_not_set = 1;
 	ele_val = ele[thisindex];
-	is_null = 0;
 	edge = 0;
 	/* this loop is needed to get the sum of weights */
 	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -305,9 +311,8 @@
 		is_worked = FLAG_GET(worked, r_nbr, c_nbr);
 		if (is_worked == 0) {
 		    ele_nbr = ele[nindex];
-		    is_null = G_is_c_null_value(&ele_nbr);
-		    edge = is_null;
-		    if (!is_null && ele_nbr <= ele_val) {
+		    edge = G_is_c_null_value(&ele_nbr);
+		    if (!edge && ele_nbr <= ele_val) {
 			if (ele_nbr < ele_val) {
 			    weight[ct_dir] =
 				mfd_pow((ele_val -
@@ -341,6 +346,7 @@
 	/* do not distribute flow along edges, this causes artifacts */
 	if (edge) {
 	    G_debug(3, "edge");
+	    FLAG_SET(worked, r, c);
 	    continue;
 	}
 
@@ -501,6 +507,7 @@
 
     workedon = 0;
 
+    /* extract streams */
     for (killer = 1; killer <= n_points; killer++) {
 	G_percent(killer, n_points, 1);
 
@@ -510,7 +517,7 @@
 	aspect = asp[thisindex];
 
 	/* do not distribute flow along edges */
-	if (aspect < 0) {
+	if (aspect <= 0) {
 	    G_debug(3, "edge");
 	    is_swale = stream[thisindex];
 	    if (is_swale) {
@@ -527,6 +534,13 @@
 		outlets[n_outlets].c = c;
 		n_outlets++;
 	    }
+	    FLAG_SET(worked, r, c);
+	    if (aspect == 0) {
+		/* can only happen with real depressions */
+		if (!have_depressions)
+		    G_fatal_error(_("Bug in stream extraction"));
+		G_debug(1, "bottom of real depression");
+	    } 
 	    continue;
 	}
 
@@ -535,6 +549,8 @@
 	    dc = c + asp_c[abs((int)aspect)];
 	}
 	else {
+	    /* can only happen with real depressions,
+	     * but should not get to here */
 	    dr = r;
 	    dc = c;
 	}
@@ -579,7 +595,7 @@
 		/* check for stream cells */
 		valued = fabs(acc[nindex]);
 		ele_nbr = ele[nindex];
-		/* if (valued >= threshold) */
+		/* check all upstream neighbours */
 		if (valued >= threshold && ct_dir != np_side &&
 		    ele_nbr > ele_val)
 		    stream_cells++;
@@ -659,6 +675,7 @@
 	}
 
 	/* update aspect */
+	/* r_max == r && c_max == c should not happen */
 	if ((r_max != dr || c_max != dc) && (r_max != r || c_max != c)) {
 	    asp[thisindex] = drain[r - r_max + 1][c - c_max + 1];
 	}
@@ -674,7 +691,7 @@
 	if (mont_exp > 0) {
 	    if (r_max == r && c_max == c)
 		G_warning
-		    ("Can't use Montgomery's method, no stream direction found");
+		    (_("Can't use Montgomery's method, no stream direction found"));
 	    else {
 		ele_nbr = ele[INDEX(r_max, c_max)];
 
@@ -720,9 +737,10 @@
 	/*********************/
 
 	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);
+		/* can't continue stream, add outlet point
+		 * r_max == r && c_max == c should not happen */
+		G_debug(1, "can't continue stream at r %d c %d", r, c);
 
 		if (n_outlets >= n_alloc_outlets) {
 		    n_alloc_outlets += stream_node_step;



More information about the grass-commit mailing list