[GRASS-SVN] r34789 - in grass/trunk/raster/r.watershed: front ram
seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Dec 8 04:11:19 EST 2008
Author: glynn
Date: 2008-12-08 04:11:19 -0500 (Mon, 08 Dec 2008)
New Revision: 34789
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_astar.h
grass/trunk/raster/r.watershed/ram/do_cum.c
grass/trunk/raster/r.watershed/ram/find_pour.c
grass/trunk/raster/r.watershed/ram/init_vars.c
grass/trunk/raster/r.watershed/ram/main.c
grass/trunk/raster/r.watershed/ram/no_stream.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_astar.h
grass/trunk/raster/r.watershed/seg/do_cum.c
grass/trunk/raster/r.watershed/seg/init_vars.c
grass/trunk/raster/r.watershed/seg/main.c
grass/trunk/raster/r.watershed/seg/no_stream.c
Log:
r.watershed MFD support from Markus Metz
Modified: grass/trunk/raster/r.watershed/front/main.c
===================================================================
--- grass/trunk/raster/r.watershed/front/main.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/front/main.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -4,7 +4,8 @@
* MODULE: front end
* AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
* Brad Douglas <rez touchofmadness.com>,
- * Hamish Bowman <hamish_b yahoo.com>
+ * Hamish Bowman <hamish_b yahoo.com>
+ * Markus Metz <markus.metz.giswork gmail.com>
* PURPOSE: Watershed determination
* COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
*
@@ -41,6 +42,8 @@
struct Option *opt14;
struct Option *opt15;
struct Option *opt16;
+ struct Option *opt17;
+ struct Flag *flag_sfd;
struct Flag *flag_flow;
struct Flag *flag_seg;
struct GModule *module;
@@ -178,13 +181,28 @@
opt15->gisprompt = "new,cell,raster";
opt15->guisection = _("Output_options");
- opt16 = G_define_option() ;
- opt16->key = "memory";
- opt16->type = TYPE_INTEGER;
- opt16->required = NO;
- opt16->answer = "300"; /* 300MB default value, please keep in sync with r.terraflow */
- opt16->description = _("Maximum memory to be used with -m flag (in MB)");
+ opt16 = G_define_option();
+ opt16->key = "convergence";
+ opt16->type = TYPE_INTEGER;
+ opt16->required = NO;
+ opt16->answer = "5";
+ opt16->label = _("Convergence factor for MFD (1-10)");
+ opt16->description =
+ _("1 = most diverging flow, 10 = most converging flow. Recommended: 5");
+ opt17 = G_define_option();
+ opt17->key = "memory";
+ opt17->type = TYPE_INTEGER;
+ opt17->required = NO;
+ opt17->answer = "300"; /* 300MB default value, please keep in sync with r.terraflow */
+ opt17->description = _("Maximum memory to be used with -m flag (in MB)");
+
+ flag_sfd = G_define_flag();
+ flag_sfd->key = 's';
+ flag_sfd->label = _("SFD (D8) flow (default is MFD)"); /* copied straight from terraflow */
+ flag_sfd->description =
+ _("SFD: single flow direction, MFD: multiple flow direction");
+
flag_flow = G_define_flag();
flag_flow->key = '4';
flag_flow->description =
@@ -240,6 +258,10 @@
else
strcat(command, "r.watershed.ram");
+ if (flag_sfd->answer) {
+ strcat(command, " -s");
+ }
+
if (flag_flow->answer)
strcat(command, " -4");
@@ -344,7 +366,14 @@
strcat(command, "\"");
}
- if (flag_seg->answer && opt16->answer) {
+ if (!flag_sfd->answer && opt16->answer) {
+ strcat(command, " conv=");
+ strcat(command, "\"");
+ strcat(command, opt16->answer);
+ strcat(command, "\"");
+ }
+
+ if (flag_seg->answer && opt17->answer) {
strcat(command, " mb=");
strcat(command, "\"");
strcat(command, opt16->answer);
@@ -356,7 +385,7 @@
ret = system(command);
- if(ret != EXIT_SUCCESS)
+ if (ret != EXIT_SUCCESS)
G_warning(_("Subprocess failed with exit code %d"), ret);
/* record map metadata/history info */
@@ -373,7 +402,8 @@
"Watershed basins", opt1->answer, flag_seg->answer);
if (opt11->answer)
write_hist(opt11->answer,
- "Watershed stream segments", opt1->answer, flag_seg->answer);
+ "Watershed stream segments", opt1->answer,
+ flag_seg->answer);
if (opt12->answer)
write_hist(opt12->answer,
"Watershed half-basins", opt1->answer, flag_seg->answer);
Modified: grass/trunk/raster/r.watershed/front/r.watershed.html
===================================================================
--- grass/trunk/raster/r.watershed/front/r.watershed.html 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/front/r.watershed.html 2008-12-08 09:11:19 UTC (rev 34789)
@@ -43,6 +43,21 @@
on disk which allows larger maps to be processes but is considerably
slower.
+<dt><em>memory</em>
+
+<dd>Maximum amount of memory in MB to be used with -m set. More memory
+speeds up the processes.
+
+<dt><em>-s</em>
+
+<dd>Use single flow direction (SFD) instead of multiple flow direction (MFD).
+
+<dt><em>convergence</em>
+
+<dd>Convergence factor for MFD. Lower values result in higher divergence,
+flow is more widely distributed. Higher values result in higher convergence,
+flow is less widely distributed, becoming more similar to SFD.
+
<dt><em>-4</em>
<dd>Allow only horizontal and vertical flow of water.
@@ -175,16 +190,19 @@
<h2>NOTES</h2>
-<em>r.watershed</em> uses an algorithm designed to minimize the impact of
-DEM data errors. This algorithm works slower than <em>r.terraflow</em> but
-provides more accurate results in areas of low slope as well as DEMs
-constructed with techniques that mistake canopy tops as the ground elevation.
-Kinner et al. (2005), for example, used SRTM and IFSAR DEMs to compare
-<em>r.watershed</em> against <em>r.terraflow</em> results in Panama.
-<em>r.terraflow</em> was unable to replicate stream locations in the larger
-valleys while <em>r.watershed</em> performed much better. Thus, if forest
-canopy exists in valleys, SRTM, IFSAR, and similar data products will cause
-major errors in <em>r.terraflow</em> stream output. Under similar conditions,
+<h4>A<sup>T</sup> least-cost search algorithm</h4>
+<em>r.watershed</em> uses an A<sup>T</sup> least-cost search algorithm
+(see <a href="#references">REFERENCES</a> section) designed to minimize the
+impact of DEM data errors. This algorithm works slower than
+<em>r.terraflow</em> but provides more accurate results in
+areas of low slope as well as DEMs constructed with techniques that
+mistake canopy tops as the ground elevation. Kinner et al. (2005), for
+example, used SRTM and IFSAR DEMs to compare <em>r.watershed</em>
+against <em>r.terraflow</em> results in Panama. <em>r.terraflow</em> was
+unable to replicate stream locations in the larger valleys while
+<em>r.watershed</em> performed much better. Thus, if forest canopy exists
+in valleys, SRTM, IFSAR, and similar data products will cause major
+errors in <em>r.terraflow</em> stream output. Under similar conditions,
<em>r.watershed</em> will generate better <b>stream</b> and <b>half_basin</b>
results. If watershed divides contain flat to low slope, <em>r.watershed</em>
will generate better basin results than <em>r.terraflow</em>.
@@ -193,53 +211,73 @@
divides contain forest canopy mixed with uncanopied areas using SRTM, IFSAR,
and similar data products, <em>r.watershed</em> will generate better basin
results than <em>r.terraflow</em>.
+The algorithm produces results similar to those obtained when running
+<em><a href="r.cost.html">r.cost</a></em> and
+<em><a href="r.drain.html">r.drain</a></em> on every cell on the map.
-<p>
-There are two versions of this program: <em>ram</em> and <em>seg</em>.
-Which is version is run depends on whether the <em>-m</em> flag is set.
+<h4>Multiple flow direction (MFD)</h4>
+<p><em>r.watershed</em> has experimental support for multiple flow
+direction (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. Although
+basins and half-basins are calculated when using MFD, it is strongly
+recommended to use SFD for basin and half-basin determination.
+<br>See example below with the North Carolina dataset for using MFD mode.
+
+<h4>In-memory mode and disk swap mode</h4>
+<p>There are two versions of this program: <em>ram</em> and <em>seg</em>.
+<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 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.
+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.
<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
out of memory and the resolution size of the current geographic region
cannot be increased, either more memory needs to be added to the computer,
-or the swap space size needs to be increased. If <em>seg</em> runs out of
+or the swap space size needs to be increased. If <em>seg</em> runs out of
memory, additional disk space needs to be freed up for the program to run.
-<p>
-Both versions use the A<sup>T</sup> least-cost search algorithm to determine
-the flow of water over the landscape (see <a href="#seealso">SEE ALSO</a>
-section).
-The algorithm produces results similar to those obtained when running
-<em><a href="r.cost.html">r.cost</a></em> and
-<em><a href="r.drain.html">r.drain</a></em> on every cell on the map.
-
-<p>
-In many situations, the elevation data will be too finely detailed for
-the amount of time or memory available. Running <em>r.watershed</em> may
-require use of a coarser resolution. To make the results more closely
+<h4>Large regions with many cells</h4>
+<p>In some situations, the region size (number of cells) may be too large for
+the amount of time or memory available. Running <em>r.watershed</em> may
+then require use of a coarser resolution. To make the results more closely
resemble the finer terrain data, create a map layer containing the
-lowest elevation values at the coarser resolution. This is done by:
+lowest elevation values at the coarser resolution. This is done by:
1) Setting the current geographic region equal to the elevation map
layer with <em>g.region</em>, and 2) Use the <em>r.neighbors</em> or
<em>r.resamp.stats</em> command to find the lowest value for an area
-equal in size to the desired resolution. For example, if the resolution
+equal in size to the desired resolution. For example, if the resolution
of the elevation data is 30 meters and the resolution of the geographic
-region for <em>r.watershed</em> will be 90 meters: use the minimum
-function for a 3 by 3 neighborhood. After changing to the resolution at
+region for <em>r.watershed</em> will be 90 meters: use the minimum
+function for a 3 by 3 neighborhood. After changing to the resolution at
which <em>r.watershed</em> will be run, <em>r.watershed</em> should be run
using the values from the <em>neighborhood</em> output map layer that
represents the minimum elevation within the region of the coarser cell.
-<p>
-The minimum size of drainage basins, defined by the <em>threshold</em>
+<h4>High-resolution elevation maps with floating point values</h4>
+<p>To get better results with high resolution elevation maps with
+floating point values, it may be necessary to multiply the original map
+with e.g. 100 to change elevation units from meters to cm, because
+<em>r.watershed</em> reads input elevation maps as integer. This allows
+control over how much detail is required to obtain useful results with a
+given dataset and region setting. Remember to adjust any <em>slope length
+and steepness (LS) factor</em> or <em>slope steepness (S) factor</em>
+output accordingly. See example below for using <em>r.mapcalc</em> to
+convert elevation from meter to millimeter.
+
+<h4>Basin threshold</h4>
+<p>The minimum size of drainage basins, defined by the <em>threshold</em>
parameter, is only relevant for those watersheds with a single stream
having at least the <em>threshold</em> of cells flowing into it.
(These watersheds are called exterior basins.)
@@ -248,6 +286,7 @@
an interior stream segment is determined by the distance between the
tributaries flowing into it.
+<h4>MASK and no data</h4>
<p>
The <em>r.watershed</em> program does not require the user to have the
current geographic region filled with elevation values. Areas without
@@ -261,6 +300,7 @@
<p>
Zero data values will be treated as elevation data (not no_data).
+<h4>Further processing of output layers</h4>
<p>
To isolate an individual river network using the output of this module,
a number of approaches may be considered.
@@ -372,6 +412,17 @@
</pre></div>
<br>
+<i>This example uses the North Carolina sample dataset.</i>
+<p>Using MFD mode on a LIDAR elevation dataset:
+<br>Convert elevation values of elev_lid972_1m from meter to millimeter
+then compare MFD mode with default medium flow convergence to SFD mode.
+<div class="code"><pre>
+ r.mapcalc "elev_lid972_1m.1mm = round(elev_lid972_1m * 1000.0)"
+ r.watershed elevation=elev_lid972_1m.1mm accumulation=elev_lid972_1m.1mm.acc drainage=elev_lid972_1m.1mm.dir convergence=5
+ r.watershed -s elevation=elev_lid972_1m.1mm accumulation=elev_lid972_1m.1mm.acc.sfd drainage=elev_lid972_1m.1mm.dir.sfd
+</pre></div>
+
+
<a name="references"></a>
<h2>REFERENCES</h2>
@@ -384,6 +435,12 @@
http://faculty.wiu.edu/CR-Ehlschlaeger2/older/IGIS/paper.html</a>
<p>
+Holmgren, P. (1994). <i>Multiple flow direction algorithms for runoff
+modelling in grid based elevation models: An empirical evaluation.</i>
+<b>Hydrological Processes</b> Vol 8(4), p.327-334.<br>
+DOI: <a href="http://dx.doi.org/10.1002/hyp.3360080405">10.1002/hyp.3360080405</a>
+
+<p>
Kinner D., H. Mitasova, R. Harmon, L. Toma, R., Stallard. (2005).
<i>GIS-based Stream Network Analysis for The Chagres River Basin,
Republic of Panama</i>. <b>The Rio Chagres: A Multidisciplinary Profile of
Modified: grass/trunk/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/Gwater.h 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/Gwater.h 2008-12-08 09:11:19 UTC (rev 34789)
@@ -40,6 +40,7 @@
extern struct Cell_head window;
+extern int mfd, c_fac;
extern int *heap_index, heap_size;
extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
extern SHORT nrows, ncols;
@@ -50,12 +51,13 @@
extern RAMSEG r_h_seg, dep_seg;
extern RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
extern POINT *astar_pts;
-extern CELL *dis, *alt, *wat, *asp, *bas, *haf, *r_h, *dep;
+extern CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
+extern DCELL *wat;
extern CELL *ril_buf;
extern int ril_fd;
extern double *s_l, *s_g, *l_s;
extern CELL one, zero;
-extern double ril_value, dzero;
+extern double ril_value, d_one, d_zero;
extern SHORT sides;
extern SHORT drain[3][3];
extern SHORT updrain[3][3];
@@ -87,11 +89,14 @@
int do_astar(void);
int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
int drop_pt(void);
+int sift_up(int, CELL);
double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
int replace(SHORT, SHORT, SHORT, SHORT);
/* do_cum.c */
int do_cum(void);
+int do_cum_mfd(void);
+double mfd_pow(double, int);
/* find_pour.c */
int find_pourpts(void);
Modified: grass/trunk/raster/r.watershed/ram/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/close_maps.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/close_maps.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -6,30 +6,84 @@
int close_maps(void)
{
- struct Colors colors;
+ struct Colors colors, logcolors;
int value, r, c, fd;
CELL *buf = NULL;
+ DCELL *dbuf = NULL;
+ struct FPRange accRange;
+ DCELL min, max;
+ DCELL clr_min, clr_max;
if (wat_flag || asp_flag || dis_flag || ls_flag || sl_flag || sg_flag)
buf = G_allocate_cell_buf();
+ dbuf = G_allocate_d_raster_buf();
G_free(alt);
if (ls_flag || sg_flag)
G_free(r_h);
if (wat_flag) {
- fd = G_open_cell_new(wat_name);
+ fd = G_open_raster_new(wat_name, DCELL_TYPE);
if (fd < 0) {
G_warning(_("unable to open new accum map layer."));
}
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++) {
- buf[c] = wat[SEG_INDEX(wat_seg, r, c)];
+ dbuf[c] = wat[SEG_INDEX(wat_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."));
+
+ /* set nice color rules: yellow, green, cyan, blue, black */
+ 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;
+ clr_max = max;
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &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);
}
}
@@ -54,6 +108,7 @@
}
G_free(asp);
+ /* visual output no longer needed */
if (dis_flag) {
fd = G_open_cell_new(dis_name);
if (fd < 0) {
Modified: grass/trunk/raster/r.watershed/ram/close_maps2.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/close_maps2.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/close_maps2.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -40,7 +40,7 @@
r = 1;
incr = 0;
while (incr >= 0) {
- G_percent(r, max, 3);
+ G_percent(r, max, 2);
for (gr = 130 + incr; gr <= 255; gr += 20) {
for (rd = 90 + incr; rd <= 255; rd += 30) {
for (bl = 90 + incr; bl <= 255; bl += 40) {
Modified: grass/trunk/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/do_astar.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -3,8 +3,6 @@
#include <grass/gis.h>
#include <grass/glocale.h>
-int sift_up(int, CELL);
-
int do_astar(void)
{
POINT *point;
@@ -96,7 +94,9 @@
}
}
G_percent(count, do_points, 3); /* finish it */
- flag_destroy(worked);
+ if (mfd == 0)
+ flag_destroy(worked);
+
flag_destroy(in_list);
G_free(heap_index);
Modified: grass/trunk/raster/r.watershed/ram/do_astar.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.h 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/do_astar.h 2008-12-08 09:11:19 UTC (rev 34789)
@@ -1,5 +1,5 @@
#ifndef __DO_ASTAR_H__
-#define __DO_ASTAR__
+#define __DO_ASTAR_H__
#define GET_PARENT(c) (int) (((c) - 2) / 3 + 1)
#define GET_CHILD(p) (int) (((p) * 3) - 1)
Modified: grass/trunk/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_cum.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/do_cum.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -9,7 +9,7 @@
CELL is_swale, value, valued;
int killer, threshold, count;
- G_message(_("SECTION 3: Accumulating Surface Flow."));
+ G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
count = 0;
if (bas_thres <= 0)
@@ -17,7 +17,7 @@
else
threshold = bas_thres;
while (first_cum != -1) {
- G_percent(count++, do_points, 3);
+ G_percent(count++, do_points, 2);
killer = first_cum;
first_cum = astar_pts[killer].nxt;
if ((dr = astar_pts[killer].downr) > -1) {
@@ -51,8 +51,258 @@
}
}
}
- G_percent(count, do_points, 3); /* finish it */
+ G_percent(count, do_points, 1); /* finish it */
G_free(astar_pts);
return 0;
}
+
+/***************************************
+ *
+ * MFD references
+ *
+ * original:
+ * Quinn, P., Beven, K., Chevallier, P., and Planchon, 0. 1991.
+ * The prediction of hillslope flow paths for distributed hydrological
+ * modelling using digital terrain models, Hydrol. Process., 5, 59-79.
+ *
+ * modified by Holmgren (1994):
+ * Holmgren, P. 1994. Multiple flow direction algorithms for runoff
+ * modelling in grid based elevation models: an empirical evaluation
+ * Hydrol. Process., 8, 327-334.
+ *
+ * implemented here:
+ * Holmgren (1994) with modifications to honour A * path in order to get
+ * out of depressions and across obstacles with gracefull flow convergence
+ * before depressions/obstacles and gracefull flow divergence after
+ * depressions/obstacles
+ *
+ * ************************************/
+
+int do_cum_mfd(void)
+{
+ int r, c, dr, dc;
+ CELL is_swale;
+ DCELL value, valued, value_sfd;
+ int killer, threshold, count;
+
+ int mfd_cells, astar_not_set;
+ double *dist_to_nbr, *weight, sum_weight, max_weight;
+ int r_nbr, c_nbr, ct_dir, np_side, in_val;
+ double dx, dy;
+ CELL ele, ele_nbr;
+ double prop;
+ int workedon;
+
+ G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
+
+ /* distances to neighbours */
+ dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
+ weight = (double *)G_malloc(sides * sizeof(double));
+
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for neighbours */
+ r_nbr = nextdr[ct_dir];
+ c_nbr = nextdc[ct_dir];
+ /* account for rare cases when ns_res != ew_res */
+ dy = ABS(r_nbr) * window.ns_res;
+ dx = ABS(c_nbr) * window.ew_res;
+ if (ct_dir < 4)
+ dist_to_nbr[ct_dir] = dx + dy;
+ else
+ dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+
+
+ }
+
+ flag_clear_all(worked);
+ workedon = 0;
+
+ count = 0;
+ if (bas_thres <= 0)
+ threshold = 60;
+ else
+ threshold = bas_thres;
+
+ while (first_cum != -1) {
+ 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;
+ value = wat[SEG_INDEX(wat_seg, r, c)];
+ value_sfd = wat[SEG_INDEX(wat_seg, dr, dc)];
+
+ /* disabled for MFD */
+ /* if (ABS(value) > threshold)
+ FLAG_SET(swale, r, c); */
+
+ /* start MFD */
+
+ /* get weights */
+ max_weight = 0;
+ sum_weight = 0;
+ np_side = -1;
+ mfd_cells = 0;
+ astar_not_set = 1;
+ ele = alt[SEG_INDEX(alt_seg, r, c)];
+ /* this loop is needed to get the sum of weights */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for neighbours */
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+ weight[ct_dir] = -1;
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
+ c_nbr < ncols) {
+ in_val = FLAG_GET(worked, r_nbr, c_nbr);
+ if (in_val == 0) {
+ ele_nbr = alt[SEG_INDEX(alt_seg, r_nbr, c_nbr)];
+ /* to consider neighbours with equal ele, change if() */
+ if (ele_nbr < ele) {
+ weight[ct_dir] =
+ mfd_pow(((ele -
+ ele_nbr) / dist_to_nbr[ct_dir]),
+ c_fac);
+ sum_weight += weight[ct_dir];
+ mfd_cells++;
+
+ if (weight[ct_dir] > max_weight)
+ max_weight = weight[ct_dir];
+
+ if (dr == r_nbr && dc == c_nbr) {
+ astar_not_set = 0;
+ }
+ }
+ }
+ if (dr == r_nbr && dc == c_nbr)
+ np_side = ct_dir;
+ }
+ }
+
+ /* 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
+ */
+
+ /* MFD, A * path not included, add to mfd_cells */
+ if (mfd_cells > 0 && astar_not_set == 1) {
+ mfd_cells++;
+ sum_weight += max_weight;
+ 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 */
+ if (mfd_cells > 1) {
+ prop = 0.0;
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
+ c_nbr < ncols && weight[ct_dir] > -0.5) {
+ in_val = FLAG_GET(worked, r_nbr, c_nbr);
+ if (in_val == 0) {
+
+ weight[ct_dir] = weight[ct_dir] / sum_weight;
+ prop += weight[ct_dir]; /* check everything sums up to 1.0 */
+
+ valued = wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)];
+ if (value > 0) {
+ if (valued > 0)
+ valued += value * weight[ct_dir];
+ else
+ valued -= value * weight[ct_dir];
+ }
+ else {
+ if (valued < 0)
+ valued += value * weight[ct_dir];
+ else
+ valued = value * weight[ct_dir] - valued;
+ }
+ wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)] = valued;
+ }
+ else if (ct_dir == np_side) {
+ workedon++; /* check for consistency with A * path */
+ }
+ }
+ }
+ if (ABS(prop - 1.0) > 5E-6f) {
+ G_warning(_("Cumulative proportion not 1.0 but %f"),
+ prop);
+ }
+ valued =
+ ABS(value_sfd) +
+ ABS(value) * weight[np_side] / sum_weight;
+ }
+
+ if (mfd_cells < 2) {
+ valued = value_sfd;
+ if (value > 0) {
+ if (valued > 0)
+ valued += value;
+ else
+ valued -= value;
+ }
+ else {
+ if (valued < 0)
+ valued += value;
+ else
+ valued = value - valued;
+ }
+ wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
+ }
+
+ /* MFD finished */
+
+ valued = ABS(valued) + 0.5;
+
+ is_swale = FLAG_GET(swale, r, c);
+ if (is_swale || (int)valued > threshold) {
+ FLAG_SET(swale, dr, dc);
+ }
+ else {
+ if (er_flag && !is_swale)
+ slope_length(r, c, dr, dc);
+ }
+ FLAG_SET(worked, r, c);
+ }
+ }
+ G_percent(count, do_points, 1); /* finish it */
+ if (workedon)
+ G_message(_("A * path already processed: %d of %d"), workedon,
+ do_points);
+
+ G_free(astar_pts);
+
+ flag_destroy(worked);
+
+ return 0;
+}
+
+double mfd_pow(double base, int exp)
+{
+ int x;
+ double result;
+
+ result = base;
+ if (exp == 1)
+ return result;
+
+ for (x = 2; x <= exp; x++) {
+ result *= base;
+ }
+ return result;
+}
Modified: grass/trunk/raster/r.watershed/ram/find_pour.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/find_pour.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/find_pour.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -8,7 +8,7 @@
basin_num = 0;
for (row = 0; row < nrows; row++) {
- G_percent(row, nrows, 3);
+ G_percent(row, nrows, 1);
northing = window.north - (row + .5) * window.ns_res;
for (col = 0; col < ncols; col++) {
value = FLAG_GET(swale, row, col);
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -1,4 +1,5 @@
#include <stdlib.h>
+#include <string.h>
#include "Gwater.h"
#include <grass/gis.h>
#include <grass/glocale.h>
@@ -15,12 +16,16 @@
ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
ob_flag = bas_flag = seg_flag = haf_flag = arm_flag = dis_flag = 0;
zero = sl_flag = sg_flag = ls_flag = er_flag = bas_thres = 0;
+ one = 1;
nxt_avail_pt = 0;
/* dep_flag = 0; */
- max_length = dzero = 0.0;
+ max_length = d_zero = 0.0;
+ d_one = 1.0;
ril_value = -1.0;
/* dep_slope = 0.0; */
sides = 8;
+ mfd = 1;
+ c_fac = 5;
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
ele_flag++;
@@ -63,9 +68,15 @@
if (sides != 4)
usage(argv[0]);
}
+ else if (sscanf(argv[r], "conv=%d", &c_fac) == 1) ;
+ else if (strcmp(argv[r], "-s") == 0)
+ mfd = 0;
else
usage(argv[0]);
}
+ if (mfd == 1 && (c_fac < 1 || c_fac > 10)) {
+ G_fatal_error("Convergence factor must be between 1 and 10.");
+ }
if ((ele_flag != 1)
||
((arm_flag == 1) &&
@@ -97,7 +108,7 @@
nrows = G_window_rows();
ncols = G_window_cols();
total_cells = nrows * ncols;
- if (max_length <= dzero)
+ if (max_length <= d_zero)
max_length = 10 * nrows * window.ns_res + 10 * ncols * window.ew_res;
if (window.ew_res < window.ns_res)
half_res = .5 * window.ew_res;
@@ -127,7 +138,8 @@
}
G_close_cell(fd);
wat =
- (CELL *) G_malloc(sizeof(CELL) * size_array(&wat_seg, nrows, ncols));
+ (DCELL *) G_malloc(sizeof(DCELL) *
+ size_array(&wat_seg, nrows, ncols));
if (run_flag) {
fd = G_open_cell_old(run_name, "");
@@ -145,10 +157,11 @@
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++)
- wat[SEG_INDEX(wat_seg, r, c)] = 1;
+ wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
}
}
- asp = (CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
+ asp =
+ (CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
if (pit_flag) {
fd = G_open_cell_old(pit_name, "");
Modified: grass/trunk/raster/r.watershed/ram/main.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/main.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/main.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -5,7 +5,7 @@
* AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
* Markus Neteler <neteler itc.it>, Roberto Flor <flor itc.it>,
* Brad Douglas <rez touchofmadness.com>,
- * Hamish Bowman <hamish_b yahoo com>,
+ * Hamish Bowman <hamish_b yahoo com>,
* Markus Metz <markus.metz.giswork gmail.com>
* PURPOSE: Watershed determination
* COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
@@ -24,6 +24,7 @@
struct Cell_head window;
+int mfd, c_fac;
int *heap_index, heap_size;
int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
SHORT nrows, ncols;
@@ -34,24 +35,28 @@
RAMSEG r_h_seg, dep_seg;
RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
POINT *astar_pts;
-CELL *dis, *alt, *wat, *asp, *bas, *haf, *r_h, *dep;
+CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
+DCELL *wat;
CELL *ril_buf;
int ril_fd;
double *s_l, *s_g, *l_s;
CELL one, zero;
-double ril_value, dzero;
+double ril_value, d_one, d_zero;
SHORT sides;
-SHORT drain[3][3] = {{ 7,6,5 },{ 8,0,4 },{ 1,2,3 }};
-SHORT updrain[3][3] = {{ 3,2,1 },{ 4,0,8 },{ 5,6,7 }};
-SHORT nextdr[8] = { 1,-1,0,0,-1,1,1,-1 };
-SHORT nextdc[8] = { 0,0,-1,1,1,-1,1,-1 };
+SHORT drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
+SHORT updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
+SHORT nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+SHORT nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
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];
-char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX], sg_name[GNAME_MAX];
-char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX], dis_name[GNAME_MAX];
+char 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], 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;
@@ -62,7 +67,12 @@
{
init_vars(argc, argv);
do_astar();
- do_cum();
+ if (mfd) {
+ do_cum_mfd();
+ }
+ else {
+ do_cum();
+ }
if (sg_flag || ls_flag) {
sg_factor();
}
Modified: grass/trunk/raster/r.watershed/ram/no_stream.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/no_stream.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/ram/no_stream.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -5,7 +5,8 @@
{
int r, rr, c, cc, uprow = 0, upcol = 0;
double slope;
- CELL downdir, max_drain, value, asp_value, hih_ele, new_ele, aspect;
+ CELL downdir, asp_value, hih_ele, new_ele, aspect;
+ DCELL value, max_drain; /* flow acc is now DCELL */
SHORT updir, riteflag, leftflag, thisdir;
while (1) {
@@ -16,9 +17,10 @@
aspect = asp[SEG_INDEX(asp_seg, r, c)];
if (aspect == drain[rr][cc]) {
value = wat[SEG_INDEX(wat_seg, r, c)];
+ /* here is the reason why MFD basins differ from SFD basins */
if (value < 0)
value = -value;
- if (value > max_drain) {
+ if ((value - max_drain) > 5E-8f) { /* floating point comparison problem workaround */
uprow = r;
upcol = c;
max_drain = value;
Modified: grass/trunk/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/Gwater.h 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/seg/Gwater.h 2008-12-08 09:11:19 UTC (rev 34789)
@@ -43,6 +43,7 @@
extern struct Cell_head window;
+extern int mfd, c_fac;
extern SSEG heap_index;
extern int heap_size;
extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
@@ -51,10 +52,11 @@
extern int bas_thres, tot_parts;
extern SSEG astar_pts;
extern BSEG worked, in_list, s_b, swale;
-extern CSEG dis, alt, wat, asp, bas, haf, r_h, dep;
+extern CSEG dis, alt, asp, bas, haf, r_h, dep;
+extern DSEG wat;
extern DSEG slp, s_l, s_g, l_s, ril;
extern CELL one, zero;
-extern double ril_value, dzero;
+extern double ril_value, d_zero, d_one;
extern SHORT sides;
extern SHORT drain[3][3];
extern SHORT updrain[3][3];
@@ -87,11 +89,14 @@
int do_astar(void);
int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
int drop_pt(void);
+int sift_up(int, CELL);
double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
int replace(SHORT, SHORT, SHORT, SHORT);
/* do_cum.c */
int do_cum(void);
+int do_cum_mfd(void);
+double mfd_pow(double, int);
/* find_pour.c */
int find_pourpts(void);
Modified: grass/trunk/raster/r.watershed/seg/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/close_maps.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/seg/close_maps.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -3,15 +3,68 @@
int close_maps(void)
{
- struct Colors colors;
+ struct Colors colors, logcolors;
int r, c;
CELL is_swale, value;
double dvalue;
+ struct FPRange accRange;
+ DCELL min, max;
+ DCELL clr_min, clr_max;
dseg_close(&slp);
cseg_close(&alt);
- if (wat_flag)
- cseg_write_cellfile(&wat, wat_name);
+ if (wat_flag) {
+ dseg_write_cellfile(&wat, wat_name);
+
+ /* set nice color rules: yellow, green, cyan, blue, black */
+ 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;
+ clr_max = max;
+ G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &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);
+
+ }
if (asp_flag) {
cseg_write_cellfile(&asp, asp_name);
G_init_colors(&colors);
@@ -19,26 +72,27 @@
G_write_colors(asp_name, this_mapset, &colors);
}
cseg_close(&asp);
+ /* visual ouput no longer needed */
if (dis_flag) {
if (bas_thres <= 0)
bas_thres = 60;
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++) {
- cseg_get(&wat, &value, r, c);
- if (value < 0) {
- value = 0;
- cseg_put(&wat, &value, r, c);
+ dseg_get(&wat, &dvalue, r, c);
+ if (dvalue < 0) {
+ dvalue = 0;
+ dseg_put(&wat, &dvalue, r, c);
}
else {
bseg_get(&swale, &is_swale, r, c);
if (is_swale) {
- value = bas_thres;
- cseg_put(&wat, &value, r, c);
+ dvalue = bas_thres;
+ dseg_put(&wat, &dvalue, r, c);
}
}
}
}
- cseg_write_cellfile(&wat, dis_name);
+ dseg_write_cellfile(&wat, dis_name);
G_init_colors(&colors);
G_make_rainbow_colors(&colors, 1, 120);
G_write_colors(dis_name, this_mapset, &colors);
Modified: grass/trunk/raster/r.watershed/seg/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_astar.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/seg/do_astar.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -5,14 +5,13 @@
#include "Gwater.h"
#include "do_astar.h"
-int sift_up(int, CELL);
-
int do_astar(void)
{
POINT point;
int doer, count;
SHORT upr, upc, r, c, ct_dir;
- CELL work_val, alt_val, alt_up, asp_up, wat_val;
+ CELL work_val, alt_val, alt_up, asp_up;
+ DCELL wat_val;
CELL in_val, drain_val;
HEAP heap_pos;
@@ -75,10 +74,10 @@
if (asp_up < -1) {
drain_val = drain[upr - r + 1][upc - c + 1];
cseg_put(&asp, &drain_val, upr, upc);
- cseg_get(&wat, &wat_val, r, c);
+ dseg_get(&wat, &wat_val, r, c);
if (wat_val > 0)
wat_val = -wat_val;
- cseg_put(&wat, &wat_val, r, c);
+ 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);
@@ -89,7 +88,10 @@
}
}
}
- bseg_close(&worked);
+
+ if (mfd == 0)
+ bseg_close(&worked);
+
bseg_close(&in_list);
seg_close(&heap_index);
@@ -288,6 +290,7 @@
heap_run = 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;
Modified: grass/trunk/raster/r.watershed/seg/do_astar.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_astar.h 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/seg/do_astar.h 2008-12-08 09:11:19 UTC (rev 34789)
@@ -1,5 +1,5 @@
#ifndef __DO_ASTAR_H__
-#define __DO_ASTAR__
+#define __DO_ASTAR_H__
#define GET_PARENT(c) ((int) (((c) - 2) / 3 + 1))
#define GET_CHILD(p) ((int) ((p) * 3 - 1))
Modified: grass/trunk/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_cum.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/seg/do_cum.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -7,11 +7,12 @@
int do_cum(void)
{
SHORT r, c, dr, dc;
- CELL is_swale, value, valued;
+ CELL is_swale;
+ DCELL value, valued;
POINT point;
int killer, threshold, count;
- G_message(_("SECTION 3: Accumulating Surface Flow."));
+ G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
count = 0;
if (bas_thres <= 0)
@@ -19,7 +20,7 @@
else
threshold = bas_thres;
while (first_cum != -1) {
- G_percent(count++, do_points, 3);
+ G_percent(count++, do_points, 2);
killer = first_cum;
seg_get(&astar_pts, (char *)&point, 0, killer);
first_cum = point.nxt;
@@ -27,10 +28,10 @@
r = point.r;
c = point.c;
dc = point.downc;
- cseg_get(&wat, &value, r, c);
+ dseg_get(&wat, &value, r, c);
if (ABS(value) >= threshold)
bseg_put(&swale, &one, r, c);
- cseg_get(&wat, &valued, dr, dc);
+ dseg_get(&wat, &valued, dr, dc);
if (value > 0) {
if (valued > 0)
valued += value;
@@ -43,7 +44,7 @@
else
valued = value - valued;
}
- cseg_put(&wat, &valued, dr, dc);
+ dseg_put(&wat, &valued, dr, dc);
bseg_get(&swale, &is_swale, r, c);
if (is_swale || ABS(valued) >= threshold) {
bseg_put(&swale, &one, dr, dc);
@@ -56,6 +57,264 @@
}
seg_close(&astar_pts);
- G_percent(count, do_points, 3); /* finish it */
+ G_percent(count, do_points, 1); /* finish it */
return 0;
}
+
+/***************************************
+ *
+ * MFD references
+ *
+ * original:
+ * Quinn, P., Beven, K., Chevallier, P., and Planchon, 0. 1991.
+ * The prediction of hillslope flow paths for distributed hydrological
+ * modelling using digital terrain models, Hydrol. Process., 5, 59-79.
+ *
+ * modified by Holmgren (1994):
+ * Holmgren, P. 1994. Multiple flow direction algorithms for runoff
+ * modelling in grid based elevation models: an empirical evaluation
+ * Hydrol. Process., 8, 327-334.
+ *
+ * implemented here:
+ * Holmgren (1994) with modifications to honour A * path in order to get
+ * out of depressions and across obstacles with gracefull flow convergence
+ * before depressions/obstacles and gracefull flow divergence after
+ * depressions/obstacles
+ *
+ * ************************************/
+
+int do_cum_mfd(void)
+{
+ int r, c, dr, dc;
+ CELL is_swale;
+ DCELL value, valued, value_sfd;
+ POINT point;
+ int killer, threshold, count;
+
+ int mfd_cells, astar_not_set;
+ double *dist_to_nbr, *weight, sum_weight, max_weight;
+ int r_nbr, c_nbr, ct_dir, np_side, in_val;
+ double dx, dy;
+ CELL ele, ele_nbr;
+ double prop;
+ int workedon;
+
+ G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
+
+ /* distances to neighbours */
+ dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
+ weight = (double *)G_malloc(sides * sizeof(double));
+
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for neighbours */
+ r_nbr = nextdr[ct_dir];
+ c_nbr = nextdc[ct_dir];
+ /* account for rare cases when ns_res != ew_res */
+ dy = ABS(r_nbr) * window.ns_res;
+ dx = ABS(c_nbr) * window.ew_res;
+ if (ct_dir < 4)
+ dist_to_nbr[ct_dir] = dx + dy;
+ else
+ dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+
+
+ }
+
+ /* reset worked, takes time... */
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ bseg_put(&worked, &zero, r, c);
+ }
+ }
+
+ workedon = 0;
+
+ count = 0;
+ if (bas_thres <= 0)
+ threshold = 60;
+ else
+ threshold = bas_thres;
+
+ while (first_cum != -1) {
+ G_percent(count++, do_points, 2);
+ 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;
+ dseg_get(&wat, &value, r, c);
+ dseg_get(&wat, &value_sfd, dr, dc);
+
+ /* disabled for MFD */
+ /* if (ABS(value) > threshold)
+ bseg_put(&swale, &one, r, c); */
+
+ /* start MFD */
+
+ /* get weights */
+ max_weight = 0;
+ sum_weight = 0;
+ np_side = -1;
+ mfd_cells = 0;
+ astar_not_set = 1;
+ cseg_get(&alt, &ele, r, c);
+ /* this loop is needed to get the sum of weights */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for neighbours */
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+ weight[ct_dir] = -1;
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
+ c_nbr < ncols) {
+ bseg_get(&worked, &in_val, r_nbr, c_nbr);
+ if (in_val == 0) {
+ cseg_get(&alt, &ele_nbr, r_nbr, c_nbr);
+ /* to consider neighbours with equal ele, change if() */
+ if (ele_nbr < ele) {
+ weight[ct_dir] =
+ mfd_pow(((ele -
+ ele_nbr) / dist_to_nbr[ct_dir]),
+ c_fac);
+ sum_weight += weight[ct_dir];
+ mfd_cells++;
+
+ if (weight[ct_dir] > max_weight)
+ max_weight = weight[ct_dir];
+
+ if (dr == r_nbr && dc == c_nbr) {
+ astar_not_set = 0;
+ }
+ }
+ }
+ if (dr == r_nbr && dc == c_nbr)
+ np_side = ct_dir;
+ }
+ }
+
+ /* 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
+ */
+
+ /* MFD, A * path not included, add to mfd_cells */
+ if (mfd_cells > 0 && astar_not_set == 1) {
+ mfd_cells++;
+ sum_weight += max_weight;
+ 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 */
+ if (mfd_cells > 1) {
+ prop = 0.0;
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+
+ /* 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, &in_val, r_nbr, c_nbr);
+ if (in_val == 0) {
+
+ weight[ct_dir] = weight[ct_dir] / sum_weight;
+ prop += weight[ct_dir]; /* check everything sums up to 1.0 */
+
+ dseg_get(&wat, &valued, r_nbr, c_nbr);
+ if (value > 0) {
+ if (valued > 0)
+ valued += value * weight[ct_dir];
+ else
+ valued -= value * weight[ct_dir];
+ }
+ else {
+ if (valued < 0)
+ valued += value * weight[ct_dir];
+ else
+ valued = value * weight[ct_dir] - valued;
+ }
+ dseg_put(&wat, &valued, r_nbr, c_nbr);
+ }
+ else if (ct_dir == np_side) {
+ workedon++; /* check for consistency with A * path */
+ }
+ }
+ }
+ if (ABS(prop - 1.0) > 5E-6f) {
+ G_warning(_("3 cumulative proportion not 1.0 but %f"),
+ prop);
+ }
+ valued =
+ ABS(value_sfd) +
+ ABS(value) * weight[np_side] / sum_weight;
+ }
+
+ if (mfd_cells < 2) {
+ valued = value_sfd;
+ if (value > 0) {
+ if (valued > 0)
+ valued += value;
+ else
+ valued -= value;
+ }
+ else {
+ if (valued < 0)
+ valued += value;
+ else
+ valued = value - valued;
+ }
+ dseg_put(&wat, &valued, dr, dc);
+ }
+
+ /* MFD finished */
+
+ valued = ABS(valued) + 0.5;
+
+ bseg_get(&swale, &is_swale, r, c);
+ if (is_swale || (int)valued > threshold) {
+ bseg_put(&swale, &one, dr, dc);
+ }
+ else {
+ if (er_flag && !is_swale)
+ slope_length(r, c, dr, dc);
+ }
+ bseg_put(&worked, &one, r, c);
+ }
+ }
+ G_percent(count, do_points, 1); /* finish it */
+ if (workedon)
+ G_message(_("A * path already processed: %d of %d"), workedon,
+ do_points);
+
+ seg_close(&astar_pts);
+
+ bseg_close(&worked);
+
+ return 0;
+}
+
+double mfd_pow(double base, int exp)
+{
+ int x;
+ double result;
+
+ result = base;
+ if (exp == 1)
+ return result;
+
+ for (x = 2; x <= exp; x++) {
+ result *= base;
+ }
+ return result;
+}
Modified: grass/trunk/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/init_vars.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/seg/init_vars.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -1,4 +1,5 @@
#include <stdlib.h>
+#include <string.h>
#include <unistd.h>
#include "Gwater.h"
#include <grass/gis.h>
@@ -14,7 +15,8 @@
/* int page_block, num_cseg; */
int max_bytes;
- CELL *buf, alt_value, wat_value, asp_value, worked_value;
+ CELL *buf, alt_value, asp_value, worked_value;
+ DCELL wat_value;
char MASK_flag;
G_gisinit(argv[0]);
@@ -23,11 +25,14 @@
zero = sl_flag = sg_flag = ls_flag = er_flag = bas_thres = 0;
nxt_avail_pt = 0;
/* dep_flag = 0; */
- max_length = dzero = 0.0;
+ max_length = d_zero = 0.0;
+ d_one = 1.0;
ril_value = -1.0;
/* dep_slope = 0.0; */
max_bytes = 0;
sides = 8;
+ mfd = 1;
+ c_fac = 5;
segs_mb = 300;
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
@@ -72,9 +77,15 @@
if (sides != 4)
usage(argv[0]);
}
+ else if (sscanf(argv[r], "conv=%d", &c_fac) == 1) ;
+ else if (strcmp(argv[r], "-s") == 0)
+ mfd = 0;
else
usage(argv[0]);
}
+ if (mfd == 1 && (c_fac < 1 || c_fac > 10)) {
+ G_fatal_error("Convergence factor must be between 1 and 10.");
+ }
if ((ele_flag != 1)
||
((arm_flag == 1) &&
@@ -107,7 +118,7 @@
G_get_set_window(&window);
nrows = G_window_rows();
ncols = G_window_cols();
- if (max_length <= dzero)
+ if (max_length <= d_zero)
max_length = 10 * nrows * window.ns_res + 10 * ncols * window.ew_res;
if (window.ew_res < window.ns_res)
half_res = .5 * window.ew_res;
@@ -158,15 +169,15 @@
cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
cseg_read_cell(&alt, ele_name, "");
cseg_read_cell(&r_h, ele_name, "");
- cseg_open(&wat, seg_rows, seg_cols, num_open_segs);
+ dseg_open(&wat, seg_rows, seg_cols, num_open_segs);
if (run_flag) {
- cseg_read_cell(&wat, run_name, "");
+ dseg_read_cell(&wat, run_name, "");
}
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++)
- if (-1 == cseg_put(&wat, &one, r, c))
+ if (-1 == dseg_put(&wat, &d_one, r, c))
exit(EXIT_FAILURE);
}
}
@@ -243,21 +254,21 @@
if (MASK_flag) {
for (r = 0; r < nrows; r++) {
- G_percent(r, nrows, 3);
+ G_percent(r, nrows, 2);
for (c = 0; c < ncols; c++) {
bseg_get(&worked, &worked_value, r, c);
if (worked_value) {
- cseg_put(&wat, &zero, r, c);
+ dseg_put(&wat, &d_zero, r, c);
}
else {
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) {
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- cseg_put(&wat, &wat_value, r, c);
+ dseg_put(&wat, &wat_value, r, c);
}
if (r == 0)
asp_value = -2;
@@ -280,10 +291,10 @@
add_pt(r, c, -1, -1, alt_value, alt_value);
asp_value = -2;
cseg_put(&asp, &asp_value, r, c);
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- cseg_put(&wat, &wat_value, r, c);
+ dseg_put(&wat, &wat_value, r, c);
}
}
else if (!bseg_get(&worked, &worked_value, r + 1, c)
@@ -292,10 +303,10 @@
add_pt(r, c, -1, -1, alt_value, alt_value);
asp_value = -6;
cseg_put(&asp, &asp_value, r, c);
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- cseg_put(&wat, &wat_value, r, c);
+ dseg_put(&wat, &wat_value, r, c);
}
}
else if (!bseg_get(&worked, &worked_value, r, c - 1)
@@ -304,10 +315,10 @@
add_pt(r, c, -1, -1, alt_value, alt_value);
asp_value = -4;
cseg_put(&asp, &asp_value, r, c);
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- cseg_put(&wat, &wat_value, r, c);
+ dseg_put(&wat, &wat_value, r, c);
}
}
else if (!bseg_get(&worked, &worked_value, r, c + 1)
@@ -316,10 +327,10 @@
add_pt(r, c, -1, -1, alt_value, alt_value);
asp_value = -8;
cseg_put(&asp, &asp_value, r, c);
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- cseg_put(&wat, &wat_value, r, c);
+ dseg_put(&wat, &wat_value, r, c);
}
}
else if (sides == 8 &&
@@ -329,10 +340,10 @@
add_pt(r, c, -1, -1, alt_value, alt_value);
asp_value = -3;
cseg_put(&asp, &asp_value, r, c);
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- cseg_put(&wat, &wat_value, r, c);
+ dseg_put(&wat, &wat_value, r, c);
}
}
else if (sides == 8 &&
@@ -342,10 +353,10 @@
add_pt(r, c, -1, -1, alt_value, alt_value);
asp_value = -1;
cseg_put(&asp, &asp_value, r, c);
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- cseg_put(&wat, &wat_value, r, c);
+ dseg_put(&wat, &wat_value, r, c);
}
}
else if (sides == 8 &&
@@ -355,10 +366,10 @@
add_pt(r, c, -1, -1, alt_value, alt_value);
asp_value = -5;
cseg_put(&asp, &asp_value, r, c);
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- cseg_put(&wat, &wat_value, r, c);
+ dseg_put(&wat, &wat_value, r, c);
}
}
else if (sides == 8 &&
@@ -368,10 +379,10 @@
add_pt(r, c, -1, -1, alt_value, alt_value);
asp_value = -7;
cseg_put(&asp, &asp_value, r, c);
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- cseg_put(&wat, &wat_value, r, c);
+ dseg_put(&wat, &wat_value, r, c);
}
}
else {
@@ -384,17 +395,17 @@
}
else {
for (r = 0; r < nrows; r++) {
- G_percent(r, nrows, 3);
+ 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);
cseg_get(&asp, &asp_value, r, c);
if (r == 0 || c == 0 || r == nrows - 1 ||
c == ncols - 1 || asp_value != 0) {
- cseg_get(&wat, &wat_value, r, c);
+ dseg_get(&wat, &wat_value, r, c);
if (wat_value > 0) {
wat_value = -wat_value;
- if (-1 == cseg_put(&wat, &wat_value, r, c))
+ if (-1 == dseg_put(&wat, &wat_value, r, c))
exit(EXIT_FAILURE);
}
if (r == 0)
@@ -419,8 +430,7 @@
}
}
}
- G_percent(r, nrows, 3); /* finish it */
+ G_percent(r, nrows, 1); /* finish it */
return 0;
}
-
Modified: grass/trunk/raster/r.watershed/seg/main.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/main.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/seg/main.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -25,6 +25,7 @@
struct Cell_head window;
+int mfd, c_fac;
SSEG heap_index;
int heap_size;
int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
@@ -33,15 +34,16 @@
int bas_thres, tot_parts;
SSEG astar_pts;
BSEG worked, in_list, s_b, swale;
-CSEG dis, alt, wat, asp, bas, haf, r_h, dep;
+CSEG dis, alt, asp, bas, haf, r_h, dep;
+DSEG wat;
DSEG slp, s_l, s_g, l_s, ril;
CELL one, zero;
-double ril_value, dzero;
+double ril_value, d_zero, d_one;
SHORT sides;
-SHORT drain[3][3] = {{ 7,6,5 },{ 8,0,4 },{ 1,2,3 }};
-SHORT updrain[3][3] = {{ 3,2,1 },{ 4,0,8 },{ 5,6,7 }};
-SHORT nextdr[8] = { 1,-1,0,0,-1,1,1,-1 };
-SHORT nextdc[8] = { 0,0,-1,1,1,-1,1,-1 };
+SHORT drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
+SHORT updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
+SHORT nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+SHORT nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
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];
@@ -65,10 +67,16 @@
{
one = 1;
zero = 0;
- dzero = 0.0;
+ d_zero = 0.0;
+ d_one = 1.0;
init_vars(argc, argv);
do_astar();
- do_cum();
+ if (mfd) {
+ do_cum_mfd();
+ }
+ else {
+ do_cum();
+ }
if (sg_flag || ls_flag) {
sg_factor();
}
Modified: grass/trunk/raster/r.watershed/seg/no_stream.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/no_stream.c 2008-12-08 06:36:09 UTC (rev 34788)
+++ grass/trunk/raster/r.watershed/seg/no_stream.c 2008-12-08 09:11:19 UTC (rev 34789)
@@ -6,7 +6,8 @@
{
int r, rr, c, cc, uprow = 0, upcol = 0;
double slope;
- CELL downdir, max_drain, value, asp_value, hih_ele, new_ele, aspect;
+ CELL downdir, asp_value, hih_ele, new_ele, aspect;
+ DCELL value, max_drain; /* flow acc is now DCELL */
SHORT updir, riteflag, leftflag, thisdir;
while (1) {
@@ -17,10 +18,11 @@
cseg_get(&asp, &aspect, r, c);
if (aspect == drain[rr][cc]) {
- cseg_get(&wat, &value, r, c);
+ dseg_get(&wat, &value, r, c);
+ /* here is the reason why MFD basins differ from SFD basins */
if (value < 0)
value = -value;
- if (value > max_drain) {
+ if ((value - max_drain) > 5E-8f) { /* floating point comparison problem workaround */
uprow = r;
upcol = c;
max_drain = value;
More information about the grass-commit
mailing list