[GRASS-SVN] r34833 - in grass/trunk/raster/r.watershed: ram seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Dec 12 08:22:44 EST 2008
Author: glynn
Date: 2008-12-12 08:22:43 -0500 (Fri, 12 Dec 2008)
New Revision: 34833
Modified:
grass/trunk/raster/r.watershed/ram/do_cum.c
grass/trunk/raster/r.watershed/ram/init_vars.c
grass/trunk/raster/r.watershed/seg/do_cum.c
grass/trunk/raster/r.watershed/seg/init_vars.c
Log:
Patch from Markus Metz (ticket #398)
Modified: grass/trunk/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_cum.c 2008-12-11 20:41:24 UTC (rev 34832)
+++ grass/trunk/raster/r.watershed/ram/do_cum.c 2008-12-12 13:22:43 UTC (rev 34833)
@@ -87,12 +87,13 @@
DCELL value, valued;
int killer, threshold, count;
- int mfd_cells, astar_not_set;
+ /* MFD */
+ int mfd_cells, stream_cells, swale_cells, astar_not_set;
double *dist_to_nbr, *weight, sum_weight, max_weight;
- int r_nbr, c_nbr, ct_dir, np_side, in_val;
+ int r_nbr, c_nbr, r_max, c_max, r_swale, c_swale, ct_dir, np_side;
double dx, dy;
- CELL ele, ele_nbr;
- double prop;
+ CELL ele, ele_nbr, aspect, is_worked;
+ double prop, max_acc, max_swale;
int workedon;
G_message(_("SECTION 3: Accumulating Surface Flow with MFD"));
@@ -135,10 +136,6 @@
value = wat[SEG_INDEX(wat_seg, r, c)];
valued = wat[SEG_INDEX(wat_seg, dr, dc)];
- /* disabled for MFD */
- /* if (ABS(value) > threshold)
- FLAG_SET(swale, r, c); */
-
/* start MFD */
/* get weights */
@@ -157,8 +154,8 @@
/* check that neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols) {
- in_val = FLAG_GET(worked, r_nbr, c_nbr);
- if (in_val == 0) {
+ is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ if (is_worked == 0) {
ele_nbr = alt[SEG_INDEX(alt_seg, r_nbr, c_nbr)];
if (ele_nbr < ele) {
weight[ct_dir] =
@@ -202,6 +199,15 @@
}
/* set flow accumulation for neighbours */
+ stream_cells = 0;
+ swale_cells = 0;
+ max_acc = -1;
+ max_swale = -1;
+ r_max = dr;
+ c_max = dc;
+ r_swale = dr;
+ c_swale = dc;
+
if (mfd_cells > 1) {
prop = 0.0;
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -212,8 +218,9 @@
/* check that neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols && weight[ct_dir] > -0.5) {
- in_val = FLAG_GET(worked, r_nbr, c_nbr);
- if (in_val == 0) {
+ is_worked = FLAG_GET(worked, r_nbr, c_nbr);
+ is_swale = FLAG_GET(swale, r_nbr, c_nbr);
+ if (is_worked == 0) {
weight[ct_dir] = weight[ct_dir] / sum_weight;
prop += weight[ct_dir]; /* check everything sums up to 1.0 */
@@ -232,10 +239,26 @@
valued = value * weight[ct_dir] - valued;
}
wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)] = valued;
+
+ if ((ABS(valued) + 0.5) >= threshold)
+ stream_cells++;
+ if (ABS(valued) > max_acc) {
+ max_acc = ABS(valued);
+ r_max = r_nbr;
+ c_max = c_nbr;
+ }
+ /* assist in stream merging */
+ if (is_swale && ABS(valued) > max_swale) {
+ max_swale = ABS(valued);
+ r_swale = r_nbr;
+ c_swale = c_nbr;
+ }
}
else if (ct_dir == np_side) {
workedon++; /* check for consistency with A * path */
}
+ if (is_swale)
+ swale_cells++;
}
}
if (ABS(prop - 1.0) > 5E-6f) {
@@ -243,6 +266,10 @@
prop);
}
valued = wat[SEG_INDEX(wat_seg, dr, dc)];
+ if (max_swale > -0.5) {
+ r_max = r_swale;
+ c_max = c_swale;
+ }
}
if (mfd_cells < 2) {
@@ -264,15 +291,28 @@
/* MFD finished */
valued = ABS(valued) + 0.5;
-
is_swale = FLAG_GET(swale, r, c);
- if (is_swale || ((int)valued) >= threshold) {
+ /* start new stream: only works with A * path */
+ if (!is_swale && (int)valued >= threshold && stream_cells < 2 &&
+ swale_cells < 1) {
FLAG_SET(swale, dr, dc);
}
+ /* continue stream */
+ if (is_swale) {
+ FLAG_SET(swale, r_max, c_max);
+ }
else {
if (er_flag && !is_swale)
- slope_length(r, c, dr, dc);
+ slope_length(r, c, r_max, c_max);
}
+ /* update asp */
+ if (dr != r_max || dc != c_max) {
+ aspect = drain[r - r_max + 1][c - c_max + 1];
+ if (asp[SEG_INDEX(asp_seg, r, c)] < 0)
+ aspect = -aspect;
+ asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+
+ }
FLAG_SET(worked, r, c);
}
}
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2008-12-11 20:41:24 UTC (rev 34832)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2008-12-12 13:22:43 UTC (rev 34833)
@@ -129,11 +129,19 @@
G_fatal_error(_("unable to open elevation map layer"));
}
+ swale = flag_create(nrows, ncols);
+ in_list = flag_create(nrows, ncols);
+ worked = flag_create(nrows, ncols);
+
for (r = 0; r < nrows; r++) {
G_get_c_raster_row(fd, buf, r);
for (c = 0; c < ncols; c++) {
index = SEG_INDEX(alt_seg, r, c);
alt[index] = r_h[index] = buf[c];
+ /* all flags need to be manually set to zero */
+ flag_unset(swale, r, c);
+ flag_unset(in_list, r, c);
+ flag_unset(worked, r, c);
}
}
G_close_cell(fd);
@@ -176,7 +184,7 @@
}
G_close_cell(fd);
}
- swale = flag_create(nrows, ncols);
+
if (ob_flag) {
fd = G_open_cell_old(ob_name, "");
if (fd < 0) {
@@ -197,8 +205,7 @@
G_fatal_error(_("unable to open rill map layer"));
}
}
- in_list = flag_create(nrows, ncols);
- worked = flag_create(nrows, ncols);
+
MASK_flag = 0;
do_points = nrows * ncols;
if (NULL != G_find_file("cell", "MASK", G_mapset())) {
@@ -237,7 +244,7 @@
sizeof(double));
}
- /* heap_index will track astar_pts in the binary min-heap */
+ /* heap_index will track astar_pts in ternary min-heap */
/* heap_index is one-based */
heap_index = (int *)G_malloc((do_points + 1) * sizeof(int));
Modified: grass/trunk/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_cum.c 2008-12-11 20:41:24 UTC (rev 34832)
+++ grass/trunk/raster/r.watershed/seg/do_cum.c 2008-12-12 13:22:43 UTC (rev 34833)
@@ -92,12 +92,13 @@
POINT point;
int killer, threshold, count;
- int mfd_cells, astar_not_set;
+ /* MFD */
+ int mfd_cells, stream_cells, swale_cells, astar_not_set;
double *dist_to_nbr, *weight, sum_weight, max_weight;
- int r_nbr, c_nbr, ct_dir, np_side, in_val;
+ int r_nbr, c_nbr, r_max, c_max, r_swale, c_swale, ct_dir, np_side;
double dx, dy;
- CELL ele, ele_nbr;
- double prop;
+ CELL ele, ele_nbr, aspect, asp_val, is_worked;
+ double prop, max_acc, max_swale;
int workedon;
G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
@@ -168,8 +169,8 @@
/* check that neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols) {
- bseg_get(&worked, &in_val, r_nbr, c_nbr);
- if (in_val == 0) {
+ bseg_get(&worked, &is_worked, r_nbr, c_nbr);
+ if (is_worked == 0) {
cseg_get(&alt, &ele_nbr, r_nbr, c_nbr);
if (ele_nbr < ele) {
weight[ct_dir] =
@@ -213,6 +214,15 @@
}
/* set flow accumulation for neighbours */
+ stream_cells = 0;
+ swale_cells = 0;
+ max_acc = -1;
+ max_swale = -1;
+ r_max = dr;
+ c_max = dc;
+ r_swale = dr;
+ c_swale = dc;
+
if (mfd_cells > 1) {
prop = 0.0;
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -223,8 +233,9 @@
/* check that neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols && weight[ct_dir] > -0.5) {
- bseg_get(&worked, &in_val, r_nbr, c_nbr);
- if (in_val == 0) {
+ bseg_get(&worked, &is_worked, r_nbr, c_nbr);
+ bseg_get(&swale, &is_swale, r_nbr, c_nbr);
+ if (is_worked == 0) {
weight[ct_dir] = weight[ct_dir] / sum_weight;
prop += weight[ct_dir]; /* check everything sums up to 1.0 */
@@ -243,6 +254,20 @@
valued = value * weight[ct_dir] - valued;
}
dseg_put(&wat, &valued, r_nbr, c_nbr);
+
+ if ((ABS(valued) + 0.5) >= threshold)
+ stream_cells++;
+ if (ABS(valued) > max_acc) {
+ max_acc = ABS(valued);
+ r_max = r_nbr;
+ c_max = c_nbr;
+ }
+ /* assist in stream merging */
+ if (is_swale && ABS(valued) > max_swale) {
+ max_swale = ABS(valued);
+ r_swale = r_nbr;
+ c_swale = c_nbr;
+ }
}
else if (ct_dir == np_side) {
workedon++; /* check for consistency with A * path */
@@ -254,6 +279,10 @@
prop);
}
dseg_get(&wat, &valued, dr, dc);
+ if (max_swale > -0.5) {
+ r_max = r_swale;
+ c_max = c_swale;
+ }
}
if (mfd_cells < 2) {
@@ -275,15 +304,29 @@
/* MFD finished */
valued = ABS(valued) + 0.5;
-
bseg_get(&swale, &is_swale, r, c);
- if (is_swale || (int)valued >= threshold) {
+ /* start new stream: only works with A * path */
+ if (!is_swale && (int)valued >= threshold && stream_cells < 2 &&
+ swale_cells < 1) {
bseg_put(&swale, &one, dr, dc);
}
+ /* continue stream */
+ if (is_swale) {
+ bseg_put(&swale, &one, r_max, c_max);
+ }
else {
if (er_flag && !is_swale)
- slope_length(r, c, dr, dc);
+ slope_length(r, c, r_max, c_max);
}
+ /* update asp */
+ if (dr != r_max || dc != c_max) {
+ aspect = drain[r - r_max + 1][c - c_max + 1];
+ cseg_get(&asp, &asp_val, r, c);
+ if (asp_val < 0)
+ aspect = -aspect;
+ cseg_put(&asp, &aspect, r, c);
+
+ }
bseg_put(&worked, &one, r, c);
}
}
Modified: grass/trunk/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/init_vars.c 2008-12-11 20:41:24 UTC (rev 34832)
+++ grass/trunk/raster/r.watershed/seg/init_vars.c 2008-12-12 13:22:43 UTC (rev 34833)
@@ -240,7 +240,7 @@
seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols,
num_open_segs, sizeof(POINT));
- /* heap_index will track astar_pts in the binary min-heap */
+ /* heap_index will track astar_pts in ternary min-heap */
/* heap_index is one-based */
seg_open(&heap_index, 1, do_points + 1, 1, seg_rows * seg_cols,
num_open_segs, sizeof(HEAP));
More information about the grass-commit
mailing list