more debug info for diagonal flow bias, improved stream extraction, comment code

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@40634 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.watershed/ram/Gwater.h

@@ -18,7 +18,6 @@
 
 #define AR_SIZE			16
 #define AR_INCR			16
-#define SHORT			int
 #define NOMASK			1
 #define MIN_SLOPE		.00001
 #define MIN_GRADIENT_DEGREES	1
@@ -35,7 +34,7 @@
 
 #define POINT       struct points
 POINT {
-    SHORT r, c; /* , downr, downc */
+    int r, c; /* , downr, downc */
     /* int nxt; */
 };
 
@@ -49,7 +48,7 @@ extern struct Cell_head window;
 extern int mfd, c_fac, abs_acc, ele_scale;
 extern int *heap_index, heap_size;
 extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
-extern SHORT nrows, ncols;
+extern int nrows, ncols;
 extern double half_res, diag, max_length, dep_slope;
 extern int bas_thres, tot_parts;
 extern CELL n_basins;
@@ -66,11 +65,11 @@ extern int ril_fd;
 extern double *s_l, *s_g, *l_s;
 extern CELL one, zero;
 extern double ril_value, d_one, d_zero;
-extern SHORT sides;
-extern SHORT drain[3][3];
-extern SHORT updrain[3][3];
-extern SHORT nextdr[8];
-extern SHORT nextdc[8];
+extern int sides;
+extern int drain[3][3];
+extern int updrain[3][3];
+extern int nextdr[8];
+extern int nextdc[8];
 extern char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
 extern char run_name[GNAME_MAX], ob_name[GNAME_MAX];
 extern char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
@@ -95,11 +94,11 @@ CELL def_basin(int, int, CELL, double, CELL);
 
 /* do_astar.c */
 int do_astar(void);
-int add_pt(SHORT, SHORT, CELL, CELL);
+int add_pt(int, int, CELL, CELL);
 int drop_pt(void);
 int sift_up(int, CELL);
-double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
-int replace(SHORT, SHORT, SHORT, SHORT);
+double get_slope(int, int, int, int, CELL, CELL);
+int replace(int, int, int, int);
 
 /* do_cum.c */
 int do_cum(void);
@@ -110,7 +109,7 @@ double mfd_pow(double, int);
 int find_pourpts(void);
 
 /* haf_side.c */
-int haf_basin_side(SHORT, SHORT, SHORT);
+int haf_basin_side(int, int, int);
 
 /* init_vars.c */
 int init_vars(int, char *[]);
@@ -130,7 +129,7 @@ int sg_factor(void);
 int len_slp_equ(double, double, double, int, int);
 
 /* slope_len.c */
-int slope_length(SHORT, SHORT, SHORT, SHORT);
+int slope_length(int, int, int, int);
 
 /* split_str.c */
 CELL split_stream(int, int, int[], int[], int, CELL, double, CELL);

raster/r.watershed/ram/close_maps2.c

@@ -29,7 +29,7 @@ int close_array_seg(void)
 	Rast_init_colors(&colors);
 	Rast_make_random_colors(&colors, 1, max);
 
-	if (max < 10000) {
+	if (max < 1000) {
 	    Rast_set_c_color((CELL) 0, 0, 0, 0, &colors);
 	    r = 1;
 	    incr = 0;
@@ -53,10 +53,6 @@ int close_array_seg(void)
 				    flag = 0;
 				    incr = -1;
 				}
-				if (r % 200 == 0)
-				    G_debug(5,
-					    "adjusting colors: r=%d\tof %d basins",
-					    r, max);
 			    }
 			}
 		    }

raster/r.watershed/ram/def_basin.c

