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