r.watershed: +SPI

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@67298 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.watershed/front/main.c

@@ -58,6 +58,7 @@ int main(int argc, char *argv[])
     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 @@ int main(int argc, char *argv[])
     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 @@ int main(int argc, char *argv[])
     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 @@ int main(int argc, char *argv[])
 	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)",

raster/r.watershed/front/r.watershed.html

@@ -113,17 +113,29 @@ geographic region. Thus, any cells with negative values cannot have
 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(&alpha; / tan(&beta;))</tt> where &alpha; is the cumulative
 upslope area draining through a point per unit contour length and
 <tt>tan(&beta;)</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>&alpha; /
 tan(&beta;) &lt; 1</tt>.
 
 <p>
+Output <b>spi</b> raster map contains stream power index SPI,
+computed as
+<tt>&alpha; * tan(&beta;)</tt> where &alpha; is the cumulative
+upslope area draining through a point per unit contour length and
+<tt>tan(&beta;)</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>&alpha; &lt; 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 &quot;aspect&quot; for each cell measured CCW from East.
 Multiplying positive values by 45 will give the direction in degrees
@@ -486,6 +498,12 @@ drainage networks from massive, radar-based elevation models with least
 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), 

raster/r.watershed/ram/Gwater.h

@@ -60,7 +60,7 @@ extern RAMSEG r_h_seg, dep_seg;
 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 char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
 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;

raster/r.watershed/ram/close_maps.c

@@ -131,7 +131,6 @@ int close_maps(void)
 
     /* TCI */
     if (tci_flag) {
-	DCELL watvalue;
 	double mean;
 
 	sum = sum_sqr = stddev = 0.0;
@@ -139,9 +138,9 @@ int close_maps(void)
 	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 @@ int close_maps(void)
 	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) {

raster/r.watershed/ram/close_maps2.c

@@ -71,7 +71,7 @@ int close_array_seg(void)
 	    }
 	    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 ... */

raster/r.watershed/ram/do_cum.c

@@ -109,7 +109,7 @@ int do_cum(void)
     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 @@ int do_cum(void)
 		    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 @@ int do_cum(void)
 	    }
 	    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 @@ int do_cum_mfd(void)
 			    if (dr == r_nbr && dc == c_nbr) {
 				astar_not_set = 0;
 			    }
+			    if (value < 0 && valued > 0)
+				wat[nbr_index] = -valued;
 			}
 		    }
 		}
@@ -401,17 +407,17 @@ int do_cum_mfd(void)
 
 			    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 @@ int do_cum_mfd(void)
 		    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 @@ int do_cum_mfd(void)
 		}
 		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;
 	    }
 	}
     }

raster/r.watershed/ram/init_vars.c

@@ -22,7 +22,8 @@ int init_vars(int argc, char *argv[])
     /* 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 @@ int init_vars(int argc, char *argv[])
 	    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 @@ int init_vars(int argc, char *argv[])
     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 @@ int init_vars(int argc, char *argv[])
 	    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);
 	}
     }

raster/r.watershed/ram/main.c

@@ -39,7 +39,7 @@ RAMSEG r_h_seg, dep_seg;
 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 @@ char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX],
     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;
 

raster/r.watershed/seg/Gwater.h

@@ -63,6 +63,12 @@ WAT_ALT {
    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 SSEG astar_pts;
 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 char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
 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;

raster/r.watershed/seg/close_maps.c

@@ -124,89 +124,175 @@ int close_maps(void)
 	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;
+
+	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"));
 
-	G_message(_("Closing TCI 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);
+
+	    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);
+	    }
 
-	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);
+
+	    /* 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 */
 
-	if (max > 0 && max > clr_max) {
+	    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);

raster/r.watershed/seg/do_cum.c

@@ -112,8 +112,9 @@ int do_cum(void)
     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 @@ int do_cum(void)
 		    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 @@ int do_cum(void)
 	    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 do_cum_mfd(void)
 {
     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 @@ int do_cum_mfd(void)
 			    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 @@ int do_cum_mfd(void)
 
 	    /* 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 @@ int do_cum_mfd(void)
 
 			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];
 
+			    if (atanb_flag) {
+				sum_contour += contour[ct_dir];
+				sca_tanb.tanb += get_slope_tci(ele, ele_nbr[ct_dir],
+				                         dist_to_nbr[ct_dir]) *
+					         weight[ct_dir];
+			    }
+
 			    if (value > 0) {
 				if (wat_nbr[ct_dir] > 0)
 				    wat_nbr[ct_dir] += value * weight[ct_dir];
@@ -473,8 +487,6 @@ int do_cum_mfd(void)
 		    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 @@ int do_cum_mfd(void)
 		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);

raster/r.watershed/seg/init_vars.c

@@ -22,6 +22,7 @@ int init_vars(int argc, char *argv[])
     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 @@ int init_vars(int argc, char *argv[])
     /* 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 @@ int init_vars(int argc, char *argv[])
 	    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 @@ int init_vars(int argc, char *argv[])
     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 @@ int init_vars(int argc, char *argv[])
     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 @@ int init_vars(int argc, char *argv[])
     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 @@ int init_vars(int argc, char *argv[])
 	    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);

raster/r.watershed/seg/main.c

@@ -39,7 +39,8 @@ SSEG astar_pts;
 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 @@ char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX],
     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;