[GRASS-SVN] r51517 - grass/trunk/raster/r.watershed/ram
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Apr 24 07:59:00 EDT 2012
Author: mmetz
Date: 2012-04-24 04:59:00 -0700 (Tue, 24 Apr 2012)
New Revision: 51517
Modified:
grass/trunk/raster/r.watershed/ram/close_maps.c
grass/trunk/raster/r.watershed/ram/do_cum.c
grass/trunk/raster/r.watershed/ram/init_vars.c
Log:
improve TCI (ram)
Modified: grass/trunk/raster/r.watershed/ram/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/close_maps.c 2012-04-24 09:36:25 UTC (rev 51516)
+++ grass/trunk/raster/r.watershed/ram/close_maps.c 2012-04-24 11:59:00 UTC (rev 51517)
@@ -18,7 +18,7 @@
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 @@
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 @@
/* 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 @@
}
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 @@
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/ram/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_cum.c 2012-04-24 09:36:25 UTC (rev 51516)
+++ grass/trunk/raster/r.watershed/ram/do_cum.c 2012-04-24 11:59:00 UTC (rev 51517)
@@ -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 @@
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 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 @@
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 @@
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 @@
* 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 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 @@
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 @@
/* 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 @@
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 @@
}
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 @@
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 @@
}
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 */
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2012-04-24 09:36:25 UTC (rev 51516)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2012-04-24 11:59:00 UTC (rev 51517)
@@ -314,10 +314,7 @@
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;
More information about the grass-commit
mailing list