GRASS 8 Programmer's Manual  8.5.0dev(2025)-d272f2f102
net_build.c
Go to the documentation of this file.
1 /*!
2  * \file lib/vector/Vlib/net_build.c
3  *
4  * \brief Vector library - related fns for vector network building
5  *
6  * Higher level functions for reading/writing/manipulating vectors.
7  *
8  * (C) 2001-2009, 2014 by the GRASS Development Team
9  *
10  * This program is free software under the GNU General Public License
11  * (>=v2). Read the file COPYING that comes with GRASS for details.
12  *
13  * \author Radim Blazek
14  * \author Stepan Turek stepan.turek seznam.cz (turns support)
15  */
16 
17 #include <grass/dbmi.h>
18 #include <grass/vector.h>
19 #include <grass/glocale.h>
20 
21 /*!
22  \brief Build network graph with turntable.
23 
24  Internal format for edge costs is integer, costs are multiplied
25  before conversion to int by 1000 and for lengths LL without geo flag by
26  1000000. The same multiplication factor is used for nodes. Costs in database
27  column may be 'integer' or 'double precision' number >= 0 or -1 for infinity
28  i.e. arc or node is closed and cannot be traversed If record in table is not
29  found for arcs, costs for arc are set to 0. If record in table is not found
30  for node, costs for node are set to 0.
31 
32  \param Map vector map
33  \param ltype line type for arcs
34  \param afield arc costs field (if 0, use length)
35  \param nfield node costs field (if 0, do not use node costs)
36  \param tfield field where turntable is attached
37  \param tucfield field with unique categories used in the turntable
38  \param afcol column with forward costs for arc
39  \param abcol column with backward costs for arc (if NULL, back costs =
40  forward costs)
41  \param ncol column with costs for nodes (if NULL, do not use
42  node costs)
43  \param geo use geodesic calculation for length (LL)
44  \param algorithm not used (in future code for algorithm)
45 
46  \return 0 on success, 1 on error
47  */
48 int Vect_net_ttb_build_graph(struct Map_info *Map, int ltype, int afield,
49  int nfield, int tfield, int tucfield,
50  const char *afcol, const char *abcol,
51  const char *ncol, int geo, int algorithm UNUSED)
52 {
53  /* TODO very long function, split into smaller ones */
54  int i, j, from, to, line, nlines, nnodes, ret, type, cat, skipped, cfound;
55  struct line_pnts *Points;
56  struct line_cats *Cats;
57  double dcost, bdcost, ll;
58  int cost, bcost;
59  dglGraph_s *gr;
60  dglInt32_t opaqueset[16] = {360000, 0, 0, 0, 0, 0, 0, 0,
61  0, 0, 0, 0, 0, 0, 0, 0};
62  struct field_info *Fi = NULL;
63  dbDriver *driver = NULL;
64  dbDriver *ttbdriver = NULL;
65  dbHandle handle;
66  dbString stmt;
67  dbColumn *Column = NULL;
68  dbCatValArray fvarr, bvarr;
69  int fctype = 0, bctype = 0, nrec, nturns;
70 
71  int ln_cat, nnode_lns, i_line, line_id, i_virt_edge;
72  struct line_cats *ln_Cats;
73  double x, y, z;
74  struct bound_box box;
75  struct boxlist *List;
76 
77  dglInt32_t dgl_cost;
78 
79  /*TODO attributes of turntable should be stored in one place */
80  const char *tcols[] = {"cat", "ln_from", "ln_to", "cost", "isec", NULL};
81  dbCatValArray tvarrs[5];
82  int tctype[5] = {0};
83  int tucfield_idx;
84 
85  int t, f;
86  int node_pt_id, turn_cat, tucfound;
87  int isec;
88 
89  /* TODO int costs -> double (waiting for dglib) */
90  G_debug(1,
91  "Vect_net_ttb_build_graph(): "
92  "ltype = %d, afield = %d, nfield = %d, tfield = %d, tucfield = %d ",
93  ltype, afield, nfield, tfield, tucfield);
94  G_debug(1, " afcol = %s, abcol = %s, ncol = %s", afcol, abcol, ncol);
95 
96  G_message(_("Building graph..."));
97 
98  Map->dgraph.line_type = ltype;
99 
100  Points = Vect_new_line_struct();
101  Cats = Vect_new_cats_struct();
102 
103  ll = 0;
104  if (G_projection() == 3)
105  ll = 1; /* LL */
106 
107  if (afcol == NULL && ll && !geo)
108  Map->dgraph.cost_multip = 1000000;
109  else
110  Map->dgraph.cost_multip = 1000;
111 
112  nlines = Vect_get_num_lines(Map);
113  nnodes = Vect_get_num_nodes(Map);
114 
115  gr = &(Map->dgraph.graph_s);
116 
117  /* Allocate space for costs, later replace by functions reading costs from
118  * graph */
119  Map->dgraph.edge_fcosts = (double *)G_malloc((nlines + 1) * sizeof(double));
120  Map->dgraph.edge_bcosts = (double *)G_malloc((nlines + 1) * sizeof(double));
121  Map->dgraph.node_costs = (double *)G_malloc((nnodes + 1) * sizeof(double));
122 
123  /* Set to -1 initially */
124  for (i = 1; i <= nlines; i++) {
125  Map->dgraph.edge_fcosts[i] = -1; /* forward */
126  Map->dgraph.edge_bcosts[i] = -1; /* backward */
127  }
128  for (i = 1; i <= nnodes; i++) {
129  Map->dgraph.node_costs[i] = 0;
130  }
131 
132  dglInitialize(gr, (dglByte_t)1, sizeof(dglInt32_t), (dglInt32_t)0,
133  opaqueset);
134 
135  if (gr == NULL)
136  G_fatal_error(_("Unable to build network graph"));
137 
138  db_init_handle(&handle);
139  db_init_string(&stmt);
140 
141  if (abcol != NULL && afcol == NULL)
142  G_fatal_error(_("Forward costs column not specified"));
143 
144  /* --- Add arcs --- */
145  /* Open db connection */
146 
147  /* Get field info */
148  if (tfield < 1)
149  G_fatal_error(_("Turntable field < 1"));
150  Fi = Vect_get_field(Map, tfield);
151  if (Fi == NULL)
152  G_fatal_error(_("Database connection not defined for layer %d"),
153  tfield);
154 
155  /* Open database */
156  ttbdriver = db_start_driver_open_database(Fi->driver, Fi->database);
157  if (ttbdriver == NULL)
158  G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
159  Fi->database, Fi->driver);
160 
161  i = 0;
162  while (tcols[i]) {
163  /* Load costs to array */
164  if (db_get_column(ttbdriver, Fi->table, tcols[i], &Column) != DB_OK)
165  G_fatal_error(_("Turntable column <%s> not found in table <%s>"),
166  tcols[i], Fi->table);
167 
168  tctype[i] = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
169  db_free_column(Column);
170 
171  if ((tctype[i] == DB_C_TYPE_INT || tctype[i] == DB_C_TYPE_DOUBLE) &&
172  !strcmp(tcols[i], "cost"))
173  ;
174  else if (tctype[i] == DB_C_TYPE_INT)
175  ;
176  else
178  _("Data type of column <%s> not supported (must be numeric)"),
179  tcols[i]);
180 
181  db_CatValArray_init(&tvarrs[i]);
182  nturns = db_select_CatValArray(ttbdriver, Fi->table, Fi->key, tcols[i],
183  NULL, &tvarrs[i]);
184  ++i;
185  }
186 
188 
189  G_debug(1, "forward costs: nrec = %d", nturns);
190 
191  /* Set node attributes */
192  G_message("Register nodes");
193  if (ncol != NULL) {
194 
195  G_debug(2, "Set nodes' costs");
196  if (nfield < 1)
197  G_fatal_error("Node field < 1");
198 
199  G_message(_("Setting node costs..."));
200 
201  Fi = Vect_get_field(Map, nfield);
202  if (Fi == NULL)
203  G_fatal_error(_("Database connection not defined for layer %d"),
204  nfield);
205 
207  if (driver == NULL)
208  G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
209  Fi->database, Fi->driver);
210 
211  /* Load costs to array */
212  if (db_get_column(driver, Fi->table, ncol, &Column) != DB_OK)
213  G_fatal_error(_("Column <%s> not found in table <%s>"), ncol,
214  Fi->table);
215 
216  fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
217  db_free_column(Column);
218 
219  if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
221  _("Data type of column <%s> not supported (must be numeric)"),
222  ncol);
223 
224  db_CatValArray_init(&fvarr);
225 
226  nrec = db_select_CatValArray(driver, Fi->table, Fi->key, ncol, NULL,
227  &fvarr);
228  G_debug(1, "node costs: nrec = %d", nrec);
230 
231  tucfield_idx = Vect_cidx_get_field_index(Map, tucfield);
232  }
233 
234  List = Vect_new_boxlist(0);
235  ln_Cats = Vect_new_cats_struct();
236 
237  G_message("Building turns graph...");
238 
239  i_virt_edge = -1;
240  for (i = 1; i <= nnodes; i++) {
241  /* TODO: what happens if we set attributes of non existing node (skipped
242  * lines, nodes without lines) */
243 
244  /* select points at node */
245  Vect_get_node_coor(Map, i, &x, &y, &z);
246  box.E = box.W = x;
247  box.N = box.S = y;
248  box.T = box.B = z;
249  Vect_select_lines_by_box(Map, &box, GV_POINT, List);
250 
251  G_debug(2, " node = %d nlines = %d", i, List->n_values);
252  cfound = 0;
253  dcost = 0;
254  tucfound = 0;
255 
256  for (j = 0; j < List->n_values; j++) {
257  line = List->id[j];
258  G_debug(2, " line (%d) = %d", j, line);
259  type = Vect_read_line(Map, NULL, Cats, line);
260  if (!(type & GV_POINT))
261  continue;
262  /* get node column costs */
263  if (ncol != NULL && !cfound &&
264  Vect_cat_get(Cats, nfield,
265  &cat)) { /* point with category of field found */
266  /* Set costs */
267  if (fctype == DB_C_TYPE_INT) {
268  ret = db_CatValArray_get_value_int(&fvarr, cat, &cost);
269  dcost = cost;
270  }
271  else { /* DB_C_TYPE_DOUBLE */
272  ret = db_CatValArray_get_value_double(&fvarr, cat, &dcost);
273  }
274  if (ret != DB_OK) {
275  G_warning(
276  _("Database record for node %d (cat = %d) not found "
277  "(cost set to 0)"),
278  i, cat);
279  }
280  cfound = 1;
281  Map->dgraph.node_costs[i] = dcost;
282  }
283 
284  /* add virtual nodes and lines, which represents the intersections
285  there are added two nodes for every intersection, which are
286  linked with the nodes (edges in primal graph). the positive node
287  - when we are going from the intersection the negative node -
288  when we are going to the intersection
289 
290  TODO There are more possible approaches in virtual nodes
291  management. We can also add and remove them dynamically as they
292  are needed for analysis when Vect_net_ttb_shortest_path is called
293  (problem of flattening graph).
294  Currently this static solution was chosen, because it cost
295  time only when graph is build. However it costs more memory
296  space. For Dijkstra algorithm this expansion should not be
297  serious problem because we can only get into positive node or go
298  from the negative node.
299 
300  */
301 
302  ret = Vect_cat_get(Cats, tucfield, &cat);
303  if (!tucfound && ret) { /* point with category of field found */
304  /* find lines which belongs to the intersection */
305  nnode_lns = Vect_get_node_n_lines(Map, i);
306 
307  for (i_line = 0; i_line < nnode_lns; i_line++) {
308 
309  line_id = Vect_get_node_line(Map, i, i_line);
310  Vect_read_line(Map, NULL, ln_Cats, abs(line_id));
311  Vect_cat_get(ln_Cats, tucfield, &ln_cat);
312 
313  if (line_id < 0)
314  ln_cat *= -1;
315  f = cat * 2;
316 
317  if (ln_cat < 0)
318  t = ln_cat * -2 + 1;
319  else
320  t = ln_cat * 2;
321 
322  G_debug(
323  5,
324  "Add arc %d for virtual node from %d to %d cost = %d",
325  i_virt_edge, f, t, 0);
326 
327  /* positive, start virtual node */
328  ret = dglAddEdge(gr, (dglInt32_t)f, (dglInt32_t)t,
329  (dglInt32_t)0, (dglInt32_t)(i_virt_edge));
330  if (ret < 0)
331  G_fatal_error(_("Cannot add network arc for virtual "
332  "node connection."));
333 
334  t = cat * 2 + 1;
335  i_virt_edge--;
336 
337  if (-ln_cat < 0)
338  f = ln_cat * 2 + 1;
339  else
340  f = ln_cat * -2;
341 
342  G_debug(
343  5,
344  "Add arc %d for virtual node from %d to %d cost = %d",
345  i_virt_edge, f, t, 0);
346 
347  /* negative, destination virtual node */
348  ret = dglAddEdge(gr, (dglInt32_t)f, (dglInt32_t)t,
349  (dglInt32_t)0, (dglInt32_t)(i_virt_edge));
350  if (ret < 0)
351  G_fatal_error(_("Cannot add network arc for virtual "
352  "node connection."));
353 
354  i_virt_edge--;
355  }
356  tucfound++;
357  }
358  else if (ret)
359  tucfound++;
360  }
361 
362  if (tucfound > 1)
363  G_warning(_("There exists more than one point of node <%d> with "
364  "unique category field <%d>.\n"
365  "The unique categories layer is not valid therefore "
366  "you will probably get incorrect results."),
367  tucfield, i);
368 
369  if (ncol != NULL && !cfound)
370  G_debug(
371  2,
372  "Category of field %d is not attached to any points in node %d"
373  "(costs set to 0)",
374  nfield, i);
375  }
376 
377  Vect_destroy_cats_struct(ln_Cats);
378 
379  for (i = 1; i <= nturns; i++) {
380  /* select points at node */
381 
382  /* TODO use cursors */
383  db_CatValArray_get_value_int(&tvarrs[0], i, &turn_cat);
384 
385  db_CatValArray_get_value_int(&tvarrs[1], i, &from);
386  db_CatValArray_get_value_int(&tvarrs[2], i, &to);
387 
388  db_CatValArray_get_value_int(&tvarrs[4], i, &isec);
389  dcost = 0.0;
390  if (ncol != NULL) {
391  /* TODO optimization do not do it for every turn in intersection
392  * again */
393  if (Vect_cidx_find_next(Map, tucfield_idx, isec, GV_POINT, 0, &type,
394  &node_pt_id) == -1) {
395  G_warning(
396  _("Unable to find point representing intersection <%d> in "
397  "unique categories field <%d>.\n"
398  "Cost for the intersection was set to 0.\n"
399  "The unique categories layer is not valid therefore you "
400  "will probably get incorrect results."),
401  isec, tucfield);
402  }
403  else {
404  Vect_read_line(Map, Points, Cats, node_pt_id);
405 
406  node_pt_id = Vect_find_node(Map, *Points->x, *Points->y,
407  *Points->z, 0.0, WITHOUT_Z);
408 
409  if (node_pt_id == 0) {
410  G_warning(
411  _("Unable to find node for point representing "
412  "intersection <%d> in unique categories field <%d>.\n"
413  "Cost for the intersection was set to 0.\n"
414  "The unique categories layer is not valid therefore "
415  "you will probably get incorrect results."),
416  isec, tucfield);
417  }
418  else {
419  G_debug(2, " node = %d", node_pt_id);
420  dcost = Map->dgraph.node_costs[node_pt_id];
421  }
422  }
423  }
424 
425  G_debug(2, "Set node's cost to %f", dcost);
426 
427  if (dcost >= 0) {
428  /* Set costs from turntable */
429  if (tctype[3] == DB_C_TYPE_INT) {
430  ret = db_CatValArray_get_value_int(&tvarrs[3], i, &cost);
431  dcost = cost;
432  }
433  else /* DB_C_TYPE_DOUBLE */
434  ret = db_CatValArray_get_value_double(&tvarrs[3], i, &dcost);
435 
436  if (ret != DB_OK) {
437  G_warning(
438  _("Database record for turn with cat = %d is not found. "
439  "(The turn was skipped."),
440  i);
441  continue;
442  }
443 
444  if (dcost >= 0) {
445 
446  if (ncol != NULL)
447  cost = (Map->dgraph.node_costs[node_pt_id] + dcost) *
449  else
450  cost = dcost * (dglInt32_t)Map->dgraph.cost_multip;
451 
452  /* dglib does not like negative id's of nodes */
453  if (from < 0)
454  f = from * -2 + 1;
455  else
456  f = from * 2;
457 
458  if (to < 0)
459  t = to * -2 + 1;
460  else
461  t = to * 2;
462 
463  G_debug(5, "Add arc/turn %d for turn from %d to %d cost = %d",
464  turn_cat, f, t, cost);
465 
466  ret = dglAddEdge(gr, (dglInt32_t)f, (dglInt32_t)t,
467  (dglInt32_t)cost, (dglInt32_t)(turn_cat));
468 
469  if (ret < 0)
471  _("Cannot add network arc representing turn."));
472  }
473  }
474  }
475 
476  Vect_destroy_boxlist(List);
477 
478  i = 0;
479  while (tcols[i]) {
480  db_CatValArray_free(&tvarrs[i]);
481  ++i;
482  }
483 
484  if (ncol != NULL) {
486  db_CatValArray_free(&fvarr);
487  }
488 
489  /* Open db connection */
490  if (afcol != NULL) {
491  /* Get field info */
492  if (afield < 1)
493  G_fatal_error(_("Arc field < 1"));
494  Fi = Vect_get_field(Map, afield);
495  if (Fi == NULL)
496  G_fatal_error(_("Database connection not defined for layer %d"),
497  afield);
498 
499  /* Open database */
501  if (driver == NULL)
502  G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
503  Fi->database, Fi->driver);
504 
505  /* Load costs to array */
506  if (db_get_column(driver, Fi->table, afcol, &Column) != DB_OK)
507  G_fatal_error(_("Column <%s> not found in table <%s>"), afcol,
508  Fi->table);
509 
510  fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
511  db_free_column(Column);
512 
513  if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
515  _("Data type of column <%s> not supported (must be numeric)"),
516  afcol);
517 
518  db_CatValArray_init(&fvarr);
519  nrec = db_select_CatValArray(driver, Fi->table, Fi->key, afcol, NULL,
520  &fvarr);
521  G_debug(1, "forward costs: nrec = %d", nrec);
522 
523  if (abcol != NULL) {
524  if (db_get_column(driver, Fi->table, abcol, &Column) != DB_OK)
525  G_fatal_error(_("Column <%s> not found in table <%s>"), abcol,
526  Fi->table);
527 
528  bctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
529  db_free_column(Column);
530 
531  if (bctype != DB_C_TYPE_INT && bctype != DB_C_TYPE_DOUBLE)
532  G_fatal_error(_("Data type of column <%s> not supported (must "
533  "be numeric)"),
534  abcol);
535 
536  db_CatValArray_init(&bvarr);
537  nrec = db_select_CatValArray(driver, Fi->table, Fi->key, abcol,
538  NULL, &bvarr);
539  G_debug(1, "backward costs: nrec = %d", nrec);
540  }
542  }
543 
544  skipped = 0;
545 
546  G_message(_("Registering arcs..."));
547 
548  for (i = 1; i <= nlines; i++) {
549  G_percent(i, nlines, 1); /* must be before any continue */
550 
551  type = Vect_read_line(Map, Points, Cats, i);
552  if (!(type & ltype & (GV_LINE | GV_BOUNDARY)))
553  continue;
554 
555  Vect_get_line_nodes(Map, i, &from, &to);
556 
557  dcost = bdcost = 0;
558 
559  cfound = Vect_cat_get(Cats, tucfield, &cat);
560  if (!cfound)
561  continue;
562 
563  if (cfound > 1)
564  G_warning(_("Line with id <%d> has more unique categories defined "
565  "in field <%d>.\n"
566  "The unique categories layer is not valid therefore "
567  "you will probably get incorrect results."),
568  i, tucfield);
569 
570  if (afcol != NULL) {
571  if (!(Vect_cat_get(Cats, afield, &cat))) {
572  G_debug(2,
573  "Category of field %d not attached to the line %d -> "
574  "cost was set to 0",
575  afield, i);
576  skipped += 2; /* Both directions */
577  }
578  else {
579  if (fctype == DB_C_TYPE_INT) {
580  ret = db_CatValArray_get_value_int(&fvarr, cat, &cost);
581  dcost = cost;
582  }
583  else { /* DB_C_TYPE_DOUBLE */
584  ret = db_CatValArray_get_value_double(&fvarr, cat, &dcost);
585  }
586  if (ret != DB_OK) {
587  G_warning(_("Database record for line %d (cat = %d, "
588  "forward/both direction(s)) not found "
589  "(cost was set to 0)"),
590  i, cat);
591  }
592 
593  if (abcol != NULL) {
594  if (bctype == DB_C_TYPE_INT) {
595  ret = db_CatValArray_get_value_int(&bvarr, cat, &bcost);
596  bdcost = bcost;
597  }
598  else { /* DB_C_TYPE_DOUBLE */
599  ret = db_CatValArray_get_value_double(&bvarr, cat,
600  &bdcost);
601  }
602  if (ret != DB_OK) {
603  G_warning(_("Database record for line %d (cat = %d, "
604  "backward direction) not found"
605  "(cost was set to 0)"),
606  i, cat);
607  }
608  }
609  else
610  bdcost = dcost;
611 
612  Vect_cat_get(Cats, tucfield, &cat);
613  }
614  }
615  else {
616  if (ll) {
617  if (geo)
618  dcost = Vect_line_geodesic_length(Points);
619  else
620  dcost = Vect_line_length(Points);
621  }
622  else
623  dcost = Vect_line_length(Points);
624 
625  bdcost = dcost;
626  }
627 
628  cost = (dglInt32_t)Map->dgraph.cost_multip * dcost;
629  dgl_cost = cost;
630 
631  cat = cat * 2;
632 
633  G_debug(5, "Setinng node %d cost: %d", cat, cost);
634  dglNodeSet_Attr(gr, dglGetNode(gr, (dglInt32_t)cat), &dgl_cost);
635 
636  Map->dgraph.edge_fcosts[i] = dcost;
637 
638  cost = (dglInt32_t)Map->dgraph.cost_multip * bdcost;
639  dgl_cost = cost;
640 
641  ++cat;
642 
643  G_debug(5, "Setinng node %d cost: %d", cat, cost);
644  dglNodeSet_Attr(gr, dglGetNode(gr, (dglInt32_t)cat), &dgl_cost);
645 
646  Map->dgraph.edge_bcosts[i] = bdcost;
647  }
648 
649  if (afcol != NULL && skipped > 0)
650  G_debug(2, "%d lines missing category of field %d skipped", skipped,
651  afield);
652 
653  if (afcol != NULL) {
655  db_CatValArray_free(&fvarr);
656 
657  if (abcol != NULL) {
658  db_CatValArray_free(&bvarr);
659  }
660  }
662 
663  Vect_destroy_line_struct(Points);
665 
666  G_message(_("Flattening the graph..."));
667  ret = dglFlatten(gr);
668  if (ret < 0)
669  G_fatal_error(_("GngFlatten error"));
670 
671  /* init SP cache */
672  /* disable to debug dglib cache */
673  dglInitializeSPCache(gr, &(Map->dgraph.spCache));
674 
675  G_message(_("Graph was built"));
676 
677  return 0;
678 }
679 
680 /*!
681  \brief Build network graph.
682 
683  Internal format for edge costs is integer, costs are multiplied
684  before conversion to int by 1000 and for lengths LL without geo flag by
685  1000000. The same multiplication factor is used for nodes. Costs in database
686  column may be 'integer' or 'double precision' number >= 0 or -1 for infinity
687  i.e. arc or node is closed and cannot be traversed If record in table is not
688  found for arcs, arc is skip. If record in table is not found for node, costs
689  for node are set to 0.
690 
691  \param Map vector map
692  \param ltype line type for arcs
693  \param afield arc costs field (if 0, use length)
694  \param nfield node costs field (if 0, do not use node costs)
695  \param afcol column with forward costs for arc
696  \param abcol column with backward costs for arc (if NULL, back costs =
697  forward costs), \param ncol column with costs for nodes (if NULL, do not use
698  node costs), \param geo use geodesic calculation for length (LL), \param
699  version graph version to create (1, 2, 3)
700 
701  \return 0 on success, 1 on error
702  */
703 int Vect_net_build_graph(struct Map_info *Map, int ltype, int afield,
704  int nfield, const char *afcol, const char *abcol,
705  const char *ncol, int geo, int version)
706 {
707  /* TODO long function, split into smaller ones */
708  int i, j, from, to, line, nlines, nnodes, ret, type, cat, skipped, cfound;
709  int dofw, dobw;
710  struct line_pnts *Points;
711  struct line_cats *Cats;
712  double dcost, bdcost, ll;
713  int cost, bcost;
714  dglGraph_s *gr;
715  dglInt32_t dgl_cost;
716  dglInt32_t opaqueset[16] = {360000, 0, 0, 0, 0, 0, 0, 0,
717  0, 0, 0, 0, 0, 0, 0, 0};
718  struct field_info *Fi = NULL;
719  dbDriver *driver = NULL;
720  dbHandle handle;
721  dbString stmt;
722  dbColumn *Column = NULL;
723  dbCatValArray fvarr, bvarr;
724  int fctype = 0, bctype = 0, nrec;
725 
726  /* TODO int costs -> double (waiting for dglib) */
727  G_debug(1, "Vect_net_build_graph(): ltype = %d, afield = %d, nfield = %d",
728  ltype, afield, nfield);
729  G_debug(1, " afcol = %s, abcol = %s, ncol = %s", afcol, abcol, ncol);
730 
731  G_message(_("Building graph..."));
732 
733  Map->dgraph.line_type = ltype;
734 
735  Points = Vect_new_line_struct();
736  Cats = Vect_new_cats_struct();
737 
738  ll = 0;
739  if (G_projection() == 3)
740  ll = 1; /* LL */
741 
742  if (afcol == NULL && ll && !geo)
743  Map->dgraph.cost_multip = 1000000;
744  else
745  Map->dgraph.cost_multip = 1000;
746 
747  nlines = Vect_get_num_lines(Map);
748  nnodes = Vect_get_num_nodes(Map);
749 
750  gr = &(Map->dgraph.graph_s);
751 
752  /* Allocate space for costs, later replace by functions reading costs from
753  * graph */
754  Map->dgraph.edge_fcosts = (double *)G_malloc((nlines + 1) * sizeof(double));
755  Map->dgraph.edge_bcosts = (double *)G_malloc((nlines + 1) * sizeof(double));
756  Map->dgraph.node_costs = (double *)G_malloc((nnodes + 1) * sizeof(double));
757  /* Set to -1 initially */
758  for (i = 1; i <= nlines; i++) {
759  Map->dgraph.edge_fcosts[i] = -1; /* forward */
760  Map->dgraph.edge_bcosts[i] = -1; /* backward */
761  }
762  for (i = 1; i <= nnodes; i++) {
763  Map->dgraph.node_costs[i] = 0;
764  }
765 
766  if (version < 1 || version > 3)
767  version = 1;
768 
769  if (ncol != NULL)
770  dglInitialize(gr, (dglByte_t)version, sizeof(dglInt32_t), (dglInt32_t)0,
771  opaqueset);
772  else
773  dglInitialize(gr, (dglByte_t)version, (dglInt32_t)0, (dglInt32_t)0,
774  opaqueset);
775 
776  if (gr == NULL)
777  G_fatal_error(_("Unable to build network graph"));
778 
779  db_init_handle(&handle);
780  db_init_string(&stmt);
781 
782  if (abcol != NULL && afcol == NULL)
783  G_fatal_error(_("Forward costs column not specified"));
784 
785  /* --- Add arcs --- */
786  /* Open db connection */
787  if (afcol != NULL) {
788  /* Get field info */
789  if (afield < 1)
790  G_fatal_error(_("Arc field < 1"));
791  Fi = Vect_get_field(Map, afield);
792  if (Fi == NULL)
793  G_fatal_error(_("Database connection not defined for layer %d"),
794  afield);
795 
796  /* Open database */
798  if (driver == NULL)
799  G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
800  Fi->database, Fi->driver);
801 
802  /* Load costs to array */
803  if (db_get_column(driver, Fi->table, afcol, &Column) != DB_OK)
804  G_fatal_error(_("Column <%s> not found in table <%s>"), afcol,
805  Fi->table);
806 
807  fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
808  db_free_column(Column);
809 
810  if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
812  _("Data type of column <%s> not supported (must be numeric)"),
813  afcol);
814 
815  db_CatValArray_init(&fvarr);
816  nrec = db_select_CatValArray(driver, Fi->table, Fi->key, afcol, NULL,
817  &fvarr);
818  G_debug(1, "forward costs: nrec = %d", nrec);
819 
820  if (abcol != NULL) {
821  if (db_get_column(driver, Fi->table, abcol, &Column) != DB_OK)
822  G_fatal_error(_("Column <%s> not found in table <%s>"), abcol,
823  Fi->table);
824 
825  bctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
826  db_free_column(Column);
827 
828  if (bctype != DB_C_TYPE_INT && bctype != DB_C_TYPE_DOUBLE)
829  G_fatal_error(_("Data type of column <%s> not supported (must "
830  "be numeric)"),
831  abcol);
832 
833  db_CatValArray_init(&bvarr);
834  nrec = db_select_CatValArray(driver, Fi->table, Fi->key, abcol,
835  NULL, &bvarr);
836  G_debug(1, "backward costs: nrec = %d", nrec);
837  }
839  }
840 
841  skipped = 0;
842 
843  G_message(_("Registering arcs..."));
844 
845  for (i = 1; i <= nlines; i++) {
846  G_percent(i, nlines, 1); /* must be before any continue */
847  dofw = dobw = 1;
848  type = Vect_read_line(Map, Points, Cats, i);
849  if (!(type & ltype & (GV_LINE | GV_BOUNDARY)))
850  continue;
851 
852  Vect_get_line_nodes(Map, i, &from, &to);
853 
854  if (afcol != NULL) {
855  if (!(Vect_cat_get(Cats, afield, &cat))) {
856  G_debug(2,
857  "Category of field %d not attached to the line %d -> "
858  "line skipped",
859  afield, i);
860  skipped += 2; /* Both directions */
861  continue;
862  }
863  else {
864  if (fctype == DB_C_TYPE_INT) {
865  ret = db_CatValArray_get_value_int(&fvarr, cat, &cost);
866  dcost = cost;
867  }
868  else { /* DB_C_TYPE_DOUBLE */
869  ret = db_CatValArray_get_value_double(&fvarr, cat, &dcost);
870  }
871  if (ret != DB_OK) {
872  G_warning(_("Database record for line %d (cat = %d, "
873  "forward/both direction(s)) not found "
874  "(forward/both direction(s) of line skipped)"),
875  i, cat);
876  dofw = 0;
877  }
878 
879  if (abcol != NULL) {
880  if (bctype == DB_C_TYPE_INT) {
881  ret = db_CatValArray_get_value_int(&bvarr, cat, &bcost);
882  bdcost = bcost;
883  }
884  else { /* DB_C_TYPE_DOUBLE */
885  ret = db_CatValArray_get_value_double(&bvarr, cat,
886  &bdcost);
887  }
888  if (ret != DB_OK) {
889  G_warning(_("Database record for line %d (cat = %d, "
890  "backward direction) not found"
891  "(direction of line skipped)"),
892  i, cat);
893  dobw = 0;
894  }
895  }
896  else {
897  if (dofw)
898  bdcost = dcost;
899  else
900  dobw = 0;
901  }
902  }
903  }
904  else {
905  if (ll) {
906  if (geo)
907  dcost = Vect_line_geodesic_length(Points);
908  else
909  dcost = Vect_line_length(Points);
910  }
911  else
912  dcost = Vect_line_length(Points);
913 
914  bdcost = dcost;
915  }
916  if (dofw && dcost != -1) {
917  cost = (dglInt32_t)Map->dgraph.cost_multip * dcost;
918  G_debug(5, "Add arc %d from %d to %d cost = %d", i, from, to, cost);
919  ret = dglAddEdge(gr, (dglInt32_t)from, (dglInt32_t)to,
920  (dglInt32_t)cost, (dglInt32_t)i);
921  Map->dgraph.edge_fcosts[i] = dcost;
922  if (ret < 0)
923  G_fatal_error("Cannot add network arc");
924  }
925 
926  G_debug(5, "bdcost = %f edge_bcosts = %f", bdcost,
927  Map->dgraph.edge_bcosts[i]);
928  if (dobw && bdcost != -1) {
929  bcost = (dglInt32_t)Map->dgraph.cost_multip * bdcost;
930  G_debug(5, "Add arc %d from %d to %d bcost = %d", -i, to, from,
931  bcost);
932  ret = dglAddEdge(gr, (dglInt32_t)to, (dglInt32_t)from,
933  (dglInt32_t)bcost, (dglInt32_t)-i);
934  Map->dgraph.edge_bcosts[i] = bdcost;
935  if (ret < 0)
936  G_fatal_error(_("Cannot add network arc"));
937  }
938  }
939 
940  if (afcol != NULL && skipped > 0)
941  G_debug(2, "%d lines missing category of field %d skipped", skipped,
942  afield);
943 
944  if (afcol != NULL) {
946  db_CatValArray_free(&fvarr);
947 
948  if (abcol != NULL) {
949  db_CatValArray_free(&bvarr);
950  }
951  }
952 
953  /* Set node attributes */
954  G_debug(2, "Register nodes");
955  if (ncol != NULL) {
956  double x, y, z;
957  struct bound_box box;
958  struct boxlist *List;
959 
960  List = Vect_new_boxlist(0);
961 
962  G_debug(2, "Set nodes' costs");
963  if (nfield < 1)
964  G_fatal_error("Node field < 1");
965 
966  G_message(_("Setting node costs..."));
967 
968  Fi = Vect_get_field(Map, nfield);
969  if (Fi == NULL)
970  G_fatal_error(_("Database connection not defined for layer %d"),
971  nfield);
972 
974  if (driver == NULL)
975  G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
976  Fi->database, Fi->driver);
977 
978  /* Load costs to array */
979  if (db_get_column(driver, Fi->table, ncol, &Column) != DB_OK)
980  G_fatal_error(_("Column <%s> not found in table <%s>"), ncol,
981  Fi->table);
982 
983  fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
984  db_free_column(Column);
985 
986  if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
988  _("Data type of column <%s> not supported (must be numeric)"),
989  ncol);
990 
991  db_CatValArray_init(&fvarr);
992  nrec = db_select_CatValArray(driver, Fi->table, Fi->key, ncol, NULL,
993  &fvarr);
994  G_debug(1, "node costs: nrec = %d", nrec);
996 
997  for (i = 1; i <= nnodes; i++) {
998  /* TODO: what happens if we set attributes of not existing node
999  * (skipped lines, nodes without lines) */
1000 
1001  /* select points at node */
1002  Vect_get_node_coor(Map, i, &x, &y, &z);
1003  box.E = box.W = x;
1004  box.N = box.S = y;
1005  box.T = box.B = z;
1006  Vect_select_lines_by_box(Map, &box, GV_POINT, List);
1007 
1008  G_debug(2, " node = %d nlines = %d", i, List->n_values);
1009  cfound = 0;
1010  dcost = 0;
1011 
1012  for (j = 0; j < List->n_values; j++) {
1013  line = List->id[j];
1014  G_debug(2, " line (%d) = %d", j, line);
1015  type = Vect_read_line(Map, NULL, Cats, line);
1016  if (!(type & GV_POINT))
1017  continue;
1018  if (Vect_cat_get(
1019  Cats, nfield,
1020  &cat)) { /* point with category of field found */
1021  /* Set costs */
1022  if (fctype == DB_C_TYPE_INT) {
1023  ret = db_CatValArray_get_value_int(&fvarr, cat, &cost);
1024  dcost = cost;
1025  }
1026  else { /* DB_C_TYPE_DOUBLE */
1027  ret = db_CatValArray_get_value_double(&fvarr, cat,
1028  &dcost);
1029  }
1030  if (ret != DB_OK) {
1031  G_warning(_("Database record for node %d (cat = %d) "
1032  "not found "
1033  "(cost set to 0)"),
1034  i, cat);
1035  }
1036  cfound = 1;
1037  break;
1038  }
1039  }
1040  if (!cfound) {
1041  G_debug(
1042  2,
1043  "Category of field %d not attached to any points in node %d"
1044  "(costs set to 0)",
1045  nfield, i);
1046  }
1047  if (dcost == -1) { /* closed */
1048  cost = -1;
1049  }
1050  else {
1051  cost = (dglInt32_t)Map->dgraph.cost_multip * dcost;
1052  }
1053 
1054  dgl_cost = cost;
1055  G_debug(3, "Set node's cost to %d", cost);
1056 
1057  dglNodeSet_Attr(gr, dglGetNode(gr, (dglInt32_t)i), &dgl_cost);
1058 
1059  Map->dgraph.node_costs[i] = dcost;
1060  }
1062  db_CatValArray_free(&fvarr);
1063 
1064  Vect_destroy_boxlist(List);
1065  }
1066 
1067  G_message(_("Flattening the graph..."));
1068  ret = dglFlatten(gr);
1069  if (ret < 0)
1070  G_fatal_error(_("GngFlatten error"));
1071 
1072  /* init SP cache */
1073  /* disable to debug dglib cache */
1074  dglInitializeSPCache(gr, &(Map->dgraph.spCache));
1075 
1076  G_message(_("Graph was built"));
1077 
1078  return 0;
1079 }
#define NULL
Definition: ccmath.h:32
Main header of GRASS DataBase Management Interface.
#define DB_C_TYPE_INT
Definition: dbmi.h:108
#define DB_C_TYPE_DOUBLE
Definition: dbmi.h:109
#define DB_OK
Definition: dbmi.h:71
void db_CatValArray_free(dbCatValArray *)
Free allocated dbCatValArray.
Definition: value.c:373
int db_CatValArray_get_value_int(dbCatValArray *, int, int *)
Find value (integer) by key.
void db_CatValArray_init(dbCatValArray *)
Initialize dbCatValArray.
Definition: value.c:361
int db_get_column(dbDriver *, const char *, const char *, dbColumn **)
Get column structure by table and column name.
int db_sqltype_to_Ctype(int)
Get C data type based on given SQL data type.
Definition: sqlCtype.c:24
int db_get_column_sqltype(dbColumn *)
Returns column sqltype for column.
dbDriver * db_start_driver_open_database(const char *, const char *)
Open driver/database connection.
Definition: db.c:28
int db_close_database_shutdown_driver(dbDriver *)
Close driver/database connection.
Definition: db.c:61
int db_select_CatValArray(dbDriver *, const char *, const char *, const char *, const char *, dbCatValArray *)
Select pairs key/value to array, values are sorted by key (must be integer)
void db_free_column(dbColumn *)
Frees column structure.
void db_init_handle(dbHandle *)
Initialize handle (i.e database/schema)
Definition: handle.c:23
void db_init_string(dbString *)
Initialize dbString.
Definition: string.c:25
int db_CatValArray_get_value_double(dbCatValArray *, int, double *)
Find value (double) by key.
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:61
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
#define G_malloc(n)
Definition: defs/gis.h:94
void G_message(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
int G_projection(void)
Query cartographic projection.
Definition: proj1.c:32
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition: line.c:77
int Vect_get_line_nodes(struct Map_info *, int, int *, int *)
Get line nodes.
Definition: level_two.c:304
int Vect_get_node_coor(struct Map_info *, int, double *, double *, double *)
Get node coordinates.
Definition: level_two.c:274
plus_t Vect_get_num_lines(struct Map_info *)
Fetch number of features (points, lines, boundaries, centroids) in vector map.
Definition: level_two.c:75
double Vect_line_length(const struct line_pnts *)
Calculate line length, 3D-length in case of 3D vector line.
Definition: line.c:575
int Vect_cidx_get_field_index(struct Map_info *, int)
Get layer index for given layer number.
Definition: Vlib/cindex.c:115
double Vect_line_geodesic_length(const struct line_pnts *)
Calculate line length.
Definition: line.c:602
int Vect_cidx_find_next(struct Map_info *, int, int, int, int, int *, int *)
Find next line/area id for given category, start_index and type_mask.
Definition: Vlib/cindex.c:327
struct boxlist * Vect_new_boxlist(int)
Creates and initializes a struct boxlist.
int Vect_cat_get(const struct line_cats *, int, int *)
Get first found category of given field.
void Vect_destroy_boxlist(struct boxlist *)
Frees all memory associated with a struct boxlist, including the struct itself.
void Vect_destroy_cats_struct(struct line_cats *)
Frees all memory associated with line_cats structure, including the struct itself.
int Vect_read_line(struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
struct field_info * Vect_get_field(struct Map_info *, int)
Get information about link to database (by layer number)
Definition: field.c:515
void Vect_destroy_field_info(struct field_info *)
Free a struct field_info and all memory associated with it.
Definition: field.c:633
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
int Vect_select_lines_by_box(struct Map_info *, const struct bound_box *, int, struct boxlist *)
Select lines with bounding boxes by box.
Definition: sindex.c:32
int Vect_get_node_n_lines(struct Map_info *, int)
Get number of lines for node.
Definition: level_two.c:381
plus_t Vect_get_num_nodes(struct Map_info *)
Get number of nodes in vector map.
Definition: level_two.c:34
int Vect_find_node(struct Map_info *, double, double, double, double, int)
Find the nearest node.
int Vect_get_node_line(struct Map_info *, int, int)
Get line id for node line index.
Definition: level_two.c:397
struct line_cats * Vect_new_cats_struct(void)
Creates and initializes line_cats structure.
#define GV_LINE
Definition: dig_defines.h:184
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:183
#define GV_BOUNDARY
Definition: dig_defines.h:185
#define WITHOUT_Z
2D/3D vector data
Definition: dig_defines.h:171
#define UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
Definition: gis.h:46
#define _(str)
Definition: glocale.h:10
int Vect_net_build_graph(struct Map_info *Map, int ltype, int afield, int nfield, const char *afcol, const char *abcol, const char *ncol, int geo, int version)
Build network graph.
Definition: net_build.c:703
int Vect_net_ttb_build_graph(struct Map_info *Map, int ltype, int afield, int nfield, int tfield, int tucfield, const char *afcol, const char *abcol, const char *ncol, int geo, int algorithm)
Build network graph with turntable.
Definition: net_build.c:48
double t
Definition: r_raster.c:39
dglSPCache_s spCache
Shortest path cache.
Definition: dig_structs.h:1216
int line_type
Line type used to build the graph.
Definition: dig_structs.h:1208
int cost_multip
Edge and node costs multiplicator.
Definition: dig_structs.h:1234
double * edge_fcosts
Forward costs used for graph.
Definition: dig_structs.h:1222
double * node_costs
Node costs used for graph.
Definition: dig_structs.h:1230
dglGraph_s graph_s
Graph structure.
Definition: dig_structs.h:1212
double * edge_bcosts
backward costs used for graph
Definition: dig_structs.h:1226
Vector map info.
Definition: dig_structs.h:1243
struct Graph_info dgraph
Graph info (built for network analysis)
Definition: dig_structs.h:1383
Bounding box.
Definition: dig_structs.h:64
double W
West.
Definition: dig_structs.h:80
double T
Top.
Definition: dig_structs.h:84
double S
South.
Definition: dig_structs.h:72
double N
North.
Definition: dig_structs.h:68
double E
East.
Definition: dig_structs.h:76
double B
Bottom.
Definition: dig_structs.h:88
List of bounding boxes with id.
Definition: dig_structs.h:1723
int * id
Array of ids.
Definition: dig_structs.h:1727
struct bound_box * box
Array of bounding boxes.
Definition: dig_structs.h:1731
int n_values
Number of items in the list.
Definition: dig_structs.h:1739
Definition: driver.h:27
Layer (old: field) information.
Definition: dig_structs.h:131
char * table
Name of DB table.
Definition: dig_structs.h:151
char * driver
Name of DB driver ('sqlite', 'dbf', ...)
Definition: dig_structs.h:143
char * database
Definition: dig_structs.h:147
char * key
Name of key column (usually 'cat')
Definition: dig_structs.h:155
Feature category info.
Definition: dig_structs.h:1677
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
double * y
Array of Y coordinates.
Definition: dig_structs.h:1659
double * x
Array of X coordinates.
Definition: dig_structs.h:1655
double * z
Array of Z coordinates.
Definition: dig_structs.h:1663
unsigned char dglByte_t
Definition: type.h:35
long dglInt32_t
Definition: type.h:36
int dglAddEdge(dglGraph_s *pGraph, dglInt32_t nHead, dglInt32_t nTail, dglInt32_t nCost, dglInt32_t nEdge)
void dglNodeSet_Attr(dglGraph_s *pGraph, dglInt32_t *pnNode, dglInt32_t *pnAttr)
int dglInitialize(dglGraph_s *pGraph, dglByte_t Version, dglInt32_t NodeAttrSize, dglInt32_t EdgeAttrSize, dglInt32_t *pOpaqueSet)
int dglInitializeSPCache(dglGraph_s *pGraph, dglSPCache_s *pCache)
int dglFlatten(dglGraph_s *pGraph)
dglInt32_t * dglGetNode(dglGraph_s *pGraph, dglInt32_t nNodeId)
#define x