GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-87d2d42d3b
geos.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/geos.c
3 
4  \brief Vector library - GEOS support
5 
6  Higher level functions for reading/writing/manipulating vectors.
7 
8  (C) 2009 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 Martin Landa <landa.martin gmail.com>
14  */
15 
16 #include <stdlib.h>
17 #include <grass/vector.h>
18 #include <grass/glocale.h>
19 
20 #ifdef HAVE_GEOS
21 
22 static GEOSGeometry *Vect__read_line_geos(struct Map_info *, long, int *);
23 static GEOSCoordSequence *V1_read_line_geos(struct Map_info *, long, int *);
24 static GEOSCoordSequence *V2_read_line_geos(struct Map_info *, int);
25 static GEOSCoordSequence *read_polygon_points(struct Map_info *, int, int *);
26 
27 /*!
28  \brief Read vector feature and stores it as GEOSGeometry instance
29 
30  Supported feature types:
31  - GV_POINT -> POINT
32  - GV_LINE -> LINESTRING
33  - GV_BOUNDARY -> LINESTRING / LINEARRING
34 
35  You should free allocated memory by GEOSGeom_destroy().
36 
37  \param Map pointer to Map_info structure
38  \param line feature id
39  \param[out] type feature type or NULL
40 
41  \return pointer to GEOSGeometry instance
42  \return empty GEOSGeometry for unsupported feature type
43  \return NULL on error
44  */
45 GEOSGeometry *Vect_read_line_geos(struct Map_info *Map, int line, int *type)
46 {
47  struct P_line *Line;
48 
49  G_debug(3, "Vect_read_line_geos(): line = %d", line);
50 
51  if (!VECT_OPEN(Map))
52  G_fatal_error("Vect_read_line_geos(): %s",
53  _("vector map is not opened"));
54 
55  if (line < 1 || line > Map->plus.n_lines)
57  _("Vect_read_line_geos(): feature id %d is not reasonable "
58  "(max features in vector map <%s>: %d)"),
59  line, Vect_get_full_name(Map), Map->plus.n_lines);
60 
61  if (Map->format != GV_FORMAT_NATIVE)
62  G_fatal_error("Vect_read_line_geos(): %s",
63  _("only native format supported"));
64 
65  Line = Map->plus.Line[line];
66  if (Line == NULL)
67  G_fatal_error("Vect_read_line_geos(): %s %d",
68  _("Attempt to read dead line"), line);
69 
70  return Vect__read_line_geos(Map, Line->offset, type);
71 }
72 
73 /*!
74  \brief Read vector area and stores it as GEOSGeometry instance (polygon)
75 
76  You should free allocated memory by GEOSGeom_destroy().
77 
78  \param Map pointer to Map_info structure
79  \param area area id
80 
81  \return pointer to GEOSGeometry instance
82  \return NULL on error
83  */
84 GEOSGeometry *Vect_read_area_geos(struct Map_info *Map, int area)
85 {
86  int i, nholes, isle;
87  GEOSGeometry *boundary, *poly, **holes;
88 
89  G_debug(3, "Vect_read_area_geos(): area = %d", area);
90 
91  boundary = GEOSGeom_createLinearRing(Vect_get_area_points_geos(Map, area));
92  if (!boundary) {
93  G_fatal_error(_("Vect_read_area_geos(): unable to read area id %d"),
94  area);
95  }
96 
97  nholes = Vect_get_area_num_isles(Map, area);
98  holes = (GEOSGeometry **)G_malloc(nholes * sizeof(GEOSGeometry *));
99  for (i = 0; i < nholes; i++) {
100  isle = Vect_get_area_isle(Map, area, i);
101  if (isle < 1) {
102  nholes--;
103  continue;
104  }
105  holes[i] =
106  GEOSGeom_createLinearRing(Vect_get_isle_points_geos(Map, isle));
107  if (!(holes[i]))
108  G_fatal_error(_("Vect_read_area_geos(): unable to read isle id %d "
109  "of area id %d"),
110  isle, area);
111  }
112 
113  poly = GEOSGeom_createPolygon(boundary, holes, nholes);
114  G_free(holes);
115 
116  return poly;
117 }
118 
119 /*!
120  \brief Create GEOSGeometry of given type from feature points.
121 
122  Supported types:
123  - GV_POINT -> POINT
124  - GV_CENTROID -> POINT
125  - GV_LINE -> LINESTRING
126  - GV_BOUNDARY -> LINEARRING
127 
128  You should free allocated memory by GEOSGeom_destroy().
129 
130  \param points pointer to line_pnts structure
131  \param type feature type (see supported types)
132  \param with_z Set to 1 if the feature is 3d, 0 otherwise
133 
134  \return pointer to GEOSGeometry instance
135  \return NULL on error
136  */
137 GEOSGeometry *Vect_line_to_geos(const struct line_pnts *points, int type,
138  int with_z)
139 {
140  int i;
141  GEOSGeometry *geom;
142  GEOSCoordSequence *pseq;
143 
144  G_debug(3, "Vect_line_to_geos(): type = %d", type);
145 
146  /* read only points / lines / boundaries */
147  if (!(type & (GV_POINT | GV_CENTROID | GV_LINES)))
148  return NULL;
149 
150  if (type == GV_POINT || type == GV_CENTROID) {
151  if (points->n_points != 1)
152  /* point is not valid */
153  return NULL;
154  }
155  else {
156  if (points->n_points < 2)
157  /* line/boundary is not valid */
158  return NULL;
159  }
160 
161  pseq = GEOSCoordSeq_create(points->n_points, with_z ? 3 : 2);
162 
163  for (i = 0; i < points->n_points; i++) {
164  GEOSCoordSeq_setX(pseq, i, points->x[i]);
165  GEOSCoordSeq_setY(pseq, i, points->y[i]);
166  if (with_z)
167  GEOSCoordSeq_setZ(pseq, i, points->z[i]);
168  }
169 
170  if (type == GV_POINT || type == GV_CENTROID)
171  geom = GEOSGeom_createPoint(pseq);
172  else if (type == GV_LINE)
173  geom = GEOSGeom_createLineString(pseq);
174  else { /* boundary */
175  geom = GEOSGeom_createLineString(pseq);
176  if (GEOSisRing(geom)) {
177  /*GEOSGeom_destroy(geom); */
178  geom = GEOSGeom_createLinearRing(pseq);
179  }
180  }
181 
182  /* GEOSCoordSeq_destroy(pseq); */
183 
184  return geom;
185 }
186 
187 /*!
188  \brief Read line from coor file
189 
190  You should free allocated memory by GEOSGeom_destroy().
191 
192  \param Map pointer to Map_info
193  \param offset line offset
194  \param[out] type feature type or NULL
195 
196  \return pointer to GEOSGeometry
197  \return NULL on error
198  \return NULL dead line
199  \return NULL end of file
200  */
201 GEOSGeometry *Vect__read_line_geos(struct Map_info *Map, long offset, int *type)
202 {
203  int ftype;
204 
205  GEOSGeometry *geom;
206  GEOSCoordSequence *pseq;
207 
208  pseq = V1_read_line_geos(Map, offset, &ftype);
209  if (!pseq)
210  G_fatal_error(_("Unable to read line offset %ld"), offset);
211 
212  if (ftype & GV_POINT) {
213  G_debug(3, " geos_type = point");
214  geom = GEOSGeom_createPoint(pseq);
215  }
216  else if (ftype & GV_LINE) {
217  G_debug(3, " geos_type = linestring");
218  geom = GEOSGeom_createLineString(pseq);
219  }
220  else { /* boundary */
221  geom = GEOSGeom_createLineString(pseq);
222  if (GEOSisRing(geom)) {
223  /* GEOSGeom_destroy(geom); */
224  geom = GEOSGeom_createLinearRing(pseq);
225  G_debug(3, " geos_type = linearring");
226  }
227  else {
228  G_debug(3, " geos_type = linestring");
229  }
230  }
231 
232  /* GEOSCoordSeq_destroy(pseq); */
233 
234  if (type)
235  *type = ftype;
236 
237  return geom;
238 }
239 
240 /*!
241  \brief Read line from coor file into GEOSCoordSequence
242 
243  You should free allocated memory by GEOSCoordSeq_destroy().
244 
245  \param Map pointer to Map_info
246  \param line line id
247 
248  \return pointer to GEOSCoordSequence
249  \return empty GEOSCoordSequence for dead line or unsuppored feature type
250  \return NULL end of file
251  */
252 GEOSCoordSequence *V2_read_line_geos(struct Map_info *Map, int line)
253 {
254  int ftype;
255  struct P_line *Line;
256 
257  G_debug(3, "V2_read_line_geos(): line = %d", line);
258 
259  Line = Map->plus.Line[line];
260 
261  if (Line == NULL)
262  G_fatal_error("V2_read_line_geos(): %s %d",
263  _("Attempt to read dead line"), line);
264 
265  return V1_read_line_geos(Map, Line->offset, &ftype);
266 }
267 
268 /*!
269  \brief Read feature from coor file into GEOSCoordSequence
270 
271  Note: Function reads only points, lines and boundaries, other
272  feature types are ignored (empty coord array is returned)!
273 
274  You should free allocated memory by GEOSCoordSeq_destroy().
275 
276  \param Map pointer to Map_info
277  \param offset line offset
278  \param[out] type feature type
279 
280  \return pointer to GEOSCoordSequence
281  \return empty GEOSCoordSequence for dead line or unsuppored feature type
282  \return NULL end of file
283  */
284 GEOSCoordSequence *V1_read_line_geos(struct Map_info *Map, long offset,
285  int *type)
286 {
287  int i, n_points;
288  int do_cats, n_cats;
289  char rhead, nc;
290  long size;
291  double *x, *y, *z;
292 
293  GEOSCoordSequence *pseq;
294 
295  G_debug(3, "V1_read_line_geos(): offset = %ld", offset);
296 
297  Map->head.last_offset = offset;
298 
299  /* reads must set in_head, but writes use default */
300  dig_set_cur_port(&(Map->head.port));
301 
302  dig_fseek(&(Map->dig_fp), offset, 0);
303 
304  if (0 >= dig__fread_port_C(&rhead, 1, &(Map->dig_fp)))
305  return NULL; /* end of file */
306 
307  if (!(rhead & 0x01)) /* dead line */
308  return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
309 
310  if (rhead & 0x02) /* categories exists */
311  do_cats =
312  1; /* do not return here let file offset moves forward to next */
313  else /* line */
314  do_cats = 0;
315 
316  rhead >>= 2;
317  *type = dig_type_from_store((int)rhead);
318 
319  /* read only points / lines / boundaries */
320  if (!(*type & (GV_POINT | GV_LINES)))
321  return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
322 
323  /* skip categories */
324  if (do_cats) {
325  if (Map->head.coor_version.minor == 1) { /* coor format 5.1 */
326  if (0 >= dig__fread_port_I(&n_cats, 1, &(Map->dig_fp)))
327  return NULL;
328  }
329  else { /* coor format 5.0 */
330  if (0 >= dig__fread_port_C(&nc, 1, &(Map->dig_fp)))
331  return NULL;
332  n_cats = (int)nc;
333  }
334  G_debug(3, " n_cats = %d", n_cats);
335 
336  if (Map->head.coor_version.minor == 1) { /* coor format 5.1 */
337  size = (2 * PORT_INT) * n_cats;
338  }
339  else { /* coor format 5.0 */
340  size = (PORT_SHORT + PORT_INT) * n_cats;
341  }
342  dig_fseek(&(Map->dig_fp), size, SEEK_CUR);
343  }
344 
345  if (*type & GV_POINTS) {
346  n_points = 1;
347  }
348  else {
349  if (0 >= dig__fread_port_I(&n_points, 1, &(Map->dig_fp)))
350  return NULL;
351  }
352 
353  G_debug(3, " n_points = %d dim = %d", n_points,
354  (Map->head.with_z) ? 3 : 2);
355 
356  pseq = GEOSCoordSeq_create(n_points, (Map->head.with_z) ? 3 : 2);
357 
358  x = (double *)G_malloc(n_points * sizeof(double));
359  y = (double *)G_malloc(n_points * sizeof(double));
360  if (Map->head.with_z)
361  z = (double *)G_malloc(n_points * sizeof(double));
362  else
363  z = NULL;
364 
365  if (0 >= dig__fread_port_D(x, n_points, &(Map->dig_fp)))
366  return NULL; /* end of file */
367 
368  if (0 >= dig__fread_port_D(y, n_points, &(Map->dig_fp)))
369  return NULL; /* end of file */
370 
371  if (Map->head.with_z) {
372  if (0 >= dig__fread_port_D(z, n_points, &(Map->dig_fp)))
373  return NULL; /* end of file */
374  }
375 
376  for (i = 0; i < n_points; i++) {
377  GEOSCoordSeq_setX(pseq, i, x[i]);
378  GEOSCoordSeq_setY(pseq, i, y[i]);
379  if (Map->head.with_z)
380  GEOSCoordSeq_setZ(pseq, i, z[i]);
381  }
382 
383  G_debug(3, " off = %ld", (long)dig_ftell(&(Map->dig_fp)));
384 
385  G_free((void *)x);
386  G_free((void *)y);
387  if (z)
388  G_free((void *)z);
389 
390  return pseq;
391 }
392 
393 /*!
394  \brief Returns the polygon array of points, i.e. outer ring (shell)
395 
396  You should free allocated memory by GEOSCoordSeq_destroy().
397 
398  See also Vect_get_area_points().
399 
400  \param Map pointer to Map_info
401  \param area area id
402 
403  \return pointer to GEOSCoordSequence
404  \return empty GEOSCoordSequence for dead area
405  \return NULL on error
406  */
408 {
409  struct Plus_head *Plus;
410  struct P_area *Area;
411 
412  G_debug(3, "Vect_get_area_points_geos(): area = %d", area);
413 
414  Plus = &(Map->plus);
415  Area = Plus->Area[area];
416 
417  if (Area == NULL) { /* dead area */
418  G_warning(_("Attempt to read points of nonexistent area id %d"), area);
419  return NULL; /* error , because we should not read dead areas */
420  }
421 
422  return read_polygon_points(Map, Area->n_lines, Area->lines);
423 }
424 
425 /*!
426  \brief Returns the polygon (isle) array of points (inner ring)
427 
428  You should free allocated memory by GEOSCoordSeq_destroy().
429 
430  See also Vect_get_isle_points().
431 
432  \param Map pointer to Map_info
433  \param isle isel id
434 
435  \return pointer to GEOSGeometry
436  \return NULL on error or dead line
437  */
439 {
440  struct Plus_head *Plus;
441  struct P_isle *Isle;
442 
443  G_debug(3, "Vect_get_isle_points_geos(): isle = %d", isle);
444 
445  Plus = &(Map->plus);
446  Isle = Plus->Isle[isle];
447 
448  return read_polygon_points(Map, Isle->n_lines, Isle->lines);
449 }
450 
451 GEOSCoordSequence *read_polygon_points(struct Map_info *Map, int n_lines,
452  int *lines)
453 {
454  int i, j, k;
455  int line, aline;
456  unsigned int n_points, n_points_shell;
457  double x, y, z;
458  int *dir;
459 
460  GEOSCoordSequence **pseq, *pseq_shell;
461 
462  G_debug(3, " n_lines = %d", n_lines);
463  pseq =
465  dir = (int *)G_malloc(n_lines * sizeof(int));
466 
467  n_points_shell = 0;
468  for (i = 0; i < n_lines; i++) {
469  line = lines[i];
470  aline = abs(line);
471  G_debug(3, " append line(%d) = %d", i, line);
472 
473  if (line > 0)
474  dir[i] = GV_FORWARD;
475  else
476  dir[i] = GV_BACKWARD;
477 
478  pseq[i] = V2_read_line_geos(Map, aline);
479  if (!(pseq[i])) {
480  G_fatal_error(_("Unable to read feature id %d"), aline);
481  }
482 
483  GEOSCoordSeq_getSize(pseq[i], &n_points);
484  G_debug(3, " line n_points = %d", n_points);
485  n_points_shell += n_points;
486  }
487 
488  /* create shell (outer ring) */
489  pseq_shell = GEOSCoordSeq_create(n_points_shell, Map->head.with_z ? 3 : 2);
490  k = 0;
491  for (i = 0; i < n_lines; i++) {
492  GEOSCoordSeq_getSize(pseq[i], &n_points);
493  if (dir[i] == GV_FORWARD) {
494  for (j = 0; j < (int)n_points; j++, k++) {
495  GEOSCoordSeq_getX(pseq[i], j, &x);
496  GEOSCoordSeq_setX(pseq_shell, k, x);
497 
498  GEOSCoordSeq_getY(pseq[i], j, &y);
499  GEOSCoordSeq_setY(pseq_shell, k, y);
500 
501  if (Map->head.with_z) {
502  GEOSCoordSeq_getY(pseq[i], j, &z);
503  GEOSCoordSeq_setZ(pseq_shell, k, z);
504  }
505  }
506  }
507  else { /* GV_BACKWARD */
508  for (j = (int)n_points - 1; j > -1; j--, k++) {
509  GEOSCoordSeq_getX(pseq[i], j, &x);
510  GEOSCoordSeq_setX(pseq_shell, k, x);
511 
512  GEOSCoordSeq_getY(pseq[i], j, &y);
513  GEOSCoordSeq_setY(pseq_shell, k, y);
514 
515  if (Map->head.with_z) {
516  GEOSCoordSeq_getY(pseq[i], j, &z);
517  GEOSCoordSeq_setZ(pseq_shell, k, z);
518  }
519  }
520  }
521  GEOSCoordSeq_destroy(pseq[i]);
522  }
523 
524  G_free((void *)pseq);
525  G_free((void *)dir);
526 
527  return pseq_shell;
528 }
529 #endif /* HAVE_GEOS */
#define NULL
Definition: ccmath.h:32
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
#define G_malloc(n)
Definition: defs/gis.h:94
int G_debug(int, const char *,...) __attribute__((format(printf
const char * Vect_get_full_name(struct Map_info *)
Get fully qualified name of vector map.
int Vect_get_area_isle(struct Map_info *, int, int)
Returns isle id for area.
int Vect_get_area_num_isles(struct Map_info *, int)
Returns number of isles for given area.
#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 VECT_OPEN(Map)
Check if vector map is open.
Definition: dig_defines.h:137
#define GV_LINES
Definition: dig_defines.h:193
#define PORT_SHORT
Definition: dig_defines.h:49
#define PORT_INT
Definition: dig_defines.h:48
#define GV_FORWARD
Line direction indicator forward/backward.
Definition: dig_defines.h:179
#define GV_POINTS
Definition: dig_defines.h:192
#define GV_BACKWARD
Definition: dig_defines.h:180
#define GV_FORMAT_NATIVE
Geometry data formats supported by lib Don't change GV_FORMAT_* values, this order is hardcoded in li...
Definition: dig_defines.h:83
int dig__fread_port_D(double *, size_t, struct gvfile *)
Read doubles from the Portable Vector Format.
Definition: portable.c:79
off_t dig_ftell(struct gvfile *file)
Get struct gvfile position.
Definition: file.c:36
int dig_set_cur_port(struct Port_info *)
Set current Port_info structure.
Definition: portable.c:996
int dig__fread_port_C(char *, size_t, struct gvfile *)
Read chars from the Portable Vector Format.
Definition: portable.c:511
int dig__fread_port_I(int *, size_t, struct gvfile *)
Read integers from the Portable Vector Format.
Definition: portable.c:345
int dig_fseek(struct gvfile *file, off_t offset, int whence)
Set struct gvfile position.
Definition: file.c:60
int dig_type_from_store(int)
Convert type from store type.
GEOSCoordSequence * Vect_get_isle_points_geos(struct Map_info *Map, int isle)
Returns the polygon (isle) array of points (inner ring)
Definition: geos.c:438
GEOSGeometry * Vect_read_line_geos(struct Map_info *Map, int line, int *type)
Read vector feature and stores it as GEOSGeometry instance.
Definition: geos.c:45
GEOSGeometry * Vect_read_area_geos(struct Map_info *Map, int area)
Read vector area and stores it as GEOSGeometry instance (polygon)
Definition: geos.c:84
GEOSCoordSequence * Vect_get_area_points_geos(struct Map_info *Map, int area)
Returns the polygon array of points, i.e. outer ring (shell)
Definition: geos.c:407
GEOSGeometry * Vect_line_to_geos(const struct line_pnts *points, int type, int with_z)
Create GEOSGeometry of given type from feature points.
Definition: geos.c:137
#define _(str)
Definition: glocale.h:10
Vector map info.
Definition: dig_structs.h:1243
struct gvfile dig_fp
GV file pointer (native format only)
Definition: dig_structs.h:1395
struct dig_head head
Header info.
Definition: dig_structs.h:1388
int format
Map format (native, ogr, postgis)
Definition: dig_structs.h:1255
struct Plus_head plus
Plus info (topology, version, ...)
Definition: dig_structs.h:1270
Area (topology) info.
Definition: dig_structs.h:1583
plus_t n_lines
Number of boundary lines.
Definition: dig_structs.h:1587
plus_t * lines
List of boundary lines.
Definition: dig_structs.h:1598
Isle (topology) info.
Definition: dig_structs.h:1623
plus_t * lines
List of boundary lines.
Definition: dig_structs.h:1638
plus_t n_lines
Number of boundary lines.
Definition: dig_structs.h:1627
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
Basic topology-related info.
Definition: dig_structs.h:769
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
struct P_area ** Area
Array of areas.
Definition: dig_structs.h:875
struct P_isle ** Isle
Array of isles.
Definition: dig_structs.h:879
int minor
Current version (minor)
Definition: dig_structs.h:275
struct Version_info coor_version
Version info for coor file.
Definition: dig_structs.h:331
off_t last_offset
Offset of last read line.
Definition: dig_structs.h:358
int with_z
2D/3D vector data
Definition: dig_structs.h:339
struct Port_info port
Portability information.
Definition: dig_structs.h:353
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
double * y
Array of Y coordinates.
Definition: dig_structs.h:1659
double * x
Array of X coordinates.
Definition: dig_structs.h:1655
int n_points
Number of points.
Definition: dig_structs.h:1667
double * z
Array of Z coordinates.
Definition: dig_structs.h:1663
struct GEOSCoordSeq_t GEOSCoordSequence
Definition: vector.h:10
struct GEOSGeom_t GEOSGeometry
Definition: vector.h:9
#define x