synchronized to ram mode

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

raster/r.watershed/seg/Gwater.h

@@ -36,7 +36,7 @@
 #define NULLFLAG         0      /* elevation is NULL */
 #define EDGEFLAG         1      /* edge cell */
 #define INLISTFLAG       2      /* in open A* list */
-#define WORKEDFLAG       3      /* in closed A* list/ accumulation done/streams done */
+#define WORKEDFLAG       3      /* in closed A* list/ accumulation done */
 #define SWALEFLAG        4      /* swale */
 #define PITFLAG          5      /* user-defined real depression */
 #define RUSLEBLOCKFLAG   6      /* is RUSLE block */

raster/r.watershed/seg/close_maps.c

@@ -72,7 +72,7 @@ int close_maps(void)
 
 	if (min < 0) {
 	    if (min < (-stddev - 1)) {
-		clr_min = min;
+		clr_min = min - 1;
 		clr_max = -stddev - 1;
 		Rast_add_d_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
 					  0, 0, &colors);
@@ -117,7 +117,7 @@ int close_maps(void)
 
 	if (max > 0 && max > stddev + 1) {
 	    clr_min = stddev + 1;
-	    clr_max = max;
+	    clr_max = max + 1;
 	    Rast_add_d_color_rule(&clr_min, 0, 0, 0, &clr_max, 0, 0, 0,
 				      &colors);
 	}

raster/r.watershed/seg/close_maps2.c

@@ -8,7 +8,8 @@ int close_array_seg(void)
     struct Colors colors;
     int incr, max, red, green, blue, rd, gr, bl, flag;
     int c, r, map_fd;
-    CELL *cellrow, value;
+    CELL *cellrow;
+    /* CELL value; */
     char cvalue;
     CSEG *theseg;
 

raster/r.watershed/seg/do_astar.c

@@ -74,18 +74,19 @@ int do_astar(void)
 	    slope[ct_dir] = alt_nbr[ct_dir] = 0;
 	    /* check that upr, upc are within region */
 	    if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
-		/* slope */
+		/* 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;
 		    slope[ct_dir] =
 			get_slope2(alt_val, alt_nbr[ct_dir],
 				   dist_to_nbr[ct_dir]);
+		}
+		if (!is_in_list) {
 		    if (ct_dir > 3 && slope[ct_dir] > 0) {
 			if (slope[nbr_ew[ct_dir]] > 0) {
 			    /* slope to ew nbr > slope to center */
@@ -103,8 +104,8 @@ int do_astar(void)
 			}
 		    }
 		}
+
 		/* put neighbour in search list if not yet in */
-		bseg_get(&bitflags, &flag_value, upr, upc);
 		if (is_in_list == 0 && skip_diag == 0) {
 		    cseg_get(&alt, &alt_up, upr, upc);
 		    /* flow direction is set here */
@@ -113,7 +114,7 @@ int do_astar(void)
 		    cseg_put(&asp, &asp_val, upr, upc);
 		}
 		/* in list, not worked. is it edge ? */
