[GRASS-SVN] r42235 - grass/trunk/raster/r.watershed/ram
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 12 05:33:34 EDT 2010
Author: mmetz
Date: 2010-05-12 05:33:34 -0400 (Wed, 12 May 2010)
New Revision: 42235
Added:
grass/trunk/raster/r.watershed/ram/do_flatarea.c
Modified:
grass/trunk/raster/r.watershed/ram/Gwater.h
grass/trunk/raster/r.watershed/ram/do_astar.c
grass/trunk/raster/r.watershed/ram/init_vars.c
grass/trunk/raster/r.watershed/ram/main.c
Log:
new option to beautify flat areas
Modified: grass/trunk/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/Gwater.h 2010-05-12 08:28:11 UTC (rev 42234)
+++ grass/trunk/raster/r.watershed/ram/Gwater.h 2010-05-12 09:33:34 UTC (rev 42235)
@@ -54,7 +54,7 @@
extern CELL n_basins;
extern OC_STACK *ocs;
extern int ocs_alloced;
-extern FLAG *worked, *in_list, *s_b, *swale;
+extern FLAG *worked, *in_list, *s_b, *swale, *flat_done;
extern RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
extern RAMSEG r_h_seg, dep_seg;
extern RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
@@ -77,7 +77,7 @@
extern char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX], thr_name[8];
extern char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX], sg_name[GNAME_MAX];
extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX], dis_name[GNAME_MAX];
-extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
+extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, flat_flag;
extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
extern char bas_flag, seg_flag, haf_flag, er_flag;
extern char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
@@ -100,6 +100,9 @@
double get_slope(int, int, int, int, CELL, CELL);
int replace(int, int, int, int);
+/* do_flatarea.c */
+int do_flatarea(int, CELL);
+
/* do_cum.c */
int do_cum(void);
int do_cum_mfd(void);
Modified: grass/trunk/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.c 2010-05-12 08:28:11 UTC (rev 42234)
+++ grass/trunk/raster/r.watershed/ram/do_astar.c 2010-05-12 09:33:34 UTC (rev 42235)
@@ -10,7 +10,7 @@
int count;
int upr, upc, r, c, ct_dir;
CELL alt_val, alt_nbr[8];
- CELL is_in_list, is_worked;
+ CELL is_in_list, is_worked, flat_is_done, nbr_flat_is_done;
int index_doer, index_up;
/* sides
* |7|1|4|
@@ -22,6 +22,7 @@
double dx, dy, dist_to_nbr[8], ew_res, ns_res;
double slope[8];
int skip_diag;
+ CELL *alt_bak;
G_message(_("SECTION 2: A* Search."));
@@ -44,6 +45,27 @@
first_astar = heap_index[1];
first_cum = do_points;
+ if (flat_flag) {
+ alt_bak =
+ (CELL *) G_malloc(sizeof(CELL) * size_array(&alt_seg, nrows, ncols));
+ flat_done = flag_create(nrows, ncols);
+ flat_is_done = 0;
+
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ index_doer = SEG_INDEX(alt_seg, r, c);
+ alt_bak[index_doer] = alt[index_doer];
+
+ flag_unset(flat_done, r, c);
+ }
+ }
+ }
+ else {
+ alt_bak = NULL;
+ flat_done = NULL;
+ flat_is_done = 1;
+ }
+
/* A* Search: search uphill, get downhill paths */
while (heap_size > 0) {
G_percent(count++, do_points, 1);
@@ -64,7 +86,9 @@
alt_val = alt[index_doer];
- FLAG_SET(worked, r, c);
+ if (flat_flag) {
+ flat_is_done = FLAG_GET(flat_done, r, c);
+ }
/* check all neighbours, breadth first search */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -78,9 +102,31 @@
is_in_list = FLAG_GET(in_list, upr, upc);
is_worked = FLAG_GET(worked, upr, upc);
skip_diag = 0;
+
+ alt_nbr[ct_dir] = alt[index_up];
+ if (flat_flag && !is_worked) {
+ alt_val = alt[index_doer];
+ if (!flat_is_done && alt_nbr[ct_dir] == alt_val) {
+ do_flatarea(index_doer, alt_val);
+ alt_nbr[ct_dir] = alt[index_up];
+ flat_is_done = 1;
+ nbr_flat_is_done = 1;
+ }
+ nbr_flat_is_done = FLAG_GET(flat_done, upr, upc);
+ if (!nbr_flat_is_done) {
+ /* use original ele values */
+ alt_val = alt_bak[index_doer];
+ alt_nbr[ct_dir] = alt_bak[index_up];
+ }
+ else {
+ /* use modified ele values */
+ alt_val = alt[index_doer];
+ alt_nbr[ct_dir] = alt[index_up];
+ }
+ }
+
/* avoid diagonal flow direction bias */
if (!is_worked) {
- alt_nbr[ct_dir] = alt[index_up];
slope[ct_dir] =
get_slope2(alt_val, alt_nbr[ct_dir],
dist_to_nbr[ct_dir]);
@@ -121,6 +167,7 @@
}
} /* end if in region */
} /* end sides */
+ FLAG_SET(worked, r, c);
}
G_percent(count, do_points, 1); /* finish it */
if (mfd == 0)
@@ -129,6 +176,17 @@
flag_destroy(in_list);
G_free(heap_index);
+ if (flat_flag) {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ index_doer = SEG_INDEX(alt_seg, r, c);
+ alt[index_doer] = alt_bak[index_doer];
+ }
+ }
+ G_free(alt_bak);
+ flag_destroy(flat_done);
+ }
+
return 0;
}
Added: grass/trunk/raster/r.watershed/ram/do_flatarea.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_flatarea.c (rev 0)
+++ grass/trunk/raster/r.watershed/ram/do_flatarea.c 2010-05-12 09:33:34 UTC (rev 42235)
@@ -0,0 +1,436 @@
+/**********************************************************************
+ *
+ * flat area beautification after Garbrecht and Martz (1997)
+ *
+ * Garbrecht, J. & Martz, L. W. (1997)
+ * The assignment of drainage direction over flat surfaces in raster
+ * digital elevation models. J. Hydrol 193, 204–213.
+ *
+ * the method is modifed for speed, only one pass is necessary to get
+ * the gradient away from higher terrain
+ *
+ *********************************************************************/
+
+
+#include <limits.h>
+#include <assert.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "Gwater.h"
+#include "do_astar.h"
+#include "rbtree.h"
+
+struct pq_node
+{
+ int idx;
+ struct pq_node *next;
+};
+
+struct pq
+{
+ struct pq_node *first, *last;
+ int size;
+};
+
+struct pq *pq_create(void)
+{
+ struct pq *q = G_malloc(sizeof(struct pq));
+
+ q->first = G_malloc(sizeof(struct pq_node));
+ q->first->next = NULL;
+ q->first->idx = -1;
+ q->last = q->first;
+ q->size = 0;
+
+ return q;
+}
+
+/* dummy end must always be allocated and empty */
+int pq_add(int idx, struct pq *q)
+{
+ assert(q->last);
+ if (q->last->idx != -1)
+ G_fatal_error("idx is %d", q->last->idx);
+
+ q->last->idx = idx;
+ if (q->last->next == NULL) {
+ struct pq_node *n = (struct pq_node *) G_malloc(sizeof(struct pq_node));
+ n->next = NULL;
+ n->idx = -1;
+ q->last->next = n;
+ q->last = q->last->next;
+ }
+ else {
+ q->last = q->last->next;
+ assert(NULL);
+ }
+
+ assert(q->last != q->last->next);
+ assert(q->first != q->last);
+ assert(q->last->idx == -1);
+ assert(q->last->next == NULL);
+ q->size++;
+
+ return 0;
+}
+
+int pq_drop(struct pq *q)
+{
+ int idx = q->first->idx;
+ struct pq_node *n = q->first;
+
+ assert(q->first != q->first->next);
+ assert(q->size);
+ q->size--;
+
+ q->first = q->first->next;
+ assert(q->first);
+ assert(q->first != q->first->next);
+ assert(n != q->first);
+ G_free(n);
+
+ return idx;
+}
+
+int pq_destroy(struct pq *q)
+{
+ struct pq_node *delme;
+
+ while (q->first) {
+ delme = q->first;
+ q->first = q->first->next;
+ G_free(delme);
+ }
+
+ G_free(q);
+
+ return 0;
+}
+
+struct orders {
+ int index, uphill, downhill;
+ char flag;
+};
+
+int cmp_orders(const void *a, const void *b)
+{
+ struct orders *oa = (struct orders *)a;
+ struct orders *ob = (struct orders *)b;
+
+ return (oa->index < ob->index ? -1 : (oa->index > ob->index));
+}
+
+int do_flatarea(int index, CELL ele)
+{
+ int upr, upc, r, c, ct_dir;
+ CELL is_in_list, is_worked, this_in_list;
+ int index_doer, index_up;
+ int n_flat_cells = 0, counter;
+ CELL ele_nbr, min_ele_diff;
+ int uphill_order, downhill_order, max_uphill_order, max_downhill_order;
+ int last_order;
+
+ struct pq *up_pq = pq_create();
+ struct pq *down_pq = pq_create();
+
+ struct orders inc_order, *order_found, *nbr_order_found;
+ struct RB_TREE *order_tree = rbtree_create(cmp_orders, sizeof(struct orders));
+
+ pq_add(index, down_pq);
+ inc_order.downhill = -1;
+ inc_order.uphill = 0;
+ inc_order.index = index;
+ inc_order.flag = 1;
+ rbtree_insert(order_tree, &inc_order);
+
+ n_flat_cells++;
+
+ min_ele_diff = INT_MAX;
+ max_uphill_order = max_downhill_order = 0;
+
+ /* get uphill start points */
+ counter = 0;
+ while (down_pq->size) {
+ if ((index_doer = pq_drop(down_pq)) == -1)
+ G_fatal_error("get start points: no more points in down queue");
+
+ seg_index_rc(alt_seg, index_doer, &r, &c);
+
+ FLAG_SET(flat_done, r, c);
+
+ /* check all neighbours, breadth first search */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for this neighbour */
+ upr = r + nextdr[ct_dir];
+ upc = c + nextdc[ct_dir];
+ /* check if r, c are within region */
+ if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
+ index_up = SEG_INDEX(alt_seg, upr, upc);
+ is_in_list = FLAG_GET(in_list, upr, upc);
+ is_worked = FLAG_GET(worked, upr, upc);
+ ele_nbr = alt[index_up];
+
+ if (ele_nbr == ele && !is_worked) {
+
+ inc_order.downhill = -1;
+ inc_order.uphill = -1;
+ inc_order.index = index_up;
+ inc_order.flag = 0;
+ /* not yet added to queue */
+ if ((order_found = rbtree_find(order_tree, &inc_order)) == NULL) {
+ n_flat_cells++;
+
+ /* add to down queue if not yet in there */
+ pq_add(index_up, down_pq);
+
+ /* add to up queue if not yet in there */
+ if (is_in_list) {
+ pq_add(index_up, up_pq);
+ /* set uphill order to 0 */
+ inc_order.uphill = 0;
+ inc_order.flag = 1;
+ counter++;
+ }
+ rbtree_insert(order_tree, &inc_order);
+ }
+ }
+ }
+ }
+ }
+
+ if (n_flat_cells < 5) {
+ /* clean up */
+ pq_destroy(up_pq);
+ pq_destroy(down_pq);
+ rbtree_destroy(order_tree);
+
+ return 1;
+ }
+
+ G_debug(2, "%d flat cells, %d cells in tree, %d start cells", n_flat_cells, (int)order_tree->count, counter);
+
+ pq_destroy(down_pq);
+ down_pq = pq_create();
+
+ /* got uphill start points, do uphill correction */
+ counter = 0;
+ uphill_order = 1;
+ while (up_pq->size) {
+ int is_in_down_queue = 0;
+
+ if ((index_doer = pq_drop(up_pq)) == -1)
+ G_fatal_error("uphill order: no more points in up queue");
+
+ seg_index_rc(alt_seg, index_doer, &r, &c);
+ this_in_list = FLAG_GET(in_list, r, c);
+
+ /* get uphill order for this point */
+ inc_order.index = index_doer;
+ if ((order_found = rbtree_find(order_tree, &inc_order)) == NULL)
+ G_fatal_error(_("flat cell escaped for uphill correction"));
+
+ last_order = uphill_order - 1;
+ uphill_order = order_found->uphill;
+
+ if (last_order > uphill_order)
+ G_warning("queue error: last uphill order %d > current uphill order %d", last_order, uphill_order);
+
+ /* debug */
+ if (uphill_order == -1)
+ G_fatal_error(_("uphill order not set"));
+
+ if (max_uphill_order < uphill_order)
+ max_uphill_order = uphill_order;
+
+ uphill_order++;
+ counter++;
+
+ /* check all neighbours, breadth first search */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for this neighbour */
+ upr = r + nextdr[ct_dir];
+ upc = c + nextdc[ct_dir];
+ /* check if r, c are within region */
+ if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
+ index_up = SEG_INDEX(alt_seg, upr, upc);
+ is_in_list = FLAG_GET(in_list, upr, upc);
+ is_worked = FLAG_GET(worked, upr, upc);
+ ele_nbr = alt[index_up];
+
+ /* all cells that are in_list should have been added previously as uphill start points */
+ if (ele_nbr == ele && !is_worked) {
+
+ inc_order.index = index_up;
+ if ((nbr_order_found = rbtree_find(order_tree, &inc_order)) == NULL) {
+ G_fatal_error(_("flat cell escaped in uphill correction"));
+ }
+
+ /* not yet added to queue */
+ if (nbr_order_found->flag != 1) {
+ if (is_in_list)
+ G_warning("cell should be in queue");
+ /* add to up queue */
+ pq_add(index_up, up_pq);
+ /* set nbr uphill order = current uphill order + 1 */
+ nbr_order_found->uphill = uphill_order;
+ nbr_order_found->flag = 1;
+ }
+ }
+ /* add focus cell to down queue */
+ if (!this_in_list && !is_in_down_queue &&
+ ele_nbr != ele && !is_in_list && !is_worked) {
+ pq_add(index_doer, down_pq);
+ /* set downhill order to 0 */
+ order_found->downhill = 0;
+ is_in_down_queue = 1;
+ }
+ if (ele_nbr > ele && min_ele_diff > ele_nbr - ele)
+ min_ele_diff = ele_nbr - ele;
+ }
+ }
+ }
+ /* debug: all flags should be set to 1 */
+
+ pq_destroy(up_pq);
+ up_pq = pq_create();
+
+ /* got downhill start points, do downhill correction */
+ downhill_order = 1;
+ while (down_pq->size) {
+ if ((index_doer = pq_drop(down_pq)) == -1)
+ G_fatal_error("downhill order: no more points in down queue");
+
+ seg_index_rc(alt_seg, index_doer, &r, &c);
+
+ /* get downhill order for this point */
+ inc_order.index = index_doer;
+ if ((order_found = rbtree_find(order_tree, &inc_order)) == NULL)
+ G_fatal_error(_("flat cell escaped for downhill correction"));
+
+ last_order = downhill_order - 1;
+ downhill_order = order_found->downhill;
+
+ if (last_order > downhill_order)
+ G_warning("queue error: last downhill order %d > current downhill order %d", last_order, downhill_order);
+
+ /* debug */
+ if (downhill_order == -1)
+ G_fatal_error(_("downhill order: downhill order not set"));
+
+ downhill_order++;
+
+ /* check all neighbours, breadth first search */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for this neighbour */
+ upr = r + nextdr[ct_dir];
+ upc = c + nextdc[ct_dir];
+ /* check if r, c are within region */
+ if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
+ index_up = SEG_INDEX(alt_seg, upr, upc);
+ is_in_list = FLAG_GET(in_list, upr, upc);
+ is_worked = FLAG_GET(worked, upr, upc);
+ ele_nbr = alt[index_up];
+
+ if (ele_nbr == ele && !is_worked) {
+
+ inc_order.index = index_up;
+ if ((nbr_order_found = rbtree_find(order_tree, &inc_order)) == NULL)
+ G_fatal_error(_("flat cell escaped in downhill correction"));
+
+ /* not yet added to queue */
+ if (nbr_order_found->downhill == -1) {
+
+ /* add to down queue */
+ pq_add(index_up, down_pq);
+ /* set nbr downhill order = current downhill order + 1 */
+ nbr_order_found->downhill = downhill_order;
+
+ if (max_downhill_order < downhill_order)
+ max_downhill_order = downhill_order;
+
+
+ /* add to up queue */
+ if (is_in_list) {
+ pq_add(index_up, up_pq);
+ /* unset flag */
+ nbr_order_found->flag = 0;
+ //nbr_order_found->downhill = 0;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /* got uphill and downhill order, adjust ele */
+
+ /* increment: ele += uphill_order + max_downhill_order - downhill_order */
+ /* decrement: ele += uphill_order - max_uphill_order - downhill_order */
+ //max_downhill_order++;
+
+ while (up_pq->size) {
+ if ((index_doer = pq_drop(up_pq)) == -1)
+ G_fatal_error("no more points in up queue");
+
+ seg_index_rc(alt_seg, index_doer, &r, &c);
+ this_in_list = FLAG_GET(in_list, r, c);
+
+ /* get uphill and downhill order for this point */
+ inc_order.index = index_doer;
+ if ((order_found = rbtree_find(order_tree, &inc_order)) == NULL)
+ G_fatal_error(_("flat cell escaped for adjustment"));
+
+ uphill_order = order_found->uphill;
+ downhill_order = order_found->downhill;
+
+ /* debug */
+ if (uphill_order == -1)
+ G_fatal_error(_("adjustment: uphill order not set"));
+ if (downhill_order == -1)
+ G_fatal_error(_("adjustment: downhill order not set"));
+
+ /* increment */
+ if (this_in_list) {
+ downhill_order = max_downhill_order;
+ uphill_order = 0;
+ }
+ alt[index_doer] += uphill_order + max_downhill_order - downhill_order;
+
+ /* check all neighbours, breadth first search */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for this neighbour */
+ upr = r + nextdr[ct_dir];
+ upc = c + nextdc[ct_dir];
+ /* check if r, c are within region */
+ if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
+ index_up = SEG_INDEX(alt_seg, upr, upc);
+ is_in_list = FLAG_GET(in_list, upr, upc);
+ is_worked = FLAG_GET(worked, upr, upc);
+ ele_nbr = alt[index_up];
+
+ if (ele_nbr == ele && !is_worked) {
+
+ inc_order.index = index_up;
+ if ((nbr_order_found = rbtree_find(order_tree, &inc_order)) == NULL)
+ G_fatal_error(_("flat cell escaped in adjustment"));
+
+ /* not yet added to queue */
+ if (nbr_order_found->flag != 0) {
+ if (is_in_list)
+ G_warning("adjustment: in_list cell should be in queue");
+ /* add to up queue */
+ pq_add(index_up, up_pq);
+ nbr_order_found->flag = 0;
+ }
+ }
+ }
+ }
+ }
+
+ /* clean up */
+ pq_destroy(up_pq);
+ pq_destroy(down_pq);
+ rbtree_destroy(order_tree);
+
+ return 0;
+}
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2010-05-12 08:28:11 UTC (rev 42234)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2010-05-12 09:33:34 UTC (rev 42235)
@@ -33,6 +33,7 @@
mfd = 1;
c_fac = 5;
abs_acc = 0;
+ flat_flag = 0;
ele_scale = 1;
for (r = 1; r < argc; r++) {
@@ -82,6 +83,8 @@
mfd = 0;
else if (strcmp(argv[r], "-a") == 0)
abs_acc = 1;
+ else if (strcmp(argv[r], "-b") == 0)
+ flat_flag = 1;
else
usage(argv[0]);
}
@@ -155,7 +158,7 @@
elebuf = Rast_allocate_buf(ele_map_type);
if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
- ele_scale = 1000; /* should be enough to do the trick */
+ ele_scale = 10000; /* should be enough to do the trick */
/* read elevation input and mark NULL/masked cells */
/* intialize accumulation and drainage direction */
Modified: grass/trunk/raster/r.watershed/ram/main.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/main.c 2010-05-12 08:28:11 UTC (rev 42234)
+++ grass/trunk/raster/r.watershed/ram/main.c 2010-05-12 09:33:34 UTC (rev 42235)
@@ -33,7 +33,7 @@
CELL n_basins;
OC_STACK *ocs;
int ocs_alloced;
-FLAG *worked, *in_list, *s_b, *swale;
+FLAG *worked, *in_list, *s_b, *swale, *flat_done;
RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
RAMSEG r_h_seg, dep_seg;
RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
@@ -59,7 +59,7 @@
sg_name[GNAME_MAX];
char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX],
dis_name[GNAME_MAX];
-char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
+char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, flat_flag;
char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
char bas_flag, seg_flag, haf_flag, er_flag;
char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
More information about the grass-commit
mailing list