1
2 /****************************************************************************
3 *
4 * MODULE: r3.flow
5 * AUTHOR(S): Anna Petrasova kratochanna <at> gmail <dot> com
6 * PURPOSE: Computes 3D flow lines and flow accumulation based on 3D raster map(s)
7 * COPYRIGHT: (C) 2014 by the GRASS Development Team
8 *
9 * This program is free software under the GNU General Public
10 * License (>=v2). Read the file COPYING that comes with GRASS
11 * for details.
12 *
13 *****************************************************************************/
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <string.h>
17 #include <math.h>
18 #include <grass/raster3d.h>
19 #include <grass/gis.h>
20 #include <grass/raster.h>
21 #include <grass/vector.h>
22 #include <grass/dbmi.h>
23 #include <grass/glocale.h>
24
25 #include "r3flow_structs.h"
26 #include "flowline.h"
27
create_table(struct Map_info * flowline_vec,struct field_info ** f_info,dbDriver ** driver,int write_scalar,int use_sampled_map)28 static void create_table(struct Map_info *flowline_vec,
29 struct field_info **f_info, dbDriver ** driver,
30 int write_scalar, int use_sampled_map)
31 {
32 dbString sql;
33 char buf[200];
34 dbDriver *drvr;
35 struct field_info *fi;
36
37 db_init_string(&sql);
38 fi = Vect_default_field_info(flowline_vec, 1, NULL, GV_1TABLE);
39 *f_info = fi;
40 Vect_map_add_dblink(flowline_vec, 1, NULL, fi->table, GV_KEY_COLUMN,
41 fi->database, fi->driver);
42 drvr = db_start_driver_open_database(fi->driver,
43 Vect_subst_var(fi->database,
44 flowline_vec));
45 if (drvr == NULL) {
46 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
47 Vect_subst_var(fi->database, flowline_vec), fi->driver);
48 }
49 db_set_error_handler_driver(drvr);
50
51 *driver = drvr;
52 sprintf(buf, "create table %s (cat integer, velocity double precision",
53 fi->table);
54 db_set_string(&sql, buf);
55 if (write_scalar)
56 db_append_string(&sql, ", input double precision");
57 if (use_sampled_map)
58 db_append_string(&sql, ", sampled double precision");
59 db_append_string(&sql, ")");
60
61 db_begin_transaction(drvr);
62 /* Create table */
63 if (db_execute_immediate(drvr, &sql) != DB_OK) {
64 G_fatal_error(_("Unable to create table: %s"), db_get_string(&sql));
65 }
66 if (db_create_index2(drvr, fi->table, fi->key) != DB_OK)
67 G_warning(_("Unable to create index for table <%s>, key <%s>"),
68 fi->table, fi->key);
69 /* Grant */
70 if (db_grant_on_table
71 (drvr, fi->table, DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK) {
72 G_fatal_error(_("Unable to grant privileges on table <%s>"),
73 fi->table);
74 }
75 }
76
check_vector_input_maps(struct Option * vector_opt,struct Option * seed_opt)77 static void check_vector_input_maps(struct Option *vector_opt,
78 struct Option *seed_opt)
79 {
80 int i;
81
82 /* Check for velocity components maps. */
83 if (vector_opt->answers != NULL) {
84 for (i = 0; i < 3; i++) {
85 if (vector_opt->answers[i] != NULL) {
86 if (NULL == G_find_raster3d(vector_opt->answers[i], ""))
87 Rast3d_fatal_error(_("3D raster map <%s> not found"),
88 vector_opt->answers[i]);
89 }
90 else {
91 Rast3d_fatal_error(_("Please provide three 3D raster maps"));
92 }
93 }
94 }
95
96 if (seed_opt->answer != NULL) {
97 if (NULL == G_find_vector2(seed_opt->answer, ""))
98 G_fatal_error(_("Vector seed map <%s> not found"),
99 seed_opt->answer);
100 }
101
102 }
103
load_input_raster3d_maps(struct Option * scalar_opt,struct Option * vector_opt,struct Gradient_info * gradient_info,RASTER3D_Region * region)104 static void load_input_raster3d_maps(struct Option *scalar_opt,
105 struct Option *vector_opt,
106 struct Gradient_info *gradient_info,
107 RASTER3D_Region * region)
108 {
109 int i;
110
111 if (scalar_opt->answer) {
112 gradient_info->scalar_map =
113 Rast3d_open_cell_old(scalar_opt->answer,
114 G_find_raster3d(scalar_opt->answer, ""),
115 region, RASTER3D_TILE_SAME_AS_FILE,
116 RASTER3D_USE_CACHE_DEFAULT);
117 if (!gradient_info->scalar_map)
118 Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
119 scalar_opt->answer);
120 gradient_info->compute_gradient = TRUE;
121 }
122 else {
123 for (i = 0; i < 3; i++) {
124 gradient_info->velocity_maps[i] =
125 Rast3d_open_cell_old(vector_opt->answers[i],
126 G_find_raster3d(vector_opt->answers[i],
127 ""), region,
128 RASTER3D_TILE_SAME_AS_FILE,
129 RASTER3D_USE_CACHE_DEFAULT);
130
131 if (!gradient_info->velocity_maps[i])
132 Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
133 vector_opt->answers[i]);
134 }
135 gradient_info->compute_gradient = FALSE;
136 }
137 }
138
init_flowaccum(RASTER3D_Region * region,RASTER3D_Map * flowacc)139 static void init_flowaccum(RASTER3D_Region * region, RASTER3D_Map * flowacc)
140 {
141 int c, r, d;
142
143 for (d = 0; d < region->depths; d++)
144 for (r = 0; r < region->rows; r++)
145 for (c = 0; c < region->cols; c++)
146 if (Rast3d_put_float(flowacc, c, r, d, 0) != 1)
147 Rast3d_fatal_error(_("init_flowaccum: error in Rast3d_put_float"));
148 }
149
main(int argc,char * argv[])150 int main(int argc, char *argv[])
151 {
152 struct Option *vector_opt, *seed_opt, *flowlines_opt, *flowacc_opt, *sampled_opt,
153 *scalar_opt, *unit_opt, *step_opt, *limit_opt, *skip_opt, *dir_opt,
154 *error_opt;
155 struct Flag *table_fl;
156 struct GModule *module;
157 RASTER3D_Region region;
158 RASTER3D_Map *flowacc, *sampled;
159 struct Integration integration;
160 struct Seed seed;
161 struct Gradient_info gradient_info;
162 struct Map_info seed_Map;
163 struct line_pnts *seed_points;
164 struct line_cats *seed_cats;
165 struct Map_info fl_map;
166 struct line_cats *fl_cats; /* for flowlines */
167 struct line_pnts *fl_points; /* for flowlines */
168 struct field_info *finfo;
169 dbDriver *driver;
170 int cat; /* cat of flowlines */
171 int if_table;
172 int i, r, c, d;
173 char *desc;
174 int n_seeds, seed_count, ltype;
175 int skip[3];
176
177 G_gisinit(argv[0]);
178 module = G_define_module();
179 G_add_keyword(_("raster3d"));
180 G_add_keyword(_("hydrology"));
181 G_add_keyword(_("voxel"));
182 module->description =
183 _("Computes 3D flow lines and 3D flow accumulation.");
184
185
186 scalar_opt = G_define_standard_option(G_OPT_R3_INPUT);
187 scalar_opt->required = NO;
188 scalar_opt->guisection = _("Input");
189
190 vector_opt = G_define_standard_option(G_OPT_R3_INPUTS);
191 vector_opt->key = "vector_field";
192 vector_opt->required = NO;
193 vector_opt->description = _("Names of three 3D raster maps describing "
194 "x, y, z components of vector field");
195 vector_opt->guisection = _("Input");
196
197 seed_opt = G_define_standard_option(G_OPT_V_INPUT);
198 seed_opt->required = NO;
199 seed_opt->key = "seed_points";
200 seed_opt->description = _("If no map is provided, "
201 "flow lines are generated "
202 "from each cell of the input 3D raster");
203 seed_opt->label = _("Name of vector map with points "
204 "from which flow lines are generated");
205 seed_opt->guisection = _("Input");
206
207 flowlines_opt = G_define_standard_option(G_OPT_V_OUTPUT);
208 flowlines_opt->key = "flowline";
209 flowlines_opt->required = NO;
210 flowlines_opt->description = _("Name for vector map of flow lines");
211 flowlines_opt->guisection = _("Output");
212
213 flowacc_opt = G_define_standard_option(G_OPT_R3_OUTPUT);
214 flowacc_opt->key = "flowaccumulation";
215 flowacc_opt->required = NO;
216 flowacc_opt->description =
217 _("Name for output flowaccumulation 3D raster");
218 flowacc_opt->guisection = _("Output");
219
220 sampled_opt = G_define_standard_option(G_OPT_R3_INPUT);
221 sampled_opt->key = "sampled";
222 sampled_opt->required = NO;
223 sampled_opt->label =
224 _("Name for 3D raster sampled by flowlines");
225 sampled_opt->description =
226 _("Values of this 3D raster will be stored "
227 "as attributes of flowlines segments");
228
229 unit_opt = G_define_option();
230 unit_opt->key = "unit";
231 unit_opt->type = TYPE_STRING;
232 unit_opt->required = NO;
233 unit_opt->answer = "cell";
234 unit_opt->options = "time,length,cell";
235 desc = NULL;
236 G_asprintf(&desc,
237 "time;%s;"
238 "length;%s;"
239 "cell;%s",
240 _("elapsed time"),
241 _("length in map units"), _("length in cells (voxels)"));
242 unit_opt->descriptions = desc;
243 unit_opt->label = _("Unit of integration step");
244 unit_opt->description = _("Default unit is cell");
245 unit_opt->guisection = _("Integration");
246
247 step_opt = G_define_option();
248 step_opt->key = "step";
249 step_opt->type = TYPE_DOUBLE;
250 step_opt->required = NO;
251 step_opt->answer = "0.25";
252 step_opt->label = _("Integration step in selected unit");
253 step_opt->description = _("Default step is 0.25 cell");
254 step_opt->guisection = _("Integration");
255
256 limit_opt = G_define_option();
257 limit_opt->key = "limit";
258 limit_opt->type = TYPE_INTEGER;
259 limit_opt->required = NO;
260 limit_opt->answer = "2000";
261 limit_opt->description = _("Maximum number of steps");
262 limit_opt->guisection = _("Integration");
263
264 error_opt = G_define_option();
265 error_opt->key = "max_error";
266 error_opt->type = TYPE_DOUBLE;
267 error_opt->required = NO;
268 error_opt->answer = "1e-5";
269 error_opt->label = _("Maximum error of integration");
270 error_opt->description = _("Influences step, increase maximum error "
271 "to allow bigger steps");
272 error_opt->guisection = _("Integration");
273
274 skip_opt = G_define_option();
275 skip_opt->key = "skip";
276 skip_opt->type = TYPE_INTEGER;
277 skip_opt->required = NO;
278 skip_opt->multiple = YES;
279 skip_opt->description =
280 _("Number of cells between flow lines in x, y and z direction");
281
282 dir_opt = G_define_option();
283 dir_opt->key = "direction";
284 dir_opt->type = TYPE_STRING;
285 dir_opt->required = NO;
286 dir_opt->multiple = NO;
287 dir_opt->options = "up,down,both";
288 dir_opt->answer = "down";
289 dir_opt->description = _("Compute flowlines upstream, "
290 "downstream or in both direction.");
291
292 table_fl = G_define_flag();
293 table_fl->key = 'a';
294 table_fl->description = _("Create and fill attribute table");
295
296 G_option_required(scalar_opt, vector_opt, NULL);
297 G_option_exclusive(scalar_opt, vector_opt, NULL);
298 G_option_required(flowlines_opt, flowacc_opt, NULL);
299 G_option_requires(seed_opt, flowlines_opt, NULL);
300 G_option_requires(table_fl, flowlines_opt, NULL);
301 G_option_requires(sampled_opt, table_fl, NULL);
302
303 if (G_parser(argc, argv))
304 exit(EXIT_FAILURE);
305
306 driver = NULL;
307 finfo = NULL;
308
309 if_table = table_fl->answer ? TRUE : FALSE;
310
311 check_vector_input_maps(vector_opt, seed_opt);
312
313 Rast3d_init_defaults();
314 Rast3d_get_window(®ion);
315
316 /* set up integration variables */
317 if (step_opt->answer) {
318 integration.step = atof(step_opt->answer);
319 integration.unit = unit_opt->answer;
320 }
321 else {
322 integration.unit = "cell";
323 integration.step = 0.25;
324 }
325 integration.max_error = atof(error_opt->answer);
326 integration.max_step = 5 * integration.step;
327 integration.min_step = integration.step / 5;
328 integration.limit = atof(limit_opt->answer);
329 if (strcmp(dir_opt->answer, "up") == 0)
330 integration.direction_type = FLOWDIR_UP;
331 else if (strcmp(dir_opt->answer, "down") == 0)
332 integration.direction_type = FLOWDIR_DOWN;
333 else
334 integration.direction_type = FLOWDIR_BOTH;
335
336
337 /* cell size is the diagonal */
338 integration.cell_size = sqrt(region.ns_res * region.ns_res +
339 region.ew_res * region.ew_res +
340 region.tb_res * region.tb_res);
341
342 /* set default skip if needed */
343 if (skip_opt->answers) {
344 for (i = 0; i < 3; i++) {
345 if (skip_opt->answers[i] != NULL) {
346 skip[i] = atoi(skip_opt->answers[i]);
347 }
348 else {
349 G_fatal_error(_("Please provide 3 integer values for skip option."));
350 }
351 }
352 }
353 else {
354 skip[0] = fmax(1, region.cols / 10);
355 skip[1] = fmax(1, region.rows / 10);
356 skip[2] = fmax(1, region.depths / 10);
357
358 }
359
360 /* open raster 3D maps of velocity components */
361 gradient_info.initialized = FALSE;
362 load_input_raster3d_maps(scalar_opt, vector_opt, &gradient_info, ®ion);
363
364
365 /* open new 3D raster map of flowacumulation */
366 if (flowacc_opt->answer) {
367 flowacc = Rast3d_open_new_opt_tile_size(flowacc_opt->answer,
368 RASTER3D_USE_CACHE_DEFAULT,
369 ®ion, FCELL_TYPE, 32);
370
371
372 if (!flowacc)
373 Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
374 flowacc_opt->answer);
375 init_flowaccum(®ion, flowacc);
376 }
377
378 /* open 3D raster map used for sampling */
379 if (sampled_opt->answer) {
380 sampled = Rast3d_open_cell_old(sampled_opt->answer,
381 G_find_raster3d(sampled_opt->answer, ""),
382 ®ion, RASTER3D_TILE_SAME_AS_FILE,
383 RASTER3D_USE_CACHE_DEFAULT);
384 if (!sampled)
385 Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
386 sampled_opt->answer);
387 }
388 else
389 sampled = NULL;
390
391 /* open new vector map of flowlines */
392 if (flowlines_opt->answer) {
393 fl_cats = Vect_new_cats_struct();
394 fl_points = Vect_new_line_struct();
395 if (Vect_open_new(&fl_map, flowlines_opt->answer, TRUE) < 0)
396 G_fatal_error(_("Unable to create vector map <%s>"),
397 flowlines_opt->answer);
398
399 Vect_hist_command(&fl_map);
400
401 if (if_table) {
402 create_table(&fl_map, &finfo, &driver,
403 gradient_info.compute_gradient, sampled ? 1 : 0);
404 }
405 }
406
407 n_seeds = 0;
408 /* open vector map of seeds */
409 if (seed_opt->answer) {
410 if (Vect_open_old2(&seed_Map, seed_opt->answer, "", "1") < 0)
411 G_fatal_error(_("Unable to open vector map <%s>"),
412 seed_opt->answer);
413 if (!Vect_is_3d(&seed_Map))
414 G_fatal_error(_("Vector map <%s> is not 3D"), seed_opt->answer);
415
416 n_seeds = Vect_get_num_primitives(&seed_Map, GV_POINT);
417 }
418 if (flowacc_opt->answer || (!seed_opt->answer && flowlines_opt->answer)) {
419 if (flowacc_opt->answer)
420 n_seeds += region.cols * region.rows * region.depths;
421 else {
422 n_seeds += ceil(region.cols / (double)skip[0]) *
423 ceil(region.rows / (double)skip[1]) *
424 ceil(region.depths / (double)skip[2]);
425 }
426 }
427 G_debug(1, "Number of seeds is %d", n_seeds);
428
429 seed_count = 0;
430 cat = 1;
431 if (seed_opt->answer) {
432
433 seed_points = Vect_new_line_struct();
434 seed_cats = Vect_new_cats_struct();
435
436 /* compute flowlines from vector seed map */
437 while (TRUE) {
438 ltype = Vect_read_next_line(&seed_Map, seed_points, seed_cats);
439 if (ltype == -1) {
440 Vect_close(&seed_Map);
441 G_fatal_error(_("Error during reading seed vector map"));
442 }
443 else if (ltype == -2) {
444 break;
445 }
446 else if (ltype == GV_POINT) {
447 seed.x = seed_points->x[0];
448 seed.y = seed_points->y[0];
449 seed.z = seed_points->z[0];
450 seed.flowline = TRUE;
451 seed.flowaccum = FALSE;
452 }
453 G_percent(seed_count, n_seeds, 1);
454 if (integration.direction_type == FLOWDIR_UP ||
455 integration.direction_type == FLOWDIR_BOTH) {
456 integration.actual_direction = FLOWDIR_UP;
457 compute_flowline(®ion, &seed, &gradient_info, flowacc, sampled,
458 &integration, &fl_map, fl_cats, fl_points,
459 &cat, if_table, finfo, driver);
460 }
461 if (integration.direction_type == FLOWDIR_DOWN ||
462 integration.direction_type == FLOWDIR_BOTH) {
463 integration.actual_direction = FLOWDIR_DOWN;
464 compute_flowline(®ion, &seed, &gradient_info, flowacc, sampled,
465 &integration, &fl_map, fl_cats, fl_points,
466 &cat, if_table, finfo, driver);
467 }
468 seed_count++;
469 }
470
471 Vect_destroy_line_struct(seed_points);
472 Vect_destroy_cats_struct(seed_cats);
473 Vect_close(&seed_Map);
474 }
475 if (flowacc_opt->answer || (!seed_opt->answer && flowlines_opt->answer)) {
476 /* compute flowlines from points on grid */
477 for (r = region.rows; r > 0; r--) {
478 for (c = 0; c < region.cols; c++) {
479 for (d = 0; d < region.depths; d++) {
480 seed.x =
481 region.west + c * region.ew_res + region.ew_res / 2;
482 seed.y =
483 region.south + r * region.ns_res - region.ns_res / 2;
484 seed.z =
485 region.bottom + d * region.tb_res + region.tb_res / 2;
486 seed.flowline = FALSE;
487 seed.flowaccum = FALSE;
488 if (flowacc_opt->answer)
489 seed.flowaccum = TRUE;
490
491 if (flowlines_opt->answer && !seed_opt->answer &&
492 (c % skip[0] == 0) && (r % skip[1] == 0) && (d % skip[2] == 0))
493 seed.flowline = TRUE;
494
495 if (seed.flowaccum || seed.flowline) {
496 G_percent(seed_count, n_seeds, 1);
497
498 if (integration.direction_type == FLOWDIR_UP ||
499 integration.direction_type == FLOWDIR_BOTH) {
500 integration.actual_direction = FLOWDIR_UP;
501 compute_flowline(®ion, &seed, &gradient_info,
502 flowacc, sampled, &integration, &fl_map,
503 fl_cats, fl_points, &cat,
504 if_table, finfo, driver);
505 }
506 if (integration.direction_type == FLOWDIR_DOWN ||
507 integration.direction_type == FLOWDIR_BOTH) {
508 integration.actual_direction = FLOWDIR_DOWN;
509 compute_flowline(®ion, &seed, &gradient_info,
510 flowacc, sampled, &integration, &fl_map,
511 fl_cats, fl_points, &cat,
512 if_table, finfo, driver);
513 }
514 seed_count++;
515 }
516 }
517 }
518 }
519 }
520 G_percent(1, 1, 1);
521 if (flowlines_opt->answer) {
522 if (if_table) {
523 db_commit_transaction(driver);
524 db_close_database_shutdown_driver(driver);
525 }
526 Vect_destroy_line_struct(fl_points);
527 Vect_destroy_cats_struct(fl_cats);
528 Vect_build(&fl_map);
529 Vect_close(&fl_map);
530 }
531
532 if (flowacc_opt->answer)
533 Rast3d_close(flowacc);
534
535
536 return EXIT_SUCCESS;
537 }
538