[GRASS-SVN] r44615 - grass-addons/grass7/raster/r.stream.extract
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Dec 16 09:01:19 EST 2010
Author: mmetz
Date: 2010-12-16 06:01:19 -0800 (Thu, 16 Dec 2010)
New Revision: 44615
Added:
grass-addons/grass7/raster/r.stream.extract/bseg.c
grass-addons/grass7/raster/r.stream.extract/cseg.c
grass-addons/grass7/raster/r.stream.extract/dseg.c
grass-addons/grass7/raster/r.stream.extract/init_search.c
grass-addons/grass7/raster/r.stream.extract/r.stream.extract.html
grass-addons/grass7/raster/r.stream.extract/seg.c
grass-addons/grass7/raster/r.stream.extract/seg.h
Removed:
grass-addons/grass7/raster/r.stream.extract/description.html
grass-addons/grass7/raster/r.stream.extract/flag.c
grass-addons/grass7/raster/r.stream.extract/rbtree.c
grass-addons/grass7/raster/r.stream.extract/rbtree.h
Modified:
grass-addons/grass7/raster/r.stream.extract/Makefile
grass-addons/grass7/raster/r.stream.extract/close.c
grass-addons/grass7/raster/r.stream.extract/del_streams.c
grass-addons/grass7/raster/r.stream.extract/do_astar.c
grass-addons/grass7/raster/r.stream.extract/flag.h
grass-addons/grass7/raster/r.stream.extract/load.c
grass-addons/grass7/raster/r.stream.extract/local_proto.h
grass-addons/grass7/raster/r.stream.extract/main.c
grass-addons/grass7/raster/r.stream.extract/streams.c
grass-addons/grass7/raster/r.stream.extract/thin.c
Log:
r.stream.extract for 7
Modified: grass-addons/grass7/raster/r.stream.extract/Makefile
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/Makefile 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/Makefile 2010-12-16 14:01:19 UTC (rev 44615)
@@ -2,8 +2,8 @@
PGM = r.stream.extract
-LIBES = $(VECTLIB) $(DBMILIB) $(GISLIB)
-DEPENDENCIES = $(VECTDEP) $(DBMIDEP) $(GISDEP)
+LIBES = $(SEGMENTLIB) $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(SEGMENTDEP) $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
EXTRA_INC = $(VECT_INC)
EXTRA_CFLAGS = $(VECT_CFLAGS)
Added: grass-addons/grass7/raster/r.stream.extract/bseg.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/bseg.c (rev 0)
+++ grass-addons/grass7/raster/r.stream.extract/bseg.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -0,0 +1,159 @@
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/glocale.h>
+#include "seg.h"
+
+int bseg_open(BSEG *bseg, int srows, int scols, int nsegs_in_memory)
+{
+ char *filename;
+ int errflag;
+ int fd;
+
+ bseg->filename = NULL;
+ bseg->fd = -1;
+ bseg->name = NULL;
+ bseg->mapset = NULL;
+
+ filename = G_tempfile();
+ if (-1 == (fd = creat(filename, 0666))) {
+ G_warning(_("bseg_open(): unable to create segment file"));
+ return -2;
+ }
+ if (0 > (errflag = segment_format(fd, Rast_window_rows(),
+ Rast_window_cols(), srows, scols,
+ sizeof(char)))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning(_("bseg_open(): could not write segment file"));
+ return -1;
+ }
+ else {
+ G_warning(_("bseg_open(): illegal configuration parameter(s)"));
+ return -3;
+ }
+ }
+ close(fd);
+ if (-1 == (fd = open(filename, 2))) {
+ unlink(filename);
+ G_warning(_("bseg_open(): unable to re-open segment file"));
+ return -4;
+ }
+ if (0 > (errflag = segment_init(&(bseg->seg), fd, nsegs_in_memory))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning(_("bseg_open(): could not read segment file"));
+ return -5;
+ }
+ else {
+ G_warning(_("bseg_open(): out of memory"));
+ return -6;
+ }
+ }
+ bseg->filename = filename;
+ bseg->fd = fd;
+ return 0;
+}
+
+int bseg_close(BSEG *bseg)
+{
+ segment_release(&(bseg->seg));
+ close(bseg->fd);
+ unlink(bseg->filename);
+ if (bseg->name) {
+ G_free(bseg->name);
+ bseg->name = NULL;
+ }
+ if (bseg->mapset) {
+ G_free(bseg->mapset);
+ bseg->mapset = NULL;
+ }
+ return 0;
+}
+
+int bseg_put(BSEG *bseg, char *value, int row, int col)
+{
+ if (segment_put(&(bseg->seg), value, row, col) < 0) {
+ G_warning(_("bseg_put(): could not write segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int bseg_put_row(BSEG *bseg, char *value, int row)
+{
+ if (segment_put_row(&(bseg->seg), value, row) < 0) {
+ G_warning(_("bseg_put_row(): could not write segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int bseg_get(BSEG *bseg, char *value, int row, int col)
+{
+ if (segment_get(&(bseg->seg), value, row, col) < 0) {
+ G_warning(_("bseg_get(): could not read segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+
+int bseg_read_raster(BSEG *bseg, char *map_name, char *mapset)
+{
+ int row, nrows;
+ int col, ncols;
+ int map_fd;
+ CELL *buffer;
+ char cbuf;
+
+ bseg->name = NULL;
+ bseg->mapset = NULL;
+
+ map_fd = Rast_open_old(map_name, mapset);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ buffer = Rast_allocate_c_buf();
+ for (row = 0; row < nrows; row++) {
+ Rast_get_c_row(map_fd, buffer, row);
+ for (col = ncols; col >= 0; col--) {
+ cbuf = (char) buffer[col];
+ bseg_put(bseg, &cbuf, row, col);
+ }
+ }
+
+ Rast_close(map_fd);
+ G_free(buffer);
+
+ bseg->name = G_store(map_name);
+ bseg->mapset = G_store(mapset);
+
+ return 0;
+}
+
+int bseg_write_raster(BSEG *bseg, char *map_name)
+{
+ int map_fd;
+ int row, nrows;
+ int col, ncols;
+ CELL *buffer;
+ char value;
+
+ map_fd = Rast_open_c_new(map_name);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ buffer = Rast_allocate_c_buf();
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 1);
+ for (col = 0; col < ncols; col++) {
+ bseg_get(bseg, &value, row, col);
+ buffer[col] = value;
+ }
+ Rast_put_row(map_fd, buffer, CELL_TYPE);
+ }
+ G_percent(row, nrows, 1); /* finish it */
+ G_free(buffer);
+ Rast_close(map_fd);
+ return 0;
+}
Property changes on: grass-addons/grass7/raster/r.stream.extract/bseg.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/raster/r.stream.extract/close.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/close.c 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/close.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,14 +1,15 @@
-#include <grass/gis.h>
+#include <grass/raster.h>
#include <grass/glocale.h>
+#include <grass/vector.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;
+ CELL stream_id, stream_nbr;
+ char aspect;
+ int next_node;
struct sstack
{
int stream_id;
@@ -17,7 +18,6 @@
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;
@@ -53,8 +53,7 @@
G_percent(i, n_outlets, 2);
r = outlets[i].r;
c = outlets[i].c;
- thisindex = INDEX(r, c);
- stream_id = stream[thisindex];
+ cseg_get(&stream, &stream_id, r, c);
if (!stream_id)
continue;
@@ -112,8 +111,11 @@
r_nbr = stream_node[stream_id].r;
c_nbr = stream_node[stream_id].c;
- draindir.pos = INDEX(r_nbr, c_nbr);
+ cseg_get(&stream, &stream_nbr, r_nbr, c_nbr);
+ if (stream_nbr <= 0)
+ G_fatal_error("stream id %d not set, top is %d, parent is %d", stream_id, top, nodestack[top - 1].stream_id);
+
Vect_cat_set(Cats, 1, stream_id);
if (stream_node[stream_id].n_trib == 0)
Vect_cat_set(Cats, 2, 0);
@@ -125,23 +127,22 @@
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];
+ bseg_get(&asp, &aspect, r_nbr, c_nbr);
+ while (aspect > 0) {
+ r_nbr = r_nbr + asp_r[(int)aspect];
+ c_nbr = c_nbr + asp_c[(int)aspect];
+
+ cseg_get(&stream, &stream_nbr, r_nbr, c_nbr);
+ if (stream_nbr <= 0)
+ G_fatal_error("stream id not set while tracing");
- 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);
+ Vect_append_point(Points, west_offset + c_nbr * ew_res,
+ north_offset - r_nbr * ns_res, 0);
+ if (stream_nbr != stream_id) {
+ /* first point of parent stream */
break;
}
- draindir.pos = INDEX(r_nbr, c_nbr);
- if (stream[INDEX(r_nbr, c_nbr)] <= 0)
- G_fatal_error(_("BUG: stream id not set"));
-
- Vect_append_point(Points, west_offset + c_nbr * ew_res,
- north_offset - r_nbr * ns_res, 0);
+ bseg_get(&asp, &aspect, r_nbr, c_nbr);
}
Vect_write_line(&Out, GV_LINE, Points, Cats);
@@ -232,8 +233,10 @@
{
int stream_fd, dir_fd, r, c, i;
CELL *cell_buf1, *cell_buf2;
- unsigned int thisindex;
struct History history;
+ char flag_value;
+ CELL stream_id;
+ char aspect;
/* cheating... */
stream_fd = dir_fd = -1;
@@ -244,54 +247,56 @@
/* write requested output rasters */
if (stream_rast) {
- stream_fd = G_open_raster_new(stream_rast, CELL_TYPE);
- cell_buf1 = G_allocate_cell_buf();
+ stream_fd = Rast_open_new(stream_rast, CELL_TYPE);
+ cell_buf1 = Rast_allocate_c_buf();
}
if (dir_rast) {
- dir_fd = G_open_raster_new(dir_rast, CELL_TYPE);
- cell_buf2 = G_allocate_cell_buf();
+ dir_fd = Rast_open_new(dir_rast, CELL_TYPE);
+ cell_buf2 = Rast_allocate_c_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 */
+ Rast_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 */
+ Rast_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 (stream_rast) {
+ cseg_get(&stream, &stream_id, r, c);
+ if (stream_id)
+ cell_buf1[c] = stream_id;
}
if (dir_rast) {
- if (!G_is_c_null_value(&ele[thisindex])) {
- cell_buf2[c] = asp[thisindex];
+ bseg_get(&bitflags, &flag_value, r, c);
+ if (!FLAG_GET(flag_value, NULLFLAG)) {
+ bseg_get(&asp, &aspect, r, c);
+ cell_buf2[c] = aspect;
}
}
}
if (stream_rast)
- G_put_raster_row(stream_fd, cell_buf1, CELL_TYPE);
+ Rast_put_row(stream_fd, cell_buf1, CELL_TYPE);
if (dir_rast)
- G_put_raster_row(dir_fd, cell_buf2, CELL_TYPE);
+ Rast_put_row(dir_fd, cell_buf2, CELL_TYPE);
}
G_percent(nrows, nrows, 2); /* finish it */
if (stream_rast) {
- G_close_cell(stream_fd);
+ Rast_close(stream_fd);
G_free(cell_buf1);
- G_short_history(stream_rast, "raster", &history);
- G_command_history(&history);
- G_write_history(stream_rast, &history);
+ Rast_short_history(stream_rast, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(stream_rast, &history);
}
if (dir_rast) {
- G_close_cell(dir_fd);
+ Rast_close(dir_fd);
G_free(cell_buf2);
- G_short_history(dir_rast, "raster", &history);
- G_command_history(&history);
- G_write_history(dir_rast, &history);
+ Rast_short_history(dir_rast, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(dir_rast, &history);
}
/* close stream vector */
@@ -301,14 +306,12 @@
}
/* 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);
Property changes on: grass-addons/grass7/raster/r.stream.extract/close.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.stream.extract/cseg.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/cseg.c (rev 0)
+++ grass-addons/grass7/raster/r.stream.extract/cseg.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -0,0 +1,154 @@
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/glocale.h>
+#include "seg.h"
+
+int cseg_open(CSEG *cseg, int srows, int scols, int nsegs_in_memory)
+{
+ char *filename;
+ int errflag;
+ int fd;
+
+ cseg->filename = NULL;
+ cseg->fd = -1;
+ cseg->name = NULL;
+ cseg->mapset = NULL;
+
+ filename = G_tempfile();
+ if (-1 == (fd = creat(filename, 0666))) {
+ G_warning(_("cseg_open(): unable to create segment file"));
+ return -2;
+ }
+ if (0 >
+ (errflag =
+ segment_format(fd, Rast_window_rows(), Rast_window_cols(), srows, scols,
+ sizeof(CELL)))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning(_("cseg_open(): could not write segment file"));
+ return -1;
+ }
+ else {
+ G_warning(_("cseg_open(): illegal configuration parameter(s)"));
+ return -3;
+ }
+ }
+ close(fd);
+ if (-1 == (fd = open(filename, 2))) {
+ unlink(filename);
+ G_warning(_("cseg_open(): unable to re-open segment file"));
+ return -4;
+ }
+ if (0 > (errflag = segment_init(&(cseg->seg), fd, nsegs_in_memory))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning(_("cseg_open(): could not read segment file"));
+ return -5;
+ }
+ else {
+ G_warning(_("cseg_open(): out of memory"));
+ return -6;
+ }
+ }
+ cseg->filename = filename;
+ cseg->fd = fd;
+ return 0;
+}
+
+int cseg_close(CSEG *cseg)
+{
+ segment_release(&(cseg->seg));
+ close(cseg->fd);
+ unlink(cseg->filename);
+ if (cseg->name) {
+ G_free(cseg->name);
+ cseg->name = NULL;
+ }
+ if (cseg->mapset) {
+ G_free(cseg->mapset);
+ cseg->mapset = NULL;
+ }
+ return 0;
+}
+
+int cseg_put(CSEG *cseg, CELL *value, int row, int col)
+{
+ if (segment_put(&(cseg->seg), value, row, col) < 0) {
+ G_warning(_("cseg_put(): could not write segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int cseg_put_row(CSEG *cseg, CELL *value, int row)
+{
+ if (segment_put_row(&(cseg->seg), value, row) < 0) {
+ G_warning(_("cseg_put_row(): could not write segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int cseg_get(CSEG *cseg, CELL *value, int row, int col)
+{
+ if (segment_get(&(cseg->seg), value, row, col) < 0) {
+ G_warning(_("cseg_get(): could not read segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int cseg_read_raster(CSEG *cseg, char *map_name, char *mapset)
+{
+ int row, nrows;
+ int map_fd;
+ CELL *buffer;
+
+ cseg->name = NULL;
+ cseg->mapset = NULL;
+
+ map_fd = Rast_open_old(map_name, mapset);
+ nrows = Rast_window_rows();
+ buffer = Rast_allocate_c_buf();
+ for (row = 0; row < nrows; row++) {
+ Rast_get_c_row(map_fd, buffer, row);
+ if (segment_put_row(&(cseg->seg), buffer, row) < 0) {
+ G_free(buffer);
+ Rast_close(map_fd);
+ G_warning(_("cseg_read_cell(): unable to segment put row for <%s> in <%s>"),
+ map_name, mapset);
+ return (-1);
+ }
+ }
+
+ Rast_close(map_fd);
+ G_free(buffer);
+
+ cseg->name = G_store(map_name);
+ cseg->mapset = G_store(mapset);
+
+ return 0;
+}
+
+int cseg_write_raster(CSEG *cseg, char *map_name)
+{
+ int map_fd;
+ int row, nrows;
+ CELL *buffer;
+
+ map_fd = Rast_open_c_new(map_name);
+ nrows = Rast_window_rows();
+ buffer = Rast_allocate_c_buf();
+ segment_flush(&(cseg->seg));
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 1);
+ segment_get_row(&(cseg->seg), buffer, row);
+ Rast_put_row(map_fd, buffer, CELL_TYPE);
+ }
+ G_percent(row, nrows, 1); /* finish it */
+ G_free(buffer);
+ Rast_close(map_fd);
+ return 0;
+}
Property changes on: grass-addons/grass7/raster/r.stream.extract/cseg.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/raster/r.stream.extract/del_streams.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/del_streams.c 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/del_streams.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,101 +1,29 @@
#include <stdlib.h>
#include <math.h>
-#include <grass/gis.h>
+#include <grass/raster.h>
#include <grass/glocale.h>
#include "local_proto.h"
-/* delete short stream segments according to threshold */
-int del_streams(int min_length)
-{
- int i;
- int n_deleted = 0;
- unsigned int thisindex, next_stream_pos = -1;
- int curr_stream, stream_id, other_trib, tmp_trib;
- int slength;
-
- G_message(_("Delete stream segments shorter than %d cells..."), min_length);
-
- /* go through all nodes */
- for (i = 1; i <= n_stream_nodes; i++) {
- G_percent(i, n_stream_nodes, 2);
-
- /* not a stream head */
- if (stream_node[i].n_trib > 0)
- continue;
-
- /* already deleted */
- thisindex = INDEX(stream_node[i].r, stream_node[i].c);
- if (stream[thisindex] == 0)
- continue;
-
- /* get length counted as n cells */
- if ((slength = seg_length(i, &next_stream_pos)) >= min_length)
- continue;
-
- stream_id = i;
-
- /* check n sibling tributaries */
- if ((curr_stream = stream[next_stream_pos]) != stream_id) {
- /* only one sibling tributary */
- if (stream_node[curr_stream].n_trib == 2) {
- if (stream_node[curr_stream].trib[0] != stream_id)
- other_trib = stream_node[curr_stream].trib[0];
- else
- other_trib = stream_node[curr_stream].trib[1];
-
- /* other trib is also stream head */
- if (stream_node[other_trib].n_trib == 0) {
- /* use shorter one */
- if (seg_length(other_trib, NULL) < slength) {
- tmp_trib = stream_id;
- stream_id = other_trib;
- other_trib = tmp_trib;
- }
- }
- del_stream_seg(stream_id);
- n_deleted++;
-
- /* update downstream IDs */
- update_stream_id(curr_stream, other_trib);
- }
- /* more than one other tributary */
- else {
- del_stream_seg(stream_id);
- n_deleted++;
- }
- }
- /* stream head is also outlet */
- else {
- del_stream_seg(stream_id);
- n_deleted++;
- }
- }
-
- G_verbose_message(_("%d stream segments deleted"), n_deleted);
-
- return n_deleted;
-}
-
/* get stream segment length */
-int seg_length(int stream_id, unsigned int *new_stream_pos)
+int seg_length(CELL stream_id, CELL *next_stream_id)
{
int r, c, r_nbr, c_nbr;
int slength = 1;
- struct ddir draindir, *founddir;
- int curr_stream;
+ CELL 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 };
+ char aspect;
r = stream_node[stream_id].r;
c = stream_node[stream_id].c;
+ if (next_stream_id)
+ *next_stream_id = stream_id;
/* get next downstream point */
- draindir.pos = INDEX(r, c);
- if (new_stream_pos)
- *new_stream_pos = draindir.pos;
- while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
- r_nbr = r + asp_r[(int)founddir->dir];
- c_nbr = c + asp_c[(int)founddir->dir];
+ bseg_get(&asp, &aspect, r, c);
+ while (aspect > 0) {
+ r_nbr = r + asp_r[(int)aspect];
+ c_nbr = c + asp_c[(int)aspect];
/* user-defined depression */
if (r_nbr == r && c_nbr == c)
@@ -104,43 +32,43 @@
if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
break;
/* next stream */
- if ((curr_stream = stream[INDEX(r_nbr, c_nbr)]) != stream_id) {
- if (new_stream_pos)
- *new_stream_pos = INDEX(r_nbr, c_nbr);
+ cseg_get(&stream, &curr_stream, r_nbr, c_nbr);
+ if (next_stream_id)
+ *next_stream_id = curr_stream;
+ if (curr_stream != stream_id)
break;
- }
slength++;
r = r_nbr;
c = c_nbr;
- if (new_stream_pos)
- *new_stream_pos = draindir.pos;
- draindir.pos = INDEX(r, c);
+ bseg_get(&asp, &aspect, r, c);
}
return slength;
}
-/* delete stream segment */
-int del_stream_seg(int stream_id)
+/* change downstream id: update or remove */
+int update_stream_id(CELL stream_id, CELL new_stream_id)
{
int i, r, c, r_nbr, c_nbr;
- unsigned int this_index;
- struct ddir draindir, *founddir;
- int curr_stream;
+ CELL new_stream = new_stream_id;
+ CELL 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 };
+ char aspect;
r = stream_node[stream_id].r;
c = stream_node[stream_id].c;
- this_index = INDEX(r, c);
- stream[this_index] = 0;
+ cseg_get(&stream, &curr_stream, r, c);
+ if (curr_stream != stream_id)
+ G_fatal_error("update downstream id: curr_stream != stream_id");
+ cseg_put(&stream, &new_stream, r, c);
curr_stream = stream_id;
/* get next downstream point */
- draindir.pos = this_index;
- while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
- r_nbr = r + asp_r[(int)founddir->dir];
- c_nbr = c + asp_c[(int)founddir->dir];
+ bseg_get(&asp, &aspect, r, c);
+ while (aspect > 0) {
+ r_nbr = r + asp_r[(int)aspect];
+ c_nbr = c + asp_c[(int)aspect];
/* user-defined depression */
if (r_nbr == r && c_nbr == c)
@@ -149,77 +77,118 @@
if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
break;
/* next stream */
- if ((curr_stream = stream[INDEX(r_nbr, c_nbr)]) != stream_id)
+ cseg_get(&stream, &curr_stream, r_nbr, c_nbr);
+ if (curr_stream != stream_id)
break;
r = r_nbr;
c = c_nbr;
- this_index = INDEX(r, c);
- stream[this_index] = 0;
- rbtree_remove(draintree, &draindir);
- draindir.pos = this_index;
+ cseg_put(&stream, &new_stream, r, c);
+ bseg_get(&asp, &aspect, r, c);
}
+
+ if (curr_stream <= 0)
+ return curr_stream;
/* update tributaries */
if (curr_stream != stream_id) {
+ int last_i = -1;
+
for (i = 0; i < stream_node[curr_stream].n_trib; i++) {
if (stream_node[curr_stream].trib[i] == stream_id) {
- stream_node[curr_stream].trib[i] = stream_node[curr_stream].trib[stream_node[curr_stream].n_trib - 1];
- stream_node[curr_stream].n_trib--;
+ if (new_stream_id) {
+ stream_node[curr_stream].trib[i] = new_stream_id;
+ }
+ else {
+ stream_node[curr_stream].n_trib--;
+ stream_node[curr_stream].trib[i] = stream_node[curr_stream].trib[stream_node[curr_stream].n_trib];
+ stream_node[curr_stream].trib[stream_node[curr_stream].n_trib] = 0;
+ }
+ last_i = i;
break;
}
}
+ for (i = 0; i < stream_node[curr_stream].n_trib; i++) {
+ if (stream_node[curr_stream].trib[i] == stream_id) {
+ G_warning("last_i %d, i %d", last_i, i);
+ G_warning("failed updating stream %d at node %d", stream_id, curr_stream);
+ }
+ }
}
- return 1;
+ return curr_stream;
}
-/* update downstream id */
-int update_stream_id(int stream_id, int new_stream_id)
+/* delete stream segments shorter than threshold */
+int del_streams(int min_length)
{
- int i, r, c, r_nbr, c_nbr;
- unsigned int this_index;
- 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 };
+ int i;
+ int n_deleted = 0;
+ CELL curr_stream, stream_id;
+ int other_trib, tmp_trib;
+ int slength;
- r = stream_node[stream_id].r;
- c = stream_node[stream_id].c;
- this_index = INDEX(r, c);
- stream[this_index] = new_stream_id;
- curr_stream = stream_id;
+ G_message(_("Delete stream segments shorter than %d cells..."), min_length);
- /* get next downstream point */
- draindir.pos = this_index;
- while ((founddir = rbtree_find(draintree, &draindir)) != NULL) {
- r_nbr = r + asp_r[(int)founddir->dir];
- c_nbr = c + asp_c[(int)founddir->dir];
+ /* TODO: proceed from stream heads to outlets
+ * -> use depth first post order traversal */
- /* user-defined depression */
- if (r_nbr == r && c_nbr == c)
- break;
- /* outside region */
- if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
- break;
- /* next stream */
- if ((curr_stream = stream[INDEX(r_nbr, c_nbr)]) != stream_id)
- break;
- r = r_nbr;
- c = c_nbr;
- this_index = INDEX(r, c);
- stream[this_index] = new_stream_id;
- draindir.pos = this_index;
- }
+ /* go through all nodes */
+ for (i = 1; i <= n_stream_nodes; i++) {
+ G_percent(i, n_stream_nodes, 2);
- /* update tributaries */
- if (curr_stream != stream_id) {
- for (i = 0; i < stream_node[curr_stream].n_trib; i++) {
- if (stream_node[curr_stream].trib[i] == stream_id) {
- stream_node[curr_stream].trib[i] = new_stream_id;
- break;
+ /* not a stream head */
+ if (stream_node[i].n_trib > 0)
+ continue;
+
+ /* already deleted */
+ cseg_get(&stream, &curr_stream, stream_node[i].r, stream_node[i].c);
+ if (curr_stream == 0)
+ continue;
+
+ /* get length counted as n cells */
+ if ((slength = seg_length(i, &curr_stream)) >= min_length)
+ continue;
+
+ stream_id = i;
+
+ /* check n sibling tributaries */
+ if (curr_stream != stream_id) {
+ /* only one sibling tributary */
+ if (stream_node[curr_stream].n_trib == 2) {
+ if (stream_node[curr_stream].trib[0] != stream_id)
+ other_trib = stream_node[curr_stream].trib[0];
+ else
+ other_trib = stream_node[curr_stream].trib[1];
+
+ /* other trib is also stream head */
+ if (stream_node[other_trib].n_trib == 0) {
+ /* use shorter one */
+ if (seg_length(other_trib, NULL) < slength) {
+ tmp_trib = stream_id;
+ stream_id = other_trib;
+ other_trib = tmp_trib;
+ }
+ }
+ update_stream_id(stream_id, 0);
+ n_deleted++;
+
+ /* update downstream IDs */
+ update_stream_id(curr_stream, other_trib);
}
+ /* more than one other tributary */
+ else {
+ update_stream_id(stream_id, 0);
+ n_deleted++;
+ }
}
+ /* stream head is also outlet */
+ else {
+ update_stream_id(stream_id, 0);
+ n_deleted++;
+ }
}
- return curr_stream;
+ G_verbose_message(_("%d stream segments deleted"), n_deleted);
+
+ return n_deleted;
}
Property changes on: grass-addons/grass7/raster/r.stream.extract/del_streams.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Deleted: grass-addons/grass7/raster/r.stream.extract/description.html
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/description.html 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/description.html 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,228 +0,0 @@
-<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 the optionally adjusted or
-weighed total contributing area in square meters or any other unit.
-<p>
-<dt><em>depression</em>
-<dd>Input map, optional: All non-NULL and non-zero cells will be regarded
-as real depressions. Streams will not be routed out of depressions. If an
-area is marked as depression but the elevation model has no depression
-at this location, streams will not stop there. If a flow accumulation map
-and a map with real depressions are provided, the flow accumulation map
-must match the depression map such that flow is not distributed out of
-the indicated depressions. It is recommended to use internally computed
-flow accumulation if a depression map is provided.
-<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_length</em>
-<dd>Minimum stream length in number of cells for first-order (head/spring)
-stream segments. All first-order stream segments shorter than
-<em>stream_length</em> will be deleted.
-
-<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 all non-NULL cells in
-input elevation. 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.
-</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 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
-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>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>
-
-<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$</i>
Modified: grass-addons/grass7/raster/r.stream.extract/do_astar.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/do_astar.c 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/do_astar.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,32 +1,27 @@
#include <stdlib.h>
#include <math.h>
-#include <grass/gis.h>
+#include <grass/raster.h>
#include <grass/glocale.h>
#include "local_proto.h"
-unsigned int get_parent(unsigned int c)
-{
- return (unsigned int)((c) - 2) / 3 + 1;
-}
+#define GET_PARENT(c) ((unsigned int)(((c) - 2) >> 2) + 1)
+#define GET_CHILD(p) ((unsigned int)((p) << 2) - 2)
-unsigned int get_child(unsigned int p)
-{
- return (unsigned int)(p) * 3 - 1;
-}
+HEAP_PNT heap_drop(void);
+int sift_up(unsigned int, HEAP_PNT);
+double get_slope(CELL, CELL, double);
-unsigned int heap_drop(void);
-
-double get_slope2(CELL, CELL, double);
-
int do_astar(void)
{
int r, c, r_nbr, c_nbr, ct_dir;
- unsigned int astp;
- int count, is_in_list, is_worked;
+ unsigned int first_cum, count;
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, ele_nbr[8];
- unsigned int thisindex, nindex;
+ WAT_ALT wa;
+ char asp_val;
+ char flag_value, is_in_list, is_worked;
+ HEAP_PNT heap_p;
/* sides
* |7|1|4|
* |2| |3|
@@ -61,31 +56,25 @@
}
ew_res = window.ew_res;
ns_res = window.ns_res;
-
+
while (heap_size > 0) {
G_percent(count++, n_points, 1);
if (count > n_points)
- G_fatal_error(_("BUG in A* Search: %d surplus points"), heap_size);
+ G_fatal_error(_("BUG in A* Search: %d surplus points"),
+ heap_size);
if (heap_size > n_points)
G_fatal_error
(_("BUG in A* Search: too many points in heap %d, should be %d"),
heap_size, n_points);
- astp = astar_pts[1];
+ heap_p = heap_drop();
- heap_drop();
+ r = heap_p.pnt.r;
+ c = heap_p.pnt.c;
- /* set flow accumulation order */
- astar_pts[first_cum] = astp;
- first_cum--;
+ ele_val = heap_p.ele;
- thisindex = astp;
- r = thisindex / ncols;
- c = thisindex - r * ncols;
-
- 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];
@@ -95,13 +84,14 @@
/* 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);
- is_worked = FLAG_GET(worked, r_nbr, c_nbr);
- nindex = INDEX(r_nbr, c_nbr);
+ bseg_get(&bitflags, &flag_value, r_nbr, c_nbr);
+ is_in_list = FLAG_GET(flag_value, INLISTFLAG);
+ is_worked = FLAG_GET(flag_value, WORKEDFLAG);
if (!is_worked) {
- ele_nbr[ct_dir] = ele[nindex];
+ seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
+ ele_nbr[ct_dir] = wa.ele;
slope[ct_dir] =
- get_slope2(ele_val, ele_nbr[ct_dir],
+ get_slope(ele_val, ele_nbr[ct_dir],
dist_to_nbr[ct_dir]);
}
/* avoid diagonal flow direction bias */
@@ -110,14 +100,14 @@
if (slope[nbr_ew[ct_dir]] > 0) {
/* slope to ew nbr > slope to center */
if (slope[ct_dir] <
- get_slope2(ele_nbr[nbr_ew[ct_dir]],
+ get_slope(ele_nbr[nbr_ew[ct_dir]],
ele_nbr[ct_dir], ew_res))
skip_me = 1;
}
if (!skip_me && slope[nbr_ns[ct_dir]] > 0) {
/* slope to ns nbr > slope to center */
if (slope[ct_dir] <
- get_slope2(ele_nbr[nbr_ns[ct_dir]],
+ get_slope(ele_nbr[nbr_ns[ct_dir]],
ele_nbr[ct_dir], ns_res))
skip_me = 1;
}
@@ -125,30 +115,47 @@
}
if (is_in_list == 0 && skip_me == 0) {
- ele_up = ele[nindex];
- asp[nindex] = drain[r_nbr - r + 1][c_nbr - c + 1];
- heap_add(r_nbr, c_nbr, ele_up, asp[nindex]);
+ ele_up = ele_nbr[ct_dir];
+ asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
+ bseg_put(&asp, &asp_val, r_nbr, c_nbr);
+ heap_add(r_nbr, c_nbr, ele_up);
+ FLAG_SET(flag_value, INLISTFLAG);
+ bseg_put(&bitflags, &flag_value, r_nbr, c_nbr);
}
else if (is_in_list && is_worked == 0) {
- /* neighbour is edge in list, not yet worked */
- if (asp[nindex] < 0) {
- asp[nindex] = drain[r_nbr - r + 1][c_nbr - c + 1];
+ if (FLAG_GET(flag_value, EDGEFLAG)) {
+ /* neighbour is edge in list, not yet worked */
+ bseg_get(&asp, &asp_val, r_nbr, c_nbr);
+ if (asp_val < 0) {
+ /* adjust flow direction for edge cell */
+ asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
+ bseg_put(&asp, &asp_val, r_nbr, c_nbr);
+ }
}
- /* neighbour is inside real depression, not yet worked */
- if (asp[nindex] == 0 && ele_val <= ele[nindex]) {
- asp[nindex] = drain[r_nbr - r + 1][c_nbr - c + 1];
+ else if (FLAG_GET(flag_value, DEPRFLAG)) {
+ G_debug(3, "real depression");
+ /* neighbour is inside real depression, not yet worked */
+ bseg_get(&asp, &asp_val, r_nbr, c_nbr);
+ if (asp_val == 0 && ele_val <= ele_nbr[ct_dir]) {
+ asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
+ bseg_put(&asp, &asp_val, r_nbr, c_nbr);
+ FLAG_UNSET(flag_value, DEPRFLAG);
+ bseg_put(&bitflags, &flag_value, r_nbr, c_nbr);
+ }
}
}
}
} /* end neighbours */
- FLAG_SET(worked, r, c);
- } /* end A* search */
+ /* add astar points to sorted list for flow accumulation and stream extraction */
+ first_cum--;
+ seg_put(&astar_pts, (char *)&heap_p.pnt, 0, first_cum);
+ bseg_get(&bitflags, &flag_value, r, c);
+ FLAG_SET(flag_value, WORKEDFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
+ } /* end A* search */
G_percent(n_points, n_points, 1); /* finish it */
- flag_destroy(in_list);
- G_free(astar_added);
-
return 1;
}
@@ -156,51 +163,40 @@
* compare function for heap
* returns 1 if point1 < point2 else 0
*/
-int heap_cmp(CELL ele1, unsigned int index1, CELL ele2, unsigned int index2)
+int heap_cmp(HEAP_PNT *a, HEAP_PNT *b)
{
- if (ele1 < ele2)
+ if (a->ele < b->ele)
return 1;
- else if (ele1 == ele2) {
- if (index1 < index2)
+ else if (a->ele == b->ele) {
+ if (a->added < b->added)
return 1;
}
-
return 0;
}
-int sift_up(unsigned int start, CELL elec)
+int sift_up(unsigned int start, HEAP_PNT child_p)
{
- unsigned int child, child_added, parent;
- CELL elep;
- unsigned int childp;
+ unsigned int parent, child;
+ HEAP_PNT heap_p;
child = start;
- child_added = astar_added[child];
- childp = astar_pts[child];
while (child > 1) {
- parent = get_parent(child);
+ parent = GET_PARENT(child);
+ seg_get(&search_heap, (char *)&heap_p, 0, parent);
- elep = ele[astar_pts[parent]];
-
- /* 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];
+ /* push parent point down if child is smaller */
+ if (heap_cmp(&child_p, &heap_p)) {
+ seg_put(&search_heap, (char *)&heap_p, 0, child);
child = parent;
}
+ /* no more sifting up, found slot for child */
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;
- }
+ /* add child to heap */
+ seg_put(&search_heap, (char *)&child_p, 0, child);
return 0;
}
@@ -209,21 +205,23 @@
* add item to heap
* returns heap_size
*/
-unsigned int heap_add(int r, int c, CELL ele, char asp)
+unsigned int heap_add(int r, int c, CELL ele)
{
- FLAG_SET(in_list, r, c);
-
+ HEAP_PNT heap_p;
+
/* add point to next free position */
heap_size++;
- astar_added[heap_size] = nxt_avail_pt;
- astar_pts[heap_size] = INDEX(r, c);
+ heap_p.added = nxt_avail_pt;
+ heap_p.ele = ele;
+ heap_p.pnt.r = r;
+ heap_p.pnt.c = c;
nxt_avail_pt++;
/* sift up: move new point towards top of heap */
- sift_up(heap_size, ele);
+ sift_up(heap_size, heap_p);
return heap_size;
}
@@ -232,62 +230,56 @@
* drop item from heap
* returns heap size
*/
-unsigned int heap_drop(void)
+HEAP_PNT heap_drop(void)
{
unsigned int child, childr, parent;
int i;
- CELL elec, eler;
+ HEAP_PNT child_p, childr_p, last_p, root_p;
+ seg_get(&search_heap, (char *)&last_p, 0, heap_size);
+ seg_get(&search_heap, (char *)&root_p, 0, 1);
+
if (heap_size == 1) {
- astar_added[1] = -1;
heap_size = 0;
- return heap_size;
+ return root_p;
}
parent = 1;
- while ((child = get_child(parent)) <= heap_size) {
+ while ((child = GET_CHILD(parent)) < heap_size) {
- elec = ele[astar_pts[child]];
+ seg_get(&search_heap, (char *)&child_p, 0, child);
if (child < heap_size) {
childr = child + 1;
- i = child + 3;
- while (childr <= heap_size && childr < i) {
- eler = ele[astar_pts[childr]];
-
- if (heap_cmp
- (eler, astar_added[childr], elec,
- astar_added[child]) == 1) {
+ i = child + 4;
+ while (childr < heap_size && childr < i) {
+ seg_get(&search_heap, (char *)&childr_p, 0, childr);
+ if (heap_cmp(&childr_p, &child_p)) {
child = childr;
- elec = eler;
+ child_p = childr_p;
}
childr++;
}
}
+ if (heap_cmp(&last_p, &child_p)) {
+ break;
+ }
+
/* move hole down */
- astar_added[parent] = astar_added[child];
- astar_pts[parent] = astar_pts[child];
+ seg_put(&search_heap, (char *)&child_p, 0, parent);
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];
+ seg_put(&search_heap, (char *)&last_p, 0, parent);
- elec = ele[astar_pts[parent]];
- /* 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;
+ return root_p;
}
-double get_slope2(CELL ele, CELL up_ele, double dist)
+double get_slope(CELL ele, CELL up_ele, double dist)
{
if (ele >= up_ele)
return 0.0;
Property changes on: grass-addons/grass7/raster/r.stream.extract/do_astar.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.stream.extract/dseg.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/dseg.c (rev 0)
+++ grass-addons/grass7/raster/r.stream.extract/dseg.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -0,0 +1,156 @@
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/glocale.h>
+#include "seg.h"
+
+int dseg_open(DSEG *dseg, int srows, int scols, int nsegs_in_memory)
+{
+ char *filename;
+ int errflag;
+ int fd;
+
+ dseg->filename = NULL;
+ dseg->fd = -1;
+ dseg->name = NULL;
+ dseg->mapset = NULL;
+
+ filename = G_tempfile();
+ if (-1 == (fd = creat(filename, 0666))) {
+ G_warning(_("dseg_open(): unable to create segment file"));
+ return -2;
+ }
+ if (0 >
+ (errflag =
+ segment_format(fd, Rast_window_rows(), Rast_window_cols(), srows, scols,
+ sizeof(DCELL)))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning(_("dseg_open(): could not write segment file"));
+ return -1;
+ }
+ else {
+ G_warning(_("dseg_open(): illegal configuration parameter(s)"));
+ return -3;
+ }
+ }
+ close(fd);
+ if (-1 == (fd = open(filename, 2))) {
+ unlink(filename);
+ G_warning(_("dseg_open(): unable to re-open segment file"));
+ return -4;
+ }
+ if (0 > (errflag = segment_init(&(dseg->seg), fd, nsegs_in_memory))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning(_("dseg_open(): could not read segment file"));
+ return -5;
+ }
+ else {
+ G_warning(_("dseg_open(): out of memory"));
+ return -6;
+ }
+ }
+ dseg->filename = filename;
+ dseg->fd = fd;
+ return 0;
+}
+
+int dseg_close(DSEG *dseg)
+{
+ segment_release(&(dseg->seg));
+ close(dseg->fd);
+ unlink(dseg->filename);
+ if (dseg->name) {
+ G_free(dseg->name);
+ dseg->name = NULL;
+ }
+ if (dseg->mapset) {
+ G_free(dseg->mapset);
+ dseg->mapset = NULL;
+ }
+ return 0;
+}
+
+int dseg_put(DSEG *dseg, DCELL *value, int row, int col)
+{
+ if (segment_put(&(dseg->seg), (DCELL *) value, row, col) < 0) {
+ G_warning(_("dseg_put(): could not write segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int dseg_put_row(DSEG *dseg, DCELL *value, int row)
+{
+ if (segment_put_row(&(dseg->seg), (DCELL *) value, row) < 0) {
+ G_warning(_("dseg_put(): could not write segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int dseg_get(DSEG *dseg, DCELL *value, int row, int col)
+{
+ if (segment_get(&(dseg->seg), (DCELL *) value, row, col) < 0) {
+ G_warning(_("dseg_get(): could not read segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int dseg_read_raster(DSEG *dseg, char *map_name, char *mapset)
+{
+ int row, nrows, ncols;
+ int map_fd;
+ DCELL *dbuffer;
+
+ dseg->name = NULL;
+ dseg->mapset = NULL;
+
+ map_fd = Rast_open_old(map_name, mapset);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ dbuffer = Rast_allocate_d_buf();
+ for (row = 0; row < nrows; row++) {
+ Rast_get_d_row(map_fd, dbuffer, row);
+ if (segment_put_row(&(dseg->seg), (DCELL *) dbuffer, row) < 0) {
+ G_free(dbuffer);
+ Rast_close(map_fd);
+ G_warning(_("dseg_read_raster(): unable to segment put row for <%s> in <%s>"),
+ map_name, mapset);
+ return (-1);
+ }
+ }
+
+ Rast_close(map_fd);
+ G_free(dbuffer);
+
+ dseg->name = G_store(map_name);
+ dseg->mapset = G_store(mapset);
+
+ return 0;
+}
+
+int dseg_write_cellfile(DSEG *dseg, char *map_name)
+{
+ int map_fd;
+ int row, nrows, ncols;
+ DCELL *dbuffer;
+
+ map_fd = Rast_open_new(map_name, DCELL_TYPE);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ dbuffer = Rast_allocate_d_buf();
+ segment_flush(&(dseg->seg));
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 1);
+ segment_get_row(&(dseg->seg), (DCELL *) dbuffer, row);
+ Rast_put_row(map_fd, dbuffer, DCELL_TYPE);
+ }
+ G_percent(row, nrows, 1); /* finish it */
+ G_free(dbuffer);
+ Rast_close(map_fd);
+ return 0;
+}
Property changes on: grass-addons/grass7/raster/r.stream.extract/dseg.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Deleted: grass-addons/grass7/raster/r.stream.extract/flag.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/flag.c 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/flag.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,49 +0,0 @@
-#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);
-}
-
-int flag_destroy(FLAG * flags)
-{
- G_free(flags->array[0]);
- G_free(flags->array);
- G_free(flags);
-
- return 0;
-}
-
-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;
-}
Modified: grass-addons/grass7/raster/r.stream.extract/flag.h
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/flag.h 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/flag.h 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,63 +1,36 @@
#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
+/* a set of routines that allow the programmer to "flag" cells in a
+ * raster map. A flag is of type unsigned char, i.e. 8 bits can be set.
+ *
+ * int flag_set(flag, bitno)
+ * sets the flag at position bitno to one.
+ *
+ * int flag_unset(flag, bitno)
+ * sets the flag at position bitno to zero.
+ *
+ * int flag_get(flag, bitno)
+ * checks if the flag is set at postion bitno.
+ *
+ * Examples:
+ * set flag at position 0: FLAG_SET(flag, 0)
+ * unset (clear) flag at position 7: FLAG_UNSET(flag, 7)
+ * check flag at position 5: is_set_at_5 = FLAG_GET(flag, 5)
*/
-#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))
+/* flag positions */
+#define NULLFLAG 0 /* elevation is NULL */
+#define EDGEFLAG 1 /* edge cell */
+#define INLISTFLAG 2 /* in open A* list */
+#define WORKEDFLAG 3 /* in closed A* list/ accumulation done */
+#define STREAMFLAG 4 /* stream */
+#define DEPRFLAG 5 /* real depression */
+#define WORKED2FLAG 6 /* extraction done */
+/* last bit is unused */
-#define FLAG_SET(flags,row,col) \
- (flags)->array[(row)][(col)>>3] |= (1<<((col) & 7))
+#define FLAG_SET(flag,bitno) ((flag) |= (1 << (bitno)))
+#define FLAG_UNSET(flag,bitno) ((flag) &= ~(1 << (bitno)))
+#define FLAG_GET(flag,bitno) ((flag) & (1 << (bitno)))
-#define FLAG_GET(flags,row,col) \
- (flags)->array[(row)][(col)>>3] & (1<<((col) & 7))
-
-/* flag.c */
-int flag_clear_all(FLAG *);
-FLAG *flag_create(int, int);
-int flag_destroy(FLAG *);
-
#endif /* __FLAG_H__ */
Property changes on: grass-addons/grass7/raster/r.stream.extract/flag.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.stream.extract/init_search.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/init_search.c (rev 0)
+++ grass-addons/grass7/raster/r.stream.extract/init_search.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -0,0 +1,123 @@
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+int init_search(int depr_fd)
+{
+ int r, c, r_nbr, c_nbr, ct_dir;
+ CELL *depr_buf, ele_value;
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ char flag_value, flag_value_nbr, is_null;
+ WAT_ALT wa;
+ char asp_value;
+ unsigned int n_depr_cells = 0;
+
+ nxt_avail_pt = heap_size = 0;
+
+ /* load edge cells and real depressions to A* heap */
+ if (depr_fd >= 0)
+ depr_buf = Rast_allocate_buf(CELL_TYPE);
+ else
+ depr_buf = NULL;
+
+ G_message(_("Initialize A* Search..."));
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 2);
+
+ if (depr_fd >= 0) {
+ Rast_get_row(depr_fd, depr_buf, r, CELL_TYPE);
+ }
+
+ for (c = 0; c < ncols; c++) {
+
+ bseg_get(&bitflags, &flag_value, r, c);
+ is_null = FLAG_GET(flag_value, NULLFLAG);
+
+ if (is_null)
+ 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;
+
+ seg_get(&watalt, (char *)&wa, r, c);
+ ele_value = wa.ele;
+ heap_add(r, c, ele_value);
+ FLAG_SET(flag_value, INLISTFLAG);
+ FLAG_SET(flag_value, EDGEFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
+ bseg_put(&asp, &asp_value, r, c);
+ continue;
+ }
+
+ /* any neighbour NULL ? */
+ asp_value = 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];
+
+ bseg_get(&bitflags, &flag_value_nbr, r_nbr, c_nbr);
+ is_null = FLAG_GET(flag_value_nbr, NULLFLAG);
+
+ if (is_null) {
+ asp_value = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+ seg_get(&watalt, (char *)&wa, r, c);
+ ele_value = wa.ele;
+ heap_add(r, c, ele_value);
+ FLAG_SET(flag_value, INLISTFLAG);
+ FLAG_SET(flag_value, EDGEFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
+ bseg_put(&asp, &asp_value, r, c);
+
+ break;
+ }
+ }
+ if (asp_value) /* some neighbour was NULL, point added to list */
+ continue;
+
+ /* real depression ? */
+ if (depr_fd >= 0) {
+ if (!Rast_is_c_null_value(&depr_buf[c]) && depr_buf[c] != 0) {
+ seg_get(&watalt, (char *)&wa, r, c);
+ ele_value = wa.ele;
+ heap_add(r, c, ele_value);
+ FLAG_SET(flag_value, INLISTFLAG);
+ FLAG_SET(flag_value, DEPRFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
+ bseg_put(&asp, &asp_value, r, c);
+ n_depr_cells++;
+ }
+ }
+ }
+ }
+ G_percent(nrows, nrows, 2); /* finish it */
+
+ if (depr_fd >= 0) {
+ Rast_close(depr_fd);
+ G_free(depr_buf);
+ }
+
+ G_debug(1, "%d edge cells", heap_size - n_depr_cells);
+ if (n_depr_cells)
+ G_debug(1, "%d cells in depressions", n_depr_cells);
+
+ return 1;
+}
Property changes on: grass-addons/grass7/raster/r.stream.extract/init_search.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/raster/r.stream.extract/load.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/load.c 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/load.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,4 +1,4 @@
-#include <grass/gis.h>
+#include <grass/raster.h>
#include <grass/glocale.h>
#include "local_proto.h"
@@ -6,16 +6,10 @@
int ele_round(double x)
{
- int n;
-
if (x >= 0.0)
- n = x + .5;
- else {
- n = -x + .5;
- n = -n;
- }
-
- return n;
+ return x + .5;
+ else
+ return x - .5;
}
/*
@@ -23,33 +17,28 @@
* gets start points for A* Search
* start points are edges
*/
-int load_maps(int ele_fd, int acc_fd, int depr_fd)
+int load_maps(int ele_fd, int acc_fd)
{
- int r, c, thisindex;
- char asp_value, *aspp;
+ int r, c;
void *ele_buf, *ptr, *acc_buf = NULL, *acc_ptr = NULL;
- CELL *loadp, ele_value;
- CELL *depr_buf;
- 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;
+ CELL ele_value, *stream_id;
+ DCELL dvalue, acc_value;
size_t ele_size, acc_size = 0;
int ele_map_type, acc_map_type = 0;
- DCELL *accp;
+ WAT_ALT *wabuf;
+ char *flag_value_buf, *aspect;
if (acc_fd < 0)
- G_message(_("Load elevation map and get start points..."));
+ G_message(_("Load elevation map..."));
else
- G_message(_("Load input maps and get start points..."));
+ G_message(_("Load input maps..."));
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);
+ ele_map_type = Rast_get_map_type(ele_fd);
+ ele_size = Rast_cell_size(ele_map_type);
+ ele_buf = Rast_allocate_buf(ele_map_type);
if (ele_buf == NULL) {
G_warning(_("Could not allocate memory"));
@@ -57,9 +46,9 @@
}
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);
+ acc_map_type = Rast_get_map_type(acc_fd);
+ acc_size = Rast_cell_size(acc_map_type);
+ acc_buf = Rast_allocate_buf(acc_map_type);
if (acc_buf == NULL) {
G_warning(_("Could not allocate memory"));
return -1;
@@ -70,73 +59,71 @@
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);
+ wabuf = G_malloc(ncols * sizeof(WAT_ALT));
+ flag_value_buf = G_malloc(ncols * sizeof(char));
+ stream_id = G_malloc(ncols * sizeof(CELL));
+ aspect = G_malloc(ncols * sizeof(char));
- loadp = ele;
- accp = acc;
- aspp = asp;
-
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;
- }
+ Rast_get_row(ele_fd, ele_buf, r, ele_map_type);
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;
- }
+ Rast_get_row(acc_fd, acc_buf, r, acc_map_type);
acc_ptr = acc_buf;
}
for (c = 0; c < ncols; c++) {
- FLAG_UNSET(worked, r, c);
- FLAG_UNSET(in_list, r, c);
+ flag_value_buf[c] = 0;
+ aspect[c] = 0;
+ stream_id[c] = 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 (Rast_is_null_value(ptr, ele_map_type)) {
+ FLAG_SET(flag_value_buf[c], NULLFLAG);
+ FLAG_SET(flag_value_buf[c], INLISTFLAG);
+ FLAG_SET(flag_value_buf[c], WORKEDFLAG);
+ FLAG_SET(flag_value_buf[c], WORKED2FLAG);
+ Rast_set_c_null_value(&ele_value, 1);
+ /* flow accumulation */
+ Rast_set_d_null_value(&acc_value, 1);
}
else {
- if (ele_map_type == CELL_TYPE) {
- *loadp = *((CELL *) ptr);
- }
- else if (ele_map_type == FCELL_TYPE) {
+ switch (ele_map_type) {
+ case CELL_TYPE:
+ ele_value = *((CELL *) ptr);
+ break;
+ case FCELL_TYPE:
dvalue = *((FCELL *) ptr);
dvalue *= ele_scale;
- *loadp = ele_round(dvalue);
- }
- else if (ele_map_type == DCELL_TYPE) {
+ ele_value = ele_round(dvalue);
+ break;
+ case DCELL_TYPE:
dvalue = *((DCELL *) ptr);
dvalue *= ele_scale;
- *loadp = ele_round(dvalue);
+ ele_value = ele_round(dvalue);
+ break;
}
if (acc_fd < 0)
- *accp = 1;
+ acc_value = 1;
else {
- if (G_is_null_value(acc_ptr, acc_map_type))
+ if (Rast_is_null_value(acc_ptr, acc_map_type))
G_fatal_error(_("Accumulation map does not match elevation map!"));
switch (acc_map_type) {
case CELL_TYPE:
- *accp = *((CELL *) acc_ptr);
+ acc_value = *((CELL *) acc_ptr);
break;
case FCELL_TYPE:
- *accp = *((FCELL *) acc_ptr);
+ acc_value = *((FCELL *) acc_ptr);
break;
case DCELL_TYPE:
- *accp = *((DCELL *) acc_ptr);
+ acc_value = *((DCELL *) acc_ptr);
break;
}
}
@@ -144,130 +131,32 @@
n_points++;
}
- loadp++;
- accp++;
+ wabuf[c].wat = acc_value;
+ wabuf[c].ele = ele_value;
ptr = G_incr_void_ptr(ptr, ele_size);
- *aspp = 0;
- aspp++;
if (acc_fd >= 0)
acc_ptr = G_incr_void_ptr(acc_ptr, acc_size);
}
+ seg_put_row(&watalt, (char *) wabuf, r);
+ bseg_put_row(&asp, aspect, r);
+ cseg_put_row(&stream, stream_id, r);
+ bseg_put_row(&bitflags, flag_value_buf, r);
}
G_percent(nrows, nrows, 1); /* finish it */
- G_close_cell(ele_fd);
+ Rast_close(ele_fd);
G_free(ele_buf);
+ G_free(wabuf);
+ G_free(flag_value_buf);
+ G_free(stream_id);
+ G_free(aspect);
if (acc_fd >= 0) {
- G_close_cell(acc_fd);
+ Rast_close(acc_fd);
G_free(acc_buf);
}
-
- astar_pts =
- (unsigned int *)G_malloc((n_points + 1) *
- sizeof(unsigned int));
-
- /* 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 and real depressions to A* heap */
- if (depr_fd >= 0)
- depr_buf = G_allocate_raster_buf(CELL_TYPE);
- else
- depr_buf = NULL;
-
- G_message(_("Set edge points..."));
- loadp = ele;
- for (r = 0; r < nrows; r++) {
- G_percent(r, nrows, 2);
-
- if (depr_fd >= 0) {
- if (G_get_raster_row(depr_fd, depr_buf, r, CELL_TYPE) < 0) {
- G_warning(_("Could not read raster map at row <%d>"), r);
- return -1;
- }
- }
-
- 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];
- asp[thisindex] = asp_value;
- heap_add(r, c, ele_value, asp_value);
- continue;
- }
-
- /* any neighbour NULL ? */
- asp_value = 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];
-
- 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];
- asp[thisindex] = asp_value;
- heap_add(r, c, ele_value, asp_value);
-
- break;
- }
- }
- if (asp_value) /* some neighbour was NULL, point added to list */
- continue;
-
- /* real depression ? */
- if (depr_fd >= 0) {
- if (!G_is_c_null_value(&depr_buf[c]) && depr_buf[c] != 0) {
- thisindex = INDEX(r, c);
- ele_value = ele[thisindex];
- asp[thisindex] = 0;
- heap_add(r, c, ele_value, 0);
- }
- }
- }
- }
- G_percent(nrows, nrows, 2); /* finish it */
-
- if (depr_fd >= 0) {
- G_close_cell(depr_fd);
- G_free(depr_buf);
- }
-
- G_debug(1, "%d edge (and depression) cells", heap_size);
+
G_debug(1, "%d non-NULL cells", n_points);
- return 1;
+ return n_points;
}
Property changes on: grass-addons/grass7/raster/r.stream.extract/load.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/raster/r.stream.extract/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/local_proto.h 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/local_proto.h 2010-12-16 14:01:19 UTC (rev 44615)
@@ -2,23 +2,31 @@
#ifndef __LOCAL_PROTO_H__
#define __LOCAL_PROTO_H__
+#include <grass/raster.h>
#include "flag.h"
-#include "rbtree.h"
+#include "seg.h"
#define INDEX(r, c) ((r) * ncols + (c))
#define MAXDEPTH 1000 /* maximum supported tree depth of stream network */
-struct ddir
-{
- int pos;
- int dir;
+#define POINT struct a_point
+POINT {
+ int r, c;
};
-struct point
-{
- int r, c;
+#define HEAP_PNT struct heap_point
+HEAP_PNT {
+ int added;
+ CELL ele;
+ POINT pnt;
};
+#define WAT_ALT struct wat_altitude
+WAT_ALT {
+ CELL ele;
+ DCELL wat;
+};
+
struct snode
{
int r, c;
@@ -27,40 +35,41 @@
int n_trib_total; /* number of all upstream stream segments */
int n_alloc; /* n allocated tributaries */
int *trib;
- double *acc;
} *stream_node;
/* global variables */
extern int nrows, ncols;
-extern unsigned int *astar_pts;
extern unsigned int n_search_points, n_points, nxt_avail_pt;
-extern unsigned int heap_size, *astar_added;
+extern unsigned int heap_size;
extern unsigned int n_stream_nodes, n_alloc_nodes;
-extern struct point *outlets;
+extern POINT *outlets;
extern unsigned int n_outlets, n_alloc_outlets;
-extern DCELL *acc;
-extern CELL *ele;
-extern char *asp;
-extern CELL *stream;
-extern FLAG *worked, *in_list;
extern char drain[3][3];
-extern unsigned int first_cum;
extern char sides;
extern int c_fac;
extern int ele_scale;
extern int have_depressions;
-extern struct RB_TREE *draintree;
+extern SSEG search_heap;
+extern SSEG astar_pts;
+extern BSEG bitflags;
+extern SSEG watalt;
+extern BSEG asp;
+extern CSEG stream;
+
/* load.c */
-int load_maps(int, int, int);
+int load_maps(int, int);
+/* init_search.c */
+int init_search(int);
+
/* do_astar.c */
int do_astar(void);
-unsigned int heap_add(int, int, CELL, char);
+unsigned int heap_add(int, int, CELL);
/* streams.c */
int do_accum(double);
-int extract_streams(double, double, int);
+int extract_streams(double, double, int, int);
/* thin.c */
int thin_streams(void);
@@ -70,9 +79,6 @@
/* del_streams.c */
int del_streams(int);
-int seg_length(int, unsigned int *);
-int del_stream_seg(int);
-int update_stream_id(int, int);
/* close.c */
int close_maps(char *, char *, char *);
Property changes on: grass-addons/grass7/raster/r.stream.extract/local_proto.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/raster/r.stream.extract/main.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/main.c 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/main.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -16,32 +16,33 @@
#include <stdlib.h>
#include <string.h>
#include <float.h>
-#include <grass/gis.h>
+#include <math.h>
+#include <grass/raster.h>
#include <grass/glocale.h>
#include "local_proto.h"
/* global variables */
int nrows, ncols;
-unsigned int *astar_pts;
unsigned int n_search_points, n_points, nxt_avail_pt;
-unsigned int heap_size, *astar_added;
+unsigned int heap_size;
unsigned int n_stream_nodes, n_alloc_nodes;
-struct point *outlets;
+POINT *outlets;
unsigned int n_outlets, n_alloc_outlets;
-DCELL *acc;
-CELL *ele;
-char *asp;
-CELL *stream;
-FLAG *worked, *in_list;
char drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
-unsigned int first_cum;
char sides;
int c_fac;
int ele_scale;
int have_depressions;
-struct RB_TREE *draintree;
+SSEG search_heap;
+SSEG astar_pts;
+BSEG bitflags;
+SSEG watalt;
+BSEG asp;
+CSEG stream;
+CELL *astar_order;
+
int main(int argc, char *argv[])
{
struct
@@ -50,6 +51,7 @@
struct Option *threshold, *d8cut;
struct Option *mont_exp;
struct Option *min_stream_length;
+ struct Option *memory;
} input;
struct
{
@@ -60,14 +62,17 @@
struct GModule *module;
int ele_fd, acc_fd, depr_fd;
double threshold, d8cut, mont_exp;
- int min_stream_length = 0;
- char *mapset;
+ int min_stream_length = 0, memory;
+ int seg_cols, seg_rows;
+ int num_open_segs, num_open_array_segs, num_seg_total;
+ const char *mapset;
G_gisinit(argv[0]);
/* Set description */
module = G_define_module();
- module->keywords = _("raster");
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("hydrology"));
module->description = _("Stream network extraction");
input.ele = G_define_standard_option(G_OPT_R_INPUT);
@@ -126,6 +131,13 @@
input.min_stream_length->description =
_("Applies only to first-order stream segments (springs/stream heads).");
+ input.memory = G_define_option();
+ input.memory->key = "memory";
+ input.memory->type = TYPE_INTEGER;
+ input.memory->required = NO;
+ input.memory->answer = "300";
+ input.memory->description = _("Maximum memory to be used in MB");
+
output.stream_rast = G_define_standard_option(G_OPT_R_OUTPUT);
output.stream_rast->key = "stream_rast";
output.stream_rast->description =
@@ -155,16 +167,16 @@
/***********************/
/* input maps exist ? */
- if (!G_find_cell(input.ele->answer, ""))
+ if (!G_find_raster(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, ""))
+ if (!G_find_raster(input.acc->answer, ""))
G_fatal_error(_("Raster map <%s> not found"), input.acc->answer);
}
if (input.depression->answer) {
- if (!G_find_cell(input.depression->answer, ""))
+ if (!G_find_raster(input.depression->answer, ""))
G_fatal_error(_("Raster map <%s> not found"), input.depression->answer);
have_depressions = 1;
}
@@ -209,6 +221,16 @@
}
else
min_stream_length = 0;
+
+ if (input.memory->answer) {
+ memory = atoi(input.memory->answer);
+ if (memory <= 0)
+ G_fatal_error(_("Memory must be positive but is %d"),
+ memory);
+ }
+ else
+ memory = 300;
+ memory = 4000;
/* Check for some output map */
if ((output.stream_rast->answer == NULL)
@@ -222,14 +244,14 @@
/*********************/
/* open input maps */
- mapset = G_find_cell2(input.ele->answer, "");
- ele_fd = G_open_cell_old(input.ele->answer, mapset);
+ mapset = G_find_raster2(input.ele->answer, "");
+ ele_fd = Rast_open_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);
+ mapset = G_find_raster2(input.acc->answer, "");
+ acc_fd = Rast_open_old(input.acc->answer, mapset);
if (acc_fd < 0)
G_fatal_error(_("Could not open input map %s"),
input.acc->answer);
@@ -238,8 +260,8 @@
acc_fd = -1;
if (input.depression->answer) {
- mapset = G_find_cell2(input.depression->answer, "");
- depr_fd = G_open_cell_old(input.depression->answer, mapset);
+ mapset = G_find_raster2(input.depression->answer, "");
+ depr_fd = Rast_open_old(input.depression->answer, mapset);
if (depr_fd < 0)
G_fatal_error(_("Could not open input map %s"),
input.depression->answer);
@@ -248,49 +270,132 @@
depr_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 */
+ nrows = Rast_window_rows();
+ ncols = Rast_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));
- asp = (char *) G_malloc(nrows * ncols * sizeof(char));
- acc = (DCELL *) G_malloc(nrows * ncols * sizeof(DCELL));
+ /* segment structures */
+ seg_rows = seg_cols = 64;
+ /* elevation + accumulation: 12 byte -> 48 KB / segment
+ * aspect: 1 byte -> 4 KB / segment
+ * stream: 4 byte -> 16 KB / segment
+ * flag: 1 byte -> 4 KB / segment
+ *
+ * Total MB / segment so far: <0.082
+ *
+ * astar_points: 9 byte -> XX KB / segment
+ * heap_points: 17 byte -> XX KB / segment
+ */
+
+ num_open_segs = memory / 0.082;
+ num_seg_total = (ncols / seg_cols + 1) * (nrows / seg_rows + 1);
+ if (num_open_segs > num_seg_total)
+ num_open_segs = num_seg_total;
+ G_verbose_message(_("%d of %d segments are kept in memory"),
+ num_open_segs, num_seg_total);
+ /* open segment files */
+ G_verbose_message(_("Create temporary files..."));
+ seg_open(&watalt, nrows, ncols, seg_rows, seg_cols, num_open_segs,
+ sizeof(WAT_ALT), 1);
+ cseg_open(&stream, seg_rows, seg_cols, num_open_segs);
+ bseg_open(&asp, seg_rows, seg_cols, num_open_segs);
+ bseg_open(&bitflags, seg_rows, seg_cols, num_open_segs);
+
/* load maps */
- if (load_maps(ele_fd, acc_fd, depr_fd) < 0)
+ if (load_maps(ele_fd, acc_fd) < 0)
G_fatal_error(_("Could not load input map(s)"));
+ else if (!n_points)
+ G_fatal_error(_("No non-NULL cells in input map(s)"));
+ G_debug(1, "open segments for A* points");
+ /* rounded down power of 2 */
+ seg_cols = (int) (pow(2, (int)(log(num_open_segs / 8.0) / log(2) + 0.5)) + 0.5);
+ if (seg_cols < 2)
+ seg_cols = 2;
+ num_open_array_segs = num_open_segs / seg_cols;
+ if (num_open_array_segs == 0)
+ num_open_array_segs = 1;
+ /* n cols in segment */
+ seg_cols *= seg_rows * seg_rows;
+ /* n segments in row */
+ num_seg_total = n_points / seg_cols;
+ if (n_points % seg_cols > 0)
+ num_seg_total++;
+ /* no need to have more segments open than exist */
+ if (num_open_array_segs > num_seg_total)
+ num_open_array_segs = num_seg_total;
+
+ if (num_open_array_segs > 4)
+ num_open_array_segs = 4;
+
+ G_debug(0, "segment size for A* points: %d", seg_cols);
+ seg_open(&astar_pts, 1, n_points, 1, seg_cols, num_open_array_segs,
+ sizeof(POINT), 1);
+
+ /* one-based d-ary search_heap with astar_pts */
+ G_debug(1, "open segments for A* search heap");
+ /* rounded down power of 2 */
+ seg_cols = (int) (pow(2, (int)(log(num_open_segs / 8.0) / log(2) + 0.5)) + 0.5);
+ if (seg_cols < 2)
+ seg_cols = 2;
+ num_open_array_segs = num_open_segs / seg_cols;
+ if (num_open_array_segs == 0)
+ num_open_array_segs = 1;
+ /* n cols in segment */
+ seg_cols *= seg_rows * seg_rows;
+ /* n segments in row */
+ num_seg_total = (n_points + 1) / seg_cols;
+ if ((n_points + 1) % seg_cols > 0)
+ num_seg_total++;
+ /* no need to have more segments open than exist */
+ if (num_open_array_segs > num_seg_total)
+ num_open_array_segs = num_seg_total;
+
+ G_debug(0, "A* search heap open segments %d, target 8, total %d",
+ num_open_array_segs, num_seg_total);
+ G_debug(0, "segment size for heap points: %d", seg_cols);
+ /* the search heap will not hold more than 5% of all points at any given time ? */
+ /* chances are good that the heap will fit into one large segment */
+ seg_open(&search_heap, 1, n_points + 1, 1, seg_cols,
+ num_open_array_segs, sizeof(HEAP_PNT), 1);
+
/********************/
/* processing */
/********************/
+ /* initialize A* search */
+ if (init_search(depr_fd) < 0)
+ G_fatal_error(_("Could not initialize search"));
+
/* sort elevation and get initial stream direction */
if (do_astar() < 0)
G_fatal_error(_("Could not sort elevation map"));
+ seg_close(&search_heap);
- /* extract streams */
if (acc_fd < 0) {
+ /* accumulate surface flow */
if (do_accum(d8cut) < 0)
G_fatal_error(_("Could not calculate flow accumulation"));
}
- stream = (CELL *) G_malloc(nrows * ncols * sizeof(CELL));
+ /* extract streams */
if (extract_streams
- (threshold, mont_exp, min_stream_length) < 0)
+ (threshold, mont_exp, min_stream_length, acc_fd < 0) < 0)
G_fatal_error(_("Could not extract streams"));
- G_free(acc);
+ seg_close(&astar_pts);
+ seg_close(&watalt);
/* thin streams */
if (thin_streams() < 0)
- G_fatal_error(_("Could not extract streams"));
+ G_fatal_error(_("Could not thin streams"));
/* delete short streams */
if (min_stream_length) {
if (del_streams(min_stream_length) < 0)
- G_fatal_error(_("Could not extract streams"));
+ G_fatal_error(_("Could not delete short stream segments"));
}
/* write output maps */
@@ -298,9 +403,9 @@
output.dir_rast->answer) < 0)
G_fatal_error(_("Could not write output maps"));
- G_free(ele);
- G_free(stream);
- G_free(asp);
+ bseg_close(&asp);
+ cseg_close(&stream);
+ bseg_close(&bitflags);
exit(EXIT_SUCCESS);
}
Property changes on: grass-addons/grass7/raster/r.stream.extract/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Copied: grass-addons/grass7/raster/r.stream.extract/r.stream.extract.html (from rev 44614, grass-addons/grass7/raster/r.stream.extract/description.html)
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/r.stream.extract.html (rev 0)
+++ grass-addons/grass7/raster/r.stream.extract/r.stream.extract.html 2010-12-16 14:01:19 UTC (rev 44615)
@@ -0,0 +1,228 @@
+<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 the optionally adjusted or
+weighed total contributing area in square meters or any other unit.
+<p>
+<dt><em>depression</em>
+<dd>Input map, optional: All non-NULL and non-zero cells will be regarded
+as real depressions. Streams will not be routed out of depressions. If an
+area is marked as depression but the elevation model has no depression
+at this location, streams will not stop there. If a flow accumulation map
+and a map with real depressions are provided, the flow accumulation map
+must match the depression map such that flow is not distributed out of
+the indicated depressions. It is recommended to use internally computed
+flow accumulation if a depression map is provided.
+<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_length</em>
+<dd>Minimum stream length in number of cells for first-order (head/spring)
+stream segments. All first-order stream segments shorter than
+<em>stream_length</em> will be deleted.
+
+<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 all non-NULL cells in
+input elevation. 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.
+</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 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
+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>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>
+
+<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$</i>
Property changes on: grass-addons/grass7/raster/r.stream.extract/r.stream.extract.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Date
Added: svn:eol-style
+ native
Deleted: grass-addons/grass7/raster/r.stream.extract/rbtree.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/rbtree.c 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/rbtree.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,550 +0,0 @@
-/*!
- * \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 */
-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 Search 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)
-{
- struct RB_NODE *it;
- struct RB_NODE *save = tree->root;
-
- /*
- Rotate away the left links so that
- we can treat this like the destruction
- of a linked list
- */
- while((it = save) != NULL) {
- if (it->link[0] == NULL) {
- /* No left links, just kill the node and move on */
- save = it->link[1];
- G_free(it->data);
- it->data = NULL;
- G_free(it);
- it = NULL;
- }
- else {
- /* Rotate away the left link and check again */
- save = it->link[0];
- it->link[0] = save->link[1];
- save->link[1] = it;
- }
- }
- G_free(tree);
- tree = NULL;
-}
-
-/* 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);
-}
Deleted: grass-addons/grass7/raster/r.stream.extract/rbtree.h
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/rbtree.h 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/rbtree.h 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,112 +0,0 @@
-/*************************************************************
- * 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/grass7/raster/r.stream.extract/seg.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/seg.c (rev 0)
+++ grass-addons/grass7/raster/r.stream.extract/seg.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -0,0 +1,112 @@
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/glocale.h>
+#include "seg.h"
+
+int
+seg_open(SSEG *sseg, int nrows, int ncols, int row_in_seg, int col_in_seg,
+ int nsegs_in_memory, int size_struct, int fill)
+{
+ char *filename;
+ int errflag;
+ int fd;
+
+ sseg->filename = NULL;
+ sseg->fd = -1;
+
+ filename = G_tempfile();
+ if (-1 == (fd = creat(filename, 0666))) {
+ G_warning(_("seg_open(): unable to create segment file"));
+ return -2;
+ }
+ if (fill)
+ errflag = segment_format(fd, nrows, ncols, row_in_seg,
+ col_in_seg, size_struct);
+ else
+ errflag = segment_format_nofill(fd, nrows, ncols, row_in_seg,
+ col_in_seg, size_struct);
+
+ if (0 > errflag) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning(_("seg_open(): could not write segment file"));
+ return -1;
+ }
+ else {
+ G_warning(_("seg_open(): illegal configuration parameter(s)"));
+ return -3;
+ }
+ }
+ close(fd);
+ if (-1 == (fd = open(filename, 2))) {
+ unlink(filename);
+ G_warning(_("seg_open(): unable to re-open segment file"));
+ return -4;
+ }
+ if (0 > (errflag = segment_init(&(sseg->seg), fd, nsegs_in_memory))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning(_("seg_open(): could not read segment file"));
+ return -5;
+ }
+ else {
+ G_warning(_("seg_open(): out of memory"));
+ return -6;
+ }
+ }
+ sseg->filename = filename;
+ sseg->fd = fd;
+ return 0;
+}
+
+int seg_close(SSEG *sseg)
+{
+ segment_release(&(sseg->seg));
+ close(sseg->fd);
+ unlink(sseg->filename);
+ return 0;
+}
+
+int seg_put(SSEG *sseg, char *value, int row, int col)
+{
+ if (segment_put(&(sseg->seg), value, row, col) < 0) {
+ G_warning(_("seg_put(): could not write segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int seg_put_row(SSEG *sseg, char *value, int row)
+{
+ if (segment_put_row(&(sseg->seg), value, row) < 0) {
+ G_warning(_("seg_put_row(): could not write segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int seg_get(SSEG *sseg, char *value, int row, int col)
+{
+ if (segment_get(&(sseg->seg), value, row, col) < 0) {
+ G_warning(_("seg_get(): could not read segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int seg_get_row(SSEG *sseg, char *value, int row)
+{
+ if (segment_get_row(&(sseg->seg), value, row) < 0) {
+ G_warning(_("seg_get_row(): could not read segment file"));
+ return -1;
+ }
+ return 0;
+}
+
+int seg_flush(SSEG *sseg)
+{
+ segment_flush(&(sseg->seg));
+ return 0;
+}
Property changes on: grass-addons/grass7/raster/r.stream.extract/seg.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.stream.extract/seg.h
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/seg.h (rev 0)
+++ grass-addons/grass7/raster/r.stream.extract/seg.h 2010-12-16 14:01:19 UTC (rev 44615)
@@ -0,0 +1,78 @@
+#ifndef __SEG_H__
+#define __SEG_H__
+
+#include <grass/raster.h>
+#include <grass/segment.h>
+
+#define CSEG struct _c_s_e_g_
+CSEG {
+ SEGMENT seg; /* segment structure */
+ int fd; /* fd for reading/writing segment file */
+ char *filename; /* name of segment file */
+ char *name; /* raster map read into segment file */
+ char *mapset;
+};
+
+#define DSEG struct _d_s_e_g_
+DSEG {
+ SEGMENT seg; /* segment structure */
+ int fd; /* fd for reading/writing segment file */
+ char *filename; /* name of segment file */
+ char *name; /* raster map read into segment file */
+ char *mapset;
+};
+
+#define BSEG struct _b_s_e_g_
+BSEG {
+ SEGMENT seg; /* segment structure */
+ int fd; /* fd for reading/writing segment file */
+ char *filename; /* name of segment file */
+ char *name; /* raster map read into segment file */
+ char *mapset;
+};
+
+#define SSEG struct _s_s_e_g_
+SSEG {
+ SEGMENT seg; /* segment structure */
+ int fd; /* fd for reading/writing segment file */
+ char *filename; /* name of segment file */
+};
+
+/* bseg.c */
+int bseg_close(BSEG *);
+int bseg_get(BSEG *, char *, int, int);
+int bseg_open(BSEG *, int, int, int);
+int bseg_put(BSEG *, char *, int, int);
+int bseg_put_row(BSEG *, char *, int);
+int bseg_read_raster(BSEG *, char *, char *);
+int bseg_write_raster(BSEG *, char *);
+
+/* cseg.c */
+int cseg_close(CSEG *);
+int cseg_get(CSEG *, CELL *, int, int);
+int cseg_open(CSEG *, int, int, int);
+int cseg_put(CSEG *, CELL *, int, int);
+int cseg_put_row(CSEG *, CELL *, int);
+int cseg_read_raster(CSEG *, char *, char *);
+int cseg_write_raster(CSEG *, char *);
+
+/* dseg.c */
+int dseg_close(DSEG *);
+int dseg_get(DSEG *, double *, int, int);
+int dseg_open(DSEG *, int, int, int);
+int dseg_put(DSEG *, double *, int, int);
+int dseg_put_row(DSEG *, double *, int);
+int dseg_read_raster(DSEG *, char *, char *);
+int dseg_write_raster(DSEG *, char *);
+
+/* seg.c */
+int seg_close(SSEG *);
+int seg_get(SSEG *, char *, int, int);
+int seg_open(SSEG *, int, int, int, int, int, int, int);
+int seg_put(SSEG *, char *, int, int);
+int seg_put_row(SSEG *, char *, int);
+int seg_get(SSEG *, char *, int, int);
+int seg_get_row(SSEG *, char *, int);
+int seg_flush(SSEG *);
+
+#endif /* __SEG_H__ */
Property changes on: grass-addons/grass7/raster/r.stream.extract/seg.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/raster/r.stream.extract/streams.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/streams.c 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/streams.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -1,28 +1,9 @@
#include <stdlib.h>
#include <math.h>
-#include <grass/gis.h>
+#include <grass/raster.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;
@@ -39,55 +20,41 @@
}
int continue_stream(CELL stream_id, int r, int c, int r_max, int c_max,
- unsigned int thisindex, int *stream_no, int min_length)
+ int *stream_no)
{
char aspect;
- int curr_stream;
- int r_nbr, c_nbr, ct_dir;
+ CELL curr_stream, stream_nbr, old_stream;
+ int r_nbr, c_nbr;
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;
+ char this_flag_value;
G_debug(3, "continue stream");
+
+ cseg_get(&stream, &curr_stream, r_max, c_max);
- /* 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)] = stream_id;
+ G_debug(3, "no confluence, just continue stream");
+ curr_stream = stream_id;
+ cseg_put(&stream, &curr_stream, r_max, c_max);
+ bseg_get(&bitflags, &this_flag_value, r_max, c_max);
+ FLAG_SET(this_flag_value, STREAMFLAG);
+ bseg_put(&bitflags, &this_flag_value, r_max, c_max);
return 0;
}
- G_debug(2, "confluence");
- /* delete short stream segments */
- /* if (min_length && stream_node[stream_id].n_trib == 0) {
- if (seg_length(stream_id, NULL) < min_length) {
- del_stream_seg(stream_id);
- return 0;
- }
- }
- */
+ G_debug(3, "confluence");
/* new confluence */
if (stream_node[curr_stream].r != r_max ||
stream_node[curr_stream].c != c_max) {
- G_debug(2, "new confluence");
+ size_t new_size;
+
+ G_debug(3, "new confluence");
/* set new stream id */
- curr_stream = stream[INDEX(r_max, c_max)] = ++(*stream_no);
+ (*stream_no)++;
/* add stream node */
if (*stream_no >= n_alloc_nodes - 1) {
n_alloc_nodes += stream_node_step;
@@ -103,7 +70,6 @@
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 */
@@ -111,75 +77,43 @@
G_warning(_("BUG: 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[*stream_no].n_alloc += 2;
+ new_size = stream_node[*stream_no].n_alloc * sizeof(int);
+ stream_node[*stream_no].trib =
+ (int *)G_realloc(stream_node[*stream_no].trib, 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);
- }
+ /* add the two tributaries */
+ G_debug(3, "add tributaries");
+ stream_node[*stream_no].trib[stream_node[*stream_no].n_trib++] =
+ curr_stream;
+ stream_node[*stream_no].trib[stream_node[*stream_no].n_trib++] =
+ stream_id;
- 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");
+ G_debug(3, "update stream IDs downstream");
r_nbr = r_max;
c_nbr = c_max;
- draindir.pos = INDEX(r_nbr, c_nbr);
+ old_stream = curr_stream;
+ curr_stream = *stream_no;
+ cseg_put(&stream, &curr_stream, r_nbr, c_nbr);
+ bseg_get(&asp, &aspect, 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(_("BUG: 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(_("BUG: stream id not set"));
- else
- stream[INDEX(r_nbr, c_nbr)] = curr_stream;
+ while (aspect > 0) {
+ r_nbr = r_nbr + asp_r[(int)aspect];
+ c_nbr = c_nbr + asp_c[(int)aspect];
+ cseg_get(&stream, &stream_nbr, r_nbr, c_nbr);
+ if (stream_nbr != old_stream)
+ aspect = -1;
+ else {
+ cseg_put(&stream, &curr_stream, r_nbr, c_nbr);
+ bseg_get(&asp, &aspect, r_nbr, c_nbr);
+ }
}
}
else {
/* stream node already existing here */
- G_debug(2, "existing confluence");
+ G_debug(3, "existing confluence");
/* add new tributary to stream node */
if (stream_node[curr_stream].n_trib >=
stream_node[curr_stream].n_alloc) {
@@ -189,17 +123,13 @@
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];
+ stream_id;
}
- G_debug(2, "%d tribs", stream_node[curr_stream].n_trib);
+ G_debug(3, "%d tribs", stream_node[curr_stream].n_trib);
if (stream_node[curr_stream].n_trib == 1)
G_warning(_("BUG: stream node %d has only 1 tributary: %d"), curr_stream,
stream_node[curr_stream].trib[0]);
@@ -208,14 +138,13 @@
}
/*
- * extracts streams for threshold
+ * accumulate surface flow
*/
int do_accum(double d8cut)
{
int r, c, dr, dc;
- CELL ele_val, ele_nbr;
- DCELL value, valued;
- int count;
+ CELL ele_val, *ele_nbr;
+ DCELL value, *wat_nbr;
struct Cell_head window;
int mfd_cells, astar_not_set;
double *dist_to_nbr, *weight, sum_weight, max_weight;
@@ -229,15 +158,19 @@
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;
+ unsigned int workedon, killer, count;
+ char *flag_nbr, this_flag_value;
+ POINT astarpoint;
+ WAT_ALT wa;
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));
+ flag_nbr = (char *)G_malloc(sides * sizeof(char));
+ wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
+ ele_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
G_get_set_window(&window);
@@ -254,38 +187,43 @@
dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
}
- /* reset worked flag */
- flag_clear_all(worked);
-
/* distribute and accumulate */
- for (killer = 1; killer <= n_points; killer++) {
+ count = 0;
+ for (killer = 0; killer < n_points; killer++) {
+
G_percent(killer, n_points, 1);
+
+ seg_get(&astar_pts, (char *)&astarpoint, 0, killer);
+ r = astarpoint.r;
+ c = astarpoint.c;
- thisindex = astar_pts[killer];
- r = thisindex / ncols;
- c = thisindex - r * ncols;
- aspect = asp[thisindex];
+ bseg_get(&asp, &aspect, r, c);
+ /* do not distribute flow along edges or out of real depressions */
+ if (aspect <= 0) {
+ bseg_get(&bitflags, &this_flag_value, r, c);
+ FLAG_UNSET(this_flag_value, WORKEDFLAG);
+ bseg_put(&bitflags, &this_flag_value, r, c);
+ continue;
+ }
+
if (aspect) {
dr = r + asp_r[abs((int)aspect)];
dc = c + asp_c[abs((int)aspect)];
}
- else {
- dr = r;
- dc = c;
- /* can only happen with real depressions */
- if (!have_depressions)
- G_fatal_error(_("Bug in stream extraction"));
- FLAG_SET(worked, r, c);
- G_debug(1, "bottom of real depression");
- continue;
- }
r_max = dr;
c_max = dc;
- value = acc[thisindex];
+ seg_get(&watalt, (char *)&wa, r, c);
+ value = wa.wat;
+ /* WORKEDFLAG has been set during A* Search
+ * reversed meaning here: 0 = done, 1 = not yet done */
+ bseg_get(&bitflags, &this_flag_value, r, c);
+ FLAG_UNSET(this_flag_value, WORKEDFLAG);
+ bseg_put(&bitflags, &this_flag_value, r, c);
+
/***************************************/
/* get weights for flow distribution */
/***************************************/
@@ -295,7 +233,7 @@
np_side = -1;
mfd_cells = 0;
astar_not_set = 1;
- ele_val = ele[thisindex];
+ ele_val = wa.ele;
edge = 0;
/* this loop is needed to get the sum of weights */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -303,22 +241,30 @@
r_nbr = r + nextdr[ct_dir];
c_nbr = c + nextdc[ct_dir];
weight[ct_dir] = -1;
+ wat_nbr[ct_dir] = 0;
+ ele_nbr[ct_dir] = 0;
+
/* 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);
+ bseg_get(&bitflags, &flag_nbr[ct_dir], r_nbr, c_nbr);
+ if ((edge = FLAG_GET(flag_nbr[ct_dir], NULLFLAG)))
+ break;
+ seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
+ wat_nbr[ct_dir] = wa.wat;
+ ele_nbr[ct_dir] = wa.ele;
- is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ /* WORKEDFLAG has been set during A* Search
+ * reversed meaning here: 0 = done, 1 = not yet done */
+ is_worked = FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG) == 0;
if (is_worked == 0) {
- ele_nbr = ele[nindex];
- edge = G_is_c_null_value(&ele_nbr);
- if (!edge && ele_nbr <= ele_val) {
- if (ele_nbr < ele_val) {
+ if (ele_nbr[ct_dir] <= ele_val) {
+ if (ele_nbr[ct_dir] < ele_val) {
weight[ct_dir] =
mfd_pow((ele_val -
- ele_nbr) / dist_to_nbr[ct_dir]);
+ ele_nbr[ct_dir]) / dist_to_nbr[ct_dir]);
}
- if (ele_nbr == ele_val) {
+ if (ele_nbr[ct_dir] == ele_val) {
weight[ct_dir] =
mfd_pow(0.5 / dist_to_nbr[ct_dir]);
}
@@ -346,7 +292,6 @@
/* do not distribute flow along edges, this causes artifacts */
if (edge) {
G_debug(3, "edge");
- FLAG_SET(worked, r, c);
continue;
}
@@ -383,18 +328,16 @@
/* 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);
+ is_worked = FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG) == 0;
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;
+ wa.wat = wat_nbr[ct_dir] + value * weight[ct_dir];
+ wa.ele = ele_nbr[ct_dir];
+ seg_put(&watalt, (char *)&wa, r_nbr, c_nbr);
}
else if (ct_dir == np_side) {
/* check for consistency with A * path */
@@ -409,17 +352,18 @@
}
/* get out of depression in SFD mode */
else {
- nindex = INDEX(dr, dc);
- valued = acc[INDEX(dr, dc)];
- valued += value;
- acc[INDEX(dr, dc)] = valued;
+ wa.wat = wat_nbr[np_side] + value;
+ wa.ele = ele_nbr[np_side];
+ seg_put(&watalt, (char *)&wa, dr, dc);
}
-
- FLAG_SET(worked, r, c);
}
+ G_percent(1, 1, 2);
G_free(dist_to_nbr);
G_free(weight);
+ G_free(wat_nbr);
+ G_free(ele_nbr);
+ G_free(flag_nbr);
return 1;
}
@@ -427,14 +371,13 @@
/*
* extracts streams for threshold, accumulation is provided
*/
-int extract_streams(double threshold, double mont_exp, int min_length)
+int extract_streams(double threshold, double mont_exp, int min_length, int internal_acc)
{
int r, c, dr, dc;
- CELL is_swale, ele_val, ele_nbr;
- CELL *streamp;
- DCELL value, valued;
+ CELL is_swale, ele_val, *ele_nbr;
+ DCELL value, valued, *wat_nbr;
struct Cell_head window;
- int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
+ int mfd_cells, stream_cells, swale_cells, astar_not_set;
double *dist_to_nbr;
double dx, dy;
int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side, max_side;
@@ -452,14 +395,14 @@
* | 2 | | 3 |
* | 5 | 0 | 6 |
*/
- unsigned int thisindex, nindex, workedon, killer;
+ unsigned int workedon, killer, count;
int stream_no = 0, stream_node_step = 1000;
double slope, diag;
+ char *flag_nbr, this_flag_value;
+ POINT astarpoint;
+ WAT_ALT wa;
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;
@@ -470,11 +413,14 @@
/* init outlet nodes */
n_alloc_outlets = stream_node_step;
outlets =
- (struct point *)G_malloc(n_alloc_outlets * sizeof(struct point));
+ (POINT *)G_malloc(n_alloc_outlets * sizeof(POINT));
n_outlets = 0;
/* distances to neighbours */
dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
+ flag_nbr = (char *)G_malloc(sides * sizeof(char));
+ wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
+ ele_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
G_get_set_window(&window);
@@ -493,44 +439,40 @@
diag = sqrt(2);
- /* reset worked flag */
- flag_clear_all(worked);
-
- /* initialize streams */
- streamp = stream;
- for (r = 0; r < nrows; r++) {
- for (c = 0; c < ncols; c++) {
- *streamp = 0;
- streamp++;
- }
- }
-
workedon = 0;
/* extract streams */
- for (killer = 1; killer <= n_points; killer++) {
+ count = 0;
+ for (killer = 0; killer < n_points; killer++) {
G_percent(killer, n_points, 1);
+
+ seg_get(&astar_pts, (char *)&astarpoint, 0, killer);
+ r = astarpoint.r;
+ c = astarpoint.c;
- thisindex = astar_pts[killer];
- r = thisindex / ncols;
- c = thisindex - r * ncols;
- aspect = asp[thisindex];
+ bseg_get(&asp, &aspect, r, c);
- FLAG_SET(worked, r, c);
+ bseg_get(&bitflags, &this_flag_value, r, c);
+ /* internal acc: SET, external acc: UNSET */
+ if (internal_acc)
+ FLAG_SET(this_flag_value, WORKEDFLAG);
+ else
+ FLAG_UNSET(this_flag_value, WORKEDFLAG);
+ bseg_put(&bitflags, &this_flag_value, r, c);
/* do not distribute flow along edges */
if (aspect <= 0) {
G_debug(3, "edge");
- is_swale = stream[thisindex];
+ is_swale = FLAG_GET(this_flag_value, STREAMFLAG);
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,
+ (POINT *)G_realloc(outlets,
n_alloc_outlets *
- sizeof(struct point));
+ sizeof(POINT));
}
outlets[n_outlets].r = r;
outlets[n_outlets].c = c;
@@ -560,7 +502,8 @@
r_nbr = r_max = dr;
c_nbr = c_max = dc;
- value = acc[thisindex];
+ seg_get(&watalt, (char *)&wa, r, c);
+ value = wa.wat;
/**********************************/
/* find main drainage direction */
@@ -572,8 +515,7 @@
stream_cells = 0;
swale_cells = 0;
astar_not_set = 1;
- ele_val = ele[thisindex];
- is_null = 0;
+ ele_val = wa.ele;
edge = 0;
flat = 1;
/* find main drainage direction */
@@ -581,33 +523,43 @@
/* get r_nbr, c_nbr for neighbours */
r_nbr = r + nextdr[ct_dir];
c_nbr = c + nextdc[ct_dir];
+ wat_nbr[ct_dir] = 0;
+ ele_nbr[ct_dir] = 0;
+ flag_nbr[ct_dir] = 0;
+
/* 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);
+ bseg_get(&bitflags, &flag_nbr[ct_dir], r_nbr, c_nbr);
+ if ((edge = FLAG_GET(flag_nbr[ct_dir], NULLFLAG)))
+ break;
+ seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
+ wat_nbr[ct_dir] = wa.wat;
+ ele_nbr[ct_dir] = wa.ele;
/* check for swale cells */
- is_swale = stream[nindex];
- if (is_swale > 0)
+ is_swale = FLAG_GET(flag_nbr[ct_dir], STREAMFLAG);
+ if (is_swale)
swale_cells++;
/* check for stream cells */
- valued = fabs(acc[nindex]);
- ele_nbr = ele[nindex];
+ valued = fabs(wat_nbr[ct_dir]);
/* check all upstream neighbours */
if (valued >= threshold && ct_dir != np_side &&
- ele_nbr > ele_val)
+ ele_nbr[ct_dir] > ele_val)
stream_cells++;
- is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ is_worked = FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG);
+ if (!internal_acc)
+ is_worked = is_worked == 0;
+
if (is_worked == 0) {
- if (ele_nbr != ele_val)
+ if (ele_nbr[ct_dir] != ele_val)
flat = 0;
- edge = G_is_c_null_value(&ele_nbr);
- if (!edge && ele_nbr <= ele_val) {
+ if (ele_nbr[ct_dir] <= ele_val) {
mfd_cells++;
@@ -624,7 +576,7 @@
}
}
}
- else if (ct_dir == np_side) {
+ else if (ct_dir == np_side && !edge) {
/* check for consistency with A * path */
workedon++;
}
@@ -635,7 +587,7 @@
break;
}
- is_swale = stream[thisindex];
+ is_swale = FLAG_GET(this_flag_value, STREAMFLAG);
/* do not continue streams along edges, these are artifacts */
if (edge) {
@@ -646,26 +598,28 @@
if (n_outlets >= n_alloc_outlets) {
n_alloc_outlets += stream_node_step;
outlets =
- (struct point *)G_realloc(outlets,
+ (POINT *)G_realloc(outlets,
n_alloc_outlets *
- sizeof(struct point));
+ sizeof(POINT));
}
outlets[n_outlets].r = r;
outlets[n_outlets].c = c;
n_outlets++;
- if (asp[thisindex] > 0) {
+ if (aspect > 0) {
aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
- asp[thisindex] = aspect;
+ bseg_put(&asp, &aspect, r, c);
}
}
continue;
}
+ if (np_side < 0)
+ G_fatal_error("np_side < 0");
+
/* 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)) {
- max_acc = fabs(acc[nindex]);
+ if (fabs(wat_nbr[np_side] >= max_acc)) {
+ max_acc = fabs(wat_nbr[np_side]);
r_max = dr;
c_max = dc;
max_side = np_side;
@@ -673,17 +627,18 @@
}
if (mfd_cells == 0) {
flat = 0;
+ r_max = dr;
+ c_max = dc;
max_side = np_side;
}
/* update aspect */
/* r_max == r && c_max == c should not happen */
if ((r_max != dr || c_max != dc) && (r_max != r || c_max != c)) {
- asp[thisindex] = drain[r - r_max + 1][c - c_max + 1];
+ aspect = drain[r - r_max + 1][c - c_max + 1];
+ bseg_put(&asp, &aspect, r, c);
}
- is_swale = stream[thisindex];
-
/**********************/
/* start new stream */
/**********************/
@@ -695,10 +650,8 @@
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[max_side]) / ele_scale;
- slope = (double)(ele_val - ele_nbr) / ele_scale;
-
if (max_side > 3)
slope /= diag;
@@ -706,10 +659,13 @@
}
}
- if (is_swale < 1 && fabs(value) >= threshold && stream_cells < 1 &&
+ if (!is_swale && fabs(value) >= threshold && stream_cells < 1 &&
swale_cells < 1 && !flat) {
G_debug(2, "start new stream");
- is_swale = stream[thisindex] = ++stream_no;
+ is_swale = ++stream_no;
+ cseg_put(&stream, &is_swale, r, c);
+ FLAG_SET(this_flag_value, STREAMFLAG);
+ bseg_put(&bitflags, &this_flag_value, r, c);
/* add stream node */
if (stream_no >= n_alloc_nodes - 1) {
n_alloc_nodes += stream_node_step;
@@ -725,7 +681,6 @@
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 */
@@ -739,38 +694,41 @@
/*********************/
if (is_swale > 0) {
+ cseg_get(&stream, &is_swale, r, c);
if (r_max == r && c_max == c) {
/* can't continue stream, add outlet point
* r_max == r && c_max == c should not happen */
- G_debug(1, "can't continue stream at r %d c %d", r, c);
+ G_debug(0, "can't continue stream at r %d c %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));
+ (POINT *)G_malloc(n_alloc_outlets *
+ sizeof(POINT));
}
outlets[n_outlets].r = r;
outlets[n_outlets].c = c;
n_outlets++;
}
else {
- continue_stream(is_swale, r, c, r_max, c_max, thisindex,
- &stream_no, min_length);
+ continue_stream(is_swale, r, c, r_max, c_max,
+ &stream_no);
}
}
}
+ G_percent(1, 1, 2);
if (workedon)
- G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
+ G_warning(_("MFD: A * path already processed when setting drainage direction: %d of %d cells"),
workedon, n_points);
- flag_destroy(worked);
G_free(dist_to_nbr);
- G_free(astar_pts);
+ G_free(wat_nbr);
+ G_free(ele_nbr);
+ G_free(flag_nbr);
- G_debug(1, "%d outlets", n_outlets);
- G_debug(1, "%d nodes", n_stream_nodes);
- G_debug(1, "%d streams", stream_no);
+ G_debug(0, "%d outlets", n_outlets);
+ G_debug(0, "%d nodes", n_stream_nodes);
+ G_debug(0, "%d streams", stream_no);
return 1;
}
Property changes on: grass-addons/grass7/raster/r.stream.extract/streams.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/raster/r.stream.extract/thin.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/thin.c 2010-12-16 13:55:32 UTC (rev 44614)
+++ grass-addons/grass7/raster/r.stream.extract/thin.c 2010-12-16 14:01:19 UTC (rev 44615)
@@ -7,58 +7,48 @@
{
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;
+ CELL curr_stream, no_stream = 0;
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 };
+ char flag_value, aspect;
r = stream_node[stream_id].r;
c = stream_node[stream_id].c;
- thisindex = INDEX(r, c);
+ cseg_get(&stream, &curr_stream, r, c);
- curr_stream = stream[thisindex];
- if (curr_stream != stream_id)
- G_fatal_error(_("BUG: 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) {
+ bseg_get(&asp, &aspect, r, c);
+ if (aspect > 0) {
/* 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)];
+ last_r = r + asp_r[(int)aspect];
+ last_c = c + asp_c[(int)aspect];
+ cseg_get(&stream, &curr_stream, 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];
+ bseg_get(&asp, &aspect, last_r, last_c);
+ while (aspect > 0) {
+ r_nbr = last_r + asp_r[(int)aspect];
+ c_nbr = last_c + asp_c[(int)aspect];
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)
+ cseg_get(&stream, &curr_stream, r_nbr, c_nbr);
+ if (curr_stream != 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);
+ cseg_put(&stream, &no_stream, last_r, last_c);
+ bseg_get(&bitflags, &flag_value, last_r, last_c);
+ FLAG_UNSET(flag_value, STREAMFLAG);
+ bseg_put(&bitflags, &flag_value, last_r, last_c);
/* update start point */
- draindir.pos = thisindex;
- founddir = rbtree_find(draintree, &draindir);
- founddir->dir = drain[r - r_nbr + 1][c - c_nbr + 1];
- asp[draindir.pos] = founddir->dir;
- last_r = r_nbr;
- last_c = c_nbr;
- draindir.pos = INDEX(last_r, last_c);
+ aspect = drain[r - r_nbr + 1][c - c_nbr + 1];
+ bseg_put(&asp, &aspect, r, c);
thinned = 1;
}
@@ -66,11 +56,10 @@
/* 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);
}
+ last_r = r_nbr;
+ last_c = c_nbr;
+ bseg_get(&asp, &aspect, last_r, last_c);
}
}
@@ -80,8 +69,8 @@
int thin_streams(void)
{
int i, j, r, c, done;
- int stream_id, next_node;
- unsigned int thisindex;
+ CELL stream_id;
+ int next_node;
struct sstack
{
int stream_id;
@@ -89,6 +78,7 @@
} *nodestack;
int top = 0, stack_step = 1000;
int n_trib_total;
+ int n_thinned = 0;
G_message(_("Thin stream segments..."));
@@ -98,8 +88,7 @@
G_percent(i, n_outlets, 2);
r = outlets[i].r;
c = outlets[i].c;
- thisindex = INDEX(r, c);
- stream_id = stream[thisindex];
+ cseg_get(&stream, &stream_id, r, c);
if (stream_id == 0)
continue;
@@ -147,8 +136,10 @@
if (thin_seg(stream_id) == 0)
G_debug(3, "segment %d not thinned", stream_id);
- else
+ else {
G_debug(3, "segment %d thinned", stream_id);
+ n_thinned++;
+ }
top--;
/* count tributaries */
@@ -174,6 +165,8 @@
G_percent(n_outlets, n_outlets, 1); /* finish it */
G_free(nodestack);
+
+ G_verbose_message(_("%d of %d stream segments were thinned"), n_thinned, n_stream_nodes);
return 1;
}
Property changes on: grass-addons/grass7/raster/r.stream.extract/thin.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list