improve TCI (seg)

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

raster/r.watershed/seg/close_maps.c

@@ -127,8 +127,9 @@ int close_maps(void)
     /* TCI */
     if (tci_flag) {
 	DCELL watvalue;
+	double mean;
 
-	G_message(_("Closing tci map"));
+	G_message(_("Closing TCI map"));
 	sum = sum_sqr = stddev = 0.0;
 	dbuf = Rast_allocate_d_buf();
 	wabuf = G_malloc(ncols * sizeof(WAT_ALT));
@@ -160,6 +161,7 @@ int close_maps(void)
 	G_free(dbuf);
 	dseg_close(&tci);
 
+	mean = sum / do_points;
 	stddev = sqrt((sum_sqr - (sum + sum / do_points)) / (do_points - 1));
 	G_debug(1, "stddev: %f", stddev);
 
@@ -174,31 +176,36 @@ int close_maps(void)
 
 	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;
+	if (min - 1 < mean - 0.5 * stddev) {
+	    clr_min = min - 1;
+	    clr_max = mean - 0.5 * stddev;
+	    Rast_add_d_color_rule(&clr_min, 255, 255, 0, &clr_max, 255,
+					  255, 0, &colors);
+	}
+
+	clr_min = mean - 0.5 * stddev;
+	clr_max = mean - 0.2 * stddev;
 	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;
+	clr_max = mean + 0.2 * stddev;
 	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;
+	clr_max = mean + 0.6 * stddev;
 	Rast_add_d_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
 				  0, 255, &colors);
 	clr_min = clr_max;
-	clr_max = max + 1.;
+	clr_max = mean + 1. * stddev;
 	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);
+	if (max > 0 && max > clr_max) {
+	    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);
     }
 

raster/r.watershed/seg/do_cum.c

@@ -4,11 +4,10 @@
 #include <grass/raster.h>
 #include <grass/glocale.h>
 
-int get_dist(double *dist_to_nbr, double *contour)
+double get_dist(double *dist_to_nbr, double *contour)
 {
     int ct_dir, r_nbr, c_nbr;
-    double dx, dy;
-    double ns_res, ew_res;
+    double dx, dy, ns_res, ew_res;
 
     if (G_projection() == PROJECTION_LL) {
 	double ew_dist1, ew_dist2, ew_dist3;
@@ -51,17 +50,17 @@ int get_dist(double *dist_to_nbr, double *contour)
 	dy = ABS(r_nbr) * ns_res;
 	dx = ABS(c_nbr) * ew_res;
 	if (ct_dir < 4)
-	    dist_to_nbr[ct_dir] = dx + dy;
+	    dist_to_nbr[ct_dir] = (dx + dy) * ele_scale;
 	else
-	    dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+	    dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy) * ele_scale;
     }
-    contour[0] = contour[1] = ew_res;
-    contour[2] = contour[3] = ns_res;
+    contour[0] = contour[1] = ew_res / 2.;
+    contour[2] = contour[3] = ns_res / 2.;
     if (sides == 8) {
-	contour[4] = contour[5] = contour[6] = contour[7] = (ew_res + ns_res) / 2.;
+	contour[4] = contour[5] = contour[6] = contour[7] = sqrt(ew_res * ns_res) / 2.;
     }
     
-    return 1;
+    return ew_res * ns_res;
 }
 
 double get_slope_tci(CELL ele, CELL down_ele, double dist)
@@ -86,7 +85,7 @@ int do_cum(void)
     WAT_ALT wa, wadown;
     ASP_FLAG af, afdown;
     double *dist_to_nbr, *contour;
-    double tci_val, tci_div; 
+    double tci_val, tci_div, cell_size;
 
     G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
 
@@ -94,7 +93,7 @@ int do_cum(void)
     dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
     contour = (double *)G_malloc(sides * sizeof(double));
 
-    get_dist(dist_to_nbr, contour);
+    cell_size = get_dist(dist_to_nbr, contour);
 
     if (bas_thres <= 0)
 	threshold = 60;
@@ -175,7 +174,7 @@ int do_cum(void)
 		tci_div = contour[np_side] * 
 		       get_slope_tci(wa.ele, wadown.ele,
 				     dist_to_nbr[np_side]);
