r3.flow: added option to sample a 3D raster by flowlines

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@61679 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster3d/r3.flow/flowline.c

@@ -37,19 +37,24 @@ static void write_segment(struct Map_info *flowline_vec,
 
 static void write_segment_db(struct field_info *finfo, dbDriver * driver,
 			     dbString * sql, const double velocity,
-			     double scalar_value, int write_scalar,
+			     double scalar_value, double sampled_map_value,
+			     int write_scalar, int use_sampled_map,
 			     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);
-
+    sprintf(buf, "insert into %s values (%d, %e", finfo->table, cat, velocity);
     db_set_string(sql, buf);
+    if (write_scalar) {
+	sprintf(buf, ", %e", scalar_value);
+	db_append_string(sql, buf);
+    }
+    if (use_sampled_map) {
+        sprintf(buf, ", %e", sampled_map_value);
+        db_append_string(sql, buf);
+    }
+    db_append_string(sql, ")");
+
 
     if (db_execute_immediate(driver, sql) != DB_OK) {
 	G_fatal_error(_("Unable to insert new record: '%s'"),
@@ -57,15 +62,14 @@ static void write_segment_db(struct field_info *finfo, dbDriver * driver,
     }
 }
 
-static double get_scalar_value(RASTER3D_Region * region,
-			       struct Gradient_info *gradient_info,
-			       double north, double east, double top)
+static double get_map_value(RASTER3D_Region * region, RASTER3D_Map *map,
+			    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,
+    Rast3d_get_value(map, col, row, depth, &val,
 		     DCELL_TYPE);
 
     return val;
@@ -86,7 +90,8 @@ static double get_scalar_value(RASTER3D_Region * region,
  */
 void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
 		      struct Gradient_info *gradient_info,
-		      RASTER3D_Map * flowacc, struct Integration *integration,
+		      RASTER3D_Map * flowacc, RASTER3D_Map * sampled_map,
+		      struct Integration *integration,
 		      struct Map_info *flowline_vec, struct line_cats *cats,
 		      struct line_pnts *points, int *cat, int if_table,
 		      struct field_info *finfo, dbDriver * driver)
@@ -100,8 +105,9 @@ void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
     int col, row, depth;
     int last_col, last_row, last_depth;
     int coor_diff;
+    DCELL scalar_value;
+    DCELL sampled_map_value;
     FCELL value;
-    double scalar_value;
     int *trav_coords;
     int size, trav_count;
     dbString sql;
@@ -114,7 +120,7 @@ void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
     last_col = last_row = last_depth = -1;
 
     size = 5;
-    scalar_value = 0;
+    value = 0;
     trav_coords = G_malloc(3 * size * sizeof(int));
 
     if (seed->flowline) {
@@ -152,11 +158,14 @@ void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
 	    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]);
+		    scalar_value = get_map_value(region, gradient_info->scalar_map,
+						 point[1], point[0], point[2]);
+		if (sampled_map)
+		    sampled_map_value = get_map_value(region, sampled_map,
+						      point[1], point[0], point[2]);
 		write_segment_db(finfo, driver, &sql, velocity, scalar_value,
-				 gradient_info->compute_gradient, *cat - 1);
+				 sampled_map_value, gradient_info->compute_gradient,
+				 sampled_map ? 1 : 0, *cat - 1);
 	    }
 	    else
 		Vect_append_point(points, point[0], point[1], point[2]);

raster3d/r3.flow/flowline.h

@@ -10,7 +10,7 @@ static const double VELOCITY_EPSILON = 1e-8;
 
 void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
 		      struct Gradient_info *gradient_info,
-		      RASTER3D_Map * flowacc,
+		      RASTER3D_Map * flowacc, RASTER3D_Map *sampled_map,
 		      struct Integration *integration,
 		      struct Map_info *flowline_vec, struct line_cats *cats,
 		      struct line_pnts *points, int *cat, int if_table,

raster3d/r3.flow/main.c

