fix for interactive mode

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

raster/r.watershed/seg/Gwater.h

@@ -5,7 +5,8 @@
 /* program to map out drainage basin structure  */
 /* this one uses the A * search algorithm       */
 /* written by Chuck Ehlschlaeger                */
-/* last modified 03/26/91                       */
+/* updated by Markus Metz                       */
+/* last modified 22/10/09                       */
 #include <math.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
@@ -62,12 +63,20 @@ WAT_ALT {
    DCELL wat;
 };
 
+#define OC_STACK struct overland_cells_stack
+OC_STACK {
+    int row, col;
+};
+
 extern struct Cell_head window;
 
 extern int mfd, c_fac, abs_acc, ele_scale;
 extern SSEG search_heap;
 extern int heap_size;
 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 double half_res, diag, max_length, dep_slope;
 extern int bas_thres, tot_parts;

raster/r.watershed/seg/close_maps2.c

@@ -20,6 +20,7 @@ int close_array_seg(void)
 	else
 	    theseg = &haf;
 	max = -9;
+	/*
 	for (r = 0; r < nrows; r++) {
 	    for (c = 0; c < ncols; c++) {
 		cseg_get(theseg, &value, r, c);
@@ -27,6 +28,8 @@ int close_array_seg(void)
 		    max = value;
 	    }
 	}
+	*/
+	max = n_basins;
 	G_debug(1, "%d basins created", max);
 	Rast_init_colors(&colors);
 	Rast_make_random_colors(&colors, 1, max);

raster/r.watershed/seg/def_basin.c