-		else if (is_in_list == 1 && is_worked == 0 &&
+		else if (is_in_list && is_worked == 0 &&
 			 FLAG_GET(flag_value, EDGEFLAG)) {
 		    cseg_get(&asp, &asp_val, upr, upc);
 		    if (asp_val < 0) {

raster/r.watershed/seg/do_cum.c

@@ -126,12 +126,12 @@ int do_cum_mfd(void)
     /* MFD */
     int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
     double *dist_to_nbr, *weight, sum_weight, max_weight;
-    int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
+    int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side, max_side;
     double dx, dy;
-    CELL ele, *ele_nbr, asp_val;
+    CELL ele, *ele_nbr, asp_val, asp_val_down;
     double prop, max_acc;
-    int edge;
-    char workedon, *flag_nbr, this_flag_value;
+    int workedon, edge, is_swale;
+    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 };
 
@@ -172,11 +172,11 @@ int do_cum_mfd(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;
 
@@ -192,12 +192,6 @@ int do_cum_mfd(void)
 	FLAG_UNSET(this_flag_value, WORKEDFLAG);
 	
 	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) {
-	    /* do not distribute flow along edges, this causes artifacts */
-	    if (FLAG_GET(this_flag_value, EDGEFLAG)) {
-		bseg_put(&bitflags, &this_flag_value, r, c);
-		continue;
-	    }
-
 	    r_max = dr;
 	    c_max = dc;
 
@@ -231,11 +225,18 @@ int do_cum_mfd(void)
 		    c_nbr < ncols) {
 
 		    bseg_get(&bitflags, &flag_nbr[ct_dir], r_nbr, c_nbr);
-		    if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
+		    seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
+		    wat_nbr[ct_dir] = wa.wat;
+		    ele_nbr[ct_dir] = wa.ele;
+
+		    /* check for swale or stream cells */
+		    is_swale = FLAG_GET(flag_nbr[ct_dir], SWALEFLAG);
+		    if (is_swale)
+			swale_cells++;
+		    if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold)
+			stream_cells++;
 
-			seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
-			wat_nbr[ct_dir] = wa.wat;
-			ele_nbr[ct_dir] = wa.ele;
+		    if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
 
 			edge = is_null = FLAG_GET(flag_nbr[ct_dir], NULLFLAG);
 			if (!is_null && ele_nbr[ct_dir] <= ele) {
@@ -270,10 +271,13 @@ int do_cum_mfd(void)
 		if (edge)
 		    break;
 	    }
-	    /* do not distribute flow along edges, this causes artifacts */
+	    /* do not continue streams along edges, this causes artifacts */
 	    if (edge) {
-		G_debug(3, "edge: should not get to here");
-		bseg_put(&bitflags, &this_flag_value, r, c);
+		is_swale = FLAG_GET(this_flag_value, SWALEFLAG);
+		if (is_swale && asp_val > 0) {
+		    asp_val = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+		    cseg_put(&asp, &asp_val, r, c);
+		}
 		continue;
 	    }
 
@@ -288,10 +292,12 @@ int do_cum_mfd(void)
 		mfd_cells++;
 		sum_weight += max_weight;
 		weight[np_side] = max_weight;
+		max_side = np_side;
 	    }
 
 	    /* set flow accumulation for neighbours */
 	    max_acc = -1;
+	    max_side = np_side;
 
 	    if (mfd_cells > 1) {
 		prop = 0.0;
@@ -332,6 +338,7 @@ int do_cum_mfd(void)
 				max_acc = ABS(wat_nbr[ct_dir]);
 				r_max = r_nbr;
 				c_max = c_nbr;
+				max_side = ct_dir;
 			    }
 			}
 			else if (ct_dir == np_side) {
@@ -340,22 +347,22 @@ int do_cum_mfd(void)
 			}
 		    }
 		}
-		if (ABS(prop - 1.0) > 5E-6f) {
-		    G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
-			      prop);
-		}
 
 		/* adjust main drainage direction to A* path if possible */
-		if (fabs(wat_nbr[np_side]) >= max_acc) {
+		/*if (fabs(wat_nbr[np_side]) >= max_acc) {
 		    max_acc = fabs(wat_nbr[ct_dir]);
 		    r_max = dr;
 		    c_max = dc;
+		} */
+
+		if (ABS(prop - 1.0) > 5E-6f) {
+		    G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
+			      prop);
 		}
 	    }
 	    /* SFD-like accumulation */
 	    else {
-		seg_get(&watalt, (char *)&wa, dr, dc);
-		valued = wa.wat;
+		valued = wat_nbr[np_side];
 		if (value > 0) {
 		    if (valued > 0)
 			valued += value;
@@ -369,18 +376,47 @@ int do_cum_mfd(void)
 			valued = value - valued;
 		}
 		wa.wat = valued;
+		wa.ele = ele_nbr[np_side];
 		seg_put(&watalt, (char *)&wa, dr, dc);
 	    }
 
