GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-7413740dd8
read_ogr.c
Go to the documentation of this file.
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(struct Map_info *, OGRGeometryH, long, struct line_pnts *);
25 static int get_line_type(struct Map_info *, long);
26 static int read_next_line_ogr(struct Map_info *, struct line_pnts *,
27  struct line_cats *, int);
28 #endif
29 
30 /*!
31  \brief Read next feature from OGR layer.
32  Skip empty features (level 1 without topology).
33 
34  This function implements sequential access.
35 
36  The action of this routine can be modified by:
37  - Vect_read_constraint_region()
38  - Vect_read_constraint_type()
39  - Vect_remove_constraints()
40 
41  \param Map pointer to Map_info structure
42  \param[out] line_p container used to store line points within
43  \param[out] line_c container used to store line categories within
44 
45  \return feature type
46  \return -2 no more features (EOF)
47  \return -1 out of memory
48  */
49 int V1_read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
50  struct line_cats *line_c)
51 {
52 #ifdef HAVE_OGR
53  return read_next_line_ogr(Map, line_p, line_c, FALSE);
54 #else
55  G_fatal_error(_("GRASS is not compiled with OGR support"));
56  return -1;
57 #endif
58 }
59 
60 /*!
61  \brief Read next feature from OGR layer on topological level.
62 
63  This function implements sequential access.
64 
65  \param Map pointer to Map_info structure
66  \param[out] line_p container used to store line points within
67  (pointer to line_pnts struct)
68  \param[out] line_c container used to store line categories within
69  (pointer to line_cats struct)
70 
71  \return feature type
72  \return -2 no more features (EOF)
73  \return -1 on failure
74  */
75 int V2_read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
76  struct line_cats *line_c)
77 {
78 #ifdef HAVE_OGR
79  int line, ret;
80  struct P_line *Line;
81  struct bound_box lbox, mbox;
82 
83  G_debug(3, "V2_read_next_line_ogr()");
84 
85  if (Map->constraint.region_flag)
86  Vect_get_constraint_box(Map, &mbox);
87 
88  while (TRUE) {
89  line = Map->next_line;
90 
91  if (Map->next_line > Map->plus.n_lines)
92  return -2; /* nothing to read */
93 
94  Map->next_line++;
95  Line = Map->plus.Line[line];
96  if (Line == NULL) { /* skip dead features */
97  continue;
98  }
99 
100  if (Map->constraint.type_flag) {
101  /* skip feature by type */
102  if (!(Line->type & Map->constraint.type))
103  continue;
104  }
105 
106  if (Line->type == GV_CENTROID) {
107  G_debug(4, "Centroid");
108 
109  if (line_p != NULL) {
110  int i, found;
111  struct bound_box box;
112  struct boxlist list;
113  struct P_topo_c *topo = (struct P_topo_c *)Line->topo;
114 
115  /* get area bbox */
116  Vect_get_area_box(Map, topo->area, &box);
117  /* search in spatial index for centroid with area bbox */
119  Vect_select_lines_by_box(Map, &box, Line->type, &list);
120 
121  found = 0;
122  for (i = 0; i < list.n_values; i++) {
123  if (list.id[i] == line) {
124  found = i;
125  break;
126  }
127  }
128 
129  Vect_reset_line(line_p);
130  Vect_append_point(line_p, list.box[found].E, list.box[found].N,
131  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  */
179 int V1_read_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
180  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 = OGR_L_GetFeature(ogr_info->layer, fid);
213  if (ogr_info->feature_cache == NULL) {
214  G_warning(_("Unable to get feature geometry, fid %ld"), fid);
215  return -1;
216  }
217  ogr_info->cache.fid = fid;
218  }
219 
220  hGeom = OGR_F_GetGeometryRef(ogr_info->feature_cache);
221  if (hGeom == NULL) {
222  G_warning(_("Unable to get feature geometry, fid %ld"), fid);
223  return -1;
224  }
225 
226  type = read_line(Map, hGeom, offset + 1, line_p);
227  }
228  else {
229  type = get_line_type(Map, fid);
230  }
231 
232  /* category */
233  if (line_c != NULL) {
234  Vect_cat_set(line_c, 1, (int)fid);
235  }
236 
237  return type;
238 #else
239  G_fatal_error(_("GRASS is not compiled with OGR support"));
240  return -1;
241 #endif
242 }
243 
244 #ifdef HAVE_OGR
245 /*!
246  \brief Recursively read feature and add all elements to points_cache and
247  types_cache.
248 
249  ftype: if > 0 use this type (because parts of Polygon are read as
250  wkbLineString)
251 
252  \param Map pointer to Map_info structure
253  \param[out] hGeom OGR geometry
254  \param ftype feature type
255 
256  \return 0 on success
257  \return 1 on error
258  */
259 int cache_feature(struct Map_info *Map, OGRGeometryH hGeom, int ftype)
260 {
261  int line, i, np, ng, tp;
262 
263  struct Format_info_ogr *ogr_info;
264 
265  OGRwkbGeometryType type;
266  OGRGeometryH hGeom2;
267 
268  G_debug(4, "cache_feature() ftype = %d", ftype);
269 
270  ogr_info = &(Map->fInfo.ogr);
271 
272  /* alloc space in lines cache */
273  line = ogr_info->cache.lines_num;
274  if (line == ogr_info->cache.lines_alloc) {
275  ogr_info->cache.lines_alloc += 1;
276  ogr_info->cache.lines = (struct line_pnts **)G_realloc(
277  (void *)ogr_info->cache.lines,
278  ogr_info->cache.lines_alloc * sizeof(struct line_pnts *));
279 
280  ogr_info->cache.lines_types =
281  (int *)G_realloc(ogr_info->cache.lines_types,
282  ogr_info->cache.lines_alloc * sizeof(int));
283 
284  for (i = ogr_info->cache.lines_num; i < ogr_info->cache.lines_alloc;
285  i++)
286  ogr_info->cache.lines[i] = Vect_new_line_struct();
287  }
288  Vect_reset_line(ogr_info->cache.lines[line]);
289 
290  type = wkbFlatten(OGR_G_GetGeometryType(hGeom));
291 
292  switch (type) {
293  case wkbPoint:
294  G_debug(4, "Point");
295  Vect_append_point(ogr_info->cache.lines[line], OGR_G_GetX(hGeom, 0),
296  OGR_G_GetY(hGeom, 0), OGR_G_GetZ(hGeom, 0));
297  ogr_info->cache.lines_types[line] = GV_POINT;
298  ogr_info->cache.lines_num++;
299  return 0;
300  break;
301 
302  case wkbLineString:
303  G_debug(4, "LineString");
304  np = OGR_G_GetPointCount(hGeom);
305  for (i = 0; i < np; i++) {
306  Vect_append_point(ogr_info->cache.lines[line], OGR_G_GetX(hGeom, i),
307  OGR_G_GetY(hGeom, i), OGR_G_GetZ(hGeom, i));
308  }
309 
310  if (ftype > 0) { /* Polygon rings */
311  ogr_info->cache.lines_types[line] = ftype;
312  }
313  else {
314  ogr_info->cache.lines_types[line] = GV_LINE;
315  }
316  ogr_info->cache.lines_num++;
317  return 0;
318  break;
319 
320  case wkbMultiPoint:
321  case wkbMultiLineString:
322  case wkbPolygon:
323  case wkbMultiPolygon:
324  case wkbGeometryCollection:
325  ng = OGR_G_GetGeometryCount(hGeom);
326  G_debug(4, "%d geoms -> next level", ng);
327  if (type == wkbPolygon) {
328  tp = GV_BOUNDARY;
329  }
330  else {
331  tp = -1;
332  }
333  for (i = 0; i < ng; i++) {
334  hGeom2 = OGR_G_GetGeometryRef(hGeom, i);
335  cache_feature(Map, hGeom2, tp);
336  }
337  return 0;
338  break;
339 
340  default:
341  G_warning(_("OGR feature type %d not supported"), type);
342  return 1;
343  break;
344  }
345 }
346 
347 int read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
348  struct line_cats *line_c, int ignore_constraint)
349 {
350  int itype;
351  struct bound_box lbox, mbox;
352  OGRFeatureH hFeature;
353  OGRGeometryH hGeom;
354 
355  struct Format_info_ogr *ogr_info;
356 
357  G_debug(3, "V1_read_next_line_ogr()");
358 
359  if (Map->constraint.region_flag && !ignore_constraint)
360  Vect_get_constraint_box(Map, &mbox);
361 
362  ogr_info = &(Map->fInfo.ogr);
363  while (TRUE) {
364  /* reset data structures */
365  if (line_p != NULL)
366  Vect_reset_line(line_p);
367  if (line_c != NULL)
368  Vect_reset_cats(line_c);
369 
370  /* read feature to cache if necessary */
371  while (ogr_info->cache.lines_next == ogr_info->cache.lines_num) {
372  hFeature = OGR_L_GetNextFeature(ogr_info->layer);
373  if (hFeature == NULL) {
374  return -2; /* nothing to read */
375  }
376 
377  hGeom = OGR_F_GetGeometryRef(hFeature);
378  if (hGeom == NULL) { /* skip feature without geometry */
379  G_warning(_("Feature without geometry. Skipped."));
380  OGR_F_Destroy(hFeature);
381  continue;
382  }
383 
384  /* cache OGR feature */
385  ogr_info->cache.fid = (int)OGR_F_GetFID(hFeature);
386  if (ogr_info->cache.fid == OGRNullFID) {
387  G_warning(_("OGR feature without ID"));
388  }
389 
390  /* cache feature */
391  ogr_info->cache.lines_num = 0;
392  cache_feature(Map, hGeom, -1);
393  G_debug(4, "%d lines read to cache", ogr_info->cache.lines_num);
394  OGR_F_Destroy(hFeature);
395 
396  /* next to be read from cache */
397  ogr_info->cache.lines_next = 0;
398  }
399 
400  /* read next part of the feature */
401  G_debug(4, "read next cached line %d", ogr_info->cache.lines_next);
402  itype = ogr_info->cache.lines_types[ogr_info->cache.lines_next];
403 
404  if (Map->constraint.type_flag && !ignore_constraint) {
405  /* skip feature by type */
406  if (!(itype & Map->constraint.type)) {
407  ogr_info->cache.lines_next++;
408  continue;
409  }
410  }
411 
412  if (Map->constraint.region_flag && !ignore_constraint) {
413  /* skip feature by region */
414  Vect_line_box(ogr_info->cache.lines[ogr_info->cache.lines_next],
415  &lbox);
416 
417  if (!Vect_box_overlap(&lbox, &mbox)) {
418  ogr_info->cache.lines_next++;
419  continue;
420  }
421  }
422 
423  /* skip feature by field - ignored */
424 
425  if (line_p != NULL)
427  line_p, ogr_info->cache.lines[ogr_info->cache.lines_next],
428  GV_FORWARD);
429 
430  if (line_c != NULL && ogr_info->cache.fid != OGRNullFID)
431  Vect_cat_set(line_c, 1, ogr_info->cache.fid);
432 
433  ogr_info->cache.lines_next++;
434  G_debug(4, "next line read, type = %d", itype);
435 
436  return itype;
437  }
438 
439  return -1; /* not reached */
440 }
441 
442 /*!
443  \brief Recursively descend to feature and read the part
444 
445  \param Map pointer to Map_info structure
446  \param hGeom OGR geometry
447  \param offset given offset
448  \param[out] Points container used to store line points within
449 
450  \return feature type
451  \return -1 on error
452  */
453 int read_line(struct Map_info *Map, OGRGeometryH hGeom, long offset,
454  struct line_pnts *Points)
455 {
456  int i, nPoints;
457  int eType, line;
458 
459  const struct Format_info_ogr *ogr_info;
460 
461  OGRGeometryH hGeom2;
462 
463  /* Read coors if hGeom is a simple element (wkbPoint,
464  * wkbLineString) otherwise descend to geometry specified by
465  * offset[offset] */
466 
467  ogr_info = &(Map->fInfo.ogr);
468 
469  eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
470  G_debug(4, "OGR geometry type: %d", eType);
471 
472  switch (eType) {
473  case wkbPoint:
474  G_debug(4, "\t->Point");
475  if (Points) {
476  Vect_append_point(Points, OGR_G_GetX(hGeom, 0),
477  OGR_G_GetY(hGeom, 0), OGR_G_GetZ(hGeom, 0));
478  }
479  return GV_POINT;
480  break;
481 
482  case wkbLineString:
483  G_debug(4, "\t->LineString");
484  if (Points) {
485  nPoints = OGR_G_GetPointCount(hGeom);
486  for (i = 0; i < nPoints; i++) {
487  Vect_append_point(Points, OGR_G_GetX(hGeom, i),
488  OGR_G_GetY(hGeom, i), OGR_G_GetZ(hGeom, i));
489  }
490  }
491  return GV_LINE;
492  break;
493 
494  case wkbPolygon:
495  case wkbMultiPoint:
496  case wkbMultiLineString:
497  case wkbMultiPolygon:
498  case wkbGeometryCollection:
499  G_debug(4, "\t->more geoms -> part %d", ogr_info->offset.array[offset]);
500  hGeom2 = OGR_G_GetGeometryRef(hGeom, ogr_info->offset.array[offset]);
501  line = read_line(Map, hGeom2, offset + 1, Points);
502  if (eType == wkbPolygon || eType == wkbMultiPolygon)
503  return GV_BOUNDARY;
504  if (eType == wkbMultiPoint)
505  return GV_POINT;
506  if (eType == wkbMultiLineString)
507  return GV_LINE;
508  return line;
509  break;
510 
511  default:
512  G_warning(_("OGR feature type '%s' not supported"),
513  OGRGeometryTypeToName(eType));
514  break;
515  }
516 
517  return -1;
518 }
519 
520 /*!
521  \brief Recursively descend to feature and read the part
522 
523  \param Map pointer to Map_info structure
524  \param hGeom OGR geometry
525  \param offset given offset
526  \param[out] Points container used to store line points within
527 
528  \return feature type
529  \return -1 on error
530  */
531 int get_line_type(struct Map_info *Map, long fid)
532 {
533  int eType;
534 
535  const struct Format_info_ogr *ogr_info;
536 
537  OGRFeatureH hFeat;
538  OGRGeometryH hGeom;
539 
540  G_debug(4, "get_line_type() fid = %ld", fid);
541 
542  ogr_info = &(Map->fInfo.ogr);
543 
544  hFeat = OGR_L_GetFeature(ogr_info->layer, fid);
545  if (hFeat == NULL)
546  return -1;
547 
548  hGeom = OGR_F_GetGeometryRef(hFeat);
549  if (hGeom == NULL)
550  return -1;
551 
552  eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
553 
554  OGR_F_Destroy(hFeat);
555 
556  G_debug(4, "OGR Geometry of type: %d", eType);
557 
558  switch (eType) {
559  case wkbPoint:
560  case wkbMultiPoint:
561  return GV_POINT;
562  break;
563 
564  case wkbLineString:
565  case wkbMultiLineString:
566  return GV_LINE;
567  break;
568 
569  case wkbPolygon:
570  case wkbMultiPolygon:
571  case wkbGeometryCollection:
572  return GV_BOUNDARY;
573  break;
574 
575  default:
576  G_warning(_("OGR feature type %d not supported"), eType);
577  break;
578  }
579 
580  return -1;
581 }
582 #endif
#define NULL
Definition: ccmath.h:32
#define G_realloc(p, n)
Definition: defs/gis.h:96
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
int Vect_reset_cats(struct line_cats *)
Reset category structure to make sure cats structure is clean to be re-used.
void Vect_line_box(const struct line_pnts *, struct bound_box *)
Get bounding box of line.
Definition: line.c:888
int Vect_cat_set(struct line_cats *, int, int)
Add new field/cat to category structure if doesn't exist yet.
int Vect_get_constraint_box(struct Map_info *, struct bound_box *)
Get constraint box.
Definition: constraint.c:79
int Vect_box_overlap(const struct bound_box *, const struct bound_box *)
Tests for overlap of two boxes.
int Vect_get_area_box(struct Map_info *, int, struct bound_box *)
Get bounding box of area.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
int Vect_select_lines_by_box(struct Map_info *, const struct bound_box *, int, struct boxlist *)
Select lines with bounding boxes by box.
Definition: sindex.c:32
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:129
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:148
int Vect_append_points(struct line_pnts *, const struct line_pnts *, int)
Appends points to the end of a line.
Definition: line.c:335
#define GV_CENTROID
Definition: dig_defines.h:186
#define GV_LINE
Definition: dig_defines.h:184
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:183
#define GV_BOUNDARY
Definition: dig_defines.h:185
#define GV_FORWARD
Line direction indicator forward/backward.
Definition: dig_defines.h:179
int dig_init_boxlist(struct boxlist *, int)
#define TRUE
Definition: gis.h:79
#define FALSE
Definition: gis.h:83
#define _(str)
Definition: glocale.h:10
int V1_read_line_ogr(struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c, off_t offset)
Read feature from OGR layer at given offset (level 1 without topology)
Definition: read_ogr.c:179
int V1_read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c)
Read next feature from OGR layer. Skip empty features (level 1 without topology).
Definition: read_ogr.c:49
int V2_read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c)
Read next feature from OGR layer on topological level.
Definition: read_ogr.c:75
int lines_next
Next line to be read from cache.
Definition: dig_structs.h:482
int * lines_types
List of line types (GV_POINT, GV_LINE, ...)
Definition: dig_structs.h:466
long fid
Feature id.
Definition: dig_structs.h:486
struct line_pnts ** lines
Lines array.
Definition: dig_structs.h:462
int lines_alloc
Number of allocated lines in cache.
Definition: dig_structs.h:474
int lines_num
Number of lines which forms current feature.
Definition: dig_structs.h:478
int * array
Offset list.
Definition: dig_structs.h:436
int array_num
Number of items in offset list.
Definition: dig_structs.h:440
Non-native format info (OGR)
Definition: dig_structs.h:505
OGRFeatureH feature_cache
Cache to avoid repeated reading (level 2)
Definition: dig_structs.h:569
OGRLayerH layer
Pointer to OGRLayer.
Definition: dig_structs.h:534
struct Format_info_offset offset
Offset list used for building pseudo-topology.
Definition: dig_structs.h:577
struct Format_info_cache cache
Lines cache for reading feature.
Definition: dig_structs.h:561
struct Format_info_ogr ogr
OGR info.
Definition: dig_structs.h:708
Vector map info.
Definition: dig_structs.h:1243
plus_t next_line
Feature id for sequential access.
Definition: dig_structs.h:1338
int type
Feature type constraint.
Definition: dig_structs.h:1359
int type_flag
Non-zero value to enable feature type constraint.
Definition: dig_structs.h:1355
struct Map_info::@11 constraint
Constraints for sequential feature access.
int region_flag
Non-zero value to enable region constraint.
Definition: dig_structs.h:1347
struct Format_info fInfo
Format info for non-native formats.
Definition: dig_structs.h:1400
struct Plus_head plus
Plus info (topology, version, ...)
Definition: dig_structs.h:1270
Vector geometry.
Definition: dig_structs.h:1553
char type
Line type.
Definition: dig_structs.h:1564
off_t offset
Offset in coor file for line.
Definition: dig_structs.h:1571
void * topo
Topology info.
Definition: dig_structs.h:1577
Centroid topology.
Definition: dig_structs.h:1514
plus_t area
Area number, negative for duplicate centroid.
Definition: dig_structs.h:1518
struct P_line ** Line
Array of vector geometries.
Definition: dig_structs.h:871
plus_t n_lines
Current number of lines.
Definition: dig_structs.h:931
Bounding box.
Definition: dig_structs.h:64
List of bounding boxes with id.
Definition: dig_structs.h:1723
Feature category info.
Definition: dig_structs.h:1677
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
Definition: manage.h:4