[GRASS-SVN] r61453 - grass-addons/grass7/raster3d/r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jul 29 16:39:47 PDT 2014
Author: annakrat
Date: 2014-07-29 16:39:47 -0700 (Tue, 29 Jul 2014)
New Revision: 61453
Modified:
grass-addons/grass7/raster3d/r3.flow/Makefile
grass-addons/grass7/raster3d/r3.flow/flowline.c
grass-addons/grass7/raster3d/r3.flow/flowline.h
grass-addons/grass7/raster3d/r3.flow/integrate.c
grass-addons/grass7/raster3d/r3.flow/integrate.h
grass-addons/grass7/raster3d/r3.flow/main.c
grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h
Log:
r3.flow: velocity and input 3d raster value added as a record in attribute table for each flowline segment (-a flag)
Modified: grass-addons/grass7/raster3d/r3.flow/Makefile
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/Makefile 2014-07-29 21:25:32 UTC (rev 61452)
+++ grass-addons/grass7/raster3d/r3.flow/Makefile 2014-07-29 23:39:47 UTC (rev 61453)
@@ -2,8 +2,8 @@
PGM=r3.flow
-LIBES = $(RASTER3DLIB) $(RASTERLIB) $(GISLIB) $(VECTORLIB)
-DEPENDENCIES = $(RASTER3DDEP) $(RASTERDEP) $(GISDEP) $(VECTORDEP)
+LIBES = $(RASTER3DLIB) $(RASTERLIB) $(GISLIB) $(VECTORLIB) $(DBMILIB)
+DEPENDENCIES = $(RASTER3DDEP) $(RASTERDEP) $(GISDEP) $(VECTORDEP) $(DBMIDEP)
EXTRA_INC = $(VECT_INC)
EXTRA_CFLAGS = $(VECT_CFLAGS)
Modified: grass-addons/grass7/raster3d/r3.flow/flowline.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/flowline.c 2014-07-29 21:25:32 UTC (rev 61452)
+++ grass-addons/grass7/raster3d/r3.flow/flowline.c 2014-07-29 23:39:47 UTC (rev 61453)
@@ -14,12 +14,63 @@
#include <grass/raster3d.h>
#include <grass/vector.h>
#include <grass/glocale.h>
+#include <grass/dbmi.h>
#include "r3flow_structs.h"
#include "integrate.h"
#include "flowline.h"
#include "voxel_traversal.h"
+static void write_segment(struct Map_info *flowline_vec,
+ struct line_pnts *points, struct line_cats *cats,
+ const double *point, int *cat)
+{
+ Vect_append_point(points, point[0], point[1], point[2]);
+
+ Vect_cat_set(cats, 1, *cat);
+ (*cat)++;
+ Vect_write_line(flowline_vec, GV_LINE, points, cats);
+ Vect_reset_line(points);
+ Vect_reset_cats(cats);
+ Vect_append_point(points, point[0], point[1], point[2]);
+}
+
+static void write_segment_db(struct field_info *finfo, dbDriver * driver,
+ dbString * sql, const double velocity,
+ double scalar_value, int write_scalar,
+ const int cat)
+{
+ char buf[200];
+
+ if (write_scalar)
+ sprintf(buf, "insert into %s values ( %d, %e, %e )", finfo->table,
+ cat, scalar_value, velocity);
+ else
+ sprintf(buf, "insert into %s values ( %d, %e)", finfo->table, cat,
+ velocity);
+
+ db_set_string(sql, buf);
+
+ if (db_execute_immediate(driver, sql) != DB_OK) {
+ G_fatal_error(_("Unable to insert new record: '%s'"),
+ db_get_string(sql));
+ }
+}
+
+static double get_scalar_value(RASTER3D_Region * region,
+ struct Gradient_info *gradient_info,
+ double north, double east, double top)
+{
+ int col, row, depth;
+ double val;
+
+ Rast3d_location2coord(region, north, east, top, &col, &row, &depth);
+ Rast3d_get_value(gradient_info->scalar_map, col, row, depth, &val,
+ DCELL_TYPE);
+
+ return val;
+}
+
/*!
\brief Computes flowline by integrating velocity field.
@@ -30,13 +81,15 @@
\param flowline_vec pointer to Map_info struct of flowline vector
\param cats pointer to line_cats struct of flowline vector
\param points pointer to line_pnts struct of flowline vector
- \param cat category of the newly created flow line
+ \param[in,out] cat starting category of the newly created flow line
+ \param if_table TRUE if attribute table should be created and filled
*/
void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
struct Gradient_info *gradient_info,
RASTER3D_Map * flowacc, struct Integration *integration,
struct Map_info *flowline_vec, struct line_cats *cats,
- struct line_pnts *points, const int cat)
+ struct line_pnts *points, int *cat, int if_table,
+ struct field_info *finfo, dbDriver * driver)
{
int i, j, count;
double delta_t;
@@ -48,8 +101,11 @@
int last_col, last_row, last_depth;
int coor_diff;
FCELL value;
+ double scalar_value;
int *trav_coords;
int size, trav_count;
+ dbString sql;
+ double velocity;
point[0] = seed->x;
point[1] = seed->y;
@@ -58,11 +114,13 @@
last_col = last_row = last_depth = -1;
size = 5;
+ scalar_value = 0;
trav_coords = G_malloc(3 * size * sizeof(int));
if (seed->flowline) {
/* append first point */
Vect_append_point(points, seed->x, seed->y, seed->z);
+ db_init_string(&sql);
}
count = 1;
while (count <= integration->limit) {
@@ -79,19 +137,30 @@
velocity_norm, integration->cell_size);
/* integrate */
- min_step = get_time_step("cell", MIN_STEP, velocity_norm,
+ min_step = get_time_step("cell", integration->min_step, velocity_norm,
integration->cell_size);
- max_step = get_time_step("cell", MAX_STEP, velocity_norm,
+ max_step = get_time_step("cell", integration->max_step, velocity_norm,
integration->cell_size);
delta_t *= (integration->actual_direction == FLOWDIR_UP ? 1 : -1);
if (rk45_integrate_next
(region, gradient_info, point, new_point,
- &delta_t, min_step, max_step) < 0)
+ &delta_t, &velocity, min_step, max_step,
+ integration->max_error) < 0)
break;
- if (seed->flowline)
- Vect_append_point(points, new_point[0], new_point[1],
- new_point[2]);
+ if (seed->flowline) {
+ if (if_table) {
+ write_segment(flowline_vec, points, cats, new_point, cat);
+ if (gradient_info->compute_gradient)
+ scalar_value = get_scalar_value(region, gradient_info,
+ point[1], point[0],
+ point[2]);
+ write_segment_db(finfo, driver, &sql, velocity, scalar_value,
+ gradient_info->compute_gradient, *cat - 1);
+ }
+ else
+ Vect_append_point(points, point[0], point[1], point[2]);
+ }
if (seed->flowaccum) {
Rast3d_location2coord(region, new_point[1], new_point[0],
new_point[2], &col, &row, &depth);
@@ -132,11 +201,13 @@
}
if (seed->flowline) {
if (points->n_points > 1) {
- Vect_cat_set(cats, 1, cat);
+ Vect_cat_set(cats, 1, *cat);
+ (*cat)++;
Vect_write_line(flowline_vec, GV_LINE, points, cats);
G_debug(1, "Flowline ended after %d steps", count - 1);
}
Vect_reset_line(points);
Vect_reset_cats(cats);
+ db_free_string(&sql);
}
}
Modified: grass-addons/grass7/raster3d/r3.flow/flowline.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/flowline.h 2014-07-29 21:25:32 UTC (rev 61452)
+++ grass-addons/grass7/raster3d/r3.flow/flowline.h 2014-07-29 23:39:47 UTC (rev 61453)
@@ -13,5 +13,6 @@
RASTER3D_Map * flowacc,
struct Integration *integration,
struct Map_info *flowline_vec, struct line_cats *cats,
- struct line_pnts *points, const int cat);
+ struct line_pnts *points, int *cat, int if_table,
+ struct field_info *finfo, dbDriver *driver);
#endif // FLOWLINE_H
Modified: grass-addons/grass7/raster3d/r3.flow/integrate.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/integrate.c 2014-07-29 21:25:32 UTC (rev 61452)
+++ grass-addons/grass7/raster3d/r3.flow/integrate.c 2014-07-29 23:39:47 UTC (rev 61453)
@@ -83,13 +83,14 @@
*/
static int rk45_next(RASTER3D_Region * region, struct Gradient_info *gradient_info,
const double *point, double *next_point,
- const double delta_t, double *error)
+ const double delta_t, double *velocity, double *error)
{
double tmp1[6][3]; /* 3 is 3 dimensions, 6 is the number of k's */
double tmp_point[3];
double vel_x, vel_y, vel_z;
double sum_tmp;
int i, j, k;
+ double vel_sq;
if (get_velocity(region, gradient_info, point[0], point[1], point[2],
&vel_x, &vel_y, &vel_z) < 0)
@@ -118,6 +119,7 @@
tmp1[i][2] = vel_z;
}
+ vel_sq = 0;
/* compute next point */
for (j = 0; j < 3; j++) {
sum_tmp = 0;
@@ -125,7 +127,10 @@
sum_tmp += C[i] * tmp1[i][j];
}
next_point[j] = point[j] + delta_t * sum_tmp;
+ vel_sq += sum_tmp * sum_tmp;
}
+ *velocity = sqrt(vel_sq);
+
if (!Rast3d_is_valid_location
(region, next_point[1], next_point[0], next_point[2]))
return -1;
@@ -158,8 +163,9 @@
*/
int rk45_integrate_next(RASTER3D_Region * region,
struct Gradient_info *gradient_info, const double *point,
- double *next_point, double *delta_t,
- const double min_step, const double max_step)
+ double *next_point, double *delta_t, double *velocity,
+ const double min_step, const double max_step,
+ const double max_error)
{
double estimated_error;
double error_ratio;
@@ -176,16 +182,17 @@
*delta_t = max_step * (*delta_t > 0 ? 1 : -1);
/* try to iteratively decrease error to less than max error */
- while (estimated_error > MAX_ERROR) {
+ while (estimated_error > max_error) {
/* compute next point and get estimated error */
if (rk45_next
- (region, gradient_info, point, next_point, *delta_t, error) == 0)
+ (region, gradient_info, point, next_point, *delta_t,
+ velocity, error) == 0)
estimated_error = norm(error[0], error[1], error[2]);
else
return -1;
/* compute new step size (empirically) */
- error_ratio = estimated_error / MAX_ERROR;
+ error_ratio = estimated_error / max_error;
if (error_ratio == 0.0)
tmp = *delta_t > 0 ? min_step : -min_step;
else if (error_ratio > 1)
@@ -211,7 +218,7 @@
/* break when the adjustment was needed (not sure why) */
if (do_break) {
if (rk45_next(region, gradient_info, point, next_point,
- *delta_t, error) < 0)
+ *delta_t, velocity, error) < 0)
return -1;
break;
}
Modified: grass-addons/grass7/raster3d/r3.flow/integrate.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/integrate.h 2014-07-29 21:25:32 UTC (rev 61452)
+++ grass-addons/grass7/raster3d/r3.flow/integrate.h 2014-07-29 23:39:47 UTC (rev 61453)
@@ -5,10 +5,6 @@
#include "r3flow_structs.h"
-static const double MAX_ERROR = 1.0e-6;
-static const double MIN_STEP = 0.01;
-static const double MAX_STEP = 1;
-
/* Cash-Karp parameters */
static const double B[5][5] = { {1. / 5, 0, 0, 0, 0},
{3. / 40, 9. / 40, 0, 0, 0},
@@ -33,6 +29,7 @@
int rk45_integrate_next(RASTER3D_Region * region,
struct Gradient_info *gradient_info, const double *point,
double *next_point, double *delta_t,
- const double min_step, const double max_step);
+ double *velocity, const double min_step,
+ const double max_step, const double max_error);
#endif // INTEGRATE_H
Modified: grass-addons/grass7/raster3d/r3.flow/main.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/main.c 2014-07-29 21:25:32 UTC (rev 61452)
+++ grass-addons/grass7/raster3d/r3.flow/main.c 2014-07-29 23:39:47 UTC (rev 61453)
@@ -19,12 +19,59 @@
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/vector.h>
+#include <grass/dbmi.h>
#include <grass/glocale.h>
#include "r3flow_structs.h"
#include "flowline.h"
+static void create_table(struct Map_info *flowline_vec,
+ struct field_info **f_info, dbDriver ** driver,
+ int write_scalar)
+{
+ dbString sql;
+ char buf[200];
+ dbDriver *drvr;
+ struct field_info *fi;
+ db_init_string(&sql);
+ fi = Vect_default_field_info(flowline_vec, 1, NULL, GV_1TABLE);
+ *f_info = fi;
+ Vect_map_add_dblink(flowline_vec, 1, NULL, fi->table, GV_KEY_COLUMN,
+ fi->database, fi->driver);
+ drvr = db_start_driver_open_database(fi->driver,
+ Vect_subst_var(fi->database,
+ flowline_vec));
+ if (drvr == NULL) {
+ G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+ Vect_subst_var(fi->database, flowline_vec), fi->driver);
+ }
+ *driver = drvr;
+ if (write_scalar)
+ sprintf(buf, "create table %s (cat integer, input double precision, "
+ "velocity double precision)", fi->table);
+ else
+ sprintf(buf,
+ "create table %s (cat integer, velocity double precision)",
+ fi->table);
+ db_set_string(&sql, buf);
+
+ db_begin_transaction(drvr);
+ /* Create table */
+ if (db_execute_immediate(drvr, &sql) != DB_OK) {
+ G_fatal_error(_("Unable to create table: %s"), db_get_string(&sql));
+ }
+ if (db_create_index2(drvr, fi->table, fi->key) != DB_OK)
+ G_warning(_("Unable to create index for table <%s>, key <%s>"),
+ fi->table, fi->key);
+ /* Grant */
+ if (db_grant_on_table
+ (drvr, fi->table, DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK) {
+ G_fatal_error(_("Unable to grant privileges on table <%s>"),
+ fi->table);
+ }
+}
+
static void check_vector_input_maps(struct Option *vector_opt,
struct Option *seed_opt)
{
@@ -101,7 +148,9 @@
int main(int argc, char *argv[])
{
struct Option *vector_opt, *seed_opt, *flowlines_opt, *flowacc_opt,
- *scalar_opt, *unit_opt, *step_opt, *limit_opt, *skip_opt, *dir_opt;
+ *scalar_opt, *unit_opt, *step_opt, *limit_opt, *skip_opt, *dir_opt,
+ *error_opt;
+ struct Flag *table_fl;
struct GModule *module;
RASTER3D_Region region;
RASTER3D_Map *flowacc;
@@ -114,7 +163,10 @@
struct Map_info fl_map;
struct line_cats *fl_cats; /* for flowlines */
struct line_pnts *fl_points; /* for flowlines */
+ struct field_info *finfo;
+ dbDriver *driver;
int cat; /* cat of flowlines */
+ int if_table;
int i, r, c, d;
char *desc;
int n_seeds, seed_count, ltype;
@@ -125,7 +177,8 @@
G_add_keyword(_("raster3d"));
G_add_keyword(_("voxel"));
G_add_keyword(_("flowline"));
- module->description = _("Computes 3D flow lines and 3D flow accumulation.");
+ module->description =
+ _("Computes 3D flow lines and 3D flow accumulation.");
scalar_opt = G_define_standard_option(G_OPT_R3_INPUT);
@@ -197,6 +250,16 @@
limit_opt->description = _("Maximum number of steps");
limit_opt->guisection = _("Integration");
+ error_opt = G_define_option();
+ error_opt->key = "max_error";
+ error_opt->type = TYPE_DOUBLE;
+ error_opt->required = NO;
+ error_opt->answer = "1e-5";
+ error_opt->label = _("Maximum error of integration");
+ error_opt->description = _("Influences step, increase maximum error "
+ "to allow bigger steps");
+ error_opt->guisection = _("Integration");
+
skip_opt = G_define_option();
skip_opt->key = "skip";
skip_opt->type = TYPE_INTEGER;
@@ -213,16 +276,26 @@
dir_opt->options = "up,down,both";
dir_opt->answer = "down";
dir_opt->description = _("Compute flowlines upstream, "
- "downstream or in both direction.");
+ "downstream or in both direction.");
+ table_fl = G_define_flag();
+ table_fl->key = 'a';
+ table_fl->description = _("Create and fill attribute table");
+
G_option_required(scalar_opt, vector_opt, NULL);
G_option_exclusive(scalar_opt, vector_opt, NULL);
G_option_required(flowlines_opt, flowacc_opt, NULL);
G_option_requires(seed_opt, flowlines_opt, NULL);
+ G_option_requires(table_fl, flowlines_opt);
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
+ driver = NULL;
+ finfo = NULL;
+
+ if_table = table_fl->answer ? TRUE : FALSE;
+
check_vector_input_maps(vector_opt, seed_opt);
Rast3d_init_defaults();
@@ -237,13 +310,16 @@
integration.unit = "cell";
integration.step = 0.25;
}
+ integration.max_error = atof(error_opt->answer);
+ integration.max_step = 5 * integration.step;
+ integration.min_step = integration.step / 5;
integration.limit = atof(limit_opt->answer);
if (strcmp(dir_opt->answer, "up") == 0)
- integration.direction_type = FLOWDIR_UP;
+ integration.direction_type = FLOWDIR_UP;
else if (strcmp(dir_opt->answer, "down") == 0)
- integration.direction_type = FLOWDIR_DOWN;
+ integration.direction_type = FLOWDIR_DOWN;
else
- integration.direction_type = FLOWDIR_BOTH;
+ integration.direction_type = FLOWDIR_BOTH;
/* cell size is the diagonal */
@@ -296,6 +372,11 @@
flowlines_opt->answer);
Vect_hist_command(&fl_map);
+
+ if (if_table) {
+ create_table(&fl_map, &finfo, &driver,
+ gradient_info.compute_gradient);
+ }
}
n_seeds = 0;
@@ -321,6 +402,7 @@
G_debug(1, "Number of seeds is %d", n_seeds);
seed_count = 0;
+ cat = 1;
if (seed_opt->answer) {
seed_points = Vect_new_line_struct();
@@ -344,18 +426,19 @@
seed.flowaccum = FALSE;
}
G_percent(seed_count, n_seeds, 1);
- cat = seed_count + 1;
- if (integration.direction_type == FLOWDIR_UP ||
+ if (integration.direction_type == FLOWDIR_UP ||
integration.direction_type == FLOWDIR_BOTH) {
integration.actual_direction = FLOWDIR_UP;
compute_flowline(®ion, &seed, &gradient_info, flowacc,
- &integration, &fl_map, fl_cats, fl_points, cat);
+ &integration, &fl_map, fl_cats, fl_points,
+ &cat, if_table, finfo, driver);
}
- if (integration.direction_type == FLOWDIR_DOWN ||
+ if (integration.direction_type == FLOWDIR_DOWN ||
integration.direction_type == FLOWDIR_BOTH) {
integration.actual_direction = FLOWDIR_DOWN;
compute_flowline(®ion, &seed, &gradient_info, flowacc,
- &integration, &fl_map, fl_cats, fl_points, cat);
+ &integration, &fl_map, fl_cats, fl_points,
+ &cat, if_table, finfo, driver);
}
seed_count++;
}
@@ -386,20 +469,22 @@
if (seed.flowaccum || seed.flowline) {
G_percent(seed_count, n_seeds, 1);
- cat = seed_count + 1;
- if (integration.direction_type == FLOWDIR_UP ||
+
+ if (integration.direction_type == FLOWDIR_UP ||
integration.direction_type == FLOWDIR_BOTH) {
integration.actual_direction = FLOWDIR_UP;
- compute_flowline(®ion, &seed, &gradient_info,
- flowacc, &integration, &fl_map,
- fl_cats, fl_points, cat);
+ compute_flowline(®ion, &seed, &gradient_info,
+ flowacc, &integration, &fl_map,
+ fl_cats, fl_points, &cat,
+ if_table, finfo, driver);
}
- if (integration.direction_type == FLOWDIR_DOWN ||
+ if (integration.direction_type == FLOWDIR_DOWN ||
integration.direction_type == FLOWDIR_BOTH) {
integration.actual_direction = FLOWDIR_DOWN;
compute_flowline(®ion, &seed, &gradient_info,
flowacc, &integration, &fl_map,
- fl_cats, fl_points, cat);
+ fl_cats, fl_points, &cat,
+ if_table, finfo, driver);
}
seed_count++;
}
@@ -407,7 +492,12 @@
}
}
}
+ G_percent(1, 1, 1);
if (flowlines_opt->answer) {
+ if (if_table) {
+ db_commit_transaction(driver);
+ db_close_database_shutdown_driver(driver);
+ }
Vect_destroy_line_struct(fl_points);
Vect_destroy_cats_struct(fl_cats);
Vect_build(&fl_map);
Modified: grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h 2014-07-29 21:25:32 UTC (rev 61452)
+++ grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h 2014-07-29 23:39:47 UTC (rev 61453)
@@ -22,6 +22,9 @@
double step;
double cell_size;
int limit;
+ double max_error;
+ double max_step;
+ double min_step;
};
struct Array
More information about the grass-commit
mailing list