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 */
Vect_read_line_geos(struct Map_info * Map,int line,int * type)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 */
Vect_read_area_geos(struct Map_info * Map,int area)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 */
Vect_line_to_geos(const struct line_pnts * points,int type,int with_z)132 GEOSGeometry *Vect_line_to_geos(const struct line_pnts *points,
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 */
Vect__read_line_geos(struct Map_info * Map,long offset,int * type)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 */
V2_read_line_geos(struct Map_info * Map,int line)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 */
V1_read_line_geos(struct Map_info * Map,long offset,int * type)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 */
Vect_get_area_points_geos(struct Map_info * Map,int area)401 GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *Map, int area)
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 */
Vect_get_isle_points_geos(struct Map_info * Map,int isle)432 GEOSCoordSequence *Vect_get_isle_points_geos(struct Map_info *Map, int isle)
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
read_polygon_points(struct Map_info * Map,int n_lines,int * lines)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 */
522