[GRASS-SVN] r42121 - in
grass/branches/releasebranch_6_4/raster/r.watershed: front ram seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 5 09:30:34 EDT 2010
Author: mmetz
Date: 2010-05-05 09:30:34 -0400 (Wed, 05 May 2010)
New Revision: 42121
Modified:
grass/branches/releasebranch_6_4/raster/r.watershed/front/description.html
grass/branches/releasebranch_6_4/raster/r.watershed/front/main.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/Gwater.h
grass/branches/releasebranch_6_4/raster/r.watershed/ram/Makefile
grass/branches/releasebranch_6_4/raster/r.watershed/ram/close_maps.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/close_maps2.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_astar.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_astar.h
grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_cum.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/find_pour.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/init_vars.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/main.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/no_stream.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/ramseg.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/sg_factor.c
grass/branches/releasebranch_6_4/raster/r.watershed/ram/slope_len.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/Gwater.h
grass/branches/releasebranch_6_4/raster/r.watershed/seg/close_maps.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/close_maps2.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_astar.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_astar.h
grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_cum.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_get.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_put.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_read.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_write.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/init_vars.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/main.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/no_stream.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/sg_factor.c
grass/branches/releasebranch_6_4/raster/r.watershed/seg/slope_len.c
Log:
replace faulty 64 version
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/front/description.html
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/front/description.html 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/front/description.html 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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,10 +24,15 @@
<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>-f</em>
+
+<dd>Use multiple flow direction (MFD) instead of single flow direction (SFD, D8).
+SFD is enabled by default.
+
<dt><em>-4</em>
<dd>Allow only horizontal and vertical flow of water.
@@ -35,15 +40,38 @@
surface flow (allows horizontal, vertical, and diagonal flow of water).
This flag will also make the drainage basins look more homogeneous.
+<dt><em>-a</em>
+
+<dd>Use positive flow accumulation even for likely underestimates. When this
+flag is not set, cells with a flow accumulation value that is likely to be
+an underestimate are converted to the negative. See below for a detailed
+description of flow accumulation output.
+
+<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>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>elevation</em>
-<dd>Input map: Elevation on which entire analysis is based.
+<dd>Input map: Elevation on which entire analysis is based. NULL (nodata)
+cells are ignored, zero and negative values are valid elevation data.
+Gaps in the elevation map that are located within the area of interest
+must be filled beforehand, e.g. with <em>r.fillnulls</em>, to avoid
+distortions.
<dt><em>depression</em>
<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. All cells that are not NULL and not zero indicate
+depressions. Water will flow into but not out of depressions.
<dt><em>flow</em>
@@ -51,7 +79,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>
@@ -64,14 +92,15 @@
<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. All cells that are not NULL and not zero indicate blocking
+terrain.
<dt><em>threshold</em>
<dd>The minimum size of an exterior watershed basin in cells, if no flow
map is input, or overland flow units when a flow map is given.
Warning: low threshold values will dramactically increase run time and
-generate difficult too read basin and half.basin results.
+generate difficult to read basin and half_basin results.
This parameter also controls the level of detail in the <em>stream</em>
segments map.
@@ -98,10 +127,10 @@
<dt><em>drainage</em>
<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
+cell measured CCW from East. Multiplying positive values by 45 will give
+the direction in degrees that the surface runoff will travel from that
+cell. The 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.
@@ -117,7 +146,8 @@
<dt><em>stream</em>
<dd>Output map: stream segments. Values correspond to the watershed
-basin values.
+basin values. Can be vectorized after thinning (<em>r.thin</em>) with
+<em>r.to.vect</em>.
<dt><em>half.basin</em>
@@ -130,38 +160,33 @@
<dt><em>visual</em>
-<dd>Output map: useful for visual display of results.
-Surface runoff accumulation with the values
-modified to provide for easy display. All negative accumulation values
-are changed to zero. All positive values above the basin threshold
-are given the value of the <em>threshold</em> parameter.
+<dd><b>DEPRECATED</b> A colortable is generated by default for the
+<em>accumulation</em> output for easy visual inspection.
<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), it is multiplied by 100.
<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 S factor is a small
+number (usually less than one), it is multiplied by 100.
</dd>
</dl>
<h2>NOTES</h2>
-<em>r.watershed</em> uses an A<sup>T</sup> least-cost search algorithm
-(see <a href="#references">REFERENCES</a> section) designed to minimize
+<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. Compared to <em>r.terraflow</em>, this
algorithm provides more accurate results in areas of low slope as well
as DEMs constructed with techniques that mistake canopy tops as the
@@ -180,52 +205,69 @@
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>
+<h4>Multiple flow direction (MFD)</h4>
+
+<em>r.watershed</em> offers two methods to calculate surface flow:
+single flow direction (SFD, D8) and multiple flow direction (MFD). With
+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. As a result, depressions and obstacles are traversed
+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.
+
+<h4>In-memory mode and disk swap mode</h4>
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.
+<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.
+version; <em>seg</em> uses the GRASS segmentation library which manages
+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
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
+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
memory, additional disk space needs to be freed up for the program to run.
+The <em>r.terraflow</em> module was specifically designed with huge
+regions in mind and may be useful here as an alternative.
-<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>
+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>
+<h4>Basin threshold</h4>
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.
@@ -235,19 +277,26 @@
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
-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).
+Gaps (NULL cells) in the elevation map that are located within the area
+of interest will heavily influence the analysis: water will
+flow into but not out of these gaps. These gaps must be filled beforehand,
+e.g. with <em>r.fillnulls</em>.
+<p>
+Zero (0) and negative 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.
@@ -261,12 +310,25 @@
starting point.
</ol>
+All individual river networks in the stream segments output can be
+identified through their ultimate outlet points. These points are all
+cells in the stream segments output with negative drainage direction.
+These points can be used as start points for <em>r.water.outlet</em> or
+<em>v.net.iso</em>.
<p>
To create <i>river mile</i> segmentation from a vectorized streams map,
try the <em>v.net.iso</em> or <em>v.lrs.segment</em> modules.
+<p>
+The stream segments output can be easily vectorized after thinning with
+<em>r.thin</em>. Each stream segment in the vector map will have the
+value of the associated basin. To isolate subbasins and streams for a
+larger basin, a MASK for the larger basin can be created with
+<em>r.water.outlet</em>. The stream segments output serves as a guide
+where to place the outlet point used as input to <em>r.water.outlet</em>.
+The basin threshold must have been sufficiently small to isolate a
+stream network and subbasins within the larger basin.
-
<h2>EXAMPLES</h2>
<i>These examples use the Spearfish sample dataset.</i>
<p>
@@ -284,7 +346,7 @@
<br>
<p>
-Set a nice color table for the accumulation map:
+Set a different color table for the accumulation map:
<div class="code"><pre>
MAP=rwater.accum
r.watershed elev=elevation.dem accum=$MAP
@@ -315,7 +377,7 @@
it to a vector output map. The accumulation cut-off, and therefore fractal
dimension, is arbitrary; in this example we use the map's mean number of
upstream catchment cells (calculated in the above example by <em>r.univar</em>)
-as the cut-off value.
+as the cut-off value. This only works with SFD, not with MFD.
<div class="code"><pre>
r.watershed elev=elevation.dem accum=rwater.accum
@@ -371,6 +433,12 @@
http://chuck.ehlschlaeger.info/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
@@ -395,8 +463,10 @@
<a href="g.region.html">g.region</a>,
<a href="r.cost.html">r.cost</a>,
<a href="r.drain.html">r.drain</a>,
+<a href="r.fillnulls.html">r.fillnulls</a>,
<a href="r.flow.html">r.flow</a>,
<!-- <a href="r.flowmd.html">r.flowmd</a>, -->
+<a href="r.mask.html">r.mask</a>,
<a href="r.neighbors.html">r.neighbors</a>,
<a href="r.param.scale.html">r.param.scale</a>,
<a href="r.resamp.interp.html">r.resamp.interp</a>,
@@ -411,7 +481,8 @@
Original version:
Charles Ehlschlaeger, U.S. Army Construction Engineering Research Laboratory
<br>
-Speedups: Markus Metz <markus.metz.giswork at gmail.com>
+Faster sorting algorithm and MFD support:
+Markus Metz <markus.metz.giswork at gmail.com>
<p>
<i>Last changed: $Date$</i>
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/front/main.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/front/main.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/front/main.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -4,9 +4,10 @@
* MODULE: front end
* AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
* Brad Douglas <rez touchofmadness.com>,
- * Hamish Bowman <hamish_b yahoo.com>
- * PURPOSE: Watershed determination
- * COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
+ * Hamish Bowman <hamish_b yahoo.com>
+ * Markus Metz <markus.metz.giswork gmail.com>
+ * PURPOSE: Hydrological analysis
+ * COPYRIGHT: (C) 1999-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -19,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[])
{
@@ -41,8 +42,11 @@
struct Option *opt14;
struct Option *opt15;
struct Option *opt16;
+ struct Option *opt17;
+ struct Flag *flag_mfd;
struct Flag *flag_flow;
struct Flag *flag_seg;
+ struct Flag *flag_abs;
struct GModule *module;
G_gisinit(argv[0]);
@@ -57,20 +61,16 @@
_("Input map: elevation on which entire analysis is based");
opt1->guisection = _("Input_options");
- opt2 = G_define_option();
+ opt2 = G_define_standard_option(G_OPT_R_INPUT);
opt2->key = "depression";
opt2->description = _("Input map: locations of real depressions");
opt2->required = NO;
- opt2->type = TYPE_STRING;
- opt2->gisprompt = "old,cell,raster";
opt2->guisection = _("Input_options");
- opt3 = G_define_option();
+ opt3 = G_define_standard_option(G_OPT_R_INPUT);
opt3->key = "flow";
opt3->description = _("Input map: amount of overland flow per cell");
opt3->required = NO;
- opt3->type = TYPE_STRING;
- opt3->gisprompt = "old,cell,raster";
opt3->guisection = _("Input_options");
opt4 = G_define_option();
@@ -82,109 +82,104 @@
opt4->gisprompt = "old,cell,raster";
opt4->guisection = _("Input_options");
- opt5 = G_define_option();
+ opt5 = G_define_standard_option(G_OPT_R_INPUT);
opt5->key = "blocking";
opt5->description =
_("Input map: terrain blocking overland surface flow, for USLE");
opt5->required = NO;
- opt5->type = TYPE_STRING;
- opt5->gisprompt = "old,cell,raster";
opt5->guisection = _("Input_options");
- opt6 = G_define_option();
- opt6->key = "threshold";
- opt6->description =
- _("Input value: minimum size of exterior watershed basin");
- opt6->required = NO;
- opt6->type = TYPE_INTEGER;
- opt6->gisprompt = "new,cell,raster";
- opt6->guisection = _("Input_options");
-
- opt7 = G_define_option();
- opt7->key = "max.slope.length";
- opt7->description =
- _("Input value: maximum length of surface flow, for USLE");
- opt7->required = NO;
- opt7->type = TYPE_DOUBLE;
- opt7->gisprompt = "new,cell,raster";
- opt7->guisection = _("Input_options");
-
- opt8 = G_define_option();
+ opt8 = G_define_standard_option(G_OPT_R_OUTPUT);
opt8->key = "accumulation";
opt8->description =
_("Output map: number of cells that drain through each cell");
opt8->required = NO;
- opt8->type = TYPE_STRING;
- opt8->gisprompt = "new,cell,raster";
opt8->guisection = _("Output_options");
- opt9 = G_define_option();
+ opt9 = G_define_standard_option(G_OPT_R_OUTPUT);
opt9->key = "drainage";
opt9->description = _("Output map: drainage direction");
opt9->required = NO;
- opt9->type = TYPE_STRING;
- opt9->gisprompt = "new,cell,raster";
opt9->guisection = _("Output_options");
- opt10 = G_define_option();
+ opt10 = G_define_standard_option(G_OPT_R_OUTPUT);
opt10->key = "basin";
opt10->description =
_("Output map: unique label for each watershed basin");
opt10->required = NO;
- opt10->type = TYPE_STRING;
- opt10->gisprompt = "new,cell,raster";
opt10->guisection = _("Output_options");
- opt11 = G_define_option();
+ opt11 = G_define_standard_option(G_OPT_R_OUTPUT);
opt11->key = "stream";
opt11->description = _("Output map: stream segments");
opt11->required = NO;
- opt11->type = TYPE_STRING;
- opt11->gisprompt = "new,cell,raster";
opt11->guisection = _("Output_options");
- opt12 = G_define_option();
+ opt12 = G_define_standard_option(G_OPT_R_OUTPUT);
opt12->key = "half.basin";
opt12->description =
_("Output map: each half-basin is given a unique value");
opt12->required = NO;
- opt12->type = TYPE_STRING;
- opt12->gisprompt = "new,cell,raster";
opt12->guisection = _("Output_options");
- opt13 = G_define_option();
+ opt13 = G_define_standard_option(G_OPT_R_OUTPUT);
opt13->key = "visual";
opt13->description =
_("Output map: useful for visual display of results");
opt13->required = NO;
- opt13->type = TYPE_STRING;
- opt13->gisprompt = "new,cell,raster";
opt13->guisection = _("Output_options");
- opt14 = G_define_option();
+ opt14 = G_define_standard_option(G_OPT_R_OUTPUT);
opt14->key = "length.slope";
opt14->description =
_("Output map: slope length and steepness (LS) factor for USLE");
opt14->required = NO;
- opt14->type = TYPE_STRING;
- opt14->gisprompt = "new,cell,raster";
opt14->guisection = _("Output_options");
- opt15 = G_define_option();
+ opt15 = G_define_standard_option(G_OPT_R_OUTPUT);
opt15->key = "slope.steepness";
opt15->description = _("Output map: slope steepness (S) factor for USLE");
opt15->required = NO;
- opt15->type = TYPE_STRING;
- 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)");
+ opt6 = G_define_option();
+ opt6->key = "threshold";
+ opt6->description =
+ _("Input value: minimum size of exterior watershed basin");
+ opt6->required = NO;
+ opt6->type = TYPE_INTEGER;
+ opt6->guisection = _("Input_options");
+ opt7 = G_define_option();
+ opt7->key = "max.slope.length";
+ opt7->description =
+ _("Input value: maximum length of surface flow, for USLE");
+ opt7->required = NO;
+ opt7->type = TYPE_DOUBLE;
+ opt7->guisection = _("Input_options");
+
+ 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_mfd = G_define_flag();
+ flag_mfd->key = 'f';
+ flag_mfd->label = _("Enable MFD flow (default is SFD (D8))");
+ flag_mfd->description =
+ _("SFD: single flow direction, MFD: multiple flow direction");
+
flag_flow = G_define_flag();
flag_flow->key = '4';
flag_flow->description =
@@ -192,9 +187,18 @@
flag_seg = G_define_flag();
flag_seg->key = 'm';
+ flag_seg->label =
+ _("Enable disk swap memory option: Operation is slow");
flag_seg->description =
- _("Enable disk swap memory option: Operation is slow");
+ _("Only needed if memory requirements exceed available RAM; see manual on how to calculate memory requirements");
+ flag_abs = G_define_flag();
+ flag_abs->key = 'a';
+ flag_abs->label =
+ _("Use positive flow accumulation even for likely underestimates");
+ flag_abs->description =
+ _("See manual for a detailed description of flow accumulation output");
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -214,18 +218,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();
@@ -248,9 +254,16 @@
strcat(command, "\"");
#endif
+ if (!flag_mfd->answer) {
+ strcat(command, " -s");
+ }
+
if (flag_flow->answer)
strcat(command, " -4");
+ if (flag_abs->answer)
+ strcat(command, " -a");
+
if (opt1->answer) {
strcat(command, " el=");
strcat(command, "\"");
@@ -352,13 +365,16 @@
strcat(command, "\"");
}
- if (flag_seg->answer && opt16->answer) {
- strcat(command, " mb=");
- strcat(command, "\"");
+ if (flag_mfd->answer && opt16->answer) {
+ strcat(command, " conv=");
strcat(command, opt16->answer);
- strcat(command, "\"");
}
+ if (flag_seg->answer && opt17->answer) {
+ strcat(command, " mb=");
+ strcat(command, opt17->answer);
+ }
+
#ifdef __MINGW32__
strcat(command, "\"");
#endif
@@ -368,45 +384,48 @@
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 */
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_mfd->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_mfd->answer);
if (opt10->answer)
write_hist(opt10->answer,
- "Watershed basins", opt1->answer, flag_seg->answer);
+ "Watershed basins", opt1->answer, flag_seg->answer,
+ flag_mfd->answer);
if (opt11->answer)
write_hist(opt11->answer,
- "Watershed stream segments", opt1->answer, flag_seg->answer);
+ "Watershed stream segments", opt1->answer,
+ flag_seg->answer, flag_mfd->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_mfd->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_mfd->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_mfd->answer);
if (opt15->answer)
write_hist(opt15->answer,
"Watershed slope steepness (S) factor",
- opt1->answer, flag_seg->answer);
+ opt1->answer, flag_seg->answer, flag_mfd->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 mfd)
{
struct History history;
@@ -416,8 +435,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", mfd ? "MFD" : "SFD (D8)");
+ 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/branches/releasebranch_6_4/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/Gwater.h 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/Gwater.h 2010-05-05 13:30:34 UTC (rev 42121)
@@ -34,65 +34,47 @@
#define POINT struct points
POINT {
- SHORT r, c, downr, downc;
+ SHORT r, c; /* , downr, downc */
int nxt;
};
-#ifdef MAIN
-#define GLOBAL
-#define DRAINVAR = {{ 7,6,5 },{ 8,0,4 },{ 1,2,3 }}
-#define UPDRAINVAR = {{ 3,2,1 },{ 4,0,8 },{ 5,6,7 }}
-#define NEXTDRVAR = { 1,-1,0,0,-1,1,1,-1 }
-#define NEXTDCVAR = { 0,0,-1,1,1,-1,1,-1 }
-#else
-#define GLOBAL extern
-#define DRAINVAR
-#define UPDRAINVAR
-#define NEXTDRVAR
-#define NEXTDCVAR
-#endif
+extern struct Cell_head window;
-GLOBAL struct Cell_head window;
+extern int mfd, c_fac, abs_acc, ele_scale;
+extern int *heap_index, heap_size;
+extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
+extern SHORT nrows, ncols;
+extern double half_res, diag, max_length, dep_slope;
+extern int bas_thres, tot_parts;
+extern FLAG *worked, *in_list, *s_b, *swale;
+extern RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
+extern RAMSEG r_h_seg, dep_seg;
+extern RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
+extern int *astar_pts;
+extern CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
+extern DCELL *wat;
+extern int ril_fd;
+extern double *s_l, *s_g, *l_s;
+extern CELL one, zero;
+extern double ril_value, d_one, d_zero;
+extern SHORT sides;
+extern SHORT drain[3][3];
+extern SHORT updrain[3][3];
+extern SHORT nextdr[8];
+extern SHORT nextdc[8];
+extern char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
+extern char run_name[GNAME_MAX], ob_name[GNAME_MAX];
+extern char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
+extern const char *this_mapset;
+extern char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX], thr_name[8];
+extern char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX], sg_name[GNAME_MAX];
+extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX], dis_name[GNAME_MAX];
+extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
+extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
+extern char bas_flag, seg_flag, haf_flag, er_flag;
+extern char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
+extern FILE *fp;
-GLOBAL int *heap_index, heap_size;
-GLOBAL int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
-GLOBAL SHORT nrows, ncols;
-GLOBAL double half_res, diag, max_length, dep_slope;
-GLOBAL int bas_thres, tot_parts;
-GLOBAL FLAG *worked, *in_list, *s_b, *swale;
-GLOBAL RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
-GLOBAL RAMSEG r_h_seg, dep_seg;
-GLOBAL RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
-GLOBAL POINT *astar_pts;
-GLOBAL CELL *dis, *alt, *wat, *asp, *bas, *haf, *r_h, *dep;
-GLOBAL CELL *ril_buf;
-GLOBAL int ril_fd;
-GLOBAL double *s_l, *s_g, *l_s;
-GLOBAL CELL one, zero;
-GLOBAL double ril_value, dzero;
-GLOBAL SHORT sides;
-GLOBAL SHORT drain[3][3] DRAINVAR;
-GLOBAL SHORT updrain[3][3] UPDRAINVAR;
-GLOBAL SHORT nextdr[8] NEXTDRVAR;
-GLOBAL SHORT nextdc[8] NEXTDCVAR;
-GLOBAL char ele_name[GNAME_MAX], *ele_mapset, pit_name[GNAME_MAX],
- *pit_mapset;
-GLOBAL char run_name[GNAME_MAX], *run_mapset, ob_name[GNAME_MAX], *ob_mapset;
-GLOBAL char ril_name[GNAME_MAX], *ril_mapset, dep_name[GNAME_MAX],
- *dep_mapset;
-GLOBAL char *this_mapset;
-GLOBAL char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX],
- thr_name[8];
-GLOBAL char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX],
- sg_name[GNAME_MAX];
-GLOBAL char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX],
- dis_name[GNAME_MAX];
-GLOBAL char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
-GLOBAL char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
-GLOBAL char bas_flag, seg_flag, haf_flag, er_flag;
-GLOBAL char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
-GLOBAL FILE *fp;
-
/* close_maps.c */
int close_maps(void);
@@ -104,13 +86,16 @@
/* 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);
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);
@@ -120,8 +105,6 @@
/* init_vars.c */
int init_vars(int, char *[]);
-int do_legal(char *);
-char *do_exist(char *);
/* no_stream.c */
int no_stream(int, int, CELL, double, CELL);
@@ -131,6 +114,7 @@
/* ramseg.c */
int size_array(int *, int, int);
+int seg_index_rc(int, int, int *, int *);
/* sg_factor.c */
int sg_factor(void);
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/ram/Makefile
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/Makefile 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/Makefile 2010-05-05 13:30:34 UTC (rev 42121)
@@ -10,4 +10,4 @@
default: etc
htmletc:
- @echo "No docs to generate."
+ @echo "No docs to generate."
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/ram/close_maps.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/close_maps.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/close_maps.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -9,30 +9,133 @@
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();
+ 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_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)];
+ if (abs_acc) {
+ G_warning("Writing out only positive flow accumulation values.");
+ G_warning("Cells with a likely underestimate for flow accumulation can no longer be identified.");
+ for (r = 0; r < nrows; r++) {
+ G_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
+ for (c = 0; c < ncols; c++) {
+ dvalue = wat[SEG_INDEX(wat_seg, r, c)];
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dvalue = ABS(dvalue);
+ dbuf[c] = dvalue;
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
+ }
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
- G_put_raster_row(fd, buf, CELL_TYPE);
}
+ 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++) {
+ 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 (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, 0, &clr_max, 0, 0,
+ 0, &colors);
+ }
+ 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) {
@@ -40,6 +143,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)];
}
@@ -54,6 +158,7 @@
}
G_free(asp);
+ /* visual output no longer needed */
if (dis_flag) {
fd = G_open_cell_new(dis_name);
if (fd < 0) {
@@ -91,16 +196,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."));
@@ -109,18 +214,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."));
@@ -131,16 +236,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/branches/releasebranch_6_4/raster/r.watershed/ram/close_maps2.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/close_maps2.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/close_maps2.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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) {
@@ -78,7 +78,7 @@
else
G_debug(1,
"Too many subbasins to reasonably check for color brightness");
- /* using the existing method ... */
+ /* using the existing stack of while/for/for/for/while loops ... */
}
/* stream segments map */
@@ -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/branches/releasebranch_6_4/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_astar.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_astar.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -3,110 +3,153 @@
#include <grass/gis.h>
#include <grass/glocale.h>
-int sift_up(int, CELL);
+double get_slope2(CELL, CELL, double);
int do_astar(void)
{
- POINT *point;
- int doer, count;
+ int count;
SHORT upr, upc, r, c, ct_dir;
- CELL alt_val, alt_up, asp_up, wat_val;
- CELL in_val, drain_val;
+ CELL alt_val, alt_nbr[8];
+ CELL is_in_list, is_worked;
int index_doer, index_up;
+ /* sides
+ * |7|1|4|
+ * |2| |3|
+ * |5|0|6|
+ */
+ int nbr_ew[8] = { 0, 1, 2, 3, 1, 0, 0, 1 };
+ int nbr_ns[8] = { 0, 1, 2, 3, 3, 2, 3, 2 };
+ double dx, dy, dist_to_nbr[8], ew_res, ns_res;
+ double slope[8];
+ int skip_diag;
-
G_message(_("SECTION 2: A * Search."));
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ upr = nextdr[ct_dir];
+ upc = nextdc[ct_dir];
+ /* account for rare cases when ns_res != ew_res */
+ dy = ABS(upr) * window.ns_res;
+ dx = ABS(upc) * 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);
+ }
+ ew_res = window.ew_res;
+ ns_res = window.ns_res;
+
count = 0;
first_astar = heap_index[1];
+ first_cum = do_points;
- /* A * Search: search uphill, get downhill path */
- while (first_astar != -1) {
+ /* A* Search: search uphill, get downhill paths */
+ while (heap_size > 0) {
G_percent(count++, do_points, 1);
/* start with point with lowest elevation, in case of equal elevation
* of following points, oldest point = point added earliest */
- /* old routine: astar_pts[first_astar] (doer = first_astar) */
- /* new routine: astar_pts[heap_index[1]] */
+ index_doer = astar_pts[1];
- doer = heap_index[1];
-
- point = &(astar_pts[doer]);
-
- /* drop astar_pts[doer] from heap */
- /* necessary to protect the current point (doer) from modification */
- /* equivalent to first_astar = point->next in old code */
drop_pt();
- /* can go, dragged on from old code */
- first_astar = heap_index[1];
+ /* add astar points to sorted list for flow accumulation */
+ astar_pts[first_cum] = index_doer;
+ first_cum--;
- /* downhill path for flow accumulation is set here */
- /* this path determines the order for flow accumulation calculation */
- point->nxt = first_cum;
- first_cum = doer;
+ seg_index_rc(alt_seg, index_doer, &r, &c);
- r = point->r;
- c = point->c;
+ G_debug(3, "A* Search: row %d, column %d, ", r, c);
- G_debug(3, "R:%2d C:%2d, ", r, c);
-
- index_doer = SEG_INDEX(alt_seg, r, c);
alt_val = alt[index_doer];
- wat_val = wat[index_doer];
FLAG_SET(worked, r, c);
- /* check all neighbours */
+ /* check neighbours */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (upr, upc) for this neighbour */
upr = r + nextdr[ct_dir];
upc = c + nextdc[ct_dir];
- index_up = SEG_INDEX(alt_seg, upr, upc);
+ slope[ct_dir] = alt_nbr[ct_dir] = 0;
/* check that r, c are within region */
if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
+ index_up = SEG_INDEX(alt_seg, upr, upc);
/* check if neighbour is in the list */
/* if not, add as new point */
- in_val = FLAG_GET(in_list, upr, upc);
- if (in_val == 0) {
- alt_up = alt[index_up];
- /* flow direction is set here */
- add_pt(upr, upc, r, c, alt_up, alt_val);
- drain_val = drain[upr - r + 1][upc - c + 1];
- asp[index_up] = drain_val;
+ is_in_list = FLAG_GET(in_list, upr, upc);
+ is_worked = FLAG_GET(worked, upr, upc);
+ skip_diag = 0;
+ /* avoid diagonal flow direction bias */
+ if (!is_worked) {
+ alt_nbr[ct_dir] = alt[index_up];
+ slope[ct_dir] =
+ get_slope2(alt_val, alt_nbr[ct_dir],
+ dist_to_nbr[ct_dir]);
+ }
+ if (!is_in_list) {
+ if (ct_dir > 3 && slope[ct_dir] > 0) {
+ if (slope[nbr_ew[ct_dir]] > 0) {
+ /* slope to ew nbr > slope to center */
+ if (slope[ct_dir] <
+ get_slope2(alt_nbr[nbr_ew[ct_dir]],
+ alt_nbr[ct_dir], ew_res))
+ skip_diag = 1;
+ }
+ if (!skip_diag && slope[nbr_ns[ct_dir]] > 0) {
+ /* slope to ns nbr > slope to center */
+ if (slope[ct_dir] <
+ get_slope2(alt_nbr[nbr_ns[ct_dir]],
+ alt_nbr[ct_dir], ns_res))
+ skip_diag = 1;
+ }
+ }
+ }
+ /* add neighbour as new point if not in the list */
+ if (is_in_list == 0 && skip_diag == 0) {
+ add_pt(upr, upc, alt_nbr[ct_dir], alt_val);
+ /* set flow direction */
+ asp[index_up] = drain[upr - r + 1][upc - c + 1];
}
- else {
- /* check if neighbour has not been worked on,
- * 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) {
- asp[index_up] = drain[upr - r + 1][upc - c + 1];
+ else if (is_in_list && is_worked == 0) {
+ /* neighbour is edge in list, not yet worked */
+ if (asp[index_up] < 0) {
+ asp[index_up] = drain[upr - r + 1][upc - c + 1];
- if (wat_val > 0)
- wat[index_doer] = -wat_val;
-
- replace(upr, upc, r, c); /* alt_up used to be */
- }
+ if (wat[index_doer] > 0)
+ wat[index_doer] = -wat[index_doer];
}
}
- }
- }
+ } /* end if in region */
+ } /* end sides */
}
- G_percent(count, do_points, 3); /* finish it */
- flag_destroy(worked);
+ G_percent(count, do_points, 1); /* finish it */
+ if (mfd == 0)
+ flag_destroy(worked);
+
flag_destroy(in_list);
G_free(heap_index);
return 0;
}
+/* compare two heap points */
+/* return 1 if a < b else 0 */
+int cmp_pnt(CELL elea, CELL eleb, int addeda, int addedb)
+{
+ if (elea < eleb)
+ return 1;
+ else if (elea == eleb) {
+ if (addeda < addedb)
+ return 1;
+ }
+ return 0;
+}
+
/* 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);
/* add point to next free position */
@@ -116,15 +159,10 @@
if (heap_size > do_points)
G_fatal_error(_("heapsize too large"));
- heap_index[heap_size] = nxt_avail_pt;
+ heap_index[heap_size] = nxt_avail_pt++;
- 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[heap_size] = SEG_INDEX(alt_seg, r, c);
- nxt_avail_pt++;
-
/* sift up: move new point towards top of heap */
sift_up(heap_size, ele);
@@ -135,9 +173,9 @@
/* new drop point routine for min heap */
int drop_pt(void)
{
- int child, childr, parent;
+ register int child, childr, parent;
CELL ele, eler;
- int i;
+ register int i;
if (heap_size == 1) {
heap_index[1] = -1;
@@ -153,38 +191,34 @@
while ((child = GET_CHILD(parent)) <= heap_size) {
/* select child with lower ele, if equal, older child
- * older child is older startpoint for flow path, important
- * chances are good that the left child will be the older child,
- * but just to make sure we test */
- ele =
- alt[SEG_INDEX
- (alt_seg, astar_pts[heap_index[child]].r,
- astar_pts[heap_index[child]].c)];
+ * older child is older startpoint for flow path, important */
+ ele = alt[astar_pts[child]];
if (child < heap_size) {
childr = child + 1;
- i = 1;
- while (childr <= heap_size && i < 3) {
- eler =
- alt[SEG_INDEX
- (alt_seg, astar_pts[heap_index[childr]].r,
- astar_pts[heap_index[childr]].c)];
- if (eler < ele) {
+ i = child + 3;
+ while (childr <= heap_size && childr < i) {
+ eler = alt[astar_pts[childr]];
+ /* get smallest child */
+ if (cmp_pnt(eler, ele, heap_index[childr], heap_index[child])) {
child = childr;
ele = eler;
- /* make sure we get the oldest child */
}
- else if (ele == eler &&
- heap_index[child] > heap_index[childr]) {
- child = childr;
- }
childr++;
- i++;
}
+ /* break if childr > last entry? that saves sifting up again
+ * OTOH, this is another comparison
+ * we have a max heap height of 20: log(INT_MAX)/log(n children per node)
+ * that would give us in the worst case 20*2 additional comparisons with 3 children
+ * the last entry will never go far up again, less than half the way
+ * so the additional comparisons for going all the way down
+ * and then a bit up again are likely less than 20*2 */
+ /* find the error in this reasoning */
}
/* move hole down */
heap_index[parent] = heap_index[child];
+ astar_pts[parent] = astar_pts[child];
parent = child;
}
@@ -192,14 +226,11 @@
/* hole is in lowest layer, move to heap end */
if (parent < heap_size) {
heap_index[parent] = heap_index[heap_size];
+ astar_pts[parent] = astar_pts[heap_size];
- ele =
- alt[SEG_INDEX
- (alt_seg, astar_pts[heap_index[parent]].r,
- astar_pts[heap_index[parent]].c)];
+ ele = alt[astar_pts[parent]];
/* sift up last swapped point, only necessary if hole moved to heap end */
sift_up(parent, ele);
-
}
/* the actual drop */
@@ -211,43 +242,33 @@
/* standard sift-up routine for d-ary min heap */
int sift_up(int start, CELL ele)
{
- int parent, parentp, child, childp;
+ register int parent, child, child_idx, child_added;
CELL elep;
child = start;
- childp = heap_index[child];
+ child_added = heap_index[child];
+ child_idx = astar_pts[child];
while (child > 1) {
parent = GET_PARENT(child);
- parentp = heap_index[parent];
- elep =
- alt[SEG_INDEX
- (alt_seg, astar_pts[parentp].r, astar_pts[parentp].c)];
- /* parent ele higher */
- if (elep > ele) {
-
+ elep = alt[astar_pts[parent]];
+ /* child smaller */
+ if (cmp_pnt(ele, elep, child_added, heap_index[parent])) {
/* push parent point down */
- heap_index[child] = parentp;
+ heap_index[child] = heap_index[parent];
+ astar_pts[child] = astar_pts[parent];
child = parent;
-
}
- /* same ele, but parent is younger */
- else if (elep == ele && parentp > childp) {
-
- /* push parent point down */
- heap_index[child] = parentp;
- child = parent;
-
- }
else
/* no more sifting up, found new slot for child */
break;
}
- /* set heap_index for child */
+ /* put point in new slot */
if (child < start) {
- heap_index[child] = childp;
+ heap_index[child] = child_added;
+ astar_pts[child] = child_idx;
}
return 0;
@@ -272,6 +293,13 @@
return (slope);
}
+double get_slope2(CELL ele, CELL up_ele, double dist)
+{
+ if (ele >= up_ele)
+ return 0.0;
+ else
+ return (double)(up_ele - ele) / dist;
+}
/* new replace */
int replace( /* ele was in there */
@@ -279,6 +307,7 @@
/* CELL ele; */
{
int now, heap_run;
+ int r2, c2;
/* find the current neighbour point and
* set flow direction to focus point */
@@ -287,9 +316,11 @@
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;
+ /* if (astar_pts[now].r == upr && astar_pts[now].c == upc) { */
+ seg_index_rc(alt_seg, astar_pts[now], &r2, &c2);
+ if (r2 == upr && c2 == upc) {
+ /* astar_pts[now].downr = r;
+ astar_pts[now].downc = c; */
return 0;
}
heap_run++;
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_astar.h
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_astar.h 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_astar.h 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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/branches/releasebranch_6_4/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_cum.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/do_cum.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -6,28 +6,36 @@
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 };
+ int this_index, down_index;
- G_message(_("SECTION 3: Accumulating Surface Flow."));
+ G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
count = 0;
if (bas_thres <= 0)
threshold = 60;
else
threshold = bas_thres;
- while (first_cum != -1) {
- G_percent(count++, do_points, 3);
- 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)];
- if (ABS(value) >= threshold)
+ for (killer = 1; killer <= do_points; killer++) {
+ G_percent(killer, do_points, 1);
+ this_index = astar_pts[killer];
+ aspect = asp[this_index];
+ seg_index_rc(alt_seg, this_index, &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) { */
+ down_index = SEG_INDEX(wat_seg, dr, dc);
+ value = wat[this_index];
+ if ((int)(ABS(value) + 0.5) >= threshold)
FLAG_SET(swale, r, c);
- valued = wat[SEG_INDEX(wat_seg, dr, dc)];
+ valued = wat[down_index];
if (value > 0) {
if (valued > 0)
valued += value;
@@ -40,9 +48,17 @@
else
valued = value - valued;
}
- wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
+ wat[down_index] = valued;
+ valued = ABS(valued) + 0.5;
is_swale = FLAG_GET(swale, r, c);
- if (is_swale || ABS(valued) >= threshold) {
+ /* update asp for depression */
+ if (is_swale && pit_flag) {
+ if (aspect > 0 && asp[down_index] == 0) {
+ aspect = -aspect;
+ asp[this_index] = aspect;
+ }
+ }
+ if (is_swale || ((int)valued) >= threshold) {
FLAG_SET(swale, dr, dc);
}
else {
@@ -51,8 +67,326 @@
}
}
}
- G_percent(count, do_points, 3); /* 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;
+ int killer, threshold, count;
+
+ /* MFD */
+ int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
+ double *dist_to_nbr, *weight, sum_weight, max_weight;
+ int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
+ double dx, dy;
+ CELL ele, ele_nbr, aspect, is_worked;
+ double prop, max_acc;
+ int workedon, edge, flat;
+ 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 };
+ int this_index, down_index, nbr_index;
+
+ 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));
+ weight = (double *)G_malloc(sides * sizeof(double));
+
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ r_nbr = nextdr[ct_dir];
+ c_nbr = nextdc[ct_dir];
+ /* account for rare cases when ns_res != ew_res */
+ dy = ABS(r_nbr) * window.ns_res;
+ dx = ABS(c_nbr) * window.ew_res;
+ if (ct_dir < 4)
+ dist_to_nbr[ct_dir] = dx + dy;
+ else
+ dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+ }
+
+ flag_clear_all(worked);
+ workedon = 0;
+
+ count = 0;
+ if (bas_thres <= 0)
+ threshold = 60;
+ else
+ threshold = bas_thres;
+
+ for (killer = 1; killer <= do_points; killer++) {
+ G_percent(killer, do_points, 1);
+ this_index = astar_pts[killer];
+ seg_index_rc(alt_seg, this_index, &r, &c);
+ aspect = asp[this_index];
+ 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[this_index];
+ down_index = SEG_INDEX(wat_seg, dr, dc);
+ valued = wat[down_index];
+
+ 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[this_index];
+ is_null = 0;
+ edge = 0;
+ flat = 1;
+ /* 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 */
+ 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) {
+
+ nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
+
+ /* check for swale or stream cells */
+ is_swale = FLAG_GET(swale, r_nbr, c_nbr);
+ if (is_swale)
+ swale_cells++;
+ valued = wat[nbr_index];
+ ele_nbr = alt[nbr_index];
+ if ((ABS(valued) + 0.5) >= threshold &&
+ ct_dir != np_side && ele_nbr > ele)
+ stream_cells++;
+
+ is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ if (is_worked == 0) {
+ is_null = G_is_c_null_value(&ele_nbr);
+ edge = is_null;
+ if (ele_nbr != ele)
+ flat = 0;
+ 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) {
+ 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;
+ }
+ else
+ edge = 1;
+ if (edge)
+ break;
+ }
+ /* do not distribute flow along edges, this causes artifacts */
+ if (edge) {
+ is_swale = FLAG_GET(swale, r, c);
+ if (is_swale && aspect > 0) {
+ aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+ asp[this_index] = aspect;
+ }
+ continue;
+ }
+
+ /* 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 > 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;
+ }
+
+ /* set flow accumulation for neighbours */
+ max_acc = -1;
+
+ 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) {
+ is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ if (is_worked == 0) {
+
+ nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
+
+ weight[ct_dir] = weight[ct_dir] / sum_weight;
+ /* check everything sums up to 1.0 */
+ prop += weight[ct_dir];
+
+ valued = wat[nbr_index];
+ 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[nbr_index] = valued;
+
+ /* get main drainage direction */
+ if (ABS(valued) >= max_acc) {
+ max_acc = ABS(valued);
+ r_max = r_nbr;
+ c_max = c_nbr;
+ }
+ }
+ else if (ct_dir == np_side) {
+ /* check for consistency with A * path */
+ workedon++;
+ }
+ }
+ }
+ if (ABS(prop - 1.0) > 5E-6f) {
+ G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
+ prop);
+ }
+ }
+
+ if (mfd_cells < 2) {
+ valued = wat[down_index];
+ if (value > 0) {
+ if (valued > 0)
+ valued += value;
+ else
+ valued -= value;
+ }
+ else {
+ if (valued < 0)
+ valued += value;
+ else
+ valued = value - valued;
+ }
+ wat[down_index] = valued;
+ }
+
+ /* update asp */
+ if (dr != r_max || dc != c_max) {
+ aspect = drain[r - r_max + 1][c - c_max + 1];
+ if (asp[this_index] < 0)
+ aspect = -aspect;
+ asp[this_index] = aspect;
+ }
+ is_swale = FLAG_GET(swale, r, c);
+ /* 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[this_index] = aspect;
+ }
+ }
+ /* start new stream */
+ value = ABS(value) + 0.5;
+ if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
+ swale_cells < 1 && !flat) {
+ FLAG_SET(swale, r, c);
+ is_swale = 1;
+ }
+ /* continue stream */
+ if (is_swale) {
+ FLAG_SET(swale, r_max, c_max);
+ }
+ else {
+ if (er_flag && !is_swale)
+ slope_length(r, c, r_max, c_max);
+ }
+ FLAG_SET(worked, r, c);
+ }
+ }
+ if (workedon)
+ 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;
+}
+
+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/branches/releasebranch_6_4/raster/r.watershed/ram/find_pour.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/find_pour.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/find_pour.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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/branches/releasebranch_6_4/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/init_vars.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/init_vars.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -1,26 +1,37 @@
#include <stdlib.h>
+#include <string.h>
#include "Gwater.h"
#include <grass/gis.h>
#include <grass/glocale.h>
+int ele_round(double);
int init_vars(int argc, char *argv[])
{
SHORT r, c;
- CELL *buf, alt_value, wat_value, asp_value;
- int fd, index;
+ CELL *buf, alt_value, asp_value, block_value;
+ DCELL dvalue, wat_value;
+ void *elebuf, *ptr;
+ int fd, index, ele_map_type;
+ size_t ele_size;
char MASK_flag;
G_gisinit(argv[0]);
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;
+ abs_acc = 0;
+ ele_scale = 1;
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
ele_flag++;
@@ -63,9 +74,17 @@
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 if (strcmp(argv[r], "-a") == 0)
+ abs_acc = 1;
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) &&
@@ -83,29 +102,8 @@
tot_parts++;
G_message(_("SECTION 1a (of %1d): Initiating Memory."), tot_parts);
this_mapset = G_mapset();
- if (asp_flag)
- do_legal(asp_name);
- if (bas_flag)
- do_legal(bas_name);
- if (seg_flag)
- do_legal(seg_name);
- if (haf_flag)
- do_legal(haf_name);
- if (sl_flag)
- do_legal(sl_name);
- if (sg_flag)
- do_legal(sg_name);
- if (ls_flag)
- do_legal(ls_name);
if (sl_flag || sg_flag || ls_flag)
er_flag = 1;
- ele_mapset = do_exist(ele_name);
- if (pit_flag)
- pit_mapset = do_exist(pit_name);
- if (ob_flag)
- ob_mapset = do_exist(ob_name);
- if (ril_flag)
- ril_mapset = do_exist(ril_name);
/* for sd factor
if (dep_flag) {
if (sscanf (dep_name, "%lf", &dep_slope) != 1) {
@@ -118,7 +116,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;
@@ -128,126 +126,162 @@
window.ns_res * window.ns_res);
if (sides == 4)
diag *= 0.5;
- 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));
+ wat =
+ (DCELL *) G_malloc(sizeof(DCELL) *
+ size_array(&wat_seg, nrows, ncols));
+ asp =
+ (CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
- fd = G_open_cell_old(ele_name, ele_mapset);
- if (fd < 0) {
- G_fatal_error(_("unable to open elevation map layer"));
+ if (er_flag) {
+ r_h =
+ (CELL *) G_malloc(sizeof(CELL) * size_array(&r_h_seg, nrows, ncols));
}
swale = flag_create(nrows, ncols);
in_list = flag_create(nrows, ncols);
worked = flag_create(nrows, ncols);
+ /* open elevation input */
+ fd = G_open_cell_old(ele_name, "");
+ if (fd < 0) {
+ G_fatal_error(_("unable to open elevation map layer"));
+ }
+
+ ele_map_type = G_get_raster_map_type(fd);
+ ele_size = G_raster_size(ele_map_type);
+ elebuf = G_allocate_raster_buf(ele_map_type);
+
+ if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
+ ele_scale = 1000; /* should be enough to do the trick */
+
+ /* read elevation input and mark NULL/masked cells */
+ /* intialize accumulation and drainage direction */
+ MASK_flag = 0;
+ do_points = nrows * ncols;
for (r = 0; r < nrows; r++) {
- G_get_c_raster_row(fd, buf, r);
+ G_get_raster_row(fd, elebuf, r, ele_map_type);
+ ptr = elebuf;
for (c = 0; c < ncols; c++) {
index = SEG_INDEX(alt_seg, r, c);
- alt[index] = 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 */
+ if (G_is_null_value(ptr, ele_map_type)) {
+ FLAG_SET(worked, r, c);
+ FLAG_SET(in_list, r, c);
+ G_set_c_null_value(&alt_value, 1);
+ G_set_d_null_value(&wat_value, 1);
+ do_points--;
+ }
+ else {
+ if (ele_map_type == CELL_TYPE) {
+ alt_value = *((CELL *)ptr);
+ }
+ else if (ele_map_type == FCELL_TYPE) {
+ dvalue = *((FCELL *)ptr);
+ dvalue *= ele_scale;
+ alt_value = ele_round(dvalue);
+ }
+ else if (ele_map_type == DCELL_TYPE) {
+ dvalue = *((DCELL *)ptr);
+ dvalue *= ele_scale;
+ alt_value = ele_round(dvalue);
+ }
+ wat_value = 1.0;
+ }
+ alt[index] = alt_value;
+ wat[index] = wat_value;
+ asp[index] = 0;
+ if (er_flag) {
+ r_h[index] = alt_value;
+ }
+ ptr = G_incr_void_ptr(ptr, ele_size);
}
}
G_close_cell(fd);
- wat =
- (CELL *) G_malloc(sizeof(CELL) * size_array(&wat_seg, nrows, ncols));
+ G_free(elebuf);
+ if (do_points < nrows * ncols)
+ MASK_flag = 1;
+ /* read flow accumulation from input map flow: amount of overland flow per cell */
+ buf = G_allocate_cell_buf();
if (run_flag) {
- run_mapset = do_exist(run_name);
- fd = G_open_cell_old(run_name, run_mapset);
+ fd = G_open_cell_old(run_name, "");
if (fd < 0) {
G_fatal_error(_("unable to open runoff map layer"));
}
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;
- }
- }
- asp = (CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
+ /* depression: drainage direction will be set to zero later */
if (pit_flag) {
- pit_mapset = do_exist(pit_name);
- fd = G_open_cell_old(pit_name, pit_mapset);
+ fd = G_open_cell_old(pit_name, "");
if (fd < 0) {
G_fatal_error(_("unable to open depression map layer"));
}
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_value)
+ asp[SEG_INDEX(asp_seg, r, c)] = 1;
}
}
G_close_cell(fd);
}
- else {
- for (r = 0; r < nrows; r++) {
- for (c = 0; c < ncols; c++)
- asp[SEG_INDEX(wat_seg, r, c)] = 0;
- }
- }
+ /* this is also creating streams... */
if (ob_flag) {
- fd = G_open_cell_old(ob_name, ob_mapset);
+ fd = G_open_cell_old(ob_name, "");
if (fd < 0) {
G_fatal_error(_("unable to open blocking map layer"));
}
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) && block_value)
FLAG_SET(swale, r, c);
}
}
G_close_cell(fd);
}
+ G_free(buf);
+
if (ril_flag) {
- ril_fd = G_open_cell_old(ril_name, ril_mapset);
+ ril_fd = G_open_cell_old(ril_name, "");
if (ril_fd < 0) {
G_fatal_error(_("unable to open rill map layer"));
}
}
- 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 =
@@ -260,7 +294,9 @@
sizeof(double));
}
- /* heap_index will track astar_pts in the binary min-heap */
+ astar_pts = (int *) G_malloc((do_points + 1) * sizeof(int));
+
+ /* 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));
@@ -278,14 +314,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;
@@ -293,15 +335,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)
@@ -309,7 +349,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)
@@ -317,7 +357,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)
@@ -325,7 +365,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)
@@ -333,7 +373,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)
@@ -341,7 +381,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)
@@ -349,7 +389,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)
@@ -357,7 +397,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)
@@ -371,7 +411,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) {
@@ -379,7 +420,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;
@@ -387,35 +433,29 @@
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;
}
-int do_legal(char *file_name)
+int ele_round(double x)
{
- if (G_legal_filename(file_name) == -1)
- G_fatal_error(_("<%s> is an illegal file name"), file_name);
+ int n;
- return 0;
-}
+ if (x >= 0.0)
+ n = x + .5;
+ else {
+ n = -x + .5;
+ n = -n;
+ }
-char *do_exist(char *file_name)
-{
- char *file_mapset = G_find_cell2(file_name, "");
-
- if (file_mapset == NULL)
- G_fatal_error(_("Raster map <%s> not found"), file_name);
-
- return (file_mapset);
+ return n;
}
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/ram/main.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/main.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/main.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -5,30 +5,73 @@
* 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
+ * PURPOSE: Hydrological analysis
+ * COPYRIGHT: (C) 1999-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
* for details.
*
*****************************************************************************/
-#define MAIN
+
#include <stdlib.h>
#include <unistd.h>
#include "Gwater.h"
#include <grass/gis.h>
#include <grass/glocale.h>
-#undef MAIN
+struct Cell_head window;
+int mfd, c_fac, abs_acc, ele_scale;
+int *heap_index, heap_size;
+int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
+SHORT nrows, ncols;
+double half_res, diag, max_length, dep_slope;
+int bas_thres, tot_parts;
+FLAG *worked, *in_list, *s_b, *swale;
+RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
+RAMSEG r_h_seg, dep_seg;
+RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
+int *astar_pts;
+CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
+DCELL *wat;
+int ril_fd;
+double *s_l, *s_g, *l_s;
+CELL one, zero;
+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 };
+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 ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
+char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
+char bas_flag, seg_flag, haf_flag, er_flag;
+char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
+FILE *fp;
+
int main(int argc, char *argv[])
{
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/branches/releasebranch_6_4/raster/r.watershed/ram/no_stream.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/no_stream.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/no_stream.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -5,23 +5,25 @@
{
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 dvalue, max_drain; /* flow acc is now DCELL */
SHORT updir, riteflag, leftflag, thisdir;
while (1) {
+ bas[SEG_INDEX(bas_seg, row, col)] = basin_num;
max_drain = -1;
for (r = row - 1, rr = 0; r <= row + 1; r++, rr++) {
for (c = col - 1, cc = 0; c <= col + 1; c++, cc++) {
if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
aspect = asp[SEG_INDEX(asp_seg, r, c)];
if (aspect == drain[rr][cc]) {
- value = wat[SEG_INDEX(wat_seg, r, c)];
- if (value < 0)
- value = -value;
- if (value > max_drain) {
+ dvalue = wat[SEG_INDEX(wat_seg, r, c)];
+ if (dvalue < 0)
+ dvalue = -dvalue;
+ if ((dvalue - max_drain) > 5E-8f) { /* floating point comparison problem workaround */
uprow = r;
upcol = c;
- max_drain = value;
+ max_drain = dvalue;
}
}
}
@@ -64,23 +66,25 @@
aspect = asp[SEG_INDEX(asp_seg, r, c)];
if (aspect == drain[rr][cc]) {
thisdir = updrain[rr][cc];
- if (haf_basin_side(updir,
- (SHORT) downdir,
- thisdir) == RITE) {
+ switch (haf_basin_side(updir,
+ (SHORT) downdir,
+ thisdir)) {
+ case RITE:
overland_cells(r, c, basin_num, basin_num,
&new_ele);
riteflag++;
- }
- else {
+ break;
+ case LEFT:
overland_cells(r, c, basin_num, basin_num - 1,
&new_ele);
leftflag++;
+ break;
}
}
}
}
}
- if (leftflag >= riteflag)
+ if (leftflag > riteflag)
haf[SEG_INDEX(haf_seg, row, col)] = basin_num - 1;
else
haf[SEG_INDEX(haf_seg, row, col)] = basin_num;
@@ -95,6 +99,7 @@
slope = MIN_SLOPE;
fprintf(fp, " %f %f\n", slope, stream_length);
}
+ haf[SEG_INDEX(haf_seg, row, col)] = basin_num;
return 0;
}
}
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/ram/ramseg.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/ramseg.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/ramseg.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -13,3 +13,16 @@
size -= (*ram_seg << RAMSEGBITS) - ncols;
return (size);
}
+
+/* get r, c from seg_index */
+int seg_index_rc(int ramseg, int seg_index, int *r, int *c)
+{
+ int seg_no, seg_remainder;
+
+ seg_no = seg_index >> DOUBLEBITS;
+ seg_remainder = seg_index - (seg_no << DOUBLEBITS);
+ *r = ((seg_no / ramseg) << RAMSEGBITS) + (seg_remainder >> RAMSEGBITS);
+ *c = ((seg_no - ((*r) >> RAMSEGBITS) * ramseg) << RAMSEGBITS) +
+ seg_remainder - (((*r) & SEGLENLESS) << RAMSEGBITS);
+ return seg_no;
+}
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/ram/sg_factor.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/sg_factor.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/sg_factor.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -2,6 +2,7 @@
#include <grass/gis.h>
#include <grass/glocale.h>
+CELL *ril_buf;
int sg_factor(void)
{
@@ -9,8 +10,11 @@
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."));
+ if (ril_flag)
+ ril_buf = G_allocate_cell_buf();
+
for (r = 0; r < nrows; r++) {
G_percent(r, nrows, 3);
if (ril_flag) {
@@ -20,7 +24,7 @@
low_elev = alt[SEG_INDEX(alt_seg, r, c)];
hih_elev = r_h[SEG_INDEX(r_h_seg, r, c)];
length = s_l[SEG_INDEX(s_l_seg, r, c)];
- height = hih_elev - low_elev;
+ height = 1.0 * (hih_elev - low_elev) / ele_scale;
if (length > max_length) {
height *= max_length / length;
length = max_length;
@@ -31,14 +35,14 @@
else
S = 16.8 * sin_theta - .50;
if (sg_flag)
- s_g[SEG_INDEX(s_g_seg, r, c)] = S;
+ s_g[SEG_INDEX(s_g_seg, r, c)] = S * 100; /* factor 100 for backwards compatibility */
if (ls_flag) {
length *= METER_TO_FOOT;
len_slp_equ(length, sin_theta, S, r, c);
}
}
}
- G_percent(r, nrows, 3); /* finish it */
+ G_percent(nrows, nrows, 1); /* finish it */
if (ril_flag) {
G_free(ril_buf);
@@ -66,7 +70,7 @@
/* rill_ratio equation from Steve Warren */
rill_ratio *= .5 + .005 * ril + .0001 * ril * ril;
s_l_exp = rill_ratio / (1 + rill_ratio);
- L = 100 * pow((slope_length / 72.6), s_l_exp);
+ L = 100 * pow((slope_length / 72.6), s_l_exp); /* factor 100 for backwards compatibility */
l_s[SEG_INDEX(l_s_seg, r, c)] = L * S;
return 0;
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/ram/slope_len.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/ram/slope_len.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/ram/slope_len.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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/branches/releasebranch_6_4/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/Gwater.h 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/Gwater.h 2010-05-05 13:30:34 UTC (rev 42121)
@@ -31,67 +31,50 @@
#define POINT struct points
POINT {
- SHORT r, c, downr, downc;
+ SHORT r, c; /* , downr, downc */
int nxt;
};
-#ifdef MAIN
-#define GLOBAL
-#define DRAINVAR = {{ 7,6,5 },{ 8,0,4 },{ 1,2,3 }}
-#define UPDRAINVAR = {{ 3,2,1 },{ 4,0,8 },{ 5,6,7 }}
-#define NEXTDRVAR = { 1,-1,0,0,-1,1,1,-1 }
-#define NEXTDCVAR = { 0,0,-1,1,1,-1,1,-1 }
-#else
-#define GLOBAL extern
-#define DRAINVAR
-#define UPDRAINVAR
-#define NEXTDRVAR
-#define NEXTDCVAR
-#endif
-
#define HEAP struct heap_item
HEAP {
int point;
CELL ele;
};
-GLOBAL struct Cell_head window;
+extern struct Cell_head window;
-GLOBAL SSEG heap_index;
-GLOBAL int heap_size;
-GLOBAL int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
-GLOBAL SHORT nrows, ncols;
-GLOBAL double half_res, diag, max_length, dep_slope;
-GLOBAL int bas_thres, tot_parts;
-GLOBAL SSEG astar_pts;
-GLOBAL BSEG worked, in_list, s_b, swale;
-GLOBAL CSEG dis, alt, wat, asp, bas, haf, r_h, dep;
-GLOBAL DSEG slp, s_l, s_g, l_s, ril;
-GLOBAL CELL one, zero;
-GLOBAL double ril_value, dzero;
-GLOBAL SHORT sides;
-GLOBAL SHORT drain[3][3] DRAINVAR;
-GLOBAL SHORT updrain[3][3] UPDRAINVAR;
-GLOBAL SHORT nextdr[8] NEXTDRVAR;
-GLOBAL SHORT nextdc[8] NEXTDCVAR;
-GLOBAL char ele_name[GNAME_MAX], *ele_mapset, pit_name[GNAME_MAX],
- *pit_mapset;
-GLOBAL char run_name[GNAME_MAX], *run_mapset, ob_name[GNAME_MAX], *ob_mapset;
-GLOBAL char ril_name[GNAME_MAX], *ril_mapset, dep_name[GNAME_MAX],
- *dep_mapset;
+extern int mfd, c_fac, abs_acc, ele_scale;
+extern SSEG heap_index;
+extern int heap_size;
+extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
+extern SHORT nrows, ncols;
+extern double half_res, diag, max_length, dep_slope;
+extern int bas_thres, tot_parts;
+extern SSEG astar_pts;
+extern BSEG worked, in_list, s_b, swale;
+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, d_zero, d_one;
+extern SHORT sides;
+extern SHORT drain[3][3];
+extern SHORT updrain[3][3];
+extern SHORT nextdr[8];
+extern SHORT nextdc[8];
+extern char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
+extern char run_name[GNAME_MAX], ob_name[GNAME_MAX];
+extern char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
-GLOBAL char *this_mapset;
-GLOBAL char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX],
- thr_name[8];
-GLOBAL char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX],
- sg_name[GNAME_MAX];
-GLOBAL char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX],
- dis_name[GNAME_MAX];
-GLOBAL char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
-GLOBAL char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
-GLOBAL char bas_flag, seg_flag, haf_flag, er_flag;
-GLOBAL char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
-GLOBAL FILE *fp;
+extern const char *this_mapset;
+extern char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX], thr_name[8];
+extern char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX], sg_name[GNAME_MAX];
+extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX], dis_name[GNAME_MAX];
+extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
+extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
+extern char bas_flag, seg_flag, haf_flag, er_flag;
+extern char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
+extern FILE *fp;
/* close_maps.c */
int close_maps(void);
@@ -104,13 +87,16 @@
/* 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);
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);
@@ -120,7 +106,6 @@
/* init_vars.c */
int init_vars(int, char *[]);
-int do_legal(char *);
char *do_exist(char *);
/* no_stream.c */
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/seg/close_maps.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/close_maps.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/close_maps.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -1,17 +1,137 @@
#include "Gwater.h"
#include <unistd.h>
+#include <grass/glocale.h>
int close_maps(void)
{
struct Colors colors;
int r, c;
- CELL is_swale, value;
- double dvalue;
+ CELL is_swale;
+ 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)
- cseg_write_cellfile(&wat, wat_name);
+ if (wat_flag) {
+ sum = sum_sqr = stddev = 0.0;
+ dbuf = G_allocate_d_raster_buf();
+ if (abs_acc) {
+ G_warning("Writing out only positive flow accumulation values.");
+ G_warning("Cells with a likely underestimate for flow accumulation can no longer be identified.");
+
+ fd = G_open_raster_new(wat_name, DCELL_TYPE);
+ if (fd < 0) {
+ G_warning(_("unable to open new accum map layer."));
+ }
+ for (r = 0; r < nrows; r++) {
+ G_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
+ for (c = 0; c < ncols; c++) {
+ dseg_get(&wat, &dvalue, r, c);
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dvalue = ABS(dvalue);
+ dbuf[c] = 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."));
+ }
+ else {
+ 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"));
+ }
+
+ 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 (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, 0, &clr_max, 0, 0, 0,
+ &colors);
+ }
+ G_write_colors(wat_name, this_mapset, &colors);
+ }
if (asp_flag) {
cseg_write_cellfile(&asp, asp_name);
G_init_colors(&colors);
@@ -19,26 +139,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);
@@ -46,7 +167,7 @@
/* error in gislib.a
G_free_colors(&colors);
*/
- cseg_close(&wat);
+ dseg_close(&wat);
if (ls_flag) {
dseg_write_cellfile(&l_s, ls_name);
dseg_close(&l_s);
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/seg/close_maps2.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/close_maps2.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/close_maps2.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -67,7 +67,7 @@
else
G_debug(1,
"Too many subbasins to reasonably check for color brightness");
- /* using the existing method ... */
+ /* using the existing stack of while/for/for/for/while loops ... */
}
/* stream segments map */
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_astar.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_astar.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_astar.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -5,21 +5,45 @@
#include "Gwater.h"
#include "do_astar.h"
-int sift_up(int, CELL);
+double get_slope2(CELL, CELL, double);
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_nbr[8], alt_up, asp_up;
+ DCELL wat_val;
CELL in_val, drain_val;
HEAP heap_pos;
+ /* sides
+ * |7|1|4|
+ * |2| |3|
+ * |5|0|6|
+ */
+ int nbr_ew[8] = { 0, 1, 2, 3, 1, 0, 0, 1 };
+ int nbr_ns[8] = { 0, 1, 2, 3, 3, 2, 3, 2 };
+ double dx, dy, dist_to_nbr[8], ew_res, ns_res;
+ double slope[8];
+ int skip_diag;
- /* double slope; */
-
G_message(_("SECTION 2: A * Search."));
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ upr = nextdr[ct_dir];
+ upc = nextdc[ct_dir];
+ /* account for rare cases when ns_res != ew_res */
+ dy = ABS(upr) * window.ns_res;
+ dx = ABS(upc) * 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);
+ }
+ ew_res = window.ew_res;
+ ns_res = window.ns_res;
+
count = 0;
seg_get(&heap_index, (char *)&heap_pos, 0, 1);
first_astar = heap_pos.point;
@@ -33,6 +57,8 @@
seg_get(&astar_pts, (char *)&point, 0, doer);
+ alt_val = heap_pos.ele;
+
/* drop astar_pts[doer] from heap */
drop_pt();
@@ -47,49 +73,71 @@
c = point.c;
bseg_put(&worked, &one, r, c);
- cseg_get(&alt, &alt_val, r, c);
/* check all neighbours, breadth first search */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (upr, upc) for this neighbour */
upr = r + nextdr[ct_dir];
upc = c + nextdc[ct_dir];
+ slope[ct_dir] = alt_nbr[ct_dir] = 0;
/* check that upr, upc are within region */
if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
- /* check if neighbour is in the list */
- /* if not, add as new point */
+ /* avoid diagonal flow direction bias */
bseg_get(&in_list, &in_val, upr, upc);
- if (in_val == 0) {
+ bseg_get(&worked, &work_val, upr, upc);
+ skip_diag = 0;
+ if (!work_val) {
cseg_get(&alt, &alt_up, upr, upc);
- add_pt(upr, upc, r, c, alt_up, alt_val);
+ alt_nbr[ct_dir] = alt_up;
+ slope[ct_dir] =
+ get_slope2(alt_val, alt_nbr[ct_dir],
+ dist_to_nbr[ct_dir]);
+ }
+ if (!in_val) {
+ if (ct_dir > 3 && slope[ct_dir] > 0) {
+ if (slope[nbr_ew[ct_dir]] > 0) {
+ /* slope to ew nbr > slope to center */
+ if (slope[ct_dir] <
+ get_slope2(alt_nbr[nbr_ew[ct_dir]],
+ alt_nbr[ct_dir], ew_res))
+ skip_diag = 1;
+ }
+ if (!skip_diag && slope[nbr_ns[ct_dir]] > 0) {
+ /* slope to ns nbr > slope to center */
+ if (slope[ct_dir] <
+ get_slope2(alt_nbr[nbr_ns[ct_dir]],
+ alt_nbr[ct_dir], ns_res))
+ skip_diag = 1;
+ }
+ }
+ }
+ /* put neighbour in search list if not yet in */
+ if (in_val == 0 && skip_diag == 0) {
+ add_pt(upr, upc, alt_nbr[ct_dir], 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 */
- bseg_get(&worked, &work_val, upr, upc);
- if (!work_val) {
- cseg_get(&asp, &asp_up, upr, upc);
- 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);
- if (wat_val > 0)
- wat_val = -wat_val;
- cseg_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); */
+ /* check if neighbour has not been worked on */
+ else if (in_val && !work_val) {
+ cseg_get(&asp, &asp_up, upr, upc);
+ 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) {
+ wat_val = -wat_val;
+ dseg_put(&wat, &wat_val, r, c);
}
}
}
}
- }
+ } /* end sides */
}
- bseg_close(&worked);
+
+ if (mfd == 0)
+ bseg_close(&worked);
+
bseg_close(&in_list);
seg_close(&heap_index);
@@ -98,15 +146,11 @@
}
/* 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;
- /* double slp_value; */
-
- /* slp_value = get_slope(r, c, downr, downc, ele, downe);
- dseg_put (&slp, &slp_value, r, c); */
bseg_put(&in_list, &one, r, c);
/* add point to next free position */
@@ -122,11 +166,11 @@
point.r = r;
point.c = c;
- 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 */
@@ -140,8 +184,8 @@
int drop_pt(void)
{
int child, childr, parent;
- int childp, childrp;
- CELL ele, eler;
+ int childp;
+ CELL ele;
int i;
HEAP heap_pos;
@@ -168,23 +212,20 @@
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) {
+ 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++;
}
}
@@ -205,7 +246,6 @@
/* sift up last swapped point, only necessary if hole moved to heap end */
sift_up(parent, ele);
-
}
/* the actual drop */
@@ -218,8 +258,7 @@
/* 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;
HEAP heap_pos;
child = start;
@@ -229,24 +268,20 @@
while (child > 1) {
parent = GET_PARENT(child);
seg_get(&heap_index, (char *)&heap_pos, 0, parent);
- 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);
child = parent;
-
}
/* 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);
child = parent;
-
}
else
/* no more sifting up, found new slot for child */
@@ -278,27 +313,32 @@
return (slope);
}
-int replace( /* ele was in there */
- SHORT upr, SHORT upc, SHORT r, SHORT c)
-/* CELL ele; */
+double get_slope2(CELL ele, CELL up_ele, double dist)
{
- int now, heap_run;
+ if (ele >= up_ele)
+ return 0.0;
+ else
+ return (double)(up_ele - ele) / dist;
+}
+
+/* no longer needed */
+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;
- 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/branches/releasebranch_6_4/raster/r.watershed/seg/do_astar.h
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_astar.h 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_astar.h 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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/branches/releasebranch_6_4/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_cum.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/do_cum.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -7,11 +7,14 @@
int do_cum(void)
{
SHORT r, c, dr, dc;
- CELL is_swale, value, valued;
+ 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."));
+ G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
count = 0;
if (bas_thres <= 0)
@@ -19,18 +22,24 @@
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;
- if ((dr = point.downr) > -1) {
- r = point.r;
- c = point.c;
- dc = point.downc;
- cseg_get(&wat, &value, r, c);
+ 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);
- cseg_get(&wat, &valued, dr, dc);
+ dseg_get(&wat, &valued, dr, dc);
if (value > 0) {
if (valued > 0)
valued += value;
@@ -43,8 +52,16 @@
else
valued = value - valued;
}
- cseg_put(&wat, &valued, dr, dc);
+ 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);
}
@@ -56,6 +73,341 @@
}
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, *wat_nbr;
+ POINT point;
+ int killer, threshold, count;
+
+ /* MFD */
+ int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
+ double *dist_to_nbr, *weight, sum_weight, max_weight;
+ int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
+ double dx, dy;
+ CELL ele, ele_nbr, asp_val, asp_val2, cvalue, *worked_nbr;
+ double prop, max_acc;
+ int workedon, edge, flat;
+ 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));
+ weight = (double *)G_malloc(sides * sizeof(double));
+
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ r_nbr = nextdr[ct_dir];
+ c_nbr = nextdc[ct_dir];
+ /* account for rare cases when ns_res != ew_res */
+ dy = ABS(r_nbr) * window.ns_res;
+ dx = ABS(c_nbr) * window.ew_res;
+ if (ct_dir < 4)
+ dist_to_nbr[ct_dir] = dx + dy;
+ else
+ dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+ }
+
+ /* reset worked, takes time... */
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ bseg_put(&worked, &zero, r, c);
+ }
+ }
+
+ worked_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
+ wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
+
+ 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;
+ 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);
+
+ 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;
+ cseg_get(&alt, &ele, r, c);
+ is_null = 0;
+ edge = 0;
+ flat = 1;
+ /* 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 */
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+ weight[ct_dir] = -1;
+ worked_nbr[ct_dir] = 1;
+ /* 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 */
+ 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;
+ cseg_get(&alt, &ele_nbr, r_nbr, c_nbr);
+ if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold &&
+ ct_dir != np_side && ele_nbr > ele)
+ stream_cells++;
+
+ bseg_get(&worked, &cvalue, r_nbr, c_nbr);
+ worked_nbr[ct_dir] = cvalue;
+ if (worked_nbr[ct_dir] == 0) {
+ is_null = G_is_c_null_value(&ele_nbr);
+ edge = is_null;
+ if (ele_nbr != ele)
+ flat = 0;
+ 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) {
+ 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;
+ }
+ else
+ edge = 1;
+ if (edge)
+ break;
+ }
+ /* do not distribute flow along edges, this causes artifacts */
+ if (edge) {
+ bseg_get(&swale, &is_swale, r, c);
+ if (is_swale && asp_val > 0) {
+ asp_val = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+ cseg_put(&asp, &asp_val, r, c);
+ }
+ continue;
+ }
+
+ /* 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 > 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;
+ }
+
+ /* set flow accumulation for neighbours */
+ max_acc = -1;
+
+ 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, &is_worked, r_nbr, c_nbr); */
+ if (worked_nbr[ct_dir] == 0) {
+
+ weight[ct_dir] = weight[ct_dir] / sum_weight;
+ /* check everything sums up to 1.0 */
+ prop += weight[ct_dir];
+
+ if (value > 0) {
+ if (wat_nbr[ct_dir] > 0)
+ wat_nbr[ct_dir] += value * weight[ct_dir];
+ else
+ wat_nbr[ct_dir] -= value * weight[ct_dir];
+ }
+ else {
+ if (wat_nbr[ct_dir] < 0)
+ wat_nbr[ct_dir] += value * weight[ct_dir];
+ else
+ 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);
+
+ /* 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;
+ }
+ }
+ else if (ct_dir == np_side) {
+ /* check for consistency with A * path */
+ workedon++;
+ }
+ }
+ }
+ if (ABS(prop - 1.0) > 5E-6f) {
+ G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
+ prop);
+ }
+ }
+
+ if (mfd_cells < 2) {
+ dseg_get(&wat, &valued, dr, dc);
+ 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);
+ }
+
+ /* 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);
+
+ }
+ /* update asp for depression */
+ bseg_get(&swale, &is_swale, r, c);
+ 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 && !flat) {
+ bseg_put(&swale, &one, r, c);
+ is_swale = 1;
+ }
+ /* continue stream */
+ if (is_swale) {
+ bseg_put(&swale, &one, r_max, c_max);
+ }
+ else {
+ if (er_flag && !is_swale)
+ slope_length(r, c, r_max, c_max);
+ }
+ 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 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;
+}
+
+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/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_get.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_get.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_get.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_put.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_put.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_put.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_read.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_read.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_read.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_write.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_write.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/dseg_write.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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/branches/releasebranch_6_4/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/init_vars.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/init_vars.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -1,22 +1,28 @@
#include <stdlib.h>
+#include <string.h>
#include <unistd.h>
#include "Gwater.h"
#include <grass/gis.h>
#include <grass/glocale.h>
+int ele_round(double);
int init_vars(int argc, char *argv[])
{
SHORT r, c;
- int fd, num_cseg_total, num_open_segs;
+ int fd, num_cseg_total, num_open_segs, n_array_segs;
int seg_rows, seg_cols;
double segs_mb;
/* int page_block, num_cseg; */
int max_bytes;
- CELL *buf, alt_value, wat_value, asp_value, worked_value;
- extern FILE *fopen();
- char MASK_flag, *do_exist();
+ CELL *buf, alt_value, asp_value, worked_value, block_value;
+ DCELL wat_value;
+ DCELL dvalue;
+ char MASK_flag;
+ void *elebuf, *ptr;
+ int ele_map_type;
+ size_t ele_size;
G_gisinit(argv[0]);
ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
@@ -24,11 +30,16 @@
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;
+ abs_acc = 0;
+ ele_scale = 1;
segs_mb = 300;
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
@@ -73,9 +84,17 @@
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 if (strcmp(argv[r], "-a") == 0)
+ abs_acc = 1;
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) &&
@@ -96,27 +115,11 @@
tot_parts);
this_mapset = G_mapset();
- if (asp_flag)
- do_legal(asp_name);
- if (bas_flag)
- do_legal(bas_name);
- if (seg_flag)
- do_legal(seg_name);
- if (haf_flag)
- do_legal(haf_name);
- if (sl_flag)
- do_legal(sl_name);
- if (sg_flag)
- do_legal(sg_name);
- if (ls_flag)
- do_legal(ls_name);
if (sl_flag || sg_flag || ls_flag)
er_flag = 1;
- ele_mapset = do_exist(ele_name);
/* for sd factor
if (dep_flag) {
if (sscanf (dep_name, "%lf", &dep_slope) != 1) {
- dep_mapset = do_exist (dep_name);
dep_flag = -1;
}
}
@@ -124,7 +127,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;
@@ -137,8 +140,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;
@@ -148,7 +150,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);
@@ -172,26 +174,124 @@
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, ele_mapset);
- cseg_read_cell(&r_h, ele_name, ele_mapset);
- cseg_open(&wat, seg_rows, seg_cols, num_open_segs);
+ cseg_read_cell(&alt, ele_name, "");
+ if (er_flag) {
+ cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
+ cseg_read_cell(&r_h, ele_name, "");
+ }
+
+ /* read elevation input and mark NULL/masked cells */
+ bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
+ bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
+ G_verbose_message("Checking for masked and NULL cells in input elevation <%s>", ele_name);
+ /* open elevation input */
+ fd = G_open_cell_old(ele_name, "");
+ if (fd < 0) {
+ G_fatal_error(_("unable to open elevation map layer"));
+ }
+
+ ele_map_type = G_get_raster_map_type(fd);
+ ele_size = G_raster_size(ele_map_type);
+ elebuf = G_allocate_raster_buf(ele_map_type);
+
+ if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
+ ele_scale = 1000; /* should be enough to do the trick */
+
+ /* read elevation input and mark NULL/masked cells */
+ MASK_flag = 0;
+ do_points = nrows * ncols;
+ for (r = 0; r < nrows; r++) {
+ G_get_raster_row(fd, elebuf, r, ele_map_type);
+ ptr = elebuf;
+ for (c = 0; c < ncols; c++) {
+
+ /* check for masked and NULL cells */
+ if (G_is_null_value(ptr, ele_map_type)) {
+ bseg_put(&worked, &one, r, c);
+ bseg_put(&in_list, &one, r, c);
+ G_set_c_null_value(&alt_value, 1);
+ do_points--;
+ }
+ else {
+ if (ele_map_type == CELL_TYPE) {
+ alt_value = *((CELL *)ptr);
+ }
+ else if (ele_map_type == FCELL_TYPE) {
+ dvalue = *((FCELL *)ptr);
+ dvalue *= ele_scale;
+ alt_value = ele_round(dvalue);
+ }
+ else if (ele_map_type == DCELL_TYPE) {
+ dvalue = *((DCELL *)ptr);
+ dvalue *= ele_scale;
+ alt_value = ele_round(dvalue);
+ }
+ }
+ cseg_put(&alt, &alt_value, r, c);
+ if (er_flag) {
+ cseg_put(&r_h, &alt_value, r, c);
+ }
+ ptr = G_incr_void_ptr(ptr, ele_size);
+ }
+ }
+ G_close_cell(fd);
+ G_free(elebuf);
+ if (do_points < nrows * ncols)
+ MASK_flag = 1;
+
+ /* initial flow accumulation */
+ dseg_open(&wat, seg_rows, seg_cols, num_open_segs);
if (run_flag) {
- run_mapset = do_exist(run_name);
- cseg_read_cell(&wat, run_name, run_mapset);
+ 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 == cseg_put(&wat, &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) {
- pit_mapset = do_exist(pit_name);
- cseg_read_cell(&asp, pit_name, pit_mapset);
+ 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) && 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++) {
@@ -202,8 +302,25 @@
}
bseg_open(&swale, seg_rows, seg_cols, num_open_segs);
if (ob_flag) {
- ob_mapset = do_exist(ob_name);
- bseg_read_cell(&swale, ob_name, ob_mapset);
+ 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) && 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++) {
@@ -212,49 +329,40 @@
}
}
if (ril_flag) {
- ril_mapset = do_exist(ril_name);
dseg_open(&ril, 1, seg_rows * seg_cols, num_open_segs);
- dseg_read_cell(&ril, ril_name, ril_mapset);
+ 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));
- /* heap_index will track astar_pts in the binary min-heap */
+ if (num_open_segs / 2 > 0)
+ n_array_segs = num_open_segs / 2;
+ else
+ n_array_segs = 1;
+
+ seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols * 2,
+ n_array_segs, 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));
+ if (seg_cols * num_open_segs * seg_rows / 10 > 0)
+ n_array_segs = seg_cols * num_open_segs * seg_rows / 10;
+ else
+ n_array_segs = 1;
+ seg_open(&heap_index, 1, do_points + 1, 1, n_array_segs,
+ 10, sizeof(HEAP));
+
G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
/* heap is empty */
@@ -264,23 +372,32 @@
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);
+ 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) {
- 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)
+ /* 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;
@@ -288,111 +405,109 @@
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);
- 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)
&& 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);
- 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)
&& 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);
- 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)
&& 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);
- 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 &&
!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);
- 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 &&
!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);
- 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 &&
!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);
- 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 &&
!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);
- 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 {
@@ -405,20 +520,29 @@
}
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);
+ 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) {
- 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)
+ /* 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;
@@ -426,12 +550,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);
@@ -440,25 +562,21 @@
}
}
}
- G_percent(r, nrows, 3); /* finish it */
+ G_percent(r, nrows, 1); /* finish it */
return 0;
}
-int do_legal(char *file_name)
+int ele_round(double x)
{
- if (G_legal_filename(file_name) == -1)
- G_fatal_error(_("<%s> is an illegal file name"), file_name);
+ int n;
- return 0;
-}
+ if (x >= 0.0)
+ n = x + .5;
+ else {
+ n = -x + .5;
+ n = -n;
+ }
-char *do_exist(char *file_name)
-{
- char *file_mapset = G_find_cell2(file_name, "");
-
- if (file_mapset == NULL)
- G_fatal_error(_("Raster map <%s> not found"), file_name);
-
- return (file_mapset);
+ return n;
}
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/seg/main.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/main.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/main.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -8,33 +8,74 @@
* Brad Douglas <rez touchofmadness.com>,
* Hamish Bowman <hamish_b yahoo.com>,
* Markus Metz <markus.metz.giswork gmail.com>
- * PURPOSE: Watershed determination using the GRASS segmentation lib
- * COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
+ * PURPOSE: Hydrological analysis using the GRASS segmentation lib
+ * COPYRIGHT: (C) 1999-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
* for details.
*
*****************************************************************************/
-#define MAIN
+
#include <stdlib.h>
#include <unistd.h>
#include "Gwater.h"
#include <grass/gis.h>
#include <grass/glocale.h>
-#undef MAIN
+struct Cell_head window;
+int mfd, c_fac, abs_acc, ele_scale;
+SSEG heap_index;
+int heap_size;
+int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
+SHORT nrows, ncols;
+double half_res, diag, max_length, dep_slope;
+int bas_thres, tot_parts;
+SSEG astar_pts;
+BSEG worked, in_list, s_b, swale;
+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, 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 };
+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 ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
+char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
+char bas_flag, seg_flag, haf_flag, er_flag;
+char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
+FILE *fp;
+
+
+
int main(int argc, char *argv[])
{
- extern FILE *fopen();
-
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/branches/releasebranch_6_4/raster/r.watershed/seg/no_stream.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/no_stream.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/no_stream.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -6,10 +6,12 @@
{
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, value;
+ DCELL dvalue, max_drain; /* flow acc is now DCELL */
SHORT updir, riteflag, leftflag, thisdir;
while (1) {
+ cseg_put(&bas, &basin_num, row, col);
max_drain = -1;
for (r = row - 1, rr = 0; r <= row + 1; r++, rr++) {
for (c = col - 1, cc = 0; c <= col + 1; c++, cc++) {
@@ -17,13 +19,13 @@
cseg_get(&asp, &aspect, r, c);
if (aspect == drain[rr][cc]) {
- cseg_get(&wat, &value, r, c);
- if (value < 0)
- value = -value;
- if (value > max_drain) {
+ dseg_get(&wat, &dvalue, r, c);
+ if (dvalue < 0)
+ dvalue = -dvalue;
+ if ((dvalue - max_drain) > 5E-8f) { /* floating point comparison problem workaround */
uprow = r;
upcol = c;
- max_drain = value;
+ max_drain = dvalue;
}
}
}
@@ -66,14 +68,15 @@
cseg_get(&asp, &aspect, r, c);
if (aspect == drain[rr][cc]) {
thisdir = updrain[rr][cc];
- if (haf_basin_side(updir,
+ switch (haf_basin_side(updir,
(SHORT) downdir,
- thisdir) == RITE) {
+ thisdir)) {
+ case RITE:
overland_cells(r, c, basin_num, basin_num,
&new_ele);
riteflag++;
- }
- else {
+ break;
+ case LEFT:
overland_cells(r, c, basin_num, basin_num - 1,
&new_ele);
leftflag++;
@@ -82,7 +85,7 @@
}
}
}
- if (leftflag >= riteflag) {
+ if (leftflag > riteflag) {
value = basin_num - 1;
cseg_put(&haf, &value, row, col);
}
@@ -99,6 +102,7 @@
slope = MIN_SLOPE;
fprintf(fp, " %f %f\n", slope, stream_length);
}
+ cseg_put(&haf, &basin_num, row, col);
return 0;
}
}
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/seg/sg_factor.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/sg_factor.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/sg_factor.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -9,16 +9,15 @@
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);
dseg_get(&s_l, &length, r, c);
cseg_get(&asp, &downer, r, c);
- height = hih_elev - low_elev;
+ height = 1.0 * (hih_elev - low_elev) / ele_scale;
if (length > max_length) {
height *= max_length / length;
length = max_length;
@@ -33,12 +32,12 @@
len_slp_equ(length, sin_theta, S, r, c);
}
if (sg_flag) {
- S *= 100.0;
+ S *= 100.0; /* factor 100 for backwards compatibility */
dseg_put(&s_g, &S, r, c);
}
}
}
- G_percent(r, nrows, 3); /* finish it */
+ G_percent(nrows, nrows, 1); /* finish it */
return 0;
}
@@ -61,7 +60,7 @@
/* rill_ratio equation from Steve Warren */
rill_ratio *= .5 + .005 * rill + .0001 * rill * rill;
s_l_exp = rill_ratio / (1 + rill_ratio);
- LS = S * 100 * pow((slope_length / 72.6), s_l_exp);
+ LS = S * 100 * pow((slope_length / 72.6), s_l_exp); /* factor 100 for backwards compatibility */
dseg_put(&l_s, &LS, r, c);
return 0;
Modified: grass/branches/releasebranch_6_4/raster/r.watershed/seg/slope_len.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.watershed/seg/slope_len.c 2010-05-05 13:25:16 UTC (rev 42120)
+++ grass/branches/releasebranch_6_4/raster/r.watershed/seg/slope_len.c 2010-05-05 13:30:34 UTC (rev 42121)
@@ -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