[GRASS-SVN] r67298 - in grass/trunk/raster/r.watershed: front ram seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Dec 21 02:06:57 PST 2015
Author: mmetz
Date: 2015-12-21 02:06:57 -0800 (Mon, 21 Dec 2015)
New Revision: 67298
Modified:
grass/trunk/raster/r.watershed/front/main.c
grass/trunk/raster/r.watershed/front/r.watershed.html
grass/trunk/raster/r.watershed/ram/Gwater.h
grass/trunk/raster/r.watershed/ram/close_maps.c
grass/trunk/raster/r.watershed/ram/close_maps2.c
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
grass/trunk/raster/r.watershed/seg/Gwater.h
grass/trunk/raster/r.watershed/seg/close_maps.c
grass/trunk/raster/r.watershed/seg/do_cum.c
grass/trunk/raster/r.watershed/seg/init_vars.c
grass/trunk/raster/r.watershed/seg/main.c
Log:
r.watershed: +SPI
Modified: grass/trunk/raster/r.watershed/front/main.c
===================================================================
--- grass/trunk/raster/r.watershed/front/main.c 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/front/main.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -58,6 +58,7 @@
struct Option *opt15;
struct Option *opt16;
struct Option *opt17;
+ struct Option *opt18;
struct Flag *flag_sfd;
struct Flag *flag_flow;
struct Flag *flag_seg;
@@ -139,6 +140,13 @@
opt17->required = NO;
opt17->guisection = _("Outputs");
+ opt18 = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt18->key = "spi";
+ opt18->label =
+ _("Stream power index a * tan(b)");
+ opt18->required = NO;
+ opt18->guisection = _("Outputs");
+
opt9 = G_define_standard_option(G_OPT_R_OUTPUT);
opt9->key = "drainage";
opt9->description = _("Name for output drainage direction raster map");
@@ -302,6 +310,7 @@
do_opt(opt7);
do_opt(opt8);
do_opt(opt17);
+ do_opt(opt18);
do_opt(opt9);
do_opt(opt10);
do_opt(opt11);
@@ -331,6 +340,10 @@
write_hist(opt17->answer,
"Watershed accumulation: topographic index ln(a / tan b)",
opt1->answer, flag_seg->answer, flag_sfd->answer);
+ if (opt18->answer)
+ write_hist(opt18->answer,
+ "Watershed accumulation: stream power index a * tan b",
+ opt1->answer, flag_seg->answer, flag_sfd->answer);
if (opt9->answer)
write_hist(opt9->answer,
"Watershed drainage direction (CCW from East divided by 45deg)",
Modified: grass/trunk/raster/r.watershed/front/r.watershed.html
===================================================================
--- grass/trunk/raster/r.watershed/front/r.watershed.html 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/front/r.watershed.html 2015-12-21 10:06:57 UTC (rev 67298)
@@ -113,17 +113,29 @@
their surface runoff and sedimentation yields calculated accurately.
<p>
-Output <b>tci</b> raster map contains topographic index TCI is
+Output <b>tci</b> raster map contains topographic index TCI,
computed as
<tt>ln(α / tan(β))</tt> where α is the cumulative
upslope area draining through a point per unit contour length and
<tt>tan(β)</tt> is the local slope angle. The TCI reflects the
tendency of water to accumulate at any point in the catchment and the
-tendency for gravitaional forces to move that water downslope (Quinn
+tendency for gravitational forces to move that water downslope (Quinn
et al. 1991). This value will be negative if <tt>α /
tan(β) < 1</tt>.
<p>
+Output <b>spi</b> raster map contains stream power index SPI,
+computed as
+<tt>α * tan(β)</tt> where α is the cumulative
+upslope area draining through a point per unit contour length and
+<tt>tan(β)</tt> is the local slope angle. The SPI reflects the
+power of ater flow at any point in the catchment and the
+tendency for gravitational forces to move that water downslope (Moore
+et al. 1991). This value will be negative if <tt>α < 0</tt>,
+i.e. for cells with possible surface runoff from outside of the current
+geographic region.
+
+<p>
Output <b>drainage</b> raster map contains drainage direction.
Provides the "aspect" for each cell measured CCW from East.
Multiplying positive values by 45 will give the direction in degrees
@@ -486,6 +498,12 @@
cost path search</i>, <b>Hydrol. Earth Syst. Sci.</b> Vol 15, 667-678.<br>
DOI: <a href="http://dx.doi.org/10.5194/hess-15-667-2011">10.5194/hess-15-667-2011</a>
+<li>
+Moore I.D., Grayson R.B., Ladson A.R. (1991). <i>Digital terrain
+modelling: a review of hydrogical, geomorphological, and biological
+applications</i>, <b>Hydrological Processes</b>, Vol 5(1), 3-30<br>
+DOI: <a href="http://dx.doi.org/10.1002/hyp.3360050103">10.1002/hyp.3360050103</a>
+
<li>Quinn P., K. Beven K., Chevallier P., Planchon O. (1991). <i>The
prediction of hillslope flow paths for distributed hydrological modelling
using Digital Elevation Models</i>, <b>Hydrological Processes</b> Vol 5(1),
Modified: grass/trunk/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/Gwater.h 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/ram/Gwater.h 2015-12-21 10:06:57 UTC (rev 67298)
@@ -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, *tci;
+extern DCELL *wat, *sca, *tanb;
extern int ril_fd;
extern double *s_l, *s_g, *l_s;
extern CELL one, zero;
@@ -76,10 +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], tci_name[GNAME_MAX];
+extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX], tci_name[GNAME_MAX], spi_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, tci_flag;
+extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag, tci_flag, spi_flag, atanb_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 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/ram/close_maps.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -131,7 +131,6 @@
/* TCI */
if (tci_flag) {
- DCELL watvalue;
double mean;
sum = sum_sqr = stddev = 0.0;
@@ -139,9 +138,9 @@
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)) {
+ if (!Rast_is_d_null_value(&tanb[SEG_INDEX(wat_seg, r, c)])) {
+ dvalue = log(sca[SEG_INDEX(wat_seg, r, c)] /
+ tanb[SEG_INDEX(wat_seg, r, c)]);
dbuf[c] = dvalue;
sum += dvalue;
sum_sqr += dvalue * dvalue;
@@ -199,6 +198,78 @@
Rast_write_colors(tci_name, this_mapset, &colors);
}
+ /* SPI */
+ if (spi_flag) {
+ double mean;
+
+ sum = sum_sqr = stddev = 0.0;
+ fd = Rast_open_new(spi_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++) {
+ if (!Rast_is_d_null_value(&tanb[SEG_INDEX(wat_seg, r, c)])) {
+ dvalue = sca[SEG_INDEX(wat_seg, r, c)] * tanb[SEG_INDEX(wat_seg, r, c)];
+ dbuf[c] = dvalue;
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
+ }
+ Rast_put_row(fd, dbuf, DCELL_TYPE);
+ }
+ Rast_close(fd);
+
+ mean = sum / do_points;
+ 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(spi_name, this_mapset, &accRange);
+ min = max = 0;
+ Rast_get_fp_range_min_max(&accRange, &min, &max);
+
+ Rast_init_colors(&colors);
+
+ 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 = 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 = 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 = mean + 1. * stddev;
+ Rast_add_d_color_rule(&clr_min, 0, 0, 255, &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(spi_name, this_mapset, &colors);
+ }
+ if (atanb_flag) {
+ G_free(sca);
+ G_free(tanb);
+ }
+
/* 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/close_maps2.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/close_maps2.c 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/ram/close_maps2.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -71,7 +71,7 @@
}
G_percent(r - 1, max, 3); /* finish it */
}
- else
+ else if (max >= 1000)
G_debug(1,
"Too many subbasins to reasonably check for color brightness");
/* using the existing stack of while/for/for/for/while loops ... */
Modified: grass/trunk/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_cum.c 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/ram/do_cum.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -109,7 +109,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, cell_size;
+ double cell_size;
G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
@@ -173,6 +173,8 @@
aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
asp[this_index] = aspect;
}
+ if (valued > 0)
+ wat[down_index] = -valued;
continue;
}
@@ -190,12 +192,14 @@
}
wat[down_index] = valued;
- /* 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]) * cell_size) / tci_div);
+ /* topographic wetness index ln(a / tan(beta)) and
+ * stream power index a * tan(beta) */
+ if (atanb_flag) {
+ sca[this_index] = fabs(wat[this_index]) *
+ (cell_size / contour[np_side]);
+ tanb[this_index] = get_slope_tci(alt[this_index],
+ alt[down_index],
+ dist_to_nbr[np_side]);
}
is_swale = FLAG_GET(swale, r, c);
@@ -357,6 +361,8 @@
if (dr == r_nbr && dc == c_nbr) {
astar_not_set = 0;
}
+ if (value < 0 && valued > 0)
+ wat[nbr_index] = -valued;
}
}
}
@@ -401,17 +407,17 @@
nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
- if (tci_flag) {
+ weight[ct_dir] = weight[ct_dir] / sum_weight;
+ /* check everything adds up to 1.0 */
+ prop += weight[ct_dir];
+
+ if (atanb_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];
-
valued = wat[nbr_index];
if (value > 0) {
if (valued > 0)
@@ -437,11 +443,9 @@
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) {
+ /* SFD-like accumulation */
+ else {
valued = wat[down_index];
if (value > 0) {
if (valued > 0)
@@ -457,16 +461,18 @@
}
wat[down_index] = valued;
- if (tci_flag) {
+ if (atanb_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]) * cell_size) /
- (sum_contour * tci_div));
+ /* topographic wetness index ln(a / tan(beta)) and
+ * stream power index a * tan(beta) */
+ if (atanb_flag) {
+ sca[this_index] = fabs(wat[this_index]) *
+ (cell_size / sum_contour);
+ tanb[this_index] = tci_div;
}
}
}
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -22,7 +22,8 @@
/* input */
ele_flag = pit_flag = run_flag = ril_flag = 0;
/* output */
- wat_flag = asp_flag = bas_flag = seg_flag = haf_flag = tci_flag = 0;
+ wat_flag = asp_flag = tci_flag = spi_flag = atanb_flag = 0;
+ bas_flag = seg_flag = haf_flag = 0;
bas_thres = 0;
/* shed, unused */
arm_flag = dis_flag = 0;
@@ -50,6 +51,8 @@
wat_flag++;
else if (sscanf(argv[r], "tci=%s", tci_name) == 1)
tci_flag++;
+ else if (sscanf(argv[r], "spi=%s", spi_name) == 1)
+ spi_flag++;
else if (sscanf(argv[r], "drainage=%s", asp_name) == 1)
asp_flag++;
else if (sscanf(argv[r], "depression=%s", pit_name) == 1)
@@ -147,11 +150,16 @@
wat =
(DCELL *) G_malloc(sizeof(DCELL) *
size_array(&wat_seg, nrows, ncols));
- if (tci_flag)
- tci = (DCELL *) G_malloc(sizeof(DCELL) *
+
+ sca = tanb = NULL;
+ atanb_flag = 0;
+ if (tci_flag || spi_flag) {
+ sca = (DCELL *) G_malloc(sizeof(DCELL) *
size_array(&wat_seg, nrows, ncols));
- else
- tci = NULL;
+ tanb = (DCELL *) G_malloc(sizeof(DCELL) *
+ size_array(&wat_seg, nrows, ncols));
+ atanb_flag = 1;
+ }
asp =
(CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
@@ -223,6 +231,10 @@
if (er_flag) {
r_h[seg_idx] = alt_value;
}
+ if (atanb_flag) {
+ Rast_set_d_null_value(&sca[seg_idx], 1);
+ Rast_set_d_null_value(&tanb[seg_idx], 1);
+ }
ptr = G_incr_void_ptr(ptr, ele_size);
}
}
Modified: grass/trunk/raster/r.watershed/ram/main.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/main.c 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/ram/main.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -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, *tci;
+DCELL *wat, *sca, *tanb;
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], tci_name[GNAME_MAX];
+char wat_name[GNAME_MAX], asp_name[GNAME_MAX], tci_name[GNAME_MAX], spi_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, tci_flag;
+char bas_flag, seg_flag, haf_flag, er_flag, tci_flag, spi_flag, atanb_flag;
char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
FILE *fp;
Modified: grass/trunk/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/Gwater.h 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/seg/Gwater.h 2015-12-21 10:06:57 UTC (rev 67298)
@@ -63,6 +63,12 @@
DCELL wat;
};
+#define A_TANB struct sca_tanb
+A_TANB {
+ DCELL sca;
+ DCELL tanb;
+};
+
#define ASP_FLAG struct aspect_flag
ASP_FLAG {
char asp;
@@ -90,7 +96,8 @@
extern BSEG s_b;
extern CSEG dis, bas, haf, r_h, dep;
extern SSEG watalt, aspflag;
-extern DSEG tci, slp, s_l, s_g, l_s, ril;
+extern DSEG slp, s_l, s_g, l_s, ril;
+extern SSEG atanb;
extern double segs_mb;
extern char zero, one;
extern double ril_value, d_zero, d_one;
@@ -106,10 +113,11 @@
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], tci_name[GNAME_MAX];
+extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX];
+extern char tci_name[GNAME_MAX], spi_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, tci_flag;
+extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag, tci_flag, spi_flag, atanb_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 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/seg/close_maps.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -124,89 +124,175 @@
Rast_write_colors(wat_name, this_mapset, &colors);
}
- /* TCI */
- if (tci_flag) {
- DCELL watvalue;
+ /* TCI, SPI */
+ if (atanb_flag) {
+ DCELL *dbuf2;
+ A_TANB sca_tanb;
double mean;
+ double sum2, sum_sqr2;
+ int fd2;
- G_message(_("Closing TCI map"));
+ if (tci_flag && spi_flag)
+ G_message(_("Closing TCI and SPI maps"));
+ else if (tci_flag && !spi_flag)
+ G_message(_("Closing TCI map"));
+ else if (!tci_flag && spi_flag)
+ G_message(_("Closing SPI map"));
+
sum = sum_sqr = stddev = 0.0;
- dbuf = Rast_allocate_d_buf();
+ sum2 = sum_sqr2 = 0.0;
wabuf = G_malloc(ncols * sizeof(WAT_ALT));
- dseg_flush(&tci);
+ seg_flush(&atanb);
if (!wat_flag)
seg_flush(&watalt);
- fd = Rast_open_new(tci_name, DCELL_TYPE);
+ fd = fd2 = -1;
+ dbuf = dbuf2 = NULL;
+ if (tci_flag) {
+ fd = Rast_open_new(tci_name, DCELL_TYPE);
+ dbuf = Rast_allocate_d_buf();
+ }
+ if (spi_flag) {
+ fd2 = Rast_open_new(spi_name, DCELL_TYPE);
+ dbuf2 = Rast_allocate_d_buf();
+ }
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);
+ if (tci_flag)
+ Rast_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
+ if (spi_flag)
+ Rast_set_d_null_value(dbuf2, ncols); /* reset row to all NULL */
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;
+ seg_get(&atanb, (char *)&sca_tanb, r, c);
+ if (!Rast_is_d_null_value(&sca_tanb.tanb)) {
+
+ if (tci_flag) {
+ dvalue = log(sca_tanb.sca / sca_tanb.tanb);
+ dbuf[c] = dvalue;
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
+ if (spi_flag) {
+ dvalue = sca_tanb.sca * sca_tanb.tanb;
+ dbuf2[c] = dvalue;
+ sum2 += dvalue;
+ sum_sqr2 += dvalue * dvalue;
+ }
}
}
- Rast_put_row(fd, dbuf, DCELL_TYPE);
+ if (tci_flag)
+ Rast_put_row(fd, dbuf, DCELL_TYPE);
+ if (spi_flag)
+ Rast_put_row(fd2, dbuf2, DCELL_TYPE);
}
G_percent(r, nrows, 1); /* finish it */
- Rast_close(fd);
G_free(wabuf);
- G_free(dbuf);
- dseg_close(&tci);
+ seg_close(&atanb);
- mean = sum / do_points;
- stddev = sqrt((sum_sqr - (sum + sum / do_points)) / (do_points - 1));
- G_debug(1, "stddev: %f", stddev);
+ if (tci_flag) {
+ Rast_close(fd);
+ G_free(dbuf);
- /* 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 */
+ mean = sum / do_points;
+ stddev = sqrt((sum_sqr - (sum + sum / do_points)) / (do_points - 1));
+ G_debug(1, "stddev: %f", stddev);
- lstddev = log(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 */
- Rast_read_fp_range(tci_name, this_mapset, &accRange);
- min = max = 0;
- Rast_get_fp_range_min_max(&accRange, &min, &max);
+ lstddev = log(stddev);
- Rast_init_colors(&colors);
+ Rast_read_fp_range(tci_name, this_mapset, &accRange);
+ min = max = 0;
+ Rast_get_fp_range_min_max(&accRange, &min, &max);
- 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);
+ Rast_init_colors(&colors);
+
+ 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 = 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 = 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 = mean + 1. * stddev;
+ Rast_add_d_color_rule(&clr_min, 0, 0, 255, &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);
}
+ if (spi_flag) {
+ Rast_close(fd2);
+ G_free(dbuf2);
- 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 = 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 = 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 = mean + 1. * stddev;
- Rast_add_d_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
- 0, &colors);
+ mean = sum2 / do_points;
+ stddev = sqrt((sum_sqr2 - (sum2 + sum2 / do_points)) / (do_points - 1));
+ G_debug(1, "stddev: %f", stddev);
- if (max > 0 && max > clr_max) {
+ /* 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(spi_name, this_mapset, &accRange);
+ min = max = 0;
+ Rast_get_fp_range_min_max(&accRange, &min, &max);
+
+ Rast_init_colors(&colors);
+
+ 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 = max + 1;
- Rast_add_d_color_rule(&clr_min, 0, 0, 0, &clr_max, 0, 0,
+ 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 = 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 = mean + 1. * stddev;
+ Rast_add_d_color_rule(&clr_min, 0, 0, 255, &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(spi_name, this_mapset, &colors);
}
- Rast_write_colors(tci_name, this_mapset, &colors);
}
seg_close(&watalt);
Modified: grass/trunk/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_cum.c 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/seg/do_cum.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -112,8 +112,9 @@
int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
WAT_ALT wa, wadown;
ASP_FLAG af, afdown;
+ A_TANB sca_tanb;
double *dist_to_nbr, *contour;
- double tci_val, tci_div, cell_size;
+ double cell_size;
G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
@@ -162,6 +163,12 @@
af.asp = -1 * drain[r - dr + 1][c - dc + 1];
}
seg_put(&aspflag, (char *)&af, r, c);
+ seg_get(&watalt, (char *)&wadown, dr, dc);
+ valued = wadown.wat;
+ if (valued > 0) {
+ wadown.wat = -valued;
+ seg_put(&watalt, (char *)&wadown, dr, dc);
+ }
continue;
}
@@ -189,13 +196,15 @@
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) * cell_size) / tci_div);
- dseg_put(&tci, &tci_val, r, c);
+ /* topographic wetness index ln(a / tan(beta)) and
+ * stream power index a * tan(beta) */
+ if (atanb_flag) {
+ sca_tanb.sca = fabs(wa.wat) *
+ (cell_size / contour[np_side]);
+
+ sca_tanb.tanb = get_slope_tci(wa.ele, wadown.ele,
+ dist_to_nbr[np_side]);
+ seg_put(&atanb, (char *)&sca_tanb, r, c);
}
/* update asp for depression */
@@ -262,10 +271,11 @@
{
int r, c, dr, dc;
DCELL value, valued, *wat_nbr;
- double tci_val, tci_div, sum_contour, cell_size;
+ double sum_contour, cell_size;
POINT point;
WAT_ALT wa;
ASP_FLAG af, afdown;
+ A_TANB sca_tanb;
GW_LARGE_INT killer;
int threshold;
@@ -384,6 +394,10 @@
if (dr == r_nbr && dc == c_nbr) {
astar_not_set = 0;
}
+ if (value < 0 && wat_nbr[ct_dir] > 0) {
+ wa.wat = -wat_nbr[ct_dir];
+ seg_put(&watalt, (char *)&wa, r_nbr, c_nbr);
+ }
}
}
}
@@ -413,7 +427,7 @@
/* set flow accumulation for neighbours */
max_val = -1;
- tci_div = sum_contour = 0.;
+ sca_tanb.tanb = sum_contour = 0.;
if (mfd_cells > 1) {
prop = 0.0;
@@ -427,17 +441,17 @@
if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
- if (tci_flag) {
+ weight[ct_dir] = weight[ct_dir] / sum_weight;
+ /* check everything adds up to 1.0 */
+ prop += weight[ct_dir];
+
+ if (atanb_flag) {
sum_contour += contour[ct_dir];
- tci_div += get_slope_tci(ele, ele_nbr[ct_dir],
+ sca_tanb.tanb += get_slope_tci(ele, ele_nbr[ct_dir],
dist_to_nbr[ct_dir]) *
- weight[ct_dir];
+ weight[ct_dir];
}
- weight[ct_dir] = weight[ct_dir] / sum_weight;
- /* check everything adds up to 1.0 */
- prop += weight[ct_dir];
-
if (value > 0) {
if (wat_nbr[ct_dir] > 0)
wat_nbr[ct_dir] += value * weight[ct_dir];
@@ -473,8 +487,6 @@
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 {
@@ -495,18 +507,18 @@
wa.ele = ele_nbr[np_side];
seg_put(&watalt, (char *)&wa, dr, dc);
- if (tci_flag) {
+ if (atanb_flag) {
sum_contour = contour[np_side];
- tci_div = get_slope_tci(ele, ele_nbr[np_side],
- dist_to_nbr[np_side]);
+ sca_tanb.tanb = 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) * cell_size) /
- (sum_contour * tci_div));
- dseg_put(&tci, &tci_val, r, c);
+ /* topographic wetness index ln(a / tan(beta)) and
+ * stream power index a * tan(beta) */
+ if (atanb_flag) {
+ sca_tanb.sca = fabs(value) * (cell_size / sum_contour);
+ seg_put(&atanb, (char *)&sca_tanb, r, c);
}
}
seg_put(&aspflag, (char *)&af, r, c);
Modified: grass/trunk/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/init_vars.c 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/seg/init_vars.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -22,6 +22,7 @@
DCELL dvalue;
WAT_ALT wa, *wabuf;
ASP_FLAG af, af_nbr, *afbuf;
+ A_TANB sca_tanb;
char MASK_flag;
void *elebuf, *ptr, *watbuf, *watptr;
int ele_map_type, wat_map_type;
@@ -32,7 +33,8 @@
/* input */
ele_flag = pit_flag = run_flag = ril_flag = 0;
/* output */
- wat_flag = asp_flag = bas_flag = seg_flag = haf_flag = tci_flag = 0;
+ wat_flag = asp_flag = tci_flag = spi_flag = atanb_flag = 0;
+ bas_flag = seg_flag = haf_flag = 0;
bas_thres = 0;
/* shed, unused */
arm_flag = dis_flag = 0;
@@ -59,6 +61,8 @@
wat_flag++;
else if (sscanf(argv[r], "tci=%s", tci_name) == 1)
tci_flag++;
+ else if (sscanf(argv[r], "spi=%s", spi_name) == 1)
+ spi_flag++;
else if (sscanf(argv[r], "drainage=%s", asp_name) == 1)
asp_flag++;
else if (sscanf(argv[r], "depression=%s", pit_name) == 1)
@@ -129,6 +133,9 @@
if (seg_flag || bas_flag || haf_flag)
tot_parts++;
+ if (tci_flag || spi_flag)
+ atanb_flag = 1;
+
G_message(n_("SECTION 1 beginning: Initiating Variables. %d section total.",
"SECTION 1 beginning: Initiating Variables. %d sections total.",
tot_parts),
@@ -184,9 +191,9 @@
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);
+ if (atanb_flag) {
+ memory_divisor += sizeof(A_TANB);
+ disk_space += sizeof(A_TANB);
}
/* RUSLE */
if (er_flag) {
@@ -258,8 +265,11 @@
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);
+ if (atanb_flag) {
+ seg_open(&atanb, nrows, ncols, seg_rows, seg_cols, num_open_segs, sizeof(A_TANB));
+ Rast_set_d_null_value(&sca_tanb.sca, 1);
+ Rast_set_d_null_value(&sca_tanb.tanb, 1);
+ }
/* open elevation input */
ele_fd = Rast_open_old(ele_name, "");
@@ -355,6 +365,9 @@
wabuf[c].wat = wat_value;
wabuf[c].ele = alt_value;
alt_value_buf[c] = alt_value;
+ if (atanb_flag) {
+ seg_put(&atanb, (char *)&sca_tanb, r, c);
+ }
ptr = G_incr_void_ptr(ptr, ele_size);
if (run_flag) {
watptr = G_incr_void_ptr(watptr, wat_size);
Modified: grass/trunk/raster/r.watershed/seg/main.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/main.c 2015-12-21 10:06:25 UTC (rev 67297)
+++ grass/trunk/raster/r.watershed/seg/main.c 2015-12-21 10:06:57 UTC (rev 67298)
@@ -39,7 +39,8 @@
BSEG s_b;
CSEG dis, alt, bas, haf, r_h, dep;
SSEG watalt, aspflag;
-DSEG tci, slp, s_l, s_g, l_s, ril;
+DSEG slp, s_l, s_g, l_s, ril;
+SSEG atanb;
double segs_mb;
char zero, one;
double ril_value, d_zero, d_one;
@@ -56,11 +57,12 @@
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], tci_name[GNAME_MAX];
+char wat_name[GNAME_MAX], asp_name[GNAME_MAX];
+char tci_name[GNAME_MAX], spi_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, tci_flag;
+char bas_flag, seg_flag, haf_flag, er_flag, tci_flag, spi_flag, atanb_flag;
char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
FILE *fp;
More information about the grass-commit
mailing list