[GRASS-SVN] r72876 - in grass-addons/grass7/raster: . r.accumulate
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 21 22:20:15 PDT 2018
Author: hcho
Date: 2018-06-21 22:20:15 -0700 (Thu, 21 Jun 2018)
New Revision: 72876
Added:
grass-addons/grass7/raster/r.accumulate/
grass-addons/grass7/raster/r.accumulate/Makefile
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
grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_comparison.png
grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_example.png
grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_r_watershed_output.png
grass-addons/grass7/raster/r.accumulate/raster.c
Modified:
grass-addons/grass7/raster/Makefile
Log:
r.accumulate: New addon for calculating weighted flow accumulation using the flow direction map only
Modified: grass-addons/grass7/raster/Makefile
===================================================================
--- grass-addons/grass7/raster/Makefile 2018-06-22 03:18:43 UTC (rev 72875)
+++ grass-addons/grass7/raster/Makefile 2018-06-22 05:20:15 UTC (rev 72876)
@@ -14,6 +14,7 @@
# r.convert \
SUBDIRS = \
+ r.accumulate \
r.agent \
r.area \
r.basin \
Added: grass-addons/grass7/raster/r.accumulate/Makefile
===================================================================
--- grass-addons/grass7/raster/r.accumulate/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.accumulate/Makefile 2018-06-22 05:20:15 UTC (rev 72876)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../../../trunk
+
+PGM = r.accumulate
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/raster/r.accumulate/accumulate.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/accumulate.c (rev 0)
+++ grass-addons/grass7/raster/r.accumulate/accumulate.c 2018-06-22 05:20:15 UTC (rev 72876)
@@ -0,0 +1,102 @@
+#include <grass/raster.h>
+#include "global.h"
+
+double
+accumulate(CELL ** dir_buf, RASTER_MAP weight_buf, RASTER_MAP acc_buf,
+ char **done, int row, int col)
+{
+ int rows = weight_buf.rows, cols = weight_buf.cols;
+ int i, j, neighbor_dir, loop_dir, has_inflow;
+ double acc_flow;
+
+ if (done[row][col])
+ return get(acc_buf, row, col);
+
+ if (weight_buf.map.v)
+ acc_flow = get(weight_buf, row, col);
+ else
+ acc_flow = 1.0;
+
+ for (i = -1; i <= 1; i++) {
+ if (row + i < 0 || row + i >= rows)
+ continue;
+ for (j = -1; j <= 1; j++) {
+ if (col + j < 0 || col + j >= cols || (i == 0 && j == 0))
+ continue;
+ has_inflow = 0;
+ neighbor_dir = dir_buf[row + i][col + j];
+ loop_dir = 0;
+ switch (i) {
+ case -1:
+ switch (j) {
+ case -1:
+ if (neighbor_dir == 315) {
+ has_inflow = 1;
+ loop_dir = 135;
+ }
+ break;
+ case 0:
+ if (neighbor_dir == 270) {
+ has_inflow = 1;
+ loop_dir = 90;
+ }
+ break;
+ case 1:
+ if (neighbor_dir == 225) {
+ has_inflow = 1;
+ loop_dir = 45;
+ }
+ break;
+ }
+ break;
+ case 0:
+ switch (j) {
+ case -1:
+ if (neighbor_dir == 360) {
+ has_inflow = 1;
+ loop_dir = 180;
+ }
+ break;
+ case 1:
+ if (neighbor_dir == 180) {
+ has_inflow = 1;
+ loop_dir = 360;
+ }
+ break;
+ }
+ break;
+ case 1:
+ switch (j) {
+ case -1:
+ if (neighbor_dir == 45) {
+ has_inflow = 1;
+ loop_dir = 225;
+ }
+ break;
+ case 0:
+ if (neighbor_dir == 90) {
+ has_inflow = 1;
+ loop_dir = 270;
+ }
+ break;
+ case 1:
+ if (neighbor_dir == 135) {
+ has_inflow = 1;
+ loop_dir = 315;
+ }
+ break;
+ }
+ break;
+ }
+ if (has_inflow && dir_buf[row][col] != loop_dir)
+ acc_flow +=
+ accumulate(dir_buf, weight_buf, acc_buf, done, row + i,
+ col + j);
+ }
+ }
+
+ set(acc_buf, row, col, acc_flow);
+ done[row][col] = 1;
+
+ return acc_flow;
+}
Added: grass-addons/grass7/raster/r.accumulate/global.h
===================================================================
--- grass-addons/grass7/raster/r.accumulate/global.h (rev 0)
+++ grass-addons/grass7/raster/r.accumulate/global.h 2018-06-22 05:20:15 UTC (rev 72876)
@@ -0,0 +1,25 @@
+#include <grass/raster.h>
+
+#define DIR_UNKNOWN 0
+#define DIR_DEG 1
+#define DIR_DEG45 2
+
+typedef struct
+{
+ RASTER_MAP_TYPE type;
+ int rows, cols;
+ union
+ {
+ void **v;
+ CELL **c;
+ FCELL **f;
+ DCELL **d;
+ } map;
+} RASTER_MAP;
+
+/* raster.c */
+void set(RASTER_MAP, int, int, double);
+double get(RASTER_MAP, int, int);
+
+/* accumulate.c */
+double accumulate(CELL **, RASTER_MAP, RASTER_MAP, char **, int, int);
Added: grass-addons/grass7/raster/r.accumulate/main.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/main.c (rev 0)
+++ grass-addons/grass7/raster/r.accumulate/main.c 2018-06-22 05:20:15 UTC (rev 72876)
@@ -0,0 +1,203 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.accumulate
+ *
+ * AUTHOR(S): Huidae Cho <grass4u gmail.com>
+ *
+ * PURPOSE: Calculates weighted flow accumulation using a flow direction
+ * map.
+ *
+ * COPYRIGHT: (C) 2018 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+int main(int argc, char *argv[])
+{
+ struct GModule *module;
+ struct
+ {
+ struct Option *dir, *format, *weight, *acc;
+ } opt;
+ char *desc;
+ char *dir_name, *weight_name, *acc_name;
+ int dir_fd, weight_fd, acc_fd;
+ double dir_format;
+ struct Range dir_range;
+ CELL dir_min, dir_max;
+ char **done;
+ CELL **dir_buf;
+ RASTER_MAP weight_buf, acc_buf;
+ int rows, cols, row, col;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("hydrology"));
+ module->description =
+ _
+ ("Calculates weighted flow accumulation using a flow direction map.");
+
+ opt.dir = G_define_standard_option(G_OPT_R_INPUT);
+ opt.dir->label = _("Name of input direction map");
+
+ opt.format = G_define_option();
+ opt.format->type = TYPE_STRING;
+ opt.format->key = "format";
+ opt.format->label = _("Format of input direction map");
+ opt.format->required = YES;
+ opt.format->options = "auto,degree,45degree";
+ opt.format->answer = "auto";
+ 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)"));
+ opt.format->descriptions = desc;
+
+ opt.weight = G_define_standard_option(G_OPT_R_INPUT);
+ opt.weight->key = "weight";
+ opt.weight->required = NO;
+ opt.weight->label = _("Name of input flow weight map");
+
+ opt.acc = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt.acc->type = TYPE_STRING;
+ opt.acc->label = _("Name for output weighted flow accumulation map");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ dir_name = opt.dir->answer;
+ weight_name = opt.weight->answer;
+ acc_name = opt.acc->answer;
+
+ dir_fd = Rast_open_old(dir_name, "");
+ if (Rast_get_map_type(dir_fd) != CELL_TYPE)
+ G_fatal_error(_("Invalid directions map <%s>"), dir_name);
+
+ /* adapted from r.path */
+ if (Rast_read_range(dir_name, "", &dir_range) < 0)
+ G_fatal_error(_("Unable to read range file"));
+ Rast_get_range_min_max(&dir_range, &dir_min, &dir_max);
+ if (dir_max <= 0)
+ G_fatal_error(_("Invalid directions map <%s>"), dir_name);
+
+ dir_format = DIR_UNKNOWN;
+ if (strcmp(opt.format->answer, "degree") == 0) {
+ if (dir_max > 360)
+ G_fatal_error(_("Directional degrees can not be > 360"));
+ dir_format = DIR_DEG;
+ }
+ else if (strcmp(opt.format->answer, "45degree") == 0) {
+ if (dir_max > 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"));
+ }
+ else if (dir_max <= 360) {
+ dir_format = DIR_DEG;
+ 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>"),
+ dir_name);
+ }
+ if (dir_format == DIR_UNKNOWN)
+ G_fatal_error(_("Invalid directions format '%s'"),
+ 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);
+
+ rows = Rast_window_rows();
+ cols = Rast_window_cols();
+
+ done = G_malloc(rows * sizeof(char *));
+ dir_buf = G_malloc(rows * sizeof(CELL *));
+ for (row = 0; row < rows; row++) {
+ done[row] = G_calloc(cols, 1);
+ dir_buf[row] = Rast_allocate_c_buf();
+ Rast_get_c_row(dir_fd, dir_buf[row], row);
+ if (dir_format == DIR_DEG45) {
+ for (col = 0; col < cols; col++)
+ dir_buf[row][col] *= 45;
+ }
+ }
+
+ if (weight_fd >= 0) {
+ weight_buf.map.v = (void **)G_malloc(rows * sizeof(void *));
+ for (row = 0; row < rows; row++) {
+ weight_buf.map.v[row] =
+ (void *)Rast_allocate_buf(weight_buf.type);
+ Rast_get_row(weight_fd, weight_buf.map.v[row], row,
+ weight_buf.type);
+ }
+ }
+ else
+ weight_buf.map.v = NULL;
+ weight_buf.rows = rows;
+ weight_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;
+
+ for (row = 0; row < rows; row++) {
+ for (col = 0; col < cols; col++)
+ accumulate(dir_buf, weight_buf, acc_buf, done, row, col);
+ }
+
+ for (row = 0; row < rows; row++)
+ Rast_put_row(acc_fd, acc_buf.map.v[row], acc_buf.type);
+
+ G_free(done);
+ for (row = 0; row < rows; row++) {
+ G_free(dir_buf[row]);
+ if (weight_fd >= 0)
+ G_free(weight_buf.map.v[row]);
+ G_free(acc_buf.map.v[row]);
+ }
+ G_free(dir_buf);
+ 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);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/grass7/raster/r.accumulate/r.accumulate.html
===================================================================
--- grass-addons/grass7/raster/r.accumulate/r.accumulate.html (rev 0)
+++ grass-addons/grass7/raster/r.accumulate/r.accumulate.html 2018-06-22 05:20:15 UTC (rev 72876)
@@ -0,0 +1,76 @@
+<h2>DESCRIPTION</h2>
+<em>r.accumulate</em> calculates weighted flow accumulatin using a flow
+direction map.
+
+<h2>NOTES</h2>
+
+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.
+
+<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)
+</pre></div>
+
+<h2>EXAMPLES</h2>
+
+These examples use the North Carolina sample dataset.
+
+<p>Create the SFD flow accumulation and drainage maps using
+<em>r.watershed</em>:
+<div class="code"><pre>
+# set computational region
+g.region raster=elevation -p
+
+# calculate positive flow accumulation and drainage directions using r.watershed
+r.watershed elevation=elevation accumulation=flow_accum drainage=drain_directions -sa
+
+# calculate flow accumulation using r.accumulate
+r.accumulate input=drain_directions output=flow_accum_new
+
+# copy color table
+r.colors map=flow_accum_new raster=flow_accum
+
+# check difference betwen flow_accum and flow_accum_new
+r.mapcalc expression="accum_diff=if(flow_accum-flow_accum_new, flow_accum-flow_accum_new, null())"
+</pre></div>
+
+<img src="r_accumulate_nc_example.png">
+
+<p>For some reaon, there are slight differences between the two output maps.
+
+<p><img src="r_accumulate_nc_r_watershed_output.png">
+
+<p>The yellow and purple cells show the difference raster map
+(<i>accum_diff</i>). The red arrows and numbers represent drainage directions
+(<i>drain_directions</i>) and flow accumulation (<i>flow_accum</i> from
+<em>r.watershed</em>), respectively. Note that some cells close to headwater
+cells are assigned 1 even if they are located downstream of other cells.
+
+<p><img src="r_accumulate_nc_comparison.png">
+
+<p>For comparison, these numbers show the new flow accumulatoin
+(<i>flow_accum_new</i> from <em>r.accumulate</em>). The same cells are properly
+accumulated from the headwater cells.
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+<a href="mailto:grass4u at gmail com">Huidae Cho</a>
+
+<p><i>Last changed: $Date: 2018-06-15 21:55:50 -0400 (Fri, 15 Jun 2018) $</i>
Added: grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_comparison.png
===================================================================
(Binary files differ)
Index: grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_comparison.png
===================================================================
--- grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_comparison.png 2018-06-22 03:18:43 UTC (rev 72875)
+++ grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_comparison.png 2018-06-22 05:20:15 UTC (rev 72876)
Property changes on: grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_comparison.png
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+image/png
\ No newline at end of property
Added: grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_example.png
===================================================================
(Binary files differ)
Index: grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_example.png
===================================================================
--- grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_example.png 2018-06-22 03:18:43 UTC (rev 72875)
+++ grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_example.png 2018-06-22 05:20:15 UTC (rev 72876)
Property changes on: grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_example.png
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+image/png
\ No newline at end of property
Added: grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_r_watershed_output.png
===================================================================
(Binary files differ)
Index: grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_r_watershed_output.png
===================================================================
--- grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_r_watershed_output.png 2018-06-22 03:18:43 UTC (rev 72875)
+++ grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_r_watershed_output.png 2018-06-22 05:20:15 UTC (rev 72876)
Property changes on: grass-addons/grass7/raster/r.accumulate/r_accumulate_nc_r_watershed_output.png
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+image/png
\ No newline at end of property
Added: grass-addons/grass7/raster/r.accumulate/raster.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/raster.c (rev 0)
+++ grass-addons/grass7/raster/r.accumulate/raster.c 2018-06-22 05:20:15 UTC (rev 72876)
@@ -0,0 +1,38 @@
+#include <grass/raster.h>
+#include "global.h"
+
+void set(RASTER_MAP buf, int row, int col, double value)
+{
+ switch (buf.type) {
+ case CELL_TYPE:
+ buf.map.c[row][col] = (CELL) value;
+ break;
+ case FCELL_TYPE:
+ buf.map.f[row][col] = (FCELL) value;
+ break;
+ case DCELL_TYPE:
+ buf.map.d[row][col] = (DCELL) value;
+ break;
+ }
+
+ return;
+}
+
+double get(RASTER_MAP buf, int row, int col)
+{
+ double value = 0;
+
+ switch (buf.type) {
+ case CELL_TYPE:
+ value = (double)buf.map.c[row][col];
+ break;
+ case FCELL_TYPE:
+ value = (double)buf.map.f[row][col];
+ break;
+ case DCELL_TYPE:
+ value = buf.map.d[row][col];
+ break;
+ }
+
+ return value;
+}
More information about the grass-commit
mailing list