1 /*!
2 \file lib/vector/Vlib/read_ogr.c
3
4 \brief Vector library - reading data (OGR format)
5
6 Higher level functions for reading/writing/manipulating vectors.
7
8 (C) 2001-2011 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, Piero Cavalieri
14 \author Martin Landa <landa.martin gmail.com>
15 */
16
17 #include <grass/vector.h>
18 #include <grass/glocale.h>
19
20 #ifdef HAVE_OGR
21 #include <ogr_api.h>
22
23 static int cache_feature(struct Map_info *, OGRGeometryH, int);
24 static int read_line(const struct Map_info *, OGRGeometryH, long,
25 struct line_pnts *);
26 static int get_line_type(const struct Map_info *, long);
27 static int read_next_line_ogr(struct Map_info *, struct line_pnts *,
28 struct line_cats *, int);
29 #endif
30
31 /*!
32 \brief Read next feature from OGR layer.
33 Skip empty features (level 1 without topology).
34
35 This function implements sequential access.
36
37 The action of this routine can be modified by:
38 - Vect_read_constraint_region()
39 - Vect_read_constraint_type()
40 - Vect_remove_constraints()
41
42 \param Map pointer to Map_info structure
43 \param[out] line_p container used to store line points within
44 \param[out] line_c container used to store line categories within
45
46 \return feature type
47 \return -2 no more features (EOF)
48 \return -1 out of memory
49 */
V1_read_next_line_ogr(struct Map_info * Map,struct line_pnts * line_p,struct line_cats * line_c)50 int V1_read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
51 struct line_cats *line_c)
52 {
53 #ifdef HAVE_OGR
54 return read_next_line_ogr(Map, line_p, line_c, FALSE);
55 #else
56 G_fatal_error(_("GRASS is not compiled with OGR support"));
57 return -1;
58 #endif
59 }
60
61 /*!
62 \brief Read next feature from OGR layer on topological level.
63
64 This function implements sequential access.
65
66 \param Map pointer to Map_info structure
67 \param[out] line_p container used to store line points within
68 (pointer to line_pnts struct)
69 \param[out] line_c container used to store line categories within
70 (pointer to line_cats struct)
71
72 \return feature type
73 \return -2 no more features (EOF)
74 \return -1 on failure
75 */
V2_read_next_line_ogr(struct Map_info * Map,struct line_pnts * line_p,struct line_cats * line_c)76 int V2_read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
77 struct line_cats *line_c)
78 {
79 #ifdef HAVE_OGR
80 int line, ret;
81 struct P_line *Line;
82 struct bound_box lbox, mbox;
83
84 G_debug(3, "V2_read_next_line_ogr()");
85
86 if (Map->constraint.region_flag)
87 Vect_get_constraint_box(Map, &mbox);
88
89 while(TRUE) {
90 line = Map->next_line;
91
92 if (Map->next_line > Map->plus.n_lines)
93 return -2; /* nothing to read */
94
95 Map->next_line++;
96 Line = Map->plus.Line[line];
97 if (Line == NULL) { /* skip dead features */
98 continue;
99 }
100
101 if (Map->constraint.type_flag) {
102 /* skip feature by type */
103 if (!(Line->type & Map->constraint.type))
104 continue;
105 }
106
107 if (Line->type == GV_CENTROID) {
108 G_debug(4, "Centroid");
109
110 if (line_p != NULL) {
111 int i, found;
112 struct bound_box box;
113 struct boxlist list;
114 struct P_topo_c *topo = (struct P_topo_c *)Line->topo;
115
116 /* get area bbox */
117 Vect_get_area_box(Map, topo->area, &box);
118 /* search in spatial index for centroid with area bbox */
119 dig_init_boxlist(&list, TRUE);
120 Vect_select_lines_by_box(Map, &box, Line->type, &list);
121
122 found = 0;
123 for (i = 0; i < list.n_values; i++) {
124 if (list.id[i] == line) {
125 found = i;
126 break;
127 }
128 }
129
130 Vect_reset_line(line_p);
131 Vect_append_point(line_p, list.box[found].E, list.box[found].N, 0.0);
132 }
133 if (line_c != NULL) {
134 /* cat = FID and offset = FID for centroid */
135 Vect_reset_cats(line_c);
136 Vect_cat_set(line_c, 1, (int) Line->offset);
137 }
138
139 ret = GV_CENTROID;
140 }
141 else {
142 ret = read_next_line_ogr(Map, line_p, line_c, TRUE);
143 }
144
145 if (line_p && Map->constraint.region_flag) {
146 /* skip feature by region */
147 Vect_line_box(line_p, &lbox);
148 if (!Vect_box_overlap(&lbox, &mbox))
149 continue;
150 }
151
152 /* skip feature by field ignored */
153
154 return ret;
155 }
156 #else
157 G_fatal_error(_("GRASS is not compiled with OGR support"));
158 return -1;
159 #endif
160 }
161
162 /*!
163 \brief Read feature from OGR layer at given offset (level 1 without topology)
164
165 This function implements random access on level 1.
166
167 \param Map pointer to Map_info structure
168 \param[out] line_p container used to store line points within
169 (pointer line_pnts struct)
170 \param[out] line_c container used to store line categories within
171 (pointer line_cats struct)
172 \param offset given offset
173
174 \return line type
175 \return 0 dead line
176 \return -2 no more features
177 \return -1 on failure
178 */
V1_read_line_ogr(struct Map_info * Map,struct line_pnts * line_p,struct line_cats * line_c,off_t offset)179 int V1_read_line_ogr(struct Map_info *Map,
180 struct line_pnts *line_p, struct line_cats *line_c, off_t offset)
181 {
182 #ifdef HAVE_OGR
183 long fid;
184 int type;
185 OGRGeometryH hGeom;
186
187 struct Format_info_ogr *ogr_info;
188
189 ogr_info = &(Map->fInfo.ogr);
190 G_debug(3, "V1_read_line_ogr(): offset = %lu offset_num = %lu",
191 (long) offset, (long) ogr_info->offset.array_num);
192
193 if (offset >= ogr_info->offset.array_num)
194 return -2; /* nothing to read */
195
196 if (line_p != NULL)
197 Vect_reset_line(line_p);
198 if (line_c != NULL)
199 Vect_reset_cats(line_c);
200
201 fid = ogr_info->offset.array[offset];
202 G_debug(4, " fid = %ld", fid);
203
204 /* coordinates */
205 if (line_p != NULL) {
206 /* read feature to cache if necessary */
207 if (ogr_info->cache.fid != fid) {
208 G_debug(4, "Read feature (fid = %ld) to cache", fid);
209 if (ogr_info->feature_cache) {
210 OGR_F_Destroy(ogr_info->feature_cache);
211 }
212 ogr_info->feature_cache =
213 OGR_L_GetFeature(ogr_info->layer, fid);
214 if (ogr_info->feature_cache == NULL) {
215 G_warning(_("Unable to get feature geometry, fid %ld"),
216 fid);
217 return -1;
218 }
219 ogr_info->cache.fid = fid;
220 }
221
222 hGeom = OGR_F_GetGeometryRef(ogr_info->feature_cache);
223 if (hGeom == NULL) {
224 G_warning(_("Unable to get feature geometry, fid %ld"),
225 fid);
226 return -1;
227 }
228
229 type = read_line(Map, hGeom, offset + 1, line_p);
230 }
231 else {
232 type = get_line_type(Map, fid);
233 }
234
235 /* category */
236 if (line_c != NULL) {
237 Vect_cat_set(line_c, 1, (int) fid);
238 }
239
240 return type;
241 #else
242 G_fatal_error(_("GRASS is not compiled with OGR support"));
243 return -1;
244 #endif
245 }
246
247 #ifdef HAVE_OGR
248 /*!
249 \brief Recursively read feature and add all elements to points_cache and types_cache.
250
251 ftype: if > 0 use this type (because parts of Polygon are read as wkbLineString)
252
253 \param Map pointer to Map_info structure
254 \param[out] hGeom OGR geometry
255 \param ftype feature type
256
257 \return 0 on success
258 \return 1 on error
259 */
cache_feature(struct Map_info * Map,OGRGeometryH hGeom,int ftype)260 int cache_feature(struct Map_info *Map, OGRGeometryH hGeom, int ftype)
261 {
262 int line, i, np, ng, tp;
263
264 struct Format_info_ogr *ogr_info;
265
266 OGRwkbGeometryType type;
267 OGRGeometryH hGeom2;
268
269 G_debug(4, "cache_feature() ftype = %d", ftype);
270
271 ogr_info = &(Map->fInfo.ogr);
272
273 /* alloc space in lines cache */
274 line = ogr_info->cache.lines_num;
275 if (line == ogr_info->cache.lines_alloc) {
276 ogr_info->cache.lines_alloc += 1;
277 ogr_info->cache.lines =
278 (struct line_pnts **)G_realloc((void *)ogr_info->cache.lines,
279 ogr_info->cache.lines_alloc *
280 sizeof(struct line_pnts *));
281
282 ogr_info->cache.lines_types =
283 (int *)G_realloc(ogr_info->cache.lines_types,
284 ogr_info->cache.lines_alloc * sizeof(int));
285
286 for (i = ogr_info->cache.lines_num; i < ogr_info->cache.lines_alloc; i++)
287 ogr_info->cache.lines[i] = Vect_new_line_struct();
288 }
289 Vect_reset_line(ogr_info->cache.lines[line]);
290
291 type = wkbFlatten(OGR_G_GetGeometryType(hGeom));
292
293 switch (type) {
294 case wkbPoint:
295 G_debug(4, "Point");
296 Vect_append_point(ogr_info->cache.lines[line],
297 OGR_G_GetX(hGeom, 0), OGR_G_GetY(hGeom, 0),
298 OGR_G_GetZ(hGeom, 0));
299 ogr_info->cache.lines_types[line] = GV_POINT;
300 ogr_info->cache.lines_num++;
301 return 0;
302 break;
303
304 case wkbLineString:
305 G_debug(4, "LineString");
306 np = OGR_G_GetPointCount(hGeom);
307 for (i = 0; i < np; i++) {
308 Vect_append_point(ogr_info->cache.lines[line],
309 OGR_G_GetX(hGeom, i), OGR_G_GetY(hGeom, i),
310 OGR_G_GetZ(hGeom, i));
311 }
312
313 if (ftype > 0) { /* Polygon rings */
314 ogr_info->cache.lines_types[line] = ftype;
315 }
316 else {
317 ogr_info->cache.lines_types[line] = GV_LINE;
318 }
319 ogr_info->cache.lines_num++;
320 return 0;
321 break;
322
323 case wkbMultiPoint:
324 case wkbMultiLineString:
325 case wkbPolygon:
326 case wkbMultiPolygon:
327 case wkbGeometryCollection:
328 ng = OGR_G_GetGeometryCount(hGeom);
329 G_debug(4, "%d geoms -> next level", ng);
330 if (type == wkbPolygon) {
331 tp = GV_BOUNDARY;
332 }
333 else {
334 tp = -1;
335 }
336 for (i = 0; i < ng; i++) {
337 hGeom2 = OGR_G_GetGeometryRef(hGeom, i);
338 cache_feature(Map, hGeom2, tp);
339 }
340 return 0;
341 break;
342
343 default:
344 G_warning(_("OGR feature type %d not supported"), type);
345 return 1;
346 break;
347 }
348 }
349
read_next_line_ogr(struct Map_info * Map,struct line_pnts * line_p,struct line_cats * line_c,int ignore_constraint)350 int read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
351 struct line_cats *line_c, int ignore_constraint)
352 {
353 int itype;
354 struct bound_box lbox, mbox;
355 OGRFeatureH hFeature;
356 OGRGeometryH hGeom;
357
358 struct Format_info_ogr *ogr_info;
359
360 G_debug(3, "V1_read_next_line_ogr()");
361
362 if (Map->constraint.region_flag && !ignore_constraint)
363 Vect_get_constraint_box(Map, &mbox);
364
365 ogr_info = &(Map->fInfo.ogr);
366 while (TRUE) {
367 /* reset data structures */
368 if (line_p != NULL)
369 Vect_reset_line(line_p);
370 if (line_c != NULL)
371 Vect_reset_cats(line_c);
372
373 /* read feature to cache if necessary */
374 while (ogr_info->cache.lines_next == ogr_info->cache.lines_num) {
375 hFeature = OGR_L_GetNextFeature(ogr_info->layer);
376 if (hFeature == NULL) {
377 return -2; /* nothing to read */
378 }
379
380 hGeom = OGR_F_GetGeometryRef(hFeature);
381 if (hGeom == NULL) { /* skip feature without geometry */
382 G_warning(_("Feature without geometry. Skipped."));
383 OGR_F_Destroy(hFeature);
384 continue;
385 }
386
387 /* cache OGR feature */
388 ogr_info->cache.fid = (int)OGR_F_GetFID(hFeature);
389 if (ogr_info->cache.fid == OGRNullFID) {
390 G_warning(_("OGR feature without ID"));
391 }
392
393 /* cache feature */
394 ogr_info->cache.lines_num = 0;
395 cache_feature(Map, hGeom, -1);
396 G_debug(4, "%d lines read to cache", ogr_info->cache.lines_num);
397 OGR_F_Destroy(hFeature);
398
399 /* next to be read from cache */
400 ogr_info->cache.lines_next = 0;
401 }
402
403 /* read next part of the feature */
404 G_debug(4, "read next cached line %d", ogr_info->cache.lines_next);
405 itype = ogr_info->cache.lines_types[ogr_info->cache.lines_next];
406
407 if (Map->constraint.type_flag && !ignore_constraint) {
408 /* skip feature by type */
409 if (!(itype & Map->constraint.type)) {
410 ogr_info->cache.lines_next++;
411 continue;
412 }
413 }
414
415 if (Map->constraint.region_flag && !ignore_constraint) {
416 /* skip feature by region */
417 Vect_line_box(ogr_info->cache.lines[ogr_info->cache.lines_next],
418 &lbox);
419
420 if (!Vect_box_overlap(&lbox, &mbox)) {
421 ogr_info->cache.lines_next++;
422 continue;
423 }
424 }
425
426 /* skip feature by field - ignored */
427
428 if (line_p != NULL)
429 Vect_append_points(line_p,
430 ogr_info->cache.lines[ogr_info->cache.lines_next], GV_FORWARD);
431
432 if (line_c != NULL && ogr_info->cache.fid != OGRNullFID)
433 Vect_cat_set(line_c, 1, ogr_info->cache.fid);
434
435 ogr_info->cache.lines_next++;
436 G_debug(4, "next line read, type = %d", itype);
437
438 return itype;
439 }
440
441 return -1; /* not reached */
442 }
443
444 /*!
445 \brief Recursively descend to feature and read the part
446
447 \param Map pointer to Map_info structure
448 \param hGeom OGR geometry
449 \param offset given offset
450 \param[out] Points container used to store line pointes within
451
452 \return feature type
453 \return -1 on error
454 */
read_line(const struct Map_info * Map,OGRGeometryH hGeom,long offset,struct line_pnts * Points)455 int read_line(const struct Map_info *Map, OGRGeometryH hGeom, long offset,
456 struct line_pnts *Points)
457 {
458 int i, nPoints;
459 int eType, line;
460
461 const struct Format_info_ogr *ogr_info;
462
463 OGRGeometryH hGeom2;
464
465 /* Read coors if hGeom is a simple element (wkbPoint,
466 * wkbLineString) otherwise descend to geometry specified by
467 * offset[offset] */
468
469 ogr_info = &(Map->fInfo.ogr);
470
471 eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
472 G_debug(4, "OGR geometry type: %d", eType);
473
474 switch (eType) {
475 case wkbPoint:
476 G_debug(4, "\t->Point");
477 if (Points) {
478 Vect_append_point(Points, OGR_G_GetX(hGeom, 0), OGR_G_GetY(hGeom, 0),
479 OGR_G_GetZ(hGeom, 0));
480 }
481 return GV_POINT;
482 break;
483
484 case wkbLineString:
485 G_debug(4, "\t->LineString");
486 if (Points) {
487 nPoints = OGR_G_GetPointCount(hGeom);
488 for (i = 0; i < nPoints; i++) {
489 Vect_append_point(Points, OGR_G_GetX(hGeom, i),
490 OGR_G_GetY(hGeom, i), OGR_G_GetZ(hGeom, i));
491 }
492 }
493 return GV_LINE;
494 break;
495
496 case wkbPolygon:
497 case wkbMultiPoint:
498 case wkbMultiLineString:
499 case wkbMultiPolygon:
500 case wkbGeometryCollection:
501 G_debug(4, "\t->more geoms -> part %d", ogr_info->offset.array[offset]);
502 hGeom2 = OGR_G_GetGeometryRef(hGeom, ogr_info->offset.array[offset]);
503 line = read_line(Map, hGeom2, offset + 1, Points);
504 if (eType == wkbPolygon || eType == wkbMultiPolygon)
505 return GV_BOUNDARY;
506 if (eType == wkbMultiPoint)
507 return GV_POINT;
508 if (eType == wkbMultiLineString)
509 return GV_LINE;
510 return line;
511 break;
512
513 default:
514 G_warning(_("OGR feature type '%s' not supported"),
515 OGRGeometryTypeToName(eType));
516 break;
517 }
518
519 return -1;
520 }
521
522 /*!
523 \brief Recursively descend to feature and read the part
524
525 \param Map pointer to Map_info structure
526 \param hGeom OGR geometry
527 \param offset given offset
528 \param[out] Points container used to store line pointes within
529
530 \return feature type
531 \return -1 on error
532 */
get_line_type(const struct Map_info * Map,long fid)533 int get_line_type(const struct Map_info *Map, long fid)
534 {
535 int eType;
536
537 const struct Format_info_ogr *ogr_info;
538
539 OGRFeatureH hFeat;
540 OGRGeometryH hGeom;
541
542 G_debug(4, "get_line_type() fid = %ld", fid);
543
544 ogr_info = &(Map->fInfo.ogr);
545
546 hFeat = OGR_L_GetFeature(ogr_info->layer, fid);
547 if (hFeat == NULL)
548 return -1;
549
550 hGeom = OGR_F_GetGeometryRef(hFeat);
551 if (hGeom == NULL)
552 return -1;
553
554 eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
555
556 OGR_F_Destroy(hFeat);
557
558 G_debug(4, "OGR Geometry of type: %d", eType);
559
560 switch (eType) {
561 case wkbPoint:
562 case wkbMultiPoint:
563 return GV_POINT;
564 break;
565
566 case wkbLineString:
567 case wkbMultiLineString:
568 return GV_LINE;
569 break;
570
571 case wkbPolygon:
572 case wkbMultiPolygon:
573 case wkbGeometryCollection:
574 return GV_BOUNDARY;
575 break;
576
577 default:
578 G_warning(_("OGR feature type %d not supported"), eType);
579 break;
580 }
581
582 return -1;
583 }
584 #endif
585