r.watershed: add TCI, remove edge artefacts for SFD (ram)

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

raster/r.watershed/ram/Gwater.h

@@ -60,7 +60,7 @@ extern RAMSEG r_h_seg, dep_seg;
 extern RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
 extern int *astar_pts;
 extern CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
-extern DCELL *wat;
+extern DCELL *wat, *tci;
 extern int ril_fd;
 extern double *s_l, *s_g, *l_s;
 extern CELL one, zero;
@@ -76,9 +76,10 @@ extern char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
 extern const char *this_mapset;
 extern char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX], thr_name[8];
 extern char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX], sg_name[GNAME_MAX];
-extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX], dis_name[GNAME_MAX];
+extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX], tci_name[GNAME_MAX];
+extern char arm_name[GNAME_MAX], dis_name[GNAME_MAX];
 extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, flat_flag;
-extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
+extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag, tci_flag;
 extern char bas_flag, seg_flag, haf_flag, er_flag;
 extern char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
 extern FILE *fp;

raster/r.watershed/ram/close_maps.c

@@ -8,7 +8,7 @@
 int close_maps(void)
 {
     struct Colors colors;
-    int value, r, c, fd;
+    int r, c, fd;
     CELL *buf = NULL;
     DCELL *dbuf = NULL;
     struct FPRange accRange;
@@ -129,6 +129,69 @@ int close_maps(void)
 	Rast_write_colors(wat_name, this_mapset, &colors);
     }
 
