[GRASS-SVN] r36273 - in grass/trunk/raster/r.watershed: front ram
seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Mar 9 04:59:34 EDT 2009
Author: mmetz
Date: 2009-03-09 04:59:34 -0400 (Mon, 09 Mar 2009)
New Revision: 36273
Modified:
grass/trunk/raster/r.watershed/front/main.c
grass/trunk/raster/r.watershed/front/r.watershed.html
grass/trunk/raster/r.watershed/ram/Gwater.h
grass/trunk/raster/r.watershed/ram/close_maps.c
grass/trunk/raster/r.watershed/ram/do_astar.c
grass/trunk/raster/r.watershed/ram/init_vars.c
grass/trunk/raster/r.watershed/ram/main.c
grass/trunk/raster/r.watershed/ram/sg_factor.c
grass/trunk/raster/r.watershed/seg/Gwater.h
grass/trunk/raster/r.watershed/seg/close_maps.c
grass/trunk/raster/r.watershed/seg/init_vars.c
grass/trunk/raster/r.watershed/seg/main.c
grass/trunk/raster/r.watershed/seg/sg_factor.c
Log:
added flag for absolute flow accumulation, some minor bugs fixed
Modified: grass/trunk/raster/r.watershed/front/main.c
===================================================================
--- grass/trunk/raster/r.watershed/front/main.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/front/main.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -6,8 +6,8 @@
* Brad Douglas <rez touchofmadness.com>,
* Hamish Bowman <hamish_b yahoo.com>
* Markus Metz <markus.metz.giswork gmail.com>
- * PURPOSE: Watershed determination
- * COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
+ * PURPOSE: Hydrological analysis
+ * COPYRIGHT: (C) 1999-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -38,7 +38,6 @@
struct Option *opt10;
struct Option *opt11;
struct Option *opt12;
- struct Option *opt13;
struct Option *opt14;
struct Option *opt15;
struct Option *opt16;
@@ -46,6 +45,7 @@
struct Flag *flag_sfd;
struct Flag *flag_flow;
struct Flag *flag_seg;
+ struct Flag *flag_abs;
struct GModule *module;
G_gisinit(argv[0]);
@@ -103,16 +103,14 @@
_("Input value: minimum size of exterior watershed basin");
opt6->required = NO;
opt6->type = TYPE_INTEGER;
- opt6->gisprompt = "new,cell,raster";
opt6->guisection = _("Input_options");
opt7 = G_define_option();
opt7->key = "max_slope_length";
opt7->description =
- _("Input value: maximum length of surface flow, for USLE");
+ _("Input value: maximum length of surface flow in map units, for USLE");
opt7->required = NO;
opt7->type = TYPE_DOUBLE;
- opt7->gisprompt = "new,cell,raster";
opt7->guisection = _("Input_options");
opt8 = G_define_option();
@@ -158,6 +156,7 @@
opt12->gisprompt = "new,cell,raster";
opt12->guisection = _("Output_options");
+/* to be removed completely. visual display of results is now done with accumulation output
opt13 = G_define_option();
opt13->key = "visual";
opt13->description =
@@ -166,7 +165,7 @@
opt13->type = TYPE_STRING;
opt13->gisprompt = "new,cell,raster";
opt13->guisection = _("Output_options");
-
+*/
opt14 = G_define_option();
opt14->key = "length_slope";
opt14->description =
@@ -213,9 +212,18 @@
flag_seg = G_define_flag();
flag_seg->key = 'm';
+ flag_seg->label =
+ _("Enable disk swap memory option: Operation is slow");
flag_seg->description =
- _("Enable disk swap memory option: Operation is slow");
+ _("Only needed if memory requirements exceed available RAM; see manual on how to calculate memory requirements");
+ flag_abs = G_define_flag();
+ flag_abs->key = 'a';
+ flag_abs->label =
+ _("Use positive flow accumulation even for likely underestimates");
+ flag_abs->description =
+ _("See manual for a detailed description of flow accumulation output");
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -228,7 +236,6 @@
&& (opt10->answer == NULL)
&& (opt11->answer == NULL)
&& (opt12->answer == NULL)
- && (opt13->answer == NULL)
&& (opt14->answer == NULL)
&& (opt15->answer == NULL)) {
G_fatal_error(_("Sorry, you must choose an output map."));
@@ -270,6 +277,9 @@
if (flag_flow->answer)
strcat(command, " -4");
+ if (flag_abs->answer)
+ strcat(command, " -a");
+
if (opt1->answer) {
strcat(command, " el=");
strcat(command, "\"");
@@ -350,12 +360,13 @@
strcat(command, "\"");
}
- if (opt13->answer) {
+/* if (opt13->answer) {
strcat(command, " di=");
strcat(command, "\"");
strcat(command, opt13->answer);
strcat(command, "\"");
}
+*/
if (opt14->answer) {
strcat(command, " LS=");
@@ -410,10 +421,10 @@
write_hist(opt12->answer,
"Watershed half-basins", opt1->answer, flag_seg->answer,
flag_sfd->answer);
- if (opt13->answer)
+ /* if (opt13->answer)
write_hist(opt13->answer,
"Watershed visualization map (filtered accumulation map)",
- opt1->answer, flag_seg->answer, flag_sfd->answer);
+ opt1->answer, flag_seg->answer, flag_sfd->answer); */
if (opt14->answer)
write_hist(opt14->answer,
"Watershed slope length and steepness (LS) factor",
Modified: grass/trunk/raster/r.watershed/front/r.watershed.html
===================================================================
--- grass/trunk/raster/r.watershed/front/r.watershed.html 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/front/r.watershed.html 2009-03-09 08:59:34 UTC (rev 36273)
@@ -40,6 +40,13 @@
surface flow (allows horizontal, vertical, and diagonal flow of water).
This flag will also make the drainage basins look more homogeneous.
+<dt><em>-a</em>
+
+<dd>Use positive flow accumulation even for likely underestimates. When this
+flag is not set, cells with a flow accumulation value that is likely to be
+an underestimate are converted to the negative. See below for a detailed
+description of flow accumulation output.
+
<dt><em>memory</em>
<dd>Maximum amount of memory in MB to be used with -m set. More memory
@@ -151,14 +158,6 @@
cells of the watershed basin are given odd values which are one less
than the value of the watershed basin.
-<dt><em>visual</em>
-
-<dd>Output map: useful for visual display of results.
-Surface runoff accumulation with the values
-modified to provide for easy display. All negative accumulation values
-are changed to zero. All positive values above the basin threshold
-are given the value of the <em>threshold</em> parameter.
-
<dt><em>length_slope</em>
<dd>Output map: slope length and steepness (LS) factor for the Revised
@@ -212,14 +211,13 @@
MFD, water flow is distributed to all neighbouring cells with lower
elevation, using slope towards neighbouring cells as a weighing factor
for proportional distribution. The A<sup>T</sup> least-cost path is
-always included. As a result, depressions and obstacles are overflown
+always included. As a result, depressions and obstacles are traversed
with a gracefull flow convergence before the overflow. The convergence
factor causes flow accumulation to converge more strongly with higher
values. The supported range is 1 to 10, recommended is a convergence
factor of 5 (Holmgren, 1994). If many small sliver basins are created
with MFD, setting the convergence factor to a higher value can reduce
the amount of small sliver basins.
-<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>.
@@ -264,17 +262,6 @@
using the values from the <em>neighborhood</em> output map layer that
represents the minimum elevation within the region of the coarser cell.
-<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
@@ -306,6 +293,16 @@
<h4>Further processing of output layers</h4>
<p>
+Problem areas, i.e. those parts of a basin with a likely underestimate of
+flow accumulation, can be easily identified with e.g.
+<div class="code"><pre>
+ r.mapcalc "problems = if(flow_acc < 0, basin, null())"
+</pre></div>
+If the region of interest contains such problem areas, and this is not
+desired, the computational region must be expanded until the catchment
+area for the region of interest is completely included.
+
+<p>
To isolate an individual river network using the output of this module,
a number of approaches may be considered.
<ol>
@@ -318,6 +315,11 @@
starting point.
</ol>
+All individual river networks in the stream segments output can be
+identified through their ultimate outlet points. These points are all
+cells in the stream segments output with negative drainage direction.
+These points can be used as start points for <em>r.water.outlet</em> or
+<em>v.net.iso</em>.
<p>
To create <i>river mile</i> segmentation from a vectorized streams map,
@@ -349,7 +351,7 @@
<br>
<p>
-Set a nice color table for the accumulation map:
+Set a different color table for the accumulation map:
<div class="code"><pre>
MAP=rwater.accum
r.watershed elev=elevation.dem accum=$MAP
@@ -380,7 +382,7 @@
it to a vector output map. The accumulation cut-off, and therefore fractal
dimension, is arbitrary; in this example we use the map's mean number of
upstream catchment cells (calculated in the above example by <em>r.univar</em>)
-as the cut-off value.
+as the cut-off value. This only works with SFD, not with MFD.
<div class="code"><pre>
r.watershed elev=elevation.dem accum=rwater.accum
@@ -424,18 +426,6 @@
</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>
@@ -496,7 +486,7 @@
Original version:
Charles Ehlschlaeger, U.S. Army Construction Engineering Research Laboratory
<br>
-Faster sorting algorithm:
+Faster sorting algorithm and MFD support:
Markus Metz <markus.metz.giswork at gmail.com>
<p>
Modified: grass/trunk/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/Gwater.h 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/ram/Gwater.h 2009-03-09 08:59:34 UTC (rev 36273)
@@ -40,7 +40,7 @@
extern struct Cell_head window;
-extern int mfd, c_fac;
+extern int mfd, c_fac, abs_acc, ele_scale;
extern int *heap_index, heap_size;
extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
extern SHORT nrows, ncols;
@@ -53,7 +53,6 @@
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;
Modified: grass/trunk/raster/r.watershed/ram/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/close_maps.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/ram/close_maps.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -30,19 +30,38 @@
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++) {
- 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;
+ if (abs_acc) {
+ G_warning("Writing out only positive flow accumulation values.");
+ G_warning("Cells with a likely underestimate for flow accumulation can no longer be identified.");
+ for (r = 0; r < nrows; r++) {
+ G_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
+ for (c = 0; c < ncols; c++) {
+ dvalue = wat[SEG_INDEX(wat_seg, r, c)];
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dvalue = ABS(dvalue);
+ dbuf[c] = dvalue;
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
}
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
- G_put_raster_row(fd, dbuf, DCELL_TYPE);
}
+ else {
+ for (r = 0; r < nrows; r++) {
+ G_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
+ for (c = 0; c < ncols; c++) {
+ dvalue = wat[SEG_INDEX(wat_seg, r, c)];
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dbuf[c] = dvalue;
+ dvalue = ABS(dvalue);
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
+ }
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
+ }
+ }
if (G_close_cell(fd) < 0)
G_warning(_("Close failed."));
Modified: grass/trunk/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/ram/do_astar.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -92,6 +92,8 @@
}
}
}
+ /* this was a lot of indentation, all in the sake of speed... */
+ /* improve code aesthetics? */
G_percent(count, do_points, 3); /* finish it */
if (mfd == 0)
flag_destroy(worked);
@@ -134,9 +136,9 @@
/* new drop point routine for min heap */
int drop_pt(void)
{
- int child, childr, parent;
+ register int child, childr, parent;
CELL ele, eler;
- int i;
+ register int i;
if (heap_size == 1) {
heap_index[1] = -1;
@@ -152,17 +154,15 @@
while ((child = GET_CHILD(parent)) <= heap_size) {
/* select child with lower ele, if equal, older child
- * older child is older startpoint for flow path, important
- * chances are good that the left child will be the older child,
- * but just to make sure we test */
+ * older child is older startpoint for flow path, important */
ele =
alt[SEG_INDEX
(alt_seg, astar_pts[heap_index[child]].r,
astar_pts[heap_index[child]].c)];
if (child < heap_size) {
childr = child + 1;
- i = 1;
- while (childr <= heap_size && i < 3) {
+ i = child + 3; /* change the number, GET_CHILD() and GET_PARENT() to play with different d-ary heaps */
+ while (childr <= heap_size && childr < i) {
eler =
alt[SEG_INDEX
(alt_seg, astar_pts[heap_index[childr]].r,
@@ -170,15 +170,23 @@
if (eler < ele) {
child = childr;
ele = eler;
- /* make sure we get the oldest child */
}
+ /* make sure we get the oldest child */
else if (ele == eler &&
heap_index[child] > heap_index[childr]) {
child = childr;
}
childr++;
- i++;
+ /* i++; */
}
+ /* break if childr > last entry? that saves sifting up again
+ * OTOH, this is another comparison
+ * we have a max heap height of 20: log(INT_MAX)/log(n children per node)
+ * that would give us in the worst case 20*2 additional comparisons with 3 children
+ * the last entry will never go far up again, less than half the way
+ * so the additional comparisons for going all the way down
+ * and then a bit up again are likely less than 20*2 */
+ /* find the error in this reasoning */
}
/* move hole down */
@@ -210,7 +218,7 @@
/* standard sift-up routine for d-ary min heap */
int sift_up(int start, CELL ele)
{
- int parent, parentp, child, childp;
+ register int parent, parentp, child, childp;
CELL elep;
child = start;
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -4,12 +4,16 @@
#include <grass/gis.h>
#include <grass/glocale.h>
+int ele_round(double);
int init_vars(int argc, char *argv[])
{
SHORT r, c;
CELL *buf, alt_value, wat_value, asp_value, block_value;
- int fd, index;
+ DCELL dvalue;
+ void *elebuf, *ptr;
+ int fd, index, ele_map_type;
+ size_t ele_size;
char MASK_flag;
G_gisinit(argv[0]);
@@ -26,6 +30,8 @@
sides = 8;
mfd = 1;
c_fac = 5;
+ abs_acc = 0;
+ ele_scale = 1;
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
ele_flag++;
@@ -71,6 +77,8 @@
else if (sscanf(argv[r], "conv=%d", &c_fac) == 1) ;
else if (strcmp(argv[r], "-s") == 0)
mfd = 0;
+ else if (strcmp(argv[r], "-a") == 0)
+ abs_acc = 1;
else
usage(argv[0]);
}
@@ -118,7 +126,7 @@
window.ns_res * window.ns_res);
if (sides == 4)
diag *= 0.5;
- buf = G_allocate_cell_buf();
+
alt =
(CELL *) G_malloc(sizeof(CELL) * size_array(&alt_seg, nrows, ncols));
if (er_flag) {
@@ -126,45 +134,79 @@
(CELL *) G_malloc(sizeof(CELL) * size_array(&r_h_seg, nrows, ncols));
}
+ swale = flag_create(nrows, ncols);
+ in_list = flag_create(nrows, ncols);
+ worked = flag_create(nrows, ncols);
+
+ /* open elevation input */
fd = G_open_cell_old(ele_name, "");
if (fd < 0) {
G_fatal_error(_("unable to open elevation map layer"));
}
- swale = flag_create(nrows, ncols);
- in_list = flag_create(nrows, ncols);
- worked = flag_create(nrows, ncols);
+ ele_map_type = G_get_raster_map_type(fd);
+ ele_size = G_raster_size(ele_map_type);
+ elebuf = G_allocate_raster_buf(ele_map_type);
+ if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
+ ele_scale = 1000; /* should be enough to do the trick */
+
+ /* read elevation input and mark NULL/masked cells */
MASK_flag = 0;
do_points = nrows * ncols;
for (r = 0; r < nrows; r++) {
- G_get_c_raster_row(fd, buf, r);
+ G_get_raster_row(fd, elebuf, r, ele_map_type);
+ ptr = elebuf;
for (c = 0; c < ncols; c++) {
index = SEG_INDEX(alt_seg, r, c);
- alt_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);
- /* check for masked and NULL cells here */
- if (G_is_c_null_value(&alt_value)) {
+
+ /* check for masked and NULL cells */
+ if (G_is_null_value(ptr, ele_map_type)) {
FLAG_SET(worked, r, c);
FLAG_SET(in_list, r, c);
+ G_set_c_null_value(&alt_value, 1);
do_points--;
}
+ else {
+ if (ele_map_type == CELL_TYPE) {
+ alt_value = *((CELL *)ptr);
+ }
+ else if (ele_map_type == FCELL_TYPE) {
+ dvalue = *((FCELL *)ptr);
+ dvalue *= ele_scale;
+ alt_value = ele_round(dvalue);
+ }
+ else if (ele_map_type == DCELL_TYPE) {
+ dvalue = *((DCELL *)ptr);
+ dvalue *= ele_scale;
+ alt_value = ele_round(dvalue);
+ }
+ }
+ alt[index] = alt_value;
+ if (er_flag) {
+ r_h[index] = alt_value;
+ }
+ ptr = G_incr_void_ptr(ptr, ele_size);
}
}
G_close_cell(fd);
+ G_free(elebuf);
if (do_points < nrows * ncols)
MASK_flag = 1;
+
+ /* initialize flow accumulation ... */
wat =
(DCELL *) G_malloc(sizeof(DCELL) *
size_array(&wat_seg, nrows, ncols));
+ buf = G_allocate_cell_buf();
if (run_flag) {
+ /* ... with input map flow: amount of overland flow per cell */
fd = G_open_cell_old(run_name, "");
if (fd < 0) {
G_fatal_error(_("unable to open runoff map layer"));
@@ -186,6 +228,7 @@
G_close_cell(fd);
}
else {
+ /* ... with 1.0 */
for (r = 0; r < nrows; r++) {
for (c = 0; c < ncols; c++) {
if (MASK_flag) {
@@ -234,6 +277,8 @@
}
G_close_cell(fd);
}
+ G_free(buf);
+
if (ril_flag) {
ril_fd = G_open_cell_old(ril_name, "");
if (ril_fd < 0) {
@@ -410,3 +455,17 @@
return 0;
}
+
+int ele_round(double x)
+{
+ int n;
+
+ if (x >= 0.0)
+ n = x + .5;
+ else {
+ n = -x + .5;
+ n = -n;
+ }
+
+ return n;
+}
Modified: grass/trunk/raster/r.watershed/ram/main.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/main.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/ram/main.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -7,8 +7,8 @@
* Brad Douglas <rez touchofmadness.com>,
* Hamish Bowman <hamish_b yahoo com>,
* Markus Metz <markus.metz.giswork gmail.com>
- * PURPOSE: Watershed determination
- * COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
+ * PURPOSE: Hydrological analysis
+ * COPYRIGHT: (C) 1999-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -24,7 +24,7 @@
struct Cell_head window;
-int mfd, c_fac;
+int mfd, c_fac, abs_acc, ele_scale;
int *heap_index, heap_size;
int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
SHORT nrows, ncols;
@@ -37,7 +37,6 @@
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;
Modified: grass/trunk/raster/r.watershed/ram/sg_factor.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/sg_factor.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/ram/sg_factor.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -2,6 +2,7 @@
#include <grass/gis.h>
#include <grass/glocale.h>
+CELL *ril_buf;
int sg_factor(void)
{
@@ -11,6 +12,9 @@
G_message(_("SECTION 4: RUSLE LS and/or S factor determination."));
+ if (ril_flag)
+ ril_buf = G_allocate_cell_buf();
+
for (r = 0; r < nrows; r++) {
G_percent(r, nrows, 3);
if (ril_flag) {
@@ -20,7 +24,7 @@
low_elev = alt[SEG_INDEX(alt_seg, r, c)];
hih_elev = r_h[SEG_INDEX(r_h_seg, r, c)];
length = s_l[SEG_INDEX(s_l_seg, r, c)];
- height = hih_elev - low_elev;
+ height = 1.0 * (hih_elev - low_elev) / ele_scale;
if (length > max_length) {
height *= max_length / length;
length = max_length;
@@ -66,7 +70,7 @@
/* rill_ratio equation from Steve Warren */
rill_ratio *= .5 + .005 * ril + .0001 * ril * ril;
s_l_exp = rill_ratio / (1 + rill_ratio);
- L = 100 * pow((slope_length / 72.6), s_l_exp);
+ L = pow((slope_length / 72.6), s_l_exp);
l_s[SEG_INDEX(l_s_seg, r, c)] = L * S;
return 0;
Modified: grass/trunk/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/Gwater.h 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/seg/Gwater.h 2009-03-09 08:59:34 UTC (rev 36273)
@@ -43,7 +43,7 @@
extern struct Cell_head window;
-extern int mfd, c_fac;
+extern int mfd, c_fac, abs_acc, ele_scale;
extern SSEG heap_index;
extern int heap_size;
extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
Modified: grass/trunk/raster/r.watershed/seg/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/close_maps.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/seg/close_maps.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -17,24 +17,50 @@
dseg_close(&slp);
cseg_close(&alt);
if (wat_flag) {
- dseg_write_cellfile(&wat, wat_name);
+ sum = sum_sqr = stddev = 0.0;
+ dbuf = G_allocate_d_raster_buf();
+ if (abs_acc) {
+ G_warning("Writing out only positive flow accumulation values.");
+ G_warning("Cells with a likely underestimate for flow accumulation can no longer be identified.");
- /* get standard deviation */
- fd = G_open_cell_old(wat_name, "");
- if (fd < 0) {
- G_fatal_error(_("unable to open flow accumulation map layer"));
+ fd = G_open_raster_new(wat_name, DCELL_TYPE);
+ if (fd < 0) {
+ G_warning(_("unable to open new accum map layer."));
+ }
+ for (r = 0; r < nrows; r++) {
+ G_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
+ for (c = 0; c < ncols; c++) {
+ dseg_get(&wat, &dvalue, r, c);
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dvalue = ABS(dvalue);
+ dbuf[c] = dvalue;
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
+ }
+ G_put_raster_row(fd, dbuf, DCELL_TYPE);
+ }
+ if (G_close_cell(fd) < 0)
+ G_warning(_("Close failed."));
}
+ else {
+ dseg_write_cellfile(&wat, wat_name);
- 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;
+ /* get standard deviation */
+ fd = G_open_cell_old(wat_name, "");
+ if (fd < 0) {
+ G_fatal_error(_("unable to open flow accumulation map layer"));
+ }
+
+ for (r = 0; r < nrows; r++) {
+ G_get_d_raster_row(fd, dbuf, r);
+ for (c = 0; c < ncols; c++) {
+ dvalue = dbuf[c];
+ if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+ dvalue = ABS(dvalue);
+ sum += dvalue;
+ sum_sqr += dvalue * dvalue;
+ }
}
}
}
Modified: grass/trunk/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/init_vars.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/seg/init_vars.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -5,6 +5,7 @@
#include <grass/gis.h>
#include <grass/glocale.h>
+int ele_round(double);
int init_vars(int argc, char *argv[])
{
@@ -17,7 +18,11 @@
int max_bytes;
CELL *buf, alt_value, asp_value, worked_value, block_value;
DCELL wat_value;
+ DCELL dvalue;
char MASK_flag;
+ void *elebuf, *ptr;
+ int ele_map_type;
+ size_t ele_size;
G_gisinit(argv[0]);
ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
@@ -33,6 +38,8 @@
sides = 8;
mfd = 1;
c_fac = 5;
+ abs_acc = 0;
+ ele_scale = 1;
segs_mb = 300;
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
@@ -80,6 +87,8 @@
else if (sscanf(argv[r], "conv=%d", &c_fac) == 1) ;
else if (strcmp(argv[r], "-s") == 0)
mfd = 0;
+ else if (strcmp(argv[r], "-a") == 0)
+ abs_acc = 1;
else
usage(argv[0]);
}
@@ -171,31 +180,63 @@
cseg_read_cell(&r_h, ele_name, "");
}
- /* NULL cells in input elevation map */
+ /* read elevation input and mark NULL/masked cells */
bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
- G_debug(1, "Checking for masked and NULL cells in input elevation <%s>.", ele_name);
- MASK_flag = 0;
- do_points = nrows * ncols;
+ G_verbose_message("Checking for masked and NULL cells in input elevation <%s>", ele_name);
+
+ /* open elevation input */
fd = G_open_cell_old(ele_name, "");
if (fd < 0) {
G_fatal_error(_("unable to open elevation map layer"));
}
- buf = G_allocate_cell_buf();
+
+ ele_map_type = G_get_raster_map_type(fd);
+ ele_size = G_raster_size(ele_map_type);
+ elebuf = G_allocate_raster_buf(ele_map_type);
+
+ if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
+ ele_scale = 1000; /* should be enough to do the trick */
+
+ /* read elevation input and mark NULL/masked cells */
+ MASK_flag = 0;
+ do_points = nrows * ncols;
for (r = 0; r < nrows; r++) {
- G_get_c_raster_row(fd, buf, r);
+ G_get_raster_row(fd, elebuf, r, ele_map_type);
+ ptr = elebuf;
for (c = 0; c < ncols; c++) {
+
/* check for masked and NULL cells */
- alt_value = buf[c];
- if (G_is_c_null_value(&alt_value)) {
+ if (G_is_null_value(ptr, ele_map_type)) {
bseg_put(&worked, &one, r, c);
bseg_put(&in_list, &one, r, c);
+ G_set_c_null_value(&alt_value, 1);
do_points--;
}
+ else {
+ if (ele_map_type == CELL_TYPE) {
+ alt_value = *((CELL *)ptr);
+ }
+ else if (ele_map_type == FCELL_TYPE) {
+ dvalue = *((FCELL *)ptr);
+ dvalue *= ele_scale;
+ alt_value = ele_round(dvalue);
+ }
+ else if (ele_map_type == DCELL_TYPE) {
+ dvalue = *((DCELL *)ptr);
+ dvalue *= ele_scale;
+ alt_value = ele_round(dvalue);
+ }
+ }
+ cseg_put(&alt, &alt_value, r, c);
+ if (er_flag) {
+ cseg_put(&r_h, &alt_value, r, c);
+ }
+ ptr = G_incr_void_ptr(ptr, ele_size);
}
}
G_close_cell(fd);
- G_free(buf);
+ G_free(elebuf);
if (do_points < nrows * ncols)
MASK_flag = 1;
@@ -515,3 +556,17 @@
return 0;
}
+
+int ele_round(double x)
+{
+ int n;
+
+ if (x >= 0.0)
+ n = x + .5;
+ else {
+ n = -x + .5;
+ n = -n;
+ }
+
+ return n;
+}
Modified: grass/trunk/raster/r.watershed/seg/main.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/main.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/seg/main.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -8,8 +8,8 @@
* Brad Douglas <rez touchofmadness.com>,
* Hamish Bowman <hamish_b yahoo.com>,
* Markus Metz <markus.metz.giswork gmail.com>
- * PURPOSE: Watershed determination using the GRASS segmentation lib
- * COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
+ * PURPOSE: Hydrological analysis using the GRASS segmentation lib
+ * COPYRIGHT: (C) 1999-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -25,7 +25,7 @@
struct Cell_head window;
-int mfd, c_fac;
+int mfd, c_fac, abs_acc, ele_scale;
SSEG heap_index;
int heap_size;
int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
Modified: grass/trunk/raster/r.watershed/seg/sg_factor.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/sg_factor.c 2009-03-09 08:54:47 UTC (rev 36272)
+++ grass/trunk/raster/r.watershed/seg/sg_factor.c 2009-03-09 08:59:34 UTC (rev 36273)
@@ -17,7 +17,7 @@
cseg_get(&r_h, &hih_elev, r, c);
dseg_get(&s_l, &length, r, c);
cseg_get(&asp, &downer, r, c);
- height = hih_elev - low_elev;
+ height = 1.0 * (hih_elev - low_elev) / ele_scale;
if (length > max_length) {
height *= max_length / length;
length = max_length;
@@ -32,7 +32,6 @@
len_slp_equ(length, sin_theta, S, r, c);
}
if (sg_flag) {
- S *= 100.0;
dseg_put(&s_g, &S, r, c);
}
}
@@ -60,7 +59,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 * pow((slope_length / 72.6), s_l_exp);
dseg_put(&l_s, &LS, r, c);
return 0;
More information about the grass-commit
mailing list