[GRASS-SVN] r51518 - grass/trunk/raster/r.watershed/seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Apr 24 08:00:01 EDT 2012
Author: mmetz
Date: 2012-04-24 05:00:01 -0700 (Tue, 24 Apr 2012)
New Revision: 51518
Modified:
grass/trunk/raster/r.watershed/seg/close_maps.c
grass/trunk/raster/r.watershed/seg/do_cum.c
Log:
improve TCI (seg)
Modified: grass/trunk/raster/r.watershed/seg/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/close_maps.c 2012-04-24 11:59:00 UTC (rev 51517)
+++ grass/trunk/raster/r.watershed/seg/close_maps.c 2012-04-24 12:00:01 UTC (rev 51518)
@@ -127,8 +127,9 @@
/* 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 @@
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 @@
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);
}
Modified: grass/trunk/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_cum.c 2012-04-24 11:59:00 UTC (rev 51517)
+++ grass/trunk/raster/r.watershed/seg/do_cum.c 2012-04-24 12:00:01 UTC (rev 51518)
@@ -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 @@
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 @@
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 @@
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 @@
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 @@
* 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 @@
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,9 +407,8 @@
/* set flow accumulation for neighbours */
max_acc = -1;
- tci_div = 0.;
+ tci_div = sum_contour = 0.;
-
if (mfd_cells > 1) {
prop = 0.0;
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -410,6 +421,13 @@
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 @@
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 @@
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 @@
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);
}
More information about the grass-commit
mailing list