-		tci_val = log(fabs(wa.wat) / tci_div);
+		tci_val = log((fabs(wa.wat) * cell_size) / tci_div);
 		dseg_put(&tci, &tci_val, r, c);
 	    }
 
@@ -220,13 +219,26 @@ int do_cum(void)
  * before depressions/obstacles and gracefull flow divergence after 
  * depressions/obstacles
  * 
+ * Topographic Convergence Index (TCI)
+ * after Quinn et al. (1991), modified and adapted for the modified 
+ * Holmgren MFD algorithm
+ * TCI: specific catchment area divided by tangens of slope
+ * specific catchment area: total catchment area divided by contour line
+ * TCI for D8:     A / (L * tanb)
+ * TCI for MFD:    A / (SUM(L_i) * (SUM(tanb_i * weight_i) / SUM(weight_i))
+ * 
+ * A: total catchment area
+ * L_i: contour length towards i_th cell
+ * tanb_i: slope = tan(b) towards i_th cell
+ * weight_i: weight for flow distribution towards i_th cell
+ *
  * ************************************/
 
 int do_cum_mfd(void)
 {
     int r, c, dr, dc;
     DCELL value, valued, *wat_nbr;
-    double tci_val, tci_div;
+    double tci_val, tci_div, sum_contour, cell_size;
     POINT point;
     WAT_ALT wa;
     ASP_FLAG af, afdown;
@@ -252,7 +264,7 @@ int do_cum_mfd(void)
     weight = (double *)G_malloc(sides * sizeof(double));
     contour = (double *)G_malloc(sides * sizeof(double));
     
-    get_dist(dist_to_nbr, contour);
+    cell_size = get_dist(dist_to_nbr, contour);
 
     flag_nbr = (char *)G_malloc(sides * sizeof(char));
     wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
@@ -395,8 +407,7 @@ int do_cum_mfd(void)
 
 	    /* set flow accumulation for neighbours */
 	    max_acc = -1;
-	    tci_div = 0.;
-
+	    tci_div = sum_contour = 0.;
 
 	    if (mfd_cells > 1) {
 		prop = 0.0;
@@ -410,6 +421,13 @@ int do_cum_mfd(void)
 
 			if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
 
+			    if (tci_flag) {
+				sum_contour += contour[ct_dir];
+				tci_div += get_slope_tci(ele, ele_nbr[ct_dir],
+				                         dist_to_nbr[ct_dir]) *
+					   weight[ct_dir];
+			    }
+
 			    weight[ct_dir] = weight[ct_dir] / sum_weight;
 			    /* check everything adds up to 1.0 */
 			    prop += weight[ct_dir];
@@ -431,11 +449,6 @@ int do_cum_mfd(void)
 			    wa.ele = ele_nbr[ct_dir];
 			    seg_put(&watalt, (char *)&wa, r_nbr, c_nbr);
 
-			    if (tci_flag)
-				tci_div += contour[ct_dir] * 
-				           get_slope_tci(ele, ele_nbr[ct_dir],
-				                         dist_to_nbr[ct_dir]);
-
 			    /* get main drainage direction */
 			    if (fabs(wat_nbr[ct_dir]) >= max_acc) {
 				max_acc = ABS(wat_nbr[ct_dir]);
@@ -461,6 +474,8 @@ int do_cum_mfd(void)
 		    G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
 			      prop);
 		}
+		if (tci_flag)
+		    tci_div /= sum_weight;
 	    }
 	    /* SFD-like accumulation */
 	    else {
@@ -481,15 +496,17 @@ int do_cum_mfd(void)
 		wa.ele = ele_nbr[np_side];
 		seg_put(&watalt, (char *)&wa, dr, dc);
 
-		if (tci_flag)
-		    tci_div = contour[np_side] * 
-			      get_slope_tci(ele, ele_nbr[np_side],
+		if (tci_flag) {
+		    sum_contour = contour[np_side];
+		    tci_div = get_slope_tci(ele, ele_nbr[np_side],
 				            dist_to_nbr[np_side]);
+		}
 	    }
 
 	    /* topographic wetness index ln(a / tan(beta)) */
 	    if (tci_flag) {
-		tci_val = log(fabs(value) / tci_div);
+		tci_val = log((fabs(value) * cell_size) /
+		              (sum_contour * tci_div));
 		dseg_put(&tci, &tci_val, r, c);
 	    }