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