@@ -5,7 +5,7 @@ CELL def_basin(int row, int col, CELL basin_num,
 {
     int r, rr, c, cc, ct, new_r[9], new_c[9];
     CELL downdir, direction, asp_value, value, new_elev;
-    SHORT oldupdir, riteflag, leftflag, thisdir;
+    int oldupdir, riteflag, leftflag, thisdir;
 
     for (;;) {
 	bas[SEG_INDEX(bas_seg, row, col)] = basin_num;
@@ -48,7 +48,7 @@ CELL def_basin(int row, int col, CELL basin_num,
 		    if (direction == drain[rr][cc]) {
 			thisdir = updrain[rr][cc];
 			switch (haf_basin_side
-				(oldupdir, (SHORT) downdir, thisdir)) {
+				(oldupdir, (int) downdir, thisdir)) {
 			case LEFT:
 			    overland_cells(r, c, basin_num, basin_num - 1,
 					   &new_elev);

raster/r.watershed/ram/do_astar.c

@@ -8,7 +8,7 @@ double get_slope2(CELL, CELL, double);
 int do_astar(void)
 {
     int count;
-    SHORT upr, upc, r, c, ct_dir;
+    int upr, upc, r, c, ct_dir;
     CELL alt_val, alt_nbr[8];
     CELL is_in_list, is_worked;
     int index_doer, index_up;
@@ -23,10 +23,10 @@ int do_astar(void)
     double slope[8];
     int skip_diag;
 
-    G_message(_("SECTION 2: A * Search."));
+    G_message(_("SECTION 2: A* Search."));
 
     for (ct_dir = 0; ct_dir < sides; ct_dir++) {
-	/* get r, c (r_nbr, c_nbr) for neighbours */
+	/* get r, c (upr, upc) for neighbours */
 	upr = nextdr[ct_dir];
 	upc = nextdc[ct_dir];
 	/* account for rare cases when ns_res != ew_res */
@@ -48,7 +48,7 @@ int do_astar(void)
     while (heap_size > 0) {
 	G_percent(count++, do_points, 1);
 
-	/* start with point with lowest elevation, in case of equal elevation
+	/* get point with lowest elevation, in case of equal elevation
 	 * of following points, oldest point = point added earliest */
 	index_doer = astar_pts[1];
 
@@ -66,13 +66,13 @@ int do_astar(void)
 
 	FLAG_SET(worked, r, c);
 
-	/* check neighbours */
+	/* check all neighbours, breadth first search */
 	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 	    /* get r, c (upr, upc) for this neighbour */
 	    upr = r + nextdr[ct_dir];
 	    upc = c + nextdc[ct_dir];
 	    slope[ct_dir] = alt_nbr[ct_dir] = 0;
-	    /* check that r, c are within region */
+	    /* check if r, c are within region */
 	    if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
 		index_up = SEG_INDEX(alt_seg, upr, upc);
 		is_in_list = FLAG_GET(in_list, upr, upc);
@@ -145,8 +145,8 @@ int cmp_pnt(CELL elea, CELL eleb, int addeda, int addedb)
     return 0;
 }
 
-/* new add point routine for min heap */
-int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
+/* add point routine for min heap */
+int add_pt(int r, int c, CELL ele, CELL downe)
 {
     FLAG_SET(in_list, r, c);
 
@@ -165,7 +165,7 @@ int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
     return 0;
 }
 
-/* new drop point routine for min heap */
+/* drop point routine for min heap */
 int drop_pt(void)
 {
     register int child, childr, parent;
@@ -182,10 +182,9 @@ int drop_pt(void)
     parent = 1;
 
     /* sift down: move hole back towards bottom of heap */
-    /* sift-down routine customised for A * Search logic */
 
     while ((child = GET_CHILD(parent)) <= heap_size) {
-	/* select child with lower ele, if equal, older child
+	/* select child with lower ele, if both are equal, older child
 	 * older child is older startpoint for flow path, important */
 	ele = alt[astar_pts[child]];
 	if (child < heap_size) {
@@ -193,21 +192,12 @@ int drop_pt(void)
 	    i = child + 3;
 	    while (childr <= heap_size && childr < i) {
 		eler = alt[astar_pts[childr]];
-		/* get smallest child */
 		if (cmp_pnt(eler, ele, heap_index[childr], heap_index[child])) {
 		    child = childr;
 		    ele = eler;
 		}
 		childr++;
 	    }
-	    /* break if childr > last entry? that saves sifting up again
-	     * OTOH, this is another comparison
-	     * we have a max heap height of 20: log(INT_MAX)/log(n children per node)
-	     * that would give us in the worst case 20*2 additional comparisons with 3 children
-	     * the last entry will never go far up again, less than half the way
-	     * so the additional comparisons for going all the way down
-	     * and then a bit up again are likely less than 20*2 */
-	    /* find the error in this reasoning */
 	}
 
 	/* move hole down */
@@ -271,7 +261,7 @@ int sift_up(int start, CELL ele)
 }
 
 double
-get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+get_slope(int r, int c, int downr, int downc, CELL ele, CELL downe)
 {
     double slope;
 
@@ -290,15 +280,15 @@ get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
 
 double get_slope2(CELL ele, CELL up_ele, double dist)
 {
-    if (ele == up_ele)
-	return 0.5 / dist;
+    if (ele >= up_ele)
+	return 0.0;
     else
 	return (double)(up_ele - ele) / dist;
 }
 
 /* replace is unused */
 int replace(			/* ele was in there */
-	       SHORT upr, SHORT upc, SHORT r, SHORT c)
+	       int upr, int upc, int r, int c)
 /* CELL ele;  */
 {
     int now, heap_run;

raster/r.watershed/ram/do_cum.c

@@ -6,11 +6,11 @@
 
 int do_cum(void)
 {
-    SHORT r, c, dr, dc;
+    int r, c, dr, dc;
     CELL is_swale, value, valued, aspect;
     int killer, threshold, count;
-    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
-    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    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 };
     int this_index, down_index;
 
     G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
@@ -110,9 +110,9 @@ int do_cum_mfd(void)
     double dx, dy;
     CELL ele, ele_nbr, aspect, is_worked;
     double prop, max_acc;
-    int workedon, edge;
-    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
-    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    int workedon, edge, flat;
+    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 };
     int this_index, down_index, nbr_index;
 
     G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
@@ -174,6 +174,7 @@ int do_cum_mfd(void)
 	    ele = alt[this_index];
 	    is_null = 0;
 	    edge = 0;
+	    flat = 1;
 	    /* this loop is needed to get the sum of weights */
 	    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 		/* get r, c (r_nbr, c_nbr) for neighbours */
@@ -184,6 +185,9 @@ int do_cum_mfd(void)
 		if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
 		    c_nbr < ncols) {
 
+		    if (dr == r_nbr && dc == c_nbr)
+			np_side = ct_dir;
+
 		    nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
 
 		    /* check for swale or stream cells */
@@ -191,12 +195,15 @@ int do_cum_mfd(void)
 		    if (is_swale)
 			swale_cells++;
 		    valued = wat[nbr_index];
-		    if ((ABS(valued) + 0.5) >= threshold)
+		    ele_nbr = alt[nbr_index];
+		    if ((ABS(valued) + 0.5) >= threshold  &&
+		        ct_dir != np_side && ele_nbr > ele)
 			stream_cells++;
 
 		    is_worked = FLAG_GET(worked, r_nbr, c_nbr);
 		    if (is_worked == 0) {
-			ele_nbr = alt[nbr_index];
+			if (ele_nbr != ele)
+			    flat = 0;
 			is_null = Rast_is_c_null_value(&ele_nbr);
 			edge = is_null;
 			if (!is_null && ele_nbr <= ele) {
@@ -223,8 +230,6 @@ int do_cum_mfd(void)
 			    }
 			}
 		    }
-		    if (dr == r_nbr && dc == c_nbr)
-			np_side = ct_dir;
 		}
 		else
 		    edge = 1;
@@ -260,7 +265,6 @@ int do_cum_mfd(void)
 	    if (mfd_cells > 1) {
 		prop = 0.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];
 		    c_nbr = c + nextdc[ct_dir];
 
@@ -273,7 +277,7 @@ int do_cum_mfd(void)
 			    nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
 
 			    weight[ct_dir] = weight[ct_dir] / sum_weight;
-			    /* check everything sums up to 1.0 */
+			    /* check everything adds up to 1.0 */
 			    prop += weight[ct_dir];
 
 			    valued = wat[nbr_index];
@@ -344,8 +348,8 @@ int do_cum_mfd(void)
 	    }
 	    /* start new stream */
 	    value = ABS(value) + 0.5;
-	    if (!is_swale && (int)value >= threshold && stream_cells < 3 &&
-		swale_cells < 1) {
+	    if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
+		swale_cells < 1 && !flat) {
 		FLAG_SET(swale, r, c);
 		is_swale = 1;
 	    }

raster/r.watershed/ram/haf_side.c

@@ -1,8 +1,8 @@
 #include "Gwater.h"
 
-int haf_basin_side(SHORT updir, SHORT downdir, SHORT thisdir)
+int haf_basin_side(int updir, int downdir, int thisdir)
 {
-    SHORT newup, newthis;
+    int newup, newthis;
 
     newup = updir - downdir;
     if (newup < 0)

raster/r.watershed/ram/init_vars.c

@@ -9,7 +9,7 @@ int ele_round(double);
 
 int init_vars(int argc, char *argv[])
 {
-    SHORT r, c;
+    int r, c;
     CELL *buf, alt_value, wat_value, asp_value, block_value;
     DCELL dvalue;
     void *elebuf, *ptr;
@@ -34,8 +34,6 @@ int init_vars(int argc, char *argv[])
     c_fac = 5;
     abs_acc = 0;
     ele_scale = 1;
-    for (r = 0; r < argc; r++)
-	G_debug(0, "%s", argv[r]);
 
     for (r = 1; r < argc; r++) {
 	if (sscanf(argv[r], "elevation=%s", ele_name) == 1)

raster/r.watershed/ram/main.c

@@ -27,7 +27,7 @@ struct Cell_head window;
 int mfd, c_fac, abs_acc, ele_scale;
 int *heap_index, heap_size;
 int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
-SHORT nrows, ncols;
+int nrows, ncols;
 double half_res, diag, max_length, dep_slope;
 int bas_thres, tot_parts;
 CELL n_basins;
@@ -44,11 +44,11 @@ int ril_fd;
 double *s_l, *s_g, *l_s;
 CELL one, zero;
 double ril_value, d_one, d_zero;
-SHORT sides;
-SHORT drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
-SHORT updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
-SHORT nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
-SHORT nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+int sides;
+int drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
+int updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
+int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
 char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
 char run_name[GNAME_MAX], ob_name[GNAME_MAX];
 char ril_name[GNAME_MAX], dep_name[GNAME_MAX];

raster/r.watershed/ram/no_stream.c

@@ -7,7 +7,7 @@ int no_stream(int row, int col, CELL basin_num,
     double slope;
     CELL downdir, asp_value, hih_ele, new_ele, aspect;
     DCELL dvalue, max_drain;	/* flow acc is now DCELL */
-    SHORT updir, riteflag, leftflag, thisdir;
+    int updir, riteflag, leftflag, thisdir;
 
     while (1) {
 	bas[SEG_INDEX(bas_seg, row, col)] = basin_num;
@@ -67,7 +67,7 @@ int no_stream(int row, int col, CELL basin_num,
 			if (aspect == drain[rr][cc]) {
 			    thisdir = updrain[rr][cc];
 			    switch (haf_basin_side
-			            (updir, (SHORT) downdir, thisdir)) {
+			            (updir, (int) downdir, thisdir)) {
 			    case LEFT:
 				overland_cells(r, c, basin_num, basin_num - 1,
 					       &new_ele);

raster/r.watershed/ram/slope_len.c

@@ -1,6 +1,6 @@
 #include "Gwater.h"
 
-int slope_length(SHORT r, SHORT c, SHORT dr, SHORT dc)
+int slope_length(int r, int c, int dr, int dc)
 {
     CELL top_alt, bot_alt, asp_value;
     double res, top_ls, bot_ls;

raster/r.watershed/ram/split_str.c

@@ -6,8 +6,8 @@ split_stream(int row, int col, int new_r[], int new_c[], int ct,
 {
     CELL downdir, old_basin, new_elev, aspect;
     double slope, easting, northing;
-    SHORT doit, ctr, updir, splitdir[9];
-    SHORT thisdir, leftflag, riteflag;
+    int doit, ctr, updir, splitdir[9];
+    int thisdir, leftflag, riteflag;
     int r, c, rr, cc;
 
     new_elev = 0;
@@ -35,7 +35,7 @@ split_stream(int row, int col, int new_r[], int new_c[], int ct,
 		    if (doit) {
 			thisdir = updrain[rr][cc];
 			switch (haf_basin_side
-				(updir, (SHORT) downdir, thisdir)) {
+				(updir, (int) downdir, thisdir)) {
 			case LEFT:
 			    overland_cells(r, c, basin_num, basin_num - 1,
 					   &new_elev);

raster/r.watershed/seg/Gwater.h

@@ -15,7 +15,6 @@
 
 #define AR_SIZE			16
 #define AR_INCR			16
-#define SHORT			int
 #define NOMASK			1
 #define MIN_SLOPE		.00001
 #define MIN_GRADIENT_DEGREES	1
@@ -45,7 +44,7 @@
 
 #define POINT       struct points
 POINT {
-    SHORT r, c; /* , downr, downc */
+    int r, c; /* , downr, downc */
     char asp;      /* drainage direction */
     char guessed;   /* accumulation will likely be an underestimate */
 };
@@ -77,7 +76,7 @@ extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 extern CELL n_basins;
 extern OC_STACK *ocs;
 extern int ocs_alloced;
-extern SHORT nrows, ncols;
+extern int nrows, ncols;
 extern double half_res, diag, max_length, dep_slope;
 extern int bas_thres, tot_parts;
 extern SSEG astar_pts;
@@ -88,11 +87,11 @@ extern DSEG slp, s_l, s_g, l_s, ril;
 extern double segs_mb;
 extern char zero, one;
 extern double ril_value, d_zero, d_one;
-extern SHORT sides;
-extern SHORT drain[3][3];
-extern SHORT updrain[3][3];
-extern SHORT nextdr[8];
-extern SHORT nextdc[8];
+extern int sides;
+extern int drain[3][3];
+extern int updrain[3][3];
+extern int nextdr[8];
+extern int nextdc[8];
 extern char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
 extern char run_name[GNAME_MAX], ob_name[GNAME_MAX];
 extern char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
@@ -140,14 +139,14 @@ CELL def_basin(int, int, CELL, double, CELL);
 
 /* do_astar.c */
 int do_astar(void);
-int add_pt(SHORT, SHORT, CELL, char, int);
+int add_pt(int, int, CELL, char, int);
 HEAP_PNT drop_pt(void);
 int sift_up(int, HEAP_PNT);
 int sift_up_mem(int, HEAP_PNT);
 int sift_down_disk(void);
 int fill_mem_heap(void);
 int cmp_pnt(HEAP_PNT *a, HEAP_PNT *b);
-double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+double get_slope(int, int, int, int, CELL, CELL);
 
 /* do_cum.c */
 int do_cum(void);
@@ -161,7 +160,7 @@ int do_stream(void);
 int find_pourpts(void);
 
 /* haf_side.c */
-int haf_basin_side(SHORT, SHORT, SHORT);
+int haf_basin_side(int, int, int);
 
 /* init_vars.c */
 int init_vars(int, char *[]);
@@ -178,7 +177,7 @@ int sg_factor(void);
 int len_slp_equ(double, double, double, int, int);
 
 /* slope_len.c */
-int slope_length(SHORT, SHORT, SHORT, SHORT);
+int slope_length(int, int, int, int);
 
 /* split_str.c */
 CELL split_stream(int, int, int[], int[], int, CELL, double, CELL);

raster/r.watershed/seg/bseg_read.c

@@ -8,7 +8,6 @@ int bseg_read_cell(BSEG * bseg, char *map_name, char *mapset)
     int row, nrows;
     int col, ncols;
     int map_fd;
-    char msg[100];
     CELL *buffer;
     char cbuf;
 

raster/r.watershed/seg/close_maps2.c

@@ -35,7 +35,7 @@ int close_array_seg(void)
 	Rast_init_colors(&colors);
 	Rast_make_random_colors(&colors, 1, max);
 
-	if (max < 10000) {
+	if (max < 1000) {
 	    Rast_set_c_color((CELL) 0, 0, 0, 0, &colors);
 	    r = 1;
 	    incr = 0;

raster/r.watershed/seg/cseg_read.c

@@ -9,7 +9,6 @@ int cseg_read_cell(CSEG * cseg, char *map_name, char *mapset)
 {
     int row, nrows;
     int map_fd;
-    char msg[100];
     CELL *buffer;
 
     cseg->name = NULL;
@@ -23,9 +22,8 @@ int cseg_read_cell(CSEG * cseg, char *map_name, char *mapset)
 	if (segment_put_row(&(cseg->seg), buffer, row) < 0) {
 	    G_free(buffer);
 	    Rast_close(map_fd);
-	    sprintf(msg, "%s(): unable to segment put row for [%s] in [%s]",
+	    G_warning("%s(): unable to segment put row for [%s] in [%s]",
 		    me, map_name, mapset);
-	    G_warning(msg);
 	    return (-1);
 	}
     }

raster/r.watershed/seg/def_basin.c

@@ -7,7 +7,7 @@ def_basin(int row, int col, CELL basin_num, double stream_length,
     int r, rr, c, cc, ct, new_r[9], new_c[9];
     CELL downdir, direction, asp_value, value, new_elev;
     char cvalue;
-    SHORT oldupdir, riteflag, leftflag, thisdir;
+    int oldupdir, riteflag, leftflag, thisdir;
 
     for (;;) {
 	cseg_put(&bas, &basin_num, row, col);
@@ -49,7 +49,7 @@ def_basin(int row, int col, CELL basin_num, double stream_length,
 		    if (direction == drain[rr][cc]) {
 			thisdir = updrain[rr][cc];
 			switch (haf_basin_side
-				(oldupdir, (SHORT) downdir, thisdir)) {
+				(oldupdir, (int) downdir, thisdir)) {
 			case LEFT:
 			    overland_cells(r, c, basin_num, basin_num - 1,
 					   &new_elev);

raster/r.watershed/seg/do_astar.c

@@ -8,7 +8,7 @@ double get_slope2(CELL, CELL, double);
 int do_astar(void)
 {
     int doer, count;
-    SHORT upr, upc, r = -1, c = -1, ct_dir;
+    int upr, upc, r = -1, c = -1, ct_dir;
     CELL alt_val, alt_nbr[8], alt_up;
     CELL asp_val;
     char flag_value, is_in_list, is_worked;
@@ -24,11 +24,12 @@ int do_astar(void)
     double dx, dy, dist_to_nbr[8], ew_res, ns_res;
     double slope[8];
     int skip_diag;
+    int count_edge = 0, count_diag = 0, count_edge_sink = 0, count_diag_sink = 0;
 
     G_message(_("SECTION 2: A* Search."));
 
     for (ct_dir = 0; ct_dir < sides; ct_dir++) {
-	/* get r, c (r_nbr, c_nbr) for neighbours */
+	/* get r, c (upr, upc) for neighbours */
 	upr = nextdr[ct_dir];
 	upc = nextdc[ct_dir];
 	/* account for rare cases when ns_res != ew_res */
@@ -51,13 +52,10 @@ int do_astar(void)
 
     doer = do_points - 1;
 
-    /* A * Search
-     * determine downhill path through all not masked cells and
-     * set drainage direction */
+    /* A* Search: search uphill, get downhill paths */
     while (heap_size > 0) {
 	G_percent(count++, do_points, 1);
 
-	/* drop next point from heap */
 	heap_p = drop_pt();
 
 	r = heap_p.pnt.r;
@@ -72,13 +70,13 @@ int do_astar(void)
 	    upr = r + nextdr[ct_dir];
 	    upc = c + nextdc[ct_dir];
 	    slope[ct_dir] = alt_nbr[ct_dir] = 0;
-	    /* check that upr, upc are within region */
+	    /* check if upr, upc are within region */
 	    if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
-		/* avoid diagonal flow direction bias */
 		bseg_get(&bitflags, &flag_value, upr, upc);
 		is_in_list = FLAG_GET(flag_value, INLISTFLAG);
 		is_worked = FLAG_GET(flag_value, WORKEDFLAG);
 		skip_diag = 0;
+		/* avoid diagonal flow direction bias */
 		if (!is_worked) {
 		    cseg_get(&alt, &alt_up, upr, upc);
 		    alt_nbr[ct_dir] = alt_up;
@@ -105,17 +103,31 @@ int do_astar(void)
 		    }
 		}
 
-		/* put neighbour in search list if not yet in */
+		/* add neighbour as new point if not in the list */
 		if (is_in_list == 0 && skip_diag == 0) {
-		    cseg_get(&alt, &alt_up, upr, upc);
-		    /* flow direction is set here */
+		    /* set flow direction */
 		    asp_val = drain[upr - r + 1][upc - c + 1];
-		    add_pt(upr, upc, alt_up, asp_val, 0);
+		    add_pt(upr, upc, alt_nbr[ct_dir], asp_val, 0);
 		    cseg_put(&asp, &asp_val, upr, upc);
+
+
+		    if (alt_nbr[ct_dir] < alt_val) {
+			if (ct_dir < 4)
+			    count_edge_sink++;
+			else
+			    count_diag_sink++;
+		    }
+		    /* includes flat areas */
+		    else {
+			if (ct_dir < 4)
+			    count_edge++;
+			else
+			    count_diag++;
+		    }
 		}
-		/* in list, not worked. is it edge ? */
 		else if (is_in_list && is_worked == 0 &&
 			 FLAG_GET(flag_value, EDGEFLAG)) {
+		    /* neighbour is edge in list, not yet worked */
 		    cseg_get(&asp, &asp_val, upr, upc);
 		    if (asp_val < 0) {
 			/* adjust flow direction for edge cell */
@@ -126,6 +138,7 @@ int do_astar(void)
 		}
 	    }
 	}
+	/* add astar points to sorted list for flow accumulation */
 	cseg_get(&asp, &asp_val, r, c);
 	heap_p.pnt.asp = asp_val;
 	seg_put(&astar_pts, (char *)&heap_p.pnt, 0, doer);
@@ -143,6 +156,11 @@ int do_astar(void)
 
     G_percent(count, do_points, 1);	/* finish it */
 
+    G_debug(1, "edge direction: %d (%.2f%%)", count_edge, (double) 100. * count_edge / (count_edge + count_diag));
+    G_debug(1, "diag direction: %d (%.2f%%)", count_diag, (double) 100. * count_diag / (count_edge + count_diag));
+    G_debug(1, "edge out of depression: %d (%.2f%%)", count_edge_sink, (double) 100. * count_edge_sink / (count_edge_sink + count_diag_sink));
+    G_debug(1, "diag out of depression: %d (%.2f%%)", count_diag_sink, (double) 100. * count_diag_sink / (count_edge_sink + count_diag_sink));
+
     return 0;
 }
 
@@ -160,7 +178,7 @@ int cmp_pnt(HEAP_PNT * a, HEAP_PNT * b)
 }
 
 /* add point routine for min heap */
-int add_pt(SHORT r, SHORT c, CELL ele, char asp, int is_edge)
+int add_pt(int r, int c, CELL ele, char asp, int is_edge)
 {
     HEAP_PNT heap_p;
     char flag_value;
@@ -204,7 +222,7 @@ HEAP_PNT drop_pt(void)
     seg_get(&search_heap, (char *)&last_p, 0, heap_size);
     seg_get(&search_heap, (char *)&root_p, 0, 1);
 
-    /* sift down */
+    /* sift down: move hole back towards bottom of heap */
     parent = 1;
     while ((child = GET_CHILD(parent)) < heap_size) {
 	/* select child with lower ele, if both are equal, older child
@@ -215,7 +233,6 @@ HEAP_PNT drop_pt(void)
 	    childr = child + 1;
 	    i = child + 4;
 	    while (childr < i && childr < heap_size) {
-		/* get smallest child */
 		seg_get(&search_heap, (char *)&childr_p, 0, childr);
 		if (cmp_pnt(&childr_p, &child_p)) {
 		    child = childr;
@@ -245,7 +262,7 @@ HEAP_PNT drop_pt(void)
 /* standard sift-up routine for d-ary min heap */
 int sift_up(int start, HEAP_PNT child_p)
 {
-    int parent, child;		/* parentp, */
+    int parent, child;
     HEAP_PNT heap_p;
 
     child = start;
@@ -271,7 +288,7 @@ int sift_up(int start, HEAP_PNT child_p)
 }
 
 double
-get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+get_slope(int r, int c, int downr, int downc, CELL ele, CELL downe)
 {
     double slope;
 
@@ -288,8 +305,8 @@ get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
 
 double get_slope2(CELL ele, CELL up_ele, double dist)
 {
-    if (ele == up_ele)
-	return 0.5 / dist;
+    if (ele >= up_ele)
+	return 0.0;
     else
 	return (double)(up_ele - ele) / dist;
 }

raster/r.watershed/seg/do_cum.c

@@ -7,14 +7,14 @@
 
 int do_cum(void)
 {
-    SHORT r, c, dr, dc;
+    int r, c, dr, dc;
     CELL asp_val, asp_val_down;
     char is_swale, this_flag_value, flag_value;
     DCELL value, valued;
     POINT point;
     int killer, threshold;
-    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
-    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    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 };
     WAT_ALT wa, wadown;
 
     G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
@@ -29,11 +29,11 @@ int do_cum(void)
 	r = point.r;
 	c = point.c;
 	asp_val = point.asp;
-	/* skip user-defined depressions */
 	if (asp_val) {
 	    dr = r + asp_r[ABS(asp_val)];
 	    dc = c + asp_c[ABS(asp_val)];
 	}
+	/* skip user-defined depressions */
 	else
 	    dr = dc = -1;
 
@@ -130,10 +130,10 @@ int do_cum_mfd(void)
     double dx, dy;
     CELL ele, *ele_nbr, asp_val, asp_val_down;
     double prop, max_acc;
-    int workedon, edge, is_swale;
+    int workedon, edge, is_swale, flat;
     char *flag_nbr, this_flag_value, flag_value;
-    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
-    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    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 };
 
     G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
     G_debug(1, "MFD convergence factor set to %d.", c_fac);
@@ -206,6 +206,7 @@ int do_cum_mfd(void)
 	    ele = wa.ele;
 	    is_null = 0;
 	    edge = 0;
+	    flat = 1;
 	    /* this loop is needed to get the sum of weights */
 	    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 		/* get r, c (r_nbr, c_nbr) for neighbours */
@@ -224,6 +225,9 @@ int do_cum_mfd(void)
 		if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
 		    c_nbr < ncols) {
 
+		    if (dr == r_nbr && dc == c_nbr)
+			np_side = ct_dir;
+
 		    bseg_get(&bitflags, &flag_nbr[ct_dir], r_nbr, c_nbr);
 		    seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
 		    wat_nbr[ct_dir] = wa.wat;
@@ -233,11 +237,15 @@ int do_cum_mfd(void)
 		    is_swale = FLAG_GET(flag_nbr[ct_dir], SWALEFLAG);
 		    if (is_swale)
 			swale_cells++;
-		    if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold)
+		    if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold &&
+		        ct_dir != np_side && ele_nbr[ct_dir] > ele)
 			stream_cells++;
 
 		    if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
 
+			if (ele_nbr[ct_dir] != ele)
+			    flat = 0;
+
 			edge = is_null = FLAG_GET(flag_nbr[ct_dir], NULLFLAG);
 			if (!is_null && ele_nbr[ct_dir] <= ele) {
 			    if (ele_nbr[ct_dir] < ele) {
@@ -263,8 +271,6 @@ int do_cum_mfd(void)
 			    }
 			}
 		    }
-		    if (dr == r_nbr && dc == c_nbr)
-			np_side = ct_dir;
 		}
 		else
 		    edge = 1;
@@ -302,7 +308,6 @@ int do_cum_mfd(void)
 	    if (mfd_cells > 1) {
 		prop = 0.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];
 		    c_nbr = c + nextdc[ct_dir];
 
@@ -313,7 +318,7 @@ int do_cum_mfd(void)
 			if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
 
 			    weight[ct_dir] = weight[ct_dir] / sum_weight;
-			    /* check everything sums up to 1.0 */
+			    /* check everything adds up to 1.0 */
 			    prop += weight[ct_dir];
 
 			    if (value > 0) {
@@ -393,9 +398,8 @@ int do_cum_mfd(void)
 	    is_swale = FLAG_GET(this_flag_value, SWALEFLAG);
 	    /* start new stream */
 	    value = ABS(value) + 0.5;
-	    /* can use stream_cells < 4 only for highres, nsres and ewres < 30 m? */
-	    if (!is_swale && (int)value >= threshold && stream_cells < 3 &&
-		swale_cells < 1) {
+	    if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
+		swale_cells < 1 && !flat) {
 		FLAG_SET(this_flag_value, SWALEFLAG);
 		is_swale = 1;
 	    }
@@ -426,9 +430,7 @@ int do_cum_mfd(void)
 	G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
 		  workedon, do_points);
 
-    if (!st_flag) {
-	seg_close(&astar_pts);
-    }
+    seg_close(&astar_pts);
     
     G_free(dist_to_nbr);
     G_free(weight);

raster/r.watershed/seg/dseg_read.c

@@ -9,7 +9,6 @@ int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
 {
     int row, nrows, ncols;
     int map_fd;
-    char msg[100];
     double *dbuffer;
 
     dseg->name = NULL;
@@ -24,9 +23,8 @@ int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
 	if (segment_put_row(&(dseg->seg), (DCELL *) dbuffer, row) < 0) {
 	    G_free(dbuffer);
 	    Rast_close(map_fd);
-	    sprintf(msg, "%s(): unable to segment put row for [%s] in [%s]",
+	    G_warning("%s(): unable to segment put row for [%s] in [%s]",
 		    me, map_name, mapset);
-	    G_warning(msg);
 	    return (-1);
 	}
     }

raster/r.watershed/seg/haf_side.c

@@ -1,8 +1,8 @@
 #include "Gwater.h"
 
-int haf_basin_side(SHORT updir, SHORT downdir, SHORT thisdir)
+int haf_basin_side(int updir, int downdir, int thisdir)
 {
-    SHORT newup, newthis;
+    int newup, newthis;
 
     newup = updir - downdir;
     if (newup < 0)

raster/r.watershed/seg/init_vars.c

@@ -9,7 +9,7 @@ int ele_round(double);
 
 int init_vars(int argc, char *argv[])
 {
-    SHORT r, c;
+    int r, c;
     int ele_fd, fd, wat_fd;
     int seg_rows, seg_cols, num_cseg_total, num_open_segs, num_open_array_segs;
 
@@ -116,13 +116,6 @@ int init_vars(int argc, char *argv[])
     /* do RUSLE */
     if (er_flag)
 	tot_parts++;
-    /* do stream extraction */
-    if (er_flag || seg_flag || bas_flag || haf_flag) {
-	st_flag = 1;
-	/* separate stream extraction only needed for MFD */
-	if (mfd)
-	    tot_parts++;
-    }
     /* define basins */
     if (seg_flag || bas_flag || haf_flag)
 	tot_parts++;
@@ -288,7 +281,7 @@ int init_vars(int argc, char *argv[])
 			else if (wat_map_type == FCELL_TYPE) {
 			    wat_value = *((FCELL *)ptr);
 			}
-			else if (ele_map_type == DCELL_TYPE) {
+			else if (wat_map_type == DCELL_TYPE) {
 			    wat_value = *((DCELL *)ptr);
 			}
 		    }
@@ -324,8 +317,7 @@ int init_vars(int argc, char *argv[])
 	G_free(watbuf);
     }
 
-    if (do_points < nrows * ncols)
-	MASK_flag = 1;
+    MASK_flag = (do_points < nrows * ncols);
     
     /* depression: drainage direction will be set to zero later */
     if (pit_flag) {
@@ -388,7 +380,11 @@ int init_vars(int argc, char *argv[])
     G_debug(1, "open segments for A* points");
     /* rounded down power of 2 */
     seg_cols = (int) (pow(2, (int)(log(num_open_segs / 8.0) / log(2) + 0.1)) + 0.1);
+    if (seg_cols < 2)
+	seg_cols = 2;
     num_open_array_segs = num_open_segs / seg_cols;
+    if (num_open_array_segs == 0)
+	num_open_array_segs = 1;
     /* n cols in segment */
     seg_cols *= seg_rows * seg_rows;
     /* n segments in row */
@@ -409,7 +405,11 @@ int init_vars(int argc, char *argv[])
     G_debug(1, "open segments for A* search heap");
     /* rounded down power of 2 */
     seg_cols = (int) (pow(2, (int)(log(num_open_segs / 8.0) / log(2) + 0.1)) + 0.1);
+    if (seg_cols < 2)
+	seg_cols = 2;
     num_open_array_segs = num_open_segs / seg_cols;
+    if (num_open_array_segs == 0)
+	num_open_array_segs = 1;
     /* n cols in segment */
     seg_cols *= seg_rows * seg_rows;
     /* n segments in row */

raster/r.watershed/seg/main.c

@@ -32,7 +32,7 @@ int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 CELL n_basins;
 OC_STACK *ocs;
 int ocs_alloced;
-SHORT nrows, ncols;
+int nrows, ncols;
 double half_res, diag, max_length, dep_slope;
 int bas_thres, tot_parts;
 SSEG astar_pts;
@@ -43,11 +43,11 @@ DSEG slp, s_l, s_g, l_s, ril;
 double segs_mb;
 char zero, one;
 double ril_value, d_zero, d_one;
-SHORT sides;
-SHORT drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
-SHORT updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
-SHORT nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
-SHORT nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+int sides;
+int drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
+int updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
+int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
 char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
 char run_name[GNAME_MAX], ob_name[GNAME_MAX];
 char ril_name[GNAME_MAX], dep_name[GNAME_MAX];

raster/r.watershed/seg/no_stream.c

@@ -8,7 +8,7 @@ no_stream(int row, int col, CELL basin_num, double stream_length,
     double slope;
     CELL downdir, asp_value, hih_ele, new_ele, aspect, value;
     DCELL dvalue, max_drain;	/* flow acc is now DCELL */
-    SHORT updir, riteflag, leftflag, thisdir;
+    int updir, riteflag, leftflag, thisdir;
     WAT_ALT wa;
 
     while (1) {
@@ -73,7 +73,7 @@ no_stream(int row, int col, CELL basin_num, double stream_length,
 			if (aspect == drain[rr][cc]) {
 			    thisdir = updrain[rr][cc];
 			    switch (haf_basin_side(updir,
-					       (SHORT) downdir,
+					       (int) downdir,
 					       thisdir)) {
 			    case RITE:
 				overland_cells(r, c, basin_num, basin_num,

raster/r.watershed/seg/slope_len.c

@@ -1,6 +1,6 @@
 #include "Gwater.h"
 
-int slope_length(SHORT r, SHORT c, SHORT dr, SHORT dc)
+int slope_length(int r, int c, int dr, int dc)
 {
     CELL top_alt, bot_alt, asp_value, ridge;
     double res, top_ls, bot_ls;

raster/r.watershed/seg/split_str.c

@@ -6,8 +6,8 @@ split_stream(int row, int col, int new_r[], int new_c[], int ct,
 {
     CELL downdir, old_basin, new_elev, aspect;
     double slope, easting, northing;
-    SHORT doit, ctr, updir, splitdir[9];
-    SHORT thisdir, leftflag, riteflag;
+    int doit, ctr, updir, splitdir[9];
+    int thisdir, leftflag, riteflag;
     int r, c, rr, cc;
 
     new_elev = 0;
@@ -35,7 +35,7 @@ split_stream(int row, int col, int new_r[], int new_c[], int ct,
 		    if (doit) {
 			thisdir = updrain[rr][cc];
 			switch (haf_basin_side
-				(updir, (SHORT) downdir, thisdir)) {
+				(updir, (int) downdir, thisdir)) {
 			case LEFT:
 			    overland_cells(r, c, basin_num, basin_num - 1,
 					   &new_elev);