@@ -27,7 +27,7 @@
 
 static void create_table(struct Map_info *flowline_vec,
 			 struct field_info **f_info, dbDriver ** driver,
-			 int write_scalar)
+			 int write_scalar, int use_sampled_map)
 {
     dbString sql;
     char buf[200];
@@ -47,14 +47,14 @@ static void create_table(struct Map_info *flowline_vec,
 		      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);
+    sprintf(buf, "create table %s (cat integer, velocity double precision",
+	    fi->table);
     db_set_string(&sql, buf);
+    if (write_scalar)
+	db_append_string(&sql, ", input double precision");
+    if (use_sampled_map)
+	db_append_string(&sql, ", sampled double precision");
+    db_append_string(&sql, ")");
 
     db_begin_transaction(drvr);
     /* Create table */
@@ -147,13 +147,13 @@ static void init_flowaccum(RASTER3D_Region * region, RASTER3D_Map * flowacc)
 
 int main(int argc, char *argv[])
 {
-    struct Option *vector_opt, *seed_opt, *flowlines_opt, *flowacc_opt,
+    struct Option *vector_opt, *seed_opt, *flowlines_opt, *flowacc_opt, *sampled_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;
+    RASTER3D_Map *flowacc, *sampled;
     struct Integration integration;
     struct Seed seed;
     struct Gradient_info gradient_info;
@@ -215,6 +215,15 @@ int main(int argc, char *argv[])
 	_("Name for output flowaccumulation 3D raster");
     flowacc_opt->guisection = _("Output");
 
+    sampled_opt = G_define_standard_option(G_OPT_R3_INPUT);
+    sampled_opt->key = "sampled";
+    sampled_opt->required = NO;
+    sampled_opt->label =
+            _("Name for 3D raster sampled by flowlines");
+    sampled_opt->description =
+            _("Values of this 3D raster will be stored "
+              "as attributes of flowlines segments");
+
     unit_opt = G_define_option();
     unit_opt->key = "unit";
     unit_opt->type = TYPE_STRING;
@@ -287,6 +296,7 @@ int main(int argc, char *argv[])
     G_option_required(flowlines_opt, flowacc_opt, NULL);
     G_option_requires(seed_opt, flowlines_opt, NULL);
     G_option_requires(table_fl, flowlines_opt);
+    G_option_requires(sampled_opt, table_fl);
 
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
@@ -363,6 +373,19 @@ int main(int argc, char *argv[])
 	init_flowaccum(&region, flowacc);
     }
 
+    /* open 3D raster map used for sampling */
+    if (sampled_opt->answer) {
+	sampled = Rast3d_open_cell_old(sampled_opt->answer,
+				       G_find_raster3d(sampled_opt->answer, ""),
+				       &region, RASTER3D_TILE_SAME_AS_FILE,
+				       RASTER3D_USE_CACHE_DEFAULT);
+	if (!sampled)
+	    Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+			       sampled_opt->answer);
+    }
+    else
+	sampled = NULL;
+
     /* open new vector map of flowlines */
     if (flowlines_opt->answer) {
 	fl_cats = Vect_new_cats_struct();
@@ -375,7 +398,7 @@ int main(int argc, char *argv[])
 
 	if (if_table) {
 	    create_table(&fl_map, &finfo, &driver,
-			 gradient_info.compute_gradient);
+			 gradient_info.compute_gradient, sampled ? 1 : 0);
 	}
     }
 
@@ -429,14 +452,14 @@ int main(int argc, char *argv[])
 	    if (integration.direction_type == FLOWDIR_UP ||
 		integration.direction_type == FLOWDIR_BOTH) {
 		integration.actual_direction = FLOWDIR_UP;
-		compute_flowline(&region, &seed, &gradient_info, flowacc,
+		compute_flowline(&region, &seed, &gradient_info, flowacc, sampled,
 				 &integration, &fl_map, fl_cats, fl_points,
 				 &cat, if_table, finfo, driver);
 	    }
 	    if (integration.direction_type == FLOWDIR_DOWN ||
 		integration.direction_type == FLOWDIR_BOTH) {
 		integration.actual_direction = FLOWDIR_DOWN;
-		compute_flowline(&region, &seed, &gradient_info, flowacc,
+		compute_flowline(&region, &seed, &gradient_info, flowacc, sampled,
 				 &integration, &fl_map, fl_cats, fl_points,
 				 &cat, if_table, finfo, driver);
 	    }
@@ -474,7 +497,7 @@ int main(int argc, char *argv[])
 			    integration.direction_type == FLOWDIR_BOTH) {
 			    integration.actual_direction = FLOWDIR_UP;
 			    compute_flowline(&region, &seed, &gradient_info,
-					     flowacc, &integration, &fl_map,
+					     flowacc, sampled, &integration, &fl_map,
 					     fl_cats, fl_points, &cat,
 					     if_table, finfo, driver);
 			}
@@ -482,7 +505,7 @@ int main(int argc, char *argv[])
 			    integration.direction_type == FLOWDIR_BOTH) {
 			    integration.actual_direction = FLOWDIR_DOWN;
 			    compute_flowline(&region, &seed, &gradient_info,
-					     flowacc, &integration, &fl_map,
+					     flowacc, sampled, &integration, &fl_map,
 					     fl_cats, fl_points, &cat,
 					     if_table, finfo, driver);
 			}

raster3d/r3.flow/r3.flow.html

@@ -33,6 +33,8 @@ is represented by one vector line with one category. With <b>a</b> flag, an attr
 and each category (record) represents one segment of a flowline, so that attributes
 specific for segments can be written. In case of <b>vector_field</b> input, only velocity is written,
 in case of <b>input</b> option, also values of the input 3D raster are written.
+Option <b>sampled</b> allows to sample (query) given 3D raster by flow lines (computed from different 3D raster) and
+write the values of the given 3D raster as attributes of the flow line segments.
 Note that using <b>a</b> flag results in longer computation time, so consider increasing
 <b>step</b> and <b>max_error</b> parameter.
 
@@ -87,6 +89,7 @@ Again, we visualize the result in 3D view and we check 'use color for thematic r
 <h2>SEE ALSO</h2>
 <em>
 <a href="http://grass.osgeo.org/grass71/manuals/r.flow.html">r.flow</a>,
+<a href="http://grass.osgeo.org/grass71/manuals/r3.gradient.html">r3.gradient</a>,
 <a href="http://grass.osgeo.org/grass71/manuals/r3.gwflow.html">r3.gwflow</a>
 </em>