[GRASS-SVN] r72888 - grass-addons/grass7/raster/r.accumulate
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jun 23 09:58:58 PDT 2018
Author: hcho
Date: 2018-06-23 09:58:58 -0700 (Sat, 23 Jun 2018)
New Revision: 72888
Modified:
grass-addons/grass7/raster/r.accumulate/accumulate.c
grass-addons/grass7/raster/r.accumulate/global.h
grass-addons/grass7/raster/r.accumulate/main.c
grass-addons/grass7/raster/r.accumulate/r.accumulate.html
Log:
r.accumulate: Add -n flag for negative flow accumulation (r.watershed's default)
Modified: grass-addons/grass7/raster/r.accumulate/accumulate.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/accumulate.c 2018-06-23 14:43:04 UTC (rev 72887)
+++ grass-addons/grass7/raster/r.accumulate/accumulate.c 2018-06-23 16:58:58 UTC (rev 72888)
@@ -1,9 +1,19 @@
+#include <math.h>
#include <grass/raster.h>
#include "global.h"
+#define NW 135
+#define N 90
+#define NE 45
+#define E 360
+#define SE 315
+#define S 270
+#define SW 225
+#define W 180
+
double
accumulate(CELL ** dir_buf, RASTER_MAP weight_buf, RASTER_MAP acc_buf,
- char **done, int row, int col)
+ char **done, char neg, int row, int col)
{
static int dir_checks[3][3][2] = {
{{SE, NW}, {S, N}, {SW, NE}},
@@ -12,10 +22,11 @@
};
int rows = weight_buf.rows, cols = weight_buf.cols;
int i, j;
+ char incomplete = 0;
double acc;
if (done[row][col])
- return get(acc_buf, row, col);
+ return fabs(get(acc_buf, row, col));
if (weight_buf.map.v)
acc = get(weight_buf, row, col);
@@ -23,21 +34,30 @@
acc = 1.0;
for (i = -1; i <= 1; i++) {
- if (row + i < 0 || row + i >= rows)
+ if (row + i < 0 || row + i >= rows) {
+ incomplete = 1;
continue;
+ }
for (j = -1; j <= 1; j++) {
- if (col + j < 0 || col + j >= cols || (i == 0 && j == 0))
+ if (i == 0 && j == 0)
continue;
+ if (col + j < 0 || col + j >= cols) {
+ incomplete = 1;
+ continue;
+ }
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])
+ dir_buf[row][col] != dir_checks[i + 1][j + 1][1]) {
acc +=
- accumulate(dir_buf, weight_buf, acc_buf, done, row + i,
- col + j);
+ accumulate(dir_buf, weight_buf, acc_buf, done, neg,
+ row + i, col + j);
+ if (done[row + i][col + j] == 2)
+ incomplete = 1;
+ }
}
}
- set(acc_buf, row, col, acc);
- done[row][col] = 1;
+ set(acc_buf, row, col, neg && incomplete ? -acc : acc);
+ done[row][col] = 1 + incomplete;
return acc;
}
Modified: grass-addons/grass7/raster/r.accumulate/global.h
===================================================================
--- grass-addons/grass7/raster/r.accumulate/global.h 2018-06-23 14:43:04 UTC (rev 72887)
+++ grass-addons/grass7/raster/r.accumulate/global.h 2018-06-23 16:58:58 UTC (rev 72888)
@@ -1,18 +1,5 @@
#include <grass/raster.h>
-#define DIR_UNKNOWN 0
-#define DIR_DEG 1
-#define DIR_DEG45 2
-
-#define NW 135
-#define N 90
-#define NE 45
-#define E 360
-#define SE 315
-#define S 270
-#define SW 225
-#define W 180
-
typedef struct
{
RASTER_MAP_TYPE type;
@@ -31,4 +18,4 @@
double get(RASTER_MAP, int, int);
/* accumulate.c */
-double accumulate(CELL **, RASTER_MAP, RASTER_MAP, char **, int, int);
+double accumulate(CELL **, RASTER_MAP, RASTER_MAP, char **, char, int, int);
Modified: grass-addons/grass7/raster/r.accumulate/main.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/main.c 2018-06-23 14:43:04 UTC (rev 72887)
+++ grass-addons/grass7/raster/r.accumulate/main.c 2018-06-23 16:58:58 UTC (rev 72888)
@@ -23,6 +23,10 @@
#include <grass/glocale.h>
#include "global.h"
+#define DIR_UNKNOWN 0
+#define DIR_DEG 1
+#define DIR_DEG45 2
+
int main(int argc, char *argv[])
{
struct GModule *module;
@@ -30,6 +34,10 @@
{
struct Option *dir, *format, *weight, *acc;
} opt;
+ struct
+ {
+ struct Flag *neg;
+ } flag;
char *desc;
char *dir_name, *weight_name, *acc_name;
int dir_fd, weight_fd, acc_fd;
@@ -36,10 +44,12 @@
double dir_format;
struct Range dir_range;
CELL dir_min, dir_max;
+ char neg;
char **done;
CELL **dir_buf;
RASTER_MAP weight_buf, acc_buf;
int rows, cols, row, col;
+ struct History hist;
G_gisinit(argv[0]);
@@ -47,8 +57,7 @@
G_add_keyword(_("raster"));
G_add_keyword(_("hydrology"));
module->description =
- _
- ("Calculates weighted flow accumulation using a flow direction map.");
+ _("Calculates weighted flow accumulation using a flow direction map.");
opt.dir = G_define_standard_option(G_OPT_R_INPUT);
opt.dir->description = _("Name of input direction map");
@@ -63,8 +72,7 @@
G_asprintf(&desc, "auto;%s;degree;%s;45degree;%s",
_("auto-detect direction format"),
_("degrees CCW from East"),
- _
- ("degrees CCW from East divided by 45 (e.g. r.watershed directions)"));
+ _("degrees CCW from East divided by 45 (e.g. r.watershed directions)"));
opt.format->descriptions = desc;
opt.weight = G_define_standard_option(G_OPT_R_INPUT);
@@ -77,6 +85,15 @@
opt.acc->description =
_("Name for output weighted flow accumulation map");
+ flag.neg = G_define_flag();
+ flag.neg->key = 'n';
+ flag.neg->label =
+ _("Use negative flow accumulation for likely underestimates");
+ flag.neg->description =
+ _("See manual for a detailed description of negative flow accumulation");
+
+ G_option_exclusive(opt.weight, flag.neg, NULL);
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -103,24 +120,20 @@
}
else if (strcmp(opt.format->answer, "45degree") == 0) {
if (dir_max > 8)
- G_fatal_error(_
- ("Directional degrees divided by 45 can not be > 8"));
+ G_fatal_error(_("Directional degrees divided by 45 can not be > 8"));
dir_format = DIR_DEG45;
}
else if (strcmp(opt.format->answer, "auto") == 0) {
if (dir_max <= 8) {
dir_format = DIR_DEG45;
- G_important_message(_
- ("Input direction format assumed to be degrees CCW from East divided by 45"));
+ G_important_message(_("Input direction format assumed to be degrees CCW from East divided by 45"));
}
else if (dir_max <= 360) {
dir_format = DIR_DEG;
- G_important_message(_
- ("Input direction format assumed to be degrees CCW from East"));
+ G_important_message(_("Input direction format assumed to be degrees CCW from East"));
}
else
- G_fatal_error(_
- ("Unable to detect format of input direction map <%s>"),
+ G_fatal_error(_("Unable to detect format of input direction map <%s>"),
dir_name);
}
if (dir_format == DIR_UNKNOWN)
@@ -140,6 +153,8 @@
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();
@@ -177,7 +192,7 @@
for (row = 0; row < rows; row++) {
for (col = 0; col < cols; col++)
- accumulate(dir_buf, weight_buf, acc_buf, done, row, col);
+ accumulate(dir_buf, weight_buf, acc_buf, done, neg, row, col);
}
for (row = 0; row < rows; row++)
@@ -200,5 +215,13 @@
Rast_close(weight_fd);
Rast_close(acc_fd);
+ Rast_put_cell_title(acc_name,
+ weight_name ? "Weighted flow accumulation" :
+ (neg ? "Flow accumulation with likely underestimates"
+ : "Flow accumulation"));
+ Rast_short_history(acc_name, "raster", &hist);
+ Rast_command_history(&hist);
+ Rast_write_history(acc_name, &hist);
+
exit(EXIT_SUCCESS);
}
Modified: grass-addons/grass7/raster/r.accumulate/r.accumulate.html
===================================================================
--- grass-addons/grass7/raster/r.accumulate/r.accumulate.html 2018-06-23 14:43:04 UTC (rev 72887)
+++ grass-addons/grass7/raster/r.accumulate/r.accumulate.html 2018-06-23 16:58:58 UTC (rev 72888)
@@ -8,16 +8,25 @@
Unlike <em>r.watershed</em>, <em>r.accumulate</em> does not require the
elevation data to calculate weigited flow accumulation. Instead, this module
only uses a flow direction map to trace and accumulate the amount of flow
-draining through each cell.
+draining through and including each cell.
+<p>With <b>-n</b> flag, the module will count the number of upstream cells plus
+one and convert it to the negative if any upstream cells are likely to receive
+flow from outside the computational region (flow direction edges). Negative
+values indentify cells with likely underestimates because not all upstream
+cells were accounted for. Since raster map <b>weight</b> may contain negative
+flow weights, <b>-n</b> flag is not compatible with <b>weight</b> option.
+Running the module twice with and without <b>-n</b> flag and <b>weight</b>
+option may be useful in this specific case.
+
<p>The module recognizes two different formats of the flow direction map:
<div class="code"><pre>
Direction encoding for neighbors of x
-
+
135 90 45 3 2 1
180 x 360 4 x 8
225 270 315 5 6 7
-
+
degrees 45 degrees
CCW from East CCW from East
(r.watershed drainage)
More information about the grass-commit
mailing list