+    /* TCI */
+    if (tci_flag) {
+	DCELL watvalue;
+
+	sum = sum_sqr = stddev = 0.0;
+	fd = Rast_open_new(tci_name, DCELL_TYPE);
+	for (r = 0; r < nrows; r++) {
+	    Rast_set_d_null_value(dbuf, ncols);	/* reset row to all NULL */
+	    for (c = 0; c < ncols; c++) {
+		dvalue = tci[SEG_INDEX(wat_seg, r, c)];
+		watvalue = wat[SEG_INDEX(wat_seg, r, c)];
+		if (!Rast_is_d_null_value(&watvalue)) {
+		    dbuf[c] = dvalue;
+		    sum += dvalue;
+		    sum_sqr += dvalue * dvalue;
+		}
+	    }
+	    Rast_put_row(fd, dbuf, DCELL_TYPE);
+	}
+	Rast_close(fd);
+
+	stddev =
+	    sqrt((sum_sqr - (sum + sum / do_points)) / (do_points - 1));
+	G_debug(1, "stddev: %f", stddev);
+
+	/* set nice color rules: yellow, green, cyan, blue, black */
+
+	lstddev = log(stddev);
+
+	Rast_read_fp_range(tci_name, this_mapset, &accRange);
+	min = max = 0;
+	Rast_get_fp_range_min_max(&accRange, &min, &max);
+
+	Rast_init_colors(&colors);
+
+	clr_min = min - 1;
+	clr_max = min + (max - min) * 0.3;
+	Rast_add_d_color_rule(&clr_min, 255, 255, 0, &clr_max, 255,
+				  255, 0, &colors);
+	clr_min = clr_max;
+	clr_max = min + (max - min) * 0.5;
+	Rast_add_d_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
+				  255, 0, &colors);
+	clr_min = clr_max;
+	clr_max = min + (max - min) * 0.6;
+	Rast_add_d_color_rule(&clr_min, 0, 255, 0, &clr_max, 0,
+				  255, 255, &colors);
+	clr_min = clr_max;
+	clr_max = min + (max - min) * 0.7;
+	Rast_add_d_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
+				  0, 255, &colors);
+	clr_min = clr_max;
+	clr_max = max + 1.;
+	Rast_add_d_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
+				  0, &colors);
+
+	clr_min = clr_max;
+	clr_max = max + 1;
+	Rast_add_d_color_rule(&clr_min, 0, 0, 0, &clr_max, 0, 0,
+				  0, &colors);
+	Rast_write_colors(tci_name, this_mapset, &colors);
+    }
+
     /* TODO: elevation == NULL -> drainage direction == NULL (wat == 0 where ele == NULL) */
     /* keep drainage direction == 0 for real depressions */
     if (asp_flag) {

raster/r.watershed/ram/do_astar.c

@@ -163,7 +163,7 @@ int do_astar(void)
 			asp[index_up] = drain[upr - r + 1][upc - c + 1];
 
 			if (wat[index_doer] > 0)
-			    wat[index_doer] = -wat[index_doer];
+			    wat[index_doer] = -1.0 * wat[index_doer];
 		    }
 		    /* neighbour is inside real depression, not yet worked */
 		    else if (asp[index_up] == 0) {
@@ -246,7 +246,7 @@ int drop_pt(void)
 
     /* sift down: move hole back towards bottom of heap */
 
-    while ((child = GET_CHILD(parent)) <= heap_size) {
+    while (GET_CHILD(child, parent) <= heap_size) {
 	/* 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]];
@@ -298,7 +298,7 @@ int sift_up(int start, CELL ele)
     child_idx = astar_pts[child];
 
     while (child > 1) {
-	parent = GET_PARENT(child);
+	GET_PARENT(parent, child);
 
 	elep = alt[astar_pts[parent]];
 	/* child smaller */

raster/r.watershed/ram/do_astar.h

@@ -1,7 +1,7 @@
 #ifndef __DO_ASTAR_H__
 #define __DO_ASTAR_H__
 
-#define GET_PARENT(c) (int) (((c) - 2) / 3 + 1)
-#define GET_CHILD(p) (int) (((p) * 3) - 1)
+#define GET_PARENT(p, c) ((p) = (int) (((c) - 2) / 3 + 1))
+#define GET_CHILD(c, p) ((c) = (int) (((p) * 3) - 1))
 
 #endif /* __DO_ASTAR_H__ */

raster/r.watershed/ram/do_cum.c

@@ -3,20 +3,95 @@
 #include <grass/raster.h>
 #include <grass/glocale.h>
 
+int get_dist(double *dist_to_nbr, double *contour)
+{
+    int ct_dir, r_nbr, c_nbr;
+    double dx, dy;
+    double ns_res, ew_res;
+
+    if (G_projection() == PROJECTION_LL) {
+	double ew_dist1, ew_dist2, ew_dist3;
+	double ns_dist1, ns_dist2, ns_dist3;
+
+	G_begin_distance_calculations();
+
+	/* EW Dist at North edge */
+	ew_dist1 = G_distance(window.east, window.north,
+	                      window.west, window.north);
+	/* EW Dist at Center */
+	ew_dist2 = G_distance(window.east, (window.north + window.south) / 2.,
+	                      window.west, (window.north + window.south) / 2.);
+	/* EW Dist at South Edge */
+	ew_dist3 = G_distance(window.east, window.south,
+	                      window.west, window.south);
+	/* NS Dist at East edge */
+	ns_dist1 = G_distance(window.east, window.north,
+	                      window.east, window.south);
+	/* NS Dist at Center */
+	ns_dist2 = G_distance((window.west + window.east) / 2., window.north,
+	                      (window.west + window.east) / 2., window.south);
+	/* NS Dist at West edge */
+	ns_dist3 = G_distance(window.west, window.north,
+	                      window.west, window.south);
+
+	ew_res = (ew_dist1 + ew_dist2 + ew_dist3) / (3 * window.cols);
+	ns_res = (ns_dist1 + ns_dist2 + ns_dist3) / (3 * window.rows);
+    }
+    else {
+	ns_res = window.ns_res;
+	ew_res = window.ew_res;
+    }
+    
+    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+	/* get r, c (r_nbr, c_nbr) for neighbours */
+	r_nbr = nextdr[ct_dir];
+	c_nbr = nextdc[ct_dir];
+	/* account for rare cases when ns_res != ew_res */
+	dy = ABS(r_nbr) * ns_res;
+	dx = ABS(c_nbr) * ew_res;
+	if (ct_dir < 4)
+	    dist_to_nbr[ct_dir] = dx + dy;
+	else
+	    dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+    }
+    contour[0] = contour[1] = ew_res;
+    contour[2] = contour[3] = ns_res;
+    if (sides == 8) {
+	contour[4] = contour[5] = contour[6] = contour[7] = (ew_res + ns_res) / 2.;
+    }
+    
+    return 1;
+}
+
+double get_slope_tci(CELL ele, CELL down_ele, double dist)
+{
+    if (down_ele >= ele)
+	return 0.5 / dist;
+    else
+	return (double)(ele - down_ele) / dist;
+}
 
 int do_cum(void)
 {
     int r, c, dr, dc;
-    CELL is_swale, aspect;
+    int r_nbr, c_nbr, ct_dir, np_side, edge;
+    CELL is_swale, aspect, ele_nbr;
     DCELL value, valued;
-    int killer, threshold, count;
+    int killer, threshold;
     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;
+    int this_index, down_index, nbr_index;
+    double *dist_to_nbr, *contour;
+    double tci_div;
 
     G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
 
-    count = 0;
+    /* distances to neighbours, contour lengths */
+    dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
+    contour = (double *)G_malloc(sides * sizeof(double));
+
+    get_dist(dist_to_nbr, contour);
+
     if (bas_thres <= 0)
 	threshold = 60;
     else
@@ -36,9 +111,44 @@ int do_cum(void)
 	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
 	    down_index = SEG_INDEX(wat_seg, dr, dc);
 	    value = wat[this_index];
-	    if ((int)(ABS(value) + 0.5) >= threshold)
+	    if (fabs(value) >= threshold)
 		FLAG_SET(swale, r, c);
 	    valued = wat[down_index];
+
+	    edge = 0;
+	    np_side = -1;
+	    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];
+
+		if (dr == r_nbr && dc == c_nbr)
+		    np_side = ct_dir;
+
+		/* check that neighbour is within region */
+		if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
+		    c_nbr < ncols) {
+
+		    nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
+		    ele_nbr = alt[nbr_index];
+		    if (Rast_is_c_null_value(&ele_nbr))
+			edge = 1;
+		}
+		else
+		    edge = 1;
+		if (edge)
+		    break;
+	    }
+	    /* do not distribute flow along edges, this causes artifacts */
+	    if (edge) {
+		is_swale = FLAG_GET(swale, r, c);
+		if (is_swale && aspect > 0) {
+		    aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+		    asp[this_index] = aspect;
+		}
+		continue;
+	    }
+
 	    if (value > 0) {
 		if (valued > 0)
 		    valued += value;
@@ -52,9 +162,17 @@ int do_cum(void)
 		    valued = value - valued;
 	    }
 	    wat[down_index] = valued;