-	    if (!st_flag) {
-		/* update asp */
-		if (dr != r_max || dc != c_max) {
-		    if (asp_val > 0) {
-			asp_val = drain[r - r_max + 1][c - c_max + 1];
-			cseg_put(&asp, &asp_val, r, c);
-		    }
+	    /* update asp */
+	    if (dr != r_max || dc != c_max) {
+		if (asp_val < 0) {
+		    asp_val = -1 * drain[r - r_max + 1][c - c_max + 1];
+		}
+		else
+		    asp_val = drain[r - r_max + 1][c - c_max + 1];
+		    
+		cseg_put(&asp, &asp_val, r, c);
+	    }
+	    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) {
+		FLAG_SET(this_flag_value, SWALEFLAG);
+		is_swale = 1;
+	    }
+	    /* update asp for depression */
+	    if (is_swale && pit_flag) {
+		cseg_get(&asp, &asp_val_down, dr, dc);
+		if (asp_val > 0 && asp_val_down == 0) {
+		    asp_val = -asp_val;
+		    cseg_put(&asp, &asp_val, r, c);
 		}
 	    }
+	    /* continue stream */
+	    if (is_swale) {
+		flag_value = flag_nbr[max_side];
+		FLAG_SET(flag_value, SWALEFLAG);
+		bseg_put(&bitflags, &flag_value, r_max, c_max);
+	    }
+	    else {
+		if (er_flag && !is_swale && !FLAG_GET(this_flag_value, RUSLEBLOCKFLAG))
+		    slope_length(r, c, r_max, c_max);
+	    }
 	}
 	bseg_put(&bitflags, &this_flag_value, r, c);
     }

raster/r.watershed/seg/do_stream.c

