[GRASS-SVN] r72896 - grass-addons/grass7/raster/r.accumulate
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Jun 24 10:12:46 PDT 2018
Author: hcho
Date: 2018-06-24 10:12:46 -0700 (Sun, 24 Jun 2018)
New Revision: 72896
Modified:
grass-addons/grass7/raster/r.accumulate/accumulate.c
grass-addons/grass7/raster/r.accumulate/main.c
Log:
r.accumulate: Fix a bug in flow weighting; Add more comments
Modified: grass-addons/grass7/raster/r.accumulate/accumulate.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/accumulate.c 2018-06-24 13:44:57 UTC (rev 72895)
+++ grass-addons/grass7/raster/r.accumulate/accumulate.c 2018-06-24 17:12:46 UTC (rev 72896)
@@ -11,9 +11,8 @@
#define SW 225
#define W 180
-double
-accumulate(CELL ** dir_buf, RASTER_MAP weight_buf, RASTER_MAP acc_buf,
- char **done, char neg, int row, int col)
+double accumulate(CELL ** dir_buf, RASTER_MAP weight_buf, RASTER_MAP acc_buf,
+ char **done, char neg, int row, int col)
{
static int dir_checks[3][3][2] = {
{{SE, NW}, {S, N}, {SW, NE}},
@@ -20,36 +19,64 @@
{{E, W}, {0, 0}, {W, E}},
{{NE, SW}, {N, S}, {NW, SE}}
};
- int rows = weight_buf.rows, cols = weight_buf.cols;
+ int rows = acc_buf.rows, cols = acc_buf.cols;
int i, j;
char incomplete = 0;
double acc;
- if (done[row][col])
- return fabs(get(acc_buf, row, col));
+ /* if the current cell has been calculated, just return its accumulation so
+ * that downstream cells can simply propagate and add it to themselves */
+ if (done[row][col]) {
+ acc = get(acc_buf, row, col);
+ /* for negative accumulation, always return its absolute value;
+ * otherwise return it as is */
+ return neg && acc < 0 ? -acc : acc;
+ }
+
+ /* if a weight map is specified (no negative accumulation is implied), use
+ * the weight value at the current cell; otherwise use 1 */
if (weight_buf.map.v)
acc = get(weight_buf, row, col);
else
acc = 1.0;
+ /* loop through all neighbor cells and see if any of them drains into the
+ * current cell (are there upstream cells?) */
for (i = -1; i <= 1; i++) {
+ /* if a neighbor cell is outside the computational region, its
+ * downstream accumulation is incomplete */
if (row + i < 0 || row + i >= rows) {
incomplete = 1;
continue;
}
+
for (j = -1; j <= 1; j++) {
+ /* skip the current cell */
if (i == 0 && j == 0)
continue;
+
+ /* if a neighbor cell is outside the computational region, its
+ * downstream accumulation is incomplete */
if (col + j < 0 || col + j >= cols) {
incomplete = 1;
continue;
}
+
+ /* if a neighbor cell drains into the current cell and the current
+ * cell doesn't flow back into the same neighbor cell (no flow
+ * loop), trace and recursively accumulate upstream cells */
if (dir_buf[row + i][col + j] == dir_checks[i + 1][j + 1][0] &&
dir_buf[row][col] != dir_checks[i + 1][j + 1][1]) {
+ /* for negative accumulation, accumulate() always returns a
+ * positive value, so acc is always positive (cell count);
+ * otherwise, acc is weighted accumulation */
acc +=
accumulate(dir_buf, weight_buf, acc_buf, done, neg,
row + i, col + j);
+
+ /* if the neighbor cell is incomplete, the current cell also
+ * becomes incomplete */
if (done[row + i][col + j] == 2)
incomplete = 1;
}
@@ -56,8 +83,15 @@
}
}
+ /* if negative accumulation is desired and the current cell is incomplete,
+ * use a negative cell count without weighting; otherwise use accumulation
+ * as is (cell count or weighted accumulation, which can be negative) */
set(acc_buf, row, col, neg && incomplete ? -acc : acc);
+
+ /* the current cell is done; 1 for no likely underestimates and 2 for
+ * likely underestimates */
done[row][col] = 1 + incomplete;
+ /* for negative accumulation, acc is already positive */
return acc;
}
Modified: grass-addons/grass7/raster/r.accumulate/main.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/main.c 2018-06-24 13:44:57 UTC (rev 72895)
+++ grass-addons/grass7/raster/r.accumulate/main.c 2018-06-24 17:12:46 UTC (rev 72896)
@@ -90,6 +90,8 @@
flag.neg->label =
_("Use negative flow accumulation for likely underestimates");
+ /* weighting doesn't support negative accumulation because weights
+ * themselves can be negative */
G_option_exclusive(opt.weight, flag.neg, NULL);
if (G_parser(argc, argv))
@@ -139,23 +141,12 @@
opt.format->answer);
/* end of r.path */
- if (weight_name) {
- weight_fd = Rast_open_old(weight_name, "");
- weight_buf.type = Rast_get_map_type(weight_fd);
- }
- else {
- weight_fd = -1;
- weight_buf.type = CELL_TYPE;
- }
-
- acc_buf.type = weight_buf.type;
- acc_fd = Rast_open_new(acc_name, acc_buf.type);
-
neg = flag.neg->answer;
rows = Rast_window_rows();
cols = Rast_window_cols();
+ /* initialize the done array and read the direction map */
done = G_malloc(rows * sizeof(char *));
dir_buf = G_malloc(rows * sizeof(CELL *));
for (row = 0; row < rows; row++) {
@@ -168,7 +159,12 @@
}
}
- if (weight_fd >= 0) {
+ /* optionally, read a weight map */
+ if (weight_name) {
+ weight_fd = Rast_open_old(weight_name, "");
+ acc_buf.type = weight_buf.type = Rast_get_map_type(weight_fd);
+ weight_buf.rows = rows;
+ weight_buf.cols = cols;
weight_buf.map.v = (void **)G_malloc(rows * sizeof(void *));
for (row = 0; row < rows; row++) {
weight_buf.map.v[row] =
@@ -177,42 +173,53 @@
weight_buf.type);
}
}
- else
+ else {
+ weight_fd = -1;
weight_buf.map.v = NULL;
- weight_buf.rows = rows;
- weight_buf.cols = cols;
+ acc_buf.type = CELL_TYPE;
+ }
+ /* create output buffer */
+ acc_buf.rows = rows;
+ acc_buf.cols = cols;
acc_buf.map.v = (void **)G_malloc(rows * sizeof(void *));
for (row = 0; row < rows; row++)
acc_buf.map.v[row] = (void *)Rast_allocate_buf(acc_buf.type);
- acc_buf.rows = rows;
- acc_buf.cols = cols;
+ /* create a new output map */
+ acc_fd = Rast_open_new(acc_name, acc_buf.type);
+
+ /* accumulate flows */
for (row = 0; row < rows; row++) {
for (col = 0; col < cols; col++)
accumulate(dir_buf, weight_buf, acc_buf, done, neg, row, col);
}
+ /* write out output buffer to the output map */
for (row = 0; row < rows; row++)
Rast_put_row(acc_fd, acc_buf.map.v[row], acc_buf.type);
- G_free(done);
+ /* close all maps */
+ Rast_close(dir_fd);
+ Rast_close(acc_fd);
+ if (weight_fd >= 0)
+ Rast_close(weight_fd);
+
+ /* free buffer memory */
for (row = 0; row < rows; row++) {
+ G_free(done[row]);
G_free(dir_buf[row]);
+ G_free(acc_buf.map.v[row]);
if (weight_fd >= 0)
G_free(weight_buf.map.v[row]);
- G_free(acc_buf.map.v[row]);
}
+ G_free(done);
G_free(dir_buf);
+ G_free(acc_buf.map.v);
if (weight_fd >= 0)
G_free(weight_buf.map.v);
- G_free(acc_buf.map.v);
- Rast_close(dir_fd);
- if (weight_fd >= 0)
- Rast_close(weight_fd);
- Rast_close(acc_fd);
-
+ /* write history */
Rast_put_cell_title(acc_name,
weight_name ? "Weighted flow accumulation" :
(neg ? "Flow accumulation with likely underestimates"
More information about the grass-commit
mailing list