[GRASS-SVN] r35412 - in grass/trunk/raster/r.watershed: front ram
seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jan 15 04:15:17 EST 2009
Author: mmetz
Date: 2009-01-15 04:15:17 -0500 (Thu, 15 Jan 2009)
New Revision: 35412
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_astar.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/sg_factor.c
grass/trunk/raster/r.watershed/ram/slope_len.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_astar.c
grass/trunk/raster/r.watershed/seg/do_cum.c
grass/trunk/raster/r.watershed/seg/dseg_get.c
grass/trunk/raster/r.watershed/seg/dseg_put.c
grass/trunk/raster/r.watershed/seg/dseg_read.c
grass/trunk/raster/r.watershed/seg/dseg_write.c
grass/trunk/raster/r.watershed/seg/init_vars.c
grass/trunk/raster/r.watershed/seg/main.c
grass/trunk/raster/r.watershed/seg/sg_factor.c
grass/trunk/raster/r.watershed/seg/slope_len.c
Log:
MFD stable, several enhancements
Modified: grass/trunk/raster/r.watershed/front/main.c
===================================================================
--- grass/trunk/raster/r.watershed/front/main.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/front/main.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -20,7 +20,7 @@
#include <grass/gis.h>
#include <grass/glocale.h>
-int write_hist(char *, char *, char *, int);
+int write_hist(char *, char *, char *, int, int);
int main(int argc, char *argv[])
{
@@ -62,7 +62,8 @@
opt2 = G_define_option();
opt2->key = "depression";
- opt2->description = _("Input map: locations of real depressions");
+ opt2->label = _("Input map: locations of real depressions");
+ opt2->description = _("All non-NULL cells are considered as real depressions");
opt2->required = NO;
opt2->type = TYPE_STRING;
opt2->gisprompt = "old,cell,raster";
@@ -87,8 +88,10 @@
opt5 = G_define_option();
opt5->key = "blocking";
+ opt5->label =
+ _("Input map: terrain blocking overland surface flow, for USLE");
opt5->description =
- _("Input map: terrain blocking overland surface flow, for USLE");
+ _("All non-NULL cells are considered as blocking terrain");
opt5->required = NO;
opt5->type = TYPE_STRING;
opt5->gisprompt = "old,cell,raster";
@@ -232,18 +235,20 @@
}
err = 0;
- /* basin and basin.thresh */
+ /* basin and basin threshold */
err += (opt10->answer != NULL && opt6->answer == NULL);
- /* stream and basin.thresh */
+ /* stream and basin threshold */
err += (opt11->answer != NULL && opt6->answer == NULL);
- /* half.basin and basin.thresh */
+ /* half_basin and basin threshold */
err += (opt12->answer != NULL && opt6->answer == NULL);
- /* slope and basin.thresh */
+ /* LS factor and basin threshold */
+ err += (opt14->answer != NULL && opt6->answer == NULL);
+ /* S factor and basin threshold */
err += (opt15->answer != NULL && opt6->answer == NULL);
if (err) {
G_message(_("Sorry, if any of the following options are set:\n"
- " basin, stream, half.basin, slope, or lS\n"
+ " basin, stream, half_basin, length_slope, or slope_steepness\n"
" you MUST provide a value for the basin "
"threshold parameter."));
G_usage();
@@ -388,39 +393,41 @@
if (opt8->answer)
write_hist(opt8->answer,
"Watershed accumulation: overland flow that traverses each cell",
- opt1->answer, flag_seg->answer);
+ opt1->answer, flag_seg->answer, flag_sfd->answer);
if (opt9->answer)
write_hist(opt9->answer,
"Watershed drainage direction (divided by 45deg)",
- opt1->answer, flag_seg->answer);
+ opt1->answer, flag_seg->answer, flag_sfd->answer);
if (opt10->answer)
write_hist(opt10->answer,
- "Watershed basins", opt1->answer, flag_seg->answer);
+ "Watershed basins", opt1->answer, flag_seg->answer,
+ flag_sfd->answer);
if (opt11->answer)
write_hist(opt11->answer,
"Watershed stream segments", opt1->answer,
- flag_seg->answer);
+ flag_seg->answer, flag_sfd->answer);
if (opt12->answer)
write_hist(opt12->answer,
- "Watershed half-basins", opt1->answer, flag_seg->answer);
+ "Watershed half-basins", opt1->answer, flag_seg->answer,
+ flag_sfd->answer);
if (opt13->answer)
write_hist(opt13->answer,
"Watershed visualization map (filtered accumulation map)",
- opt1->answer, flag_seg->answer);
+ opt1->answer, flag_seg->answer, flag_sfd->answer);
if (opt14->answer)
write_hist(opt14->answer,
"Watershed slope length and steepness (LS) factor",
- opt1->answer, flag_seg->answer);
+ opt1->answer, flag_seg->answer, flag_sfd->answer);
if (opt15->answer)
write_hist(opt15->answer,
"Watershed slope steepness (S) factor",
- opt1->answer, flag_seg->answer);
+ opt1->answer, flag_seg->answer, flag_sfd->answer);
exit(ret);
}
/* record map history info */
-int write_hist(char *map_name, char *title, char *source_name, int mode)
+int write_hist(char *map_name, char *title, char *source_name, int mode, int sfd)
{
struct History history;
@@ -430,8 +437,10 @@
strncpy(history.datsrc_1, source_name, RECORD_LEN);
history.datsrc_1[RECORD_LEN - 1] = '\0'; /* strncpy() doesn't null terminate if maxfill */
sprintf(history.edhist[0],
- "Processing mode: %s", mode ? "Segmented" : "All in RAM");
- history.edlinecnt = 1;
+ "Processing mode: %s", sfd ? "SFD (D8)" : "MFD");
+ sprintf(history.edhist[1],
+ "Memory mode: %s", mode ? "Segmented" : "All in RAM");
+ history.edlinecnt = 2;
G_command_history(&history);
return G_write_history(map_name, &history);
Modified: grass/trunk/raster/r.watershed/front/r.watershed.html
===================================================================
--- grass/trunk/raster/r.watershed/front/r.watershed.html 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/front/r.watershed.html 2009-01-15 09:15:17 UTC (rev 35412)
@@ -1,7 +1,7 @@
<h2>DESCRIPTION</h2>
<em>r.watershed</em> generates a set of maps indicating:
-1) the location of watershed basins, and
+1) flow accumulation, drainage direction, the location of streams and watershed basins, and
2) the LS and S factors of the Revised Universal Soil Loss Equation (RUSLE).
<p>
@@ -24,13 +24,13 @@
<dd>Without this flag set, the entire analysis is run in memory
maintained by the operating system. This can be limiting, but is
-relatively fast. Setting the flag causes the program to manage memory
+very fast. Setting the flag causes the program to manage memory
on disk which allows larger maps to be processed but is considerably
slower.
<dt><em>-s</em>
-<dd>Use single flow direction (SFD) instead of multiple flow direction (MFD).
+<dd>Use single flow direction (SFD, D8) instead of multiple flow direction (MFD).
MFD is enabled by default.
<dt><em>-4</em>
@@ -59,7 +59,8 @@
<dd>Input map: Map layer of actual depressions or sinkholes in the
landscape that are large enough to slow and store surface runoff from
-a storm event. Any non-zero values indicate depressions.
+a storm event. Any non-NULL (not nodata) values indicate depressions.
+Water will flow into depressions, but not out of depressions.
<dt><em>flow</em>
@@ -67,7 +68,7 @@
amount of overland flow units that each cell will contribute to the
watershed basin model. Overland flow units represent the amount of
overland flow each cell contributes to surface flow. If omitted, a
-value of one (1) is assumed. The algorithm is D8 flow accumulation.
+value of one (1) is assumed.
<dt><em>disturbed_land</em>
@@ -80,7 +81,7 @@
<dd>Input map: terrain that will block overland surface flow. Terrain
that will block overland surface flow and restart the slope length
-for the RUSLE. Any non-zero values indicate blocking terrain.
+for the RUSLE. Any non-NULL (not nodata) values indicate blocking terrain.
<dt><em>threshold</em>
@@ -116,8 +117,8 @@
<dd>Output map: drainage direction. Provides the "aspect" for each
cell. Multiplying positive values by 45 will give the direction in
degrees that the surface runoff will travel from that cell. The
-value -1 indicates that the cell is a depression area (defined by
-the depression input map). Other negative values indicate that
+value 0 (zero) indicates that the cell is a depression area (defined by
+the depression input map). Negative values indicate that
surface runoff is leaving the boundaries of the current geographic
region. The absolute value of these negative cells indicates the
direction of flow.
@@ -154,22 +155,19 @@
<dt><em>length_slope</em>
-<dd>Output map: slope length and steepness (LS) factor. Contains the LS
-factor for the Revised Universal Soil Loss Equation. Equations taken
-from <em>Revised Universal Soil Loss Equation for Western Rangelands</em>
-(Weltz et al. 1987).
-Since the LS factor is a small number, it is multiplied by 100 for the
-GRASS output map.
+<dd>Output map: slope length and steepness (LS) factor for the Revised
+Universal Soil Loss Equation (RUSLE). Equations taken from <em>Revised
+Universal Soil Loss Equation for Western Rangelands</em>
+(Weltz et al. 1987). Since the LS factor is a small number (usually less
+than one), the GRASS output map is of type DCELL.
<dt><em>slope_steepness</em>
-<dd>Output map: slope steepness (S) factor for RUSLE.
-Contains the revised S factor for the Universal Soil
-Loss Equation. Equations taken from article entitled
+<dd>Output map: slope steepness (S) factor for the Universal Soil
+Loss Equation (RUSLE). Equations taken from article entitled
<em>Revised Slope Steepness Factor for the Universal Soil
-Loss Equation</em> (McCool et al. 1987). Since the S factor
-is a small number (usually less than one), it is multiplied
-by 100 for the GRASS output map layer.
+Loss Equation</em> (McCool et al. 1987). Since the LS factor is a small
+number (usually less than one), the GRASS output map is of type DCELL.
</dd>
</dl>
@@ -208,13 +206,13 @@
MFD, water flow is distributed to all neighbouring cells with lower
elevation, using slope towards neighbouring cells as a weighing factor
for proportional distribution. The A<sup>T</sup> least-cost path is
-always included and assigned the maxmimum observed weighing factor.
-As a result, depressions and obstacles are overflown with a gracefull
-flow convergence before the overflow. The convergence factor causes flow
-accumulation to converge more strongly with higher values. The supported
-range is 1 to 10, recommended is a convergence factor of 5. If many
-small sliver basins are created with MFD, setting the convergence factor
-to a higher value can reduce the amount of small sliver basins.
+always included. As a result, depressions and obstacles are overflown
+with a gracefull flow convergence before the overflow. The convergence
+factor causes flow accumulation to converge more strongly with higher
+values. The supported range is 1 to 10, recommended is a convergence
+factor of 5 (Holmgren, 1994). If many small sliver basins are created
+with MFD, setting the convergence factor to a higher value can reduce
+the amount of small sliver basins.
<br>See example below with the North Carolina dataset for using MFD mode.
<h4>In-memory mode and disk swap mode</h4>
@@ -222,12 +220,17 @@
<em>ram</em> is used by default, <em>seg</em> can be used by setting
the <em>-m</em> flag.
<br>
+The <em>ram</em> version requires a maximum of 31 MB of RAM for 1 million
+cells. Together with the amount of system memory (RAM) available, this
+value can be used to estimate whether the current region can be
+processed with the <em>ram</em> version.
+<br>
The <em>ram</em> version uses virtual memory managed by the operating
system to store all the data structures and is faster than the <em>seg</em>
version; <em>seg</em> uses the GRASS segmentation library which manages
-data in disk files. Thus <em>seg</em> uses much less system memory (RAM)
-allowing other processes to operate on the same CPU, even when the
-current geographic region is huge.
+data in disk files. <em>seg</em> uses only as much system memory (RAM) as
+specified with the <em>memory</em> option, allowing other processes to
+operate on the same system, even when the current geographic region is huge.
<br>
Due to memory requirements of both programs, it is quite easy to run out of
memory when working with huge map regions. If the <em>ram</em> version runs
@@ -280,15 +283,16 @@
<p>
The <em>r.watershed</em> program does not require the user to have the
current geographic region filled with elevation values. Areas without
-elevation data MUST be masked out, by creating a raster map (or raster
-reclassification) named <tt>MASK</tt>. Areas
-masked out will be treated as if they are off the edge of the region.
-MASKs will reduce the memory necessary to run the program. Masking out
-unimportant areas can significantly reduce processing time if the watersheds
-of interest occupy a small percentage of the overall area.
+elevation data (masked or NULL cells) are ignored. It is NOT necessary to
+create a raster map (or raster reclassification) named <tt>MASK</tt> for
+NULL cells. Areas without elevation data will be treated as if they are
+off the edge of the region. Such areas will reduce the memory necessary
+to run the program. Masking out unimportant areas can significantly
+reduce processing time if the watersheds of interest occupy a small
+percentage of the overall area.
<p>
-Zero data values will be treated as elevation data (not no_data).
+Zero (0) and negative data values will be treated as elevation data (not no_data).
<h4>Further processing of output layers</h4>
<p>
Modified: grass/trunk/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/Gwater.h 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/ram/Gwater.h 2009-01-15 09:15:17 UTC (rev 35412)
@@ -34,7 +34,7 @@
#define POINT struct points
POINT {
- SHORT r, c, downr, downc;
+ SHORT r, c; /* , downr, downc */
int nxt;
};
@@ -87,7 +87,7 @@
/* do_astar.c */
int do_astar(void);
-int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int add_pt(SHORT, SHORT, CELL, CELL);
int drop_pt(void);
int sift_up(int, CELL);
double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
Modified: grass/trunk/raster/r.watershed/ram/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/close_maps.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/ram/close_maps.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -6,21 +6,24 @@
int close_maps(void)
{
- struct Colors colors, logcolors;
+ struct Colors colors;
int value, r, c, fd;
CELL *buf = NULL;
DCELL *dbuf = NULL;
struct FPRange accRange;
DCELL min, max;
DCELL clr_min, clr_max;
+ DCELL sum, sum_sqr, stddev, lstddev, dvalue;
- if (wat_flag || asp_flag || dis_flag || ls_flag || sl_flag || sg_flag)
+ if (asp_flag || dis_flag)
buf = G_allocate_cell_buf();
- dbuf = G_allocate_d_raster_buf();
+ if (wat_flag || ls_flag || sl_flag || sg_flag)
+ dbuf = G_allocate_d_raster_buf();
G_free(alt);
if (ls_flag || sg_flag)
G_free(r_h);
+ sum = sum_sqr = stddev = 0.0;
if (wat_flag) {
fd = G_open_raster_new(wat_name, DCELL_TYPE);
if (fd < 0) {
@@ -28,65 +31,92 @@
}
else {
for (r = 0; r < nrows; r++) {
+ G_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
for (c = 0; c < ncols; c++) {
- dbuf[c] = wat[SEG_INDEX(wat_seg, r, c)];
+ dvalue = wat[SEG_INDEX(wat_seg, r, c)];
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dbuf[c] = dvalue;
+ dvalue = ABS(dvalue);
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
}
G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
if (G_close_cell(fd) < 0)
G_warning(_("Close failed."));
+ 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);
+
G_read_fp_range(wat_name, this_mapset, &accRange);
min = max = 0;
G_get_fp_range_min_max(&accRange, &min, &max);
G_init_colors(&colors);
- if (ABS(min) > max) {
- clr_min = 1;
- clr_max = 0.25 * max;
- G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
- 255, 0, &colors);
- clr_min = 0.25 * max;
- clr_max = 0.5 * max;
- G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0,
+
+ if (min < 0) {
+ if (min < (-stddev - 1)) {
+ clr_min = min;
+ clr_max = -stddev - 1;
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
+ 0, 0, &colors);
+ }
+ clr_min = -stddev - 1.;
+ clr_max = -1. * exp(lstddev * 0.75);
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
+ 0, 255, &colors);
+ clr_min = clr_max;
+ clr_max = -1. * exp(lstddev * 0.5);
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0,
255, 255, &colors);
- clr_min = 0.5 * max;
- clr_max = 0.75 * max;
+ clr_min = clr_max;
+ clr_max = -1. * exp(lstddev * 0.35);
G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
- 0, 255, &colors);
- clr_min = 0.75 * max;
+ 255, 0, &colors);
+ clr_min = clr_max;
+ clr_max = -1.;
+ G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 255,
+ 255, 0, &colors);
+ }
+ clr_min = -1.;
+ clr_max = 1.;
+ G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 255,
+ 255, 0, &colors);
+ clr_min = 1;
+ clr_max = exp(lstddev * 0.35);
+ G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
+ 255, 0, &colors);
+ clr_min = clr_max;
+ clr_max = exp(lstddev * 0.5);
+ G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0,
+ 255, 255, &colors);
+ clr_min = clr_max;
+ clr_max = exp(lstddev * 0.75);
+ G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
+ 0, 255, &colors);
+ clr_min = clr_max;
+ clr_max = stddev + 1.;
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
+ 0, &colors);
+
+ if (max > 0 && max > stddev + 1) {
+ clr_min = stddev + 1;
clr_max = max;
- G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0, 0,
0, &colors);
- max = -max;
}
- else {
- min = ABS(min);
- clr_min = 1;
- clr_max = 0.25 * min;
- G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
- 255, 0, &colors);
- clr_min = 0.25 * min;
- clr_max = 0.5 * min;
- G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0,
- 255, 255, &colors);
- clr_min = 0.5 * min;
- clr_max = 0.75 * min;
- G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
- 0, 255, &colors);
- clr_min = 0.75 * min;
- clr_max = min;
- G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
- 0, &colors);
- }
- G_abs_log_colors(&logcolors, &colors, 5);
- G_add_d_raster_color_rule(&min, 0, 0, 0, &max, 0, 0, 0,
- &logcolors);
- G_write_colors(wat_name, this_mapset, &logcolors);
+ G_write_colors(wat_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) {
fd = G_open_cell_new(asp_name);
if (fd < 0) {
@@ -94,6 +124,7 @@
}
else {
for (r = 0; r < nrows; r++) {
+ G_set_c_null_value(buf, ncols); /* reset row to all NULL */
for (c = 0; c < ncols; c++) {
buf[c] = asp[SEG_INDEX(asp_seg, r, c)];
}
@@ -146,16 +177,16 @@
G_free(wat);
if (ls_flag) {
- fd = G_open_cell_new(ls_name);
+ fd = G_open_raster_new(ls_name, DCELL_TYPE);
if (fd < 0) {
- G_warning(_("unable to open new L map layer."));
+ G_warning(_("unable to open new LS factor map layer."));
}
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++) {
- buf[c] = l_s[SEG_INDEX(l_s_seg, r, c)] + .5;
+ dbuf[c] = l_s[SEG_INDEX(l_s_seg, r, c)];
}
- G_put_raster_row(fd, buf, CELL_TYPE);
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
if (G_close_cell(fd) < 0)
G_warning(_("Close failed."));
@@ -164,18 +195,18 @@
}
if (sl_flag) {
- fd = G_open_cell_new(sl_name);
+ fd = G_open_raster_new(sl_name, DCELL_TYPE);
if (fd < 0) {
G_warning(_("unable to open new slope length map layer."));
}
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++) {
- buf[c] = s_l[SEG_INDEX(s_l_seg, r, c)] + .5;
- if (buf[c] > max_length)
- buf[c] = max_length + .5;
+ dbuf[c] = s_l[SEG_INDEX(s_l_seg, r, c)];
+ if (dbuf[c] > max_length)
+ dbuf[c] = max_length;
}
- G_put_raster_row(fd, buf, CELL_TYPE);
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
if (G_close_cell(fd) < 0)
G_warning(_("Close failed."));
@@ -186,16 +217,16 @@
G_free(s_l);
if (sg_flag) {
- fd = G_open_cell_new(sg_name);
+ fd = G_open_raster_new(sg_name, DCELL_TYPE);
if (fd < 0) {
- G_warning(_("unable to open new S map layer."));
+ G_warning(_("unable to open new S factor map layer."));
}
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++) {
- buf[c] = s_g[SEG_INDEX(s_g_seg, r, c)] * 100 + .5;
+ dbuf[c] = s_g[SEG_INDEX(s_g_seg, r, c)];
}
- G_put_raster_row(fd, buf, CELL_TYPE);
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
if (G_close_cell(fd) < 0)
G_warning(_("Close failed."));
Modified: grass/trunk/raster/r.watershed/ram/close_maps2.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/close_maps2.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/ram/close_maps2.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -112,7 +112,7 @@
G_write_colors(bas_name, this_mapset, &colors);
}
- /* half.basins map */
+ /* half_basins map */
if (haf_flag) {
map_fd = G_open_cell_new(haf_name);
for (r = 0; r < nrows; r++) {
Modified: grass/trunk/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/ram/do_astar.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -36,7 +36,7 @@
/* equivalent to first_astar = point->next in old code */
drop_pt();
- /* can go, dragged on from old code */
+ /* can go, dragged on from old code: replace first_astar with heap_index[1] in line 22 */
first_astar = heap_index[1];
/* downhill path for flow accumulation is set here */
@@ -51,7 +51,6 @@
index_doer = SEG_INDEX(alt_seg, r, c);
alt_val = alt[index_doer];
- wat_val = wat[index_doer];
FLAG_SET(worked, r, c);
@@ -69,24 +68,24 @@
if (in_val == 0) {
alt_up = alt[index_up];
/* flow direction is set here */
- add_pt(upr, upc, r, c, alt_up, alt_val);
+ add_pt(upr, upc, alt_up, alt_val);
drain_val = drain[upr - r + 1][upc - c + 1];
asp[index_up] = drain_val;
-
}
else {
- /* check if neighbour has not been worked on,
- * update values for asp and wat */
+ /* check if neighbour has been worked on,
+ * if not, update values for asp and wat */
in_val = FLAG_GET(worked, upr, upc);
if (in_val == 0) {
asp_up = asp[index_up];
- if (asp_up < -1) {
+ if (asp_up < 0) {
asp[index_up] = drain[upr - r + 1][upc - c + 1];
+ wat_val = wat[index_doer];
if (wat_val > 0)
wat[index_doer] = -wat_val;
- replace(upr, upc, r, c); /* alt_up used to be */
+ /* replace(upr, upc, r, c); */ /* alt_up used to be */
}
}
}
@@ -104,7 +103,7 @@
}
/* new add point routine for min heap */
-int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
{
FLAG_SET(in_list, r, c);
@@ -120,8 +119,8 @@
astar_pts[nxt_avail_pt].r = r;
astar_pts[nxt_avail_pt].c = c;
- astar_pts[nxt_avail_pt].downr = downr;
- astar_pts[nxt_avail_pt].downc = downc;
+/* astar_pts[nxt_avail_pt].downr = downr;
+ astar_pts[nxt_avail_pt].downc = downc; */
nxt_avail_pt++;
@@ -288,8 +287,8 @@
while (heap_run <= heap_size) {
now = heap_index[heap_run];
if (astar_pts[now].r == upr && astar_pts[now].c == upc) {
- astar_pts[now].downr = r;
- astar_pts[now].downc = c;
+ /*astar_pts[now].downr = r;
+ astar_pts[now].downc = c; */
return 0;
}
heap_run++;
Modified: grass/trunk/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_cum.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/ram/do_cum.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -6,8 +6,10 @@
int do_cum(void)
{
SHORT r, c, dr, dc;
- CELL is_swale, value, valued;
+ CELL is_swale, value, valued, aspect;
int killer, threshold, count;
+ SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
@@ -20,12 +22,18 @@
G_percent(count++, do_points, 2);
killer = first_cum;
first_cum = astar_pts[killer].nxt;
- if ((dr = astar_pts[killer].downr) > -1) {
- r = astar_pts[killer].r;
- c = astar_pts[killer].c;
- dc = astar_pts[killer].downc;
+ r = astar_pts[killer].r;
+ c = astar_pts[killer].c;
+ aspect = asp[SEG_INDEX(asp_seg, r, c)];
+ if (aspect) {
+ dr = r + asp_r[ABS(aspect)];
+ dc = c + asp_c[ABS(aspect)];
+ }
+ else
+ dr = dc = -1;
+ if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
value = wat[SEG_INDEX(wat_seg, r, c)];
- if (ABS(value) >= threshold)
+ if ((int)(ABS(value) + 0.5) >= threshold)
FLAG_SET(swale, r, c);
valued = wat[SEG_INDEX(wat_seg, dr, dc)];
if (value > 0) {
@@ -43,6 +51,13 @@
wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
valued = ABS(valued) + 0.5;
is_swale = FLAG_GET(swale, r, c);
+ /* update asp for depression */
+ if (is_swale && pit_flag) {
+ if (aspect > 0 && asp[SEG_INDEX(asp_seg, dr, dc)] == 0) {
+ aspect = -aspect;
+ asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ }
+ }
if (is_swale || ((int)valued) >= threshold) {
FLAG_SET(swale, dr, dc);
}
@@ -86,17 +101,19 @@
CELL is_swale;
DCELL value, valued;
int killer, threshold, count;
+ SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
/* MFD */
- int mfd_cells, stream_cells, swale_cells, astar_not_set;
+ 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, r_swale, c_swale, ct_dir, np_side;
+ 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, max_swale;
+ double prop, max_acc;
int workedon;
- G_message(_("SECTION 3: Accumulating Surface Flow with MFD"));
+ G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
G_debug(1, "MFD convergence factor set to %d.", c_fac);
/* distances to neighbours */
@@ -129,22 +146,32 @@
G_percent(count++, do_points, 2);
killer = first_cum;
first_cum = astar_pts[killer].nxt;
- if ((dr = astar_pts[killer].downr) > -1) {
- r = astar_pts[killer].r;
- c = astar_pts[killer].c;
- dc = astar_pts[killer].downc;
+ r = astar_pts[killer].r;
+ c = astar_pts[killer].c;
+ aspect = asp[SEG_INDEX(asp_seg, r, c)];
+ if (aspect) {
+ dr = r + asp_r[ABS(aspect)];
+ dc = c + asp_c[ABS(aspect)];
+ }
+ else
+ dr = dc = -1;
+ if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
value = wat[SEG_INDEX(wat_seg, r, c)];
valued = wat[SEG_INDEX(wat_seg, dr, dc)];
- /* start MFD */
+ r_max = dr;
+ c_max = dc;
/* get weights */
max_weight = 0;
sum_weight = 0;
np_side = -1;
mfd_cells = 0;
+ stream_cells = 0;
+ swale_cells = 0;
astar_not_set = 1;
ele = alt[SEG_INDEX(alt_seg, r, c)];
+ is_null = 0;
/* this loop is needed to get the sum of weights */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (r_nbr, c_nbr) for neighbours */
@@ -154,19 +181,37 @@
/* check that neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols) {
+
+ /* check for swale or stream cells */
+ is_swale = FLAG_GET(swale, r_nbr, c_nbr);
+ if (is_swale)
+ swale_cells++;
+ valued = wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)];
+ if ((ABS(valued) + 0.5) >= threshold)
+ stream_cells++;
+
is_worked = FLAG_GET(worked, r_nbr, c_nbr);
if (is_worked == 0) {
ele_nbr = alt[SEG_INDEX(alt_seg, r_nbr, c_nbr)];
- if (ele_nbr < ele) {
- weight[ct_dir] =
- mfd_pow(((ele -
- ele_nbr) / dist_to_nbr[ct_dir]),
- c_fac);
+ is_null = G_is_c_null_value(&ele_nbr);
+ if (!is_null && ele_nbr <= ele) {
+ if (ele_nbr < ele) {
+ weight[ct_dir] =
+ mfd_pow(((ele -
+ ele_nbr) / dist_to_nbr[ct_dir]),
+ c_fac);
+ }
+ if (ele_nbr == ele) {
+ weight[ct_dir] =
+ mfd_pow((0.5 / dist_to_nbr[ct_dir]),
+ c_fac);
+ }
sum_weight += weight[ct_dir];
mfd_cells++;
- if (weight[ct_dir] > max_weight)
+ if (weight[ct_dir] > max_weight) {
max_weight = weight[ct_dir];
+ }
if (dr == r_nbr && dc == c_nbr) {
astar_not_set = 0;
@@ -181,7 +226,6 @@
/* honour A * path
* mfd_cells == 0: fine, SFD along A * path
* mfd_cells == 1 && astar_not_set == 0: fine, SFD along A * path
- * mfd_cells > 1 && astar_not_set == 0: MFD, set A * path to max_weight
* mfd_cells > 0 && astar_not_set == 1: A * path not included, add to mfd_cells
*/
@@ -192,21 +236,8 @@
weight[np_side] = max_weight;
}
- /* MFD, A * path included, set A * path to max_weight */
- if (mfd_cells > 1 && astar_not_set == 0) {
- sum_weight = sum_weight - weight[np_side] + max_weight;
- weight[np_side] = max_weight;
- }
-
/* set flow accumulation for neighbours */
- stream_cells = 0;
- swale_cells = 0;
max_acc = -1;
- max_swale = -1;
- r_max = dr;
- c_max = dc;
- r_swale = dr;
- c_swale = dc;
if (mfd_cells > 1) {
prop = 0.0;
@@ -219,11 +250,11 @@
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols && weight[ct_dir] > -0.5) {
is_worked = FLAG_GET(worked, r_nbr, c_nbr);
- is_swale = FLAG_GET(swale, r_nbr, c_nbr);
if (is_worked == 0) {
weight[ct_dir] = weight[ct_dir] / sum_weight;
- prop += weight[ct_dir]; /* check everything sums up to 1.0 */
+ /* check everything sums up to 1.0 */
+ prop += weight[ct_dir];
valued = wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)];
if (value > 0) {
@@ -240,39 +271,27 @@
}
wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)] = valued;
- if ((ABS(valued) + 0.5) >= threshold)
- stream_cells++;
- if (ABS(valued) > max_acc) {
+ /* get main drainage direction */
+ if (ABS(valued) >= max_acc) {
max_acc = ABS(valued);
r_max = r_nbr;
c_max = c_nbr;
}
- /* assist in stream merging */
- if (is_swale && ABS(valued) > max_swale) {
- max_swale = ABS(valued);
- r_swale = r_nbr;
- c_swale = c_nbr;
- }
}
else if (ct_dir == np_side) {
- workedon++; /* check for consistency with A * path */
+ /* check for consistency with A * path */
+ workedon++;
}
- if (is_swale)
- swale_cells++;
}
}
if (ABS(prop - 1.0) > 5E-6f) {
G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
prop);
}
- valued = wat[SEG_INDEX(wat_seg, dr, dc)];
- if (max_swale > -0.5) {
- r_max = r_swale;
- c_max = c_swale;
- }
}
if (mfd_cells < 2) {
+ valued = wat[SEG_INDEX(wat_seg, dr, dc)];
if (value > 0) {
if (valued > 0)
valued += value;
@@ -288,14 +307,27 @@
wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
}
- /* MFD finished */
-
- valued = ABS(valued) + 0.5;
+ /* update asp */
+ if (dr != r_max || dc != c_max) {
+ aspect = drain[r - r_max + 1][c - c_max + 1];
+ if (asp[SEG_INDEX(asp_seg, r, c)] < 0)
+ aspect = -aspect;
+ asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ }
is_swale = FLAG_GET(swale, r, c);
- /* start new stream: only works with A * path */
- if (!is_swale && (int)valued >= threshold && stream_cells < 2 &&
+ /* update asp for depression */
+ if (is_swale && pit_flag) {
+ if (aspect > 0 && asp[SEG_INDEX(asp_seg, r_max, c_max)] == 0) {
+ aspect = -aspect;
+ asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ }
+ }
+ /* start new stream */
+ value = ABS(value) + 0.5;
+ if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
swale_cells < 1) {
- FLAG_SET(swale, dr, dc);
+ FLAG_SET(swale, r, c);
+ is_swale = 1;
}
/* continue stream */
if (is_swale) {
@@ -305,26 +337,21 @@
if (er_flag && !is_swale)
slope_length(r, c, r_max, c_max);
}
- /* update asp */
- if (dr != r_max || dc != c_max) {
- aspect = drain[r - r_max + 1][c - c_max + 1];
- if (asp[SEG_INDEX(asp_seg, r, c)] < 0)
- aspect = -aspect;
- asp[SEG_INDEX(asp_seg, r, c)] = aspect;
-
- }
FLAG_SET(worked, r, c);
}
}
G_percent(count, do_points, 1); /* finish it */
if (workedon)
- G_warning(_("MFD: A * path already processed when distributing flow: %d of %d"),
+ G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
workedon, do_points);
G_free(astar_pts);
flag_destroy(worked);
+ G_free(dist_to_nbr);
+ G_free(weight);
+
return 0;
}
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -8,7 +8,7 @@
int init_vars(int argc, char *argv[])
{
SHORT r, c;
- CELL *buf, alt_value, wat_value, asp_value;
+ CELL *buf, alt_value, wat_value, asp_value, block_value;
int fd, index;
char MASK_flag;
@@ -121,8 +121,10 @@
buf = G_allocate_cell_buf();
alt =
(CELL *) G_malloc(sizeof(CELL) * size_array(&alt_seg, nrows, ncols));
- r_h =
- (CELL *) G_malloc(sizeof(CELL) * size_array(&r_h_seg, nrows, ncols));
+ if (er_flag) {
+ r_h =
+ (CELL *) G_malloc(sizeof(CELL) * size_array(&r_h_seg, nrows, ncols));
+ }
fd = G_open_cell_old(ele_name, "");
if (fd < 0) {
@@ -133,18 +135,31 @@
in_list = flag_create(nrows, ncols);
worked = flag_create(nrows, ncols);
+ MASK_flag = 0;
+ do_points = nrows * ncols;
for (r = 0; r < nrows; r++) {
G_get_c_raster_row(fd, buf, r);
for (c = 0; c < ncols; c++) {
index = SEG_INDEX(alt_seg, r, c);
- alt[index] = r_h[index] = buf[c];
+ alt_value = alt[index] = buf[c];
+ if (er_flag) {
+ r_h[index] = buf[c];
+ }
/* all flags need to be manually set to zero */
- flag_unset(swale, r, c);
- flag_unset(in_list, r, c);
- flag_unset(worked, r, c);
+ flag_unset(swale, r, c);
+ flag_unset(in_list, r, c);
+ flag_unset(worked, r, c);
+ /* check for masked and NULL cells here */
+ if (G_is_c_null_value(&alt_value)) {
+ FLAG_SET(worked, r, c);
+ FLAG_SET(in_list, r, c);
+ do_points--;
+ }
}
}
G_close_cell(fd);
+ if (do_points < nrows * ncols)
+ MASK_flag = 1;
wat =
(DCELL *) G_malloc(sizeof(DCELL) *
size_array(&wat_seg, nrows, ncols));
@@ -157,20 +172,36 @@
for (r = 0; r < nrows; r++) {
G_get_c_raster_row(fd, buf, r);
for (c = 0; c < ncols; c++) {
- wat[SEG_INDEX(wat_seg, r, c)] = buf[c];
+ if (MASK_flag) {
+ index = FLAG_GET(worked, r, c);
+ if (!index)
+ wat[SEG_INDEX(wat_seg, r, c)] = buf[c];
+ else
+ wat[SEG_INDEX(wat_seg, r, c)] = 0.0;
+ }
+ else
+ wat[SEG_INDEX(wat_seg, r, c)] = buf[c];
}
}
G_close_cell(fd);
}
else {
for (r = 0; r < nrows; r++) {
- for (c = 0; c < ncols; c++)
- wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
+ for (c = 0; c < ncols; c++) {
+ if (MASK_flag) {
+ index = FLAG_GET(worked, r, c);
+ if (!index)
+ wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
+ }
+ else
+ wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
+ }
}
}
asp =
(CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
+ /* depression: drainage direction will be set to zero later */
if (pit_flag) {
fd = G_open_cell_old(pit_name, "");
if (fd < 0) {
@@ -179,12 +210,15 @@
for (r = 0; r < nrows; r++) {
G_get_c_raster_row(fd, buf, r);
for (c = 0; c < ncols; c++) {
- asp[SEG_INDEX(asp_seg, r, c)] = buf[c];
+ asp_value = buf[c];
+ if (!G_is_c_null_value(&asp_value))
+ asp[SEG_INDEX(asp_seg, r, c)] = 1;
}
}
G_close_cell(fd);
}
+ /* this is also creating streams... */
if (ob_flag) {
fd = G_open_cell_old(ob_name, "");
if (fd < 0) {
@@ -193,7 +227,8 @@
for (r = 0; r < nrows; r++) {
G_get_c_raster_row(fd, buf, r);
for (c = 0; c < ncols; c++) {
- if (buf[c])
+ block_value = buf[c];
+ if (!G_is_c_null_value(&block_value))
FLAG_SET(swale, r, c);
}
}
@@ -206,32 +241,12 @@
}
}
- MASK_flag = 0;
- do_points = nrows * ncols;
- if (NULL != G_find_file("cell", "MASK", G_mapset())) {
- MASK_flag = 1;
- if ((fd = G_open_cell_old("MASK", G_mapset())) < 0) {
- G_fatal_error(_("unable to open MASK"));
- }
- else {
- for (r = 0; r < nrows; r++) {
- G_get_c_raster_row_nomask(fd, buf, r);
- for (c = 0; c < ncols; c++) {
- if (!buf[c]) {
- FLAG_SET(worked, r, c);
- FLAG_SET(in_list, r, c);
- do_points--;
- }
- }
- }
- G_close_cell(fd);
- }
+ /* RUSLE: LS and/or S factor */
+ if (er_flag) {
+ s_l =
+ (double *)G_malloc(size_array(&s_l_seg, nrows, ncols) *
+ sizeof(double));
}
- s_l =
- (double *)G_malloc(size_array(&s_l_seg, nrows, ncols) *
- sizeof(double));
- /* astar_pts = (POINT *) G_malloc (nrows * ncols * sizeof (POINT)); */
- astar_pts = (POINT *) G_malloc(do_points * sizeof(POINT));
if (sg_flag) {
s_g =
@@ -244,6 +259,8 @@
sizeof(double));
}
+ astar_pts = (POINT *) G_malloc(do_points * sizeof(POINT));
+
/* heap_index will track astar_pts in ternary min-heap */
/* heap_index is one-based */
heap_index = (int *)G_malloc((do_points + 1) * sizeof(int));
@@ -262,14 +279,20 @@
wat[SEG_INDEX(wat_seg, r, c)] = 0;
}
else {
- s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
+ if (er_flag)
+ s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
asp_value = asp[SEG_INDEX(asp_seg, r, c)];
if (r == 0 || c == 0 || r == nrows - 1 ||
c == ncols - 1 || asp_value != 0) {
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
- if (r == 0)
+ /* set depression */
+ if (asp_value) {
+ asp_value = 0;
+ wat[SEG_INDEX(wat_seg, r, c)] = ABS(wat_value);
+ }
+ else if (r == 0)
asp_value = -2;
else if (c == 0)
asp_value = -4;
@@ -277,15 +300,13 @@
asp_value = -6;
else if (c == ncols - 1)
asp_value = -8;
- else
- asp_value = -1;
asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
}
else if (FLAG_GET(worked, r - 1, c)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -2;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -293,7 +314,7 @@
}
else if (FLAG_GET(worked, r + 1, c)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -6;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -301,7 +322,7 @@
}
else if (FLAG_GET(worked, r, c - 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -4;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -309,7 +330,7 @@
}
else if (FLAG_GET(worked, r, c + 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -8;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -317,7 +338,7 @@
}
else if (sides == 8 && FLAG_GET(worked, r - 1, c - 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -3;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -325,7 +346,7 @@
}
else if (sides == 8 && FLAG_GET(worked, r - 1, c + 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -1;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -333,7 +354,7 @@
}
else if (sides == 8 && FLAG_GET(worked, r + 1, c - 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -5;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -341,7 +362,7 @@
}
else if (sides == 8 && FLAG_GET(worked, r + 1, c + 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -7;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -355,7 +376,8 @@
for (r = 0; r < nrows; r++) {
G_percent(r, nrows, 3);
for (c = 0; c < ncols; c++) {
- s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
+ if (er_flag)
+ s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
asp_value = asp[SEG_INDEX(asp_seg, r, c)];
if (r == 0 || c == 0 || r == nrows - 1 ||
c == ncols - 1 || asp_value != 0) {
@@ -363,7 +385,12 @@
if (wat_value > 0) {
wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
}
- if (r == 0)
+ /* set depression */
+ if (asp_value) {
+ asp_value = 0;
+ wat[SEG_INDEX(wat_seg, r, c)] = ABS(wat_value);
+ }
+ else if (r == 0)
asp_value = -2;
else if (c == 0)
asp_value = -4;
@@ -371,17 +398,15 @@
asp_value = -6;
else if (c == ncols - 1)
asp_value = -8;
- else
- asp_value = -1;
asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
}
}
}
}
- G_percent(r, nrows, 3); /* finish it */
+ G_percent(r, nrows, 1); /* finish it */
return 0;
}
Modified: grass/trunk/raster/r.watershed/ram/sg_factor.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/sg_factor.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/ram/sg_factor.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -9,7 +9,7 @@
CELL low_elev, hih_elev;
double height, length, S, sin_theta;
- G_message(_("SECTION 4: Length Slope determination."));
+ G_message(_("SECTION 4: RUSLE LS and/or S factor determination."));
for (r = 0; r < nrows; r++) {
G_percent(r, nrows, 3);
@@ -38,7 +38,7 @@
}
}
}
- G_percent(r, nrows, 3); /* finish it */
+ G_percent(nrows, nrows, 1); /* finish it */
if (ril_flag) {
G_free(ril_buf);
Modified: grass/trunk/raster/r.watershed/ram/slope_len.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/slope_len.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/ram/slope_len.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -20,7 +20,7 @@
if (asp_value == 2 || asp_value == 6)
res = window.ns_res;
else /* asp_value == 4, 8, -2, -4, -6, or -8 */
- res = diag;
+ res = diag; /* how can res be diag with sides == 4??? */
}
else { /* c == dc */
Modified: grass/trunk/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/Gwater.h 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/Gwater.h 2009-01-15 09:15:17 UTC (rev 35412)
@@ -31,7 +31,7 @@
#define POINT struct points
POINT {
- SHORT r, c, downr, downc;
+ SHORT r, c; /* , downr, downc */
int nxt;
};
@@ -87,7 +87,7 @@
/* do_astar.c */
int do_astar(void);
-int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int add_pt(SHORT, SHORT, CELL, CELL);
int drop_pt(void);
int sift_up(int, CELL);
double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
Modified: grass/trunk/raster/r.watershed/seg/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/close_maps.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/close_maps.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -1,69 +1,110 @@
#include "Gwater.h"
#include <unistd.h>
+#include <grass/glocale.h>
int close_maps(void)
{
- struct Colors colors, logcolors;
+ struct Colors colors;
int r, c;
CELL is_swale;
- double dvalue;
+ DCELL *dbuf = NULL;
+ int fd;
struct FPRange accRange;
DCELL min, max;
DCELL clr_min, clr_max;
+ DCELL sum, sum_sqr, stddev, lstddev, dvalue;
dseg_close(&slp);
cseg_close(&alt);
if (wat_flag) {
dseg_write_cellfile(&wat, wat_name);
+ /* get standard deviation */
+ fd = G_open_cell_old(wat_name, "");
+ if (fd < 0) {
+ G_fatal_error(_("unable to open flow accumulation map layer"));
+ }
+
+ sum = sum_sqr = stddev = 0.0;
+ dbuf = G_allocate_d_raster_buf();
+ for (r = 0; r < nrows; r++) {
+ G_get_d_raster_row(fd, dbuf, r);
+ for (c = 0; c < ncols; c++) {
+ dvalue = dbuf[c];
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dvalue = ABS(dvalue);
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
+ }
+ }
+
+ 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);
+
G_read_fp_range(wat_name, this_mapset, &accRange);
min = max = 0;
G_get_fp_range_min_max(&accRange, &min, &max);
G_init_colors(&colors);
- if (ABS(min) > max) {
- clr_min = 1;
- clr_max = 0.25 * max;
- G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0, 255,
- 0, &colors);
- clr_min = 0.25 * max;
- clr_max = 0.5 * max;
- G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0, 255,
- 255, &colors);
- clr_min = 0.5 * max;
- clr_max = 0.75 * max;
- G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0, 0,
- 255, &colors);
- clr_min = 0.75 * max;
+
+ if (min < 0) {
+ if (min < (-stddev - 1)) {
+ clr_min = min;
+ clr_max = -stddev - 1;
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
+ 0, 0, &colors);
+ }
+ clr_min = -stddev - 1.;
+ clr_max = -1. * exp(lstddev * 0.75);
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
+ 0, 255, &colors);
+ clr_min = clr_max;
+ clr_max = -1. * exp(lstddev * 0.5);
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0,
+ 255, 255, &colors);
+ clr_min = clr_max;
+ clr_max = -1. * exp(lstddev * 0.35);
+ G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
+ 255, 0, &colors);
+ clr_min = clr_max;
+ clr_max = -1.;
+ G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 255,
+ 255, 0, &colors);
+ }
+ clr_min = -1.;
+ clr_max = 1.;
+ G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 255,
+ 255, 0, &colors);
+ clr_min = 1;
+ clr_max = exp(lstddev * 0.35);
+ G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
+ 255, 0, &colors);
+ clr_min = clr_max;
+ clr_max = exp(lstddev * 0.5);
+ G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0,
+ 255, 255, &colors);
+ clr_min = clr_max;
+ clr_max = exp(lstddev * 0.75);
+ G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
+ 0, 255, &colors);
+ clr_min = clr_max;
+ clr_max = stddev + 1.;
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
+ 0, &colors);
+
+ if (max > 0 && max > stddev + 1) {
+ clr_min = stddev + 1;
clr_max = max;
- G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0, 0,
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0, 0, 0,
&colors);
- max = -max;
}
- else {
- min = ABS(min);
- clr_min = 1;
- clr_max = 0.25 * min;
- G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0, 255,
- 0, &colors);
- clr_min = 0.25 * min;
- clr_max = 0.5 * min;
- G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0, 255,
- 255, &colors);
- clr_min = 0.5 * min;
- clr_max = 0.75 * min;
- G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0, 0,
- 255, &colors);
- clr_min = 0.75 * min;
- clr_max = min;
- G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0, 0,
- &colors);
- }
- G_abs_log_colors(&logcolors, &colors, 5);
- G_add_d_raster_color_rule(&min, 0, 0, 0, &max, 0, 0, 0, &logcolors);
- G_write_colors(wat_name, this_mapset, &logcolors);
-
+ G_write_colors(wat_name, this_mapset, &colors);
}
if (asp_flag) {
cseg_write_cellfile(&asp, asp_name);
Modified: grass/trunk/raster/r.watershed/seg/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_astar.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/do_astar.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -46,7 +46,8 @@
c = point.c;
bseg_put(&worked, &one, r, c);
- cseg_get(&alt, &alt_val, r, c);
+ /* cseg_get(&alt, &alt_val, r, c); */
+ alt_val = heap_pos.ele;
/* check all neighbours, breadth first search */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -60,26 +61,27 @@
bseg_get(&in_list, &in_val, upr, upc);
if (in_val == 0) {
cseg_get(&alt, &alt_up, upr, upc);
- add_pt(upr, upc, r, c, alt_up, alt_val);
+ add_pt(upr, upc, alt_up, alt_val);
/* flow direction is set here */
drain_val = drain[upr - r + 1][upc - c + 1];
cseg_put(&asp, &drain_val, upr, upc);
}
else {
/* check if neighbour has not been worked on,
- * update values for asp, wat and slp */
+ * update values for asp and wat */
bseg_get(&worked, &work_val, upr, upc);
if (!work_val) {
cseg_get(&asp, &asp_up, upr, upc);
- if (asp_up < -1) {
+ if (asp_up < 0) {
drain_val = drain[upr - r + 1][upc - c + 1];
cseg_put(&asp, &drain_val, upr, upc);
dseg_get(&wat, &wat_val, r, c);
- if (wat_val > 0)
+ if (wat_val > 0) {
wat_val = -wat_val;
- dseg_put(&wat, &wat_val, r, c);
- cseg_get(&alt, &alt_up, upr, upc);
- replace(upr, upc, r, c); /* alt_up used to be */
+ dseg_put(&wat, &wat_val, r, c);
+ }
+ /* cseg_get(&alt, &alt_up, upr, upc); */
+ /* replace(upr, upc, r, c); */ /* alt_up used to be */
/* slope = get_slope (upr, upc, r, c, alt_up, alt_val);
dseg_put (&slp, &slope, upr, upc); */
}
@@ -100,7 +102,7 @@
}
/* new add point routine for min heap */
-int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
{
POINT point;
HEAP heap_pos;
@@ -124,11 +126,13 @@
point.r = r;
point.c = c;
- point.downr = downr;
- point.downc = downc;
+/* point.downr = downr;
+ point.downc = downc; */
seg_put(&astar_pts, (char *)&point, 0, nxt_avail_pt);
+ /* cseg_put(&pnt_index, &nxt_avail_pt, r, c); */
+
nxt_avail_pt++;
/* sift up: move new point towards top of heap */
@@ -142,8 +146,8 @@
int drop_pt(void)
{
int child, childr, parent;
- int childp, childrp;
- CELL ele, eler;
+ int childp; /* , childrp */
+ CELL ele; /* , eler */
int i;
HEAP heap_pos;
@@ -170,23 +174,23 @@
ele = heap_pos.ele;
if (child < heap_size) {
childr = child + 1;
- i = 1;
- while (childr <= heap_size && i < 3) {
+ i = child + 3;
+ while (childr <= heap_size && childr < i) {
seg_get(&heap_index, (char *)&heap_pos, 0, childr);
- childrp = heap_pos.point;
- eler = heap_pos.ele;
- if (eler < ele) {
+ /* childrp = heap_pos.point;
+ eler = heap_pos.ele; */
+ if (heap_pos.ele < ele) {
child = childr;
- childp = childrp;
- ele = eler;
+ childp = heap_pos.point;
+ ele = heap_pos.ele;
}
/* make sure we get the oldest child */
- else if (ele == eler && childp > childrp) {
+ else if (ele == heap_pos.ele && childp > heap_pos.point) {
child = childr;
- childp = childrp;
+ childp = heap_pos.point;
}
childr++;
- i++;
+ /* i++; */
}
}
@@ -220,8 +224,8 @@
/* standard sift-up routine for d-ary min heap */
int sift_up(int start, CELL ele)
{
- int parent, parentp, child, childp;
- CELL elep;
+ int parent, child, childp; /* parentp, */
+ /* CELL elep; */
HEAP heap_pos;
child = start;
@@ -231,11 +235,11 @@
while (child > 1) {
parent = GET_PARENT(child);
seg_get(&heap_index, (char *)&heap_pos, 0, parent);
- parentp = heap_pos.point;
- elep = heap_pos.ele;
+ /* parentp = heap_pos.point;
+ elep = heap_pos.ele; */
/* parent ele higher */
- if (elep > ele) {
+ if (heap_pos.ele > ele) {
/* push parent point down */
seg_put(&heap_index, (char *)&heap_pos, 0, child);
@@ -243,7 +247,7 @@
}
/* same ele, but parent is younger */
- else if (elep == ele && parentp > childp) {
+ else if (heap_pos.ele == ele && heap_pos.point > childp) {
/* push parent point down */
seg_put(&heap_index, (char *)&heap_pos, 0, child);
@@ -280,28 +284,23 @@
return (slope);
}
-int replace( /* ele was in there */
- SHORT upr, SHORT upc, SHORT r, SHORT c)
-/* CELL ele; */
-{
- int now, heap_run;
+int replace(SHORT upr, SHORT upc, SHORT r, SHORT c)
+{ /* ele was in there */
+ /* CELL ele; */
+ int now;
POINT point;
- HEAP heap_pos;
- heap_run = 0;
+ now = 0;
- /* this is dumb, works for ram, for seg make new index pt_index[row][col] */
- while (heap_run <= heap_size) {
- seg_get(&heap_index, (char *)&heap_pos, 0, heap_run);
- now = heap_pos.point;
- seg_get(&astar_pts, (char *)&point, 0, now);
- if (point.r == upr && point.c == upc) {
- point.downr = r;
- point.downc = c;
- seg_put(&astar_pts, (char *)&point, 0, now);
- return 0;
- }
- heap_run++;;
+ /* cseg_get(&pnt_index, &now, upr, upc); */
+ seg_get(&astar_pts, (char *)&point, 0, now);
+ if (point.r != upr || point.c != upc) {
+ G_warning("pnt_index incorrect!");
+ return 1;
}
+/* point.downr = r;
+ point.downc = c; */
+ seg_put(&astar_pts, (char *)&point, 0, now);
+
return 0;
}
Modified: grass/trunk/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_cum.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/do_cum.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -7,13 +7,14 @@
int do_cum(void)
{
SHORT r, c, dr, dc;
- CELL is_swale;
+ CELL is_swale, asp_val, asp_val_down;
DCELL value, valued;
POINT point;
int killer, threshold, count;
+ SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
- G_debug(1, "MFD convergence factor set to %d.", c_fac);
count = 0;
if (bas_thres <= 0)
@@ -25,10 +26,16 @@
killer = first_cum;
seg_get(&astar_pts, (char *)&point, 0, killer);
first_cum = point.nxt;
- if ((dr = point.downr) > -1) {
- r = point.r;
- c = point.c;
- dc = point.downc;
+ r = point.r;
+ c = point.c;
+ cseg_get(&asp, &asp_val, r, c);
+ if (asp_val) {
+ dr = r + asp_r[ABS(asp_val)];
+ dc = c + asp_c[ABS(asp_val)];
+ }
+ else
+ dr = dc = -1;
+ if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = point.downr) > -1) { */
dseg_get(&wat, &value, r, c);
if (ABS(value) >= threshold)
bseg_put(&swale, &one, r, c);
@@ -47,6 +54,14 @@
}
dseg_put(&wat, &valued, dr, dc);
bseg_get(&swale, &is_swale, r, c);
+ /* update asp for depression */
+ if (is_swale && pit_flag) {
+ cseg_get(&asp, &asp_val_down, dr, dc);
+ if (asp_val > 0 && asp_val_down == 0) {
+ asp_val = -asp_val;
+ cseg_put(&asp, &asp_val, r, c);
+ }
+ }
if (is_swale || ABS(valued) >= threshold) {
bseg_put(&swale, &one, dr, dc);
}
@@ -88,20 +103,23 @@
{
int r, c, dr, dc;
CELL is_swale;
- DCELL value, valued;
+ DCELL value, valued, *wat_nbr;
POINT point;
int killer, threshold, count;
/* MFD */
- int mfd_cells, stream_cells, swale_cells, astar_not_set;
+ 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, r_swale, c_swale, ct_dir, np_side;
+ int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
double dx, dy;
- CELL ele, ele_nbr, aspect, asp_val, is_worked;
- double prop, max_acc, max_swale;
+ CELL ele, ele_nbr, asp_val, asp_val2, cvalue, *worked_nbr;
+ double prop, max_acc;
int workedon;
+ SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
+ G_debug(1, "MFD convergence factor set to %d.", c_fac);
/* distances to neighbours */
dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
@@ -127,6 +145,9 @@
}
}
+ worked_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
+ wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
+
workedon = 0;
count = 0;
@@ -140,26 +161,32 @@
killer = first_cum;
seg_get(&astar_pts, (char *)&point, 0, killer);
first_cum = point.nxt;
- if ((dr = point.downr) > -1) {
- r = point.r;
- c = point.c;
- dc = point.downc;
+ r = point.r;
+ c = point.c;
+ cseg_get(&asp, &asp_val, r, c);
+ if (asp_val) {
+ dr = r + asp_r[ABS(asp_val)];
+ dc = c + asp_c[ABS(asp_val)];
+ }
+ else
+ dr = dc = -1;
+ if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = point.downr) > -1) { */
+ /* dc = point.downc; */
dseg_get(&wat, &value, r, c);
- dseg_get(&wat, &valued, dr, dc);
- /* disabled for MFD */
- /* if (ABS(value) > threshold)
- bseg_put(&swale, &one, r, c); */
+ r_max = dr;
+ c_max = dc;
- /* start MFD */
-
/* get weights */
max_weight = 0;
sum_weight = 0;
np_side = -1;
mfd_cells = 0;
+ stream_cells = 0;
+ swale_cells = 0;
astar_not_set = 1;
cseg_get(&alt, &ele, r, c);
+ is_null = 0;
/* this loop is needed to get the sum of weights */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (r_nbr, c_nbr) for neighbours */
@@ -169,19 +196,39 @@
/* check that neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols) {
- bseg_get(&worked, &is_worked, r_nbr, c_nbr);
- if (is_worked == 0) {
+
+ /* check for swale or stream cells */
+ bseg_get(&swale, &is_swale, r_nbr, c_nbr);
+ if (is_swale)
+ swale_cells++;
+ dseg_get(&wat, &valued, r_nbr, c_nbr);
+ wat_nbr[ct_dir] = valued;
+ if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold)
+ stream_cells++;
+
+ bseg_get(&worked, &cvalue, r_nbr, c_nbr);
+ worked_nbr[ct_dir] = cvalue;
+ if (worked_nbr[ct_dir] == 0) {
cseg_get(&alt, &ele_nbr, r_nbr, c_nbr);
- if (ele_nbr < ele) {
- weight[ct_dir] =
- mfd_pow(((ele -
- ele_nbr) / dist_to_nbr[ct_dir]),
- c_fac);
+ is_null = G_is_c_null_value(&ele_nbr);
+ if (!is_null && ele_nbr <= ele) {
+ if (ele_nbr < ele) {
+ weight[ct_dir] =
+ mfd_pow(((ele -
+ ele_nbr) / dist_to_nbr[ct_dir]),
+ c_fac);
+ }
+ if (ele_nbr == ele) {
+ weight[ct_dir] =
+ mfd_pow((0.5 / dist_to_nbr[ct_dir]),
+ c_fac);
+ }
sum_weight += weight[ct_dir];
mfd_cells++;
- if (weight[ct_dir] > max_weight)
+ if (weight[ct_dir] > max_weight) {
max_weight = weight[ct_dir];
+ }
if (dr == r_nbr && dc == c_nbr) {
astar_not_set = 0;
@@ -196,7 +243,6 @@
/* honour A * path
* mfd_cells == 0: fine, SFD along A * path
* mfd_cells == 1 && astar_not_set == 0: fine, SFD along A * path
- * mfd_cells > 1 && astar_not_set == 0: MFD, set A * path to max_weight
* mfd_cells > 0 && astar_not_set == 1: A * path not included, add to mfd_cells
*/
@@ -207,21 +253,8 @@
weight[np_side] = max_weight;
}
- /* MFD, A * path included, set A * path to max_weight */
- if (mfd_cells > 1 && astar_not_set == 0) {
- sum_weight = sum_weight - weight[np_side] + max_weight;
- weight[np_side] = max_weight;
- }
-
/* set flow accumulation for neighbours */
- stream_cells = 0;
- swale_cells = 0;
max_acc = -1;
- max_swale = -1;
- r_max = dr;
- c_max = dc;
- r_swale = dr;
- c_swale = dc;
if (mfd_cells > 1) {
prop = 0.0;
@@ -233,44 +266,38 @@
/* check that neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols && weight[ct_dir] > -0.5) {
- bseg_get(&worked, &is_worked, r_nbr, c_nbr);
- bseg_get(&swale, &is_swale, r_nbr, c_nbr);
- if (is_worked == 0) {
+ /* bseg_get(&worked, &is_worked, r_nbr, c_nbr); */
+ if (worked_nbr[ct_dir] == 0) {
weight[ct_dir] = weight[ct_dir] / sum_weight;
- prop += weight[ct_dir]; /* check everything sums up to 1.0 */
+ /* check everything sums up to 1.0 */
+ prop += weight[ct_dir];
- dseg_get(&wat, &valued, r_nbr, c_nbr);
if (value > 0) {
- if (valued > 0)
- valued += value * weight[ct_dir];
+ if (wat_nbr[ct_dir] > 0)
+ wat_nbr[ct_dir] += value * weight[ct_dir];
else
- valued -= value * weight[ct_dir];
+ wat_nbr[ct_dir] -= value * weight[ct_dir];
}
else {
- if (valued < 0)
- valued += value * weight[ct_dir];
+ if (wat_nbr[ct_dir] < 0)
+ wat_nbr[ct_dir] += value * weight[ct_dir];
else
- valued = value * weight[ct_dir] - valued;
+ wat_nbr[ct_dir] = value * weight[ct_dir] - wat_nbr[ct_dir];
}
+ valued = wat_nbr[ct_dir];
dseg_put(&wat, &valued, r_nbr, c_nbr);
- if ((ABS(valued) + 0.5) >= threshold)
- stream_cells++;
- if (ABS(valued) > max_acc) {
- max_acc = ABS(valued);
+ /* get main drainage direction */
+ if (ABS(wat_nbr[ct_dir]) >= max_acc) {
+ max_acc = ABS(wat_nbr[ct_dir]);
r_max = r_nbr;
c_max = c_nbr;
}
- /* assist in stream merging */
- if (is_swale && ABS(valued) > max_swale) {
- max_swale = ABS(valued);
- r_swale = r_nbr;
- c_swale = c_nbr;
- }
}
else if (ct_dir == np_side) {
- workedon++; /* check for consistency with A * path */
+ /* check for consistency with A * path */
+ workedon++;
}
}
}
@@ -278,14 +305,10 @@
G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
prop);
}
- dseg_get(&wat, &valued, dr, dc);
- if (max_swale > -0.5) {
- r_max = r_swale;
- c_max = c_swale;
- }
}
if (mfd_cells < 2) {
+ dseg_get(&wat, &valued, dr, dc);
if (value > 0) {
if (valued > 0)
valued += value;
@@ -301,14 +324,30 @@
dseg_put(&wat, &valued, dr, dc);
}
- /* MFD finished */
+ /* update asp */
+ if (dr != r_max || dc != c_max) {
+ asp_val2 = drain[r - r_max + 1][c - c_max + 1];
+ /* cseg_get(&asp, &asp_val, r, c); */
+ if (asp_val < 0)
+ asp_val2 = -asp_val2;
+ cseg_put(&asp, &asp_val2, r, c);
- valued = ABS(valued) + 0.5;
+ }
+ /* update asp for depression */
bseg_get(&swale, &is_swale, r, c);
- /* start new stream: only works with A * path */
- if (!is_swale && (int)valued >= threshold && stream_cells < 2 &&
+ if (is_swale && pit_flag) {
+ cseg_get(&asp, &asp_val2, r_max, c_max);
+ if (asp_val > 0 && asp_val2 == 0) {
+ asp_val = -asp_val;
+ cseg_put(&asp, &asp_val, r, c);
+ }
+ }
+ /* start new stream */
+ value = ABS(value) + 0.5;
+ if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
swale_cells < 1) {
- bseg_put(&swale, &one, dr, dc);
+ bseg_put(&swale, &one, r, c);
+ is_swale = 1;
}
/* continue stream */
if (is_swale) {
@@ -318,26 +357,22 @@
if (er_flag && !is_swale)
slope_length(r, c, r_max, c_max);
}
- /* update asp */
- if (dr != r_max || dc != c_max) {
- aspect = drain[r - r_max + 1][c - c_max + 1];
- cseg_get(&asp, &asp_val, r, c);
- if (asp_val < 0)
- aspect = -aspect;
- cseg_put(&asp, &aspect, r, c);
-
- }
bseg_put(&worked, &one, r, c);
}
}
G_percent(count, do_points, 1); /* finish it */
if (workedon)
- G_warning(_("MFD: A * path already processed when distributing flow: %d of %d"),
+ G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
workedon, do_points);
seg_close(&astar_pts);
bseg_close(&worked);
+
+ G_free(dist_to_nbr);
+ G_free(weight);
+ G_free(wat_nbr);
+ G_free(worked_nbr);
return 0;
}
Modified: grass/trunk/raster/r.watershed/seg/dseg_get.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/dseg_get.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/dseg_get.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -4,7 +4,7 @@
int dseg_get(DSEG * dseg, double *value, int row, int col)
{
- if (segment_get(&(dseg->seg), (CELL *) value, row, col) < 0) {
+ if (segment_get(&(dseg->seg), (DCELL *) value, row, col) < 0) {
G_warning("dseg_get(): could not read segment file");
return -1;
}
Modified: grass/trunk/raster/r.watershed/seg/dseg_put.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/dseg_put.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/dseg_put.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -4,7 +4,7 @@
int dseg_put(DSEG * dseg, double *value, int row, int col)
{
- if (segment_put(&(dseg->seg), (CELL *) value, row, col) < 0) {
+ if (segment_put(&(dseg->seg), (DCELL *) value, row, col) < 0) {
G_warning("dseg_put(): could not write segment file");
return -1;
}
Modified: grass/trunk/raster/r.watershed/seg/dseg_read.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/dseg_read.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/dseg_read.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -6,10 +6,9 @@
int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
{
- int row, col, nrows, ncols;
+ int row, nrows, ncols;
int map_fd;
char msg[100];
- CELL *buffer;
double *dbuffer;
dseg->name = NULL;
@@ -24,11 +23,9 @@
}
nrows = G_window_rows();
ncols = G_window_cols();
- buffer = G_allocate_cell_buf();
- dbuffer = (double *)G_malloc(ncols * sizeof(double));
+ dbuffer = G_allocate_d_raster_buf();
for (row = 0; row < nrows; row++) {
- if (G_get_c_raster_row(map_fd, buffer, row) < 0) {
- G_free(buffer);
+ if (G_get_d_raster_row(map_fd, dbuffer, row) < 0) {
G_free(dbuffer);
G_close_cell(map_fd);
sprintf(msg, "%s(): unable to read file [%s] in [%s], %d %d",
@@ -36,11 +33,7 @@
G_warning(msg);
return -2;
}
- for (col = ncols - 1; col >= 0; col--) {
- dbuffer[col] = (double)buffer[col];
- }
- if (segment_put_row(&(dseg->seg), (CELL *) dbuffer, row) < 0) {
- G_free(buffer);
+ if (segment_put_row(&(dseg->seg), (DCELL *) dbuffer, row) < 0) {
G_free(dbuffer);
G_close_cell(map_fd);
sprintf(msg, "%s(): unable to segment put row for [%s] in [%s]",
@@ -51,7 +44,6 @@
}
G_close_cell(map_fd);
- G_free(buffer);
G_free(dbuffer);
dseg->name = G_store(map_name);
Modified: grass/trunk/raster/r.watershed/seg/dseg_write.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/dseg_write.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/dseg_write.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -7,27 +7,21 @@
int dseg_write_cellfile(DSEG * dseg, char *map_name)
{
int map_fd;
- int row, col, nrows, ncols;
- CELL *buffer;
+ int row, nrows, ncols;
double *dbuffer;
- map_fd = G_open_cell_new(map_name);
+ map_fd = G_open_raster_new(map_name, DCELL_TYPE);
if (map_fd < 0) {
G_warning("%s(): unable to open new map layer [%s]", me, map_name);
return -1;
}
nrows = G_window_rows();
ncols = G_window_cols();
- buffer = G_allocate_cell_buf();
- dbuffer = (double *)G_malloc(ncols * sizeof(double));
+ dbuffer = G_allocate_d_raster_buf();
segment_flush(&(dseg->seg));
for (row = 0; row < nrows; row++) {
- segment_get_row(&(dseg->seg), (CELL *) dbuffer, row);
- for (col = ncols - 1; col >= 0; col--) {
- buffer[col] = (CELL) (dbuffer[col] + 0.5);
- }
- if (G_put_raster_row(map_fd, buffer, CELL_TYPE) < 0) {
- G_free(buffer);
+ segment_get_row(&(dseg->seg), (DCELL *) dbuffer, row);
+ if (G_put_raster_row(map_fd, dbuffer, DCELL_TYPE) < 0) {
G_free(dbuffer);
G_unopen_cell(map_fd);
G_warning("%s(): unable to write new map layer [%s], row %d",
@@ -35,7 +29,6 @@
return -2;
}
}
- G_free(buffer);
G_free(dbuffer);
G_close_cell(map_fd);
return 0;
Modified: grass/trunk/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/init_vars.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/init_vars.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -15,7 +15,7 @@
/* int page_block, num_cseg; */
int max_bytes;
- CELL *buf, alt_value, asp_value, worked_value;
+ CELL *buf, alt_value, asp_value, worked_value, block_value;
DCELL wat_value;
char MASK_flag;
@@ -131,8 +131,7 @@
/* segment parameters: one size fits all. Fine tune? */
/* Segment rows and cols: 200 */
- /* 1 segment open for all rasters: 2.86 MB */
- /* num_open_segs = segs_mb / 2.86 */
+ /* 1 segment open for all rasters: 1.34 MB */
seg_rows = SROW;
seg_cols = SCOL;
@@ -142,7 +141,7 @@
G_warning(_("Maximum memory to be used was smaller than 3 MB, set to default = 300 MB."));
}
- num_open_segs = segs_mb / 2.86;
+ num_open_segs = segs_mb / 1.34;
G_debug(1, "segs MB: %.0f", segs_mb);
G_debug(1, "region rows: %d", nrows);
@@ -166,24 +165,92 @@
G_debug(1, " open segments after adjusting:\t%d", num_open_segs);
cseg_open(&alt, seg_rows, seg_cols, num_open_segs);
- cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
cseg_read_cell(&alt, ele_name, "");
- cseg_read_cell(&r_h, ele_name, "");
+ if (er_flag) {
+ cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
+ cseg_read_cell(&r_h, ele_name, "");
+ }
+
+ /* NULL cells in input elevation map */
+ bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
+ bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
+ G_debug(1, "Checking for masked and NULL cells in input elevation <%s>.", ele_name);
+ MASK_flag = 0;
+ do_points = nrows * ncols;
+ fd = G_open_cell_old(ele_name, "");
+ if (fd < 0) {
+ G_fatal_error(_("unable to open elevation map layer"));
+ }
+ buf = G_allocate_cell_buf();
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ /* check for masked and NULL cells */
+ alt_value = buf[c];
+ if (G_is_c_null_value(&alt_value)) {
+ bseg_put(&worked, &one, r, c);
+ bseg_put(&in_list, &one, r, c);
+ do_points--;
+ }
+ }
+ }
+ G_close_cell(fd);
+ G_free(buf);
+ if (do_points < nrows * ncols)
+ MASK_flag = 1;
+
+ /* initial flow accumulation */
dseg_open(&wat, seg_rows, seg_cols, num_open_segs);
-
if (run_flag) {
dseg_read_cell(&wat, run_name, "");
+ if (MASK_flag) {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ bseg_get(&worked, &worked_value, r, c);
+ if (worked_value)
+ dseg_put(&wat, &d_zero, r, c);
+ }
+ }
+ }
}
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++)
- if (-1 == dseg_put(&wat, &d_one, r, c))
- exit(EXIT_FAILURE);
+ if (MASK_flag) {
+ bseg_get(&worked, &worked_value, r, c);
+ if (worked_value)
+ dseg_put(&wat, &d_zero, r, c);
+ else
+ dseg_put(&wat, &d_one, r, c);
+ }
+ else {
+ if (-1 == dseg_put(&wat, &d_one, r, c))
+ exit(EXIT_FAILURE);
+ }
}
}
cseg_open(&asp, seg_rows, seg_cols, num_open_segs);
+ /* depression: drainage direction will be set to zero later */
if (pit_flag) {
- cseg_read_cell(&asp, pit_name, "");
+ fd = G_open_cell_old(pit_name, "");
+ if (fd < 0) {
+ G_fatal_error(_("unable to open depression map layer"));
+ }
+ buf = G_allocate_cell_buf();
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ asp_value = buf[c];
+ if (!G_is_c_null_value(&asp_value)) {
+ cseg_put(&asp, &one, r, c);
+ }
+ else {
+ cseg_put(&asp, &zero, r, c);
+ }
+ }
+ }
+ G_close_cell(fd);
+ G_free(buf);
}
else {
for (r = 0; r < nrows; r++) {
@@ -194,7 +261,25 @@
}
bseg_open(&swale, seg_rows, seg_cols, num_open_segs);
if (ob_flag) {
- bseg_read_cell(&swale, ob_name, "");
+ fd = G_open_cell_old(ob_name, "");
+ if (fd < 0) {
+ G_fatal_error(_("unable to open blocking map layer"));
+ }
+ buf = G_allocate_cell_buf();
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ block_value = buf[c];
+ if (!G_is_c_null_value(&block_value)) {
+ bseg_put(&swale, &one, r, c);
+ }
+ else {
+ bseg_put(&swale, &zero, r, c);
+ }
+ }
+ }
+ G_close_cell(fd);
+ G_free(buf);
}
else {
for (r = 0; r < nrows; r++) {
@@ -206,44 +291,26 @@
dseg_open(&ril, 1, seg_rows * seg_cols, num_open_segs);
dseg_read_cell(&ril, ril_name, "");
}
- bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
- bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
- MASK_flag = 0;
- do_points = nrows * ncols;
- if (NULL != G_find_file("cell", "MASK", G_mapset())) {
- MASK_flag = 1;
- if ((fd = G_open_cell_old("MASK", G_mapset())) < 0) {
- G_fatal_error(_("Unable to open MASK"));
- }
- else {
- buf = G_allocate_cell_buf();
- for (r = 0; r < nrows; r++) {
- G_get_c_raster_row_nomask(fd, buf, r);
- for (c = 0; c < ncols; c++) {
- if (!buf[c]) {
- do_points--;
- bseg_put(&worked, &one, r, c);
- bseg_put(&in_list, &one, r, c);
- }
- }
- }
- G_close_cell(fd);
- G_free(buf);
- }
+
+ /* dseg_open(&slp, SROW, SCOL, num_open_segs); */
+
+ /* RUSLE: LS and/or S factor */
+
+ if (er_flag) {
+ dseg_open(&s_l, seg_rows, seg_cols, num_open_segs);
}
- /* dseg_open(&slp, SROW, SCOL, num_open_segs); */
- 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);
if (ls_flag)
dseg_open(&l_s, 1, seg_rows * seg_cols, num_open_segs);
- seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols,
- num_open_segs, sizeof(POINT));
+ seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols * 2,
+ num_open_segs / 2, sizeof(POINT));
+
/* heap_index will track astar_pts in ternary min-heap */
/* heap_index is one-based */
- seg_open(&heap_index, 1, do_points + 1, 1, seg_rows * seg_cols,
- num_open_segs, sizeof(HEAP));
+ seg_open(&heap_index, 1, do_points + 1, 1, seg_cols * num_open_segs * seg_rows / 10,
+ 10, sizeof(HEAP));
G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
@@ -261,7 +328,8 @@
dseg_put(&wat, &d_zero, r, c);
}
else {
- dseg_put(&s_l, &half_res, r, c);
+ if (er_flag)
+ dseg_put(&s_l, &half_res, r, c);
cseg_get(&asp, &asp_value, r, c);
if (r == 0 || c == 0 || r == nrows - 1 ||
c == ncols - 1 || asp_value != 0) {
@@ -270,7 +338,15 @@
wat_value = -wat_value;
dseg_put(&wat, &wat_value, r, c);
}
- if (r == 0)
+ /* set depression */
+ if (asp_value) {
+ asp_value = 0;
+ if (wat_value < 0) {
+ wat_value = -wat_value;
+ dseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else if (r == 0)
asp_value = -2;
else if (c == 0)
asp_value = -4;
@@ -278,17 +354,15 @@
asp_value = -6;
else if (c == ncols - 1)
asp_value = -8;
- else
- asp_value = -1;
if (-1 == cseg_put(&asp, &asp_value, r, c))
exit(EXIT_FAILURE);
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
}
else if (!bseg_get(&worked, &worked_value, r - 1, c)
&& worked_value != 0) {
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp_value = -2;
cseg_put(&asp, &asp_value, r, c);
dseg_get(&wat, &wat_value, r, c);
@@ -300,7 +374,7 @@
else if (!bseg_get(&worked, &worked_value, r + 1, c)
&& worked_value != 0) {
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp_value = -6;
cseg_put(&asp, &asp_value, r, c);
dseg_get(&wat, &wat_value, r, c);
@@ -312,7 +386,7 @@
else if (!bseg_get(&worked, &worked_value, r, c - 1)
&& worked_value != 0) {
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp_value = -4;
cseg_put(&asp, &asp_value, r, c);
dseg_get(&wat, &wat_value, r, c);
@@ -324,7 +398,7 @@
else if (!bseg_get(&worked, &worked_value, r, c + 1)
&& worked_value != 0) {
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp_value = -8;
cseg_put(&asp, &asp_value, r, c);
dseg_get(&wat, &wat_value, r, c);
@@ -337,7 +411,7 @@
!bseg_get(&worked, &worked_value, r - 1, c - 1)
&& worked_value != 0) {
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp_value = -3;
cseg_put(&asp, &asp_value, r, c);
dseg_get(&wat, &wat_value, r, c);
@@ -350,7 +424,7 @@
!bseg_get(&worked, &worked_value, r - 1, c + 1)
&& worked_value != 0) {
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp_value = -1;
cseg_put(&asp, &asp_value, r, c);
dseg_get(&wat, &wat_value, r, c);
@@ -363,7 +437,7 @@
!bseg_get(&worked, &worked_value, r + 1, c - 1)
&& worked_value != 0) {
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp_value = -5;
cseg_put(&asp, &asp_value, r, c);
dseg_get(&wat, &wat_value, r, c);
@@ -376,7 +450,7 @@
!bseg_get(&worked, &worked_value, r + 1, c + 1)
&& worked_value != 0) {
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp_value = -7;
cseg_put(&asp, &asp_value, r, c);
dseg_get(&wat, &wat_value, r, c);
@@ -398,7 +472,8 @@
G_percent(r, nrows, 2);
for (c = 0; c < ncols; c++) {
bseg_put(&worked, &zero, r, c);
- dseg_put(&s_l, &half_res, r, c);
+ if (er_flag)
+ dseg_put(&s_l, &half_res, r, c);
cseg_get(&asp, &asp_value, r, c);
if (r == 0 || c == 0 || r == nrows - 1 ||
c == ncols - 1 || asp_value != 0) {
@@ -408,7 +483,15 @@
if (-1 == dseg_put(&wat, &wat_value, r, c))
exit(EXIT_FAILURE);
}
- if (r == 0)
+ /* set depression */
+ if (asp_value) {
+ asp_value = 0;
+ if (wat_value < 0) {
+ wat_value = -wat_value;
+ dseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else if (r == 0)
asp_value = -2;
else if (c == 0)
asp_value = -4;
@@ -416,12 +499,10 @@
asp_value = -6;
else if (c == ncols - 1)
asp_value = -8;
- else
- asp_value = -1;
if (-1 == cseg_put(&asp, &asp_value, r, c))
exit(EXIT_FAILURE);
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
}
else {
bseg_put(&in_list, &zero, r, c);
Modified: grass/trunk/raster/r.watershed/seg/main.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/main.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/main.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -47,7 +47,6 @@
char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
char run_name[GNAME_MAX], ob_name[GNAME_MAX];
char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
-
const char *this_mapset;
char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX],
thr_name[8];
Modified: grass/trunk/raster/r.watershed/seg/sg_factor.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/sg_factor.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/sg_factor.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -9,10 +9,9 @@
CELL downer, low_elev, hih_elev;
double height, length, S, sin_theta;
- G_message(_("SECTION 4: Length Slope determination."));
+ G_message(_("SECTION 4: RUSLE LS and/or S factor determination."));
for (r = nrows - 1; r >= 0; r--) {
-/* FIXME: G_percent params count backwards */
- G_percent(r, nrows, 3);
+ G_percent(nrows - r, nrows, 3);
for (c = ncols - 1; c >= 0; c--) {
cseg_get(&alt, &low_elev, r, c);
cseg_get(&r_h, &hih_elev, r, c);
@@ -38,7 +37,7 @@
}
}
}
- G_percent(r, nrows, 3); /* finish it */
+ G_percent(nrows, nrows, 1); /* finish it */
return 0;
}
Modified: grass/trunk/raster/r.watershed/seg/slope_len.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/slope_len.c 2009-01-14 20:52:38 UTC (rev 35411)
+++ grass/trunk/raster/r.watershed/seg/slope_len.c 2009-01-15 09:15:17 UTC (rev 35412)
@@ -20,7 +20,7 @@
if (asp_value == 2 || asp_value == 6)
res = window.ns_res;
else /* asp_value == 4, 8, -2, -4, -6, or -8 */
- res = diag;
+ res = diag; /* how can res be diag with sides == 4??? */
}
else { /* c == dc */
More information about the grass-commit
mailing list