@@ -1,224 +0,0 @@
-#include "Gwater.h"
-#include <unistd.h>
-#include <grass/gis.h>
-#include <grass/glocale.h>
-
-int do_stream(void)
-{
-    int r, c, dr, dc;
-    char is_swale;
-    DCELL value, *wat_nbr;
-    WAT_ALT wa;
-    POINT point;
-    int killer, threshold, count;
-
-    /* Breadth First Search */
-    int stream_cells, swale_cells;
-    int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
-    CELL ele, *ele_nbr, asp_val, asp_val_down;
-    double max_acc;
-    int edge;
-    char *flag_nbr, this_flag_value, flag_value;
-    int workedon;
-    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 };
-
-    G_message(_("SECTION 4: Extracting Streams."));
-
-    flag_nbr = (char *)G_malloc(sides * sizeof(char));
-    wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
-    ele_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
-
-    workedon = 0;
-
-    count = 0;
-    if (bas_thres <= 0)
-	threshold = 60;
-    else
-	threshold = bas_thres;
-
-    for (killer = 0; killer < do_points; killer++) {
-	G_percent(count++, do_points, 1);
-	seg_get(&astar_pts, (char *)&point, 0, killer);
-	r = point.r;
-	c = point.c;
-	/* here is original A* direction */
-	asp_val = point.asp;
-	/* do_stream() is only executed for MFD and if stream segments are needed
-	 * therefore it is ok to fetch asp from point.asp
-	 * if for some reason do_cum_mfd updates asp even if followed by
-	 * do_stream, the following line needs to be uncommented and
-	 * the above line can be commented out */
-	/* cseg_get(&asp, &asp_val, r, c); */
-	if (asp_val) {
-	    dr = r + asp_r[ABS(asp_val)];
-	    dc = c + asp_c[ABS(asp_val)];
-	}
-	else
-	    dr = dc = -1;
-	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) {
-	    bseg_get(&bitflags, &this_flag_value, r, c);
-	    /* do not continue streams along edges, this causes artifacts */
-	    if (FLAG_GET(this_flag_value, EDGEFLAG)) {
-		is_swale = FLAG_GET(this_flag_value, SWALEFLAG);
-		if (is_swale && asp_val > 0) {
-
-		    /* find first neighbour that is NULL or outside region */
-		    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
-			r_nbr = r + nextdr[ct_dir];
-			c_nbr = c + nextdc[ct_dir];
-
-			/* check if neighbour is within region */
-			if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
-			    c_nbr < ncols) {
-
-			    bseg_get(&bitflags, &flag_value, r_nbr, c_nbr);
-			    if (FLAG_GET(flag_value, NULLFLAG)) {
-				asp_val = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
-				break;
-			    }
-			}
-			else {
-			    asp_val = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
-			    break;
-			}
-		    }
-		    cseg_put(&asp, &asp_val, r, c);
-		}
-		continue;
-	    }
-
-	    seg_get(&watalt, (char *)&wa, r, c);
-	    value = wa.wat;
-	    if (point.guessed && value > 0)
-		value = -value;
-
-	    r_max = dr;
-	    c_max = dc;
-
-	    np_side = -1;
-	    stream_cells = 0;
-	    swale_cells = 0;
-	    max_acc = -1;
-	    ele = wa.ele;
-	    edge = 0;
-	    /* visit all neighbours */
-	    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];
-
-		wat_nbr[ct_dir] = 0;
-		ele_nbr[ct_dir] = 0;
-		FLAG_SET(flag_nbr[ct_dir], WORKEDFLAG);
-
-		if (dr == r_nbr && dc == c_nbr)
-		    np_side = ct_dir;
-
-		/* check if neighbour is within region */
-		if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
-		    c_nbr < ncols) {
-
-		    /* check for swale or stream cells */
-		    bseg_get(&bitflags, &flag_nbr[ct_dir], r_nbr, c_nbr);
-		    is_swale = FLAG_GET(flag_nbr[ct_dir], SWALEFLAG);
-		    if (is_swale)
-			swale_cells++;
-		    seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
-		    wat_nbr[ct_dir] = wa.wat;
-
-		    if (fabs(wat_nbr[ct_dir]) >= threshold)
-			stream_cells++;
-
-		    if (!FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
-
-			ele_nbr[ct_dir] = wa.ele;
-			
-			edge = FLAG_GET(flag_nbr[ct_dir], NULLFLAG);
-			if (!edge && ele_nbr[ct_dir] <= ele) {
-
-			    /* set main drainage direction */
-			    if (fabs(wat_nbr[ct_dir]) >= max_acc) {
-				max_acc = fabs(wat_nbr[ct_dir]);
-				r_max = r_nbr;
-				c_max = c_nbr;
-			    }
-			}
-		    }
-		    else if (ct_dir == np_side && !edge) {
-			/* check for consistency with main drainage direction */
-			workedon++;
-		    }
-		}
-		else
-		    edge = 1;
-		if (edge)
-		    break;
-	    }
-	    /* do not continue streams along edges, this causes artifacts */
-	    if (edge) {
-		is_swale = FLAG_GET(this_flag_value, SWALEFLAG);
-		if (is_swale && asp_val > 0) {
-		    asp_val = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
-		    cseg_put(&asp, &asp_val, r, c);
-		}
-		continue;
-	    }
-
-	    /* adjust main drainage direction to A* path if possible */
-	    if (fabs(wat_nbr[np_side]) >= max_acc) {
-		max_acc = fabs(wat_nbr[ct_dir]);
-		r_max = dr;
-		c_max = dc;
-	    }
-	    /* update asp */
-	    if (dr != r_max || dc != c_max) {
-		if (asp_val > 0) {
-		    asp_val = drain[r - r_max + 1][c - c_max + 1];
-		    cseg_put(&asp, &asp_val, r, c);
-		}
-	    }
-	    is_swale = FLAG_GET(this_flag_value, SWALEFLAG);
-	    /* start new stream */
-	    value = fabs(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) {
-		FLAG_SET(this_flag_value, SWALEFLAG);
-		is_swale = 1;
-	    }
-	    /* update asp for depression */
-	    if (is_swale && pit_flag) {
-		cseg_get(&asp, &asp_val_down, dr, dc);
-		if (asp_val > 0 && asp_val_down == 0) {
-		    asp_val = -asp_val;
-		    cseg_put(&asp, &asp_val, r, c);
-		}
-	    }
-	    /* continue stream */
-	    if (is_swale) {
-		bseg_get(&bitflags, &flag_value, r_max, c_max);
-		FLAG_SET(flag_value, SWALEFLAG);
-		bseg_put(&bitflags, &flag_value, r_max, c_max);
-	    }
-	    else {
-		if (er_flag && !is_swale && !FLAG_GET(this_flag_value, RUSLEBLOCKFLAG))
-		    slope_length(r, c, r_max, c_max);
-	    }
-	    FLAG_SET(this_flag_value, WORKEDFLAG);
-	    bseg_put(&bitflags, &this_flag_value, r, c);
-	}
-    }
-    G_percent(count, do_points, 1);	/* finish it */
-    if (workedon)
-	G_warning(_("MFD: A * path already processed when extracting streams: %d of %d cells"),
-		  workedon, do_points);
-
-    seg_close(&astar_pts);
-    
-    G_free(wat_nbr);
-    G_free(ele_nbr);
-    G_free(flag_nbr);
-
-    return 0;
-}

