GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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(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 */
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 */
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 */
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 */
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 
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 */
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 */
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
int Vect_get_area_box(const struct Map_info *, int, struct bound_box *)
Get bounding box of area.
int Vect_append_points(struct line_pnts *, const struct line_pnts *, int)
Appends points to the end of a line.
Definition: line.c:337
int * array
Offset list.
Definition: dig_structs.h:446
#define TRUE
Definition: gis.h:59
Bounding box.
Definition: dig_structs.h:65
#define GV_FORWARD
Line direction indicator forward/backward.
Definition: dig_defines.h:178
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
plus_t area
Area number, negative for duplicate centroid.
Definition: dig_structs.h:1537
int * id
Array of ids.
Definition: dig_structs.h:1755
off_t offset
Offset in coor file for line.
Definition: dig_structs.h:1593
Vector geometry.
Definition: dig_structs.h:1574
struct P_line ** Line
Array of vector geometries.
Definition: dig_structs.h:887
int region_flag
Non-zero value to enable region constraint.
Definition: dig_structs.h:1362
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:76
struct Format_info_cache cache
Lines cache for reading feature.
Definition: dig_structs.h:573
struct Format_info_offset offset
Offset list used for building pseudo-topology.
Definition: dig_structs.h:589
long fid
Feature id.
Definition: dig_structs.h:497
#define GV_CENTROID
Definition: dig_defines.h:185
struct Format_info fInfo
Format info for non-native formats.
Definition: dig_structs.h:1415
double E
East.
Definition: dig_structs.h:78
void Vect_line_box(const struct line_pnts *, struct bound_box *)
Get bounding box of line.
Definition: line.c:922
#define NULL
Definition: ccmath.h:32
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 type
Feature type constraint.
Definition: dig_structs.h:1374
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
struct bound_box * box
Array of bounding boxes.
Definition: dig_structs.h:1759
int lines_num
Number of lines which forms current feature.
Definition: dig_structs.h:489
plus_t n_lines
Current number of lines.
Definition: dig_structs.h:947
Feature category info.
Definition: dig_structs.h:1702
#define GV_LINE
Definition: dig_defines.h:183
double N
North.
Definition: dig_structs.h:70
char type
Line type.
Definition: dig_structs.h:1586
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
Centroid topology.
Definition: dig_structs.h:1532
plus_t next_line
Feature id for sequential access.
Definition: dig_structs.h:1353
struct Map_info::@11 constraint
Constraints for sequential feature access.
#define FALSE
Definition: gis.h:63
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
void * topo
Topology info.
Definition: dig_structs.h:1599
int lines_alloc
Number of allocated lines in cache.
Definition: dig_structs.h:485
int Vect_get_constraint_box(const struct Map_info *, struct bound_box *)
Get constraint box.
Definition: constraint.c:79
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:149
#define GV_BOUNDARY
Definition: dig_defines.h:184
struct Plus_head plus
Plus info (topology, version, ...)
Definition: dig_structs.h:1286
Non-native format info (OGR)
Definition: dig_structs.h:516
struct Format_info_ogr ogr
OGR info.
Definition: dig_structs.h:722
int Vect_reset_cats(struct line_cats *)
Reset category structure to make sure cats structure is clean to be re-used.
int dig_init_boxlist(struct boxlist *, int)
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:50
int type_flag
Non-zero value to enable feature type constraint.
Definition: dig_structs.h:1370
int array_num
Number of items in offset list.
Definition: dig_structs.h:450
int * lines_types
List of line types (GV_POINT, GV_LINE, ...)
Definition: dig_structs.h:477
Vector map info.
Definition: dig_structs.h:1259
OGRFeatureH feature_cache
Cache to avoid repeated reading (level 2)
Definition: dig_structs.h:581
List of bounding boxes with id.
Definition: dig_structs.h:1750
void G_warning(const char *,...) __attribute__((format(printf
int lines_next
Next line to be read from cache.
Definition: dig_structs.h:493
struct line_pnts ** lines
Lines array.
Definition: dig_structs.h:473
#define G_realloc(p, n)
Definition: defs/gis.h:114
#define _(str)
Definition: glocale.h:10
OGRLayerH layer
Pointer to OGRLayer.
Definition: dig_structs.h:546
int Vect_box_overlap(const struct bound_box *, const struct bound_box *)
Tests for overlap of two boxes.
int G_debug(int, const char *,...) __attribute__((format(printf
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:130
int Vect_cat_set(struct line_cats *, int, int)
Add new field/cat to category structure if doesn&#39;t exist yet.