[GRASS-SVN] r43387 - grass-addons/raster/r.stream.extract
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Sep 1 14:36:31 EDT 2010
Author: mmetz
Date: 2010-09-01 18:36:31 +0000 (Wed, 01 Sep 2010)
New Revision: 43387
Modified:
grass-addons/raster/r.stream.extract/description.html
grass-addons/raster/r.stream.extract/load.c
grass-addons/raster/r.stream.extract/local_proto.h
grass-addons/raster/r.stream.extract/main.c
grass-addons/raster/r.stream.extract/streams.c
Log:
remove weight option and update manual
Modified: grass-addons/raster/r.stream.extract/description.html
===================================================================
--- grass-addons/raster/r.stream.extract/description.html 2010-09-01 13:18:36 UTC (rev 43386)
+++ grass-addons/raster/r.stream.extract/description.html 2010-09-01 18:36:31 UTC (rev 43387)
@@ -26,21 +26,9 @@
should be identical to the <em>accumulation map</em>. Flow direction is
first calculated from <em>elevation</em> and then adjusted to
<em>accumulation</em>. It is not necessary to provide <em>accumulation</em>
-as the number of cells, it can also be (adjusted) total contributing
-area in square meters or any other unit.
+as the number of cells, it can also be the optionally adjusted or
+weighed total contributing area in square meters or any other unit.
<p>
-<dt><em>weight</em>
-<dd>Input map, optional: Map used as <em>weight</em> for accumulation
-values. If local accumulation multiplied by local weight reaches or
-exceeds treshold, a new stream is initiated. If both
-<em>accumulation</em> and <em>weight</em> are given, memory can be safed
-by multiplying accumulation with <em>weight</em> first using
-<a href="r.mapcalc.html">r.mapcalc</a>, and then only giving the new
-accumulation map as input, <em>weight</em> is now already built in. This
-option allows e.g. to decrease the number of streams in dry areas and
-increase the number of streams in wet areas by setting <em>weight</em>
-to smaller than 1 in dry areas and larger than 1 in wet areas.
-<p>
<dt><em>threshold</em>
<dd>Required: <em>threshold</em> for stream initiation by overland flow:
the minumum (optionally modifed) flow accumulation value that will initiate
@@ -106,12 +94,33 @@
MFD (FD8) after Holmgren 1994, as for
<a href="r.watershed.html">r.watershed</a>. The <em>threshold</em>
option determines the number of streams and detail of stream networks.
-Whenever (optionally weighed) flow accumulation reaches
-<em>threshold</em>, a new stream is started and traced downstream to its
-outlet point. As for <a href="r.watershed.html">r.watershed</a>,
-flow accumulation is calculated as the number of cells draining through
-a cell.
+Whenever flow accumulation reaches <em>threshold</em>, a new stream is
+started and traced downstream to its outlet point. As for
+<a href="r.watershed.html">r.watershed</a>, flow accumulation is
+calculated as the number of cells draining through a cell.
+<h4>Weighed flow accumulation</h4>
+Flow accumulation can be calculated first, e.g. with
+<a href="r.watershed.html">r.watershed</a>, and then modified before
+using it as input for <em>r.stream.extract</em>. In its general form, a
+weighed accumulation map is generated by first creating a weighing map
+and then multiplying the accumulation map with the weighing map using
+<a href="r.mapcalc.html">r.mapcalc</a>. It is highly recommended to
+evaluate the weighed flow accumulation map first, before using it as
+input for <em>r.stream.extract</em>.
+<p>
+This allows e.g. to decrease the number of streams in dry areas and
+increase the number of streams in wet areas by setting <em>weight</em>
+to smaller than 1 in dry areas and larger than 1 in wet areas.
+<p>
+Another possibility is to restrict channel initiation to valleys
+determined from terrain morphology. Valleys can be determined with
+<a href="r.param.scale.html">r.param.scale</a> param=crosc
+(cross-sectional or tangential curvature). Curvature values < 0
+indicate concave features, i.e. valleys. The size of the processing
+window determines whether narrow or broad valleys will be identified
+(See example below).
+
<h4>Stream output</h4>
The <em>stream_rast</em> output raster and vector contains stream
segments with unique IDs. Note that these IDs are different from the IDs
@@ -120,11 +129,67 @@
segment, at confluences and at stream network outlet locations.
<p>
+<h2>EXAMPLE</h2>
+This example is based on the elevation map <em>elevation.10m</em> in the
+sample dataset spearfish60 and uses valleys determined with
+<a href="r.param.scale.html">r.param.scale</a> to weigh an accumulation
+map produced with <a href="r.watershed.html">r.watershed</a>.
+
+<pre>
+# set region
+g.region -p rast=elevation.10m at PERMANENT
+
+# calculate flow accumulation
+r.watershed ele=elevation.10m at PERMANENT acc=elevation.10m.acc -f
+
+# curvature to get narrow valleys
+r.param.scale input="elevation.10m at PERMANENT" output="tangential_curv_5" size=5 param=crosc
+
+# curvature to get a bit broader valleys
+r.param.scale input="elevation.10m at PERMANENT" output="tangential_curv_7" size=7 param=crosc
+
+# curvature to get broad valleys
+r.param.scale input="elevation.10m at PERMANENT" output="tangential_curv_11" size=11 param=crosc
+
+# create weight map
+r.mapcalc "weight = if(tangential_curv_5 < 0, -100 * tangential_curv_5, \
+ if(tangential_curv_7 < 0, -100 * tangential_curv_7, \
+ if(tangential_curv_11 < 0, -100 * tangential_curv_11, 0.000001)))"
+
+# weigh accumulation map
+r.mapcalc "elevation.10m.acc.weighed = elevation.10m.acc * weight"
+
+# copy color table from original accumulation map
+r.colors map=elevation.10m.acc.weighed raster=elevation.10m.acc
+</pre>
+
+Display both the original and the weighed accumulation map.
+<br>
+Compare them and proceed if the weighed accumulation map makes sense.
+
+<pre>
+# extract streams
+r.stream.extract elevation=elevation.10m at PERMANENT \
+ accumulation=elevation.10m.acc.weighed \
+ threshold=1000 \
+ stream_rast=elevation.10m.streams
+
+# extract streams using the original accumulation map
+r.stream.extract elevation=elevation.10m at PERMANENT \
+ accumulation=elevation.10m.acc \
+ threshold=1000 \
+ stream_rast=elevation.10m.streams.noweight
+</pre>
+
+Now display both stream maps and decide which one is more realistic.
+
<h2>SEE ALSO</h2>
<em>
<a href="r.watershed.html">r.watershed</a>,
<a href="r.terraflow.html">r.terraflow</a>,
+<a href="r.param.scale.html">r.param.scale</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
<a href="r.thin.html">r.thin</a>,
<a href="r.to.vect.html">r.to.vect</a>
</em>
Modified: grass-addons/raster/r.stream.extract/load.c
===================================================================
--- grass-addons/raster/r.stream.extract/load.c 2010-09-01 13:18:36 UTC (rev 43386)
+++ grass-addons/raster/r.stream.extract/load.c 2010-09-01 18:36:31 UTC (rev 43387)
@@ -23,23 +23,22 @@
* gets start points for A* Search
* start points are edges
*/
-int load_maps(int ele_fd, int acc_fd, int weight_fd)
+int load_maps(int ele_fd, int acc_fd)
{
int r, c, thisindex;
char asp_value, *aspp;
- void *ele_buf, *ptr, *acc_buf = NULL, *acc_ptr = NULL, *weight_buf =
- NULL, *weight_ptr = NULL;
+ void *ele_buf, *ptr, *acc_buf = NULL, *acc_ptr = NULL;
CELL *loadp, ele_value;
DCELL dvalue;
int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
int r_nbr, c_nbr, ct_dir;
int is_worked;
- size_t ele_size, acc_size = 0, weight_size = 0;
- int ele_map_type, acc_map_type = 0, weight_map_type = 0;
- DCELL *accp, *weightp;
+ size_t ele_size, acc_size = 0;
+ int ele_map_type, acc_map_type = 0;
+ DCELL *accp;
- if (acc_fd < 0 && weight_fd < 0)
+ if (acc_fd < 0)
G_message(_("load elevation map and get start points"));
else
G_message(_("load input maps and get start points"));
@@ -66,16 +65,6 @@
}
}
- if (weight_fd >= 0) {
- weight_map_type = G_get_raster_map_type(weight_fd);
- weight_size = G_raster_size(weight_map_type);
- weight_buf = G_allocate_raster_buf(weight_map_type);
- if (weight_buf == NULL) {
- G_warning(_("could not allocate memory"));
- return -1;
- }
- }
-
ele_scale = 1;
if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
ele_scale = 1000; /* should be enough to do the trick */
@@ -85,7 +74,6 @@
loadp = ele;
accp = acc;
- weightp = accweight;
aspp = asp;
G_debug(1, "start loading %d rows, %d cols", nrows, ncols);
@@ -107,15 +95,6 @@
acc_ptr = acc_buf;
}
- if (weight_fd >= 0) {
- if (G_get_raster_row(weight_fd, weight_buf, r, weight_map_type) <
- 0) {
- G_warning(_("could not read raster maps at row <%d>"), r);
- return -1;
- }
- weight_ptr = weight_buf;
- }
-
for (c = 0; c < ncols; c++) {
FLAG_UNSET(worked, r, c);
@@ -127,8 +106,6 @@
FLAG_SET(in_list, r, c);
G_set_c_null_value(loadp, 1);
*accp = 0;
- if (weight_fd >= 0)
- *weightp = 0;
}
else {
if (ele_map_type == CELL_TYPE) {
@@ -157,17 +134,6 @@
*accp = *((DCELL *) acc_ptr);
}
}
- if (weight_fd >= 0) {
- if (weight_map_type == CELL_TYPE) {
- *weightp = *((CELL *) weight_ptr);
- }
- else if (weight_map_type == FCELL_TYPE) {
- *weightp = *((FCELL *) weight_ptr);
- }
- else if (weight_map_type == DCELL_TYPE) {
- *weightp = *((DCELL *) weight_ptr);
- }
- }
n_points++;
}
@@ -179,10 +145,6 @@
aspp++;
if (acc_fd >= 0)
acc_ptr = G_incr_void_ptr(acc_ptr, acc_size);
- if (weight_fd >= 0) {
- weight_ptr = G_incr_void_ptr(weight_ptr, weight_size);
- weightp++;
- }
}
}
G_percent(nrows, nrows, 1); /* finish it */
@@ -195,11 +157,6 @@
G_free(acc_buf);
}
- if (weight_fd >= 0) {
- G_close_cell(weight_fd);
- G_free(weight_buf);
- }
-
astar_pts =
(unsigned int *)G_malloc((n_points + 1) *
sizeof(unsigned int));
Modified: grass-addons/raster/r.stream.extract/local_proto.h
===================================================================
--- grass-addons/raster/r.stream.extract/local_proto.h 2010-09-01 13:18:36 UTC (rev 43386)
+++ grass-addons/raster/r.stream.extract/local_proto.h 2010-09-01 18:36:31 UTC (rev 43387)
@@ -37,7 +37,7 @@
unsigned int n_stream_nodes, n_alloc_nodes;
struct point *outlets;
unsigned int n_outlets, n_alloc_outlets;
-DCELL *acc, *accweight;
+DCELL *acc;
CELL *ele;
char *asp;
CELL *stream;
@@ -50,7 +50,7 @@
struct RB_TREE *draintree;
/* load.c */
-int load_maps(int, int, int);
+int load_maps(int, int);
/* do_astar.c */
int do_astar(void);
@@ -58,11 +58,14 @@
/* streams.c */
int do_accum(double);
-int extract_streams(double, double, int, int);
+int extract_streams(double, double, int);
/* thin.c */
int thin_streams(void);
+/* basins.c */
+int basin_borders(void);
+
/* del_streams.c */
int del_streams(int);
int seg_length(int, unsigned int *);
Modified: grass-addons/raster/r.stream.extract/main.c
===================================================================
--- grass-addons/raster/r.stream.extract/main.c 2010-09-01 13:18:36 UTC (rev 43386)
+++ grass-addons/raster/r.stream.extract/main.c 2010-09-01 18:36:31 UTC (rev 43387)
@@ -27,7 +27,7 @@
{
struct
{
- struct Option *ele, *acc, *weight;
+ struct Option *ele, *acc;
struct Option *threshold, *d8cut;
struct Option *mont_exp;
struct Option *min_stream_length;
@@ -39,7 +39,7 @@
struct Option *dir_rast;
} output;
struct GModule *module;
- int ele_fd, acc_fd, weight_fd;
+ int ele_fd, acc_fd;
double threshold, d8cut, mont_exp;
int min_stream_length = 0;
char *mapset;
@@ -63,13 +63,6 @@
input.acc->description =
_("Stream extraction will use provided accumulation instead of calculating it anew");
- input.weight = G_define_standard_option(G_OPT_R_INPUT);
- input.weight->key = "weight";
- input.weight->label = _("Weight map for accumulation");
- input.weight->required = NO;
- input.weight->description =
- _("Map used as weight for flow accumulation when initiating streams");
-
input.threshold = G_define_option();
input.threshold->key = "threshold";
input.threshold->label = _("Minimum flow accumulation for streams");
@@ -144,12 +137,6 @@
G_fatal_error(_("Raster map <%s> not found"), input.acc->answer);
}
- if (input.weight->answer) {
- if (!G_find_cell(input.weight->answer, ""))
- G_fatal_error(_("Raster map <%s> not found"),
- input.weight->answer);
- }
-
/* threshold makes sense */
threshold = atof(input.threshold->answer);
if (threshold <= 0)
@@ -191,7 +178,8 @@
/* Check for some output map */
if ((output.stream_rast->answer == NULL)
- && (output.stream_vect->answer == NULL)) {
+ && (output.stream_vect->answer == NULL)
+ && (output.dir_rast->answer == NULL)) {
G_fatal_error(_("Sorry, you must choose at least one output map."));
}
@@ -215,16 +203,6 @@
else
acc_fd = -1;
- if (input.weight->answer) {
- mapset = G_find_cell2(input.weight->answer, "");
- weight_fd = G_open_cell_old(input.weight->answer, mapset);
- if (weight_fd < 0)
- G_fatal_error(_("could not open input map %s"),
- input.weight->answer);
- }
- else
- weight_fd = -1;
-
/* set global variables */
nrows = G_window_rows();
ncols = G_window_cols();
@@ -235,13 +213,9 @@
ele = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
asp = (char *) G_malloc(nrows * ncols * sizeof(char));
acc = (DCELL *) G_malloc(nrows * ncols * sizeof(DCELL));
- if (input.weight->answer)
- accweight = (DCELL *) G_malloc(nrows * ncols * sizeof(DCELL));
- else
- accweight = NULL;
/* load maps */
- if (load_maps(ele_fd, acc_fd, weight_fd) < 0)
+ if (load_maps(ele_fd, acc_fd) < 0)
G_fatal_error(_("could not load input map(s)"));
/********************/
@@ -260,12 +234,10 @@
stream = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
if (extract_streams
- (threshold, mont_exp, (input.weight->answer != NULL), min_stream_length) < 0)
+ (threshold, mont_exp, min_stream_length) < 0)
G_fatal_error(_("could not extract streams"));
G_free(acc);
- if (input.weight->answer)
- G_free(accweight);
/* thin streams */
if (thin_streams() < 0)
Modified: grass-addons/raster/r.stream.extract/streams.c
===================================================================
--- grass-addons/raster/r.stream.extract/streams.c 2010-09-01 13:18:36 UTC (rev 43386)
+++ grass-addons/raster/r.stream.extract/streams.c 2010-09-01 18:36:31 UTC (rev 43387)
@@ -421,7 +421,7 @@
/*
* extracts streams for threshold, accumulation is provided
*/
-int extract_streams(double threshold, double mont_exp, int use_weight, int min_length)
+int extract_streams(double threshold, double mont_exp, int min_length)
{
int r, c, dr, dc;
CELL is_swale, ele_val, ele_nbr;
@@ -646,7 +646,7 @@
/* set main drainage direction to A* path if possible */
if (mfd_cells > 0 && max_side != np_side) {
nindex = INDEX(dr, dc);
- if (fabs(acc[nindex]) >= max_acc) {
+ if (fabs(acc[nindex] >= max_acc)) {
max_acc = fabs(acc[nindex]);
r_max = dr;
c_max = dc;
@@ -669,10 +669,6 @@
/* start new stream */
/**********************/
- /* weight map has precedence over Montgomery */
- if (use_weight) {
- value *= accweight[thisindex];
- }
/* Montgomery's stream initiation acc * (tan(slope))^mont_exp */
/* uses whatever unit is accumulation */
if (mont_exp > 0) {
More information about the grass-commit
mailing list