raster/r.watershed/seg/init_vars.c

@@ -399,7 +399,7 @@ 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)));
+    seg_cols = (int) (pow(2, (int)(log(num_open_segs / 8.0) / log(2) + 0.1)) + 0.1);
     num_open_array_segs = num_open_segs / seg_cols;
     /* n cols in segment */
     seg_cols *= seg_rows * seg_rows;
@@ -420,7 +420,7 @@ int init_vars(int argc, char *argv[])
     /* one-based d-ary search_heap with astar_pts */
     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)));
+    seg_cols = (int) (pow(2, (int)(log(num_open_segs / 8.0) / log(2) + 0.1)) + 0.1);
     num_open_array_segs = num_open_segs / seg_cols;
     /* n cols in segment */
     seg_cols *= seg_rows * seg_rows;

raster/r.watershed/seg/main.c

@@ -82,9 +82,6 @@ int main(int argc, char *argv[])
     else {
 	do_cum();
     }
-    if (st_flag && mfd) {
-	do_stream();
-    }
     if (sg_flag || ls_flag) {
 	sg_factor();
     }

raster/r.watershed/seg/no_stream.c

@@ -12,6 +12,7 @@ no_stream(int row, int col, CELL basin_num, double stream_length,
     WAT_ALT wa;
 
     while (1) {
+	cseg_put(&bas, &basin_num, row, col);
 	max_drain = -1;
 	for (r = row - 1, rr = 0; r <= row + 1; r++, rr++) {
 	    for (c = col - 1, cc = 0; c <= col + 1; c++, cc++) {
@@ -71,23 +72,25 @@ no_stream(int row, int col, CELL basin_num, double stream_length,
 			cseg_get(&asp, &aspect, r, c);
 			if (aspect == drain[rr][cc]) {
 			    thisdir = updrain[rr][cc];
-			    if (haf_basin_side(updir,
+			    switch (haf_basin_side(updir,
 					       (SHORT) downdir,
-					       thisdir) == RITE) {
+					       thisdir)) {
+			    case RITE:
 				overland_cells(r, c, basin_num, basin_num,
 					       &new_ele);
 				riteflag++;
-			    }
-			    else {
+				break;
+			    case LEFT:
 				overland_cells(r, c, basin_num, basin_num - 1,
 					       &new_ele);
 				leftflag++;
+				break;
 			    }
 			}
 		    }
 		}
 	    }
-	    if (leftflag >= riteflag) {
+	    if (leftflag > riteflag) {
 		value = basin_num - 1;
 		cseg_put(&haf, &value, row, col);
 	    }
@@ -104,6 +107,7 @@ no_stream(int row, int col, CELL basin_num, double stream_length,
 		    slope = MIN_SLOPE;
 		fprintf(fp, " %f %f\n", slope, stream_length);
 	    }
+	    cseg_put(&haf, &basin_num, row, col);
 	    return 0;
 	}
     }

raster/r.watershed/seg/over_cells.c

@@ -40,7 +40,6 @@ overland_cells(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
 {
     int r, rr, c, cc, next_r, next_c;
     CELL value;
-    OC_STACK cur;
     int top = 0;
 
     /* put root on stack */
@@ -53,9 +52,8 @@ overland_cells(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
 
     while (top) {
 	/* assign basin numbers */
-	cur = ocs[--top];
-	next_r = cur.row;
-	next_c = cur.col;
+	next_r = ocs[--top].row;
+	next_c = ocs[top].col;
 
 	/* add all neighbours pouring into current cell to stack */
 	for (r = next_r - 1, rr = 0; r <= next_r + 1; r++, rr++) {