[GRASS-SVN] r39281 - grass-addons/raster/r.stream.extract
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Sep 22 10:42:16 EDT 2009
Author: mmetz
Date: 2009-09-22 10:42:16 -0400 (Tue, 22 Sep 2009)
New Revision: 39281
Added:
grass-addons/raster/r.stream.extract/Makefile
grass-addons/raster/r.stream.extract/close.c
grass-addons/raster/r.stream.extract/description.html
grass-addons/raster/r.stream.extract/do_astar.c
grass-addons/raster/r.stream.extract/flag.h
grass-addons/raster/r.stream.extract/flag_clr_all.c
grass-addons/raster/r.stream.extract/flag_create.c
grass-addons/raster/r.stream.extract/flag_destroy.c
grass-addons/raster/r.stream.extract/flag_get.c
grass-addons/raster/r.stream.extract/flag_set.c
grass-addons/raster/r.stream.extract/flag_unset.c
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/rbtree.c
grass-addons/raster/r.stream.extract/rbtree.h
grass-addons/raster/r.stream.extract/streams.c
grass-addons/raster/r.stream.extract/thin.c
Log:
new r.stream.extract
Added: grass-addons/raster/r.stream.extract/Makefile
===================================================================
--- grass-addons/raster/r.stream.extract/Makefile (rev 0)
+++ grass-addons/raster/r.stream.extract/Makefile 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,12 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.stream.extract
+
+LIBES = $(VECTLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/raster/r.stream.extract/close.c
===================================================================
--- grass-addons/raster/r.stream.extract/close.c (rev 0)
+++ grass-addons/raster/r.stream.extract/close.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,320 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/dbmi.h>
+#include <grass/Vect.h>
+#include "local_proto.h"
+
+int close_streamvect(char *stream_vect)
+{
+ int i, r, c, r_nbr, c_nbr, done;
+ int stream_id, next_node;
+ unsigned int thisindex;
+ struct sstack
+ {
+ int stream_id;
+ int next_trib;
+ } *nodestack;
+ int top = 0, stack_step = 1000;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ struct ddir draindir, *founddir;
+ struct Map_info Out;
+ static struct line_pnts *Points;
+ struct line_cats *Cats;
+ dbDriver *driver;
+ dbHandle handle;
+ dbString table_name, dbsql, valstr;
+ struct field_info *Fi;
+ char *cat_col_name = "cat", buf[2000];
+ struct Cell_head window;
+ double north_offset, west_offset, ns_res, ew_res;
+ int next_cat;
+
+ G_message(_("write vector maps"));
+
+ if (0 > Vect_open_new(&Out, stream_vect, 0)) {
+ G_fatal_error(_("Unable to create vector map <%s>"), stream_vect);
+ }
+
+ nodestack = (struct sstack *)G_malloc(stack_step * sizeof(struct sstack));
+
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ G_get_set_window(&window);
+ ns_res = window.ns_res;
+ ew_res = window.ew_res;
+ north_offset = window.north - 0.5 * ns_res;
+ west_offset = window.west + 0.5 * ew_res;
+
+ next_cat = n_stream_nodes + 1;
+
+ for (i = 0; i < n_outlets; i++, next_cat++) {
+ G_percent(i, n_outlets, 2);
+ r = outlets[i].r;
+ c = outlets[i].c;
+ thisindex = INDEX(r, c);
+ stream_id = stream[thisindex];
+
+ Vect_reset_line(Points);
+ Vect_reset_cats(Cats);
+
+ /* outlet */
+ Vect_cat_set(Cats, 1, stream_id);
+ Vect_cat_set(Cats, 2, 2);
+ Vect_append_point(Points, west_offset + c * ew_res,
+ north_offset - r * ns_res, 0);
+ Vect_write_line(&Out, GV_POINT, Points, Cats);
+
+ /* add root node to stack */
+ G_debug(3, "add root node");
+ top = 0;
+ nodestack[top].stream_id = stream_id;
+ nodestack[top].next_trib = 0;
+
+ /* depth first post order traversal */
+ G_debug(3, "traverse");
+ while (top >= 0) {
+
+ done = 1;
+ stream_id = nodestack[top].stream_id;
+ G_debug(3, "stream_id %d", stream_id);
+ if (nodestack[top].next_trib < stream_node[stream_id].n_trib) {
+ /* add to stack */
+ next_node =
+ stream_node[stream_id].trib[nodestack[top].next_trib];
+ G_debug(3, "add to stack: next %d, trib %d, n trib %d",
+ next_node, nodestack[top].next_trib,
+ stream_node[stream_id].n_trib);
+ nodestack[top].next_trib++;
+ top++;
+ if (top >= stack_step) {
+ /* need more space */
+ stack_step += 1000;
+ nodestack =
+ (struct sstack *)G_realloc(nodestack,
+ stack_step *
+ sizeof(struct sstack));
+ }
+ nodestack[top].next_trib = 0;
+ nodestack[top].stream_id = next_node;
+ done = 0;
+ G_debug(3, "go further down");
+ }
+ if (done) {
+ G_debug(3, "write stream segment");
+
+ Vect_reset_line(Points);
+ Vect_reset_cats(Cats);
+
+ r_nbr = stream_node[stream_id].r;
+ c_nbr = stream_node[stream_id].c;
+ draindir.pos = INDEX(r_nbr, c_nbr);
+
+ Vect_cat_set(Cats, 1, stream_id);
+ if (stream_node[stream_id].n_trib == 0)
+ Vect_cat_set(Cats, 2, 0);
+ else
+ Vect_cat_set(Cats, 2, 1);
+
+ Vect_append_point(Points, west_offset + c_nbr * ew_res,
+ north_offset - r_nbr * ns_res, 0);
+
+ Vect_write_line(&Out, GV_POINT, Points, Cats);
+
+ while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
+ r_nbr = r_nbr + asp_r[(int)founddir->dir];
+ c_nbr = c_nbr + asp_c[(int)founddir->dir];
+
+ if (stream[founddir->pos] != stream_id) {
+ /* append first point of parent stream */
+ Vect_append_point(Points,
+ west_offset + c_nbr * ew_res,
+ north_offset - r_nbr * ns_res, 0);
+ break;
+ }
+ draindir.pos = INDEX(r_nbr, c_nbr);
+ if (stream[INDEX(r_nbr, c_nbr)] <= 0)
+ G_fatal_error("stream id not set");
+
+ Vect_append_point(Points, west_offset + c_nbr * ew_res,
+ north_offset - r_nbr * ns_res, 0);
+ }
+
+ Vect_write_line(&Out, GV_LINE, Points, Cats);
+
+ top--;
+ }
+ }
+ }
+ G_percent(n_outlets, n_outlets, 1); /* finish it */
+
+ G_message(_("write vector attribute table"));
+
+ /* Prepeare strings for use in db_* calls */
+ db_init_string(&dbsql);
+ db_init_string(&valstr);
+ db_init_string(&table_name);
+ db_init_handle(&handle);
+
+ /* Preparing database for use */
+ /* Create database for new vector map */
+ Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
+ driver = db_start_driver_open_database(Fi->driver, Fi->database);
+ if (driver == NULL) {
+ G_fatal_error(_("Unable to start driver <%s>"), Fi->driver);
+ }
+
+ G_debug(1, "table: %s", Fi->table);
+ G_debug(1, "driver: %s", Fi->driver);
+ G_debug(1, "database: %s", Fi->database);
+
+ sprintf(buf,
+ "create table %s (%s integer, stream_type varchar(20), type_code integer)",
+ Fi->table, cat_col_name);
+ db_set_string(&dbsql, buf);
+
+ if (db_execute_immediate(driver, &dbsql) != DB_OK) {
+ db_close_database(driver);
+ db_shutdown_driver(driver);
+ G_fatal_error(_("Cannot create table: %s"), db_get_string(&dbsql));
+ }
+
+ if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
+ G_warning(_("Cannot create index"));
+
+ if (db_grant_on_table(driver, Fi->table, DB_PRIV_SELECT,
+ DB_GROUP | DB_PUBLIC) != DB_OK)
+ G_fatal_error(_("Cannot grant privileges on table %s"), Fi->table);
+
+ db_begin_transaction(driver);
+
+ /* stream nodes */
+ for (i = 1; i <= n_stream_nodes; i++) {
+
+ sprintf(buf, "insert into %s values ( %d, \'%s\', %d )",
+ Fi->table, i,
+ (stream_node[i].n_trib > 0 ? "start" : "intermediate"),
+ (stream_node[i].n_trib > 0));
+
+ db_set_string(&dbsql, buf);
+
+ if (db_execute_immediate(driver, &dbsql) != DB_OK) {
+ db_close_database(driver);
+ db_shutdown_driver(driver);
+ G_fatal_error(_("Cannot insert new row: %s"),
+ db_get_string(&dbsql));
+ }
+ }
+
+ G_message(_("close vector"));
+
+ db_commit_transaction(driver);
+ db_close_database_shutdown_driver(driver);
+
+ Vect_map_add_dblink(&Out, 1, NULL, Fi->table,
+ cat_col_name, Fi->database, Fi->driver);
+
+ G_message(_("close vector"));
+
+ Vect_hist_command(&Out);
+ Vect_build(&Out);
+ Vect_close(&Out);
+
+ G_free(nodestack);
+
+ return 1;
+}
+
+
+int close_maps(char *stream_rast, char *stream_vect, char *dir_rast)
+{
+ int stream_fd, dir_fd, r, c, i;
+ CELL *cell_buf1, *cell_buf2;
+ unsigned int thisindex;
+ struct History history;
+ struct ddir draindir, *founddir;
+
+ /* cheating... */
+ stream_fd = dir_fd = -1;
+ cell_buf1 = cell_buf2 = NULL;
+
+ G_message(_("write raster maps"));
+
+ /* write requested output rasters */
+ if (stream_rast) {
+ stream_fd = G_open_raster_new(stream_rast, CELL_TYPE);
+ cell_buf1 = G_allocate_cell_buf();
+ }
+ if (dir_rast) {
+ dir_fd = G_open_raster_new(dir_rast, CELL_TYPE);
+ cell_buf2 = G_allocate_cell_buf();
+ }
+
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 2);
+ if (stream_rast)
+ G_set_c_null_value(cell_buf1, ncols); /* reset row to all NULL */
+ if (dir_rast)
+ G_set_c_null_value(cell_buf2, ncols); /* reset row to all NULL */
+
+ for (c = 0; c < ncols; c++) {
+ thisindex = INDEX(r, c);
+ if (stream[thisindex] > 0) {
+ if (stream_rast)
+ cell_buf1[c] = stream[thisindex];
+ if (dir_rast) {
+ draindir.pos = thisindex;
+ if ((founddir =
+ rbtree_find(draintree, &draindir)) != NULL) {
+ cell_buf2[c] = founddir->dir;
+ }
+ else {
+ cell_buf2[c] = 0;
+ }
+ }
+ }
+ }
+ if (stream_rast)
+ G_put_raster_row(stream_fd, cell_buf1, CELL_TYPE);
+ if (dir_rast)
+ G_put_raster_row(dir_fd, cell_buf2, CELL_TYPE);
+ }
+ G_percent(nrows, nrows, 2); /* finish it */
+
+ if (stream_rast) {
+ G_close_cell(stream_fd);
+ G_free(cell_buf1);
+ G_short_history(stream_rast, "raster", &history);
+ G_command_history(&history);
+ G_write_history(stream_rast, &history);
+ }
+ if (dir_rast) {
+ G_close_cell(dir_fd);
+ G_free(cell_buf2);
+ G_short_history(dir_rast, "raster", &history);
+ G_command_history(&history);
+ G_write_history(dir_rast, &history);
+ }
+
+ /* close stream vector */
+ if (stream_vect) {
+ if (close_streamvect(stream_vect) < 0)
+ G_fatal_error(_("Unable to write vector map <%s>"), stream_vect);
+ }
+
+ /* rearranging desk chairs on the Titanic... */
+ rbtree_destroy(draintree);
+ G_free(outlets);
+
+ /* free stream nodes */
+ for (i = 1; i <= n_stream_nodes; i++) {
+ if (stream_node[i].n_alloc > 0) {
+ G_free(stream_node[i].trib);
+ G_free(stream_node[i].acc);
+ }
+ }
+ G_free(stream_node);
+
+ return 1;
+}
Added: grass-addons/raster/r.stream.extract/description.html
===================================================================
--- grass-addons/raster/r.stream.extract/description.html (rev 0)
+++ grass-addons/raster/r.stream.extract/description.html 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,151 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.stream.extract</em> extracts streams in both raster and vector
+format from a required input <em>elevation</em> map and optional input
+<em>accumulation map</em>.
+
+<h2>OPTIONS</h2>
+
+<dl>
+<dt><em>elevation</em>
+
+<dd>Input map, required: 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.
+<p>
+<dt><em>accumulation</em>
+<dd>Input map, optional: Accumulation values of the provided
+<em>accumulation</em> map are used and not calculated from the input
+<em>elevation</em> map. If <em>accumulation</em> is given,
+<em>elevation</em> must be exactly the same map used to calculate
+<em>accumulation</em>. If <em>accumulation</em> was calculated with
+<a href="r.terraflow.html">r.terraflow</a>, the filled elevation output
+of r.terraflow must be used. Further on, the current region's resolution
+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.
+<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
+a new stream. If Montgomery's method for channel initiation is used, the
+cell value of the accumulation input map is multiplied by
+(tan(local slope))<sup>mexp</sup> and then compared to <em>threshold</em>.
+<p>
+<dt><em>d8cut</em>
+<dd>Minimum amount of overland flow (accumulation) when SFD (D8) will be
+used instead of MFD (FD8) to calculate flow accumulation. Only applies
+if no accumulation map is provided. Setting to 0 disables MFD completely.
+<p>
+<dt><em>mexp</em>
+<dd>Use the method of Montgomery and Foufoula-Georgiou (1993) to
+initiate a stream with exponent <em>mexp</em>. The cell value of the
+accumulation input map is multiplied by (tan(local slope))<sup>mexp</sup>
+and then compared to <em>threshold</em>. If threshold is reached or
+exceeded, a new stream is initiated. The default value 0 disables
+Montgomery. Montgomery and Foufoula-Georgiou (1993) generally recommend
+to use 2.0 as exponent. <em>mexp</em> values closer to 0 will produce
+streams more similar to streams extracted with Montgomery disabled.
+Larger <em>mexp</em> values decrease the number of streams in flat areas
+and increase the number of streams in steep areas. If <em>weight</em> is
+given, the weight is applied first.
+
+<p>
+<dt><em>stream_rast</em>
+<dd>Output raster map with extracted streams. Cell values encode unique
+ID for each stream segment.
+<p>
+<dt><em>stream_vect</em>
+<dd>Output vector map with extracted stream segments and points. Points
+are written at the start location of each stream segment and at the
+outlet of a stream network. In layer 1, categories are unique IDs,
+identical to the cell value of the raster output. The attribute table
+for layer 1 holds information about the type of stream segment: start
+segment, or intermediate segment with tributaries. Columns are cat int,
+stream_type varchar(), type_code int. The encoding for type_code is 0 =
+start, 1 = intermediate. In layer 2, categories are identical to
+type_code in layer 1 with additional category 2 = outlet for outlet
+points. Points with category 1 = intermediate in layer 2 are at the
+location of confluences.
+<p>
+<dt><em>direction</em>
+<dd>Output raster map with flow direction for extracted streams. Flow
+direction is of D8 type with a range of 1 to 8. Multiplying values with
+45 gives degrees CCW from East. Flow direction was adjusted during
+thinning, taking shortcuts and skipping cells that were eliminated by
+the thinning procedure. Non-stream cells are set to NULL. A full,
+corrected flow direction map can be created by patching the
+<em>direction</em> output map with the flow direction map of r.watershed.
+</dl>
+
+<h2>NOTES</h2>
+
+<h4>Stream extraction</h4>
+If no accumulation input map is provided, flow accumulation is
+determined with a hydrological anaylsis similar to
+<a href="r.watershed.html">r.watershed</a>. The algorithm is
+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.
+
+<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
+assigned by <a href="r.watershed.html">r.watershed</a>. The vector
+output also contains points at the location of the start of a stream
+segment, at confluences and at stream network outlet locations.
+<p>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.terraflow.html">r.terraflow</a>,
+<a href="r.thin.html">r.thin</a>,
+<a href="r.to.vect.html">r.to.vect</a>
+</em>
+
+<h2>REFERENCES</h2>
+Ehlschlaeger, C. (1989). <i>Using the A<sup>T</sup> Search Algorithm
+to Develop Hydrologic Models from Digital Elevation Data</i>,
+<b>Proceedings of International Geographic Information Systems (IGIS)
+Symposium '89</b>, pp 275-281 (Baltimore, MD, 18-19 March 1989).<br>
+URL: <a href="http://faculty.wiu.edu/CR-Ehlschlaeger2/older/IGIS/paper.html">
+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), pp 327-334.<br>
+DOI: <a href="http://dx.doi.org/10.1002/hyp.3360080405">10.1002/hyp.3360080405</a>
+
+<p>
+Montgomery, D.R., Foufoula-Georgiou, E. (1993). <i>Channel network source
+representation using digital elevation models.</i>
+<b>Water Resources Research</b> Vol 29(12), pp 3925-3934.
+
+<h2>AUTHOR</h2>
+Markus Metz
+
+<p><i>Last changed: $Date: 2008-05-16 21:09:06 +0200 (Fri, 16 May 2008) $</i>
Added: grass-addons/raster/r.stream.extract/do_astar.c
===================================================================
--- grass-addons/raster/r.stream.extract/do_astar.c (rev 0)
+++ grass-addons/raster/r.stream.extract/do_astar.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,229 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+unsigned int get_parent(unsigned int c)
+{
+ return (unsigned int)((c) - 2) / 3 + 1;
+}
+
+unsigned int get_child(unsigned int p)
+{
+ return (unsigned int)(p) * 3 - 1;
+}
+
+unsigned int heap_drop(void);
+
+int do_astar(void)
+{
+ int r, c, r_nbr, c_nbr, ct_dir;
+ struct ast_point astp;
+ int count, is_in_list;
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ CELL ele_val, ele_up;
+ char asp_val;
+ unsigned int thisindex, nindex;
+
+ count = 0;
+
+ first_cum = n_points;
+
+ G_message(_("A* Search..."));
+
+ while (heap_size > 0) {
+ G_percent(count++, n_points, 1);
+ if (count > n_points)
+ G_fatal_error("broken A* Search, %d surplus points", heap_size);
+
+ if (heap_size > n_points)
+ G_fatal_error
+ ("broken A* Search, too many points in heap %d, should be %d",
+ heap_size, n_points);
+
+ astp = astar_pts[1];
+
+ heap_drop();
+
+ /* set flow accumulation order */
+ astar_pts[first_cum] = astp;
+ first_cum--;
+
+ r = astp.r;
+ c = astp.c;
+
+ thisindex = INDEX(r, c);
+ ele_val = ele[thisindex];
+
+ 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) {
+
+ is_in_list = FLAG_GET(in_list, r_nbr, c_nbr);
+
+ if (is_in_list == 0) {
+ nindex = INDEX(r_nbr, c_nbr);
+ ele_up = ele[nindex];
+ asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
+ heap_add(r_nbr, c_nbr, ele_up, asp_val);
+ FLAG_SET(in_list, r_nbr, c_nbr);
+ }
+ }
+ } /* end neighbours */
+ FLAG_SET(worked, r, c);
+ } /* end A* search */
+
+ G_percent(n_points, n_points, 1); /* finish it */
+
+ flag_destroy(in_list);
+ G_free(astar_added);
+
+ return 1;
+}
+
+/*
+ * compare function for heap
+ * returns 1 if point1 < point2 else 0
+ */
+int heap_cmp(CELL ele1, unsigned int index1, CELL ele2, unsigned int index2)
+{
+ if (ele1 < ele2)
+ return 1;
+ else if (ele1 == ele2) {
+ if (index1 < index2)
+ return 1;
+ }
+
+ return 0;
+}
+
+int sift_up(unsigned int start, CELL elec)
+{
+ unsigned int child, child_added, parent, nindex;
+ CELL elep;
+ struct ast_point childp;
+
+ child = start;
+ child_added = astar_added[child];
+ childp = astar_pts[child];
+
+ while (child > 1) {
+ parent = get_parent(child);
+
+ nindex = INDEX(astar_pts[parent].r, astar_pts[parent].c);
+
+ elep = ele[nindex];
+
+ /* child < parent */
+ if (heap_cmp(elec, child_added, elep, astar_added[parent]) == 1) {
+ /* push parent point down */
+ astar_added[child] = astar_added[parent];
+ astar_pts[child] = astar_pts[parent];
+ child = parent;
+ }
+ else
+ /* no more sifting up, found new slot for child */
+ break;
+ }
+
+ /* set heap_index for child */
+ if (child < start) {
+ astar_added[child] = child_added;
+ astar_pts[child] = childp;
+ return 1;
+ }
+
+ return 0;
+}
+
+/*
+ * add item to heap
+ * returns heap_size
+ */
+unsigned int heap_add(int r, int c, CELL ele, char asp)
+{
+ FLAG_SET(in_list, r, c);
+
+ /* add point to next free position */
+
+ heap_size++;
+
+ if (heap_size > n_points)
+ G_fatal_error(_("heapsize too large"));
+
+ astar_added[heap_size] = nxt_avail_pt;
+ astar_pts[heap_size].r = r;
+ astar_pts[heap_size].c = c;
+ astar_pts[heap_size].asp = asp;
+
+ nxt_avail_pt++;
+
+ /* sift up: move new point towards top of heap */
+
+ sift_up(heap_size, ele);
+
+ return heap_size;
+}
+
+/*
+ * drop item from heap
+ * returns heap size
+ */
+unsigned int heap_drop(void)
+{
+ unsigned int child, childr, parent;
+ int i;
+ CELL elec, eler;
+
+ if (heap_size == 1) {
+ astar_added[1] = -1;
+ heap_size = 0;
+ return heap_size;
+ }
+
+ parent = 1;
+ while ((child = get_child(parent)) <= heap_size) {
+
+ elec = ele[INDEX(astar_pts[child].r, astar_pts[child].c)];
+
+ if (child < heap_size) {
+ childr = child + 1;
+ 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 = ele[INDEX(astar_pts[childr].r, astar_pts[childr].c)];
+
+ if (heap_cmp
+ (eler, astar_added[childr], elec,
+ astar_added[child]) == 1) {
+ child = childr;
+ elec = eler;
+ }
+ childr++;
+ }
+ }
+
+ /* move hole down */
+ astar_added[parent] = astar_added[child];
+ astar_pts[parent] = astar_pts[child];
+ parent = child;
+ }
+
+ /* hole is in lowest layer, move to heap end */
+ if (parent < heap_size) {
+ astar_added[parent] = astar_added[heap_size];
+ astar_pts[parent] = astar_pts[heap_size];
+
+ elec = ele[INDEX(astar_pts[parent].r, astar_pts[parent].c)];
+ /* sift up last swapped point, only necessary if hole moved to heap end */
+ sift_up(parent, elec);
+ }
+
+ /* the actual drop */
+ heap_size--;
+
+ return heap_size;
+}
Added: grass-addons/raster/r.stream.extract/flag.h
===================================================================
--- grass-addons/raster/r.stream.extract/flag.h (rev 0)
+++ grass-addons/raster/r.stream.extract/flag.h 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,77 @@
+#ifndef __FLAG_H__
+#define __FLAG_H__
+
+
+/* flag.[ch] is a set of routines which will set up an array of bits
+ ** that allow the programmer to "flag" cells in a raster map.
+ **
+ ** FLAG *
+ ** flag_create(nrows,ncols)
+ ** int nrows, ncols;
+ ** opens the structure flag.
+ ** The flag structure will be a two dimensional array of bits the
+ ** size of nrows by ncols. Will initalize flags to zero (unset).
+ **
+ ** flag_destroy(flags)
+ ** FLAG *flags;
+ ** closes flags and gives the memory back to the system.
+ **
+ ** flag_clear_all(flags)
+ ** FLAG *flags;
+ ** sets all values in flags to zero.
+ **
+ ** flag_unset(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ ** sets the value of (row, col) in flags to zero.
+ **
+ ** flag_set(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ ** will set the value of (row, col) in flags to one.
+ **
+ ** int
+ ** flag_get(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ ** returns the value in flags that is at (row, col).
+ **
+ ** idea by Michael Shapiro
+ ** code by Chuck Ehlschlaeger
+ ** April 03, 1989
+ */
+#define FLAG struct _flagsss_
+FLAG {
+ int nrows, ncols, leng;
+ unsigned char **array;
+};
+
+#define FLAG_UNSET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] &= ~(1<<((col) & 7))
+
+#define FLAG_SET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] |= (1<<((col) & 7))
+
+#define FLAG_GET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] & (1<<((col) & 7))
+
+/* flag_clr_all.c */
+int flag_clear_all(FLAG *);
+
+/* flag_create.c */
+FLAG *flag_create(int, int);
+
+/* flag_destroy.c */
+int flag_destroy(FLAG *);
+
+/* flag_get.c */
+int flag_get(FLAG *, int, int);
+
+/* flag_set.c */
+int flag_set(FLAG *, int, int);
+
+/* flag_unset.c */
+int flag_unset(FLAG *, int, int);
+
+
+#endif /* __FLAG_H__ */
Added: grass-addons/raster/r.stream.extract/flag_clr_all.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_clr_all.c (rev 0)
+++ grass-addons/raster/r.stream.extract/flag_clr_all.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,14 @@
+#include "flag.h"
+
+int flag_clear_all(FLAG * flags)
+{
+ register int r, c;
+
+ for (r = 0; r < flags->nrows; r++) {
+ for (c = 0; c < flags->leng; c++) {
+ flags->array[r][c] = 0;
+ }
+ }
+
+ return 0;
+}
Added: grass-addons/raster/r.stream.extract/flag_create.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_create.c (rev 0)
+++ grass-addons/raster/r.stream.extract/flag_create.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,27 @@
+#include <grass/gis.h>
+#include "flag.h"
+
+FLAG *flag_create(int nrows, int ncols)
+{
+ unsigned char *temp;
+ FLAG *new_flag;
+ register int i;
+
+ new_flag = (FLAG *) G_malloc(sizeof(FLAG));
+ new_flag->nrows = nrows;
+ new_flag->ncols = ncols;
+ new_flag->leng = (ncols + 7) / 8;
+ new_flag->array =
+ (unsigned char **)G_malloc(nrows * sizeof(unsigned char *));
+
+ temp =
+ (unsigned char *)G_malloc(nrows * new_flag->leng *
+ sizeof(unsigned char));
+
+ for (i = 0; i < nrows; i++) {
+ new_flag->array[i] = temp;
+ temp += new_flag->leng;
+ }
+
+ return (new_flag);
+}
Added: grass-addons/raster/r.stream.extract/flag_destroy.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_destroy.c (rev 0)
+++ grass-addons/raster/r.stream.extract/flag_destroy.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,11 @@
+#include <grass/gis.h>
+#include "flag.h"
+
+int flag_destroy(FLAG * flags)
+{
+ G_free(flags->array[0]);
+ G_free(flags->array);
+ G_free(flags);
+
+ return 0;
+}
Added: grass-addons/raster/r.stream.extract/flag_get.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_get.c (rev 0)
+++ grass-addons/raster/r.stream.extract/flag_get.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,6 @@
+#include "flag.h"
+
+int flag_get(FLAG * flags, int row, int col)
+{
+ return (flags->array[row][col >> 3] & (1 << (col & 7)));
+}
Added: grass-addons/raster/r.stream.extract/flag_set.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_set.c (rev 0)
+++ grass-addons/raster/r.stream.extract/flag_set.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,8 @@
+#include "flag.h"
+
+int flag_set(FLAG * flags, int row, int col)
+{
+ flags->array[row][col >> 3] |= (1 << (col & 7));
+
+ return 0;
+}
Added: grass-addons/raster/r.stream.extract/flag_unset.c
===================================================================
--- grass-addons/raster/r.stream.extract/flag_unset.c (rev 0)
+++ grass-addons/raster/r.stream.extract/flag_unset.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,8 @@
+#include "flag.h"
+
+int flag_unset(FLAG * flags, int row, int col)
+{
+ flags->array[row][col >> 3] &= ~(1 << (col & 7));
+
+ return 0;
+}
Added: grass-addons/raster/r.stream.extract/load.c
===================================================================
--- grass-addons/raster/r.stream.extract/load.c (rev 0)
+++ grass-addons/raster/r.stream.extract/load.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,283 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+/* need elevation map, do A* search on elevation like for r.watershed */
+
+int ele_round(double x)
+{
+ int n;
+
+ if (x >= 0.0)
+ n = x + .5;
+ else {
+ n = -x + .5;
+ n = -n;
+ }
+
+ return n;
+}
+
+/*
+ * loads elevation and optional flow accumulation map to memory and
+ * gets start points for A* Search
+ * start points are edges
+ */
+int load_maps(int ele_fd, int acc_fd, int weight_fd)
+{
+ int r, c, thisindex;
+ char asp_value;
+ void *ele_buf, *ptr, *acc_buf = NULL, *acc_ptr = NULL, *weight_buf =
+ NULL, *weight_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;
+ CELL *streamp;
+ DCELL *accp, *weightp;
+
+ if (acc_fd < 0 && weight_fd < 0)
+ G_message(_("load elevation map and get start points"));
+ else
+ G_message(_("load input maps and get start points"));
+
+ n_search_points = n_points = 0;
+
+ G_debug(1, "get buffers");
+ ele_map_type = G_get_raster_map_type(ele_fd);
+ ele_size = G_raster_size(ele_map_type);
+ ele_buf = G_allocate_raster_buf(ele_map_type);
+
+ if (ele_buf == NULL) {
+ G_warning(_("could not allocate memory"));
+ return -1;
+ }
+
+ if (acc_fd >= 0) {
+ acc_map_type = G_get_raster_map_type(acc_fd);
+ acc_size = G_raster_size(acc_map_type);
+ acc_buf = G_allocate_raster_buf(acc_map_type);
+ if (acc_buf == NULL) {
+ G_warning(_("could not allocate memory"));
+ return -1;
+ }
+ }
+
+ 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 */
+
+ worked = flag_create(nrows, ncols);
+ in_list = flag_create(nrows, ncols);
+
+ loadp = ele;
+ streamp = stream;
+ accp = acc;
+ weightp = accweight;
+
+ G_debug(1, "start loading %d rows, %d cols", nrows, ncols);
+ for (r = 0; r < nrows; r++) {
+
+ G_percent(r, nrows, 2);
+
+ if (G_get_raster_row(ele_fd, ele_buf, r, ele_map_type) < 0) {
+ G_warning(_("could not read raster maps at row <%d>"), r);
+ return -1;
+ }
+ ptr = ele_buf;
+
+ if (acc_fd >= 0) {
+ if (G_get_raster_row(acc_fd, acc_buf, r, acc_map_type) < 0) {
+ G_warning(_("could not read raster maps at row <%d>"), r);
+ return -1;
+ }
+ 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);
+ FLAG_UNSET(in_list, r, c);
+
+ *streamp = 0;
+
+ /* 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(loadp, 1);
+ *accp = 0;
+ if (weight_fd >= 0)
+ *weightp = 0;
+ }
+ else {
+ if (ele_map_type == CELL_TYPE) {
+ *loadp = *((CELL *) ptr);
+ }
+ else if (ele_map_type == FCELL_TYPE) {
+ dvalue = *((FCELL *) ptr);
+ dvalue *= ele_scale;
+ *loadp = ele_round(dvalue);
+ }
+ else if (ele_map_type == DCELL_TYPE) {
+ dvalue = *((DCELL *) ptr);
+ dvalue *= ele_scale;
+ *loadp = ele_round(dvalue);
+ }
+ if (acc_fd < 0)
+ *accp = 1;
+ else {
+ if (acc_map_type == CELL_TYPE) {
+ *accp = *((CELL *) acc_ptr);
+ }
+ else if (acc_map_type == FCELL_TYPE) {
+ *accp = *((FCELL *) acc_ptr);
+ }
+ else if (acc_map_type == DCELL_TYPE) {
+ *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++;
+ }
+
+ loadp++;
+ accp++;
+ streamp++;
+ ptr = G_incr_void_ptr(ptr, ele_size);
+ 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 */
+
+ G_close_cell(ele_fd);
+ G_free(ele_buf);
+
+ if (acc_fd >= 0) {
+ G_close_cell(acc_fd);
+ G_free(acc_buf);
+ }
+
+ if (weight_fd >= 0) {
+ G_close_cell(weight_fd);
+ G_free(weight_buf);
+ }
+
+ astar_pts =
+ (struct ast_point *)G_malloc((n_points + 1) *
+ sizeof(struct ast_point));
+
+ /* astar_heap will track astar_pts in ternary min-heap */
+ /* astar_heap is one-based */
+ astar_added =
+ (unsigned int *)G_malloc((n_points + 1) * sizeof(unsigned int));
+
+ nxt_avail_pt = heap_size = 0;
+
+ /* load edge cells to A* heap */
+ G_message(_("set edge points"));
+ loadp = ele;
+ for (r = 0; r < nrows; r++) {
+
+ G_percent(r, nrows, 2);
+ for (c = 0; c < ncols; c++) {
+
+ is_worked = FLAG_GET(worked, r, c);
+
+ if (is_worked)
+ continue;
+
+ if (r == 0 || r == nrows - 1 || c == 0 || c == ncols - 1) {
+
+ asp_value = 0;
+ if (r == 0 && c == 0)
+ asp_value = -7;
+ else if (r == 0 && c == ncols - 1)
+ asp_value = -5;
+ else if (r == nrows - 1 && c == 0)
+ asp_value = -1;
+ else if (r == nrows - 1 && c == ncols - 1)
+ asp_value = -3;
+ else if (r == 0)
+ asp_value = -2;
+ else if (c == 0)
+ asp_value = -4;
+ else if (r == nrows - 1)
+ asp_value = -6;
+ else if (c == ncols - 1)
+ asp_value = -8;
+
+ thisindex = INDEX(r, c);
+ ele_value = ele[thisindex];
+ heap_add(r, c, ele_value, asp_value);
+ FLAG_SET(in_list, r, c);
+ continue;
+ }
+
+ /* any neighbour NULL ? */
+ 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];
+
+ is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+
+ if (is_worked) {
+ asp_value = drain[r - r_nbr + 1][c - c_nbr + 1];
+ thisindex = INDEX(r, c);
+ ele_value = ele[thisindex];
+ heap_add(r, c, ele_value, asp_value);
+ FLAG_SET(in_list, r, c);
+
+ break;
+ }
+ }
+
+ }
+ }
+ G_percent(nrows, nrows, 2); /* finish it */
+
+ G_debug(1, "%d edge cells", heap_size);
+ G_debug(1, "%d non-NULL cells", n_points);
+
+ return 1;
+}
Added: grass-addons/raster/r.stream.extract/local_proto.h
===================================================================
--- grass-addons/raster/r.stream.extract/local_proto.h (rev 0)
+++ grass-addons/raster/r.stream.extract/local_proto.h 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,75 @@
+
+#ifndef __LOCAL_PROTO_H__
+#define __LOCAL_PROTO_H__
+
+#include "flag.h"
+#include "rbtree.h"
+
+#define INDEX(r, c) ((r) * ncols + (c))
+#define MAXDEPTH 1000 /* maximum supported tree depth of stream network */
+
+struct ddir
+{
+ int pos;
+ int dir;
+};
+
+struct ast_point
+{
+ int r, c;
+ char asp;
+};
+
+struct point
+{
+ int r, c;
+};
+
+struct snode
+{
+ int r, c;
+ int id;
+ int n_trib; /* number of tributaries */
+ int n_trib_total; /* number of all upstream stream segments */
+ int n_alloc; /* n allocated tributaries */
+ int *trib;
+ double *acc;
+} *stream_node;
+
+int nrows, ncols;
+struct ast_point *astar_pts;
+unsigned int n_search_points, n_points, nxt_avail_pt;
+unsigned int heap_size, *astar_added;
+unsigned int n_stream_nodes, n_alloc_nodes;
+struct point *outlets;
+unsigned int n_outlets, n_alloc_outlets;
+DCELL *acc, *accweight;
+CELL *ele;
+CELL *stream;
+int *strahler, *horton; /* strahler and horton order */
+FLAG *worked, *in_list;
+extern char drain[3][3];
+unsigned int first_cum;
+char sides;
+int c_fac;
+int ele_scale;
+struct RB_TREE *draintree;
+
+/* load.c */
+int load_maps(int, int, int);
+
+/* do_astar.c */
+int do_astar(void);
+unsigned int heap_add(int, int, CELL, char);
+
+/* thin.c */
+int thin_streams(void);
+
+/* streams.c */
+int do_accum(double);
+int extract_streams(double, double, int);
+
+/* close.c */
+int close_maps(char *, char *, char *);
+
+#endif /* __LOCAL_PROTO_H__ */
Added: grass-addons/raster/r.stream.extract/main.c
===================================================================
--- grass-addons/raster/r.stream.extract/main.c (rev 0)
+++ grass-addons/raster/r.stream.extract/main.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,267 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.stream.extract
+ * AUTHOR(S): Markus Metz <markus.metz.giswork gmail.com>
+ * PURPOSE: Hydrological analysis
+ * Extracts stream networks from accumulation raster with
+ * given threshold
+ * COPYRIGHT: (C) 1999-2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+char drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
+
+
+int main(int argc, char *argv[])
+{
+ struct
+ {
+ struct Option *ele, *acc, *weight;
+ struct Option *threshold, *d8cut;
+ struct Option *mont_exp;
+ } input;
+ struct
+ {
+ struct Option *stream_rast;
+ struct Option *stream_vect;
+ struct Option *dir_rast;
+ } output;
+ struct GModule *module;
+ int ele_fd, acc_fd, weight_fd;
+ double threshold, d8cut, mont_exp;
+ char *mapset;
+
+ G_gisinit(argv[0]);
+
+ /* Set description */
+ module = G_define_module();
+ module->keywords = _("raster");
+ module->description = _("Stream network extraction");
+
+ input.ele = G_define_standard_option(G_OPT_R_INPUT);
+ input.ele->key = "elevation";
+ input.ele->label = _("Elevation map");
+ input.ele->description = _("Elevation on which entire analysis is based");
+
+ input.acc = G_define_standard_option(G_OPT_R_INPUT);
+ input.acc->key = "accumulation";
+ input.acc->label = _("Accumulation map");
+ input.acc->required = NO;
+ 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");
+ input.threshold->description = _("Must be > 0");
+ input.threshold->required = YES;
+ input.threshold->type = TYPE_DOUBLE;
+
+ input.d8cut = G_define_option();
+ input.d8cut->key = "d8cut";
+ input.d8cut->label = _("Use SFD above this threshold");
+ input.d8cut->description =
+ _("If accumulation is larger than d8cut, SFD is used instead of MFD."
+ " Applies only if no accumulation map is given.");
+ input.d8cut->required = NO;
+ input.d8cut->answer = "infinity";
+ input.d8cut->type = TYPE_DOUBLE;
+
+ input.mont_exp = G_define_option();
+ input.mont_exp->key = "mexp";
+ input.mont_exp->type = TYPE_DOUBLE;
+ input.mont_exp->required = NO;
+ input.mont_exp->answer = "0";
+ input.mont_exp->label =
+ _("Montgomery exponent for slope, disabled with 0");
+ input.mont_exp->description =
+ _("Montgomery: accumulation is multiplied with pow(slope,mexp) and then compared with threshold.");
+
+ output.stream_rast = G_define_standard_option(G_OPT_R_OUTPUT);
+ output.stream_rast->key = "stream_rast";
+ output.stream_rast->description =
+ _("Output raster map with unique stream ids");
+ output.stream_rast->required = NO;
+ output.stream_rast->guisection = _("Output options");
+
+ output.stream_vect = G_define_standard_option(G_OPT_V_OUTPUT);
+ output.stream_vect->key = "stream_vect";
+ output.stream_vect->description =
+ _("Output vector with unique stream ids");
+ output.stream_vect->required = NO;
+ output.stream_vect->guisection = _("Output options");
+
+ output.dir_rast = G_define_standard_option(G_OPT_R_OUTPUT);
+ output.dir_rast->key = "direction";
+ output.dir_rast->description =
+ _("Output raster map with flow direction for streams");
+ output.dir_rast->required = NO;
+ output.dir_rast->guisection = _("Output options");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /***********************/
+ /* check options */
+
+ /***********************/
+
+ /* input maps exist ? */
+ if (!G_find_cell(input.ele->answer, ""))
+ G_fatal_error(_("Raster map <%s> not found"), input.ele->answer);
+
+ if (input.acc->answer) {
+ if (!G_find_cell(input.acc->answer, ""))
+ 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)
+ G_fatal_error(_("Threshold must be > 0 but is %f"), threshold);
+
+ /* d8cut */
+ if (strcmp(input.d8cut->answer, "infinity") == 0) {
+ d8cut = DBL_MAX;
+ }
+ else {
+ d8cut = atof(input.d8cut->answer);
+ if (d8cut < 0)
+ G_fatal_error(_("d8cut must be positive or zero but is %f"),
+ d8cut);
+ }
+
+ /* Montgomery stream initiation */
+ if (input.mont_exp->answer) {
+ mont_exp = atof(input.mont_exp->answer);
+ if (mont_exp < 0)
+ G_fatal_error(_("Montgomery exponent must be positive or zero but is %f"),
+ mont_exp);
+ if (mont_exp > 3)
+ G_warning(_("Montgomery exponent is %f, recommended range is 0.0 - 3.0"),
+ mont_exp);
+ }
+ else
+ mont_exp = 0;
+
+ /* Check for some output map */
+ if ((output.stream_rast->answer == NULL)
+ && (output.stream_vect->answer == NULL)) {
+ G_fatal_error(_("Sorry, you must choose at least one output map."));
+ }
+
+ /*********************/
+ /* preparation */
+
+ /*********************/
+
+ /* open input maps */
+ mapset = G_find_cell2(input.ele->answer, "");
+ ele_fd = G_open_cell_old(input.ele->answer, mapset);
+ if (ele_fd < 0)
+ G_fatal_error(_("could not open input map %s"), input.ele->answer);
+
+ if (input.acc->answer) {
+ mapset = G_find_cell2(input.acc->answer, "");
+ acc_fd = G_open_cell_old(input.acc->answer, mapset);
+ if (acc_fd < 0)
+ G_fatal_error(_("could not open input map %s"),
+ input.acc->answer);
+ }
+ 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();
+ sides = 8; /* not a user option */
+ c_fac = 5; /* not a user option, MFD covergence factor 5 gives best results */
+
+ /* allocate memory */
+ ele = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
+ acc = (DCELL *) G_malloc(nrows * ncols * sizeof(DCELL));
+ stream = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
+ 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)
+ G_fatal_error(_("could not load input map(s)"));
+
+ /********************/
+ /* processing */
+
+ /********************/
+
+ /* sort elevation and get initial stream direction */
+ if (do_astar() < 0)
+ G_fatal_error(_("could not sort elevation map"));
+
+ /* extract streams */
+ if (acc_fd < 0) {
+ if (do_accum(d8cut) < 0)
+ G_fatal_error(_("could not calculate flow accumulation"));
+ if (extract_streams
+ (threshold, mont_exp, (input.weight->answer != NULL)) < 0)
+ G_fatal_error(_("could not extract streams"));
+ }
+ else {
+ if (extract_streams
+ (threshold, mont_exp, (input.weight->answer != NULL)) < 0)
+ G_fatal_error(_("could not extract streams"));
+ }
+
+ G_free(ele);
+ G_free(acc);
+ if (input.weight->answer)
+ G_free(accweight);
+
+ /* thin streams */
+ if (thin_streams() < 0)
+ G_fatal_error(_("could not extract streams"));
+
+ /* write output maps */
+ if (close_maps(output.stream_rast->answer, output.stream_vect->answer,
+ output.dir_rast->answer) < 0)
+ G_fatal_error(_("could not write output maps"));
+
+ G_free(stream);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/raster/r.stream.extract/rbtree.c
===================================================================
--- grass-addons/raster/r.stream.extract/rbtree.c (rev 0)
+++ grass-addons/raster/r.stream.extract/rbtree.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,536 @@
+/*!
+ * \file rbtree.c
+ *
+ * \brief binary search tree
+ *
+ * Generic balanced binary search tree (Red Black Tree) implementation
+ *
+ * (C) 2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public License
+ * (>=v2). Read the file COPYING that comes with GRASS for details.
+ *
+ * \author Original author Julienne Walker 2003, 2008
+ * GRASS implementation Markus Metz, 2009
+ */
+
+/* balanced binary search tree implementation
+ *
+ * this one is a Red Black Tree, the bare version, no parent pointers, no threads
+ * The core code comes from Julienne Walker's tutorials on binary search trees
+ * original license: public domain
+ * http://eternallyconfuzzled.com/tuts/datastructures/jsw_tut_rbtree.aspx
+ * some ideas come from libavl (GPL >= 2)
+ * I could have used some off-the-shelf solution, but that's boring
+ *
+ * Red Black Trees are used to maintain a data structure with
+ * search, insertion and deletion in O(log N) time
+ */
+
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "rbtree.h"
+
+/* internal functions */
+void rbtree_destroy2(struct RB_NODE *);
+struct RB_NODE *rbtree_single(struct RB_NODE *, int);
+struct RB_NODE *rbtree_double(struct RB_NODE *, int);
+void *rbtree_first(struct RB_TRAV *);
+void *rbtree_next(struct RB_TRAV *);
+struct RB_NODE *rbtree_make_node(size_t, void *);
+int is_red(struct RB_NODE *);
+
+
+/* create new tree and initialize
+ * returns pointer to new tree, NULL for memory allocation error
+ */
+struct RB_TREE *rbtree_create(rb_compare_fn *compare, size_t rb_datasize)
+{
+ struct RB_TREE *tree = G_malloc(sizeof(*tree));
+
+ if (tree == NULL) {
+ G_warning("RB tree: Out of memory!");
+ return NULL;
+ }
+
+ assert(compare);
+
+ tree->datasize = rb_datasize;
+ tree->rb_compare = compare;
+ tree->count = 0;
+ tree->root = NULL;
+
+ return tree;
+}
+
+/* add an item to a tree
+ * non-recursive top-down insertion
+ * the algorithm does not allow duplicates and also does not warn about a duplicate
+ * returns 1 on success, 0 on failure
+ */
+int rbtree_insert(struct RB_TREE *tree, void *data)
+{
+ assert(tree && data);
+
+ if (tree->root == NULL) {
+ /* create a new root node for tree */
+ tree->root = rbtree_make_node(tree->datasize, data);
+ if (tree->root == NULL)
+ return 0;
+ }
+ else {
+ struct RB_NODE head = {0}; /* False tree root */
+
+ struct RB_NODE *g, *t; /* Grandparent & parent */
+ struct RB_NODE *p, *q; /* Iterator & parent */
+ int dir = 0, last = 0;
+
+ /* Set up helpers */
+ t = &head;
+ g = p = NULL;
+ q = t->link[1] = tree->root;
+
+ /* Search down the tree */
+ for ( ; ; ) {
+ if (q == NULL) {
+ /* Insert new node at the bottom */
+ p->link[dir] = q = rbtree_make_node(tree->datasize, data);
+ if (q == NULL)
+ return 0;
+ }
+ else if (is_red(q->link[0]) && is_red(q->link[1])) {
+ /* Color flip */
+ q->red = 1;
+ q->link[0]->red = 0;
+ q->link[1]->red = 0;
+ }
+
+ /* Fix red violation */
+ if (is_red(q) && is_red(p)) {
+ int dir2 = t->link[1] == g;
+
+ if (q == p->link[last])
+ t->link[dir2] = rbtree_single(g, !last);
+ else
+ t->link[dir2] = rbtree_double(g, !last);
+ }
+
+ last = dir;
+ dir = tree->rb_compare(q->data, data);
+
+ /* Stop if found. This check also disallows duplicates in the tree */
+ if (dir == 0)
+ break;
+
+ dir = dir < 0;
+
+ /* Move the helpers down */
+ if (g != NULL)
+ t = g;
+
+ g = p, p = q;
+ q = q->link[dir];
+ }
+
+ /* Update root */
+ tree->root = head.link[1];
+ }
+
+ /* Make root black */
+ tree->root->red = 0;
+
+ tree->count++;
+
+ return 1;
+}
+
+/* remove an item from a tree that matches given data
+ * non-recursive top-down removal
+ * returns 1 on successful removal
+ * returns 0 if data item was not found
+ */
+int rbtree_remove(struct RB_TREE *tree, const void *data)
+{
+ struct RB_NODE head = {0}; /* False tree root */
+ struct RB_NODE *q, *p, *g; /* Helpers */
+ struct RB_NODE *f = NULL; /* Found item */
+ int dir = 1, removed = 0;
+
+ assert(tree && data);
+
+ if (tree->root == NULL) {
+ return 0; /* empty tree, nothing to remove */
+ }
+
+ /* Set up helpers */
+ q = &head;
+ g = p = NULL;
+ q->link[1] = tree->root;
+
+ /* Search and push a red down */
+ while (q->link[dir] != NULL) {
+ int last = dir;
+
+ /* Update helpers */
+ g = p, p = q;
+ q = q->link[dir];
+ dir = tree->rb_compare(q->data, data);
+
+ /* Save found node */
+ if (dir == 0)
+ f = q;
+
+ dir = dir < 0;
+
+ /* Push the red node down */
+ if (!is_red(q) && !is_red(q->link[dir])) {
+ if (is_red(q->link[!dir]))
+ p = p->link[last] = rbtree_single(q, dir);
+ else if (!is_red(q->link[!dir])) {
+ struct RB_NODE *s = p->link[!last];
+
+ if (s != NULL) {
+ if (!is_red(s->link[!last]) &&
+ !is_red(s->link[last])) {
+ /* Color flip */
+ p->red = 0;
+ s->red = 1;
+ q->red = 1;
+ }
+ else {
+ int dir2 = g->link[1] == p;
+
+ if (is_red(s->link[last]))
+ g->link[dir2] = rbtree_double(p, last);
+ else if (is_red(s->link[!last]))
+ g->link[dir2] = rbtree_single(p, last);
+
+ /* Ensure correct coloring */
+ q->red = g->link[dir2]->red = 1;
+ g->link[dir2]->link[0]->red = 0;
+ g->link[dir2]->link[1]->red = 0;
+ }
+ }
+ }
+ }
+ }
+
+ /* Replace and remove if found */
+ if (f != NULL) {
+ G_free(f->data);
+ f->data = q->data;
+ p->link[p->link[1] == q] = q->link[q->link[0] == NULL];
+ G_free(q);
+ tree->count--;
+ removed = 1;
+ }
+ else
+ G_debug(2, "RB tree: data not found in search tree");
+
+ /* Update root and make it black */
+ tree->root = head.link[1];
+ if ( tree->root != NULL)
+ tree->root->red = 0;
+
+ return removed;
+}
+
+/* find data item in tree
+ * returns pointer to data item if found else NULL
+ */
+void *rbtree_find(struct RB_TREE *tree, const void *data)
+{
+ struct RB_NODE *curr_node = tree->root;
+ int cmp = 0;
+
+ assert(tree && data);
+
+ while (curr_node != NULL) {
+ cmp = tree->rb_compare(curr_node->data, data);
+ if (cmp == 0)
+ return curr_node->data; /* found */
+ else {
+ curr_node = curr_node->link[cmp < 0];
+ }
+ }
+ return NULL;
+}
+
+/* initialize tree traversal
+ * (re-)sets trav structure
+ * returns 0
+ */
+int rbtree_init_trav(struct RB_TRAV *trav, struct RB_TREE *tree)
+{
+ assert(trav && tree);
+
+ trav->tree = tree;
+ trav->curr_node = tree->root;
+ trav->first = 1;
+ trav->top = 0;
+
+ return 0;
+}
+
+/* traverse the tree in ascending order
+ * useful to get all items in the tree non-recursively
+ * struct RB_TRAV *trav needs to be initialized first
+ * returns pointer to data, NULL when finished
+ */
+void *rbtree_traverse(struct RB_TRAV *trav)
+{
+ assert(trav);
+
+ if (trav->curr_node == NULL) {
+ if (trav->first)
+ G_debug(1, "RB tree: empty tree");
+ else
+ G_debug(1, "RB tree: finished traversing");
+
+ return NULL;
+ }
+
+ if (!trav->first)
+ return rbtree_next(trav);
+ else {
+ trav->first = 0;
+ return rbtree_first(trav);
+ }
+}
+
+/* find start point to traverse the tree in ascending order
+ * useful to get a selection of items in the tree
+ * magnitudes faster than traversing the whole tree
+ * may return first item that's smaller or first item that's larger
+ * struct RB_TRAV *trav needs to be initialized first
+ * returns pointer to data, NULL when finished
+ */
+void *rbtree_traverse_start(struct RB_TRAV *trav, const void *data)
+{
+ int dir = 0;
+
+ assert(trav && data);
+
+ if (trav->curr_node == NULL) {
+ if (trav->first)
+ G_warning("RB tree: empty tree");
+ else
+ G_warning("RB tree: finished traversing");
+
+ return NULL;
+ }
+
+ if (!trav->first)
+ return rbtree_next(trav);
+
+ /* else first time, get start node */
+
+ trav->first = 0;
+ trav->top = 0;
+
+ while (trav->curr_node != NULL) {
+ dir = trav->tree->rb_compare(trav->curr_node->data, data);
+ /* exact match, great! */
+ if (dir == 0)
+ return trav->curr_node->data;
+ else {
+ dir = dir < 0;
+ /* end of branch, also reached if
+ * smallest item is larger than search template or
+ * largest item is smaller than search template */
+ if (trav->curr_node->link[dir] == NULL)
+ return trav->curr_node->data;
+
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[dir];
+ }
+ }
+
+ return NULL; /* should not happen */
+}
+
+/* two functions needed to fully traverse the tree: initialize and continue
+ * useful to get all items in the tree non-recursively
+ * this one here uses a stack
+ * parent pointers or threads would also be possible
+ * but these would need to be added to RB_NODE
+ * -> more memory needed for standard operations
+ */
+
+/* start traversing the tree
+ * returns pointer to smallest data item
+ */
+void *rbtree_first(struct RB_TRAV *trav)
+{
+ /* get smallest item */
+ while (trav->curr_node->link[0] != NULL) {
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[0];
+ }
+
+ return trav->curr_node->data; /* return smallest item */
+}
+
+/* continue traversing the tree in ascending order
+ * returns pointer to data item, NULL when finished
+ */
+void *rbtree_next(struct RB_TRAV *trav)
+{
+ if (trav->curr_node->link[1] != NULL) {
+ /* something on the right side: larger item */
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[1];
+
+ /* go down, find smallest item in this branch */
+ while (trav->curr_node->link[0] != NULL) {
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[0];
+ }
+ }
+ else {
+ /* at smallest item in this branch, go back up */
+ struct RB_NODE *last;
+ do {
+ if (trav->top == 0) {
+ trav->curr_node = NULL;
+ break;
+ }
+ last = trav->curr_node;
+ trav->curr_node = trav->up[--trav->top];
+ } while (last == trav->curr_node->link[1]);
+ }
+
+ if (trav->curr_node != NULL) {
+ return trav->curr_node->data;
+ }
+ else
+ return NULL; /* finished traversing */
+}
+
+/* destroy the tree */
+void rbtree_destroy(struct RB_TREE *tree) {
+ rbtree_destroy2(tree->root);
+ G_free(tree);
+}
+
+void rbtree_destroy2(struct RB_NODE *root)
+{
+ if (root != NULL) {
+ rbtree_destroy2(root->link[0]);
+ rbtree_destroy2(root->link[1]);
+ G_free(root->data);
+ G_free(root);
+ }
+}
+
+/* used for debugging: check for errors in tree structure */
+int rbtree_debug(struct RB_TREE *tree, struct RB_NODE *root)
+{
+ int lh, rh;
+
+ if (root == NULL)
+ return 1;
+ else {
+ struct RB_NODE *ln = root->link[0];
+ struct RB_NODE *rn = root->link[1];
+ int lcmp = 0, rcmp = 0;
+
+ /* Consecutive red links */
+ if (is_red(root)) {
+ if (is_red(ln) || is_red(rn)) {
+ G_warning("Red Black Tree debugging: Red violation");
+ return 0;
+ }
+ }
+
+ lh = rbtree_debug(tree, ln);
+ rh = rbtree_debug(tree, rn);
+
+ if (ln) {
+ lcmp = tree->rb_compare(ln->data, root->data);
+ }
+
+ if (rn) {
+ rcmp = tree->rb_compare(rn->data, root->data);
+ }
+
+ /* Invalid binary search tree:
+ * left node >= parent or right node <= parent */
+ if ((ln != NULL && lcmp > -1)
+ || (rn != NULL && rcmp < 1)) {
+ G_warning("Red Black Tree debugging: Binary tree violation" );
+ return 0;
+ }
+
+ /* Black height mismatch */
+ if (lh != 0 && rh != 0 && lh != rh) {
+ G_warning("Red Black Tree debugging: Black violation");
+ return 0;
+ }
+
+ /* Only count black links */
+ if (lh != 0 && rh != 0)
+ return is_red(root) ? lh : lh + 1;
+ else
+ return 0;
+ }
+}
+
+/*******************************************************
+ * *
+ * internal functions for Red Black Tree maintenance *
+ * *
+ *******************************************************/
+
+/* add a new node to the tree */
+struct RB_NODE *rbtree_make_node(size_t datasize, void *data)
+{
+ struct RB_NODE *new_node = G_malloc(sizeof(*new_node));
+
+ if (new_node == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+
+ new_node->data = G_malloc(datasize);
+ if (new_node->data == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+
+ memcpy(new_node->data, data, datasize);
+ new_node->red = 1; /* 1 is red, 0 is black */
+ new_node->link[0] = NULL;
+ new_node->link[1] = NULL;
+
+ return new_node;
+}
+
+/* check for red violation */
+int is_red(struct RB_NODE *root)
+{
+ if (root)
+ return root->red == 1;
+
+ return 0;
+}
+
+/* single rotation */
+struct RB_NODE *rbtree_single(struct RB_NODE *root, int dir)
+{
+ struct RB_NODE *newroot = root->link[!dir];
+
+ root->link[!dir] = newroot->link[dir];
+ newroot->link[dir] = root;
+
+ root->red = 1;
+ newroot->red = 0;
+
+ return newroot;
+}
+
+/* double rotation */
+struct RB_NODE *rbtree_double(struct RB_NODE *root, int dir)
+{
+ root->link[!dir] = rbtree_single(root->link[!dir], !dir);
+ return rbtree_single(root, dir);
+}
Added: grass-addons/raster/r.stream.extract/rbtree.h
===================================================================
--- grass-addons/raster/r.stream.extract/rbtree.h (rev 0)
+++ grass-addons/raster/r.stream.extract/rbtree.h 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,112 @@
+/*************************************************************
+ * USAGE *
+ *************************************************************
+ *
+ * NOTE: duplicates are not supported
+ *
+ * custom compare function
+ * extern int my_compare_fn(const void *, const void *);
+ * int my_compare_fn(const void *a, const void *b) {
+ * if ((mydatastruct *) a < (mydatastruct *) b)
+ * return -1;
+ * else if ((mydatastruct *) a > (mydatastruct *) b)
+ * return 1;
+ * else if ((mydatastruct *) a == (mydatastruct *) b)
+ * return 0;
+ * }
+ *
+ * create and initialize tree:
+ * struct RB_TREE *mytree = rbtree_create(my_compare_fn, item_size);
+ *
+ * insert items to tree:
+ * struct mydatastruct data = <some data>;
+ * if (rbtree_insert(mytree, &data) == 0)
+ * G_warning("could not insert data");
+ *
+ * find item in tree:
+ * struct mydatastruct data = <some data>;
+ * if (rbtree_find(mytree, &data) == 0)
+ * G_message("data not found");
+ *
+ * delete item from tree:
+ * struct mydatastruct data = <some data>;
+ * if (rbtree_remove(mytree, &data) == 0)
+ * G_warning("could not find data in tree");
+ *
+ * traverse tree (get all items in tree in ascending order):
+ * struct RB_TRAV trav;
+ * rbtree_init_trav(&trav, tree);
+ * while ((data = rbtree_traverse(&trav)) != NULL) {
+ * if (my_compare_fn(data, threshold_data) == 0) break;
+ * <do something with data>;
+ * }
+ *
+ * get a selection of items: all data > data1 and < data2
+ * start in tree where data is last smaller or first larger compared to data1
+ * struct RB_TRAV trav;
+ * rbtree_init_trav(&trav, tree);
+ * data = rbtree_traverse_start(&trav, &data1);
+ * <do something with data>;
+ * while ((data = rbtree_traverse(&trav)) != NULL) {
+ * if (data > data2) break;
+ * <do something with data>;
+ * }
+ *
+ * destroy tree:
+ * rbtree_destroy(mytree);
+ *
+ * debug the whole tree with
+ * rbtree_debug(mytree, mytree->root);
+ *
+ *************************************************************/
+
+#include <stddef.h>
+
+/* maximum RB Tree height */
+#define RBTREE_MAX_HEIGHT 64 /* should be more than enough */
+
+/* routine to compare data items
+ * return -1 if rb_a < rb_b
+ * return 0 if rb_a == rb_b
+ * return 1 if rb_a > rb_b
+ */
+typedef int rb_compare_fn(const void *rb_a, const void *rb_b);
+
+struct RB_NODE
+{
+ unsigned char red; /* 0 = black, 1 = red */
+ void *data; /* any kind of data */
+ struct RB_NODE *link[2]; /* link to children: link[0] for smaller, link[1] for larger */
+};
+
+struct RB_TREE
+{
+ struct RB_NODE *root; /* root node */
+ size_t datasize; /* item size */
+ size_t count; /* number of items in tree. */
+ rb_compare_fn *rb_compare; /* function to compare data */
+};
+
+struct RB_TRAV
+{
+ struct RB_TREE *tree; /* tree being traversed */
+ struct RB_NODE *curr_node; /* current node */
+ struct RB_NODE *up[RBTREE_MAX_HEIGHT]; /* stack of parent nodes */
+ int top; /* index for stack */
+ int first; /* little helper flag */
+};
+
+/* tree functions */
+struct RB_TREE *rbtree_create(rb_compare_fn *, size_t);
+void rbtree_destroy(struct RB_TREE *);
+int rbtree_insert(struct RB_TREE *, void *);
+int rbtree_remove(struct RB_TREE *, const void *);
+void *rbtree_find(struct RB_TREE *, const void *);
+
+/* tree traversal functions */
+int rbtree_init_trav(struct RB_TRAV *, struct RB_TREE *);
+void* rbtree_traverse(struct RB_TRAV *);
+void *rbtree_traverse_start(struct RB_TRAV *, const void *);
+
+/* debug tree from given node downwards */
+int rbtree_debug(struct RB_TREE *, struct RB_NODE *);
Added: grass-addons/raster/r.stream.extract/streams.c
===================================================================
--- grass-addons/raster/r.stream.extract/streams.c (rev 0)
+++ grass-addons/raster/r.stream.extract/streams.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,753 @@
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+/*
+ * compare function for search tree
+ * returns 1 if a > b
+ * returns -1 if a < b
+ * returns 0 if a == b
+ */
+int draindir_compare(const void *itema, const void *itemb)
+{
+ struct ddir *a = (struct ddir *)itema;
+ struct ddir *b = (struct ddir *)itemb;
+
+ if (a->pos > b->pos)
+ return 1;
+ else if (a->pos < b->pos)
+ return -1;
+
+ return 0;
+}
+
+double mfd_pow(double base)
+{
+ int x;
+ double result;
+
+ result = base;
+ if (c_fac == 1)
+ return result;
+
+ for (x = 2; x <= c_fac; x++) {
+ result *= base;
+ }
+ return result;
+}
+
+int continue_stream(CELL is_swale, int r, int c, int r_max, int c_max,
+ unsigned int thisindex, int *stream_no)
+{
+ char aspect;
+ int curr_stream;
+ int r_nbr, c_nbr, ct_dir;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ int stream_node_step = 1000;
+ struct ddir draindir, *founddir;
+
+ G_debug(3, "continue stream");
+
+ /* set drainage direction */
+ aspect = drain[r - r_max + 1][c - c_max + 1];
+
+ /* add to search tree */
+ draindir.pos = thisindex;
+ draindir.dir = aspect;
+ rbtree_insert(draintree, &draindir);
+
+ curr_stream = stream[INDEX(r_max, c_max)];
+ if (curr_stream < 0)
+ curr_stream = 0;
+ /* confluence */
+ if (curr_stream <= 0) {
+ /* no confluence, just continue */
+ G_debug(2, "no confluence, just continue stream");
+ stream[INDEX(r_max, c_max)] = is_swale;
+ }
+ else {
+ G_debug(2, "confluence");
+ /* new confluence */
+ if (stream_node[curr_stream].r != r_max ||
+ stream_node[curr_stream].c != c_max) {
+ G_debug(2, "new confluence");
+ /* set new stream id */
+ curr_stream = stream[INDEX(r_max, c_max)] = ++(*stream_no);
+ /* add stream node */
+ if (*stream_no >= n_alloc_nodes - 1) {
+ n_alloc_nodes += stream_node_step;
+ stream_node =
+ (struct snode *)G_realloc(stream_node,
+ n_alloc_nodes *
+ sizeof(struct snode));
+ }
+ stream_node[*stream_no].r = r_max;
+ stream_node[*stream_no].c = c_max;
+ stream_node[*stream_no].id = *stream_no;
+ stream_node[*stream_no].n_trib = 0;
+ stream_node[*stream_no].n_trib_total = 0;
+ stream_node[*stream_no].n_alloc = 0;
+ stream_node[*stream_no].trib = NULL;
+ stream_node[*stream_no].acc = NULL;
+ n_stream_nodes++;
+
+ /* debug */
+ if (n_stream_nodes != *stream_no)
+ G_warning(_("stream_no %d and n_stream_nodes %d out of sync"),
+ *stream_no, n_stream_nodes);
+
+ /* add all tributaries */
+ G_debug(2, "add all tributaries");
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r_nbr, c_nbr for neighbours */
+ r_nbr = r_max + nextdr[ct_dir];
+ c_nbr = c_max + nextdc[ct_dir];
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
+ c_nbr < ncols) {
+ draindir.pos = INDEX(r_nbr, c_nbr);
+ if ((founddir =
+ rbtree_find(draintree, &draindir)) != NULL) {
+ if (r_nbr + asp_r[(int)founddir->dir] == r_max &&
+ c_nbr + asp_c[(int)founddir->dir] == c_max) {
+
+ /* add tributary to stream node */
+ if (stream_node[curr_stream].n_trib >=
+ stream_node[curr_stream].n_alloc) {
+ size_t new_size;
+
+ stream_node[curr_stream].n_alloc += 2;
+ new_size =
+ stream_node[curr_stream].n_alloc *
+ sizeof(int);
+ stream_node[curr_stream].trib =
+ (int *)G_realloc(stream_node[curr_stream].
+ trib, new_size);
+ new_size =
+ stream_node[curr_stream].n_alloc *
+ sizeof(double);
+ stream_node[curr_stream].acc =
+ (double *)
+ G_realloc(stream_node[curr_stream].acc,
+ new_size);
+ }
+
+ stream_node[curr_stream].
+ trib[stream_node[curr_stream].n_trib] =
+ stream[draindir.pos];
+ stream_node[curr_stream].
+ acc[stream_node[curr_stream].n_trib++] =
+ acc[draindir.pos];
+ }
+ }
+ }
+ }
+
+ /* update stream IDs downstream */
+ G_debug(2, "update stream IDs downstream");
+ r_nbr = r_max;
+ c_nbr = c_max;
+ draindir.pos = INDEX(r_nbr, c_nbr);
+
+ while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
+ if (asp_r[(int)founddir->dir] == 0 &&
+ asp_c[(int)founddir->dir] == 0)
+ G_fatal_error("no valid stream direction");
+ r_nbr = r_nbr + asp_r[(int)founddir->dir];
+ c_nbr = c_nbr + asp_c[(int)founddir->dir];
+ draindir.pos = INDEX(r_nbr, c_nbr);
+ if (stream[INDEX(r_nbr, c_nbr)] <= 0)
+ G_fatal_error("stream id not set");
+ else
+ stream[INDEX(r_nbr, c_nbr)] = curr_stream;
+ }
+ }
+ else {
+ /* stream node already existing here */
+ G_debug(2, "existing confluence");
+ /* add new tributary to stream node */
+ if (stream_node[curr_stream].n_trib >=
+ stream_node[curr_stream].n_alloc) {
+ size_t new_size;
+
+ stream_node[curr_stream].n_alloc += 2;
+ new_size = stream_node[curr_stream].n_alloc * sizeof(int);
+ stream_node[curr_stream].trib =
+ (int *)G_realloc(stream_node[curr_stream].trib, new_size);
+ new_size = stream_node[curr_stream].n_alloc * sizeof(double);
+ stream_node[curr_stream].acc =
+ (double *)G_realloc(stream_node[curr_stream].acc,
+ new_size);
+ }
+
+ stream_node[curr_stream].trib[stream_node[curr_stream].n_trib++] =
+ stream[thisindex];
+ }
+ /* end new confluence */
+
+ G_debug(2, "%d tribs", stream_node[curr_stream].n_trib);
+ if (stream_node[curr_stream].n_trib == 1)
+ G_warning("stream node %d only 1 trib: %d", curr_stream,
+ stream_node[curr_stream].trib[0]);
+ }
+
+ return 1;
+}
+
+/*
+ * extracts streams for threshold
+ */
+int do_accum(double d8cut)
+{
+ int r, c, dr, dc;
+ CELL ele_val, ele_nbr;
+ DCELL value, valued;
+ int count;
+ struct Cell_head window;
+ int mfd_cells, astar_not_set, is_null;
+ double *dist_to_nbr, *weight, sum_weight, max_weight;
+ double dx, dy;
+ int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
+ int is_worked;
+ char aspect;
+ double max_acc, prop;
+ int edge;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ unsigned int thisindex, nindex, workedon, killer;
+
+ G_message(_("calculate flow accumulation..."));
+
+ count = 0;
+
+ /* distances to neighbours */
+ dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
+ weight = (double *)G_malloc(sides * sizeof(double));
+
+ G_get_set_window(&window);
+
+ 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 flag */
+ flag_clear_all(worked);
+
+ for (killer = 1; killer <= n_points; killer++) {
+ G_percent(killer, n_points, 1);
+
+ r = astar_pts[killer].r;
+ c = astar_pts[killer].c;
+ aspect = astar_pts[killer].asp;
+
+ thisindex = INDEX(r, c);
+
+ /* do not distribute flow along edges */
+ if (aspect < 0) {
+ G_debug(3, "edge");
+ continue;
+ }
+
+ if (aspect) {
+ dr = r + asp_r[abs((int)aspect)];
+ dc = c + asp_c[abs((int)aspect)];
+ }
+ else {
+ dr = r;
+ dc = c;
+ }
+
+ r_max = dr;
+ c_max = dc;
+
+ value = acc[thisindex];
+
+ /***************************************/
+ /* get weights for flow distribution */
+
+ /***************************************/
+
+ max_weight = 0;
+ sum_weight = 0;
+ np_side = -1;
+ mfd_cells = 0;
+ astar_not_set = 1;
+ ele_val = ele[thisindex];
+ 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_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) {
+
+ nindex = INDEX(r_nbr, c_nbr);
+
+ is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ if (is_worked == 0) {
+ ele_nbr = ele[nindex];
+ is_null = G_is_c_null_value(&ele_nbr);
+ edge = is_null;
+ if (!is_null && ele_nbr <= ele_val) {
+ if (ele_nbr < ele_val) {
+ weight[ct_dir] =
+ mfd_pow((ele_val -
+ ele_nbr) / dist_to_nbr[ct_dir]);
+ }
+ if (ele_nbr == ele_val) {
+ weight[ct_dir] =
+ mfd_pow(0.5 / dist_to_nbr[ct_dir]);
+ }
+ 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) {
+ G_debug(3, "edge");
+ 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
+ */
+
+ /************************************/
+ /* distribute and accumulate flow */
+
+ /************************************/
+
+ /* 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;
+ }
+
+ /* use SFD (D8) if d8cut threshold exceeded */
+ if (fabs(value) > d8cut)
+ mfd_cells = 0;
+
+ 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];
+
+ nindex = INDEX(r_nbr, c_nbr);
+
+ valued = acc[nindex];
+ valued += value * weight[ct_dir];
+ acc[nindex] = valued;
+ }
+ else if (ct_dir == np_side) {
+ /* check for consistency with A * path */
+ workedon++;
+ }
+ }
+ }
+ if (fabs(prop - 1.0) > 5E-6f) {
+ G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
+ prop);
+ }
+ }
+ /* get out of depression in SFD mode */
+ else {
+ nindex = INDEX(dr, dc);
+ valued = acc[INDEX(dr, dc)];
+ valued += value;
+ acc[INDEX(dr, dc)] = valued;
+ }
+
+ FLAG_SET(worked, r, c);
+ }
+
+ G_free(dist_to_nbr);
+ G_free(weight);
+
+ return 1;
+}
+
+/*
+ * extracts streams for threshold, accumulation is provided
+ */
+int extract_streams(double threshold, double mont_exp, int use_weight)
+{
+ int r, c, dr, dc;
+ CELL is_swale, ele_val, ele_nbr;
+ DCELL value, valued;
+ struct Cell_head window;
+ int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
+ double *dist_to_nbr;
+ double dx, dy;
+ int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side, max_side;
+ int is_worked;
+ char aspect;
+ double max_acc;
+ int edge, flat;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ /* sides */
+ /*
+ * | 7 | 1 | 4 |
+ * | 2 | | 3 |
+ * | 5 | 0 | 6 |
+ */
+ unsigned int thisindex, nindex, workedon, killer;
+ int stream_no = 0, stream_node_step = 1000;
+ double slope, diag;
+
+ G_message(_("extract streams..."));
+
+ /* init BST for drainage direction */
+ draintree = rbtree_create(draindir_compare, sizeof(struct ddir));
+
+ /* init stream nodes */
+ n_alloc_nodes = stream_node_step;
+ stream_node =
+ (struct snode *)G_malloc(n_alloc_nodes * sizeof(struct snode));
+ n_stream_nodes = 0;
+
+ /* init outlet nodes */
+ n_alloc_outlets = stream_node_step;
+ outlets =
+ (struct point *)G_malloc(n_alloc_outlets * sizeof(struct point));
+ n_outlets = 0;
+
+ /* distances to neighbours */
+ dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
+
+ G_get_set_window(&window);
+
+ 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);
+ }
+
+ diag = sqrt(2);
+
+ /* reset worked flag */
+ flag_clear_all(worked);
+
+ workedon = 0;
+
+ for (killer = 1; killer <= n_points; killer++) {
+ G_percent(killer, n_points, 1);
+
+ r = astar_pts[killer].r;
+ c = astar_pts[killer].c;
+ aspect = astar_pts[killer].asp;
+
+ thisindex = INDEX(r, c);
+
+ /* do not distribute flow along edges */
+ if (aspect < 0) {
+ G_debug(3, "edge");
+ is_swale = stream[thisindex];
+ if (is_swale) {
+ G_debug(2, "edge outlet");
+ /* add outlet point */
+ if (n_outlets >= n_alloc_outlets) {
+ n_alloc_outlets += stream_node_step;
+ outlets =
+ (struct point *)G_realloc(outlets,
+ n_alloc_outlets *
+ sizeof(struct point));
+ }
+ outlets[n_outlets].r = r;
+ outlets[n_outlets].c = c;
+ n_outlets++;
+ }
+ continue;
+ }
+
+ if (aspect) {
+ dr = r + asp_r[abs((int)aspect)];
+ dc = c + asp_c[abs((int)aspect)];
+ }
+ else {
+ dr = r;
+ dc = c;
+ }
+
+ r_max = dr;
+ c_max = dc;
+
+ value = acc[thisindex];
+
+ /************************************/
+ /* find main drainage direction */
+
+ /************************************/
+
+ max_acc = -1;
+ max_side = np_side = -1;
+ mfd_cells = 0;
+ stream_cells = 0;
+ swale_cells = 0;
+ astar_not_set = 1;
+ ele_val = ele[thisindex];
+ is_null = 0;
+ edge = 0;
+ flat = 1;
+ /* this loop is needed to get the sum of weights */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r_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) {
+
+ if (dr == r_nbr && dc == c_nbr)
+ np_side = ct_dir;
+
+ nindex = INDEX(r_nbr, c_nbr);
+
+ /* check for swale cells */
+ is_swale = stream[nindex];
+ if (is_swale > 0)
+ swale_cells++;
+ /* check for stream cells */
+ valued = fabs(acc[nindex]);
+ if (valued >= threshold)
+ stream_cells++;
+
+ is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ if (is_worked == 0) {
+ ele_nbr = ele[nindex];
+ if (ele_nbr != ele_val)
+ flat = 0;
+ edge = G_is_c_null_value(&ele_nbr);
+ if (!edge && ele_nbr <= ele_val) {
+
+ mfd_cells++;
+
+ /* set main drainage direction */
+ if (valued >= max_acc) {
+ max_acc = valued;
+ r_max = r_nbr;
+ c_max = c_nbr;
+ max_side = ct_dir;
+ }
+
+ if (dr == r_nbr && dc == c_nbr) {
+ astar_not_set = 0;
+ }
+ }
+ }
+ else if (ct_dir == np_side) {
+ /* check for consistency with A * path */
+ workedon++;
+ }
+ }
+ else
+ edge = 1;
+ if (edge)
+ break;
+ }
+
+ is_swale = stream[thisindex];
+
+ /* do not distribute flow along edges, this causes artifacts */
+ if (edge) {
+ G_debug(3, "edge");
+ if (is_swale) {
+ G_debug(2, "edge outlet");
+ /* add outlet point */
+ if (n_outlets >= n_alloc_outlets) {
+ n_alloc_outlets += stream_node_step;
+ outlets =
+ (struct point *)G_realloc(outlets,
+ n_alloc_outlets *
+ sizeof(struct point));
+ }
+ outlets[n_outlets].r = r;
+ outlets[n_outlets].c = c;
+ n_outlets++;
+ }
+ 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++;
+ /* get main drainage direction */
+ nindex = INDEX(dr, dc);
+ if (fabs(acc[nindex]) >= max_acc) {
+ max_acc = fabs(acc[nindex]);
+ r_max = dr;
+ c_max = dc;
+ max_side = np_side;
+ }
+ }
+ else if (mfd_cells == 0) {
+ flat = 0;
+ max_side = np_side;
+ }
+
+ is_swale = stream[thisindex];
+
+ /**********************/
+ /* 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) {
+ if (r_max == r && c_max == c)
+ G_warning
+ ("Can't use Montgomery's method, no stream direction found");
+ else {
+
+ ele_nbr = ele[INDEX(r_max, c_max)];
+
+ slope = (double)(ele_val - ele_nbr) / ele_scale;
+
+ if (max_side > 3)
+ slope /= diag;
+
+ value *= pow(fabs(slope), mont_exp);
+
+ }
+ }
+
+ if (is_swale < 1 && flat == 0 && fabs(value) >= threshold &&
+ stream_cells < 4 && swale_cells < 1) {
+ G_debug(2, "start new stream");
+ is_swale = stream[thisindex] = ++stream_no;
+ /* add stream node */
+ if (stream_no >= n_alloc_nodes - 1) {
+ n_alloc_nodes += stream_node_step;
+ stream_node =
+ (struct snode *)G_realloc(stream_node,
+ n_alloc_nodes *
+ sizeof(struct snode));
+ }
+ stream_node[stream_no].r = r;
+ stream_node[stream_no].c = c;
+ stream_node[stream_no].id = stream_no;
+ stream_node[stream_no].n_trib = 0;
+ stream_node[stream_no].n_trib_total = 0;
+ stream_node[stream_no].n_alloc = 0;
+ stream_node[stream_no].trib = NULL;
+ stream_node[stream_no].acc = NULL;
+ n_stream_nodes++;
+
+ /* debug */
+ if (n_stream_nodes != stream_no)
+ G_warning(_("stream_no %d and n_stream_nodes %d out of sync"),
+ stream_no, n_stream_nodes);
+ }
+
+ /*********************/
+ /* continue stream */
+
+ /*********************/
+
+ /* can't continue stream, add outlet point */
+ if (is_swale > 0 && r_max == r && c_max == c) {
+ G_debug(2, "can't continue stream at %d %d", r, c);
+
+ if (n_outlets >= n_alloc_outlets) {
+ n_alloc_outlets += stream_node_step;
+ outlets =
+ (struct point *)G_malloc(n_alloc_outlets *
+ sizeof(struct point));
+ }
+ outlets[n_outlets].r = r;
+ outlets[n_outlets].c = c;
+ n_outlets++;
+ }
+ else if (is_swale > 0) {
+ continue_stream(is_swale, r, c, r_max, c_max, thisindex,
+ &stream_no);
+ }
+
+ FLAG_SET(worked, r, c);
+ }
+ if (workedon)
+ G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
+ workedon, n_points);
+
+ flag_destroy(worked);
+
+ G_free(dist_to_nbr);
+
+ G_debug(1, "%d outlets", n_outlets);
+ G_debug(1, "%d nodes", n_stream_nodes);
+ G_debug(1, "%d streams", stream_no);
+
+ return 1;
+}
Added: grass-addons/raster/r.stream.extract/thin.c
===================================================================
--- grass-addons/raster/r.stream.extract/thin.c (rev 0)
+++ grass-addons/raster/r.stream.extract/thin.c 2009-09-22 14:42:16 UTC (rev 39281)
@@ -0,0 +1,175 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+int thin_seg(int stream_id)
+{
+ int thinned = 0;
+ int r, c, r_nbr, c_nbr, last_r, last_c;
+ unsigned int thisindex, lastindex;
+ struct ddir draindir, *founddir;
+ int curr_stream;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ r = stream_node[stream_id].r;
+ c = stream_node[stream_id].c;
+
+ thisindex = INDEX(r, c);
+
+ curr_stream = stream[thisindex];
+ if (curr_stream != stream_id)
+ G_fatal_error(_("stream node and stream not identical: stream id %d, stream node id %d, stream %d"),
+ stream_id, stream_node[stream_id].id, curr_stream);
+
+ draindir.pos = thisindex;
+ if ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
+ /* get downstream point */
+ last_r = r + asp_r[(int)founddir->dir];
+ last_c = c + asp_c[(int)founddir->dir];
+ curr_stream = stream[INDEX(last_r, last_c)];
+
+ if (curr_stream != stream_id)
+ return thinned;
+
+ /* get next downstream point */
+ draindir.pos = INDEX(last_r, last_c);
+ while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
+ r_nbr = last_r + asp_r[(int)founddir->dir];
+ c_nbr = last_c + asp_c[(int)founddir->dir];
+
+ if (r_nbr == last_r && c_nbr == last_c)
+ return thinned;
+ if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
+ return thinned;
+ if ((curr_stream = stream[INDEX(r_nbr, c_nbr)]) != stream_id)
+ return thinned;
+ if (abs(r_nbr - r) < 2 && abs(c_nbr - c) < 2) {
+ /* eliminate last point */
+ lastindex = INDEX(last_r, last_c);
+ stream[lastindex] = 0;
+ draindir.pos = lastindex;
+ rbtree_remove(draintree, &draindir);
+ /* update start point */
+ draindir.pos = thisindex;
+ founddir = rbtree_find(draintree, &draindir);
+ founddir->dir = drain[r - r_nbr + 1][c - c_nbr + 1];
+ last_r = r_nbr;
+ last_c = c_nbr;
+ draindir.pos = INDEX(last_r, last_c);
+
+ thinned = 1;
+ }
+ else {
+ /* nothing to eliminate, continue from last point */
+ r = last_r;
+ c = last_c;
+ last_r = r_nbr;
+ last_c = c_nbr;
+ thisindex = INDEX(r, c);
+ draindir.pos = INDEX(last_r, last_c);
+ }
+ }
+ }
+
+ return thinned;
+}
+
+int thin_streams(void)
+{
+ int i, j, r, c, done;
+ int stream_id, next_node;
+ unsigned int thisindex;
+ struct sstack
+ {
+ int stream_id;
+ int next_trib;
+ } *nodestack;
+ int top = 0, stack_step = 1000;
+ int n_trib_total;
+
+ G_message(_("thin stream segments"));
+
+ nodestack = (struct sstack *)G_malloc(stack_step * sizeof(struct sstack));
+
+ for (i = 0; i < n_outlets; i++) {
+ G_percent(i, n_outlets, 2);
+ r = outlets[i].r;
+ c = outlets[i].c;
+ thisindex = INDEX(r, c);
+ stream_id = stream[thisindex];
+
+ /* add root node to stack */
+ G_debug(2, "add root node");
+ top = 0;
+ nodestack[top].stream_id = stream_id;
+ nodestack[top].next_trib = 0;
+
+ /* depth first post order traversal */
+ G_debug(2, "traverse");
+ while (top >= 0) {
+
+ done = 1;
+ stream_id = nodestack[top].stream_id;
+ G_debug(3, "stream_id %d, top %d", stream_id, top);
+ if (nodestack[top].next_trib < stream_node[stream_id].n_trib) {
+ /* add to stack */
+ G_debug(3, "get next node");
+ next_node =
+ stream_node[stream_id].trib[nodestack[top].next_trib];
+ G_debug(3, "add to stack: next %d, trib %d, n trib %d",
+ next_node, nodestack[top].next_trib,
+ stream_node[stream_id].n_trib);
+ nodestack[top].next_trib++;
+ top++;
+ if (top >= stack_step) {
+ /* need more space */
+ stack_step += 1000;
+ nodestack =
+ (struct sstack *)G_realloc(nodestack,
+ stack_step *
+ sizeof(struct sstack));
+ }
+
+ nodestack[top].next_trib = 0;
+ nodestack[top].stream_id = next_node;
+ done = 0;
+ G_debug(3, "go further down");
+ }
+ if (done) {
+ /* thin stream segment */
+ G_debug(3, "thin stream segment %d", stream_id);
+
+ if (thin_seg(stream_id) == 0)
+ G_debug(3, "segment %d not thinned", stream_id);
+ else
+ G_debug(3, "segment %d thinned", stream_id);
+
+ top--;
+ /* count tributaries */
+ if (top >= 0) {
+ n_trib_total = 0;
+ stream_id = nodestack[top].stream_id;
+ for (j = 0; j < stream_node[stream_id].n_trib; j++) {
+ /* intermediate */
+ if (stream_node[stream_node[stream_id].trib[j]].
+ n_trib > 0)
+ n_trib_total +=
+ stream_node[stream_node[stream_id].trib[j]].
+ n_trib_total;
+ /* start */
+ else
+ n_trib_total++;
+ }
+ stream_node[stream_id].n_trib_total = n_trib_total;
+ }
+ }
+ }
+ }
+ G_percent(n_outlets, n_outlets, 1); /* finish it */
+
+ G_free(nodestack);
+
+ return 1;
+}
More information about the grass-commit
mailing list