@@ -11,9 +11,6 @@ def_basin(int row, int col, CELL basin_num, double stream_length,
 
     for (;;) {
 	cseg_put(&bas, &basin_num, row, col);
-	cseg_get(&asp, &asp_value, row, col);
-	if (asp_value < 0)
-	    asp_value = -asp_value;
 	ct = 0;
 	for (r = row - 1, rr = 0; rr < 3; r++, rr++) {
 	    for (c = col - 1, cc = 0; cc < 3; c++, cc++) {
@@ -75,28 +72,33 @@ def_basin(int row, int col, CELL basin_num, double stream_length,
 	else {
 	    cseg_put(&haf, &basin_num, row, col);
 	}
-	if (sides == 8) {
-	    if (new_r[1] != row && new_c[1] != col)
-		stream_length += diag;
-	    else if (new_r[1] != row)
-		stream_length += window.ns_res;
-	    else
-		stream_length += window.ew_res;
-	}
-	else {			/* sides == 4 */
-
-	    if (asp_value == 2 || asp_value == 6) {
-		if (new_r[1] != row)
+	if (arm_flag) {
+	    if (sides == 8) {
+		if (new_r[1] != row && new_c[1] != col)
+		    stream_length += diag;
+		else if (new_r[1] != row)
 		    stream_length += window.ns_res;
 		else
-		    stream_length += diag;
+		    stream_length += window.ew_res;
 	    }
-	    else {		/* asp_value == 4, 8 */
+	    else {			/* sides == 4 */
 
-		if (new_c[1] != col)
-		    stream_length += window.ew_res;
-		else
-		    stream_length += diag;
+		cseg_get(&asp, &asp_value, row, col);
+		if (asp_value < 0)
+		    asp_value = -asp_value;
+		if (asp_value == 2 || asp_value == 6) {
+		    if (new_r[1] != row)
+			stream_length += window.ns_res;
+		    else
+			stream_length += diag;
+		}
+		else {		/* asp_value == 4, 8 */
+
+		    if (new_c[1] != col)
+			stream_length += window.ew_res;
+		    else
+			stream_length += diag;
+		}
 	    }
 	}
 	row = new_r[1];

raster/r.watershed/seg/do_cum.c

@@ -29,6 +29,7 @@ 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)];
@@ -171,6 +172,7 @@ 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)];

raster/r.watershed/seg/do_stream.c

@@ -42,7 +42,14 @@ int do_stream(void)
 	seg_get(&astar_pts, (char *)&point, 0, killer);
 	r = point.r;
 	c = point.c;
-	cseg_get(&asp, &asp_val, r, 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)];
@@ -174,7 +181,8 @@ int do_stream(void)
 	    is_swale = FLAG_GET(this_flag_value, SWALEFLAG);
 	    /* start new stream */
 	    value = fabs(value) + 0.5;
-	    if (!is_swale && (int)value >= threshold && stream_cells < 4 &&
+	    /* 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;

raster/r.watershed/seg/find_pour.c

@@ -7,36 +7,34 @@ int find_pourpts(void)
     CELL old_elev, basin_num, value;
     char is_swale;
 
+    ocs_alloced = 2 * bas_thres;
+    ocs = (OC_STACK *)G_malloc(ocs_alloced * sizeof(OC_STACK));
+
     basin_num = 0;
+    stream_length = old_elev = 0;
     for (row = 0; row < nrows; row++) {
 	G_percent(row, nrows, 1);
 	northing = window.north - (row + .5) * window.ns_res;
 	for (col = 0; col < ncols; col++) {
-	    /* cseg_get (&wat, &value, row, col);
-	       if (value < 0)
-	       {
-	       value = -value;
-	       } */
 	    cseg_get(&asp, &value, row, col);
 	    bseg_get(&bitflags, &is_swale, row, col);
 	    is_swale = FLAG_GET(is_swale, SWALEFLAG);
-	    /* bseg_get(&swale, &is_swale, row, col); */
 	    if (value < 0 && is_swale > 0) {
 		basin_num += 2;
-		cseg_get(&alt, &old_elev, row, col);
 		if (arm_flag) {
 		    easting = window.west + (col + .5) * window.ew_res;
 		    fprintf(fp, "%5d drains into %5d at %3d %3d %.3f %.3f",
 			    (int)basin_num, 0, row, col, easting, northing);
-		}
-		if (col == 0 || col == ncols - 1) {
-		    stream_length = .5 * window.ew_res;
-		}
-		else if (row == 0 || row == nrows - 1) {
-		    stream_length = .5 * window.ns_res;
-		}
-		else {
-		    stream_length = 0.0;
+		    if (col == 0 || col == ncols - 1) {
+			stream_length = .5 * window.ew_res;
+		    }
+		    else if (row == 0 || row == nrows - 1) {
+			stream_length = .5 * window.ns_res;
+		    }
+		    else {
+			stream_length = 0.0;
+		    }
+		    cseg_get(&alt, &old_elev, row, col);
 		}
 		basin_num =
 		    def_basin(row, col, basin_num, stream_length, old_elev);
@@ -44,6 +42,8 @@ int find_pourpts(void)
 	}
     }
     G_percent(nrows, nrows, 1);	/* finish it */
+    n_basins = basin_num;
+    G_free(ocs);
 
     return 0;
 }

raster/r.watershed/seg/init_vars.c

@@ -153,18 +153,22 @@ int init_vars(int argc, char *argv[])
 	diag *= 0.5;
 
     /* segment parameters: one size fits all. Fine tune? */
-    /* Segment rows and cols: 128 */
-    /* 1 segment open for all data structures: 0.122 MB */
+    /* Segment rows and cols: 64 */
+    /* 1 segment open for all data structures used in A* Search: 0.18 MB */
+    /* is really 0.22 MB but search heap holds max 5% of all points
+     * i.e. we will probably never need all open segments for search heap
+     */
 
     seg_rows = SROW;
     seg_cols = SCOL;
 
     if (segs_mb < 3.0) {
 	segs_mb = 300;
-	G_warning(_("Maximum memory to be used was smaller than 3 MB, set to default = 300 MB."));
+	G_warning(_("Maximum memory to be used was smaller than 3 MB,"
+	            " set to default = 300 MB."));
     }
 
-    num_open_segs = segs_mb / 0.122;
+    num_open_segs = segs_mb / 0.18;
 
     G_debug(1, "segs MB: %.0f", segs_mb);
     G_debug(1, "region rows: %d", nrows);
@@ -187,6 +191,12 @@ int init_vars(int argc, char *argv[])
 	num_open_segs = num_cseg_total;
     G_debug(1, "  open segments after adjusting:\t%d", num_open_segs);
 
+    if (num_cseg_total * 0.18 < 1024.0)
+	G_verbose_message(_("Will need at most %.2f MB of disk space"), num_cseg_total * 0.18);
+    else
+	G_verbose_message(_("Will need at most %.2f GB (%.2f MB) of disk space"),
+	           (num_cseg_total * 0.18) / 1024.0, num_cseg_total * 0.18);
+
     if (er_flag) {
 	cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
 	cseg_read_cell(&r_h, ele_name, "");
@@ -348,48 +358,49 @@ int init_vars(int argc, char *argv[])
 	G_free(buf);
     }
 
-    if (ob_flag) {
-	fd = Rast_open_old(ob_name, "");
-	if (fd < 0) {
-	    G_fatal_error(_("unable to open blocking map layer"));
-	}
-	buf = Rast_allocate_c_buf();
-	for (r = 0; r < nrows; r++) {
-	    G_percent(r, nrows, 1);
-	    Rast_get_c_row(fd, buf, r);
-	    for (c = 0; c < ncols; c++) {
-		block_value = buf[c];
-		if (!Rast_is_c_null_value(&block_value) && block_value) {
-		    bseg_get(&bitflags, &flag_value, r, c);
-		    FLAG_SET(flag_value, RUSLEBLOCKFLAG);
-		    bseg_put(&bitflags, &flag_value, r, c);
+    /* do RUSLE */
+    if (er_flag) {
+	if (ob_flag) {
+	    fd = Rast_open_old(ob_name, "");
+	    if (fd < 0) {
+		G_fatal_error(_("unable to open blocking map layer"));
+	    }
+	    buf = Rast_allocate_c_buf();
+	    for (r = 0; r < nrows; r++) {
+		G_percent(r, nrows, 1);
+		Rast_get_c_row(fd, buf, r);
+		for (c = 0; c < ncols; c++) {
+		    block_value = buf[c];
+		    if (!Rast_is_c_null_value(&block_value) && block_value) {
+			bseg_get(&bitflags, &flag_value, r, c);
+			FLAG_SET(flag_value, RUSLEBLOCKFLAG);
+			bseg_put(&bitflags, &flag_value, r, c);
+		    }
 		}
 	    }
+	    G_percent(nrows, nrows, 1);    /* finish it */
+	    Rast_close(fd);
+	    G_free(buf);
 	}
-	G_percent(nrows, nrows, 1);    /* finish it */
-	Rast_close(fd);
-	G_free(buf);
-    }
 
-    if (ril_flag) {
-	dseg_open(&ril, 1, seg_rows * seg_cols, num_open_segs);
-	dseg_read_cell(&ril, ril_name, "");
-    }
-    
-    /* dseg_open(&slp, SROW, SCOL, num_open_segs); */
+	if (ril_flag) {
+	    dseg_open(&ril, 1, seg_rows * seg_cols, num_open_segs);
+	    dseg_read_cell(&ril, ril_name, "");
+	}
+	
+	/* dseg_open(&slp, SROW, SCOL, num_open_segs); */
 
-    /* RUSLE: LS and/or S factor */
-    if (er_flag) {
-	dseg_open(&s_l, seg_rows, seg_cols, num_open_segs);
+	    dseg_open(&s_l, seg_rows, seg_cols, num_open_segs);
+	if (sg_flag)
+	    dseg_open(&s_g, 1, seg_rows * seg_cols, num_open_segs);
+	if (ls_flag)
+	    dseg_open(&l_s, 1, seg_rows * seg_cols, num_open_segs);
     }
-    if (sg_flag)
-	dseg_open(&s_g, 1, seg_rows * seg_cols, num_open_segs);
-    if (ls_flag)
-	dseg_open(&l_s, 1, seg_rows * seg_cols, num_open_segs);
 
     G_debug(1, "open segments for A* points");
-    /* next smaller power of 2 */
-    seg_cols = (int) pow(2, (int)(log(num_open_segs / 8) / log(2)));
+    /* rounded down power of 2 */
+    seg_cols = (int) pow(2, (int)(log(num_open_segs / 8.0) / log(2)));
+    num_open_array_segs = num_open_segs / seg_cols;
     /* n cols in segment */
     seg_cols *= seg_rows * seg_rows;
     /* n segments in row */
@@ -397,20 +408,20 @@ int init_vars(int argc, char *argv[])
     if (do_points % seg_cols > 0)
 	num_cseg_total++;
     /* no need to have more segments open than exist */
-    if ((num_open_segs * seg_rows * seg_rows) / (seg_cols * 2) > num_cseg_total)
+    if (num_open_array_segs > num_cseg_total)
 	num_open_array_segs = num_cseg_total;
-    else
-	num_open_array_segs = (num_open_segs * seg_rows * seg_rows) / (seg_cols * 2);
-    
-    G_debug(2, "A* points open segments %d, target 4", num_open_array_segs);
+
+    if (num_open_array_segs > 4)
+	num_open_array_segs = 4;
     
     seg_open(&astar_pts, 1, do_points, 1, seg_cols, num_open_array_segs,
 	     sizeof(POINT));
 
     /* one-based d-ary search_heap with astar_pts */
     G_debug(1, "open segments for A* search heap");
-    /* next smaller power of 2 */
-    seg_cols = (int) pow(2, (int)(log(num_open_segs / 16) / log(2)));
+    /* rounded down power of 2 */
+    seg_cols = (int) pow(2, (int)(log(num_open_segs / 8.0) / log(2)));
+    num_open_array_segs = num_open_segs / seg_cols;
     /* n cols in segment */
     seg_cols *= seg_rows * seg_rows;
     /* n segments in row */
@@ -418,12 +429,12 @@ int init_vars(int argc, char *argv[])
     if ((do_points + 1) % seg_cols > 0)
 	num_cseg_total++;
     /* no need to have more segments open than exist */
-    if ((num_open_segs * seg_rows * seg_rows) / (seg_cols * 2) > num_cseg_total)
+    if (num_open_array_segs > num_cseg_total)
 	num_open_array_segs = num_cseg_total;
-    else
-	num_open_array_segs = (num_open_segs * seg_rows * seg_rows) / (seg_cols * 2);
 
-    G_debug(2, "A* search heap open segments %d, target 8", num_open_array_segs);
+    G_debug(1, "A* search heap open segments %d, target 8", num_open_array_segs);
+    /* the search heap will not hold more than 5% of all points at any given time ? */
+    /* chances are good that the heap will fit into one large segment */
     seg_open(&search_heap, 1, do_points + 1, 1, seg_cols,
 	     num_open_array_segs, sizeof(HEAP_PNT));
 

raster/r.watershed/seg/main.c

@@ -29,6 +29,9 @@ int mfd, c_fac, abs_acc, ele_scale;
 SSEG search_heap;
 int heap_size;
 int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
+CELL n_basins;
+OC_STACK *ocs;
+int ocs_alloced;
 SHORT nrows, ncols;
 double half_res, diag, max_length, dep_slope;
 int bas_thres, tot_parts;
@@ -94,7 +97,7 @@ int main(int argc, char *argv[])
 	if (arm_flag) {
 	    fp = fopen(arm_name, "w");
 	}
-	num_open_segs = segs_mb / 0.122;
+	num_open_segs = segs_mb / 0.2;
 	if (num_open_segs > (ncols / SCOL + 1) * (nrows / SROW + 1)) {
 	    num_open_segs = (ncols / SCOL + 1) * (nrows / SROW + 1);
 	}

raster/r.watershed/seg/no_stream.c

@@ -23,7 +23,7 @@ no_stream(int row, int col, CELL basin_num, double stream_length,
 			dvalue = wa.wat;
 			if (dvalue < 0)
 			    dvalue = -dvalue;
-			if ((dvalue - max_drain) > 5E-8f) {
+			if (dvalue > max_drain) {
 			    uprow = r;
 			    upcol = c;
 			    max_drain = dvalue;
@@ -37,29 +37,31 @@ no_stream(int row, int col, CELL basin_num, double stream_length,
 	    cseg_get(&asp, &downdir, row, col);
 	    if (downdir < 0)
 		downdir = -downdir;
-	    if (sides == 8) {
-		if (uprow != row && upcol != col)
-		    stream_length += diag;
-		else if (uprow != row)
-		    stream_length += window.ns_res;
-		else
-		    stream_length += window.ew_res;
-	    }
-	    else {		/* sides == 4 */
-
-		cseg_get(&asp, &asp_value, uprow, upcol);
-		if (downdir == 2 || downdir == 6) {
-		    if (asp_value == 2 || asp_value == 6)
+	    if (arm_flag) {
+		if (sides == 8) {
+		    if (uprow != row && upcol != col)
+			stream_length += diag;
+		    else if (uprow != row)
 			stream_length += window.ns_res;
 		    else
-			stream_length += diag;
+			stream_length += window.ew_res;
 		}
-		else {		/* downdir == 4,8 */
+		else {		/* sides == 4 */
 
-		    if (asp_value == 4 || asp_value == 8)
-			stream_length += window.ew_res;
-		    else
-			stream_length += diag;
+		    cseg_get(&asp, &asp_value, uprow, upcol);
+		    if (downdir == 2 || downdir == 6) {
+			if (asp_value == 2 || asp_value == 6)
+			    stream_length += window.ns_res;
+			else
+			    stream_length += diag;
+		    }
+		    else {		/* downdir == 4,8 */
+
+			if (asp_value == 4 || asp_value == 8)
+			    stream_length += window.ew_res;
+			else
+			    stream_length += diag;
+		    }
 		}
 	    }
 	    riteflag = leftflag = 0;

raster/r.watershed/seg/over_cells.c

@@ -2,7 +2,7 @@
 #define BIGNEG	-9999999
 
 int
-overland_cells(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
+overland_cells_recursive(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
 {
     int r, rr, c, cc;
     CELL new_ele, new_max_ele, value;
@@ -15,25 +15,80 @@ overland_cells(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
 	    if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
 		cseg_get(&asp, &value, r, c);
 		if (value == drain[rr][cc]) {
-		    if (r != row && c != col) {
-			overland_cells(r, c, basin_num, haf_num, &new_ele);
-		    }
-		    else if (r != row) {
-			overland_cells(r, c, basin_num, haf_num, &new_ele);
-		    }
-		    else {
-			overland_cells(r, c, basin_num, haf_num, &new_ele);
-		    }
+		    overland_cells(r, c, basin_num, haf_num, &new_ele);
 		}
 	    }
 	}
     }
-    if (new_max_ele == BIGNEG) {
-	cseg_get(&alt, hih_ele, row, col);
+    /*
+    if (arm_flag) {
+	if (new_max_ele == BIGNEG) {
+	    cseg_get(&alt, hih_ele, row, col);
+	}
+	else {
+	    *hih_ele = new_max_ele;
+	}
     }
-    else {
-	*hih_ele = new_max_ele;
+    */
+
+    return 0;
+}
+
+/* non-recursive version */
+int
+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 */
+    ocs[top].row = row;
+    ocs[top].col = col;
+    cseg_put(&bas, &basin_num, row, col);
+    cseg_put(&haf, &haf_num, row, col);
+
+    top++;
+
+    while (top) {
+	/* assign basin numbers */
+	cur = ocs[--top];
+	next_r = cur.row;
+	next_c = cur.col;
+
+	/* add all neighbours pouring into current cell to stack */
+	for (r = next_r - 1, rr = 0; r <= next_r + 1; r++, rr++) {
+	    for (c = next_c - 1, cc = 0; c <= next_c + 1; c++, cc++) {
+		if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+		    cseg_get(&asp, &value, r, c);
+		    if (value == drain[rr][cc]) {
+			if (top >= ocs_alloced) {
+			    ocs_alloced += bas_thres;
+			    ocs = (OC_STACK *)G_realloc(ocs, ocs_alloced * sizeof(OC_STACK));
+			}
+			ocs[top].row = r;
+			ocs[top].col = c;
+			cseg_put(&bas, &basin_num, r, c);
+			cseg_put(&haf, &haf_num, r, c);
+			top++;
+		    }
+		}
+	    }
+	}
+	
+    }
+
+    /*
+    if (arm_flag) {
+	if (new_max_ele == BIGNEG) {
+	    cseg_get(&alt, hih_ele, row, col);
+	}
+	else {
+	    *hih_ele = new_max_ele;
+	}
     }
+    */
 
     return 0;
 }

raster/r.watershed/seg/split_str.c

@@ -10,6 +10,8 @@ split_stream(int row, int col, int new_r[], int new_c[], int ct,
     SHORT thisdir, leftflag, riteflag;
     int r, c, rr, cc;
 
+    new_elev = 0;
+
     for (ctr = 1; ctr <= ct; ctr++)
 	splitdir[ctr] = drain[row - new_r[ctr] + 1][col - new_c[ctr] + 1];
     updir = splitdir[1];
@@ -58,16 +60,17 @@ split_stream(int row, int col, int new_r[], int new_c[], int ct,
 	cseg_put(&haf, &basin_num, row, col);
     }
     old_basin = basin_num;
-    cseg_get(&alt, &new_elev, row, col);
-    if ((slope = (new_elev - old_elev) / stream_length) < MIN_SLOPE)
-	slope = MIN_SLOPE;
-    if (arm_flag)
+    if (arm_flag) {
+	cseg_get(&alt, &new_elev, row, col);
+	if ((slope = (new_elev - old_elev) / stream_length) < MIN_SLOPE)
+	    slope = MIN_SLOPE;
 	fprintf(fp, " %f %f\n", slope, stream_length);
+    }
     for (r = 1; r <= ct; r++) {
 	basin_num += 2;
-	easting = window.west + (new_c[r] + .5) * window.ew_res;
-	northing = window.north - (new_r[r] + .5) * window.ns_res;
 	if (arm_flag) {
+	    easting = window.west + (new_c[r] + .5) * window.ew_res;
+	    northing = window.north - (new_r[r] + .5) * window.ns_res;
 	    fprintf(fp, "%5d drains into %5d at %3d %3d %.3f %.3f",
 		    (int)basin_num, old_basin, new_r[r], new_c[r], easting,
 		    northing);