[GRASS-SVN] r35905 - in
grass/branches/develbranch_6/raster/r.watershed: front ram seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Feb 17 07:16:36 EST 2009
Author: mmetz
Date: 2009-02-17 07:16:36 -0500 (Tue, 17 Feb 2009)
New Revision: 35905
Modified:
grass/branches/develbranch_6/raster/r.watershed/front/description.html
grass/branches/develbranch_6/raster/r.watershed/front/main.c
grass/branches/develbranch_6/raster/r.watershed/ram/Gwater.h
grass/branches/develbranch_6/raster/r.watershed/ram/close_maps.c
grass/branches/develbranch_6/raster/r.watershed/ram/close_maps2.c
grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c
grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.h
grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c
grass/branches/develbranch_6/raster/r.watershed/ram/find_pour.c
grass/branches/develbranch_6/raster/r.watershed/ram/init_vars.c
grass/branches/develbranch_6/raster/r.watershed/ram/main.c
grass/branches/develbranch_6/raster/r.watershed/ram/no_stream.c
grass/branches/develbranch_6/raster/r.watershed/ram/sg_factor.c
grass/branches/develbranch_6/raster/r.watershed/ram/slope_len.c
grass/branches/develbranch_6/raster/r.watershed/seg/Gwater.h
grass/branches/develbranch_6/raster/r.watershed/seg/close_maps.c
grass/branches/develbranch_6/raster/r.watershed/seg/close_maps2.c
grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c
grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.h
grass/branches/develbranch_6/raster/r.watershed/seg/do_cum.c
grass/branches/develbranch_6/raster/r.watershed/seg/dseg_get.c
grass/branches/develbranch_6/raster/r.watershed/seg/dseg_put.c
grass/branches/develbranch_6/raster/r.watershed/seg/dseg_read.c
grass/branches/develbranch_6/raster/r.watershed/seg/dseg_write.c
grass/branches/develbranch_6/raster/r.watershed/seg/init_vars.c
grass/branches/develbranch_6/raster/r.watershed/seg/main.c
grass/branches/develbranch_6/raster/r.watershed/seg/no_stream.c
grass/branches/develbranch_6/raster/r.watershed/seg/sg_factor.c
grass/branches/develbranch_6/raster/r.watershed/seg/slope_len.c
Log:
backporting MFD with full backwards compatibility
Modified: grass/branches/develbranch_6/raster/r.watershed/front/description.html
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/front/description.html 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/front/description.html 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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,31 @@
surface flow (allows horizontal, vertical, and diagonal flow of water).
This flag will also make the drainage basins look more homogeneous.
+<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 +72,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 +85,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 +120,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 +139,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>
@@ -138,30 +161,28 @@
<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 +201,81 @@
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 overflown
+with a gracefull flow convergence before the overflow. The convergence
+factor causes flow accumulation to converge more strongly with higher
+values. The supported range is 1 to 10, recommended is a convergence
+factor of 5 (Holmgren, 1994). If many small sliver basins are created
+with MFD, setting the convergence factor to a higher value can reduce
+the amount of small sliver basins.
+<br>See example below with the North Carolina dataset for using MFD mode.
+
+<h4>In-memory mode and disk swap mode</h4>
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>High-resolution elevation maps with floating point values</h4>
+To get better results with high resolution elevation maps with
+floating point values, it may be necessary to multiply the original map
+with e.g. 100 to change elevation units from meters to cm, because
+<em>r.watershed</em> reads input elevation maps as integer. This allows
+control over how much detail is required to obtain useful results with a
+given dataset and region setting. Remember to adjust any <em>slope length
+and steepness (LS) factor</em> or <em>slope steepness (S) factor</em>
+output accordingly. See example below for using <em>r.mapcalc</em> to
+convert elevation from meter to millimeter.
+
+<h4>Basin threshold</h4>
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 +285,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.
@@ -265,8 +322,16 @@
<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>
@@ -359,6 +424,18 @@
</pre></div>
<br>
+<i>This example uses the North Carolina sample dataset.</i>
+<p>
+Using MFD mode on a LIDAR elevation dataset:
+<br>Convert elevation values of elev_lid972_1m from meter to millimeter
+then compare MFD mode with default medium flow convergence to SFD mode.
+<div class="code"><pre>
+ r.mapcalc "elev_lid972_1m.1mm = round(elev_lid972_1m * 1000.0)"
+ r.watershed elevation=elev_lid972_1m.1mm accumulation=elev_lid972_1m.1mm.acc drainage=elev_lid972_1m.1mm.dir convergence=5
+ r.watershed -s elevation=elev_lid972_1m.1mm accumulation=elev_lid972_1m.1mm.acc.sfd drainage=elev_lid972_1m.1mm.dir.sfd
+</pre></div>
+
+
<a name="references"></a>
<h2>REFERENCES</h2>
@@ -371,6 +448,12 @@
http://faculty.wiu.edu/CR-Ehlschlaeger2/older/IGIS/paper.html</a>
<p>
+Holmgren, P. (1994). <i>Multiple flow direction algorithms for runoff
+modelling in grid based elevation models: An empirical evaluation.</i>
+<b>Hydrological Processes</b> Vol 8(4), p.327-334.<br>
+DOI: <a href="http://dx.doi.org/10.1002/hyp.3360080405">10.1002/hyp.3360080405</a>
+
+<p>
Kinner D., H. Mitasova, R. Harmon, L. Toma, R., Stallard. (2005).
<i>GIS-based Stream Network Analysis for The Chagres River Basin,
Republic of Panama</i>. <b>The Rio Chagres: A Multidisciplinary Profile of
@@ -395,8 +478,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 +496,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:
+Markus Metz <markus.metz.giswork at gmail.com>
<p>
<i>Last changed: $Date$</i>
Modified: grass/branches/develbranch_6/raster/r.watershed/front/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/front/main.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/front/main.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -4,7 +4,8 @@
* MODULE: front end
* AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
* Brad Douglas <rez touchofmadness.com>,
- * Hamish Bowman <hamish_b yahoo.com>
+ * Hamish Bowman <hamish_b yahoo.com>
+ * Markus Metz <markus.metz.giswork gmail.com>
* PURPOSE: Watershed determination
* COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
*
@@ -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,6 +42,8 @@
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 GModule *module;
@@ -59,7 +62,8 @@
opt2 = G_define_option();
opt2->key = "depression";
- opt2->description = _("Input map: locations of real depressions");
+ opt2->label = _("Input map: locations of real depressions");
+ opt2->description = _("All non-NULL and non-zero cells are considered as real depressions");
opt2->required = NO;
opt2->type = TYPE_STRING;
opt2->gisprompt = "old,cell,raster";
@@ -84,8 +88,10 @@
opt5 = G_define_option();
opt5->key = "blocking";
+ opt5->label =
+ _("Input map: terrain blocking overland surface flow, for USLE");
opt5->description =
- _("Input map: terrain blocking overland surface flow, for USLE");
+ _("All non-NULL and non-zero cells are considered as blocking terrain");
opt5->required = NO;
opt5->type = TYPE_STRING;
opt5->gisprompt = "old,cell,raster";
@@ -178,13 +184,28 @@
opt15->gisprompt = "new,cell,raster";
opt15->guisection = _("Output_options");
- opt16 = G_define_option() ;
- opt16->key = "memory";
- opt16->type = TYPE_INTEGER;
- opt16->required = NO;
- opt16->answer = "300"; /* 300MB default value, please keep in sync with r.terraflow */
- opt16->description = _("Maximum memory to be used with -m flag (in MB)");
+ opt16 = G_define_option();
+ opt16->key = "convergence";
+ opt16->type = TYPE_INTEGER;
+ opt16->required = NO;
+ opt16->answer = "5";
+ opt16->label = _("Convergence factor for MFD (1-10)");
+ opt16->description =
+ _("1 = most diverging flow, 10 = most converging flow. Recommended: 5");
+ opt17 = G_define_option();
+ opt17->key = "memory";
+ opt17->type = TYPE_INTEGER;
+ opt17->required = NO;
+ opt17->answer = "300"; /* 300MB default value, please keep in sync with r.terraflow */
+ opt17->description = _("Maximum memory to be used with -m flag (in MB)");
+
+ flag_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 =
@@ -214,18 +235,20 @@
}
err = 0;
- /* basin and basin.thresh */
+ /* basin and basin threshold */
err += (opt10->answer != NULL && opt6->answer == NULL);
- /* stream and basin.thresh */
+ /* stream and basin threshold */
err += (opt11->answer != NULL && opt6->answer == NULL);
- /* half.basin and basin.thresh */
+ /* half.basin and basin threshold */
err += (opt12->answer != NULL && opt6->answer == NULL);
- /* slope and basin.thresh */
+ /* LS factor and basin threshold */
+ err += (opt14->answer != NULL && opt6->answer == NULL);
+ /* S factor and basin threshold */
err += (opt15->answer != NULL && opt6->answer == NULL);
if (err) {
G_message(_("Sorry, if any of the following options are set:\n"
- " basin, stream, half.basin, slope, or lS\n"
+ " basin, stream, half.basin, length.slope, or slope.steepness\n"
" you MUST provide a value for the basin "
"threshold parameter."));
G_usage();
@@ -240,6 +263,10 @@
else
strcat(command, "r.watershed.ram");
+ if (!flag_mfd->answer) {
+ strcat(command, " -s");
+ }
+
if (flag_flow->answer)
strcat(command, " -4");
@@ -344,57 +371,63 @@
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);
+ }
+
G_debug(1, "Mode: %s", flag_seg->answer ? "Segmented" : "All in RAM");
G_debug(1, "Running: %s", command);
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;
@@ -404,8 +437,10 @@
strncpy(history.datsrc_1, source_name, RECORD_LEN);
history.datsrc_1[RECORD_LEN - 1] = '\0'; /* strncpy() doesn't null terminate if maxfill */
sprintf(history.edhist[0],
- "Processing mode: %s", mode ? "Segmented" : "All in RAM");
- history.edlinecnt = 1;
+ "Processing mode: %s", 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/develbranch_6/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/Gwater.h 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/Gwater.h 2009-02-17 12:16:36 UTC (rev 35905)
@@ -34,65 +34,48 @@
#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;
+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 POINT *astar_pts;
+extern CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
+extern DCELL *wat;
+extern CELL *ril_buf;
+extern int ril_fd;
+extern double *s_l, *s_g, *l_s;
+extern CELL one, zero;
+extern double ril_value, 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 +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,8 +106,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);
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/close_maps.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/close_maps.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/close_maps.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -9,30 +9,114 @@
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++) {
+ G_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
for (c = 0; c < ncols; c++) {
- buf[c] = wat[SEG_INDEX(wat_seg, r, c)];
+ dvalue = wat[SEG_INDEX(wat_seg, r, c)];
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dbuf[c] = dvalue;
+ dvalue = ABS(dvalue);
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
}
- G_put_raster_row(fd, buf, CELL_TYPE);
+ 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 +124,7 @@
}
else {
for (r = 0; r < nrows; r++) {
+ G_set_c_null_value(buf, ncols); /* reset row to all NULL */
for (c = 0; c < ncols; c++) {
buf[c] = asp[SEG_INDEX(asp_seg, r, c)];
}
@@ -54,6 +139,7 @@
}
G_free(asp);
+ /* visual output no longer needed */
if (dis_flag) {
fd = G_open_cell_new(dis_name);
if (fd < 0) {
@@ -91,16 +177,16 @@
G_free(wat);
if (ls_flag) {
- fd = G_open_cell_new(ls_name);
+ fd = G_open_raster_new(ls_name, DCELL_TYPE);
if (fd < 0) {
- G_warning(_("unable to open new L map layer."));
+ G_warning(_("unable to open new LS factor map layer."));
}
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++) {
- buf[c] = l_s[SEG_INDEX(l_s_seg, r, c)] + .5;
+ dbuf[c] = l_s[SEG_INDEX(l_s_seg, r, c)];
}
- G_put_raster_row(fd, buf, CELL_TYPE);
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
if (G_close_cell(fd) < 0)
G_warning(_("Close failed."));
@@ -109,18 +195,18 @@
}
if (sl_flag) {
- fd = G_open_cell_new(sl_name);
+ fd = G_open_raster_new(sl_name, DCELL_TYPE);
if (fd < 0) {
G_warning(_("unable to open new slope length map layer."));
}
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++) {
- buf[c] = s_l[SEG_INDEX(s_l_seg, r, c)] + .5;
- if (buf[c] > max_length)
- buf[c] = max_length + .5;
+ dbuf[c] = s_l[SEG_INDEX(s_l_seg, r, c)];
+ if (dbuf[c] > max_length)
+ dbuf[c] = max_length;
}
- G_put_raster_row(fd, buf, CELL_TYPE);
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
if (G_close_cell(fd) < 0)
G_warning(_("Close failed."));
@@ -131,16 +217,16 @@
G_free(s_l);
if (sg_flag) {
- fd = G_open_cell_new(sg_name);
+ fd = G_open_raster_new(sg_name, DCELL_TYPE);
if (fd < 0) {
- G_warning(_("unable to open new S map layer."));
+ G_warning(_("unable to open new S factor map layer."));
}
else {
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++) {
- buf[c] = s_g[SEG_INDEX(s_g_seg, r, c)] * 100 + .5;
+ dbuf[c] = s_g[SEG_INDEX(s_g_seg, r, c)];
}
- G_put_raster_row(fd, buf, CELL_TYPE);
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
if (G_close_cell(fd) < 0)
G_warning(_("Close failed."));
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/close_maps2.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/close_maps2.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/close_maps2.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -3,8 +3,6 @@
#include <grass/gis.h>
#include <grass/glocale.h>
-int sift_up(int, CELL);
-
int do_astar(void)
{
POINT *point;
@@ -38,7 +36,7 @@
/* equivalent to first_astar = point->next in old code */
drop_pt();
- /* can go, dragged on from old code */
+ /* can go, dragged on from old code: replace first_astar with heap_index[1] in line 22 */
first_astar = heap_index[1];
/* downhill path for flow accumulation is set here */
@@ -53,7 +51,6 @@
index_doer = SEG_INDEX(alt_seg, r, c);
alt_val = alt[index_doer];
- wat_val = wat[index_doer];
FLAG_SET(worked, r, c);
@@ -71,24 +68,24 @@
if (in_val == 0) {
alt_up = alt[index_up];
/* flow direction is set here */
- add_pt(upr, upc, r, c, alt_up, alt_val);
+ add_pt(upr, upc, alt_up, alt_val);
drain_val = drain[upr - r + 1][upc - c + 1];
asp[index_up] = drain_val;
-
}
else {
- /* check if neighbour has not been worked on,
- * update values for asp and wat */
+ /* check if neighbour has been worked on,
+ * if not, update values for asp and wat */
in_val = FLAG_GET(worked, upr, upc);
if (in_val == 0) {
asp_up = asp[index_up];
- if (asp_up < -1) {
+ if (asp_up < 0) {
asp[index_up] = drain[upr - r + 1][upc - c + 1];
+ wat_val = wat[index_doer];
if (wat_val > 0)
wat[index_doer] = -wat_val;
- replace(upr, upc, r, c); /* alt_up used to be */
+ /* replace(upr, upc, r, c); */ /* alt_up used to be */
}
}
}
@@ -96,7 +93,9 @@
}
}
G_percent(count, do_points, 3); /* finish it */
- flag_destroy(worked);
+ if (mfd == 0)
+ flag_destroy(worked);
+
flag_destroy(in_list);
G_free(heap_index);
@@ -104,7 +103,7 @@
}
/* new add point routine for min heap */
-int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
{
FLAG_SET(in_list, r, c);
@@ -120,8 +119,8 @@
astar_pts[nxt_avail_pt].r = r;
astar_pts[nxt_avail_pt].c = c;
- astar_pts[nxt_avail_pt].downr = downr;
- astar_pts[nxt_avail_pt].downc = downc;
+/* astar_pts[nxt_avail_pt].downr = downr;
+ astar_pts[nxt_avail_pt].downc = downc; */
nxt_avail_pt++;
@@ -273,7 +272,7 @@
}
-/* new replace */
+/* no longer needed */
int replace( /* ele was in there */
SHORT upr, SHORT upc, SHORT r, SHORT c)
/* CELL ele; */
@@ -288,8 +287,8 @@
while (heap_run <= heap_size) {
now = heap_index[heap_run];
if (astar_pts[now].r == upr && astar_pts[now].c == upc) {
- astar_pts[now].downr = r;
- astar_pts[now].downc = c;
+ /*astar_pts[now].downr = r;
+ astar_pts[now].downc = c; */
return 0;
}
heap_run++;
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.h
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.h 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.h 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -6,10 +6,12 @@
int do_cum(void)
{
SHORT r, c, dr, dc;
- CELL is_swale, value, valued;
+ CELL is_swale, value, valued, aspect;
int killer, threshold, count;
+ SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
- G_message(_("SECTION 3: Accumulating Surface Flow."));
+ G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
count = 0;
if (bas_thres <= 0)
@@ -17,15 +19,21 @@
else
threshold = bas_thres;
while (first_cum != -1) {
- G_percent(count++, do_points, 3);
+ G_percent(count++, do_points, 2);
killer = first_cum;
first_cum = astar_pts[killer].nxt;
- if ((dr = astar_pts[killer].downr) > -1) {
- r = astar_pts[killer].r;
- c = astar_pts[killer].c;
- dc = astar_pts[killer].downc;
+ r = astar_pts[killer].r;
+ c = astar_pts[killer].c;
+ aspect = asp[SEG_INDEX(asp_seg, r, c)];
+ if (aspect) {
+ dr = r + asp_r[ABS(aspect)];
+ dc = c + asp_c[ABS(aspect)];
+ }
+ else
+ dr = dc = -1;
+ if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
value = wat[SEG_INDEX(wat_seg, r, c)];
- if (ABS(value) >= threshold)
+ if ((int)(ABS(value) + 0.5) >= threshold)
FLAG_SET(swale, r, c);
valued = wat[SEG_INDEX(wat_seg, dr, dc)];
if (value > 0) {
@@ -41,8 +49,16 @@
valued = value - valued;
}
wat[SEG_INDEX(wat_seg, dr, dc)] = 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[SEG_INDEX(asp_seg, dr, dc)] == 0) {
+ aspect = -aspect;
+ asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ }
+ }
+ if (is_swale || ((int)valued) >= threshold) {
FLAG_SET(swale, dr, dc);
}
else {
@@ -51,8 +67,320 @@
}
}
}
- G_percent(count, do_points, 3); /* finish it */
+ G_percent(count, do_points, 1); /* finish it */
G_free(astar_pts);
return 0;
}
+
+/***************************************
+ *
+ * MFD references
+ *
+ * original:
+ * Quinn, P., Beven, K., Chevallier, P., and Planchon, 0. 1991.
+ * The prediction of hillslope flow paths for distributed hydrological
+ * modelling using digital terrain models, Hydrol. Process., 5, 59-79.
+ *
+ * modified by Holmgren (1994):
+ * Holmgren, P. 1994. Multiple flow direction algorithms for runoff
+ * modelling in grid based elevation models: an empirical evaluation
+ * Hydrol. Process., 8, 327-334.
+ *
+ * implemented here:
+ * Holmgren (1994) with modifications to honour A * path in order to get
+ * out of depressions and across obstacles with gracefull flow convergence
+ * before depressions/obstacles and gracefull flow divergence after
+ * depressions/obstacles
+ *
+ * ************************************/
+
+int do_cum_mfd(void)
+{
+ int r, c, dr, dc;
+ CELL is_swale;
+ DCELL value, valued;
+ 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;
+ 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);
+ }
+
+ flag_clear_all(worked);
+ workedon = 0;
+
+ count = 0;
+ if (bas_thres <= 0)
+ threshold = 60;
+ else
+ threshold = bas_thres;
+
+ while (first_cum != -1) {
+ G_percent(count++, do_points, 2);
+ killer = first_cum;
+ first_cum = astar_pts[killer].nxt;
+ r = astar_pts[killer].r;
+ c = astar_pts[killer].c;
+ aspect = asp[SEG_INDEX(asp_seg, r, c)];
+ if (aspect) {
+ dr = r + asp_r[ABS(aspect)];
+ dc = c + asp_c[ABS(aspect)];
+ }
+ else
+ dr = dc = -1;
+ if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
+ value = wat[SEG_INDEX(wat_seg, r, c)];
+ valued = wat[SEG_INDEX(wat_seg, dr, dc)];
+
+ r_max = dr;
+ c_max = dc;
+
+ /* get weights */
+ max_weight = 0;
+ sum_weight = 0;
+ np_side = -1;
+ mfd_cells = 0;
+ stream_cells = 0;
+ swale_cells = 0;
+ astar_not_set = 1;
+ ele = alt[SEG_INDEX(alt_seg, r, c)];
+ is_null = 0;
+ edge = 0;
+ /* this loop is needed to get the sum of weights */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ 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) {
+
+ /* check for swale or stream cells */
+ is_swale = FLAG_GET(swale, r_nbr, c_nbr);
+ if (is_swale)
+ swale_cells++;
+ valued = wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)];
+ if ((ABS(valued) + 0.5) >= threshold)
+ stream_cells++;
+
+ is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ if (is_worked == 0) {
+ ele_nbr = alt[SEG_INDEX(alt_seg, r_nbr, c_nbr)];
+ is_null = G_is_c_null_value(&ele_nbr);
+ edge = is_null;
+ 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[SEG_INDEX(asp_seg, r, c)] = 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) {
+
+ weight[ct_dir] = weight[ct_dir] / sum_weight;
+ /* check everything sums up to 1.0 */
+ prop += weight[ct_dir];
+
+ valued = wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)];
+ if (value > 0) {
+ if (valued > 0)
+ valued += value * weight[ct_dir];
+ else
+ valued -= value * weight[ct_dir];
+ }
+ else {
+ if (valued < 0)
+ valued += value * weight[ct_dir];
+ else
+ valued = value * weight[ct_dir] - valued;
+ }
+ wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)] = valued;
+
+ /* 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[SEG_INDEX(wat_seg, dr, dc)];
+ if (value > 0) {
+ if (valued > 0)
+ valued += value;
+ else
+ valued -= value;
+ }
+ else {
+ if (valued < 0)
+ valued += value;
+ else
+ valued = value - valued;
+ }
+ wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
+ }
+
+ /* update asp */
+ if (dr != r_max || dc != c_max) {
+ aspect = drain[r - r_max + 1][c - c_max + 1];
+ if (asp[SEG_INDEX(asp_seg, r, c)] < 0)
+ aspect = -aspect;
+ asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ }
+ is_swale = FLAG_GET(swale, r, c);
+ /* update asp for depression */
+ if (is_swale && pit_flag) {
+ if (aspect > 0 && asp[SEG_INDEX(asp_seg, r_max, c_max)] == 0) {
+ aspect = -aspect;
+ asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ }
+ }
+ /* start new stream */
+ value = ABS(value) + 0.5;
+ if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
+ swale_cells < 1) {
+ FLAG_SET(swale, 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);
+ }
+ }
+ 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);
+
+ 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/develbranch_6/raster/r.watershed/ram/find_pour.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/find_pour.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/find_pour.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/init_vars.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/init_vars.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -1,4 +1,5 @@
#include <stdlib.h>
+#include <string.h>
#include "Gwater.h"
#include <grass/gis.h>
#include <grass/glocale.h>
@@ -7,7 +8,7 @@
int init_vars(int argc, char *argv[])
{
SHORT r, c;
- CELL *buf, alt_value, wat_value, asp_value;
+ CELL *buf, alt_value, wat_value, asp_value, block_value;
int fd, index;
char MASK_flag;
@@ -15,12 +16,16 @@
ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
ob_flag = bas_flag = seg_flag = haf_flag = arm_flag = dis_flag = 0;
zero = sl_flag = sg_flag = ls_flag = er_flag = bas_thres = 0;
+ one = 1;
nxt_avail_pt = 0;
/* dep_flag = 0; */
- max_length = dzero = 0.0;
+ max_length = d_zero = 0.0;
+ d_one = 1.0;
ril_value = -1.0;
/* dep_slope = 0.0; */
sides = 8;
+ mfd = 1;
+ c_fac = 5;
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
ele_flag++;
@@ -63,9 +68,15 @@
if (sides != 4)
usage(argv[0]);
}
+ else if (sscanf(argv[r], "conv=%d", &c_fac) == 1) ;
+ else if (strcmp(argv[r], "-s") == 0)
+ mfd = 0;
else
usage(argv[0]);
}
+ if (mfd == 1 && (c_fac < 1 || c_fac > 10)) {
+ G_fatal_error("Convergence factor must be between 1 and 10.");
+ }
if ((ele_flag != 1)
||
((arm_flag == 1) &&
@@ -83,29 +94,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 +108,7 @@
nrows = G_window_rows();
ncols = G_window_cols();
total_cells = nrows * ncols;
- if (max_length <= dzero)
+ if (max_length <= d_zero)
max_length = 10 * nrows * window.ns_res + 10 * ncols * window.ew_res;
if (window.ew_res < window.ns_res)
half_res = .5 * window.ew_res;
@@ -131,10 +121,12 @@
buf = G_allocate_cell_buf();
alt =
(CELL *) G_malloc(sizeof(CELL) * size_array(&alt_seg, nrows, ncols));
- r_h =
- (CELL *) G_malloc(sizeof(CELL) * size_array(&r_h_seg, nrows, ncols));
+ if (er_flag) {
+ r_h =
+ (CELL *) G_malloc(sizeof(CELL) * size_array(&r_h_seg, nrows, ncols));
+ }
- fd = G_open_cell_old(ele_name, ele_mapset);
+ fd = G_open_cell_old(ele_name, "");
if (fd < 0) {
G_fatal_error(_("unable to open elevation map layer"));
}
@@ -143,105 +135,118 @@
in_list = flag_create(nrows, ncols);
worked = flag_create(nrows, ncols);
+ MASK_flag = 0;
+ do_points = nrows * ncols;
for (r = 0; r < nrows; r++) {
G_get_c_raster_row(fd, buf, r);
for (c = 0; c < ncols; c++) {
index = SEG_INDEX(alt_seg, r, c);
- alt[index] = r_h[index] = buf[c];
+ alt_value = alt[index] = buf[c];
+ if (er_flag) {
+ r_h[index] = buf[c];
+ }
/* all flags need to be manually set to zero */
- flag_unset(swale, r, c);
- flag_unset(in_list, r, c);
- flag_unset(worked, r, c);
+ flag_unset(swale, r, c);
+ flag_unset(in_list, r, c);
+ flag_unset(worked, r, c);
+ /* check for masked and NULL cells here */
+ if (G_is_c_null_value(&alt_value)) {
+ FLAG_SET(worked, r, c);
+ FLAG_SET(in_list, r, c);
+ do_points--;
+ }
}
}
G_close_cell(fd);
+ if (do_points < nrows * ncols)
+ MASK_flag = 1;
wat =
- (CELL *) G_malloc(sizeof(CELL) * size_array(&wat_seg, nrows, ncols));
+ (DCELL *) G_malloc(sizeof(DCELL) *
+ size_array(&wat_seg, nrows, ncols));
if (run_flag) {
- 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;
+ for (c = 0; c < ncols; c++) {
+ if (MASK_flag) {
+ index = FLAG_GET(worked, r, c);
+ if (!index)
+ wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
+ }
+ else
+ wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
+ }
}
}
- asp = (CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
+ 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);
}
+ /* 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);
}
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 =
@@ -254,7 +259,9 @@
sizeof(double));
}
- /* heap_index will track astar_pts in the binary min-heap */
+ astar_pts = (POINT *) G_malloc(do_points * sizeof(POINT));
+
+ /* heap_index will track astar_pts in ternary min-heap */
/* heap_index is one-based */
heap_index = (int *)G_malloc((do_points + 1) * sizeof(int));
@@ -272,14 +279,20 @@
wat[SEG_INDEX(wat_seg, r, c)] = 0;
}
else {
- s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
+ if (er_flag)
+ s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
asp_value = asp[SEG_INDEX(asp_seg, r, c)];
if (r == 0 || c == 0 || r == nrows - 1 ||
c == ncols - 1 || asp_value != 0) {
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
- if (r == 0)
+ /* set depression */
+ if (asp_value) {
+ asp_value = 0;
+ wat[SEG_INDEX(wat_seg, r, c)] = ABS(wat_value);
+ }
+ else if (r == 0)
asp_value = -2;
else if (c == 0)
asp_value = -4;
@@ -287,15 +300,13 @@
asp_value = -6;
else if (c == ncols - 1)
asp_value = -8;
- else
- asp_value = -1;
asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
}
else if (FLAG_GET(worked, r - 1, c)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -2;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -303,7 +314,7 @@
}
else if (FLAG_GET(worked, r + 1, c)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -6;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -311,7 +322,7 @@
}
else if (FLAG_GET(worked, r, c - 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -4;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -319,7 +330,7 @@
}
else if (FLAG_GET(worked, r, c + 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -8;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -327,7 +338,7 @@
}
else if (sides == 8 && FLAG_GET(worked, r - 1, c - 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -3;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -335,7 +346,7 @@
}
else if (sides == 8 && FLAG_GET(worked, r - 1, c + 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -1;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -343,7 +354,7 @@
}
else if (sides == 8 && FLAG_GET(worked, r + 1, c - 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -5;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -351,7 +362,7 @@
}
else if (sides == 8 && FLAG_GET(worked, r + 1, c + 1)) {
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
asp[SEG_INDEX(asp_seg, r, c)] = -7;
wat_value = wat[SEG_INDEX(wat_seg, r, c)];
if (wat_value > 0)
@@ -365,7 +376,8 @@
for (r = 0; r < nrows; r++) {
G_percent(r, nrows, 3);
for (c = 0; c < ncols; c++) {
- s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
+ if (er_flag)
+ s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
asp_value = asp[SEG_INDEX(asp_seg, r, c)];
if (r == 0 || c == 0 || r == nrows - 1 ||
c == ncols - 1 || asp_value != 0) {
@@ -373,7 +385,12 @@
if (wat_value > 0) {
wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
}
- if (r == 0)
+ /* set depression */
+ if (asp_value) {
+ asp_value = 0;
+ wat[SEG_INDEX(wat_seg, r, c)] = ABS(wat_value);
+ }
+ else if (r == 0)
asp_value = -2;
else if (c == 0)
asp_value = -4;
@@ -381,35 +398,15 @@
asp_value = -6;
else if (c == ncols - 1)
asp_value = -8;
- else
- asp_value = -1;
asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
alt_value = alt[SEG_INDEX(alt_seg, r, c)];
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
}
}
}
}
- G_percent(r, nrows, 3); /* finish it */
+ G_percent(r, nrows, 1); /* finish it */
return 0;
}
-
-int do_legal(char *file_name)
-{
- if (G_legal_filename(file_name) == -1)
- G_fatal_error(_("<%s> is an illegal file name"), file_name);
-
- return 0;
-}
-
-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);
-}
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/main.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/main.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -5,7 +5,7 @@
* AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
* Markus Neteler <neteler itc.it>, Roberto Flor <flor itc.it>,
* Brad Douglas <rez touchofmadness.com>,
- * Hamish Bowman <hamish_b yahoo com>,
+ * Hamish Bowman <hamish_b yahoo com>,
* Markus Metz <markus.metz.giswork gmail.com>
* PURPOSE: Watershed determination
* COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
@@ -15,20 +15,64 @@
* 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;
+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;
+POINT *astar_pts;
+CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
+DCELL *wat;
+CELL *ril_buf;
+int ril_fd;
+double *s_l, *s_g, *l_s;
+CELL one, zero;
+double ril_value, 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/develbranch_6/raster/r.watershed/ram/no_stream.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/no_stream.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/no_stream.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -5,7 +5,8 @@
{
int r, rr, c, cc, uprow = 0, upcol = 0;
double slope;
- CELL downdir, max_drain, value, asp_value, hih_ele, new_ele, aspect;
+ CELL downdir, asp_value, hih_ele, new_ele, aspect;
+ DCELL dvalue, max_drain; /* flow acc is now DCELL */
SHORT updir, riteflag, leftflag, thisdir;
while (1) {
@@ -15,13 +16,13 @@
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;
}
}
}
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/sg_factor.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/sg_factor.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/sg_factor.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -9,7 +9,7 @@
CELL low_elev, hih_elev;
double height, length, S, sin_theta;
- G_message(_("SECTION 4: Length Slope determination."));
+ G_message(_("SECTION 4: RUSLE LS and/or S factor determination."));
for (r = 0; r < nrows; r++) {
G_percent(r, nrows, 3);
@@ -31,14 +31,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 +66,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/develbranch_6/raster/r.watershed/ram/slope_len.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/slope_len.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/slope_len.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/Gwater.h 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/Gwater.h 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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;
+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/develbranch_6/raster/r.watershed/seg/close_maps.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/close_maps.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/close_maps.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -1,17 +1,111 @@
#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) {
+ dseg_write_cellfile(&wat, wat_name);
+
+ /* get standard deviation */
+ fd = G_open_cell_old(wat_name, "");
+ if (fd < 0) {
+ G_fatal_error(_("unable to open flow accumulation map layer"));
+ }
+
+ sum = sum_sqr = stddev = 0.0;
+ dbuf = G_allocate_d_raster_buf();
+ for (r = 0; r < nrows; r++) {
+ G_get_d_raster_row(fd, dbuf, r);
+ for (c = 0; c < ncols; c++) {
+ dvalue = dbuf[c];
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dvalue = ABS(dvalue);
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
+ }
+ }
+
+ stddev = sqrt((sum_sqr - (sum + sum / do_points)) / (do_points - 1));
+ G_debug(1, "stddev: %f", stddev);
+
+ /* set nice color rules: yellow, green, cyan, blue, black */
+ /* start with white to get more detail? NULL cells are white by default, may be confusing */
+
+ lstddev = log(stddev);
+
+ G_read_fp_range(wat_name, this_mapset, &accRange);
+ min = max = 0;
+ G_get_fp_range_min_max(&accRange, &min, &max);
+
+ G_init_colors(&colors);
+
+ if (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 +113,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 +141,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/develbranch_6/raster/r.watershed/seg/close_maps2.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/close_maps2.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/close_maps2.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/seg/do_astar.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -5,14 +5,13 @@
#include "Gwater.h"
#include "do_astar.h"
-int sift_up(int, CELL);
-
int do_astar(void)
{
POINT point;
int doer, count;
SHORT upr, upc, r, c, ct_dir;
- CELL work_val, alt_val, alt_up, asp_up, wat_val;
+ CELL work_val, alt_val, alt_up, asp_up;
+ DCELL wat_val;
CELL in_val, drain_val;
HEAP heap_pos;
@@ -47,7 +46,8 @@
c = point.c;
bseg_put(&worked, &one, r, c);
- cseg_get(&alt, &alt_val, r, c);
+ /* cseg_get(&alt, &alt_val, r, c); */
+ alt_val = heap_pos.ele;
/* check all neighbours, breadth first search */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -61,26 +61,27 @@
bseg_get(&in_list, &in_val, upr, upc);
if (in_val == 0) {
cseg_get(&alt, &alt_up, upr, upc);
- add_pt(upr, upc, r, c, alt_up, alt_val);
+ add_pt(upr, upc, alt_up, alt_val);
/* flow direction is set here */
drain_val = drain[upr - r + 1][upc - c + 1];
cseg_put(&asp, &drain_val, upr, upc);
}
else {
/* check if neighbour has not been worked on,
- * update values for asp, wat and slp */
+ * update values for asp and wat */
bseg_get(&worked, &work_val, upr, upc);
if (!work_val) {
cseg_get(&asp, &asp_up, upr, upc);
- if (asp_up < -1) {
+ if (asp_up < 0) {
drain_val = drain[upr - r + 1][upc - c + 1];
cseg_put(&asp, &drain_val, upr, upc);
- cseg_get(&wat, &wat_val, r, c);
- if (wat_val > 0)
+ dseg_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 */
+ dseg_put(&wat, &wat_val, r, c);
+ }
+ /* cseg_get(&alt, &alt_up, upr, upc); */
+ /* replace(upr, upc, r, c); */ /* alt_up used to be */
/* slope = get_slope (upr, upc, r, c, alt_up, alt_val);
dseg_put (&slp, &slope, upr, upc); */
}
@@ -89,7 +90,10 @@
}
}
}
- bseg_close(&worked);
+
+ if (mfd == 0)
+ bseg_close(&worked);
+
bseg_close(&in_list);
seg_close(&heap_index);
@@ -98,7 +102,7 @@
}
/* new add point routine for min heap */
-int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
{
POINT point;
HEAP heap_pos;
@@ -122,11 +126,13 @@
point.r = r;
point.c = c;
- point.downr = downr;
- point.downc = downc;
+/* point.downr = downr;
+ point.downc = downc; */
seg_put(&astar_pts, (char *)&point, 0, nxt_avail_pt);
+ /* cseg_put(&pnt_index, &nxt_avail_pt, r, c); */
+
nxt_avail_pt++;
/* sift up: move new point towards top of heap */
@@ -140,8 +146,8 @@
int drop_pt(void)
{
int child, childr, parent;
- int childp, childrp;
- CELL ele, eler;
+ int childp; /* , childrp */
+ CELL ele; /* , eler */
int i;
HEAP heap_pos;
@@ -168,23 +174,23 @@
ele = heap_pos.ele;
if (child < heap_size) {
childr = child + 1;
- i = 1;
- while (childr <= heap_size && i < 3) {
+ i = child + 3;
+ while (childr <= heap_size && childr < i) {
seg_get(&heap_index, (char *)&heap_pos, 0, childr);
- childrp = heap_pos.point;
- eler = heap_pos.ele;
- if (eler < ele) {
+ /* childrp = heap_pos.point;
+ eler = heap_pos.ele; */
+ if (heap_pos.ele < ele) {
child = childr;
- childp = childrp;
- ele = eler;
+ childp = heap_pos.point;
+ ele = heap_pos.ele;
}
/* make sure we get the oldest child */
- else if (ele == eler && childp > childrp) {
+ else if (ele == heap_pos.ele && childp > heap_pos.point) {
child = childr;
- childp = childrp;
+ childp = heap_pos.point;
}
childr++;
- i++;
+ /* i++; */
}
}
@@ -218,8 +224,8 @@
/* standard sift-up routine for d-ary min heap */
int sift_up(int start, CELL ele)
{
- int parent, parentp, child, childp;
- CELL elep;
+ int parent, child, childp; /* parentp, */
+ /* CELL elep; */
HEAP heap_pos;
child = start;
@@ -229,11 +235,11 @@
while (child > 1) {
parent = GET_PARENT(child);
seg_get(&heap_index, (char *)&heap_pos, 0, parent);
- parentp = heap_pos.point;
- elep = heap_pos.ele;
+ /* parentp = heap_pos.point;
+ elep = heap_pos.ele; */
/* parent ele higher */
- if (elep > ele) {
+ if (heap_pos.ele > ele) {
/* push parent point down */
seg_put(&heap_index, (char *)&heap_pos, 0, child);
@@ -241,7 +247,7 @@
}
/* same ele, but parent is younger */
- else if (elep == ele && parentp > childp) {
+ else if (heap_pos.ele == ele && heap_pos.point > childp) {
/* push parent point down */
seg_put(&heap_index, (char *)&heap_pos, 0, child);
@@ -278,27 +284,24 @@
return (slope);
}
-int replace( /* ele was in there */
- SHORT upr, SHORT upc, SHORT r, SHORT c)
-/* CELL ele; */
-{
- int now, heap_run;
+/* 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/develbranch_6/raster/r.watershed/seg/do_astar.h
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.h 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.h 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/do_cum.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/do_cum.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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,336 @@
}
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;
+ 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;
+ /* 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) {
+
+ /* check for swale or stream cells */
+ bseg_get(&swale, &is_swale, r_nbr, c_nbr);
+ if (is_swale)
+ swale_cells++;
+ dseg_get(&wat, &valued, r_nbr, c_nbr);
+ wat_nbr[ct_dir] = valued;
+ if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold)
+ stream_cells++;
+
+ bseg_get(&worked, &cvalue, r_nbr, c_nbr);
+ worked_nbr[ct_dir] = cvalue;
+ if (worked_nbr[ct_dir] == 0) {
+ cseg_get(&alt, &ele_nbr, r_nbr, c_nbr);
+ is_null = G_is_c_null_value(&ele_nbr);
+ edge = is_null;
+ 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) {
+ 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/develbranch_6/raster/r.watershed/seg/dseg_get.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/dseg_get.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/dseg_get.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/seg/dseg_put.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/dseg_put.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/dseg_put.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/seg/dseg_read.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/dseg_read.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/dseg_read.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/seg/dseg_write.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/dseg_write.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/dseg_write.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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/develbranch_6/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/init_vars.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/init_vars.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -1,4 +1,5 @@
#include <stdlib.h>
+#include <string.h>
#include <unistd.h>
#include "Gwater.h"
#include <grass/gis.h>
@@ -14,9 +15,9 @@
/* 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;
+ char MASK_flag;
G_gisinit(argv[0]);
ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
@@ -24,11 +25,14 @@
zero = sl_flag = sg_flag = ls_flag = er_flag = bas_thres = 0;
nxt_avail_pt = 0;
/* dep_flag = 0; */
- max_length = dzero = 0.0;
+ max_length = d_zero = 0.0;
+ d_one = 1.0;
ril_value = -1.0;
/* dep_slope = 0.0; */
max_bytes = 0;
sides = 8;
+ mfd = 1;
+ c_fac = 5;
segs_mb = 300;
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
@@ -73,9 +77,15 @@
if (sides != 4)
usage(argv[0]);
}
+ else if (sscanf(argv[r], "conv=%d", &c_fac) == 1) ;
+ else if (strcmp(argv[r], "-s") == 0)
+ mfd = 0;
else
usage(argv[0]);
}
+ if (mfd == 1 && (c_fac < 1 || c_fac > 10)) {
+ G_fatal_error("Convergence factor must be between 1 and 10.");
+ }
if ((ele_flag != 1)
||
((arm_flag == 1) &&
@@ -96,27 +106,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 +118,7 @@
G_get_set_window(&window);
nrows = G_window_rows();
ncols = G_window_cols();
- if (max_length <= dzero)
+ if (max_length <= d_zero)
max_length = 10 * nrows * window.ns_res + 10 * ncols * window.ew_res;
if (window.ew_res < window.ns_res)
half_res = .5 * window.ew_res;
@@ -137,8 +131,7 @@
/* segment parameters: one size fits all. Fine tune? */
/* Segment rows and cols: 200 */
- /* 1 segment open for all rasters: 2.86 MB */
- /* num_open_segs = segs_mb / 2.86 */
+ /* 1 segment open for all rasters: 1.34 MB */
seg_rows = SROW;
seg_cols = SCOL;
@@ -148,7 +141,7 @@
G_warning(_("Maximum memory to be used was smaller than 3 MB, set to default = 300 MB."));
}
- num_open_segs = segs_mb / 2.86;
+ num_open_segs = segs_mb / 1.34;
G_debug(1, "segs MB: %.0f", segs_mb);
G_debug(1, "region rows: %d", nrows);
@@ -172,26 +165,92 @@
G_debug(1, " open segments after adjusting:\t%d", num_open_segs);
cseg_open(&alt, seg_rows, seg_cols, num_open_segs);
- cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
- cseg_read_cell(&alt, ele_name, 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, "");
+ }
+
+ /* NULL cells in input elevation map */
+ bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
+ bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
+ G_debug(1, "Checking for masked and NULL cells in input elevation <%s>.", ele_name);
+ MASK_flag = 0;
+ do_points = nrows * ncols;
+ fd = G_open_cell_old(ele_name, "");
+ if (fd < 0) {
+ G_fatal_error(_("unable to open elevation map layer"));
+ }
+ buf = G_allocate_cell_buf();
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ /* check for masked and NULL cells */
+ alt_value = buf[c];
+ if (G_is_c_null_value(&alt_value)) {
+ bseg_put(&worked, &one, r, c);
+ bseg_put(&in_list, &one, r, c);
+ do_points--;
+ }
+ }
+ }
+ G_close_cell(fd);
+ G_free(buf);
+ if (do_points < nrows * ncols)
+ MASK_flag = 1;
+
+ /* initial flow accumulation */
+ dseg_open(&wat, seg_rows, seg_cols, num_open_segs);
if (run_flag) {
- 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 +261,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,48 +288,29 @@
}
}
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 */
+ seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols * 2,
+ num_open_segs / 2, sizeof(POINT));
+
+ /* heap_index will track astar_pts in ternary min-heap */
/* heap_index is one-based */
- seg_open(&heap_index, 1, do_points + 1, 1, seg_rows * seg_cols,
- num_open_segs, sizeof(HEAP));
+ seg_open(&heap_index, 1, do_points + 1, 1, seg_cols * num_open_segs * seg_rows / 10,
+ 10, sizeof(HEAP));
G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
@@ -264,23 +321,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 +354,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 +469,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 +499,10 @@
asp_value = -6;
else if (c == ncols - 1)
asp_value = -8;
- else
- asp_value = -1;
if (-1 == cseg_put(&asp, &asp_value, r, c))
exit(EXIT_FAILURE);
cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, -1, -1, alt_value, alt_value);
+ add_pt(r, c, alt_value, alt_value);
}
else {
bseg_put(&in_list, &zero, r, c);
@@ -440,25 +511,7 @@
}
}
}
- G_percent(r, nrows, 3); /* finish it */
+ G_percent(r, nrows, 1); /* finish it */
return 0;
}
-
-int do_legal(char *file_name)
-{
- if (G_legal_filename(file_name) == -1)
- G_fatal_error(_("<%s> is an illegal file name"), file_name);
-
- return 0;
-}
-
-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);
-}
Modified: grass/branches/develbranch_6/raster/r.watershed/seg/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/main.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/main.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -16,25 +16,66 @@
* 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;
+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/develbranch_6/raster/r.watershed/seg/no_stream.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/no_stream.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/no_stream.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -6,7 +6,8 @@
{
int r, rr, c, cc, uprow = 0, upcol = 0;
double slope;
- CELL downdir, max_drain, value, asp_value, hih_ele, new_ele, aspect;
+ CELL downdir, asp_value, hih_ele, new_ele, aspect, value;
+ DCELL dvalue, max_drain; /* flow acc is now DCELL */
SHORT updir, riteflag, leftflag, thisdir;
while (1) {
@@ -17,13 +18,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;
}
}
}
Modified: grass/branches/develbranch_6/raster/r.watershed/seg/sg_factor.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/sg_factor.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/sg_factor.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -9,10 +9,9 @@
CELL downer, low_elev, hih_elev;
double height, length, S, sin_theta;
- G_message(_("SECTION 4: Length Slope determination."));
+ G_message(_("SECTION 4: RUSLE LS and/or S factor determination."));
for (r = nrows - 1; r >= 0; r--) {
-/* FIXME: G_percent params count backwards */
- G_percent(r, nrows, 3);
+ G_percent(nrows - r, nrows, 3);
for (c = ncols - 1; c >= 0; c--) {
cseg_get(&alt, &low_elev, r, c);
cseg_get(&r_h, &hih_elev, r, c);
@@ -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/develbranch_6/raster/r.watershed/seg/slope_len.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/slope_len.c 2009-02-16 14:56:00 UTC (rev 35904)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/slope_len.c 2009-02-17 12:16:36 UTC (rev 35905)
@@ -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