[GRASS-SVN] r50109 - grass-addons/grass7/raster/r.hydrodem
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jan 9 12:41:57 EST 2012
Author: mmetz
Date: 2012-01-09 09:41:56 -0800 (Mon, 09 Jan 2012)
New Revision: 50109
Added:
grass-addons/grass7/raster/r.hydrodem/hydro_con.c
grass-addons/grass7/raster/r.hydrodem/r.hydrodem.html
Removed:
grass-addons/grass7/raster/r.hydrodem/del_streams.c
grass-addons/grass7/raster/r.hydrodem/r.stream.extract.html
grass-addons/grass7/raster/r.hydrodem/streams.c
grass-addons/grass7/raster/r.hydrodem/thin.c
Modified:
grass-addons/grass7/raster/r.hydrodem/Makefile
grass-addons/grass7/raster/r.hydrodem/close.c
grass-addons/grass7/raster/r.hydrodem/do_astar.c
grass-addons/grass7/raster/r.hydrodem/init_search.c
grass-addons/grass7/raster/r.hydrodem/load.c
grass-addons/grass7/raster/r.hydrodem/local_proto.h
grass-addons/grass7/raster/r.hydrodem/main.c
Log:
new module for sink removal
Modified: grass-addons/grass7/raster/r.hydrodem/Makefile
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/Makefile 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/Makefile 2012-01-09 17:41:56 UTC (rev 50109)
@@ -1,11 +1,9 @@
MODULE_TOPDIR = ../..
-PGM = r.stream.extract
+PGM = r.hydrodem
-LIBES = $(SEGMENTLIB) $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB)
-DEPENDENCIES = $(SEGMENTDEP) $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
-EXTRA_INC = $(VECT_INC)
-EXTRA_CFLAGS = $(VECT_CFLAGS)
+LIBES = $(SEGMENTLIB) $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(SEGMENTDEP) $(RASTERDEP) $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Property changes on: grass-addons/grass7/raster/r.hydrodem/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/raster/r.hydrodem/close.c
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/close.c 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/close.c 2012-01-09 17:41:56 UTC (rev 50109)
@@ -1,322 +1,57 @@
+#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/glocale.h>
-#include <grass/vector.h>
-#include <grass/dbmi.h>
#include "local_proto.h"
-int close_streamvect(char *stream_vect)
+int close_map(char *ele_rast, int ele_map_type)
{
- int i, r, c, r_nbr, c_nbr, done;
- CELL stream_id, stream_nbr;
- char aspect;
- int next_node;
- struct sstack
- {
- int stream_id;
- int next_trib;
- } *nodestack;
- int top = 0, stack_step = 1000;
- int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
- int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
- struct Map_info Out;
- static struct line_pnts *Points;
- struct line_cats *Cats;
- dbDriver *driver;
- dbHandle handle;
- dbString table_name, dbsql, valstr;
- struct field_info *Fi;
- char *cat_col_name = "cat", buf[2000];
- struct Cell_head window;
- double north_offset, west_offset, ns_res, ew_res;
- int next_cat;
+ int ele_fd, r, c;
+ void *ele_buf, *ele_ptr;
+ struct History history;
+ size_t ele_size;
+ CELL ele_val;
- G_message(_("Write vector map <%s>..."), stream_vect);
+ G_message(_("Write conditioned DEM"));
- if (0 > Vect_open_new(&Out, stream_vect, 0)) {
- G_fatal_error(_("Unable to create vector map <%s>"), stream_vect);
- }
+ ele_fd = Rast_open_new(ele_rast, ele_map_type);
+ ele_size = Rast_cell_size(ele_map_type);
+ ele_buf = Rast_allocate_buf(ele_map_type);
- nodestack = (struct sstack *)G_malloc(stack_step * sizeof(struct sstack));
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 2);
- Points = Vect_new_line_struct();
- Cats = Vect_new_cats_struct();
+ Rast_set_null_value(ele_buf, ncols, ele_map_type); /* reset row to all NULL */
- G_get_set_window(&window);
- ns_res = window.ns_res;
- ew_res = window.ew_res;
- north_offset = window.north - 0.5 * ns_res;
- west_offset = window.west + 0.5 * ew_res;
+ ele_ptr = ele_buf;
- next_cat = n_stream_nodes + 1;
+ for (c = 0; c < ncols; c++) {
+ cseg_get(&ele, &ele_val, r, c);
- for (i = 0; i < n_outlets; i++, next_cat++) {
- G_percent(i, n_outlets, 2);
- r = outlets[i].r;
- c = outlets[i].c;
- cseg_get(&stream, &stream_id, r, c);
+ if (!Rast_is_c_null_value(&ele_val)) {
- if (!stream_id)
- continue;
-
- Vect_reset_line(Points);
- Vect_reset_cats(Cats);
-
- /* outlet */
- Vect_cat_set(Cats, 1, stream_id);
- Vect_cat_set(Cats, 2, 2);
- Vect_append_point(Points, west_offset + c * ew_res,
- north_offset - r * ns_res, 0);
- Vect_write_line(&Out, GV_POINT, Points, Cats);
-
- /* add root node to stack */
- G_debug(3, "add root node");
- top = 0;
- nodestack[top].stream_id = stream_id;
- nodestack[top].next_trib = 0;
-
- /* depth first post order traversal */
- G_debug(3, "traverse");
- while (top >= 0) {
-
- done = 1;
- stream_id = nodestack[top].stream_id;
- G_debug(3, "stream_id %d", stream_id);
- if (nodestack[top].next_trib < stream_node[stream_id].n_trib) {
- /* add to stack */
- next_node =
- stream_node[stream_id].trib[nodestack[top].next_trib];
- G_debug(3, "add to stack: next %d, trib %d, n trib %d",
- next_node, nodestack[top].next_trib,
- stream_node[stream_id].n_trib);
- nodestack[top].next_trib++;
- top++;
- if (top >= stack_step) {
- /* need more space */
- stack_step += 1000;
- nodestack =
- (struct sstack *)G_realloc(nodestack,
- stack_step *
- sizeof(struct sstack));
+ switch (ele_map_type) {
+ case CELL_TYPE:
+ *((CELL *) ele_ptr) = ele_val;
+ break;
+ case FCELL_TYPE:
+ *((FCELL *) ele_ptr) = (FCELL) ele_val / ele_scale;
+ break;
+ case DCELL_TYPE:
+ *((DCELL *) ele_ptr) = (DCELL) ele_val / ele_scale;
+ break;
}
- nodestack[top].next_trib = 0;
- nodestack[top].stream_id = next_node;
- done = 0;
- G_debug(3, "go further down");
}
- if (done) {
- G_debug(3, "write stream segment");
-
- Vect_reset_line(Points);
- Vect_reset_cats(Cats);
-
- r_nbr = stream_node[stream_id].r;
- c_nbr = stream_node[stream_id].c;
-
- 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);
- else
- Vect_cat_set(Cats, 2, 1);
-
- Vect_append_point(Points, west_offset + c_nbr * ew_res,
- north_offset - r_nbr * ns_res, 0);
-
- Vect_write_line(&Out, GV_POINT, Points, Cats);
-
- 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");
-
- 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;
- }
- bseg_get(&asp, &aspect, r_nbr, c_nbr);
- }
-
- Vect_write_line(&Out, GV_LINE, Points, Cats);
-
- top--;
- }
+ ele_ptr = G_incr_void_ptr(ele_ptr, ele_size);
}
+ Rast_put_row(ele_fd, ele_buf, ele_map_type);
}
- G_percent(n_outlets, n_outlets, 1); /* finish it */
-
- G_message(_("Write vector attribute table"));
-
- /* Prepeare strings for use in db_* calls */
- db_init_string(&dbsql);
- db_init_string(&valstr);
- db_init_string(&table_name);
- db_init_handle(&handle);
-
- /* Preparing database for use */
- /* Create database for new vector map */
- Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
- driver = db_start_driver_open_database(Fi->driver,
- Vect_subst_var(Fi->database,
- &Out));
- if (driver == NULL) {
- G_fatal_error(_("Unable to start driver <%s>"), Fi->driver);
- }
-
- G_debug(1, "table: %s", Fi->table);
- G_debug(1, "driver: %s", Fi->driver);
- G_debug(1, "database: %s", Fi->database);
-
- sprintf(buf,
- "create table %s (%s integer, stream_type varchar(20), type_code integer)",
- Fi->table, cat_col_name);
- db_set_string(&dbsql, buf);
-
- if (db_execute_immediate(driver, &dbsql) != DB_OK) {
- db_close_database(driver);
- db_shutdown_driver(driver);
- G_fatal_error(_("Cannot create table: %s"), db_get_string(&dbsql));
- }
-
- if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
- G_warning(_("Cannot create index"));
-
- if (db_grant_on_table(driver, Fi->table, DB_PRIV_SELECT,
- DB_GROUP | DB_PUBLIC) != DB_OK)
- G_fatal_error(_("Cannot grant privileges on table %s"), Fi->table);
-
- db_begin_transaction(driver);
-
- /* stream nodes */
- for (i = 1; i <= n_stream_nodes; i++) {
-
- sprintf(buf, "insert into %s values ( %d, \'%s\', %d )",
- Fi->table, i,
- (stream_node[i].n_trib > 0 ? "intermediate" : "start"),
- (stream_node[i].n_trib > 0));
-
- db_set_string(&dbsql, buf);
-
- if (db_execute_immediate(driver, &dbsql) != DB_OK) {
- db_close_database(driver);
- db_shutdown_driver(driver);
- G_fatal_error(_("Cannot insert new row: %s"),
- db_get_string(&dbsql));
- }
- }
-
- db_commit_transaction(driver);
- db_close_database_shutdown_driver(driver);
-
- Vect_map_add_dblink(&Out, 1, NULL, Fi->table,
- cat_col_name, Fi->database, Fi->driver);
-
- G_debug(1, "close vector");
-
- Vect_hist_command(&Out);
- Vect_build(&Out);
- Vect_close(&Out);
-
- G_free(nodestack);
-
- return 1;
-}
-
-
-int close_maps(char *stream_rast, char *stream_vect, char *dir_rast)
-{
- int stream_fd, dir_fd, r, c, i;
- CELL *cell_buf1, *cell_buf2;
- struct History history;
- char flag_value;
- CELL stream_id;
- char aspect;
-
- /* cheating... */
- stream_fd = dir_fd = -1;
- cell_buf1 = cell_buf2 = NULL;
-
- G_message(_("Write raster %s"),
- (stream_rast != NULL) + (dir_rast != NULL) > 1 ? "maps" : "map");
-
- /* write requested output rasters */
- if (stream_rast) {
- stream_fd = Rast_open_new(stream_rast, CELL_TYPE);
- cell_buf1 = Rast_allocate_c_buf();
- }
- if (dir_rast) {
- 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)
- Rast_set_c_null_value(cell_buf1, ncols); /* reset row to all NULL */
- if (dir_rast)
- Rast_set_c_null_value(cell_buf2, ncols); /* reset row to all NULL */
-
- for (c = 0; c < ncols; c++) {
- if (stream_rast) {
- cseg_get(&stream, &stream_id, r, c);
- if (stream_id)
- cell_buf1[c] = stream_id;
- }
- if (dir_rast) {
- 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)
- Rast_put_row(stream_fd, cell_buf1, CELL_TYPE);
- if (dir_rast)
- Rast_put_row(dir_fd, cell_buf2, CELL_TYPE);
- }
G_percent(nrows, nrows, 2); /* finish it */
- if (stream_rast) {
- Rast_close(stream_fd);
- G_free(cell_buf1);
- Rast_short_history(stream_rast, "raster", &history);
- Rast_command_history(&history);
- Rast_write_history(stream_rast, &history);
- }
- if (dir_rast) {
- Rast_close(dir_fd);
- G_free(cell_buf2);
- Rast_short_history(dir_rast, "raster", &history);
- Rast_command_history(&history);
- Rast_write_history(dir_rast, &history);
- }
+ Rast_close(ele_fd);
+ G_free(ele_buf);
+ Rast_short_history(ele_rast, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(ele_rast, &history);
- /* close stream vector */
- if (stream_vect) {
- if (close_streamvect(stream_vect) < 0)
- G_fatal_error(_("Unable to write vector map <%s>"), stream_vect);
- }
-
- /* rearranging desk chairs on the Titanic... */
- 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);
-
return 1;
}
Deleted: grass-addons/grass7/raster/r.hydrodem/del_streams.c
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/del_streams.c 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/del_streams.c 2012-01-09 17:41:56 UTC (rev 50109)
@@ -1,194 +0,0 @@
-#include <stdlib.h>
-#include <math.h>
-#include <grass/raster.h>
-#include <grass/glocale.h>
-#include "local_proto.h"
-
-/* get stream segment length */
-int seg_length(CELL stream_id, CELL *next_stream_id)
-{
- int r, c, r_nbr, c_nbr;
- int slength = 1;
- 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 */
- 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)
- break;
- /* outside region */
- if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
- break;
- /* next stream */
- 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;
- bseg_get(&asp, &aspect, r, c);
- }
-
- return slength;
-}
-
-/* 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;
- 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;
- 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 */
- 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)
- break;
- /* outside region */
- if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
- break;
- /* next stream */
- cseg_get(&stream, &curr_stream, r_nbr, c_nbr);
- if (curr_stream != stream_id)
- break;
- r = r_nbr;
- c = c_nbr;
- 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) {
- 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 curr_stream;
-}
-
-/* delete stream segments shorter than threshold */
-int del_streams(int min_length)
-{
- int i;
- int n_deleted = 0;
- CELL curr_stream, stream_id;
- int other_trib, tmp_trib;
- int slength;
-
- G_message(_("Delete stream segments shorter than %d cells..."), min_length);
-
- /* TODO: proceed from stream heads to outlets
- * -> use depth first post order traversal */
-
- /* 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 */
- 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++;
- }
- }
-
- G_verbose_message(_("%d stream segments deleted"), n_deleted);
-
- return n_deleted;
-}
Modified: grass-addons/grass7/raster/r.hydrodem/do_astar.c
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/do_astar.c 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/do_astar.c 2012-01-09 17:41:56 UTC (rev 50109)
@@ -1,43 +1,53 @@
#include <stdlib.h>
#include <math.h>
-#include <grass/raster.h>
+#include <grass/gis.h>
#include <grass/glocale.h>
#include "local_proto.h"
#define GET_PARENT(c) ((unsigned int)(((c) - 2) >> 2) + 1)
#define GET_CHILD(p) ((unsigned int)((p) << 2) - 2)
-HEAP_PNT heap_drop(void);
-int sift_up(unsigned int, HEAP_PNT);
+#define HEAP_CMP(a, b) (((a)->ele < (b)->ele) ? 1 : \
+ (((a)->ele > (b)->ele) ? 0 : \
+ (((a)->added < (b)->added) ? 1 : 0)))
+
+
+struct heap_point heap_drop(void);
+
double get_slope(CELL, CELL, double);
int do_astar(void)
{
int r, c, r_nbr, c_nbr, ct_dir;
- unsigned int first_cum, count;
+ int 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];
- WAT_ALT wa;
- char asp_val;
+ 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 };
+ CELL ele_val, ele_down, ele_nbr[8];
+ char is_bottom;
+ char asp_val, asp_val_this;
char flag_value, is_in_list, is_worked;
- HEAP_PNT heap_p;
+ struct heap_point heap_p;
/* sides
* |7|1|4|
* |2| |3|
* |5|0|6|
*/
- int nbr_ew[8] = { 0, 1, 2, 3, 1, 0, 0, 1 };
- int nbr_ns[8] = { 0, 1, 2, 3, 3, 2, 3, 2 };
+ int nbr_ew[8] = { 0, 1, 2, 3, 1, 0, 0, 1};
+ int nbr_ns[8] = { 0, 1, 2, 3, 3, 2, 3, 2};
double dx, dy, dist_to_nbr[8], ew_res, ns_res;
double slope[8];
+ int skip_diag;
struct Cell_head window;
- int skip_me;
count = 0;
first_cum = n_points;
+ sinks = first_sink = NULL;
+ n_sinks = 0;
+
G_message(_("A* Search..."));
Rast_get_window(&window);
@@ -56,7 +66,7 @@
}
ew_res = window.ew_res;
ns_res = window.ns_res;
-
+
while (heap_size > 0) {
G_percent(count++, n_points, 1);
if (count > n_points)
@@ -70,17 +80,31 @@
heap_p = heap_drop();
- r = heap_p.pnt.r;
- c = heap_p.pnt.c;
+ /* flow accumulation order is not needed */
+ r = heap_p.r;
+ c = heap_p.c;
- ele_val = heap_p.ele;
+ first_cum--;
+ bseg_get(&bitflags, &flag_value, r, c);
+ FLAG_SET(flag_value, WORKEDFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
+ ele_val = ele_down = heap_p.ele;
+ is_bottom = 0;
+ bseg_get(&draindir, &asp_val_this, r, c);
+ if (asp_val_this > 0 && !FLAG_GET(flag_value, EDGEFLAG)) {
+ r_nbr = r + asp_r[(int)asp_val_this];
+ c_nbr = c + asp_c[(int)asp_val_this];
+ cseg_get(&ele, &ele_down, r_nbr, c_nbr);
+ if (ele_down > ele_val)
+ is_bottom = 1;
+ }
+
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];
slope[ct_dir] = ele_nbr[ct_dir] = 0;
- skip_me = 0;
/* check that neighbour is within region */
if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
@@ -89,73 +113,105 @@
bseg_get(&bitflags, &flag_value, r_nbr, c_nbr);
is_in_list = FLAG_GET(flag_value, INLISTFLAG);
is_worked = FLAG_GET(flag_value, WORKEDFLAG);
+ skip_diag = 0;
+
+ /* avoid diagonal flow direction bias */
if (!is_worked) {
- seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
- ele_nbr[ct_dir] = wa.ele;
+ cseg_get(&ele, &ele_nbr[ct_dir], r_nbr, c_nbr);
slope[ct_dir] = get_slope(ele_val, ele_nbr[ct_dir],
- dist_to_nbr[ct_dir]);
+ dist_to_nbr[ct_dir]);
}
- /* avoid diagonal flow direction bias */
+
if (!is_in_list) {
if (ct_dir > 3 && slope[ct_dir] > 0) {
if (slope[nbr_ew[ct_dir]] > 0) {
/* slope to ew nbr > slope to center */
- if (slope[ct_dir] <
- get_slope(ele_nbr[nbr_ew[ct_dir]],
- ele_nbr[ct_dir], ew_res))
- skip_me = 1;
+ if (slope[ct_dir] < get_slope(ele_nbr[nbr_ew[ct_dir]],
+ ele_nbr[ct_dir], ew_res))
+ skip_diag = 1;
}
- if (!skip_me && slope[nbr_ns[ct_dir]] > 0) {
+ if (!skip_diag && slope[nbr_ns[ct_dir]] > 0) {
/* slope to ns nbr > slope to center */
- if (slope[ct_dir] <
- get_slope(ele_nbr[nbr_ns[ct_dir]],
- ele_nbr[ct_dir], ns_res))
- skip_me = 1;
+ if (slope[ct_dir] < get_slope(ele_nbr[nbr_ns[ct_dir]],
+ ele_nbr[ct_dir], ns_res))
+ skip_diag = 1;
}
}
}
- if (is_in_list == 0 && skip_me == 0) {
- ele_up = ele_nbr[ct_dir];
+ if (is_in_list == 0 && skip_diag == 0) {
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);
+ heap_add(r_nbr, c_nbr, ele_nbr[ct_dir], asp_val, flag_value);
+
+ if (ele_nbr[ct_dir] < ele_val)
+ is_bottom = 0;
}
else if (is_in_list && is_worked == 0) {
if (FLAG_GET(flag_value, EDGEFLAG)) {
- /* neighbour is edge in list, not yet worked */
- bseg_get(&asp, &asp_val, r_nbr, c_nbr);
+ bseg_get(&draindir, &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);
+ /* update edge cell ? no */
+ /* asp_val = drain[r_nbr - r + 1][c_nbr - c + 1]; */
+ /* check if this causes trouble */
+ bseg_put(&draindir, &asp_val, r_nbr, c_nbr);
+
+ if (ele_nbr[ct_dir] < ele_val)
+ is_bottom = 0;
}
}
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);
+ bseg_get(&draindir, &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);
+ bseg_put(&draindir, &asp_val, r_nbr, c_nbr);
FLAG_UNSET(flag_value, DEPRFLAG);
bseg_put(&bitflags, &flag_value, r_nbr, c_nbr);
}
}
}
} /* end neighbours */
- /* 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);
+
+
+ if (is_bottom) {
+ /* add sink bottom */
+ /* process upstream */
+ if (1) {
+ if (first_sink) {
+ sinks->next =
+ (struct sink_list *)
+ G_malloc(sizeof(struct sink_list));
+ sinks = sinks->next;
+ }
+ /* first sink */
+ else {
+ first_sink =
+ (struct sink_list *)
+ G_malloc(sizeof(struct sink_list));
+ sinks = first_sink;
+ }
+ sinks->next = NULL;
+ }
+ /* process downstream */
+ if (0) {
+ sinks =
+ (struct sink_list *)
+ G_malloc(sizeof(struct sink_list));
+ sinks->next = first_sink;
+ first_sink = sinks;
+ }
+ sinks->r = r;
+ sinks->c = c;
+ n_sinks++;
+ }
} /* end A* search */
G_percent(n_points, n_points, 1); /* finish it */
+ if (first_cum)
+ G_warning(_("processed points mismatch of %u"), first_cum);
+
return 1;
}
@@ -163,7 +219,8 @@
* compare function for heap
* returns 1 if point1 < point2 else 0
*/
-int heap_cmp(HEAP_PNT *a, HEAP_PNT *b)
+
+static int heap_cmp(struct heap_point *a, struct heap_point *b)
{
if (a->ele < b->ele)
return 1;
@@ -174,15 +231,16 @@
return 0;
}
-int sift_up(unsigned int start, HEAP_PNT child_p)
+int sift_up(unsigned int start, struct heap_point child_p)
{
unsigned int parent, child;
- HEAP_PNT heap_p;
+ struct heap_point heap_p;
child = start;
while (child > 1) {
parent = GET_PARENT(child);
+
seg_get(&search_heap, (char *)&heap_p, 0, parent);
/* push parent point down if child is smaller */
@@ -191,7 +249,7 @@
child = parent;
}
else
- /* no more sifting up, found slot for child */
+ /* no more sifting up, found new slot for child */
break;
}
@@ -205,22 +263,30 @@
* add item to heap
* returns heap_size
*/
-unsigned int heap_add(int r, int c, CELL ele)
+unsigned int heap_add(int r, int c, CELL ele, char asp, char flag_value)
{
- HEAP_PNT heap_p;
-
+ struct heap_point heap_p;
+
/* add point to next free position */
heap_size++;
+ if (heap_size > n_points)
+ G_fatal_error(_("Heapsize too large"));
+
+ heap_p.r = r;
+ heap_p.c = c;
+ heap_p.ele = ele;
heap_p.added = nxt_avail_pt;
- heap_p.ele = ele;
- heap_p.pnt.r = r;
- heap_p.pnt.c = c;
+ bseg_put(&draindir, &asp, r, c);
+ FLAG_SET(flag_value, INLISTFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
+
nxt_avail_pt++;
/* sift up: move new point towards top of heap */
+
sift_up(heap_size, heap_p);
return heap_size;
@@ -230,11 +296,11 @@
* drop item from heap
* returns heap size
*/
-HEAP_PNT heap_drop(void)
+struct heap_point heap_drop(void)
{
unsigned int child, childr, parent;
int i;
- HEAP_PNT child_p, childr_p, last_p, root_p;
+ struct heap_point 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);
@@ -245,14 +311,14 @@
}
parent = 1;
- while ((child = GET_CHILD(parent)) < heap_size) {
+ while ((child = GET_CHILD(parent)) <= heap_size) {
seg_get(&search_heap, (char *)&child_p, 0, child);
if (child < heap_size) {
childr = child + 1;
i = child + 4;
- while (childr < heap_size && childr < i) {
+ while (childr <= heap_size && childr < i) {
seg_get(&search_heap, (char *)&childr_p, 0, childr);
if (heap_cmp(&childr_p, &child_p)) {
child = childr;
@@ -261,7 +327,6 @@
childr++;
}
}
-
if (heap_cmp(&last_p, &child_p)) {
break;
}
Added: grass-addons/grass7/raster/r.hydrodem/hydro_con.c
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/hydro_con.c (rev 0)
+++ grass-addons/grass7/raster/r.hydrodem/hydro_con.c 2012-01-09 17:41:56 UTC (rev 50109)
@@ -0,0 +1,859 @@
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+struct ns
+{
+ int r, c;
+ int next_side;
+ int n_ngbrs;
+};
+
+static struct ns *n_stack;
+static unsigned int stack_alloc;
+
+int count_fill(int peak_r, int peak_c, CELL peak_ele, int *n_splits, CELL *next_ele)
+{
+ int r, c, r_nbr, c_nbr, ct_dir;
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ 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 };
+ CELL ele_val, ele_nbr;
+ int n_to_fill;
+ int top, done;
+ char drain_val;
+
+ /* go upstream from spill point */
+ r = peak_r;
+ c = peak_c;
+ ele_val = peak_ele;
+ *n_splits = 0;
+ n_to_fill = 0;
+ /* post-order traversal */
+ /* add spill point as root to stack */
+ top = 0;
+ n_stack[top].r = peak_r;
+ n_stack[top].c = peak_c;
+ n_stack[top].next_side = 0;
+ n_stack[top].n_ngbrs = 0;
+ while (top >= 0) {
+ done = 1;
+ r = n_stack[top].r;
+ c = n_stack[top].c;
+ for (ct_dir = n_stack[top].next_side; ct_dir < sides; ct_dir++) {
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 && c_nbr < ncols) {
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+ if (drain_val > 0) {
+ /* contributing cell */
+ if (r_nbr + asp_r[(int)drain_val] == r &&
+ c_nbr + asp_c[(int)drain_val] == c) {
+
+ /* contributing neighbour is lower than spill point, add to stack */
+ if (peak_ele >= ele_nbr) {
+
+ if (peak_ele > ele_nbr)
+ n_to_fill++;
+ n_stack[top].next_side = ct_dir + 1;
+ n_stack[top].n_ngbrs++;
+ top++;
+ if (top >= stack_alloc) {
+ stack_alloc += nrows + ncols;
+ n_stack =
+ (struct ns *)G_realloc(n_stack,
+ stack_alloc *
+ sizeof(struct ns));
+ }
+ n_stack[top].r = r_nbr;
+ n_stack[top].c = c_nbr;
+ n_stack[top].next_side = 0;
+ n_stack[top].n_ngbrs = 0;
+ done = 0;
+ break;
+ }
+ else if (next_ele && peak_ele < ele_nbr) {
+ if (*next_ele > ele_nbr)
+ *next_ele = ele_nbr;
+ }
+ }
+ } /* end contributing */
+ } /* end in region */
+ } /* end sides */
+
+ if (done) {
+ *n_splits += (n_stack[top].n_ngbrs > 1);
+ n_stack[top].next_side = sides;
+
+ top--;
+ }
+ }
+ return n_to_fill;
+}
+
+/* detect flat area:
+ * count neighbouring cells with same ele
+ * exclude contributing cells and cell this one is contributing to */
+int is_flat(int r, int c, CELL this_ele)
+{
+ int r_nbr, c_nbr, ct_dir;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ CELL ele_nbr;
+ char drain_val, drain_val_nbr;
+ int counter = 0;
+
+ bseg_get(&draindir, &drain_val, r, c);
+
+ if (drain_val < 1)
+ return 0;
+
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 && c_nbr < ncols) {
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ bseg_get(&draindir, &drain_val_nbr, r_nbr, c_nbr);
+ if (drain_val_nbr > 0) {
+ /* not a contributing cell */
+ if (r_nbr + asp_r[(int)drain_val_nbr] != r ||
+ c_nbr + asp_c[(int)drain_val_nbr] != c) {
+ /* does not contribute to this cell */
+ if (r + asp_r[(int)drain_val] != r_nbr ||
+ c + asp_c[(int)drain_val] != c_nbr) {
+ if (ele_nbr == this_ele)
+ counter++;
+ }
+ }
+ }
+ }
+ }
+
+ return (counter > 0);
+}
+
+int fill_sink(int peak_r, int peak_c, CELL peak_ele)
+{
+ int r, c, r_nbr, c_nbr, ct_dir;
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ 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 };
+ CELL ele_nbr;
+ int n_to_fill = 0;
+ int top, done;
+ char drain_val;
+
+ /* post-order traversal */
+ /* add spill point as root to stack */
+ top = 0;
+ n_stack[top].r = peak_r;
+ n_stack[top].c = peak_c;
+ n_stack[top].next_side = 0;
+ n_stack[top].n_ngbrs = 0;
+ while (top >= 0) {
+ done = 1;
+ r = n_stack[top].r;
+ c = n_stack[top].c;
+ for (ct_dir = n_stack[top].next_side; ct_dir < sides; ct_dir++) {
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 && c_nbr < ncols) {
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+ if (drain_val > 0) {
+ /* contributing cell */
+ if (r_nbr + asp_r[(int)drain_val] == r &&
+ c_nbr + asp_c[(int)drain_val] == c) {
+
+ /* contributing neighbour is lower than spill point, add to stack */
+ if (peak_ele >= ele_nbr) {
+
+ n_to_fill++;
+ n_stack[top].next_side = ct_dir + 1;
+ n_stack[top].n_ngbrs++;
+ top++;
+ if (top >= stack_alloc) {
+ stack_alloc += nrows + ncols;
+ n_stack =
+ (struct ns *)G_realloc(n_stack,
+ stack_alloc *
+ sizeof(struct ns));
+ }
+ n_stack[top].r = r_nbr;
+ n_stack[top].c = c_nbr;
+ n_stack[top].next_side = 0;
+ n_stack[top].n_ngbrs = 0;
+ done = 0;
+ break;
+ }
+ }
+ } /* end contributing */
+ } /* end in region */
+ } /* end sides */
+
+ if (done) {
+ cseg_put(&ele, &peak_ele, r, c);
+ n_stack[top].next_side = sides;
+ top--;
+ }
+ }
+
+ return n_to_fill;
+}
+
+/* carve channel */
+int carve(int bottom_r, int bottom_c, CELL bottom_ele, int peak_r, int peak_c)
+{
+ int r, c, r_nbr, c_nbr, carved = 0;
+ char drain_val;
+ int ct_dir;
+ CELL ele_val, ele_nbr;
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ 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 top, done;
+
+ /* carve upstream from spill point --> */
+ /* post-order traversal */
+ /* add spill point as root to stack */
+ top = 0;
+ n_stack[top].r = peak_r;
+ n_stack[top].c = peak_c;
+ n_stack[top].next_side = 0;
+ n_stack[top].n_ngbrs = 0;
+ while (top >= 0) {
+ done = 1;
+ r = n_stack[top].r;
+ c = n_stack[top].c;
+ cseg_get(&ele, &ele_val, r, c);
+ for (ct_dir = n_stack[top].next_side; ct_dir < sides; ct_dir++) {
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 && c_nbr < ncols) {
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+ if (drain_val > 0) {
+ /* contributing cell */
+ if (r_nbr + asp_r[(int) drain_val] == r &&
+ c_nbr + asp_c[(int) drain_val] == c) {
+
+ /* contributing neighbour is lower than current point, add to stack */
+ if (ele_val > ele_nbr && ele_nbr >= bottom_ele) {
+
+ n_stack[top].next_side = ct_dir + 1;
+ n_stack[top].n_ngbrs++;
+ top++;
+ if (top >= stack_alloc) {
+ stack_alloc += nrows + ncols;
+ n_stack =
+ (struct ns *)G_realloc(n_stack,
+ stack_alloc *
+ sizeof(struct ns));
+ }
+ n_stack[top].r = r_nbr;
+ n_stack[top].c = c_nbr;
+ n_stack[top].next_side = 0;
+ n_stack[top].n_ngbrs = 0;
+ done = 0;
+ break;
+ }
+ }
+ } /* end contributing */
+ } /* end in region */
+ } /* end sides */
+
+ if (done) {
+ /* lower all cells to bottom ele that have lower lying contributing neighbours */
+ if (n_stack[top].n_ngbrs > 0)
+ cseg_put(&ele, &bottom_ele, r, c);
+ n_stack[top].next_side = sides;
+ top--;
+ }
+ }
+ /* <-- carve upstream from spill point */
+
+ /* carve downstream from sink bottom */
+ bseg_get(&draindir, &drain_val, bottom_r, bottom_c);
+ if (drain_val < 1) {
+ G_warning(_("Can't carve downstream from r %d, c %d"), bottom_r,
+ bottom_c);
+ return 0;
+ }
+
+ r_nbr = bottom_r + asp_r[(int) drain_val];
+ c_nbr = bottom_c + asp_c[(int) drain_val];
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+
+ /* go downstream up to peak and down to bottom again */
+ while (ele_nbr >= bottom_ele) {
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+ cseg_put(&ele, &bottom_ele, r_nbr, c_nbr);
+ carved++;
+
+ if (drain_val > 0) {
+ r_nbr = r_nbr + asp_r[(int) drain_val];
+ c_nbr = c_nbr + asp_c[(int) drain_val];
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ }
+ else {
+ G_debug(2, "%d carved to edge", carved);
+ return carved;
+ }
+ }
+
+ return carved;
+}
+
+/* remove sinks */
+int hydro_con(void)
+{
+ int counter = 1;
+ int r, c, r_nbr, c_nbr, peak_r, peak_c, noedge, force_carve;
+ CELL ele_val, ele_nbr, ele_last, peak_ele, bottom_ele;
+ int down_uphill, down_downhill, n_cells;
+ struct sink_list *last_sink;
+ 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 drain_val, flag_value;
+ int n_splits, n_to_fill;
+ int n_filled, n_carved, n_just_a_bit, n_processed;
+ int carve_to_edge, skipme;
+
+ stack_alloc = nrows + ncols;
+
+ n_stack = (struct ns *)G_malloc(stack_alloc * sizeof(struct ns));
+
+ G_message(_("Processing %d sinks"), n_sinks);
+
+ n_filled = n_carved = n_just_a_bit = n_processed = 0;
+
+ while (first_sink) {
+ G_percent(counter++, n_sinks, 1);
+
+ r = first_sink->r;
+ c = first_sink->c;
+
+ cseg_get(&ele, &ele_val, r, c);
+ bottom_ele = ele_val;
+
+ skipme = 0;
+
+ G_debug(1, "get spill point");
+
+ /* go downstream, get spill point */
+ peak_r = r;
+ peak_c = c;
+ bseg_get(&draindir, &drain_val, r, c);
+
+ r_nbr = r + asp_r[(int) drain_val];
+ c_nbr = c + asp_c[(int) drain_val];
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+
+ down_uphill = down_downhill = 0;
+ noedge = 1;
+ ele_last = bottom_ele;
+ n_cells = 0;
+ while (ele_nbr >= ele_last) {
+ n_cells++;
+ if (ele_nbr > ele_last) {
+ down_uphill = n_cells;
+ peak_r = r_nbr;
+ peak_c = c_nbr;
+ }
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+
+ if (drain_val > 0) {
+ r_nbr = r_nbr + asp_r[(int) drain_val];
+ c_nbr = c_nbr + asp_c[(int) drain_val];
+
+ ele_last = ele_nbr;
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ }
+ else {
+ break;
+ }
+ }
+
+ cseg_get(&ele, &peak_ele, peak_r, peak_c);
+
+ bseg_get(&bitflags, &flag_value, peak_r, peak_c);
+ noedge = (FLAG_GET(flag_value, EDGEFLAG) == 0);
+ r_nbr = peak_r;
+ c_nbr = peak_c;
+ ele_nbr = peak_ele;
+
+ /* check */
+ if (down_uphill == 0) {
+ G_debug(2, "sink already processed");
+ n_processed++;
+
+ last_sink = first_sink;
+ first_sink = first_sink->next;
+ G_free(last_sink);
+ continue;
+ }
+
+ carve_to_edge = 0;
+ if (noedge) {
+ G_debug(1, "go downstream from spill point");
+ /* go downstream from spill point until we reach bottom ele again */
+ while (ele_nbr >= bottom_ele) {
+ if (ele_nbr > bottom_ele)
+ down_downhill++;
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+
+ /* detect flat area ? */
+
+ if (drain_val > 0) {
+ r_nbr = r_nbr + asp_r[(int) drain_val];
+ c_nbr = c_nbr + asp_c[(int) drain_val];
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ }
+ else {
+ carve_to_edge = 1;
+ break;
+ }
+ }
+ }
+ G_debug(1, "%d cells downstream from spill point", down_downhill);
+
+ /* count number of cells that would be filled from current spill point upstream */
+ G_debug(1, "count fill");
+ n_to_fill = count_fill(peak_r, peak_c, peak_ele, &n_splits, NULL);
+
+ skipme = (!do_all && n_to_fill <= size_max);
+
+ nat_thresh = 10;
+
+ /* shallow inland pans:
+ * - bottom is close to center of the depression
+ * - long channel to carve downstream from spill point
+ * - shallowness */
+ if (n_to_fill > nat_thresh && peak_ele - bottom_ele < 10 && sqrt(n_to_fill) / 4.0 < down_uphill)
+ G_debug(1, "shallow inland pan?");
+
+
+ /*************************/
+ /* set correction method */
+ /*************************/
+
+ /* force carving if edge not past spill point ? */
+ force_carve = (noedge == 0);
+ force_carve = 0;
+ if (!force_carve && (do_all || n_mod_max >= 4) &&
+ down_uphill + down_downhill <= 4 &&
+ down_uphill + down_downhill <= n_to_fill)
+ force_carve = 0;
+
+ /************************************/
+ /* apply selected correction method */
+ /************************************/
+
+ /* force carve channel to edge if edge not past spill point */
+ if (force_carve) {
+ G_debug(1, "force carve small channel or to edge");
+ carve(first_sink->r, first_sink->c, bottom_ele, peak_r, peak_c);
+ n_carved++;
+ }
+ /* impact reduction algorithm */
+ else if (!skipme) {
+ int next_r, next_c;
+ int n_new_splits, is_edge = 0, this_is_flat;
+ CELL min_bottom_ele, last_ele, next_bottom_ele;
+ int least_impact, li_max;
+ int least_impact_r, least_impact_c;
+ CELL least_impact_ele;
+
+ /* initialize least impact */
+ li_max = n_to_fill + down_uphill + down_downhill;
+ least_impact = down_uphill + down_downhill;
+ least_impact_r = first_sink->r;
+ least_impact_c = first_sink->c;
+ least_impact_ele = bottom_ele;
+
+ /* IRA --> */
+
+ G_debug(1, "impact reduction algorithm");
+
+ /* start with sink bottom, proceed to spill point */
+ /* fill incrementally, for each step check if IRA condition fulfilled */
+ /* proceed until condition no longer met, then use least impact solution */
+ /* the determined new sink bottom can be
+ * - the original sink bottom: no filling, only carving
+ * - the spill point: no carving, only filling
+ * - any cell in between the two */
+
+ next_r = first_sink->r;
+ next_c = first_sink->c;
+ min_bottom_ele = bottom_ele;
+
+ G_debug(1, "find new fill point");
+
+ /* include no filling, full carving, and full filling, no carving */
+ while (bottom_ele <= peak_ele) {
+
+ /* get n_to_fill, down_uphill, down_downhill */
+
+ down_uphill = 0;
+ down_downhill = 0;
+
+ /* count cells to fill */
+ next_bottom_ele = peak_ele;
+ n_to_fill = count_fill(next_r, next_c, bottom_ele,
+ &n_new_splits, &next_bottom_ele);
+
+ /* count channel length from spill point */
+ /* go downstream from spill point until we reach bottom ele again */
+ G_debug(2, "count channel length from spill point");
+ r_nbr = peak_r;
+ c_nbr = peak_c;
+ last_ele = ele_nbr = peak_ele;
+ min_bottom_ele = bottom_ele;
+ this_is_flat = 0;
+ while (ele_nbr >= bottom_ele && this_is_flat < 3) {
+ if (ele_nbr > bottom_ele) {
+ down_downhill++;
+ /* catch min 2 consecutive cells in flat area with same ele */
+ if (is_flat(r_nbr, c_nbr, ele_nbr)) {
+ if (ele_nbr == last_ele) {
+ this_is_flat++;
+ if (this_is_flat > 1 && ele_nbr > min_bottom_ele)
+ min_bottom_ele = bottom_ele + 1;
+ }
+ else
+ this_is_flat = 1;
+ }
+ if (next_bottom_ele > ele_nbr)
+ next_bottom_ele = ele_nbr;
+ }
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+
+ if (drain_val > 0) {
+ r_nbr = r_nbr + asp_r[(int) drain_val];
+ c_nbr = c_nbr + asp_c[(int) drain_val];
+
+ last_ele = ele_nbr;
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ }
+ else {
+ G_debug(2, "reached edge past spill point");
+ break;
+ }
+ }
+
+ /* carving ok */
+ if (this_is_flat < 3) {
+ min_bottom_ele = bottom_ele;
+ /* count cells to carve from sink bottom to spill point */
+ ele_last = ele_nbr = bottom_ele;
+ r_nbr = next_r;
+ c_nbr = next_c;
+ n_cells = 0;
+ while (ele_nbr >= ele_last) {
+ n_cells++;
+ if (ele_nbr > bottom_ele) {
+ down_uphill = n_cells;
+ }
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+
+ if (drain_val > 0) {
+ r_nbr = r_nbr + asp_r[(int) drain_val];
+ c_nbr = c_nbr + asp_c[(int) drain_val];
+
+ ele_last = ele_nbr;
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ }
+ else {
+ break;
+ }
+ }
+ /* remember least impact */
+ if (least_impact > n_to_fill + down_uphill + down_downhill) {
+ least_impact = n_to_fill + down_uphill + down_downhill;
+ least_impact_r = next_r;
+ least_impact_c = next_c;
+ least_impact_ele = bottom_ele;
+ }
+ }
+ /* carving not ok, need to fill more */
+ else {
+ G_debug(2, "carving not ok");
+ least_impact = li_max;
+ least_impact_r = peak_r;
+ least_impact_c = peak_c;
+ least_impact_ele = peak_ele;
+ min_bottom_ele = bottom_ele + 1;
+ }
+
+ G_debug(2,
+ "channel length to spill point %d, channel length from spill point %d, cells to fill %d",
+ down_uphill, down_downhill, n_to_fill);
+
+ /* check IRA condition */
+ /* preference for channel carving, relaxed stream split condition */
+ /* continue filling ? */
+ /* remains a mystery why this condition results in less modifications... */
+ /* sqrt(n_to_fill) < (down_uphill + down_downhill) / 2.0 ||
+ * n_new_splits <= n_splits_max */
+ if (sqrt(n_to_fill) < (down_uphill + down_downhill) / 1.0 ||
+ n_new_splits <= 4 ||
+ bottom_ele < min_bottom_ele ||
+ least_impact == li_max) {
+ /* continue with IRA */
+ G_debug(2, "IRA conditions fulfilled");
+ }
+ /* IRA condition no longer fulfilled */
+ else {
+ G_debug(2, "IRA conditions not met");
+ G_debug(2, "n to fill: %d", n_to_fill);
+ G_debug(2, "down_uphill: %d", down_uphill);
+ G_debug(2, "down_downhill: %d", down_downhill);
+ break;
+ }
+
+ if (bottom_ele >= peak_ele)
+ break;
+
+ /* below spill point, get next higher cell downstream of current bottom */
+ G_debug(2, "get next fill point candidate");
+
+ r_nbr = next_r;
+ c_nbr = next_c;
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+ if (drain_val > 0) {
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ /* go to next higher cell downstream */
+ while (ele_nbr == bottom_ele || ele_nbr < min_bottom_ele) {
+ bseg_get(&draindir, &drain_val, r_nbr, c_nbr);
+
+ if (drain_val > 0) {
+ r_nbr = r_nbr + asp_r[(int) drain_val];
+ c_nbr = c_nbr + asp_c[(int) drain_val];
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+ }
+ else {
+ is_edge = 1;
+ break;
+ }
+ }
+ /* got next downstream cell that's higher than sink bottom */
+ G_debug(2, "got next downstream cell");
+ }
+ else
+ is_edge = 1;
+
+ if (is_edge) {
+ G_warning(_("IRA reached edge !!!"));
+ break;
+ }
+
+ /* increase ele or go to next point */
+ /*
+ if (ele_nbr > bottom_ele + 1 && bottom_ele + 1 >= min_bottom_ele) {
+ G_debug(2, "increase ele, keep point");
+ bottom_ele = bottom_ele + 1;
+ }
+ */
+ if (ele_nbr > next_bottom_ele && next_bottom_ele >= min_bottom_ele) {
+ G_debug(2, "increase ele, keep point");
+ bottom_ele = next_bottom_ele;
+ }
+ else {
+ next_r = r_nbr;
+ next_c = c_nbr;
+ bottom_ele = ele_nbr;
+ }
+ }
+ /* <-- IRA */
+
+ if (least_impact == li_max)
+ G_warning(_("IRA error"));
+
+ if (do_all || least_impact <= n_mod_max) {
+
+ /* (partially) fill the sink */
+ n_to_fill = 0;
+ cseg_get(&ele, &ele_val, first_sink->r, first_sink->c);
+ if (least_impact_ele != ele_val) {
+ G_debug(1, "IRA sink filling");
+
+ cseg_get(&ele, &ele_val, least_impact_r, least_impact_c);
+ if (ele_val != least_impact_ele) {
+ G_debug(1, "changing ele from %d to %d", ele_val, least_impact_ele);
+ cseg_put(&ele, &least_impact_ele, least_impact_r, least_impact_c);
+ }
+ n_to_fill = fill_sink(least_impact_r, least_impact_c, least_impact_ele);
+ first_sink->r = least_impact_r;
+ first_sink->c = least_impact_c;
+ if (least_impact_ele < peak_ele)
+ n_just_a_bit++;
+ else
+ n_filled++;
+ }
+ else {
+ n_carved++;
+ G_debug(1, "IRA at bottom");
+ }
+
+ /* carve if not completely filled */
+ down_downhill = 0;
+ if (least_impact_ele < peak_ele) {
+ G_debug(2, "IRA carving");
+ down_downhill = carve(first_sink->r, first_sink->c,
+ least_impact_ele, peak_r, peak_c);
+ }
+ }
+ }
+
+ G_debug(1, "get next sink");
+
+ last_sink = first_sink;
+ first_sink = first_sink->next;
+ G_free(last_sink);
+ }
+
+ G_verbose_message
+ ("------------------------------------------------------");
+ G_verbose_message(_("%d sinks processed"), n_sinks);
+ G_verbose_message(_("%d sinks within larger sinks"), n_processed);
+ G_verbose_message(_("%d sinks completely filled"), n_filled);
+ G_verbose_message(_("%d sinks partially filled and carved with IRA"), n_just_a_bit);
+ G_verbose_message(_("%d channels carved"), n_carved);
+ G_verbose_message
+ ("------------------------------------------------------");
+
+ return 1;
+}
+
+int one_cell_extrema(int do_peaks, int do_pits, int all)
+{
+ int r, c, r_nbr, c_nbr;
+ int ct_dir;
+ int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+ int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+ int skipme;
+ CELL ele_min, ele_max, ele_this, ele_nbr;
+ unsigned int n_pits = 0;
+ unsigned int n_peaks = 0;
+ unsigned int counter = 0;
+ double sumofsquares, mean, stddev, upperci, lowerci;
+ int n_valid;
+
+ G_message(_("Remove one cell extremas"));
+
+ /* not for edges */
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 2);
+
+ for (c = 0; c < ncols; c++) {
+
+ skipme = 0;
+ cseg_get(&ele, &ele_this, r, c);
+
+ if (Rast_is_c_null_value(&ele_this))
+ continue;
+
+ ele_min = INT_MAX;
+ ele_max = INT_MIN;
+ sumofsquares = mean = n_valid = 0;
+
+ /* TODO: better use 5x5 neighbourhood ? */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
+ /* check that neighbour is within region */
+ if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
+ c_nbr < ncols) {
+
+ cseg_get(&ele, &ele_nbr, r_nbr, c_nbr);
+
+ if (Rast_is_c_null_value(&ele_nbr)) {
+ skipme = 1;
+ break;
+ }
+
+ n_valid++;
+ sumofsquares += ele_nbr * ele_nbr;
+ mean += ele_nbr;
+
+ if (ele_min > ele_nbr) {
+ ele_min = ele_nbr;
+ }
+ if (ele_max < ele_nbr) {
+ ele_max = ele_nbr;
+ }
+ }
+ else {
+ skipme = 1;
+ break;
+ }
+ }
+
+ if (skipme)
+ continue;
+
+ counter++;
+
+ if (all) {
+ upperci = ele_this - 1;
+ lowerci = ele_this + 1;
+ }
+ else {
+ /* remove extremas only if outside 95% CI = mean +- 1.96 * SD */
+ mean /= n_valid;
+ stddev = sqrt(sumofsquares / n_valid - mean * mean);
+ upperci = mean + 1.96 * stddev;
+ lowerci = mean - 1.96 * stddev;
+ }
+
+ /* one cell pit, lower than all surrounding neighbours */
+ if (do_pits) {
+ if (ele_this < ele_min && ele_this < lowerci) {
+ ele_this = ele_min;
+ cseg_put(&ele, &ele_this, r, c);
+ n_pits++;
+ }
+ }
+
+ /* one cell peak, higher than all surrounding neighbours */
+ if (do_peaks) {
+ if (ele_this > ele_max && ele_this > upperci) {
+ ele_this = ele_max;
+ cseg_put(&ele, &ele_this, r, c);
+ n_peaks++;
+ }
+ }
+ }
+ }
+ G_percent(nrows, nrows, 1); /* finish it */
+
+ G_verbose_message("%u cells checked", counter);
+ G_verbose_message("%u one-cell peaks removed", n_peaks);
+ G_verbose_message("%u one-cell pits removed", n_pits);
+
+ return 1;
+}
Property changes on: grass-addons/grass7/raster/r.hydrodem/hydro_con.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/raster/r.hydrodem/init_search.c
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/init_search.c 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/init_search.c 2012-01-09 17:41:56 UTC (rev 50109)
@@ -5,28 +5,39 @@
int init_search(int depr_fd)
{
int r, c, r_nbr, c_nbr, ct_dir;
- CELL *depr_buf, ele_value;
+ void *depr_buf, *depr_ptr;
+ CELL 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;
+ unsigned int n_null_cells = nrows * ncols - n_points;
+ int depr_map_type, depr_size;
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
+ G_message(_("Set edge points"));
+
+ if (depr_fd > -1) {
+ depr_map_type = Rast_get_map_type(depr_fd);
+ depr_size = Rast_cell_size(depr_map_type);
+ depr_buf = Rast_allocate_buf(depr_map_type);
+ }
+ else {
depr_buf = NULL;
+ depr_map_type = 0;
+ depr_size = 0;
+ }
+ depr_ptr = 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);
+ Rast_get_row(depr_fd, depr_buf, r, depr_map_type);
+ depr_ptr = depr_buf;
}
for (c = 0; c < ncols; c++) {
@@ -34,10 +45,15 @@
bseg_get(&bitflags, &flag_value, r, c);
is_null = FLAG_GET(flag_value, NULLFLAG);
- if (is_null)
+ if (is_null) {
+ if (depr_fd > -1)
+ depr_ptr = G_incr_void_ptr(depr_ptr, depr_size);
+
continue;
+ }
asp_value = 0;
+
if (r == 0 || r == nrows - 1 || c == 0 || c == ncols - 1) {
if (r == 0 && c == 0)
@@ -57,60 +73,62 @@
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);
+ cseg_get(&ele, &ele_value, r, c);
FLAG_SET(flag_value, EDGEFLAG);
- bseg_put(&bitflags, &flag_value, r, c);
- bseg_put(&asp, &asp_value, r, c);
+ heap_add(r, c, ele_value, asp_value, flag_value);
+
+ if (depr_fd > -1)
+ depr_ptr = G_incr_void_ptr(depr_ptr, depr_size);
+
continue;
}
- /* any neighbour NULL ? */
- for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r, c (r_nbr, c_nbr) for neighbours */
- r_nbr = r + nextdr[ct_dir];
- c_nbr = c + nextdc[ct_dir];
+ if (n_null_cells > 0) {
+ /* any neighbour NULL ? */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ r_nbr = r + nextdr[ct_dir];
+ c_nbr = c + nextdc[ct_dir];
- bseg_get(&bitflags, &flag_value_nbr, r_nbr, c_nbr);
- is_null = FLAG_GET(flag_value_nbr, NULLFLAG);
+ 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);
+ if (is_null) {
+ asp_value = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+ cseg_get(&ele, &ele_value, r, c);
+ FLAG_SET(flag_value, EDGEFLAG);
+ heap_add(r, c, ele_value, asp_value, flag_value);
- break;
+ break;
+ }
}
}
- if (asp_value) /* some neighbour was NULL, point added to list */
+
+ if (asp_value) { /* some neighbour was NULL, point added to list */
+
+ if (depr_fd > -1)
+ depr_ptr = G_incr_void_ptr(depr_ptr, depr_size);
+
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);
+ if (depr_fd > -1) {
+ DCELL depr_val = Rast_get_d_value(depr_ptr, depr_map_type);
+
+ if (!Rast_is_null_value(depr_ptr, depr_map_type) && depr_val != 0) {
+ cseg_get(&ele, &ele_value, r, c);
FLAG_SET(flag_value, DEPRFLAG);
- bseg_put(&bitflags, &flag_value, r, c);
- bseg_put(&asp, &asp_value, r, c);
+ heap_add(r, c, ele_value, asp_value, flag_value);
n_depr_cells++;
}
+ depr_ptr = G_incr_void_ptr(depr_ptr, depr_size);
}
}
}
G_percent(nrows, nrows, 2); /* finish it */
- if (depr_fd >= 0) {
- Rast_close(depr_fd);
+ if (depr_fd > -1) {
G_free(depr_buf);
}
Modified: grass-addons/grass7/raster/r.hydrodem/load.c
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/load.c 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/load.c 2012-01-09 17:41:56 UTC (rev 50109)
@@ -1,3 +1,4 @@
+#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/glocale.h>
#include "local_proto.h"
@@ -2,36 +3,25 @@
-/* need elevation map, do A* search on elevation like for r.watershed */
-
int ele_round(double x)
{
- if (x >= 0.0)
- return x + .5;
- else
- return x - .5;
+ return (x > 0.0 ? x + .5 : x - .5);
}
/*
- * loads elevation and optional flow accumulation map to memory and
- * gets start points for A* Search
+ * loads elevation map to segment file and gets start points for A* Search,
* start points are edges
*/
-int load_maps(int ele_fd, int acc_fd)
+int load_map(int ele_fd, int depr_fd)
{
int r, c;
- void *ele_buf, *ptr, *acc_buf = NULL, *acc_ptr = NULL;
- CELL ele_value, *stream_id;
- DCELL dvalue, acc_value;
- size_t ele_size, acc_size = 0;
- int ele_map_type, acc_map_type = 0;
- WAT_ALT *wabuf;
- char *flag_value_buf, *aspect;
+ char flag_value, asp_value = 0;
+ void *ele_buf, *ptr;
+ CELL ele_value;
+ DCELL dvalue;
+ int ele_size;
+ int ele_map_type;
- if (acc_fd < 0)
- G_message(_("Load elevation map..."));
- else
- G_message(_("Load input maps..."));
+ G_message(_("Load elevation map"));
n_search_points = n_points = 0;
- G_debug(1, "get buffers");
ele_map_type = Rast_get_map_type(ele_fd);
@@ -45,25 +35,10 @@
return -1;
}
- if (acc_fd >= 0) {
- 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;
- }
- }
-
ele_scale = 1;
if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
ele_scale = 1000; /* should be enough to do the trick */
- 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));
-
G_debug(1, "start loading %d rows, %d cols", nrows, ncols);
for (r = 0; r < nrows; r++) {
@@ -72,26 +47,17 @@
Rast_get_row(ele_fd, ele_buf, r, ele_map_type);
ptr = ele_buf;
- if (acc_fd >= 0) {
- Rast_get_row(acc_fd, acc_buf, r, acc_map_type);
- acc_ptr = acc_buf;
- }
-
for (c = 0; c < ncols; c++) {
- flag_value_buf[c] = 0;
- aspect[c] = 0;
- stream_id[c] = 0;
+ flag_value = 0;
/* check for masked and NULL cells */
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);
+ FLAG_SET(flag_value, NULLFLAG);
+ FLAG_SET(flag_value, INLISTFLAG);
+ FLAG_SET(flag_value, WORKEDFLAG);
+ FLAG_SET(flag_value, WORKED2FLAG);
Rast_set_c_null_value(&ele_value, 1);
- /* flow accumulation */
- Rast_set_d_null_value(&acc_value, 1);
}
else {
switch (ele_map_type) {
@@ -109,54 +75,23 @@
ele_value = ele_round(dvalue);
break;
}
- if (acc_fd < 0)
- acc_value = 1;
- else {
- 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:
- acc_value = *((CELL *) acc_ptr);
- break;
- case FCELL_TYPE:
- acc_value = *((FCELL *) acc_ptr);
- break;
- case DCELL_TYPE:
- acc_value = *((DCELL *) acc_ptr);
- break;
- }
- }
-
n_points++;
}
- wabuf[c].wat = acc_value;
- wabuf[c].ele = ele_value;
+ cseg_put(&ele, &ele_value, r, c);
+ bseg_put(&draindir, &asp_value, r, c);
+ bseg_put(&bitflags, &flag_value, r, c);
+
ptr = G_incr_void_ptr(ptr, ele_size);
- 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 */
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) {
- Rast_close(acc_fd);
- G_free(acc_buf);
- }
-
G_debug(1, "%d non-NULL cells", n_points);
- return n_points;
+ return 1;
}
Modified: grass-addons/grass7/raster/r.hydrodem/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/local_proto.h 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/local_proto.h 2012-01-09 17:41:56 UTC (rev 50109)
@@ -2,39 +2,38 @@
#ifndef __LOCAL_PROTO_H__
#define __LOCAL_PROTO_H__
+#include <grass/gis.h>
#include <grass/raster.h>
-#include "flag.h"
#include "seg.h"
+#include "flag.h"
#define INDEX(r, c) ((r) * ncols + (c))
#define MAXDEPTH 1000 /* maximum supported tree depth of stream network */
-#define POINT struct a_point
-POINT {
+struct ddir
+{
+ int pos;
+ int dir;
+};
+
+struct point
+{
int r, c;
};
-#define HEAP_PNT struct heap_point
-HEAP_PNT {
+struct heap_point {
unsigned int added;
CELL ele;
- POINT pnt;
+ int r, c;
};
-#define WAT_ALT struct wat_altitude
-WAT_ALT {
- CELL ele;
- DCELL wat;
+struct sink_list
+{
+ int r, c;
+ struct sink_list *next;
};
-/* global variables */
-#ifdef MAIN
-# define GLOBAL
-#else
-# define GLOBAL extern
-#endif
-
-GLOBAL struct snode
+struct snode
{
int r, c;
int id;
@@ -42,51 +41,48 @@
int n_trib_total; /* number of all upstream stream segments */
int n_alloc; /* n allocated tributaries */
int *trib;
+ double *acc;
} *stream_node;
-GLOBAL int nrows, ncols;
-GLOBAL unsigned int n_search_points, n_points, nxt_avail_pt;
-GLOBAL unsigned int heap_size;
-GLOBAL unsigned int n_stream_nodes, n_alloc_nodes;
-GLOBAL POINT *outlets;
-GLOBAL unsigned int n_outlets, n_alloc_outlets;
-GLOBAL char drain[3][3];
-GLOBAL char sides;
-GLOBAL int c_fac;
-GLOBAL int ele_scale;
-GLOBAL int have_depressions;
+extern int nrows, ncols;
+extern unsigned int n_search_points, n_points, nxt_avail_pt;
+extern unsigned int heap_size;
+extern unsigned int n_sinks;
+extern int n_mod_max, size_max;
+extern int do_all, keep_nat, nat_thresh;
+extern unsigned int n_stream_nodes, n_alloc_nodes;
+extern struct point *outlets;
+extern struct sink_list *sinks, *first_sink;
+extern unsigned int n_outlets, n_alloc_outlets;
+extern char drain[3][3];
+extern unsigned int first_cum;
+extern char sides;
+extern int c_fac;
+extern int ele_scale;
+extern struct RB_TREE *draintree;
-GLOBAL SSEG search_heap;
-GLOBAL SSEG astar_pts;
-GLOBAL BSEG bitflags;
-GLOBAL SSEG watalt;
-GLOBAL BSEG asp;
-GLOBAL CSEG stream;
+extern SSEG search_heap;
+extern SSEG astar_pts;
+extern BSEG bitflags;
+extern CSEG ele;
+extern BSEG draindir;
+extern CSEG stream;
/* load.c */
-int load_maps(int, int);
+int load_map(int, int);
/* init_search.c */
int init_search(int);
/* do_astar.c */
int do_astar(void);
-unsigned int heap_add(int, int, CELL);
+unsigned int heap_add(int, int, CELL, char, char);
-/* streams.c */
-int do_accum(double);
-int extract_streams(double, double, int, int);
+/* hydro_con.c */
+int hydro_con(void);
+int one_cell_extrema(int, int, int);
-/* thin.c */
-int thin_streams(void);
-
-/* basins.c */
-int basin_borders(void);
-
-/* del_streams.c */
-int del_streams(int);
-
/* close.c */
-int close_maps(char *, char *, char *);
+int close_map(char *, int);
#endif /* __LOCAL_PROTO_H__ */
Modified: grass-addons/grass7/raster/r.hydrodem/main.c
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/main.c 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/main.c 2012-01-09 17:41:56 UTC (rev 50109)
@@ -1,11 +1,10 @@
/****************************************************************************
*
- * MODULE: r.stream.extract
+ * MODULE: r.hydrodem
* AUTHOR(S): Markus Metz <markus.metz.giswork gmail.com>
* PURPOSE: Hydrological analysis
- * Extracts stream networks from accumulation raster with
- * given threshold
+ * DEM hydrological conditioning based on A* Search
* COPYRIGHT: (C) 1999-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public
@@ -16,59 +15,59 @@
#include <stdlib.h>
#include <string.h>
#include <float.h>
-#include <math.h>
+#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/glocale.h>
-#define MAIN
#include "local_proto.h"
-/* global variables */
int nrows, ncols;
unsigned int n_search_points, n_points, nxt_avail_pt;
unsigned int heap_size;
+unsigned int n_sinks;
+int n_mod_max, size_max;
+int do_all, keep_nat, nat_thresh;
unsigned int n_stream_nodes, n_alloc_nodes;
-POINT *outlets;
+struct point *outlets;
+struct sink_list *sinks, *first_sink;
unsigned int n_outlets, n_alloc_outlets;
+
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 ele;
+BSEG draindir;
CSEG stream;
-CELL *astar_order;
int main(int argc, char *argv[])
{
struct
{
- struct Option *ele, *acc, *depression;
- struct Option *threshold, *d8cut;
- struct Option *mont_exp;
- struct Option *min_stream_length;
- struct Option *memory;
+ struct Option *ele, *depr, *memory;
} input;
struct
{
- struct Option *stream_rast;
- struct Option *stream_vect;
- struct Option *dir_rast;
+ struct Option *ele_hydro;
+ struct Option *mod_max;
+ struct Option *size_max;
+ struct Flag *do_all;
} output;
struct GModule *module;
- int ele_fd, acc_fd, depr_fd;
- double threshold, d8cut, mont_exp;
- int min_stream_length = 0, memory;
+ int ele_fd, ele_map_type, depr_fd;
+ int memory;
int seg_cols, seg_rows;
double seg2kb;
int num_open_segs, num_open_array_segs, num_seg_total;
double memory_divisor, heap_mem, disk_space;
const char *mapset;
+ struct Colors colors;
G_gisinit(argv[0]);
@@ -76,64 +75,21 @@
module = G_define_module();
G_add_keyword(_("raster"));
G_add_keyword(_("hydrology"));
- module->description = _("Stream network extraction");
+ module->description = _("Hydrological conditioning, sink removal");
input.ele = G_define_standard_option(G_OPT_R_INPUT);
- input.ele->key = "elevation";
+ input.ele->key = "input";
input.ele->label = _("Elevation map");
- input.ele->description = _("Elevation on which entire analysis is based");
+ input.ele->description =
+ _("Elevation map to be hydrologically corrected");
- input.acc = G_define_standard_option(G_OPT_R_INPUT);
- input.acc->key = "accumulation";
- input.acc->label = _("Accumulation map");
- input.acc->required = NO;
- input.acc->description =
- _("Stream extraction will use provided accumulation instead of calculating it anew");
+ input.depr = G_define_standard_option(G_OPT_R_INPUT);
+ input.depr->key = "depression";
+ input.depr->required = NO;
+ input.depr->label = _("Depression map");
+ input.depr->description =
+ _("Map indicating real depressions that must not be modified");
- input.depression = G_define_standard_option(G_OPT_R_INPUT);
- input.depression->key = "depression";
- input.depression->label = _("Map with real depressions");
- input.depression->required = NO;
- input.depression->description =
- _("Streams will not be routed out of real depressions");
-
- input.threshold = G_define_option();
- input.threshold->key = "threshold";
- input.threshold->label = _("Minimum flow accumulation for streams");
- input.threshold->description = _("Must be > 0");
- input.threshold->required = YES;
- input.threshold->type = TYPE_DOUBLE;
-
- input.d8cut = G_define_option();
- input.d8cut->key = "d8cut";
- input.d8cut->label = _("Use SFD above this threshold");
- input.d8cut->description =
- _("If accumulation is larger than d8cut, SFD is used instead of MFD."
- " Applies only if no accumulation map is given.");
- input.d8cut->required = NO;
- input.d8cut->answer = "infinity";
- input.d8cut->type = TYPE_DOUBLE;
-
- input.mont_exp = G_define_option();
- input.mont_exp->key = "mexp";
- input.mont_exp->type = TYPE_DOUBLE;
- input.mont_exp->required = NO;
- input.mont_exp->answer = "0";
- input.mont_exp->label =
- _("Montgomery exponent for slope, disabled with 0");
- input.mont_exp->description =
- _("Montgomery: accumulation is multiplied with pow(slope,mexp) and then compared with threshold.");
-
- input.min_stream_length = G_define_option();
- input.min_stream_length->key = "stream_length";
- input.min_stream_length->type = TYPE_INTEGER;
- input.min_stream_length->required = NO;
- input.min_stream_length->answer = "0";
- input.min_stream_length->label =
- _("Delete stream segments shorter than stream_length cells.");
- 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;
@@ -141,27 +97,36 @@
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 =
- _("Output raster map with unique stream ids");
- output.stream_rast->required = NO;
- output.stream_rast->guisection = _("Output options");
+ output.ele_hydro = G_define_standard_option(G_OPT_R_OUTPUT);
+ output.ele_hydro->key = "output";
+ output.ele_hydro->description =
+ _("Name of hydrologically conditioned raster map");
+ output.ele_hydro->required = YES;
- output.stream_vect = G_define_standard_option(G_OPT_V_OUTPUT);
- output.stream_vect->key = "stream_vect";
- output.stream_vect->description =
- _("Output vector with unique stream ids");
- output.stream_vect->required = NO;
- output.stream_vect->guisection = _("Output options");
+ output.mod_max = G_define_option();
+ output.mod_max->key = "mod";
+ output.mod_max->description =
+ (_("Only remove sinks requiring not more than <mod> cell modifications."));
+ output.mod_max->type = TYPE_INTEGER;
+ output.mod_max->answer = "4";
+ output.mod_max->required = YES;
- output.dir_rast = G_define_standard_option(G_OPT_R_OUTPUT);
- output.dir_rast->key = "direction";
- output.dir_rast->description =
- _("Output raster map with flow direction");
- output.dir_rast->required = NO;
- output.dir_rast->guisection = _("Output options");
+ output.size_max = G_define_option();
+ output.size_max->key = "size";
+ output.size_max->description =
+ (_("Only remove sinks not larger than <size> cells."));
+ output.size_max->type = TYPE_INTEGER;
+ output.size_max->answer = "4";
+ output.size_max->required = YES;
+ output.do_all = G_define_flag();
+ output.do_all->key = 'a';
+ output.do_all->label = (_("Remove all sinks."));
+ output.do_all->description =
+ (_("By default only minor corrections are done to the DEM and "
+ "the result will not be 100% hydrologically correct."
+ "Use this flag to override default."));
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -172,59 +137,21 @@
/* input maps exist ? */
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_raster(input.acc->answer, ""))
- G_fatal_error(_("Raster map <%s> not found"), input.acc->answer);
+
+ if (input.depr->answer) {
+ if (!G_find_raster(input.depr->answer, ""))
+ G_fatal_error(_("Raster map <%s> not found"),
+ input.depr->answer);
}
- if (input.depression->answer) {
- if (!G_find_raster(input.depression->answer, ""))
- G_fatal_error(_("Raster map <%s> not found"), input.depression->answer);
- have_depressions = 1;
- }
- else
- have_depressions = 0;
+ if ((n_mod_max = atoi(output.mod_max->answer)) <= 0)
+ G_fatal_error(_("'%s' must be a positive integer"),
+ output.mod_max->key);
- /* threshold makes sense */
- threshold = atof(input.threshold->answer);
- if (threshold <= 0)
- G_fatal_error(_("Threshold must be > 0 but is %f"), threshold);
+ if ((size_max = atoi(output.size_max->answer)) <= 0)
+ G_fatal_error(_("'%s' must be a positive integer"),
+ output.size_max->key);
- /* d8cut */
- if (strcmp(input.d8cut->answer, "infinity") == 0) {
- d8cut = DBL_MAX;
- }
- else {
- d8cut = atof(input.d8cut->answer);
- if (d8cut < 0)
- G_fatal_error(_("d8cut must be positive or zero but is %f"),
- d8cut);
- }
-
- /* Montgomery stream initiation */
- if (input.mont_exp->answer) {
- mont_exp = atof(input.mont_exp->answer);
- if (mont_exp < 0)
- G_fatal_error(_("Montgomery exponent must be positive or zero but is %f"),
- mont_exp);
- if (mont_exp > 3)
- G_warning(_("Montgomery exponent is %f, recommended range is 0.0 - 3.0"),
- mont_exp);
- }
- else
- mont_exp = 0;
-
- /* Minimum stream segment length */
- if (input.min_stream_length->answer) {
- min_stream_length = atoi(input.min_stream_length->answer);
- if (min_stream_length < 0)
- G_fatal_error(_("Minimum stream length must be positive or zero but is %d"),
- min_stream_length);
- }
- else
- min_stream_length = 0;
-
if (input.memory->answer) {
memory = atoi(input.memory->answer);
if (memory <= 0)
@@ -234,11 +161,13 @@
else
memory = 300;
- /* Check for some output map */
- if ((output.stream_rast->answer == NULL)
- && (output.stream_vect->answer == NULL)
- && (output.dir_rast->answer == NULL)) {
- G_fatal_error(_("Sorry, you must choose at least one output map."));
+ do_all = output.do_all->answer;
+
+ if (!do_all) {
+ G_verbose_message(_("All sinks with max %d cells to be modified will be removed"),
+ n_mod_max);
+ G_verbose_message(_("All sinks not larger than %d cells will be removed."),
+ size_max);
}
/*********************/
@@ -251,22 +180,12 @@
if (ele_fd < 0)
G_fatal_error(_("Could not open input map %s"), input.ele->answer);
- if (input.acc->answer) {
- 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);
- }
- else
- acc_fd = -1;
-
- if (input.depression->answer) {
- mapset = G_find_raster2(input.depression->answer, "");
- depr_fd = Rast_open_old(input.depression->answer, mapset);
+ if (input.depr->answer) {
+ mapset = G_find_raster2(input.depr->answer, "");
+ depr_fd = Rast_open_old(input.depr->answer, mapset);
if (depr_fd < 0)
G_fatal_error(_("Could not open input map %s"),
- input.depression->answer);
+ input.depr->answer);
}
else
depr_fd = -1;
@@ -274,111 +193,77 @@
/* set global variables */
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 */
+ sides = 8; /* not a user option */
+ ele_map_type = Rast_get_map_type(ele_fd);
+
/* segment structures */
seg_rows = seg_cols = 64;
seg2kb = seg_rows * seg_cols / 1024.;
- /* 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.07
- *
- * astar_points: 8 byte -> 32 KB / segment
- * heap_points: 16 byte -> 64 KB / segment
- *
- * Total MB / segment: 0.16
- */
/* balance segment files */
- /* elevation + accumulation: * 2 */
- memory_divisor = sizeof(WAT_ALT) * 2;
- disk_space = sizeof(WAT_ALT);
- /* aspect: as is */
- memory_divisor += sizeof(char);
- disk_space += sizeof(char);
- /* stream ids: / 2 */
- memory_divisor += sizeof(int) / 2.;
- disk_space += sizeof(int);
+ /* elevation: * 2 */
+ memory_divisor = seg2kb * sizeof(CELL) * 2;
+ /* drainage direction: as is */
+ memory_divisor += seg2kb * sizeof(char);
/* flags: * 4 */
- memory_divisor += sizeof(char) * 4;
- disk_space += sizeof(char);
- /* astar_points: / 16 */
- /* ideally only a few but large segments */
- memory_divisor += sizeof(POINT) / 16.;
- disk_space += sizeof(POINT);
+ memory_divisor += seg2kb * sizeof(char) * 4;
/* heap points: / 4 */
- memory_divisor += sizeof(HEAP_PNT) / 4.;
- disk_space += sizeof(HEAP_PNT);
+ memory_divisor += seg2kb * sizeof(struct heap_point) / 4.;
/* KB -> MB */
- memory_divisor *= seg2kb / 1024.;
- disk_space *= seg2kb / 1024.;
+ memory_divisor /= 1024.;
num_open_segs = memory / memory_divisor;
- heap_mem = num_open_segs * seg2kb * sizeof(HEAP_PNT) /
+ /* heap_mem is in MB */
+ heap_mem = num_open_segs * seg2kb * sizeof(struct heap_point) /
(4. * 1024.);
num_seg_total = (ncols / seg_cols + 1) * (nrows / seg_rows + 1);
if (num_open_segs > num_seg_total) {
heap_mem += (num_open_segs - num_seg_total) * memory_divisor;
heap_mem -= (num_open_segs - num_seg_total) * seg2kb *
- sizeof(HEAP_PNT) / (4. * 1024.);
+ sizeof(struct heap_point) / (4. * 1024.);
num_open_segs = num_seg_total;
}
if (num_open_segs < 16) {
num_open_segs = 16;
- heap_mem = num_open_segs * seg2kb * sizeof(HEAP_PNT) /
+ heap_mem = num_open_segs * seg2kb * sizeof(struct heap_point) /
(4. * 1024.);
}
+ disk_space = (1. * sizeof(CELL) + 2 * sizeof(char) +
+ sizeof(struct heap_point));
+ disk_space *= (num_seg_total * seg2kb / 1024.); /* KB -> MB */
+
G_verbose_message(_("%.2f of data are kept in memory"),
100. * num_open_segs / num_seg_total);
- disk_space *= num_seg_total;
- if (disk_space < 1024.0)
- G_verbose_message(_("Will need up to %.2f MB of disk space"), disk_space);
- else
- G_verbose_message(_("Will need up to %.2f GB (%.0f MB) of disk space"),
- disk_space / 1024.0, disk_space);
+ G_verbose_message(_("Will need up to %.2f MB of disk space"), disk_space);
/* open segment files */
G_verbose_message(_("Create temporary files..."));
- seg_open(&watalt, nrows, ncols, seg_rows, seg_cols, num_open_segs * 2,
- sizeof(WAT_ALT), 1);
+ cseg_open(&ele, seg_rows, seg_cols, num_open_segs * 2);
if (num_open_segs * 2 > num_seg_total)
- heap_mem += (num_open_segs * 2 - num_seg_total) * seg2kb *
- sizeof(WAT_ALT) / 1024.;
- cseg_open(&stream, seg_rows, seg_cols, num_open_segs / 2.);
- bseg_open(&asp, seg_rows, seg_cols, num_open_segs);
+ heap_mem += (num_open_segs * 2 - num_seg_total) * seg2kb * sizeof(CELL) / 1024.;
+ bseg_open(&draindir, seg_rows, seg_cols, num_open_segs);
bseg_open(&bitflags, seg_rows, seg_cols, num_open_segs * 4);
if (num_open_segs * 4 > num_seg_total)
heap_mem += (num_open_segs * 4 - num_seg_total) * seg2kb / 1024.;
- /* load maps */
- 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");
- /* columns per segment */
- seg_cols = seg_rows * seg_rows;
- 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 */
- num_open_array_segs = num_open_segs / 16.;
- if (num_open_array_segs > num_seg_total)
- num_open_array_segs = num_seg_total;
- if (num_open_array_segs < 1)
- num_open_array_segs = 1;
+ /* load map */
+ if (load_map(ele_fd, depr_fd) < 0) {
+ cseg_close(&ele);
+ bseg_close(&draindir);
+ bseg_close(&bitflags);
+ G_fatal_error(_("Could not load input map"));
+ }
- G_debug(1, "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);
+ if (n_points == 0) {
+ cseg_close(&ele);
+ bseg_close(&draindir);
+ bseg_close(&bitflags);
+ G_fatal_error(_("No non-NULL cells loaded from input map"));
+ }
- /* one-based d-ary search_heap with astar_pts */
+ /* one-based d-ary search_heap */
G_debug(1, "open segments for A* search heap");
/* allowed memory for search heap in MB */
@@ -390,7 +275,8 @@
if (n_points % seg_cols > 0)
num_seg_total++;
/* no need to have more segments open than exist */
- num_open_array_segs = (1 << 20) * heap_mem / (seg_cols * sizeof(HEAP_PNT));
+ num_open_array_segs = (1 << 20) * heap_mem /
+ (seg_cols * sizeof(struct heap_point));
if (num_open_array_segs > num_seg_total)
num_open_array_segs = num_seg_total;
if (num_open_array_segs < 2)
@@ -402,53 +288,60 @@
/* 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);
+ num_open_array_segs, sizeof(struct heap_point), 1);
/********************/
/* processing */
/********************/
+ /* remove one cell extrema */
+ one_cell_extrema(1, 1, 0);
+
/* initialize A* search */
- if (init_search(depr_fd) < 0)
+ if (init_search(depr_fd) < 0) {
+ seg_close(&search_heap);
+ cseg_close(&ele);
+ bseg_close(&draindir);
+ bseg_close(&bitflags);
G_fatal_error(_("Could not initialize search"));
+ }
+ if (depr_fd >= 0) {
+ Rast_close(depr_fd);
+ }
+
/* sort elevation and get initial stream direction */
- if (do_astar() < 0)
+ if (do_astar() < 0) {
+ seg_close(&search_heap);
+ cseg_close(&ele);
+ bseg_close(&draindir);
+ bseg_close(&bitflags);
G_fatal_error(_("Could not sort elevation map"));
+ }
seg_close(&search_heap);
- if (acc_fd < 0) {
- /* accumulate surface flow */
- if (do_accum(d8cut) < 0)
- G_fatal_error(_("Could not calculate flow accumulation"));
+ /* hydrological corrections */
+ if (hydro_con() < 0) {
+ cseg_close(&ele);
+ bseg_close(&draindir);
+ bseg_close(&bitflags);
+ G_fatal_error(_("Could not apply hydrological conditioning"));
}
- /* extract streams */
- if (extract_streams
- (threshold, mont_exp, min_stream_length, acc_fd < 0) < 0)
- G_fatal_error(_("Could not extract streams"));
-
- seg_close(&astar_pts);
- seg_close(&watalt);
-
- /* thin streams */
- if (thin_streams() < 0)
- 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 delete short stream segments"));
+ /* write output maps */
+ if (close_map(output.ele_hydro->answer, ele_map_type) < 0) {
+ cseg_close(&ele);
+ bseg_close(&draindir);
+ bseg_close(&bitflags);
+ G_fatal_error(_("Could not write output map"));
}
- /* write output maps */
- if (close_maps(output.stream_rast->answer, output.stream_vect->answer,
- output.dir_rast->answer) < 0)
- G_fatal_error(_("Could not write output maps"));
-
- bseg_close(&asp);
- cseg_close(&stream);
+ cseg_close(&ele);
+ bseg_close(&draindir);
bseg_close(&bitflags);
+ Rast_read_colors(input.ele->answer, mapset, &colors);
+ Rast_write_colors(output.ele_hydro->answer, G_mapset(), &colors);
+
exit(EXIT_SUCCESS);
}
Copied: grass-addons/grass7/raster/r.hydrodem/r.hydrodem.html (from rev 50108, grass-addons/grass7/raster/r.hydrodem/r.stream.extract.html)
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/r.hydrodem.html (rev 0)
+++ grass-addons/grass7/raster/r.hydrodem/r.hydrodem.html 2012-01-09 17:41:56 UTC (rev 50109)
@@ -0,0 +1,84 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.hydrodem</em> applies hydrological conditioning (sink removal) to
+a required input <em>elevation</em> map. If the conditioned elevation
+map is going to be used as input elevation for <em>r.watershed</em>,
+only small sinks should be removed and the amount of modifications
+restricted with the <b>mod</b> option. For other modules such as
+<em>r.terraflow</em> or third-party software, full sink removal is
+recommended.
+
+<h2>OPTIONS</h2>
+
+<dl>
+<dt><b>input</b>
+
+<dd>Input map, required: Digital elevation model to be corrected. Gaps
+in the elevation map that are located within the area of interest should
+be filled beforehand, e.g. with <em>r.fillnulls</em> or
+<em>r.resamp.bspline</em>, to avoid distortions.
+<p>
+<dt><b>output</b>
+<dd>Output map, required: Hydrologically conditioned digital elevation
+model. By default, only minor modifications are done and not all sinks
+are removed.
+<p>
+<dt><b>size</b>
+<dd>All sinks of up to <b>size</b> cells will be removed.
+Default is 4, if in doubt, decrease and not increase.
+<p>
+<dt><b>mod</b>
+<dd>All sinks will be removed that require not more than <b>mod</b>
+cells to be modifed. Often, rather large sinks can be removed by
+carving through only few cells. Default is 4, if in doubt, increase
+and not decrease.
+<p>
+<dt><b>-a</b>
+<dd><b>Not recommended if input for <em>r.watershed</em> is generated.</b>
+<br>
+With the <b>-a</b> flag set, all sinks will be removed using an impact
+reduction approach based on Lindsay & Creed (2005). The output will be a
+depression-less digital elevation model, suitable for e.g.
+<em>r.terraflow</em> or other hydrological analyses that require a
+depression-less DEM as input.
+</dl>
+
+<h2>NOTES</h2>
+
+This module is designed for <em>r.watershed</em> with the purpose to
+slightly denoise a DEM prior to analysis. First, all one-cell peaks and
+pits are removed, then the actual hydrological corrections are applied.
+In most cases, the removal of one-cell extrema could already be
+sufficient to improve <em>r.watershed</em> results in difficult terrain,
+particularly nearly flat areas.
+<p>
+The impact reduction algorithm used by <em>r.hydrodem</em> is based on
+Lindsay & Creed (2005), with some additional checks for
+hydrological consistency. With complete sink removal, results of
+<em>r.terraflow</em> are very similar to results of <em>r.watershed</em>.
+<p>
+<em>r.hydrodem</em> uses the same method to determine drainage directions
+like <em>r.watershed</em>.
+<p>
+
+<h2>REFERENCES</h2>
+Lindsay, J. B., and Creed, I. F. 2005. Removal of artifact depressions
+from digital elevation models: towards a minimum impact approach.
+Hydrological Processes 19, 3113-3126.
+DOI: <a href=http://dx.doi.org/10.1002/hyp.5835>10.1002/hyp.5835</a>
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.terraflow.html">r.terraflow</a>,
+<a href="r.fill.dir.html">r.fill.dir</a>,
+<a href="r.fillnulls.html">r.fillnulls</a>,
+<a href="r.resamp.bspline.html">r.resamp.bspline</a>
+</em>
+
+<h2>AUTHOR</h2>
+Markus Metz
+
+<p><i>Last changed: $Date$</i>
Deleted: grass-addons/grass7/raster/r.hydrodem/r.stream.extract.html
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/r.stream.extract.html 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/r.stream.extract.html 2012-01-09 17:41:56 UTC (rev 50109)
@@ -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>
Deleted: grass-addons/grass7/raster/r.hydrodem/streams.c
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/streams.c 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/streams.c 2012-01-09 17:41:56 UTC (rev 50109)
@@ -1,730 +0,0 @@
-#include <stdlib.h>
-#include <math.h>
-#include <grass/raster.h>
-#include <grass/glocale.h>
-#include "local_proto.h"
-
-double mfd_pow(double base)
-{
- int i;
- double result = base;
-
- for (i = 2; i <= c_fac; i++) {
- result *= base;
- }
- return result;
-}
-
-int continue_stream(CELL stream_id, int r, int c, int r_max, int c_max,
- int *stream_no)
-{
- char aspect;
- 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 stream_node_step = 1000;
- char this_flag_value;
-
- G_debug(3, "continue stream");
-
- cseg_get(&stream, &curr_stream, r_max, c_max);
-
- if (curr_stream <= 0) {
- /* no confluence, just continue */
- 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(3, "confluence");
-
- /* new confluence */
- if (stream_node[curr_stream].r != r_max ||
- stream_node[curr_stream].c != c_max) {
- size_t new_size;
-
- G_debug(3, "new confluence");
- /* set new stream id */
- (*stream_no)++;
- /* add stream node */
- if (*stream_no >= n_alloc_nodes - 1) {
- n_alloc_nodes += stream_node_step;
- stream_node =
- (struct snode *)G_realloc(stream_node,
- n_alloc_nodes *
- sizeof(struct snode));
- }
- stream_node[*stream_no].r = r_max;
- stream_node[*stream_no].c = c_max;
- stream_node[*stream_no].id = *stream_no;
- stream_node[*stream_no].n_trib = 0;
- stream_node[*stream_no].n_trib_total = 0;
- stream_node[*stream_no].n_alloc = 0;
- stream_node[*stream_no].trib = NULL;
- n_stream_nodes++;
-
- /* debug */
- if (n_stream_nodes != *stream_no)
- G_warning(_("BUG: stream_no %d and n_stream_nodes %d out of sync"),
- *stream_no, n_stream_nodes);
-
-
- 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);
-
- /* 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;
-
- /* update stream IDs downstream */
- G_debug(3, "update stream IDs downstream");
- r_nbr = r_max;
- c_nbr = c_max;
- 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 (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(3, "existing confluence");
- /* add new tributary to stream node */
- if (stream_node[curr_stream].n_trib >=
- stream_node[curr_stream].n_alloc) {
- size_t new_size;
-
- stream_node[curr_stream].n_alloc += 2;
- new_size = stream_node[curr_stream].n_alloc * sizeof(int);
- stream_node[curr_stream].trib =
- (int *)G_realloc(stream_node[curr_stream].trib, new_size);
- }
-
- stream_node[curr_stream].trib[stream_node[curr_stream].n_trib++] =
- stream_id;
- }
-
- 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]);
-
- return 1;
-}
-
-/*
- * accumulate surface flow
- */
-int do_accum(double d8cut)
-{
- int r, c, dr, dc;
- 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;
- double dx, dy;
- int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
- int is_worked;
- char aspect;
- double max_acc, prop;
- int edge;
- int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
- int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
- int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
- int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
- unsigned int workedon, killer, count;
- char *flag_nbr, this_flag_value;
- POINT astarpoint;
- WAT_ALT wa;
-
- G_message(_("Calculate flow accumulation..."));
-
- /* 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);
-
- for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r, c (r_nbr, c_nbr) for neighbours */
- r_nbr = nextdr[ct_dir];
- c_nbr = nextdc[ct_dir];
- /* account for rare cases when ns_res != ew_res */
- dy = abs(r_nbr) * window.ns_res;
- dx = abs(c_nbr) * window.ew_res;
- if (ct_dir < 4)
- dist_to_nbr[ct_dir] = dx + dy;
- else
- dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
- }
-
- /* distribute and accumulate */
- 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;
-
- 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)];
- }
-
- r_max = dr;
- c_max = dc;
-
- 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 */
- /***************************************/
-
- max_weight = 0;
- sum_weight = 0;
- np_side = -1;
- mfd_cells = 0;
- astar_not_set = 1;
- 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++) {
- /* get r_nbr, c_nbr for neighbours */
- 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) {
-
- 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;
-
- /* 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) {
- if (ele_nbr[ct_dir] <= ele_val) {
- if (ele_nbr[ct_dir] < ele_val) {
- weight[ct_dir] =
- mfd_pow((ele_val -
- ele_nbr[ct_dir]) / dist_to_nbr[ct_dir]);
- }
- if (ele_nbr[ct_dir] == ele_val) {
- weight[ct_dir] =
- mfd_pow(0.5 / dist_to_nbr[ct_dir]);
- }
- sum_weight += weight[ct_dir];
- mfd_cells++;
-
- if (weight[ct_dir] > max_weight) {
- max_weight = weight[ct_dir];
- }
-
- if (dr == r_nbr && dc == c_nbr) {
- astar_not_set = 0;
- }
- }
- }
- if (dr == r_nbr && dc == c_nbr)
- np_side = ct_dir;
- }
- else
- edge = 1;
- if (edge)
- break;
- }
-
- /* do not distribute flow along edges, this causes artifacts */
- if (edge) {
- G_debug(3, "edge");
- continue;
- }
-
- /* honour A * path
- * mfd_cells == 0: fine, SFD along A * path
- * mfd_cells == 1 && astar_not_set == 0: fine, SFD along A * path
- * mfd_cells > 0 && astar_not_set == 1: A * path not included, add to mfd_cells
- */
-
- /************************************/
- /* distribute and accumulate flow */
- /************************************/
-
- /* MFD, A * path not included, add to mfd_cells */
- if (mfd_cells > 0 && astar_not_set == 1) {
- mfd_cells++;
- sum_weight += max_weight;
- weight[np_side] = max_weight;
- }
-
- /* use SFD (D8) if d8cut threshold exceeded */
- if (fabs(value) > d8cut)
- mfd_cells = 0;
-
- max_acc = -1;
-
- if (mfd_cells > 1) {
- prop = 0.0;
- for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r, c (r_nbr, c_nbr) for neighbours */
- r_nbr = r + nextdr[ct_dir];
- c_nbr = c + nextdc[ct_dir];
-
- /* check that neighbour is within region */
- if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
- c_nbr < ncols && weight[ct_dir] > -0.5) {
- is_worked = FLAG_GET(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];
-
- 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 */
- workedon++;
- }
- }
- }
- if (fabs(prop - 1.0) > 5E-6f) {
- G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
- prop);
- }
- }
- /* get out of depression in SFD mode */
- else {
- wa.wat = wat_nbr[np_side] + value;
- wa.ele = ele_nbr[np_side];
- seg_put(&watalt, (char *)&wa, dr, dc);
- }
- }
- 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;
-}
-
-/*
- * extracts streams for threshold, accumulation is provided
- */
-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;
- DCELL value, valued, *wat_nbr;
- struct Cell_head window;
- 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;
- int is_worked;
- char aspect;
- double max_acc;
- int edge, flat;
- int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
- int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
- int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
- int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
- /* sides */
- /*
- * | 7 | 1 | 4 |
- * | 2 | | 3 |
- * | 5 | 0 | 6 |
- */
- unsigned int 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 stream nodes */
- n_alloc_nodes = stream_node_step;
- stream_node =
- (struct snode *)G_malloc(n_alloc_nodes * sizeof(struct snode));
- n_stream_nodes = 0;
-
- /* init outlet nodes */
- n_alloc_outlets = stream_node_step;
- outlets =
- (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);
-
- for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r, c (r_nbr, c_nbr) for neighbours */
- r_nbr = nextdr[ct_dir];
- c_nbr = nextdc[ct_dir];
- /* account for rare cases when ns_res != ew_res */
- dy = abs(r_nbr) * window.ns_res;
- dx = abs(c_nbr) * window.ew_res;
- if (ct_dir < 4)
- dist_to_nbr[ct_dir] = dx + dy;
- else
- dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
- }
-
- diag = sqrt(2);
-
- workedon = 0;
-
- /* extract streams */
- 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;
-
- bseg_get(&asp, &aspect, 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 = 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 =
- (POINT *)G_realloc(outlets,
- n_alloc_outlets *
- sizeof(POINT));
- }
- outlets[n_outlets].r = r;
- outlets[n_outlets].c = c;
- n_outlets++;
- }
-
- if (aspect == 0) {
- /* can only happen with real depressions */
- if (!have_depressions)
- G_fatal_error(_("Bug in stream extraction"));
- G_debug(2, "bottom of real depression");
- }
- continue;
- }
-
- if (aspect) {
- dr = r + asp_r[abs((int)aspect)];
- dc = c + asp_c[abs((int)aspect)];
- }
- else {
- /* can only happen with real depressions,
- * but should not get to here */
- dr = r;
- dc = c;
- }
-
- r_nbr = r_max = dr;
- c_nbr = c_max = dc;
-
- seg_get(&watalt, (char *)&wa, r, c);
- value = wa.wat;
-
- /**********************************/
- /* find main drainage direction */
- /**********************************/
-
- max_acc = -1;
- max_side = np_side = -1;
- mfd_cells = 0;
- stream_cells = 0;
- swale_cells = 0;
- astar_not_set = 1;
- ele_val = wa.ele;
- edge = 0;
- flat = 1;
- /* find main drainage direction */
- for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r_nbr, c_nbr for neighbours */
- r_nbr = r + nextdr[ct_dir];
- c_nbr = c + nextdc[ct_dir];
- 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;
-
- 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 = FLAG_GET(flag_nbr[ct_dir], STREAMFLAG);
- if (is_swale)
- swale_cells++;
-
- /* check for stream cells */
- valued = fabs(wat_nbr[ct_dir]);
- /* check all upstream neighbours */
- if (valued >= threshold && ct_dir != np_side &&
- ele_nbr[ct_dir] > ele_val)
- stream_cells++;
-
- is_worked = FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG);
- if (!internal_acc)
- is_worked = is_worked == 0;
-
- if (is_worked == 0) {
- if (ele_nbr[ct_dir] != ele_val)
- flat = 0;
- if (ele_nbr[ct_dir] <= ele_val) {
-
- mfd_cells++;
-
- /* set main drainage direction */
- if (valued >= max_acc) {
- max_acc = valued;
- r_max = r_nbr;
- c_max = c_nbr;
- max_side = ct_dir;
- }
-
- if (dr == r_nbr && dc == c_nbr) {
- astar_not_set = 0;
- }
- }
- }
- else if (ct_dir == np_side && !edge) {
- /* check for consistency with A * path */
- workedon++;
- }
- }
- else
- edge = 1;
- if (edge)
- break;
- }
-
- is_swale = FLAG_GET(this_flag_value, STREAMFLAG);
-
- /* do not continue streams along edges, these are artifacts */
- if (edge) {
- G_debug(3, "edge");
- if (is_swale) {
- G_debug(2, "edge outlet");
- /* add outlet point */
- if (n_outlets >= n_alloc_outlets) {
- n_alloc_outlets += stream_node_step;
- outlets =
- (POINT *)G_realloc(outlets,
- n_alloc_outlets *
- sizeof(POINT));
- }
- outlets[n_outlets].r = r;
- outlets[n_outlets].c = c;
- n_outlets++;
- if (aspect > 0) {
- aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
- 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) {
- 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;
- }
- }
- 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)) {
- aspect = drain[r - r_max + 1][c - c_max + 1];
- bseg_put(&asp, &aspect, r, c);
- }
-
- /**********************/
- /* start new stream */
- /**********************/
-
- /* Montgomery's stream initiation acc * (tan(slope))^mont_exp */
- /* uses whatever unit is accumulation */
- if (mont_exp > 0) {
- if (r_max == r && c_max == c)
- G_warning
- (_("Can't use Montgomery's method, no stream direction found"));
- else {
- slope = (double)(ele_val - ele_nbr[max_side]) / ele_scale;
-
- if (max_side > 3)
- slope /= diag;
-
- value *= pow(fabs(slope), mont_exp);
- }
- }
-
- if (!is_swale && fabs(value) >= threshold && stream_cells < 1 &&
- swale_cells < 1 && !flat) {
- G_debug(2, "start new stream");
- 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;
- stream_node =
- (struct snode *)G_realloc(stream_node,
- n_alloc_nodes *
- sizeof(struct snode));
- }
- stream_node[stream_no].r = r;
- stream_node[stream_no].c = c;
- stream_node[stream_no].id = stream_no;
- stream_node[stream_no].n_trib = 0;
- stream_node[stream_no].n_trib_total = 0;
- stream_node[stream_no].n_alloc = 0;
- stream_node[stream_no].trib = NULL;
- n_stream_nodes++;
-
- /* debug */
- if (n_stream_nodes != stream_no)
- G_warning(_("BUG: stream_no %d and n_stream_nodes %d out of sync"),
- stream_no, n_stream_nodes);
- }
-
- /*********************/
- /* continue stream */
- /*********************/
-
- 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);
-
- if (n_outlets >= n_alloc_outlets) {
- n_alloc_outlets += stream_node_step;
- outlets =
- (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,
- &stream_no);
- }
- }
- }
- G_percent(1, 1, 2);
- if (workedon)
- G_warning(_("MFD: A * path already processed when setting drainage direction: %d of %d cells"),
- workedon, n_points);
-
- G_free(dist_to_nbr);
- 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);
-
- return 1;
-}
Deleted: grass-addons/grass7/raster/r.hydrodem/thin.c
===================================================================
--- grass-addons/grass7/raster/r.hydrodem/thin.c 2012-01-09 17:37:16 UTC (rev 50108)
+++ grass-addons/grass7/raster/r.hydrodem/thin.c 2012-01-09 17:41:56 UTC (rev 50109)
@@ -1,172 +0,0 @@
-#include <stdlib.h>
-#include <grass/gis.h>
-#include <grass/glocale.h>
-#include "local_proto.h"
-
-int thin_seg(int stream_id)
-{
- int thinned = 0;
- int r, c, r_nbr, c_nbr, last_r, last_c;
- 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;
-
- cseg_get(&stream, &curr_stream, r, c);
-
- bseg_get(&asp, &aspect, r, c);
- if (aspect > 0) {
- /* get downstream point */
- 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 */
- 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;
- 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 */
- 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 */
- aspect = drain[r - r_nbr + 1][c - c_nbr + 1];
- bseg_put(&asp, &aspect, r, c);
-
- thinned = 1;
- }
- else {
- /* nothing to eliminate, continue from last point */
- r = last_r;
- c = last_c;
- }
- last_r = r_nbr;
- last_c = c_nbr;
- bseg_get(&asp, &aspect, last_r, last_c);
- }
- }
-
- return thinned;
-}
-
-int thin_streams(void)
-{
- int i, j, r, c, done;
- CELL stream_id;
- int next_node;
- struct sstack
- {
- int stream_id;
- int next_trib;
- } *nodestack;
- int top = 0, stack_step = 1000;
- int n_trib_total;
- int n_thinned = 0;
-
- G_message(_("Thin stream segments..."));
-
- nodestack = (struct sstack *)G_malloc(stack_step * sizeof(struct sstack));
-
- for (i = 0; i < n_outlets; i++) {
- G_percent(i, n_outlets, 2);
- r = outlets[i].r;
- c = outlets[i].c;
- cseg_get(&stream, &stream_id, r, c);
-
- if (stream_id == 0)
- continue;
-
- /* add root node to stack */
- G_debug(2, "add root node");
- top = 0;
- nodestack[top].stream_id = stream_id;
- nodestack[top].next_trib = 0;
-
- /* depth first post order traversal */
- G_debug(2, "traverse");
- while (top >= 0) {
-
- done = 1;
- stream_id = nodestack[top].stream_id;
- G_debug(3, "stream_id %d, top %d", stream_id, top);
- if (nodestack[top].next_trib < stream_node[stream_id].n_trib) {
- /* add to stack */
- G_debug(3, "get next node");
- next_node =
- stream_node[stream_id].trib[nodestack[top].next_trib];
- G_debug(3, "add to stack: next %d, trib %d, n trib %d",
- next_node, nodestack[top].next_trib,
- stream_node[stream_id].n_trib);
- nodestack[top].next_trib++;
- top++;
- if (top >= stack_step) {
- /* need more space */
- stack_step += 1000;
- nodestack =
- (struct sstack *)G_realloc(nodestack,
- stack_step *
- sizeof(struct sstack));
- }
-
- nodestack[top].next_trib = 0;
- nodestack[top].stream_id = next_node;
- done = 0;
- G_debug(3, "go further down");
- }
- if (done) {
- /* thin stream segment */
- G_debug(3, "thin stream segment %d", stream_id);
-
- if (thin_seg(stream_id) == 0)
- G_debug(3, "segment %d not thinned", stream_id);
- else {
- G_debug(3, "segment %d thinned", stream_id);
- n_thinned++;
- }
-
- top--;
- /* count tributaries */
- if (top >= 0) {
- n_trib_total = 0;
- stream_id = nodestack[top].stream_id;
- for (j = 0; j < stream_node[stream_id].n_trib; j++) {
- /* intermediate */
- if (stream_node[stream_node[stream_id].trib[j]].
- n_trib > 0)
- n_trib_total +=
- stream_node[stream_node[stream_id].trib[j]].
- n_trib_total;
- /* start */
- else
- n_trib_total++;
- }
- stream_node[stream_id].n_trib_total = n_trib_total;
- }
- }
- }
- }
- G_percent(n_outlets, n_outlets, 1); /* finish it */
-
- G_free(nodestack);
-
- G_verbose_message(_("%d of %d stream segments were thinned"), n_thinned, n_stream_nodes);
-
- return 1;
-}
More information about the grass-commit
mailing list