[GRASS-SVN] r51490 - grass/trunk/raster/r.watershed/seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Apr 23 08:44:53 EDT 2012
Author: mmetz
Date: 2012-04-23 05:44:53 -0700 (Mon, 23 Apr 2012)
New Revision: 51490
Modified:
grass/trunk/raster/r.watershed/seg/Gwater.h
grass/trunk/raster/r.watershed/seg/close_maps.c
grass/trunk/raster/r.watershed/seg/cseg.h
grass/trunk/raster/r.watershed/seg/do_cum.c
grass/trunk/raster/r.watershed/seg/dseg_get.c
grass/trunk/raster/r.watershed/seg/init_vars.c
grass/trunk/raster/r.watershed/seg/main.c
Log:
r.watershed: add TCI, remove edge artefacts for SFD (seg)
Modified: grass/trunk/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/Gwater.h 2012-04-23 12:44:35 UTC (rev 51489)
+++ grass/trunk/raster/r.watershed/seg/Gwater.h 2012-04-23 12:44:53 UTC (rev 51490)
@@ -90,7 +90,7 @@
extern BSEG s_b;
extern CSEG dis, bas, haf, r_h, dep;
extern SSEG watalt, aspflag;
-extern DSEG slp, s_l, s_g, l_s, ril;
+extern DSEG tci, slp, s_l, s_g, l_s, ril;
extern double segs_mb;
extern char zero, one;
extern double ril_value, d_zero, d_one;
@@ -106,9 +106,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;
-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/seg/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/close_maps.c 2012-04-23 12:44:35 UTC (rev 51489)
+++ grass/trunk/raster/r.watershed/seg/close_maps.c 2012-04-23 12:44:53 UTC (rev 51490)
@@ -123,6 +123,85 @@
}
Rast_write_colors(wat_name, this_mapset, &colors);
}
+
+ /* TCI */
+ if (tci_flag) {
+ DCELL watvalue;
+
+ G_message(_("Closing tci map"));
+ sum = sum_sqr = stddev = 0.0;
+ dbuf = Rast_allocate_d_buf();
+ wabuf = G_malloc(ncols * sizeof(WAT_ALT));
+ dseg_flush(&tci);
+ if (!wat_flag)
+ seg_flush(&watalt);
+
+ fd = Rast_open_new(tci_name, DCELL_TYPE);
+
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 1);
+ Rast_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
+ seg_get_row(&watalt, (char *)wabuf, r);
+ for (c = 0; c < ncols; c++) {
+ dseg_get(&tci, &dvalue, r, c);
+ watvalue = wabuf[c].wat;
+ if (!Rast_is_d_null_value(&watvalue)) {
+ dbuf[c] = dvalue;
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
+ }
+ Rast_put_row(fd, dbuf, DCELL_TYPE);
+ }
+ G_percent(r, nrows, 1); /* finish it */
+
+ Rast_close(fd);
+ G_free(wabuf);
+ G_free(dbuf);
+ dseg_close(&tci);
+
+ 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 */
+ /* start with white to get more detail? NULL cells are white by default, may be confusing */
+
+ 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);
+ }
+
seg_close(&watalt);
if (asp_flag) {
G_message(_("Closing flow direction map"));
Modified: grass/trunk/raster/r.watershed/seg/cseg.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/cseg.h 2012-04-23 12:44:35 UTC (rev 51489)
+++ grass/trunk/raster/r.watershed/seg/cseg.h 2012-04-23 12:44:53 UTC (rev 51490)
@@ -81,6 +81,7 @@
/* dseg_get.c */
int dseg_get(DSEG *, double *, GW_LARGE_INT, GW_LARGE_INT);
+int dseg_flush(DSEG *);
/* dseg_open.c */
int dseg_open(DSEG *, int, int, int);
Modified: grass/trunk/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_cum.c 2012-04-23 12:44:35 UTC (rev 51489)
+++ grass/trunk/raster/r.watershed/seg/do_cum.c 2012-04-23 12:44:53 UTC (rev 51490)
@@ -4,10 +4,78 @@
#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;
+ int r_nbr, c_nbr, ct_dir, np_side;
char is_swale;
DCELL value, valued;
POINT point;
@@ -17,9 +85,17 @@
int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
WAT_ALT wa, wadown;
ASP_FLAG af, afdown;
+ double *dist_to_nbr, *contour;
+ double tci_val, tci_div;
G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
+ /* 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
@@ -41,7 +117,35 @@
FLAG_UNSET(af.flag, WORKEDFLAG);
if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) {
- /* TODO: do not distribute flow along edges, this causes artifacts */
+ np_side = -1;
+ r_nbr = dr;
+ c_nbr = dc;
+ 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) {
+
+ seg_get(&aspflag, (char *)&afdown, r_nbr, c_nbr);
+ if (FLAG_GET(afdown.flag, NULLFLAG))
+ break;
+ }
+ }
+ /* do not distribute flow along edges, this causes artifacts */
+ if (FLAG_GET(af.flag, EDGEFLAG)) {
+ if (FLAG_GET(af.flag, SWALEFLAG) && af.asp > 0) {
+ af.asp = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+ }
+ seg_put(&aspflag, (char *)&af, r, c);
+ continue;
+ }
+
seg_get(&watalt, (char *)&wa, r, c);
value = wa.wat;
is_swale = FLAG_GET(af.flag, SWALEFLAG);
@@ -65,6 +169,16 @@
}
wadown.wat = valued;
seg_put(&watalt, (char *)&wadown, dr, dc);
+
+ /* topographic wetness index ln(a / tan(beta)) */
+ if (tci_flag) {
+ 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);
+ dseg_put(&tci, &tci_val, r, c);
+ }
+
/* update asp for depression */
if (is_swale || fabs(valued) >= threshold) {
seg_get(&aspflag, (char *)&afdown, dr, dc);
@@ -112,6 +226,7 @@
{
int r, c, dr, dc;
DCELL value, valued, *wat_nbr;
+ double tci_val, tci_div;
POINT point;
WAT_ALT wa;
ASP_FLAG af, afdown;
@@ -120,9 +235,8 @@
/* MFD */
int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
- double *dist_to_nbr, *weight, sum_weight, max_weight;
- int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side, max_side;
- double dx, dy;
+ double *dist_to_nbr, *contour, *weight, sum_weight, max_weight;
+ int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
CELL ele, *ele_nbr;
double prop, max_acc;
int workedon, edge, is_swale, flat;
@@ -136,20 +250,10 @@
/* distances to neighbours */
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_nbr = (char *)G_malloc(sides * sizeof(char));
wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
ele_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
@@ -287,13 +391,13 @@
mfd_cells++;
sum_weight += max_weight;
weight[np_side] = max_weight;
- max_side = np_side;
}
/* set flow accumulation for neighbours */
max_acc = -1;
- max_side = np_side;
+ tci_div = 0.;
+
if (mfd_cells > 1) {
prop = 0.0;
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -327,12 +431,16 @@
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]);
r_max = r_nbr;
c_max = c_nbr;
- max_side = ct_dir;
}
}
else if (ct_dir == np_side) {
@@ -372,8 +480,19 @@
wa.wat = valued;
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],
+ dist_to_nbr[np_side]);
}
+ /* topographic wetness index ln(a / tan(beta)) */
+ if (tci_flag) {
+ tci_val = log(fabs(value) / tci_div);
+ dseg_put(&tci, &tci_val, r, c);
+ }
+
/* update asp */
if (dr != r_max || dc != c_max) {
if (af.asp < 0) {
Modified: grass/trunk/raster/r.watershed/seg/dseg_get.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/dseg_get.c 2012-04-23 12:44:35 UTC (rev 51489)
+++ grass/trunk/raster/r.watershed/seg/dseg_get.c 2012-04-23 12:44:53 UTC (rev 51490)
@@ -10,3 +10,9 @@
}
return 0;
}
+
+int dseg_flush(DSEG * dseg)
+{
+ segment_flush(&(dseg->seg));
+ return 0;
+}
Modified: grass/trunk/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/init_vars.c 2012-04-23 12:44:35 UTC (rev 51489)
+++ grass/trunk/raster/r.watershed/seg/init_vars.c 2012-04-23 12:44:53 UTC (rev 51490)
@@ -29,9 +29,15 @@
int ct_dir, r_nbr, c_nbr;
G_gisinit(argv[0]);
- ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
- st_flag = bas_flag = seg_flag = haf_flag = arm_flag = dis_flag = 0;
- ob_flag = 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;
nxt_avail_pt = 0;
/* dep_flag = 0; */
max_length = d_zero = 0.0;
@@ -51,6 +57,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)
@@ -173,6 +181,11 @@
/* heap points: / 4 */
memory_divisor += sizeof(HEAP_PNT) / 4.;
disk_space += sizeof(HEAP_PNT);
+ /* TCI: as is */
+ if (tci_flag) {
+ memory_divisor += sizeof(double);
+ disk_space += sizeof(double);
+ }
/* RUSLE */
if (er_flag) {
/* r_h */
@@ -243,6 +256,9 @@
seg_open(&watalt, nrows, ncols, seg_rows, seg_cols, num_open_segs * 2, sizeof(WAT_ALT));
seg_open(&aspflag, nrows, ncols, seg_rows, seg_cols, num_open_segs * 4, sizeof(ASP_FLAG));
+ if (tci_flag)
+ dseg_open(&tci, seg_rows, seg_cols, num_open_segs);
+
/* open elevation input */
ele_fd = Rast_open_old(ele_name, "");
@@ -385,7 +401,7 @@
}
if (ril_flag) {
- dseg_open(&ril, 1, seg_rows * seg_cols, num_open_segs);
+ dseg_open(&ril, seg_rows, seg_cols, num_open_segs);
dseg_read_cell(&ril, ril_name, "");
}
@@ -393,9 +409,9 @@
dseg_open(&s_l, seg_rows, seg_cols, num_open_segs);
if (sg_flag)
- dseg_open(&s_g, 1, seg_rows * seg_cols, num_open_segs);
+ dseg_open(&s_g, seg_rows, seg_cols, num_open_segs);
if (ls_flag)
- dseg_open(&l_s, 1, seg_rows * seg_cols, num_open_segs);
+ dseg_open(&l_s, seg_rows, seg_cols, num_open_segs);
}
G_debug(1, "open segments for A* points");
Modified: grass/trunk/raster/r.watershed/seg/main.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/main.c 2012-04-23 12:44:35 UTC (rev 51489)
+++ grass/trunk/raster/r.watershed/seg/main.c 2012-04-23 12:44:53 UTC (rev 51490)
@@ -39,7 +39,7 @@
BSEG s_b;
CSEG dis, alt, bas, haf, r_h, dep;
SSEG watalt, aspflag;
-DSEG slp, s_l, s_g, l_s, ril;
+DSEG tci, slp, s_l, s_g, l_s, ril;
double segs_mb;
char zero, one;
double ril_value, d_zero, d_one;
@@ -56,11 +56,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;
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