[GRASS-SVN] r72898 - grass-addons/grass7/raster/r.accumulate
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 25 02:19:07 PDT 2018
Author: hcho
Date: 2018-06-25 02:19:07 -0700 (Mon, 25 Jun 2018)
New Revision: 72898
Modified:
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/raster.c
Log:
r.accumulate: Implement stream delineation (stream option); Rename input/output options to direction/accumulation (multiple output options)
Modified: grass-addons/grass7/raster/r.accumulate/Makefile
===================================================================
--- grass-addons/grass7/raster/r.accumulate/Makefile 2018-06-24 17:17:46 UTC (rev 72897)
+++ grass-addons/grass7/raster/r.accumulate/Makefile 2018-06-25 09:19:07 UTC (rev 72898)
@@ -2,8 +2,8 @@
PGM = r.accumulate
-LIBES = $(RASTERLIB) $(GISLIB)
-DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+LIBES = $(GISLIB) $(RASTERLIB) $(VECTORLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP) $(VECTORDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Modified: grass-addons/grass7/raster/r.accumulate/accumulate.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/accumulate.c 2018-06-24 17:17:46 UTC (rev 72897)
+++ grass-addons/grass7/raster/r.accumulate/accumulate.c 2018-06-25 09:19:07 UTC (rev 72898)
@@ -1,44 +1,30 @@
-#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, char neg, int row, int col)
+double accumulate(struct cell_map *dir_buf, struct raster_map *weight_buf,
+ struct raster_map *accum_buf, char **done, char neg,
+ int row, int col)
{
- static int dir_checks[3][3][2] = {
- {{SE, NW}, {S, N}, {SW, NE}},
- {{E, W}, {0, 0}, {W, E}},
- {{NE, SW}, {N, S}, {NW, SE}}
- };
- int rows = acc_buf.rows, cols = acc_buf.cols;
+ int rows = dir_buf->rows, cols = dir_buf->cols;
int i, j;
char incomplete = 0;
- double acc;
+ double accum;
/* 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);
+ accum = get(accum_buf, row, col);
/* for negative accumulation, always return its absolute value;
* otherwise return it as is */
- return neg && acc < 0 ? -acc : acc;
+ return neg && accum < 0 ? -accum : accum;
}
/* 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);
+ if (weight_buf->map.v)
+ accum = get(weight_buf, row, col);
else
- acc = 1.0;
+ accum = 1.0;
/* loop through all neighbor cells and see if any of them drains into the
* current cell (are there upstream cells?) */
@@ -65,13 +51,13 @@
/* 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]) {
+ if (dir_buf->c[row + i][col + j] == dir_checks[i + 1][j + 1][0] &&
+ dir_buf->c[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,
+ * positive value, so accum is always positive (cell count);
+ * otherwise, accum is weighted accumulation */
+ accum +=
+ accumulate(dir_buf, weight_buf, accum_buf, done, neg,
row + i, col + j);
/* if the neighbor cell is incomplete, the current cell also
@@ -85,12 +71,12 @@
/* 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);
+ set(accum_buf, row, col, neg && incomplete ? -accum : accum);
/* 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;
+ /* for negative accumulation, accum is already positive */
+ return accum;
}
Modified: grass-addons/grass7/raster/r.accumulate/global.h
===================================================================
--- grass-addons/grass7/raster/r.accumulate/global.h 2018-06-24 17:17:46 UTC (rev 72897)
+++ grass-addons/grass7/raster/r.accumulate/global.h 2018-06-25 09:19:07 UTC (rev 72898)
@@ -1,7 +1,26 @@
+#ifndef _GLOBAL_H_
+#define _GLOBAL_H_
+
#include <grass/raster.h>
+#include <grass/vector.h>
-typedef struct
+#define NE 1
+#define N 2
+#define NW 3
+#define W 4
+#define SW 5
+#define S 6
+#define SE 7
+#define E 8
+
+struct cell_map
{
+ int rows, cols;
+ CELL **c;
+};
+
+struct raster_map
+{
RASTER_MAP_TYPE type;
int rows, cols;
union
@@ -11,11 +30,46 @@
FCELL **f;
DCELL **d;
} map;
-} RASTER_MAP;
+};
+struct point_list
+{
+ int nalloc, n;
+ double *x, *y;
+};
+
+#ifdef _MAIN_C_
+#define GLOBAL
+#else
+#define GLOBAL extern
+#endif
+
+GLOBAL int dir_checks[3][3][2]
+#ifdef _MAIN_C_
+ = {
+ {{SE, NW}, {S, N}, {SW, NE}},
+ {{E, W}, {0, 0}, {W, E}},
+ {{NE, SW}, {N, S}, {NW, SE}}
+}
+#endif
+;
+
/* raster.c */
-void set(RASTER_MAP, int, int, double);
-double get(RASTER_MAP, int, int);
+void set(struct raster_map *, int, int, double);
+double get(struct raster_map *, int, int);
+/* point_list.c */
+void init_point_list(struct point_list *);
+void reset_point_list(struct point_list *);
+void free_point_list(struct point_list *);
+void add_point(struct point_list *, double, double);
+
/* accumulate.c */
-double accumulate(CELL **, RASTER_MAP, RASTER_MAP, char **, char, int, int);
+double accumulate(struct cell_map *, struct raster_map *, struct raster_map *,
+ char **, char, int, int);
+
+/* delineate_streams.c */
+void delineate_streams(struct Map_info *, double, struct cell_map *,
+ struct raster_map *);
+
+#endif
Modified: grass-addons/grass7/raster/r.accumulate/main.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/main.c 2018-06-24 17:17:46 UTC (rev 72897)
+++ grass-addons/grass7/raster/r.accumulate/main.c 2018-06-25 09:19:07 UTC (rev 72898)
@@ -5,8 +5,8 @@
*
* AUTHOR(S): Huidae Cho <grass4u gmail.com>
*
- * PURPOSE: Calculates weighted flow accumulation using a flow direction
- * map.
+ * PURPOSE: Calculates weighted flow accumulation and delineates stream
+ * networks using a flow direction map.
*
* COPYRIGHT: (C) 2018 by the GRASS Development Team
*
@@ -16,10 +16,13 @@
*
*****************************************************************************/
+#define _MAIN_C_
+
#include <stdlib.h>
#include <string.h>
#include <grass/gis.h>
#include <grass/raster.h>
+#include <grass/vector.h>
#include <grass/glocale.h>
#include "global.h"
@@ -32,7 +35,12 @@
struct GModule *module;
struct
{
- struct Option *dir, *format, *weight, *acc;
+ struct Option *dir;
+ struct Option *format;
+ struct Option *weight;
+ struct Option *accum;
+ struct Option *thresh;
+ struct Option *stream;
} opt;
struct
{
@@ -39,17 +47,18 @@
struct Flag *neg;
} flag;
char *desc;
- char *dir_name, *weight_name, *acc_name;
- int dir_fd, weight_fd, acc_fd;
- double dir_format;
+ char *dir_name, *weight_name, *accum_name, *stream_name;
+ int dir_fd, weight_fd, accum_fd;
+ double dir_format, thresh;
struct Range dir_range;
CELL dir_min, dir_max;
char neg;
char **done;
- CELL **dir_buf;
- RASTER_MAP weight_buf, acc_buf;
+ struct cell_map dir_buf;
+ struct raster_map weight_buf, accum_buf;
int rows, cols, row, col;
struct History hist;
+ struct Map_info Map;
G_gisinit(argv[0]);
@@ -57,9 +66,10 @@
G_add_keyword(_("raster"));
G_add_keyword(_("hydrology"));
module->description =
- _("Calculates weighted flow accumulation using a flow direction map.");
+ _("Calculates weighted flow accumulation and delineates stream networks using a flow direction map.");
opt.dir = G_define_standard_option(G_OPT_R_INPUT);
+ opt.dir->key = "direction";
opt.dir->description = _("Name of input direction map");
opt.format = G_define_option();
@@ -80,11 +90,23 @@
opt.weight->required = NO;
opt.weight->description = _("Name of input flow weight map");
- opt.acc = G_define_standard_option(G_OPT_R_OUTPUT);
- opt.acc->type = TYPE_STRING;
- opt.acc->description =
+ opt.accum = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt.accum->key = "accumulation";
+ opt.accum->required = NO;
+ opt.accum->type = TYPE_STRING;
+ opt.accum->description =
_("Name for output weighted flow accumulation map");
+ opt.thresh = G_define_option();
+ opt.thresh->key = "threshold";
+ opt.thresh->type = TYPE_DOUBLE;
+ opt.thresh->label = _("Minimum flow accumulation for streams");
+
+ opt.stream = G_define_standard_option(G_OPT_V_OUTPUT);
+ opt.stream->key = "stream";
+ opt.stream->required = NO;
+ opt.stream->description = _("Name for output stream vector map");
+
flag.neg = G_define_flag();
flag.neg->key = 'n';
flag.neg->label =
@@ -93,6 +115,8 @@
/* weighting doesn't support negative accumulation because weights
* themselves can be negative */
G_option_exclusive(opt.weight, flag.neg, NULL);
+ G_option_required(opt.accum, opt.stream, NULL);
+ G_option_collective(opt.thresh, opt.stream, NULL);
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -99,7 +123,8 @@
dir_name = opt.dir->answer;
weight_name = opt.weight->answer;
- acc_name = opt.acc->answer;
+ accum_name = opt.accum->answer;
+ stream_name = opt.stream->answer;
dir_fd = Rast_open_old(dir_name, "");
if (Rast_get_map_type(dir_fd) != CELL_TYPE)
@@ -141,6 +166,7 @@
opt.format->answer);
/* end of r.path */
+ thresh = opt.thresh->answer ? atof(opt.thresh->answer) : 0.0;
neg = flag.neg->answer;
rows = Rast_window_rows();
@@ -148,21 +174,24 @@
/* initialize the done array and read the direction map */
done = G_malloc(rows * sizeof(char *));
- dir_buf = G_malloc(rows * sizeof(CELL *));
+ dir_buf.rows = rows;
+ dir_buf.cols = cols;
+ dir_buf.c = 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) {
+ dir_buf.c[row] = Rast_allocate_c_buf();
+ Rast_get_c_row(dir_fd, dir_buf.c[row], row);
+ if (dir_format == DIR_DEG) {
for (col = 0; col < cols; col++)
- dir_buf[row][col] *= 45;
+ dir_buf.c[row][col] /= 45.0;
}
}
+ Rast_close(dir_fd);
/* 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);
+ accum_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 *));
@@ -172,61 +201,80 @@
Rast_get_row(weight_fd, weight_buf.map.v[row], row,
weight_buf.type);
}
+ Rast_close(weight_fd);
}
else {
weight_fd = -1;
weight_buf.map.v = NULL;
- acc_buf.type = CELL_TYPE;
+ accum_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 *));
+ /* create accumulation buffer */
+ accum_buf.rows = rows;
+ accum_buf.cols = cols;
+ accum_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);
+ accum_buf.map.v[row] = (void *)Rast_allocate_buf(accum_buf.type);
- /* create a new output map */
- acc_fd = Rast_open_new(acc_name, acc_buf.type);
+ /* create a new accumulation map if requested */
+ accum_fd = accum_name ? Rast_open_new(accum_name, accum_buf.type) : -1;
/* 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);
+ accumulate(&dir_buf, &weight_buf, &accum_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);
+ /* write out buffer to the accumulatoin map */
+ if (accum_fd >= 0) {
+ for (row = 0; row < rows; row++)
+ Rast_put_row(accum_fd, accum_buf.map.v[row], accum_buf.type);
+ Rast_close(accum_fd);
- /* close all maps */
- Rast_close(dir_fd);
- Rast_close(acc_fd);
- if (weight_fd >= 0)
- Rast_close(weight_fd);
+ /* write history */
+ Rast_put_cell_title(accum_name,
+ weight_name ? "Weighted flow accumulation" :
+ (neg ?
+ "Flow accumulation with likely underestimates" :
+ "Flow accumulation"));
+ Rast_short_history(accum_name, "raster", &hist);
+ Rast_command_history(&hist);
+ Rast_write_history(accum_name, &hist);
+ }
/* 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(done);
- G_free(dir_buf);
- G_free(acc_buf.map.v);
if (weight_fd >= 0)
G_free(weight_buf.map.v);
- /* write history */
- 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);
+ /* delineate stream networks */
+ if (stream_name) {
+ if (Vect_open_new(&Map, stream_name, 0) < 0)
+ G_fatal_error(_("Unable to create vector map <%s>"), stream_name);
+ Vect_set_map_name(&Map, "Stream network");
+ Vect_hist_command(&Map);
+ delineate_streams(&Map, thresh, &dir_buf, &accum_buf);
+
+ if (!Vect_build(&Map))
+ G_warning(_("Unable to build topology for vector map <%s>"),
+ stream_name);
+ Vect_close(&Map);
+ }
+
+ /* free buffer memory */
+ for (row = 0; row < rows; row++) {
+ G_free(dir_buf.c[row]);
+ G_free(accum_buf.map.v[row]);
+ }
+ G_free(dir_buf.c);
+ G_free(accum_buf.map.v);
+
exit(EXIT_SUCCESS);
}
Modified: grass-addons/grass7/raster/r.accumulate/r.accumulate.html
===================================================================
--- grass-addons/grass7/raster/r.accumulate/r.accumulate.html 2018-06-24 17:17:46 UTC (rev 72897)
+++ grass-addons/grass7/raster/r.accumulate/r.accumulate.html 2018-06-25 09:19:07 UTC (rev 72898)
@@ -1,7 +1,7 @@
<h2>DESCRIPTION</h2>
-<em>r.accumulate</em> calculates weighted flow accumulatin using a flow
-direction map.
+<em>r.accumulate</em> calculates weighted flow accumulatin and delineates
+stream networks using a flow direction map.
<h2>NOTES</h2>
@@ -36,17 +36,19 @@
These examples use the North Carolina sample dataset.
-<p>Create the SFD flow accumulation and drainage maps using
-<em>r.watershed</em>:
+<h3>Flow accumulation</h3>
+
+<p>Calculate flow accumulation using <em>r.watershed</em> and
+<em>r.accumulate</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
+r.watershed elevation=elevation accumulation=flow_accum drainage=drain_directions -s -a
# calculate flow accumulation using r.accumulate
-r.accumulate input=drain_directions output=flow_accum_new
+r.accumulate direction=drain_directions accumulation=flow_accum_new
# copy color table
r.colors map=flow_accum_new raster=flow_accum
@@ -73,6 +75,36 @@
(<i>flow_accum_new</i> from <em>r.accumulate</em>). The same cells are properly
accumulated from the headwater cells.
+<h3>Stream network delineation</h3>
+
+<p>Calculate flow accumulation and delineate stream networks at once:
+<div class="code"><pre>
+# use r.accumulate to create flow_accum_new and streams_new at once
+r.accumulate direction=drain_directions accumulation=flow_accum_new threshold=50000 stream=streams_new
+
+# or delineate stream networks only without creating an accumulation map
+r.accumulate direction=drain_directions threshold=50000 stream=streams_new_only
+
+# use r.stream.extract, elevation, and flow_accum_new to delineate stream networks
+r.stream.extract elevation=elevation accumulation=flow_accum_new threshold=50000 stream_vector=streams direction=drain_directions_2
+</pre></div>
+
+<img src="r_accumulate_nc_stream_example.png">
+
+<p>In this example, <em>r.accumulate</em> and <em>r.stream.extract</em>
+produced slightly different stream networks where the flow turns 90 degrees or
+there are multiple downstream cells with the same elevation. The following
+figure shows an example where the <i>streams</i> (red from
+<em>r.stream.extract</em>) and <i>streams_new</i> (blue from
+<em>r.accumulate</em>) vector maps present different stream paths. The numbers
+show cell elevations (<i>elevation</i> map), and the yellow
+(<i>drain_directions</i> map) and green (<i>drain_directions_2</i> map) arrows
+represent flow directions generated by <em>r.watershed</em> with <b>-s</b> flag
+and <em>r.stream.extract</em>, respectively. Note that the two flow direction
+maps are different.
+
+<p><img src="r_accumulate_nc_stream_comparison.png">
+
<h2>SEE ALSO</h2>
<em>
Modified: grass-addons/grass7/raster/r.accumulate/raster.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/raster.c 2018-06-24 17:17:46 UTC (rev 72897)
+++ grass-addons/grass7/raster/r.accumulate/raster.c 2018-06-25 09:19:07 UTC (rev 72898)
@@ -1,17 +1,17 @@
#include <grass/raster.h>
#include "global.h"
-void set(RASTER_MAP buf, int row, int col, double value)
+void set(struct raster_map *buf, int row, int col, double value)
{
- switch (buf.type) {
+ switch (buf->type) {
case CELL_TYPE:
- buf.map.c[row][col] = (CELL) value;
+ buf->map.c[row][col] = (CELL) value;
break;
case FCELL_TYPE:
- buf.map.f[row][col] = (FCELL) value;
+ buf->map.f[row][col] = (FCELL) value;
break;
case DCELL_TYPE:
- buf.map.d[row][col] = (DCELL) value;
+ buf->map.d[row][col] = (DCELL) value;
break;
}
@@ -18,19 +18,19 @@
return;
}
-double get(RASTER_MAP buf, int row, int col)
+double get(struct raster_map *buf, int row, int col)
{
double value = 0;
- switch (buf.type) {
+ switch (buf->type) {
case CELL_TYPE:
- value = (double)buf.map.c[row][col];
+ value = (double)buf->map.c[row][col];
break;
case FCELL_TYPE:
- value = (double)buf.map.f[row][col];
+ value = (double)buf->map.f[row][col];
break;
case DCELL_TYPE:
- value = buf.map.d[row][col];
+ value = buf->map.d[row][col];
break;
}
More information about the grass-commit
mailing list