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