improve TCI (ram)

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

raster/r.watershed/ram/close_maps.c

@@ -18,7 +18,7 @@ int close_maps(void)
 
     if (asp_flag || dis_flag)
 	buf = Rast_allocate_c_buf();
-    if (wat_flag || ls_flag || sl_flag || sg_flag)
+    if (wat_flag || ls_flag || sl_flag || sg_flag || tci_flag)
 	dbuf = Rast_allocate_d_buf();
     G_free(alt);
     if (ls_flag || sg_flag)
@@ -49,7 +49,7 @@ int close_maps(void)
 		Rast_set_d_null_value(dbuf, ncols);	/* reset row to all NULL */
 		for (c = 0; c < ncols; c++) {
 		    dvalue = wat[SEG_INDEX(wat_seg, r, c)];
-		    if (Rast_is_d_null_value(&dvalue) == 0 && dvalue) {
+		    if (!Rast_is_d_null_value(&dvalue)) {
 			dbuf[c] = dvalue;
 			dvalue = ABS(dvalue);
 			sum += dvalue;
@@ -132,6 +132,7 @@ int close_maps(void)
     /* TCI */
     if (tci_flag) {
 	DCELL watvalue;
+	double mean;
 
 	sum = sum_sqr = stddev = 0.0;
 	fd = Rast_open_new(tci_name, DCELL_TYPE);
@@ -150,6 +151,7 @@ int close_maps(void)
 	}
 	Rast_close(fd);
 
+	mean = sum / do_points;
 	stddev =
 	    sqrt((sum_sqr - (sum + sum / do_points)) / (do_points - 1));
 	G_debug(1, "stddev: %f", stddev);
@@ -164,31 +166,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/ram/do_cum.c

@@ -3,11 +3,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;
@@ -50,17 +49,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)
@@ -82,7 +81,7 @@ int do_cum(void)
     int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
     int this_index, down_index, nbr_index;
     double *dist_to_nbr, *contour;
-    double tci_div;
+    double tci_div, cell_size;
 
     G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
 
@@ -90,7 +89,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;
@@ -168,7 +167,7 @@ int do_cum(void)
 		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);
+		tci[this_index] = log((fabs(wat[this_index]) * cell_size) / tci_div);
 	    }
 
 	    is_swale = FLAG_GET(swale, r, c);
@@ -206,12 +205,18 @@ 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)
+ * 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 neighbor
- * tanb_i: slope = tan(b) towards i-th neighbor
+ * 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
  *
  * ************************************/
 
@@ -219,7 +224,7 @@ int do_cum_mfd(void)
 {
     int r, c, dr, dc;
     CELL is_swale;
-    DCELL value, valued, tci_div;
+    DCELL value, valued, tci_div, sum_contour, cell_size;
     int killer, threshold, count;
 
     /* MFD */
@@ -241,7 +246,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_clear_all(worked);
     workedon = 0;
@@ -370,7 +375,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;
@@ -386,6 +391,13 @@ int do_cum_mfd(void)
 
 			    nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
 
+			    if (tci_flag) {
+				sum_contour += contour[ct_dir];
+				tci_div += get_slope_tci(ele, alt[nbr_index],
+				                         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];
@@ -405,11 +417,6 @@ int do_cum_mfd(void)
 			    }
 			    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) {
 				max_acc = ABS(valued);
@@ -427,6 +434,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;
 	    }
 
 	    if (mfd_cells < 2) {
@@ -445,14 +454,16 @@ int do_cum_mfd(void)
 		}
 		wat[down_index] = valued;
 
-		if (tci_flag)
-		    tci_div = contour[np_side] * 
-			      get_slope_tci(ele, alt[down_index],
+		if (tci_flag) {
+		    sum_contour = contour[np_side];
+		    tci_div = 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);
+		tci[this_index] = log((fabs(wat[this_index]) * cell_size) /
+		                      (sum_contour * tci_div));
 	    }
 
 	    /* update asp */

raster/r.watershed/ram/init_vars.c

@@ -314,10 +314,7 @@ int init_vars(int argc, char *argv[])
 	    Rast_get_c_row(fd, buf, r);
 	for (c = 0; c < ncols; c++) {
 	    seg_idx = SEG_INDEX(wat_seg, r, c);
-	    if (FLAG_GET(worked, r, c)) {
-		wat[seg_idx] = 0;
-	    }
-	    else {
+	    if (!(FLAG_GET(worked, r, c))) {
 		asp_value = asp[seg_idx];
 		if (er_flag)
 		    s_l[seg_idx] = half_res;