[GRASS-SVN] r51489 - grass/trunk/raster/r.watershed/ram
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Apr 23 08:44:35 EDT 2012
Author: mmetz
Date: 2012-04-23 05:44:35 -0700 (Mon, 23 Apr 2012)
New Revision: 51489
Modified:
grass/trunk/raster/r.watershed/ram/Gwater.h
grass/trunk/raster/r.watershed/ram/close_maps.c
grass/trunk/raster/r.watershed/ram/do_astar.c
grass/trunk/raster/r.watershed/ram/do_astar.h
grass/trunk/raster/r.watershed/ram/do_cum.c
grass/trunk/raster/r.watershed/ram/init_vars.c
grass/trunk/raster/r.watershed/ram/main.c
Log:
r.watershed: add TCI, remove edge artefacts for SFD (ram)
Modified: grass/trunk/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/Gwater.h 2012-04-23 11:30:28 UTC (rev 51488)
+++ grass/trunk/raster/r.watershed/ram/Gwater.h 2012-04-23 12:44:35 UTC (rev 51489)
@@ -60,7 +60,7 @@
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 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;
Modified: grass/trunk/raster/r.watershed/ram/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/close_maps.c 2012-04-23 11:30:28 UTC (rev 51488)
+++ grass/trunk/raster/r.watershed/ram/close_maps.c 2012-04-23 12:44:35 UTC (rev 51489)
@@ -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 @@
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) {
Modified: grass/trunk/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.c 2012-04-23 11:30:28 UTC (rev 51488)
+++ grass/trunk/raster/r.watershed/ram/do_astar.c 2012-04-23 12:44:35 UTC (rev 51489)
@@ -163,7 +163,7 @@
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 @@
/* 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 @@
child_idx = astar_pts[child];
while (child > 1) {
- parent = GET_PARENT(child);
+ GET_PARENT(parent, child);
elep = alt[astar_pts[parent]];
/* child smaller */
Modified: grass/trunk/raster/r.watershed/ram/do_astar.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.h 2012-04-23 11:30:28 UTC (rev 51488)
+++ grass/trunk/raster/r.watershed/ram/do_astar.h 2012-04-23 12:44:35 UTC (rev 51489)
@@ -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__ */
Modified: grass/trunk/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_cum.c 2012-04-23 11:30:28 UTC (rev 51488)
+++ grass/trunk/raster/r.watershed/ram/do_cum.c 2012-04-23 12:44:35 UTC (rev 51489)
@@ -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 @@
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 @@
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,47 +206,43 @@
* 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));
+ contour = (double *)G_malloc(sides * sizeof(double));
+
+ get_dist(dist_to_nbr, contour);
- 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);
- }
-
flag_clear_all(worked);
workedon = 0;
@@ -153,7 +267,6 @@
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 @@
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 @@
/* set flow accumulation for neighbours */
max_acc = -1;
+ tci_div = 0.;
if (mfd_cells > 1) {
prop = 0.0;
@@ -289,6 +404,11 @@
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,7 +444,16 @@
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 */
if (dr != r_max || dc != c_max) {
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2012-04-23 11:30:28 UTC (rev 51488)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2012-04-23 12:44:35 UTC (rev 51489)
@@ -19,9 +19,16 @@
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 @@
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 @@
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));
Modified: grass/trunk/raster/r.watershed/ram/main.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/main.c 2012-04-23 11:30:28 UTC (rev 51488)
+++ grass/trunk/raster/r.watershed/ram/main.c 2012-04-23 12:44:35 UTC (rev 51489)
@@ -39,7 +39,7 @@
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 @@
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;
More information about the grass-commit
mailing list