-	    valued = ABS(valued) + 0.5;
+
+	    /* topographic wetness index ln(a / tan(beta)) */
+	    if (tci_flag) {
+		tci_div = contour[np_side] * 
+		       get_slope_tci(alt[this_index], alt[down_index],
+				     dist_to_nbr[np_side]);
+		tci[this_index] = log(fabs(wat[this_index]) / tci_div);
+	    }
+
 	    is_swale = FLAG_GET(swale, r, c);
-	    if (is_swale || ((int)valued) >= threshold) {
+	    if (is_swale || fabs(valued) >= threshold) {
 		FLAG_SET(swale, dr, dc);
 	    }
 	    else {
@@ -88,46 +206,42 @@ int do_cum(void)
  * before depressions/obstacles and gracefull flow divergence after 
  * depressions/obstacles
  * 
+ * Topograohic Wetness Index (TCI)
+ * TCI = SUM( (A / L_i) * weight / (tanb_i * weight))
+ * TCI = A / SUM (L_i / tanb_i)
+ * A: total catchment area
+ * L_i: contour length towards i-th neighbor
+ * tanb_i: slope = tan(b) towards i-th neighbor
+ *
  * ************************************/
 
 int do_cum_mfd(void)
 {
     int r, c, dr, dc;
     CELL is_swale;
-    DCELL value, valued;
+    DCELL value, valued, tci_div;
     int killer, threshold, count;
 
     /* MFD */
     int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
-    double *dist_to_nbr, *weight, sum_weight, max_weight;
+    double *dist_to_nbr, *contour, *weight, sum_weight, max_weight;
     int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
-    double dx, dy;
     CELL ele, ele_nbr, aspect, is_worked;
     double prop, max_acc;
     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."));
     G_debug(1, "MFD convergence factor set to %d.", c_fac);
 
-    /* distances to neighbours */
+    /* distances to neighbours, weights, contour lengths */
     dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
     weight = (double *)G_malloc(sides * sizeof(double));
-
-    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
-	/* get r, c (r_nbr, c_nbr) for neighbours */
-	r_nbr = nextdr[ct_dir];
-	c_nbr = nextdc[ct_dir];
-	/* account for rare cases when ns_res != ew_res */
-	dy = ABS(r_nbr) * window.ns_res;
-	dx = ABS(c_nbr) * window.ew_res;
-	if (ct_dir < 4)
-	    dist_to_nbr[ct_dir] = dx + dy;
-	else
-	    dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
-    }
+    contour = (double *)G_malloc(sides * sizeof(double));
+    
+    get_dist(dist_to_nbr, contour);
 
     flag_clear_all(worked);
     workedon = 0;
@@ -153,7 +267,6 @@ int do_cum_mfd(void)
 	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
 	    value = wat[this_index];
 	    down_index = SEG_INDEX(wat_seg, dr, dc);
-	    valued = wat[down_index];
 
 	    r_max = dr;
 	    c_max = dc;
@@ -176,13 +289,14 @@ int do_cum_mfd(void)
 		r_nbr = r + nextdr[ct_dir];
 		c_nbr = c + nextdc[ct_dir];
 		weight[ct_dir] = -1;
+
+		if (dr == r_nbr && dc == c_nbr)
+		    np_side = ct_dir;
+
 		/* check that neighbour is within region */
 		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 */
@@ -256,6 +370,7 @@ int do_cum_mfd(void)
 
 	    /* set flow accumulation for neighbours */
 	    max_acc = -1;
+	    tci_div = 0.;
 
 	    if (mfd_cells > 1) {
 		prop = 0.0;
@@ -289,6 +404,11 @@ int do_cum_mfd(void)
 				    valued = value * weight[ct_dir] - valued;
 			    }
 			    wat[nbr_index] = valued;
+			    
+			    if (tci_flag)
+				tci_div += contour[ct_dir] * 
+				           get_slope_tci(ele, alt[nbr_index],
+				                         dist_to_nbr[ct_dir]);
 
 			    /* get main drainage direction */
 			    if (ABS(valued) >= max_acc) {
@@ -324,6 +444,15 @@ int do_cum_mfd(void)
 			valued = value - valued;
 		}
 		wat[down_index] = valued;
+
+		if (tci_flag)
+		    tci_div = contour[np_side] * 
+			      get_slope_tci(ele, alt[down_index],
+				            dist_to_nbr[np_side]);
+	    }
+	    /* topographic wetness index ln(a / tan(beta)) */
+	    if (tci_flag) {
+		tci[this_index] = log(fabs(wat[this_index]) / tci_div);
 	    }
 
 	    /* update asp */

raster/r.watershed/ram/init_vars.c

@@ -19,9 +19,16 @@ int init_vars(int argc, char *argv[])
     int seg_idx;
 
     G_gisinit(argv[0]);
-    ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
-    ob_flag = bas_flag = seg_flag = haf_flag = arm_flag = dis_flag = 0;
-    zero = sl_flag = sg_flag = ls_flag = er_flag = bas_thres = 0;
+    /* input */
+    ele_flag = pit_flag = run_flag = ril_flag = 0;
+    /* output */
+    wat_flag = asp_flag = bas_flag = seg_flag = haf_flag = tci_flag = 0;
+    bas_thres = 0;
+    /* shed, unused */
+    arm_flag = dis_flag = 0;
+    /* RUSLE */
+    ob_flag = st_flag = sl_flag = sg_flag = ls_flag = er_flag = 0;
+    zero = 0;
     one = 1;
     nxt_avail_pt = 0;
     /* dep_flag = 0; */
@@ -41,6 +48,8 @@ int init_vars(int argc, char *argv[])
 	    ele_flag++;
 	else if (sscanf(argv[r], "accumulation=%s", wat_name) == 1)
 	    wat_flag++;
+	else if (sscanf(argv[r], "tci=%s", tci_name) == 1)
+	    tci_flag++;
 	else if (sscanf(argv[r], "drainage=%s", asp_name) == 1)
 	    asp_flag++;
 	else if (sscanf(argv[r], "depression=%s", pit_name) == 1)
@@ -138,6 +147,12 @@ int init_vars(int argc, char *argv[])
     wat =
 	(DCELL *) G_malloc(sizeof(DCELL) *
 			   size_array(&wat_seg, nrows, ncols));
+    if (tci_flag)
+	tci = (DCELL *) G_malloc(sizeof(DCELL) *
+			         size_array(&wat_seg, nrows, ncols));
+    else
+	tci = NULL;
+
     asp =
 	(CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
 

raster/r.watershed/ram/main.c

@@ -39,7 +39,7 @@ RAMSEG r_h_seg, dep_seg;
 RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
 int *astar_pts;
 CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
-DCELL *wat;
+DCELL *wat, *tci;
 int ril_fd;
 double *s_l, *s_g, *l_s;
 CELL one, zero;
@@ -57,11 +57,11 @@ char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX],
     thr_name[8];
 char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX],
     sg_name[GNAME_MAX];
-char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX],
-    dis_name[GNAME_MAX];
+char wat_name[GNAME_MAX], asp_name[GNAME_MAX], tci_name[GNAME_MAX];
+char arm_name[GNAME_MAX], dis_name[GNAME_MAX];
 char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, flat_flag;
 char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
-char bas_flag, seg_flag, haf_flag, er_flag;
+char bas_flag, seg_flag, haf_flag, er_flag, tci_flag;
 char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
 FILE *fp;