1
2 /****************************************************************
3 *
4 * MODULE: v.profile
5 *
6 * AUTHOR(S): Maris Nartiss <maris.gis@gmail.com>
7 * with hints from v.out.ascii, v.buffer, v.what
8 * and other GRASS GIS modules
9 *
10 * PURPOSE: Output vector point/line values along sampling line
11 *
12 * COPYRIGHT: (C) 2008, 2017 by 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 * TODO: Attach a centroid to buffer with tolerance value;
20 * Ability to have "interrupted" profiling line - with holes,
21 * that are not counted into whole profile length;
22 * Implement area sampling by printing out boundary crossing place?
23 * String quoting is unoptimal:
24 * * What if delimiter equals string quote character?
25 * * Quotes within strings are not escaped
26 * * What if user wants to have different quote symbol or no quotes at all?
27 *
28 ****************************************************************/
29
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <grass/config.h>
33 #include <grass/gis.h>
34 #include <grass/vector.h>
35 #include <grass/dbmi.h>
36 #include <grass/glocale.h>
37
38 #include "local_proto.h"
39
40 #if HAVE_GEOS
41 #include <float.h>
42
43 /* A copy of ring2pts from v.buffer geos.c
44 * It is a terrible approach to copy/pasta functions around,
45 * but all this GEOS stuff is just a temporary solution before native
46 * buffers are fixed. (Will happen at some point after next ten years
47 * or so as there's nothing more permanent than a temporary solution.)
48 * 2017-11-19
49 */
ring2pts(const GEOSGeometry * geom,struct line_pnts * Points)50 static int ring2pts(const GEOSGeometry *geom, struct line_pnts *Points)
51 {
52 int i, ncoords;
53 double x, y, z;
54 const GEOSCoordSequence *seq = NULL;
55
56 G_debug(3, "ring2pts()");
57
58 Vect_reset_line(Points);
59 if (!geom) {
60 G_warning(_("Invalid GEOS geometry!"));
61 return 0;
62 }
63 z = 0.0;
64 ncoords = GEOSGetNumCoordinates(geom);
65 if (!ncoords) {
66 G_warning(_("No coordinates in GEOS geometry (can be ok for negative distance)!"));
67 return 0;
68 }
69 seq = GEOSGeom_getCoordSeq(geom);
70 for (i = 0; i < ncoords; i++) {
71 GEOSCoordSeq_getX(seq, i, &x);
72 GEOSCoordSeq_getY(seq, i, &y);
73 if (x != x || x > DBL_MAX || x < -DBL_MAX)
74 G_fatal_error(_("Invalid x coordinate %f"), x);
75 if (y != y || y > DBL_MAX || y < -DBL_MAX)
76 G_fatal_error(_("Invalid y coordinate %f"), y);
77 Vect_append_point(Points, x, y, z);
78 }
79
80 return 1;
81 }
82
83 /* Helper for converting multipoligons to GRASS poligons */
add_poly(const GEOSGeometry * OGeom,struct line_pnts * Buffer)84 static void add_poly(const GEOSGeometry *OGeom, struct line_pnts *Buffer) {
85 const GEOSGeometry *geom2;
86 static struct line_pnts *gPoints;
87 int i, nrings;
88
89 gPoints = Vect_new_line_struct();
90
91 geom2 = GEOSGetExteriorRing(OGeom);
92 if (!ring2pts(geom2, gPoints)) {
93 G_fatal_error(_("Corrupt GEOS geometry"));
94 }
95
96 Vect_append_points(Buffer, gPoints, GV_FORWARD);
97 Vect_reset_line(gPoints);
98
99 nrings = GEOSGetNumInteriorRings(OGeom);
100
101 for (i = 0; i < nrings; i++) {
102 geom2 = GEOSGetInteriorRingN(OGeom, i);
103 if (!ring2pts(geom2, gPoints)) {
104 G_fatal_error(_("Corrupt GEOS geometry"));
105 }
106 Vect_append_points(Buffer, gPoints, GV_FORWARD);
107 Vect_reset_line(gPoints);
108 }
109 }
110 #endif
111
112 static int compdist(const void *, const void *);
113
114 Result *resultset;
115
main(int argc,char * argv[])116 int main(int argc, char *argv[])
117 {
118 struct Map_info In, Out, Pro;
119 struct line_pnts *Points, *Profil, *Buffer, *Ipoints;
120 struct line_cats *Cats;
121 struct ilist *Catlist;
122 FILE *ascii;
123
124 int i, dp, type, otype, id, nrows, ncats, col, more, open3d,
125 layer, pro_layer, *cats, c, field_index, cat;
126 size_t j, rescount;
127
128 /* GCC */
129 int ncols = 0;
130 double xval, yval, bufsize;
131 const char *mapset, *pro_mapset;
132 char sql[200], *fs;
133
134 struct GModule *module;
135 struct Option *old_map, *new_map, *coords_opt, *buffer_opt, *delim_opt,
136 *dp_opt, *layer_opt, *where_opt, *inline_map, *inline_where,
137 *inline_layer, *type_opt, *file_opt;
138 struct Flag *no_column_flag, *no_z_flag;
139
140 struct field_info *Fi, *Fpro;
141 dbDriver *driver;
142 dbHandle handle;
143 dbCursor cursor;
144 dbTable *table;
145 dbColumn *column;
146 dbString table_name, dbsql, valstr;
147
148 /* initialize GIS environment */
149 G_gisinit(argv[0]);
150
151 /* initialize module */
152 module = G_define_module();
153 G_add_keyword(_("vector"));
154 G_add_keyword(_("profile"));
155 G_add_keyword(_("transect"));
156 module->description = _("Vector map profiling tool");
157
158 /* Map to be profiled */
159 old_map = G_define_standard_option(G_OPT_V_INPUT);
160 old_map->required = YES;
161
162 type_opt = G_define_standard_option(G_OPT_V_TYPE);
163 type_opt->options = "point,line";
164 type_opt->answer = "point,line";
165 type_opt->guisection = _("Selection");
166
167 where_opt = G_define_standard_option(G_OPT_DB_WHERE);
168 where_opt->guisection = _("Selection");
169
170 layer_opt = G_define_standard_option(G_OPT_V_FIELD);
171 layer_opt->answer = "1";
172 layer_opt->description = _("Use features only from specified layer");
173 layer_opt->guisection = _("Selection");
174
175 /* Text output details */
176 file_opt = G_define_option();
177 file_opt->key = "output";
178 file_opt->type = TYPE_STRING;
179 file_opt->required = NO;
180 file_opt->multiple = NO;
181 file_opt->gisprompt = "new_file,file,output";
182 file_opt->answer = "-";
183 file_opt->description = _("Path to output text file or - for stdout");
184 file_opt->guisection = _("Format");
185
186 delim_opt = G_define_standard_option(G_OPT_F_SEP);
187 delim_opt->guisection = _("Format");
188
189 dp_opt = G_define_option();
190 dp_opt->key = "dp";
191 dp_opt->type = TYPE_INTEGER;
192 dp_opt->required = NO;
193 dp_opt->options = "0-32";
194 dp_opt->answer = "2";
195 dp_opt->description = _("Number of significant digits");
196 dp_opt->guisection = _("Format");
197
198 /* Profiling tolerance */
199 buffer_opt = G_define_option();
200 buffer_opt->key = "buffer";
201 buffer_opt->type = TYPE_DOUBLE;
202 buffer_opt->required = YES;
203 buffer_opt->answer = "10";
204 buffer_opt->label = _("Buffer (tolerance) for points in map units");
205 buffer_opt->description = _("How far points can be from sampling line");
206
207 /* Storing tolerance area (buffer) is useful to examine which points match */
208 new_map = G_define_option();
209 new_map->key = "map_output";
210 new_map->type = TYPE_STRING;
211 new_map->key_desc = "name";
212 new_map->required = NO;
213 new_map->multiple = NO;
214 new_map->gisprompt = "new,vector,vector";
215 new_map->label = _("Name for profile line and buffer output map");
216 new_map->description =
217 _("Profile line and buffer around it will be written");
218 new_map->guisection = _("Output");
219
220 no_column_flag = G_define_flag();
221 no_column_flag->key = 'c';
222 no_column_flag->description = _("Do not print column names");
223 no_column_flag->guisection = _("Output");
224
225 no_z_flag = G_define_flag();
226 no_z_flag->key = 'z';
227 no_z_flag->label = _("Do not print 3D vector data (z values)");
228 no_z_flag->description = _("Only affects 3D vectors");
229 no_z_flag->guisection = _("Output");
230
231 /* Two ways of defining profiling line:
232 * - by list of coordinates */
233 coords_opt = G_define_standard_option(G_OPT_M_COORDS);
234 coords_opt->multiple = YES;
235 coords_opt->label = _("Coordinates for profiling line nodes");
236 coords_opt->description = _("Specify profiling line vertexes and nodes");
237 coords_opt->guisection = _("Profiling line");
238
239 /* - or profiling line can be taken form other vector map */
240 inline_map = G_define_option();
241 inline_map->key = "profile_map";
242 inline_map->type = TYPE_STRING;
243 inline_map->key_desc = "name";
244 inline_map->required = NO;
245 inline_map->multiple = NO;
246 inline_map->gisprompt = "old,vector,vector";
247 inline_map->label = _("Profiling line map");
248 inline_map->description = _("Vector map containing profiling line");
249 inline_map->guisection = _("Profiling line");
250
251 inline_where = G_define_option();
252 inline_where->key = "profile_where";
253 inline_where->type = TYPE_STRING;
254 inline_where->key_desc = "sql_query";
255 inline_where->required = NO;
256 inline_where->multiple = NO;
257 inline_where->label = _("WHERE conditions for input profile line map");
258 inline_where->description =
259 _("Use to select only one line from profiling line map");
260 inline_where->guisection = _("Profiling line");
261
262 inline_layer = G_define_option();
263 inline_layer->key = "profile_layer";
264 inline_layer->type = TYPE_INTEGER;
265 inline_layer->required = NO;
266 inline_layer->answer = "1";
267 inline_layer->description = _("Profiling line map layer");
268 inline_layer->guisection = _("Profiling line");
269
270 if (G_parser(argc, argv))
271 exit(EXIT_FAILURE);
272
273 #if HAVE_GEOS
274 #if (GEOS_VERSION_MAJOR < 3 || (GEOS_VERSION_MAJOR >= 3 && GEOS_VERSION_MINOR < 3))
275 G_fatal_error("This module requires GEOS >= 3.3");
276 #endif
277 initGEOS(G_message, G_fatal_error);
278 #else
279 G_fatal_error("GRASS native buffering functions are known to return incorrect results.\n"
280 "Till those errors are fixed, this module requires GRASS to be compiled with GEOS support.");
281 #endif
282
283 /* Start with simple input validation and then move to more complex ones */
284 otype = Vect_option_to_types(type_opt);
285 layer = atoi(layer_opt->answer);
286 pro_layer = atoi(inline_layer->answer);
287 if (layer < 1 || pro_layer < 1)
288 G_fatal_error(_("Layer 0 not supported"));
289
290 /* The precision of the output */
291 if (dp_opt->answer) {
292 if (sscanf(dp_opt->answer, "%d", &dp) != 1)
293 G_fatal_error(_("Failed to interpreter 'dp' parameter as an integer"));
294 }
295
296 /* Get buffer size */
297 bufsize = fabs(atof(buffer_opt->answer));
298 if (!(bufsize > 0))
299 G_fatal_error(_("Tolerance value must be greater than 0"));
300
301 /* If new map name is provided, it has to be useable */
302 if (new_map->answer != NULL)
303 if (Vect_legal_filename(new_map->answer) < 1)
304 G_fatal_error(_("<%s> is not a valid vector map name"),
305 new_map->answer);
306
307 /* inline_where has no use if inline_map has been not provided */
308 if (inline_where->answer != NULL && inline_map->answer == NULL)
309 G_fatal_error(_("No input profile map name provided, but WHERE conditions for it have been set"));
310
311 /* Currently only one profile input method is supported */
312 if (inline_map->answer != NULL && coords_opt->answer != NULL)
313 G_fatal_error(_("Profile input coordinates and vector map are provided. "
314 "Please provide only one of them"));
315
316 if (inline_map->answer == NULL && coords_opt->answer == NULL)
317 G_fatal_error(_("No profile input coordinates nor vector map are provided. "
318 "Please provide one of them"));
319
320 /* Where to put module output */
321 if (file_opt->answer) {
322 if (strcmp(file_opt->answer, "-") == 0) {
323 ascii = stdout;
324 }
325 else {
326 ascii = fopen(file_opt->answer, "w");
327 }
328
329 if (ascii == NULL) {
330 G_fatal_error(_("Unable to open file <%s>"), file_opt->answer);
331 }
332 }
333 else {
334 ascii = stdout;
335 }
336
337 /* Create and initialize struct's where to store points/lines and categories */
338 Points = Vect_new_line_struct();
339 Profil = Vect_new_line_struct();
340 Buffer = Vect_new_line_struct();
341 Ipoints = Vect_new_line_struct();
342 Cats = Vect_new_cats_struct();
343
344 /* Construct profile line from user supplied points */
345 if (inline_map->answer == NULL) {
346 for (i = 0; coords_opt->answers[i] != NULL; i += 2) {
347 xval = atof(coords_opt->answers[i]);
348 yval = atof(coords_opt->answers[i + 1]);
349 Vect_append_point(Profil, xval, yval, 0);
350 }
351
352 /* Line is built from two coordinates */
353 if (i == 2)
354 G_fatal_error(_("At least profile start and end coordinates are required!"));
355 }
356 else {
357 /* Check provided profile map name validity */
358 if ((pro_mapset = G_find_vector2(inline_map->answer, "")) == NULL)
359 G_fatal_error(_("Vector map <%s> not found"), inline_map->answer);
360 }
361
362 if ((mapset = G_find_vector2(old_map->answer, "")) == NULL)
363 G_fatal_error(_("Vector map <%s> not found"), old_map->answer);
364
365 if (Vect_set_open_level(2))
366 G_fatal_error(_("Unable to set predetermined vector open level"));
367
368 /* Open existing vector map for reading */
369 if (Vect_open_old(&In, old_map->answer, mapset) < 1)
370 G_fatal_error(_("Unable to open vector map <%s>"), old_map->answer);
371
372 /* Process input as 3D only if it's required */
373 if (!no_z_flag->answer && Vect_is_3d(&In))
374 open3d = WITH_Z;
375 else
376 open3d = WITHOUT_Z;
377
378 /* the field separator */
379 fs = G_option_to_separator(delim_opt);
380
381 /* Let's get vector layers db connections information */
382 Fi = Vect_get_field(&In, layer);
383 if (!Fi && where_opt->answer != NULL) {
384 Vect_close(&In);
385 G_fatal_error(_("No database connection defined for map <%s> layer %d, "
386 "but WHERE condition is provided"), old_map->answer,
387 layer);
388 }
389
390 /* Get profile line from an existing vector map */
391 if (inline_map->answer != NULL) {
392 /* If we get here, pro_mapset is inicialized */
393 if (1 > Vect_open_old(&Pro, inline_map->answer, pro_mapset))
394 G_fatal_error(_("Unable to open vector map <%s>"),
395 inline_map->answer);
396 if (inline_where->answer != NULL) {
397 Fpro = Vect_get_field(&Pro, pro_layer);
398 if (!Fpro) {
399 Vect_close(&In);
400 Vect_close(&Pro);
401 G_fatal_error(_("No database connection defined for map <%s> layer %d, "
402 "but WHERE condition is provided"),
403 inline_map->answer, pro_layer);
404 }
405 /* Prepeare strings for use in db_* calls */
406 db_init_string(&dbsql);
407 db_init_string(&valstr);
408 db_init_string(&table_name);
409 db_init_handle(&handle);
410 G_debug(1,
411 "Field number:%d; Name:<%s>; Driver:<%s>; Database:<%s>; Table:<%s>; Key:<%s>",
412 Fpro->number, Fpro->name, Fpro->driver, Fpro->database,
413 Fpro->table, Fpro->key);
414
415 /* Prepearing database for use */
416 driver = db_start_driver(Fpro->driver);
417 if (driver == NULL) {
418 Vect_close(&In);
419 Vect_close(&Pro);
420 G_fatal_error(_("Unable to start driver <%s>"), Fpro->driver);
421 }
422 db_set_handle(&handle, Fpro->database, NULL);
423 if (db_open_database(driver, &handle) != DB_OK) {
424 Vect_close(&In);
425 Vect_close(&Pro);
426 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
427 Fpro->database, Fpro->driver);
428 }
429 db_set_string(&table_name, Fpro->table);
430 if (db_describe_table(driver, &table_name, &table) != DB_OK) {
431 Vect_close(&In);
432 Vect_close(&Pro);
433 G_fatal_error(_("Unable to describe table <%s>"),
434 Fpro->table);
435 }
436 ncols = db_get_table_number_of_columns(table);
437
438 ncats =
439 db_select_int(driver, Fpro->table, Fpro->key,
440 inline_where->answer, &cats);
441 if (ncats < 1) {
442 Vect_close(&In);
443 Vect_close(&Pro);
444 G_fatal_error(_("No features match Your query"));
445 }
446 if (ncats > 1) {
447 Vect_close(&In);
448 Vect_close(&Pro);
449 G_fatal_error(_("Your query matches more than one record in input profiling map. "
450 "Currently it's not supported. Enhance WHERE conditions to get only one line."));
451 }
452 if (!(Catlist = Vect_new_list())) {
453 Vect_close(&In);
454 Vect_close(&Pro);
455 G_fatal_error(_("Error while initialising line list"));
456 }
457 /* Get all features matching specified CAT value */
458 Vect_cidx_find_all(&Pro, pro_layer, GV_LINE, cats[0], Catlist);
459 if (Catlist->n_values > 1) {
460 Vect_close(&In);
461 Vect_close(&Pro);
462 G_fatal_error(_("Your query matches more than one record in input profiling map. "
463 "Currently it's not supported. Enhance WHERE conditions to get only one line."));
464 }
465 if (Vect_read_line(&Pro, Profil, NULL, Catlist->value[0]) !=
466 GV_LINE) {
467 G_fatal_error(_("Error while reading vector feature from profile line map"));
468 }
469 }
470 else {
471 /* WHERE not provided -> assuming profiling line map contains only one line */
472 c = 0;
473 while ((type = Vect_read_next_line(&Pro, Profil, NULL)) > 0) {
474 if (type & GV_LINE)
475 c++;
476 }
477 /* Currently we support only one SINGLE input line */
478 if (c > 1) {
479 Vect_close(&In);
480 Vect_close(&Pro);
481 G_fatal_error(_("Your input profile map contains more than one line. "
482 "Currently it's not supported. Provide WHERE conditions to get only one line."));
483 }
484 }
485 }
486
487 /* Create a buffer around profile line for point sampling
488 Tolerance is calculated in such way that buffer will have flat end and no cap. */
489 /* Native buffering is known to fail.
490 Vect_line_buffer(Profil, bufsize, 1 - (bufsize * cos((2 * M_PI) / 2)),
491 Buffer);
492 */
493 #ifdef HAVE_GEOS
494 /* Code lifted from v.buffer geos.c (with modifications) */
495 GEOSGeometry *IGeom;
496 GEOSGeometry *OGeom = NULL;
497 const GEOSGeometry *geom2 = NULL;
498
499 IGeom = Vect_line_to_geos(Profil, GV_LINE, 0);
500 if (!IGeom) {
501 G_fatal_error(_("Failed to convert GRASS line to GEOS line"));
502 }
503
504 GEOSBufferParams* geos_params = GEOSBufferParams_create();
505 GEOSBufferParams_setEndCapStyle(geos_params, GEOSBUF_CAP_FLAT);
506 OGeom = GEOSBufferWithParams(IGeom, geos_params, bufsize);
507 GEOSBufferParams_destroy(geos_params);
508 if (!OGeom) {
509 G_fatal_error(_("Buffering failed"));
510 }
511
512 if (GEOSGeomTypeId(OGeom) == GEOS_MULTIPOLYGON) {
513 int ngeoms = GEOSGetNumGeometries(OGeom);
514 for (i = 0; i < ngeoms; i++) {
515 geom2 = GEOSGetGeometryN(OGeom, i);
516 add_poly(geom2, Buffer);
517 }
518 }
519 else {
520 add_poly(OGeom, Buffer);
521 }
522
523 if (IGeom)
524 GEOSGeom_destroy(IGeom);
525 if (OGeom)
526 GEOSGeom_destroy(OGeom);
527 finishGEOS();
528 #endif
529
530 Vect_cat_set(Cats, 1, 1);
531
532 /* Should we store used buffer for later examination? */
533 if (new_map->answer != NULL) {
534 /* Open new vector for reading/writing */
535 if (0 > Vect_open_new(&Out, new_map->answer, WITHOUT_Z)) {
536 Vect_close(&In);
537 G_fatal_error(_("Unable to create vector map <%s>"),
538 new_map->answer);
539 }
540
541 /* Write profile line and it's buffer into new vector map */
542 Vect_write_line(&Out, GV_LINE, Profil, Cats);
543 /* No category for boundary */
544 Vect_reset_cats(Cats);
545 Vect_write_line(&Out, GV_BOUNDARY, Buffer, Cats);
546 }
547
548 /* Here starts processing of input map */
549 rescount = 0;
550 resultset = NULL;
551
552 /* If input vector has a database connection... */
553 if (Fi != NULL) {
554 field_index = Vect_cidx_get_field_index(&In, layer);
555 if (field_index < 0) {
556 G_fatal_error(_("Vector map <%s> does not have cat's defined on layer %d"),
557 old_map->answer, layer);
558 }
559
560 /* Prepeare strings for use in db_* calls */
561 db_init_string(&dbsql);
562 db_init_string(&valstr);
563 db_init_string(&table_name);
564 db_init_handle(&handle);
565
566 G_debug(1,
567 "Field number:%d; Name:<%s>; Driver:<%s>; Database:<%s>; Table:<%s>; Key:<%s>",
568 Fi->number, Fi->name, Fi->driver, Fi->database, Fi->table,
569 Fi->key);
570
571 /* Prepearing database for use */
572 driver = db_start_driver(Fi->driver);
573 if (driver == NULL) {
574 Vect_close(&In);
575 G_fatal_error(_("Unable to start driver <%s>"), Fi->driver);
576 }
577 db_set_handle(&handle, Fi->database, NULL);
578 if (db_open_database(driver, &handle) != DB_OK) {
579 Vect_close(&In);
580 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
581 Fi->database, Fi->driver);
582 }
583 db_set_string(&table_name, Fi->table);
584 if (db_describe_table(driver, &table_name, &table) != DB_OK) {
585 Vect_close(&In);
586 G_fatal_error(_("Unable to describe table <%s>"), Fi->table);
587 }
588 ncols = db_get_table_number_of_columns(table);
589
590 /* Create vector feature list for future processing by applying SQL WHERE conditions... */
591 if (where_opt->answer != NULL) {
592 ncats =
593 db_select_int(driver, Fi->table, Fi->key, where_opt->answer,
594 &cats);
595 if (ncats < 1)
596 G_fatal_error(_("No features match Your query"));
597 for (i = 0; i < ncats; i++) {
598 c = Vect_cidx_find_next(&In, field_index, cats[i],
599 otype, 0, &type, &id);
600 /* Crunch over all points/lines, that match specified CAT */
601 while (c >= 0) {
602 c++;
603 if (type & otype) {
604 switch (Vect_read_line(&In, Points, Cats, id)) {
605 case GV_POINT:
606 Vect_cat_get(Cats, layer, &cat);
607 proc_point(Points, Profil, Buffer, cat,
608 &rescount, open3d);
609 break;
610 case GV_LINE:
611 Vect_reset_line(Ipoints);
612 if (Vect_line_get_intersections
613 (Profil, Points, Ipoints, open3d) > 0) {
614 Vect_cat_get(Cats, layer, &cat);
615 proc_line(Ipoints, Profil, cat, &rescount,
616 open3d);
617 }
618 break;
619 }
620 }
621 else
622 G_fatal_error
623 ("Error in Vect_cidx_find_next function! Report a bug.");
624 c = Vect_cidx_find_next(&In, field_index, cats[i], otype,
625 c, &type, &id);
626 }
627 }
628 }
629 }
630
631 /* Process all lines IF no database exists or WHERE was not provided.
632 Read in single line and get its type */
633 if (Fi == NULL || (where_opt->answer == NULL && Fi != NULL)) {
634 while ((type = Vect_read_next_line(&In, Points, Cats)) > 0) {
635 if (type & GV_POINT) {
636 Vect_cat_get(Cats, layer, &cat);
637 proc_point(Points, Profil, Buffer, cat, &rescount, open3d);
638 }
639 if (type & GV_LINE) {
640 Vect_reset_line(Ipoints);
641 if (Vect_line_get_intersections
642 (Profil, Points, Ipoints, open3d) > 0)
643 Vect_cat_get(Cats, layer, &cat);
644 proc_line(Ipoints, Profil, cat, &rescount, open3d);
645 }
646 }
647 }
648 /* We don't need input vector anymore */
649 Vect_close(&In);
650 G_debug(1, "There are %zu features matching profile line", rescount);
651
652 /* Sort results by distance, cat */
653 qsort(&resultset[0], rescount, sizeof(Result), compdist);
654
655 /* Print out column names */
656 if (!no_column_flag->answer) {
657 fprintf(ascii, "Number%sDistance", fs);
658 if (open3d == WITH_Z)
659 fprintf(ascii, "%sZ", fs);
660 if (Fi != NULL) {
661 /* ncols are initialized here from previous Fi != NULL if block */
662 for (col = 0; col < ncols; col++) {
663 column = db_get_table_column(table, col);
664 fprintf(ascii, "%s%s", fs, db_get_column_name(column));
665 }
666 }
667 fprintf(ascii, "\n");
668 }
669
670 /* Print out result */
671 for (j = 0; j < rescount; j++) {
672 fprintf(ascii, "%zu%s%.*f", j + 1, fs, dp, resultset[j].distance);
673 if (open3d == WITH_Z)
674 fprintf(ascii, "%s%.*f", fs, dp, resultset[j].z);
675 if (Fi != NULL) {
676 sprintf(sql, "select * from %s where %s=%d", Fi->table, Fi->key,
677 resultset[j].cat);
678 G_debug(2, "SQL: \"%s\"", sql);
679 db_set_string(&dbsql, sql);
680 /* driver IS initialized here in case if Fi != NULL */
681 if (db_open_select_cursor(driver, &dbsql, &cursor, DB_SEQUENTIAL)
682 != DB_OK)
683 G_warning(_("Unabale to get attribute data for cat %d"),
684 resultset[j].cat);
685 else {
686 nrows = db_get_num_rows(&cursor);
687 G_debug(1, "Result count: %d", nrows);
688 table = db_get_cursor_table(&cursor);
689
690 if (nrows > 0) {
691 if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK) {
692 G_warning(_("Error while retreiving database record for cat %d"),
693 resultset[j].cat);
694 }
695 else {
696 for (col = 0; col < ncols; col++) {
697 /* Column description retreiving is fast, as they live in provided table structure */
698 column = db_get_table_column(table, col);
699 db_convert_column_value_to_string(column,
700 &valstr);
701 type = db_get_column_sqltype(column);
702
703 /* Those values should be quoted */
704 if (type == DB_SQL_TYPE_CHARACTER ||
705 type == DB_SQL_TYPE_DATE ||
706 type == DB_SQL_TYPE_TIME ||
707 type == DB_SQL_TYPE_TIMESTAMP ||
708 type == DB_SQL_TYPE_INTERVAL ||
709 type == DB_SQL_TYPE_TEXT ||
710 type == DB_SQL_TYPE_SERIAL)
711 fprintf(ascii, "%s\"%s\"", fs,
712 db_get_string(&valstr));
713 else
714 fprintf(ascii, "%s%s", fs,
715 db_get_string(&valstr));
716 }
717 }
718 }
719 else {
720 for (col = 0; col < ncols; col++) {
721 fprintf(ascii, "%s", fs);
722 }
723 }
724 }
725 db_close_cursor(&cursor);
726 }
727 /* Terminate attribute data line and flush data to provided output (file/stdout) */
728 fprintf(ascii, "\n");
729 if (fflush(ascii))
730 G_fatal_error(_("Can not write data portion to provided output"));
731 }
732
733 /* Build topology for vector map and close them */
734 if (new_map->answer != NULL) {
735 Vect_build(&Out);
736 Vect_close(&Out);
737 }
738
739 if (Fi != NULL)
740 db_close_database_shutdown_driver(driver);
741
742 if (ascii != NULL)
743 fclose(ascii);
744
745 /* Don't forget to report to caller successful end of data processing :) */
746 exit(EXIT_SUCCESS);
747 }
748
749 /* Qsort distance comparison function */
compdist(const void * d1,const void * d2)750 static int compdist(const void *d1, const void *d2)
751 {
752 Result *r1 = (Result *) d1, *r2 = (Result *) d2;
753
754 G_debug(5, "Comparing %f with %f", r1->distance, r2->distance);
755
756 if (r1->distance == r2->distance) {
757 if (r1->cat < r2->cat)
758 return -1;
759 else
760 return (r1->cat > r2->cat);
761 }
762 if (r1->distance < r2->distance)
763 return -1;
764 else
765 return (r1->distance > r2->distance);
766 }
767