1
2 /****************************************************************
3 *
4 * MODULE: v.net.distance
5 *
6 * AUTHOR(S): Daniel Bundala
7 * Markus Metz
8 *
9 * PURPOSE: Computes shortest distance via the network between
10 * two given sets of features.
11 *
12 * COPYRIGHT: (C) 2009-2010, 2012 by Daniel Bundala, and the GRASS Development Team
13 *
14 * This program is free software under the
15 * GNU General Public License (>=v2).
16 * Read the file COPYING that comes with GRASS
17 * for details.
18 *
19 ****************************************************************/
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <grass/gis.h>
24 #include <grass/vector.h>
25 #include <grass/glocale.h>
26 #include <grass/dbmi.h>
27 #include <grass/neta.h>
28
main(int argc,char * argv[])29 int main(int argc, char *argv[])
30 {
31 struct Map_info In, Out;
32 static struct line_pnts *Points, *PPoints;
33 struct line_cats *Cats, *TCats;
34 struct ilist *slist;
35 struct GModule *module; /* GRASS module for parsing arguments */
36 struct Option *map_in, *map_out;
37 struct Option *catf_opt, *fieldf_opt, *wheref_opt;
38 struct Option *catt_opt, *fieldt_opt, *wheret_opt, *typet_opt;
39 struct Option *afield_opt, *nfield_opt, *abcol, *afcol, *ncol, *atype_opt;
40 struct Flag *geo_f, *segments_f;
41 int with_z, geo, segments;
42 int atype, ttype;
43 struct varray *varrayf, *varrayt;
44 int flayer, tlayer;
45 int afield, nfield;
46 dglGraph_s *graph;
47 struct ilist *nodest;
48 int i, j, nnodes, nlines;
49 int *dst, *nodes_to_features;
50 int from_nr; /* 'from' features not reachable */
51 dglInt32_t **nxt;
52 struct line_cats **on_path;
53 char *segdir;
54 char buf[2000];
55
56 /* Attribute table */
57 dbString sql;
58 dbDriver *driver;
59 struct field_info *Fi;
60
61 /* initialize GIS environment */
62 G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
63
64 /* initialize module */
65 module = G_define_module();
66 G_add_keyword(_("vector"));
67 G_add_keyword(_("network"));
68 G_add_keyword(_("shortest path"));
69 module->label = _("Computes shortest distance via the network between "
70 "the given sets of features.");
71 module->description =
72 _("Finds the shortest paths from each 'from' point to the nearest 'to' feature "
73 "and various information about this relation are uploaded to the attribute table.");
74
75 /* Define the different options as defined in gis.h */
76 map_in = G_define_standard_option(G_OPT_V_INPUT);
77 map_out = G_define_standard_option(G_OPT_V_OUTPUT);
78
79 afield_opt = G_define_standard_option(G_OPT_V_FIELD);
80 afield_opt->key = "arc_layer";
81 afield_opt->answer = "1";
82 afield_opt->label = _("Arc layer");
83 afield_opt->guisection = _("Cost");
84
85 atype_opt = G_define_standard_option(G_OPT_V_TYPE);
86 atype_opt->key = "arc_type";
87 atype_opt->options = "line,boundary";
88 atype_opt->answer = "line,boundary";
89 atype_opt->label = _("Arc type");
90 atype_opt->guisection = _("Cost");
91
92 nfield_opt = G_define_standard_option(G_OPT_V_FIELD);
93 nfield_opt->key = "node_layer";
94 nfield_opt->answer = "2";
95 nfield_opt->label = _("Node layer");
96 nfield_opt->guisection = _("Cost");
97
98 fieldf_opt = G_define_standard_option(G_OPT_V_FIELD);
99 fieldf_opt->key = "from_layer";
100 fieldf_opt->label = _("From layer number or name");
101 fieldf_opt->guisection = _("From");
102
103 catf_opt = G_define_standard_option(G_OPT_V_CATS);
104 catf_opt->key = "from_cats";
105 catf_opt->label = _("From category values");
106 catf_opt->guisection = _("From");
107
108 wheref_opt = G_define_standard_option(G_OPT_DB_WHERE);
109 wheref_opt->key = "from_where";
110 wheref_opt->label =
111 _("From WHERE conditions of SQL statement without 'where' keyword");
112 wheref_opt->guisection = _("From");
113
114 fieldt_opt = G_define_standard_option(G_OPT_V_FIELD);
115 fieldt_opt->key = "to_layer";
116 fieldt_opt->description = _("To layer number or name");
117 fieldt_opt->guisection = _("To");
118
119 typet_opt = G_define_standard_option(G_OPT_V_TYPE);
120 typet_opt->key = "to_type";
121 typet_opt->options = "point,line,boundary";
122 typet_opt->answer = "point";
123 typet_opt->description = _("To feature type");
124 typet_opt->guisection = _("To");
125
126 catt_opt = G_define_standard_option(G_OPT_V_CATS);
127 catt_opt->key = "to_cats";
128 catt_opt->label = _("To category values");
129 catt_opt->guisection = _("To");
130
131 wheret_opt = G_define_standard_option(G_OPT_DB_WHERE);
132 wheret_opt->key = "to_where";
133 wheret_opt->label =
134 _("To WHERE conditions of SQL statement without 'where' keyword");
135 wheret_opt->guisection = _("To");
136
137 afcol = G_define_standard_option(G_OPT_DB_COLUMN);
138 afcol->key = "arc_column";
139 afcol->required = NO;
140 afcol->description =
141 _("Arc forward/both direction(s) cost column (number)");
142 afcol->guisection = _("Cost");
143
144 abcol = G_define_standard_option(G_OPT_DB_COLUMN);
145 abcol->key = "arc_backward_column";
146 abcol->required = NO;
147 abcol->description = _("Arc backward direction cost column (number)");
148 abcol->guisection = _("Cost");
149
150 ncol = G_define_standard_option(G_OPT_DB_COLUMN);
151 ncol->key = "node_column";
152 ncol->required = NO;
153 ncol->description = _("Node cost column (number)");
154 ncol->guisection = _("Cost");
155
156 geo_f = G_define_flag();
157 geo_f->key = 'g';
158 geo_f->description =
159 _("Use geodesic calculation for longitude-latitude locations");
160
161 segments_f = G_define_flag();
162 #if 0
163 /* use this to sync with v.net.path */
164 segments_f->key = 's';
165 segments_f->description = _("Write output as original input segments, "
166 "not each path as one line.");
167 #else
168 segments_f->key = 'l';
169 segments_f->description = _("Write each output path as one line, "
170 "not as original input segments.");
171 #endif
172
173 /* options and flags parser */
174 if (G_parser(argc, argv))
175 exit(EXIT_FAILURE);
176
177 atype = Vect_option_to_types(atype_opt);
178 ttype = Vect_option_to_types(typet_opt);
179
180 Points = Vect_new_line_struct();
181 PPoints = Vect_new_line_struct();
182 Cats = Vect_new_cats_struct();
183 TCats = Vect_new_cats_struct();
184 slist = G_new_ilist();
185
186 Vect_check_input_output_name(map_in->answer, map_out->answer,
187 G_FATAL_EXIT);
188
189 Vect_set_open_level(2);
190
191 if (1 > Vect_open_old(&In, map_in->answer, ""))
192 G_fatal_error(_("Unable to open vector map <%s>"), map_in->answer);
193
194 with_z = Vect_is_3d(&In);
195
196 if (0 > Vect_open_new(&Out, map_out->answer, with_z)) {
197 Vect_close(&In);
198 G_fatal_error(_("Unable to create vector map <%s>"), map_out->answer);
199 }
200
201
202 if (geo_f->answer) {
203 geo = 1;
204 if (G_projection() != PROJECTION_LL)
205 G_warning(_("The current projection is not longitude-latitude"));
206 }
207 else
208 geo = 0;
209
210 #if 0
211 /* use this to sync with v.net.path */
212 segments = segments_f->answer;
213 #else
214 segments = !segments_f->answer;
215 #endif
216
217 nnodes = Vect_get_num_nodes(&In);
218 nlines = Vect_get_num_lines(&In);
219
220 dst = (int *)G_calloc(nnodes + 1, sizeof(int));
221 nxt = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
222 nodes_to_features = (int *)G_calloc(nnodes + 1, sizeof(int));
223 on_path =
224 (struct line_cats **)G_calloc(nlines + 1, sizeof(struct line_cats *));
225 segdir = (char *)G_calloc(nlines + 1, sizeof(char));
226
227 if (!dst || !nxt || !nodes_to_features || !on_path || !segdir)
228 G_fatal_error(_("Out of memory"));
229
230 for (i = 1; i <= nlines; i++) {
231 on_path[i] = Vect_new_cats_struct();
232 segdir[i] = 0;
233 }
234
235 /*initialise varrays and nodes list appropriatelly */
236 afield = Vect_get_field_number(&In, afield_opt->answer);
237 nfield = Vect_get_field_number(&In, nfield_opt->answer);
238
239 flayer = atoi(fieldf_opt->answer);
240 tlayer = atoi(fieldt_opt->answer);
241
242 if (NetA_initialise_varray(&In, flayer, GV_POINT, wheref_opt->answer,
243 catf_opt->answer, &varrayf) <= 0) {
244 G_fatal_error(_("No 'from' features selected. "
245 "Please check options '%s', '%s', '%s'."),
246 fieldf_opt->key, wheref_opt->key, catf_opt->key);
247 }
248
249 if (NetA_initialise_varray(&In, tlayer, ttype, wheret_opt->answer,
250 catt_opt->answer, &varrayt) <= 0) {
251 G_fatal_error(_("No 'to' features selected. "
252 "Please check options '%s', '%s', '%s'."),
253 fieldt_opt->key, wheret_opt->key, catt_opt->key);
254 }
255
256 nodest = Vect_new_list();
257 NetA_varray_to_nodes(&In, varrayt, nodest, nodes_to_features);
258
259 if (nodest->n_values == 0)
260 G_fatal_error(_("No 'to' features"));
261
262 if (0 != Vect_net_build_graph(&In, atype, afield, nfield, afcol->answer, abcol->answer,
263 ncol->answer, geo, 2))
264 G_fatal_error(_("Unable to build graph for vector map <%s>"), Vect_get_full_name(&In));
265
266 graph = Vect_net_get_graph(&In);
267
268 G_message(_("Distances to 'to' features ..."));
269
270 NetA_distance_to_points(graph, nodest, dst, nxt);
271
272 /* Create table */
273 Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
274 Vect_map_add_dblink(&Out, 1, NULL, Fi->table, GV_KEY_COLUMN, Fi->database,
275 Fi->driver);
276 db_init_string(&sql);
277 driver = db_start_driver_open_database(Fi->driver, Fi->database);
278 if (driver == NULL)
279 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
280 Fi->database, Fi->driver);
281 db_set_error_handler_driver(driver);
282
283 sprintf(buf,
284 "create table %s ( cat integer, tcat integer, dist double precision)",
285 Fi->table);
286
287 db_set_string(&sql, buf);
288 G_debug(2, "%s", db_get_string(&sql));
289
290 if (db_execute_immediate(driver, &sql) != DB_OK) {
291 G_fatal_error(_("Unable to create table: '%s'"), db_get_string(&sql));
292 }
293
294 if (db_create_index2(driver, Fi->table, GV_KEY_COLUMN) != DB_OK)
295 G_warning(_("Cannot create index"));
296
297 if (db_grant_on_table
298 (driver, Fi->table, DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK)
299 G_fatal_error(_("Cannot grant privileges on table <%s>"), Fi->table);
300
301 db_begin_transaction(driver);
302
303 Vect_copy_head_data(&In, &Out);
304 Vect_hist_copy(&In, &Out);
305 Vect_hist_command(&Out);
306
307 G_message(_("Tracing paths from 'from' features ..."));
308 from_nr = 0;
309 for (i = 1; i <= nlines; i++) {
310 if (varrayf->c[i]) {
311 int type = Vect_read_line(&In, Points, Cats, i);
312 int node, tcat, cat;
313 double cost;
314 dglInt32_t *vertex, vertex_id;
315
316 if (!Vect_cat_get(Cats, flayer, &cat))
317 continue;
318
319 if (type & GV_POINTS) {
320 node = Vect_find_node(&In, Points->x[0], Points->y[0], Points->z[0], 0, 0);
321 }
322 else {
323 Vect_get_line_nodes(&In, i, &node, NULL);
324 }
325 if (node < 1)
326 continue;
327 if (dst[node] < 0) {
328 /* unreachable */
329 from_nr++;
330 continue;
331 }
332 cost = dst[node] / (double)In.dgraph.cost_multip;
333 vertex = dglGetNode(graph, node);
334 vertex_id = node;
335 slist->n_values = 0;
336 while (nxt[vertex_id] != NULL) {
337 int edge_id;
338
339 edge_id = (int) dglEdgeGet_Id(graph, nxt[vertex_id]);
340 if (segments) {
341 Vect_cat_set(on_path[abs(edge_id)], 1, cat);
342 if (edge_id < 0) {
343 segdir[abs(edge_id)] = 1;
344 }
345 }
346 else
347 G_ilist_add(slist, edge_id);
348
349 vertex = dglEdgeGet_Tail(graph, nxt[vertex_id]);
350 vertex_id = dglNodeGet_Id(graph, vertex);
351 }
352 G_debug(3, "read line %d, vertex id %d", nodes_to_features[vertex_id], (int)vertex_id);
353 Vect_read_line(&In, NULL, TCats, nodes_to_features[vertex_id]);
354 if (!Vect_cat_get(TCats, tlayer, &tcat))
355 continue;
356
357 Vect_write_line(&Out, type, Points, Cats);
358 sprintf(buf, "insert into %s values (%d, %d, %f)", Fi->table, cat,
359 tcat, cost);
360 db_set_string(&sql, buf);
361 G_debug(3, "%s", db_get_string(&sql));
362 if (db_execute_immediate(driver, &sql) != DB_OK) {
363 G_fatal_error(_("Cannot insert new record: %s"),
364 db_get_string(&sql));
365 };
366
367 if (!segments) {
368 Vect_reset_line(PPoints);
369 for (j = 0; j < slist->n_values; j++) {
370 Vect_read_line(&In, Points, NULL, abs(slist->value[j]));
371 if (slist->value[j] > 0)
372 Vect_append_points(PPoints, Points,
373 GV_FORWARD);
374 else
375 Vect_append_points(PPoints, Points,
376 GV_BACKWARD);
377 PPoints->n_points--;
378 }
379 PPoints->n_points++;
380 Vect_reset_cats(Cats);
381 Vect_cat_set(Cats, 1, cat);
382 Vect_write_line(&Out, GV_LINE, PPoints, Cats);
383 }
384
385 }
386 }
387
388 if (segments) {
389 for (i = 1; i <= nlines; i++) {
390 if (on_path[i]->n_cats > 0) {
391 int type;
392
393 if (segdir[i]) {
394 type = Vect_read_line(&In, PPoints, NULL, i);
395 Vect_reset_line(Points);
396 Vect_append_points(Points, PPoints, GV_BACKWARD);
397 }
398 else
399 type = Vect_read_line(&In, Points, NULL, i);
400
401 Vect_write_line(&Out, type, Points, on_path[i]);
402 }
403 }
404 }
405
406 db_commit_transaction(driver);
407 db_close_database_shutdown_driver(driver);
408
409 Vect_build(&Out);
410
411 Vect_close(&In);
412 Vect_close(&Out);
413
414 for (i = 1; i <= nlines; i++)
415 Vect_destroy_cats_struct(on_path[i]);
416 G_free(on_path);
417 G_free(nodes_to_features);
418 G_free(dst);
419 G_free(nxt);
420 G_free(segdir);
421
422 if (from_nr)
423 G_warning(n_("%d 'from' feature was not reachable",
424 "%d 'from' features were not reachable",
425 from_nr), from_nr);
426
427 exit(EXIT_SUCCESS);
428 }
429