[GRASS-SVN] r72906 - grass-addons/grass7/raster/r.accumulate
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 25 20:08:04 PDT 2018
Author: hcho
Date: 2018-06-25 20:08:04 -0700 (Mon, 25 Jun 2018)
New Revision: 72906
Added:
grass-addons/grass7/raster/r.accumulate/calculate_lfp.c
grass-addons/grass7/raster/r.accumulate/line_list.c
Modified:
grass-addons/grass7/raster/r.accumulate/Makefile
grass-addons/grass7/raster/r.accumulate/delineate_streams.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: Implement longest flow path using flow directions (no dependency on r.stream.distance addon unlike r.lfp)
Modified: grass-addons/grass7/raster/r.accumulate/Makefile
===================================================================
--- grass-addons/grass7/raster/r.accumulate/Makefile 2018-06-25 21:28:53 UTC (rev 72905)
+++ grass-addons/grass7/raster/r.accumulate/Makefile 2018-06-26 03:08:04 UTC (rev 72906)
@@ -2,8 +2,8 @@
PGM = r.accumulate
-LIBES = $(RASTERLIB) $(VECTORLIB) $(GISLIB)
-DEPENDENCIES = $(RASTERDEP) $(VECTORDEP) $(GISDEP)
+LIBES = $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
EXTRA_INC = $(VECT_INC)
EXTRA_CFLAGS = $(VECT_CFLAGS)
Added: grass-addons/grass7/raster/r.accumulate/calculate_lfp.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/calculate_lfp.c (rev 0)
+++ grass-addons/grass7/raster/r.accumulate/calculate_lfp.c 2018-06-26 03:08:04 UTC (rev 72906)
@@ -0,0 +1,250 @@
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+struct neighbor_accum
+{
+ int row;
+ int col;
+ int accum;
+};
+
+static void add_table(struct Map_info *, char *, dbDriver **,
+ struct field_info **);
+static int trace_up(struct cell_map *, struct raster_map *,
+ struct Cell_head *, int, int, struct point_list *,
+ struct line_list *);
+static int compare_neighbor_accum(const void *, const void *);
+static int compare_line(const void *, const void *);
+
+void calculate_lfp(struct Map_info *Map, struct cell_map *dir_buf,
+ struct raster_map *accum_buf, int *id, char *id_colname,
+ struct point_list *outlet_pl)
+{
+ struct Cell_head window;
+ int rows = accum_buf->rows, cols = accum_buf->cols;
+ struct point_list pl;
+ struct line_list ll;
+ struct line_cats *Cats;
+ int i, cat;
+ dbDriver *driver = NULL;
+ struct field_info *Fi;
+ dbString sql;
+
+ if (id_colname) {
+ add_table(Map, id_colname, &driver, &Fi);
+ db_init_string(&sql);
+ }
+
+ G_get_set_window(&window);
+
+ init_point_list(&pl);
+ init_line_list(&ll);
+
+ Cats = Vect_new_cats_struct();
+
+ /* loop through all outlets and find the longest flow path for each */
+ cat = 1;
+ for (i = 0; i < outlet_pl->n; i++) {
+ int row = (int)Rast_northing_to_row(outlet_pl->y[i], &window);
+ int col = (int)Rast_easting_to_col(outlet_pl->x[i], &window);
+ int n;
+
+ /* if the outlet is outside the computational region, skip */
+ if (row < 0 || row >= rows || col < 0 || col >= cols) {
+ G_warning(_("Skip outlet (%f, %f) outside the current region"),
+ outlet_pl->x[i], outlet_pl->y[i]);
+ continue;
+ }
+
+ /* trace up flow accumulation */
+ reset_point_list(&pl);
+ reset_line_list(&ll);
+
+ trace_up(dir_buf, accum_buf, &window, row, col, &pl, &ll);
+
+ /* sort lines by length in descending order */
+ qsort(ll.lines, ll.n, sizeof(struct line *), compare_line);
+
+ if (!ll.n)
+ G_fatal_error(_("Failed to calculate the longest flow path for outlet (%f, %f)"),
+ outlet_pl->x[i], outlet_pl->y[i]);
+
+ /* write out the longest flow path */
+ for (n = 0; n < ll.n && ll.lines[n]->length == ll.lines[0]->length;
+ n++) {
+ Vect_reset_cats(Cats);
+
+ Vect_cat_set(Cats, 1, cat);
+ Vect_write_line(Map, GV_LINE, ll.lines[n]->Points, Cats);
+
+ if (driver) {
+ char *buf;
+
+ G_asprintf(&buf, "insert into %s (%s, %s) values (%d, %d)",
+ Fi->table, Fi->key, id_colname, cat, id[i]);
+ db_set_string(&sql, buf);
+
+ if (db_execute_immediate(driver, &sql) != DB_OK)
+ G_fatal_error(_("Unable to create table: %s"),
+ db_get_string(&sql));
+ db_free_string(&sql);
+ }
+
+ cat++;
+ }
+ }
+
+ free_point_list(&pl);
+ free_line_list(&ll);
+
+ Vect_destroy_cats_struct(Cats);
+
+ if (driver) {
+ db_commit_transaction(driver);
+ db_close_database_shutdown_driver(driver);
+ }
+}
+
+static void add_table(struct Map_info *Map, char *id_colname,
+ dbDriver ** pdriver, struct field_info **pFi)
+{
+ dbDriver *driver;
+ struct field_info *Fi;
+ char *buf;
+ dbString sql;
+
+ Fi = Vect_default_field_info(Map, 1, NULL, GV_1TABLE);
+
+ driver =
+ db_start_driver_open_database(Fi->driver,
+ Vect_subst_var(Fi->database, Map));
+ db_set_error_handler_driver(driver);
+ db_begin_transaction(driver);
+
+ if (!driver)
+ G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+ Fi->database, Fi->driver);
+
+ G_asprintf(&buf, "create table %s (%s integer, %s integer)", Fi->table,
+ Fi->key, id_colname);
+ db_init_string(&sql);
+ db_set_string(&sql, buf);
+
+ if (db_execute_immediate(driver, &sql) != DB_OK)
+ G_fatal_error(_("Unable to create table: %s"), db_get_string(&sql));
+ db_free_string(&sql);
+
+ if (db_grant_on_table
+ (driver, Fi->table, DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK)
+ G_fatal_error(_("Unable to grant privileges on table <%s>"),
+ Fi->table);
+
+ if (Vect_map_add_dblink
+ (Map, 1, NULL, Fi->table, GV_KEY_COLUMN, Fi->database, Fi->driver))
+ G_fatal_error(_("Unable to add database link for vector map <%s>"),
+ Vect_get_full_name(Map));
+
+ *pdriver = driver;
+ *pFi = Fi;
+}
+
+static int trace_up(struct cell_map *dir_buf, struct raster_map *accum_buf,
+ struct Cell_head *window, int row, int col,
+ struct point_list *pl, struct line_list *ll)
+{
+ int rows = dir_buf->rows, cols = dir_buf->cols;
+ double x, y;
+ int i, j, nup;
+ struct neighbor_accum up_accum[8];
+
+ /* if the current cell is outside the computational region, stop tracing */
+ if (row < 0 || row >= rows || col < 0 || col >= cols)
+ return 1;
+
+ /* add the current cell */
+ x = Rast_col_to_easting(col + 0.5, window);
+ y = Rast_row_to_northing(row + 0.5, window);
+ add_point(pl, x, y);
+
+ /* if the current accumulation is 1 (headwater), stop tracing */
+ if (get(accum_buf, row, col) == 1)
+ return 1;
+
+ nup = 0;
+ for (i = -1; i <= 1; i++) {
+ /* skip edge cells */
+ if (row + i < 0 || row + i >= rows)
+ continue;
+
+ for (j = -1; j <= 1; j++) {
+ /* skip the current and edge cells */
+ if ((i == 0 && j == 0) || col + j < 0 || col + j >= cols)
+ continue;
+
+ /* if a neighbor cell flows into the current cell, store its
+ * accumulation in the accum array */
+ if (dir_buf->c[row + i][col + j] == dir_checks[i + 1][j + 1][0]) {
+ up_accum[nup].row = row + i;
+ up_accum[nup].col = col + j;
+ up_accum[nup++].accum = get(accum_buf, row + i, col + j);
+ }
+ }
+ }
+
+ if (!nup)
+ G_fatal_error(_("No upstream cells found for a non-headwater cell at (%f, %f)"),
+ x, y);
+
+ /* sort upstream cells by accumulation in descending order */
+ qsort(up_accum, nup, sizeof(struct neighbor_accum),
+ compare_neighbor_accum);
+
+ /* trace up the largest accumulation */
+ for (i = 0; i < nup && up_accum[i].accum == up_accum[0].accum; i++) {
+ /* store pl->n to come back later */
+ int split_pl_n = pl->n;
+
+ /* if tracing is successful, store the line and its length */
+ if (trace_up
+ (dir_buf, accum_buf, window, up_accum[i].row, up_accum[i].col, pl,
+ ll) && pl->n > 0) {
+ struct line *line = (struct line *)G_malloc(sizeof(struct line));
+
+ line->Points = Vect_new_line_struct();
+ Vect_copy_xyz_to_pnts(line->Points, pl->x, pl->y, NULL, pl->n);
+ Vect_line_reverse(line->Points);
+ line->length = Vect_line_length(line->Points);
+ add_line(ll, line);
+ }
+
+ /* keep tracing if the next upstream cell has the same largest
+ * accumulation */
+ pl->n = split_pl_n;
+ }
+
+ return 0;
+}
+
+static int compare_neighbor_accum(const void *a, const void *b)
+{
+ const struct neighbor_accum *ap = (const struct neighbor_accum *)a;
+ const struct neighbor_accum *bp = (const struct neighbor_accum *)b;
+
+ /* descending by accum */
+ return bp->accum - ap->accum;
+}
+
+static int compare_line(const void *a, const void *b)
+{
+ const struct line **ap = (const struct line **)a;
+ const struct line **bp = (const struct line **)b;
+
+ /* descending by length */
+ if ((*ap)->length < (*bp)->length)
+ return 1;
+ if ((*ap)->length > (*bp)->length)
+ return -1;
+ return 0;
+}
Modified: grass-addons/grass7/raster/r.accumulate/delineate_streams.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/delineate_streams.c 2018-06-25 21:28:53 UTC (rev 72905)
+++ grass-addons/grass7/raster/r.accumulate/delineate_streams.c 2018-06-26 03:08:04 UTC (rev 72906)
@@ -11,16 +11,15 @@
struct Cell_head window;
int rows = accum_buf->rows, cols = accum_buf->cols;
int row, col;
- int i, j;
- char has_thresh_inflow;
struct point_list pl;
struct line_pnts *Points;
struct line_cats *Cats;
- int stream_id = 0;
+ int stream_cat = 0;
G_get_set_window(&window);
init_point_list(&pl);
+
Points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
@@ -27,6 +26,9 @@
/* loop through all cells to find headwater cells */
for (row = 0; row < rows; row++) {
for (col = 0; col < cols; col++) {
+ int i, j;
+ char has_thresh_inflow = 0;
+
/* if the current cell is less than the threshold, skip */
if (get(accum_buf, row, col) < thresh)
continue;
@@ -33,8 +35,6 @@
/* the current cell is greater than the threshold; check if it is
* headwater (no upstream cells greater than the threshold) */
- has_thresh_inflow = 0;
-
for (i = -1; i <= 1; i++) {
/* skip edge cells */
if (row + i < 0 || row + i >= rows)
@@ -66,7 +66,7 @@
Vect_reset_cats(Cats);
Vect_copy_xyz_to_pnts(Points, pl.x, pl.y, NULL, pl.n);
- Vect_cat_set(Cats, 1, ++stream_id);
+ Vect_cat_set(Cats, 1, ++stream_cat);
Vect_write_line(Map, GV_LINE, Points, Cats);
}
}
@@ -74,6 +74,7 @@
}
free_point_list(&pl);
+
Vect_destroy_line_struct(Points);
Vect_destroy_cats_struct(Cats);
}
Modified: grass-addons/grass7/raster/r.accumulate/global.h
===================================================================
--- grass-addons/grass7/raster/r.accumulate/global.h 2018-06-25 21:28:53 UTC (rev 72905)
+++ grass-addons/grass7/raster/r.accumulate/global.h 2018-06-26 03:08:04 UTC (rev 72906)
@@ -38,6 +38,18 @@
double *x, *y;
};
+struct line
+{
+ struct line_pnts *Points;
+ double length;
+};
+
+struct line_list
+{
+ int nalloc, n;
+ struct line **lines;
+};
+
#ifdef _MAIN_C_
#define GLOBAL
#else
@@ -64,6 +76,12 @@
void free_point_list(struct point_list *);
void add_point(struct point_list *, double, double);
+/* line_list.c */
+void init_line_list(struct line_list *);
+void reset_line_list(struct line_list *);
+void free_line_list(struct line_list *);
+void add_line(struct line_list *, struct line *);
+
/* accumulate.c */
void accumulate(struct cell_map *, struct raster_map *, struct raster_map *,
char **, char);
@@ -72,4 +90,8 @@
void delineate_streams(struct Map_info *, double, struct cell_map *,
struct raster_map *);
+/* calculate_lfp.c */
+void calculate_lfp(struct Map_info *, struct cell_map *, struct raster_map *,
+ int *, char *, struct point_list *);
+
#endif
Added: grass-addons/grass7/raster/r.accumulate/line_list.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/line_list.c (rev 0)
+++ grass-addons/grass7/raster/r.accumulate/line_list.c 2018-06-26 03:08:04 UTC (rev 72906)
@@ -0,0 +1,43 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+#define LINES_INCREMENT 1024
+
+void init_line_list(struct line_list *ll)
+{
+ ll->nalloc = ll->n = 0;
+ ll->lines = NULL;
+}
+
+void reset_line_list(struct line_list *ll)
+{
+ free_line_list(ll);
+}
+
+void free_line_list(struct line_list *ll)
+{
+ if (ll->lines) {
+ int i;
+
+ for (i = 0; i < ll->n; i++)
+ Vect_destroy_line_struct(ll->lines[i]->Points);
+ G_free(ll->lines);
+ }
+ init_line_list(ll);
+}
+
+/* adapted from r.path */
+void add_line(struct line_list *ll, struct line *l)
+{
+ if (ll->n == ll->nalloc) {
+ ll->nalloc += LINES_INCREMENT;
+ ll->lines =
+ (struct line **)G_realloc(ll->lines,
+ ll->nalloc * sizeof(struct line *));
+ if (!ll->lines)
+ G_fatal_error(_("Unable to increase line list"));
+ }
+ ll->lines[ll->n] = l;
+ ll->n++;
+}
Modified: grass-addons/grass7/raster/r.accumulate/main.c
===================================================================
--- grass-addons/grass7/raster/r.accumulate/main.c 2018-06-25 21:28:53 UTC (rev 72905)
+++ grass-addons/grass7/raster/r.accumulate/main.c 2018-06-26 03:08:04 UTC (rev 72906)
@@ -5,8 +5,8 @@
*
* AUTHOR(S): Huidae Cho <grass4u gmail.com>
*
- * PURPOSE: Calculates weighted flow accumulation and delineates stream
- * networks using a flow direction map.
+ * PURPOSE: Calculates weighted flow accumulation, stream networks, and
+ * the longest flow path using a flow direction map.
*
* COPYRIGHT: (C) 2018 by the GRASS Development Team
*
@@ -41,6 +41,10 @@
struct Option *accum;
struct Option *thresh;
struct Option *stream;
+ struct Option *lfp;
+ struct Option *coords;
+ struct Option *id_column;
+ struct Option *id;
} opt;
struct
{
@@ -47,7 +51,7 @@
struct Flag *neg;
} flag;
char *desc;
- char *dir_name, *weight_name, *accum_name, *stream_name;
+ char *dir_name, *weight_name, *accum_name, *stream_name, *lfp_name;
int dir_fd;
double dir_format, thresh;
struct Range dir_range;
@@ -57,8 +61,10 @@
struct cell_map dir_buf;
struct raster_map weight_buf, accum_buf;
int rows, cols, row, col;
- struct History hist;
struct Map_info Map;
+ struct point_list outlet_pl;
+ int *id;
+ char *id_colname;
G_gisinit(argv[0]);
@@ -66,7 +72,7 @@
G_add_keyword(_("raster"));
G_add_keyword(_("hydrology"));
module->description =
- _("Calculates weighted flow accumulation and delineates stream networks using a flow direction map.");
+ _("Calculates weighted flow accumulation, stream networks, and the longest flow path using a flow direction map.");
opt.dir = G_define_standard_option(G_OPT_R_INPUT);
opt.dir->key = "direction";
@@ -107,6 +113,29 @@
opt.stream->required = NO;
opt.stream->description = _("Name for output stream vector map");
+ opt.lfp = G_define_standard_option(G_OPT_V_OUTPUT);
+ opt.lfp->key = "longest_flow_path";
+ opt.lfp->required = NO;
+ opt.lfp->description = _("Name for output longest flow path vector map");
+
+ opt.coords = G_define_standard_option(G_OPT_M_COORDS);
+ opt.coords->multiple = YES;
+ opt.coords->description =
+ _("Coordinates of longest flow path outlet point");
+
+ opt.id_column = G_define_standard_option(G_OPT_DB_COLUMN);
+ opt.id_column->key = "id_column";
+ opt.id_column->description =
+ _("Name for output longest flow path ID column");
+
+ opt.id = G_define_option();
+ opt.id->key = "id";
+ opt.id->type = TYPE_INTEGER;
+ opt.id->multiple = YES;
+ opt.id->description = _("ID for longest flow path");
+
+
+
flag.neg = G_define_flag();
flag.neg->key = 'n';
flag.neg->label =
@@ -113,10 +142,14 @@
_("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);
- G_option_required(opt.accum, opt.stream, NULL);
+ * themselves can be negative; the longest flow path requires positive
+ * non-weighted accumulation */
+ G_option_exclusive(opt.weight, opt.lfp, flag.neg, NULL);
+ G_option_required(opt.accum, opt.stream, opt.lfp, NULL);
G_option_collective(opt.thresh, opt.stream, NULL);
+ G_option_collective(opt.lfp, opt.coords, NULL);
+ G_option_collective(opt.id_column, opt.id, NULL);
+ G_option_requires(opt.id, opt.coords, NULL);
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -125,6 +158,7 @@
weight_name = opt.weight->answer;
accum_name = opt.accum->answer;
stream_name = opt.stream->answer;
+ lfp_name = opt.lfp->answer;
dir_fd = Rast_open_old(dir_name, "");
if (Rast_get_map_type(dir_fd) != CELL_TYPE)
@@ -166,6 +200,35 @@
opt.format->answer);
/* end of r.path */
+ /* read id and outlet */
+ id = NULL;
+ init_point_list(&outlet_pl);
+ if (opt.coords->answers) {
+ int i = 0;
+
+ while (opt.coords->answers[i]) {
+ add_point(&outlet_pl, atof(opt.coords->answers[i]),
+ atof(opt.coords->answers[i + 1]));
+ i += 2;
+ }
+ if (opt.id->answers) {
+ i = 0;
+ while (opt.id->answers[i])
+ i++;
+ if (i < outlet_pl.n)
+ G_fatal_error(_("Too few longest flow path IDs specified"));
+ if (i > outlet_pl.n)
+ G_fatal_error(_("Too many longest flow path IDs specified"));
+ id = (int *)G_malloc(outlet_pl.n * sizeof(int));
+ i = 0;
+ while (opt.id->answers[i]) {
+ id[i] = atoi(opt.id->answers[i]);
+ i++;
+ }
+ }
+ }
+ id_colname = opt.id_column->answer;
+
thresh = opt.thresh->answer ? atof(opt.thresh->answer) : 0.0;
neg = flag.neg->answer;
@@ -222,6 +285,7 @@
/* write out buffer to the accumulatoin map if requested */
if (accum_name) {
int accum_fd = Rast_open_new(accum_name, accum_buf.type);
+ struct History hist;
for (row = 0; row < rows; row++)
Rast_put_row(accum_fd, accum_buf.map.v[row], accum_buf.type);
@@ -263,6 +327,21 @@
Vect_close(&Map);
}
+ /* calculate the longest flow path */
+ if (lfp_name) {
+ if (Vect_open_new(&Map, lfp_name, 0) < 0)
+ G_fatal_error(_("Unable to create vector map <%s>"), lfp_name);
+ Vect_set_map_name(&Map, "Longest flow path");
+ Vect_hist_command(&Map);
+
+ calculate_lfp(&Map, &dir_buf, &accum_buf, id, id_colname, &outlet_pl);
+
+ if (!Vect_build(&Map))
+ G_warning(_("Unable to build topology for vector map <%s>"),
+ lfp_name);
+ Vect_close(&Map);
+ }
+
/* free buffer memory */
for (row = 0; row < rows; row++) {
G_free(dir_buf.c[row]);
Modified: grass-addons/grass7/raster/r.accumulate/r.accumulate.html
===================================================================
--- grass-addons/grass7/raster/r.accumulate/r.accumulate.html 2018-06-25 21:28:53 UTC (rev 72905)
+++ grass-addons/grass7/raster/r.accumulate/r.accumulate.html 2018-06-26 03:08:04 UTC (rev 72906)
@@ -1,7 +1,7 @@
<h2>DESCRIPTION</h2>
-<em>r.accumulate</em> calculates weighted flow accumulatin and delineates
-stream networks using a flow direction map.
+<em>r.accumulate</em> calculates weighted flow accumulatin, stream networks,
+and the longest flow path using a flow direction map.
<h2>NOTES</h2>
More information about the grass-commit
mailing list