1 /******************************************************************************
2 * $Id$
3 *
4 * Project: MapServer
5 * Purpose: MapServer-GEOS integration.
6 * Author: Steve Lime and the MapServer team.
7 *
8 ******************************************************************************
9 * Copyright (c) 1996-2005 Regents of the University of Minnesota.
10 *
11 * Permission is hereby granted, free of charge, to any person obtaining a
12 * copy of this software and associated documentation files (the "Software"),
13 * to deal in the Software without restriction, including without limitation
14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 * and/or sell copies of the Software, and to permit persons to whom the
16 * Software is furnished to do so, subject to the following conditions:
17 *
18 * The above copyright notice and this permission notice shall be included in
19 * all copies of this Software or works derived from this Software.
20 *
21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 * DEALINGS IN THE SOFTWARE.
28 *****************************************************************************/
29
30 #include "mapserver.h"
31
32 #ifdef USE_GEOS
33
34 // To avoid accidental use of non reentrant GEOS API.
35 // (check only effective in GEOS >= 3.5)
36 #define GEOS_USE_ONLY_R_API
37
38 #include <geos_c.h>
39
40 /*
41 ** Error handling...
42 */
msGEOSError(const char * format,...)43 static void msGEOSError(const char *format, ...)
44 {
45 va_list args;
46
47 va_start (args, format);
48 msSetError(MS_GEOSERR, format, "msGEOSError()", args); /* just pass along to MapServer error handling */
49 va_end(args);
50
51 return;
52 }
53
msGEOSNotice(const char * fmt,...)54 static void msGEOSNotice(const char *fmt, ...)
55 {
56 return; /* do nothing with notices at this point */
57 }
58
59 #ifndef USE_THREAD
60 static GEOSContextHandle_t geos_handle;
61
msGetGeosContextHandle()62 static inline GEOSContextHandle_t msGetGeosContextHandle()
63 {
64 return geos_handle;
65 }
66
67 #else
68 #include "mapthread.h"
69 typedef struct geos_thread_info {
70 struct geos_thread_info *next;
71 void* thread_id;
72 GEOSContextHandle_t geos_handle;
73 } geos_thread_info_t;
74
75 static geos_thread_info_t *geos_list = NULL;
76
msGetGeosContextHandle()77 static GEOSContextHandle_t msGetGeosContextHandle()
78 {
79 geos_thread_info_t *link;
80 GEOSContextHandle_t ret_obj;
81 void* thread_id;
82
83 msAcquireLock( TLOCK_GEOS );
84
85 thread_id = msGetThreadId();
86
87 /* find link for this thread */
88
89 for( link = geos_list;
90 link != NULL && link->thread_id != thread_id
91 && link->next != NULL && link->next->thread_id != thread_id;
92 link = link->next ) {}
93
94 /* If the target thread link is already at the head of the list were ok */
95 if( geos_list != NULL && geos_list->thread_id == thread_id ) {
96 }
97
98 /* We don't have one ... initialize one. */
99 else if( link == NULL || link->next == NULL ) {
100 geos_thread_info_t *new_link;
101 new_link = (geos_thread_info_t *) malloc(sizeof(geos_thread_info_t));
102 new_link->next = geos_list;
103 new_link->thread_id = thread_id;
104 new_link->geos_handle = initGEOS_r(msGEOSNotice, msGEOSError);
105
106 geos_list = new_link;
107 }
108
109 /* If the link is not already at the head of the list, promote it */
110 else if( link != NULL && link->next != NULL ) {
111 geos_thread_info_t *target = link->next;
112
113 link->next = link->next->next;
114 target->next = geos_list;
115 geos_list = target;
116 }
117
118 ret_obj = geos_list->geos_handle;
119
120 msReleaseLock( TLOCK_GEOS );
121
122 return ret_obj;
123 }
124 #endif
125
126 /*
127 ** Setup/Cleanup wrappers
128 */
msGEOSSetup()129 void msGEOSSetup()
130 {
131 #ifndef USE_THREAD
132 geos_handle = initGEOS_r(msGEOSNotice, msGEOSError);
133 #endif
134 }
135
msGEOSCleanup()136 void msGEOSCleanup()
137 {
138 #ifndef USE_THREAD
139 finishGEOS_r(geos_handle);
140 geos_handle = NULL;
141 #else
142 geos_thread_info_t *link;
143 msAcquireLock( TLOCK_GEOS );
144 for( link = geos_list; link != NULL;) {
145 geos_thread_info_t *cur = link;
146 link = link->next;
147 finishGEOS_r(cur->geos_handle);
148 free(cur);
149 }
150 geos_list = NULL;
151 msReleaseLock( TLOCK_GEOS );
152 #endif
153 }
154
155 /*
156 ** Translation functions
157 */
msGEOSShape2Geometry_point(pointObj * point)158 static GEOSGeom msGEOSShape2Geometry_point(pointObj *point)
159 {
160 GEOSCoordSeq coords;
161 GEOSGeom g;
162 GEOSContextHandle_t handle = msGetGeosContextHandle();
163
164 if(!point) return NULL;
165
166 coords = GEOSCoordSeq_create_r(handle,1, 2); /* todo handle z's */
167 if(!coords) return NULL;
168
169 GEOSCoordSeq_setX_r(handle,coords, 0, point->x);
170 GEOSCoordSeq_setY_r(handle,coords, 0, point->y);
171 /* GEOSCoordSeq_setY(coords, 0, point->z); */
172
173 g = GEOSGeom_createPoint_r(handle,coords); /* g owns the coordinate in coords */
174
175 return g;
176 }
177
msGEOSShape2Geometry_multipoint(lineObj * multipoint)178 static GEOSGeom msGEOSShape2Geometry_multipoint(lineObj *multipoint)
179 {
180 int i;
181 GEOSGeom g;
182 GEOSGeom *points;
183 GEOSContextHandle_t handle = msGetGeosContextHandle();
184
185 if(!multipoint) return NULL;
186
187 points = malloc(multipoint->numpoints*sizeof(GEOSGeom));
188 if(!points) return NULL;
189
190 for(i=0; i<multipoint->numpoints; i++)
191 points[i] = msGEOSShape2Geometry_point(&(multipoint->point[i]));
192
193 g = GEOSGeom_createCollection_r(handle,GEOS_MULTIPOINT, points, multipoint->numpoints);
194
195 free(points);
196
197 return g;
198 }
199
msGEOSShape2Geometry_line(lineObj * line)200 static GEOSGeom msGEOSShape2Geometry_line(lineObj *line)
201 {
202 int i;
203 GEOSGeom g;
204 GEOSCoordSeq coords;
205 GEOSContextHandle_t handle = msGetGeosContextHandle();
206
207 if(!line) return NULL;
208
209 coords = GEOSCoordSeq_create_r(handle,line->numpoints, 2); /* todo handle z's */
210 if(!coords) return NULL;
211
212 for(i=0; i<line->numpoints; i++) {
213 GEOSCoordSeq_setX_r(handle,coords, i, line->point[i].x);
214 GEOSCoordSeq_setY_r(handle,coords, i, line->point[i].y);
215 /* GEOSCoordSeq_setZ(coords, i, line->point[i].z); */
216 }
217
218 g = GEOSGeom_createLineString_r(handle,coords); /* g owns the coordinates in coords */
219
220 return g;
221 }
222
msGEOSShape2Geometry_multiline(shapeObj * multiline)223 static GEOSGeom msGEOSShape2Geometry_multiline(shapeObj *multiline)
224 {
225 int i;
226 GEOSGeom g;
227 GEOSGeom *lines;
228 GEOSContextHandle_t handle = msGetGeosContextHandle();
229
230 if(!multiline) return NULL;
231
232 lines = malloc(multiline->numlines*sizeof(GEOSGeom));
233 if(!lines) return NULL;
234
235 for(i=0; i<multiline->numlines; i++)
236 lines[i] = msGEOSShape2Geometry_line(&(multiline->line[i]));
237
238 g = GEOSGeom_createCollection_r(handle,GEOS_MULTILINESTRING, lines, multiline->numlines);
239
240 free(lines);
241
242 return g;
243 }
244
msGEOSShape2Geometry_simplepolygon(shapeObj * shape,int r,int * outerList)245 static GEOSGeom msGEOSShape2Geometry_simplepolygon(shapeObj *shape, int r, int *outerList)
246 {
247 int i, j, k;
248 GEOSCoordSeq coords;
249 GEOSGeom g;
250 GEOSGeom outerRing;
251 GEOSGeom *innerRings=NULL;
252 int numInnerRings=0, *innerList;
253 GEOSContextHandle_t handle = msGetGeosContextHandle();
254
255 if(!shape || !outerList) return NULL;
256
257 /* build the outer shell */
258 coords = GEOSCoordSeq_create_r(handle,shape->line[r].numpoints, 2); /* todo handle z's */
259 if(!coords) return NULL;
260
261 for(i=0; i<shape->line[r].numpoints; i++) {
262 GEOSCoordSeq_setX_r(handle,coords, i, shape->line[r].point[i].x);
263 GEOSCoordSeq_setY_r(handle,coords, i, shape->line[r].point[i].y);
264 /* GEOSCoordSeq_setZ(coords, i, shape->line[r].point[i].z); */
265 }
266
267 outerRing = GEOSGeom_createLinearRing_r(handle,coords); /* outerRing owns the coordinates in coords */
268
269 /* build the holes */
270 innerList = msGetInnerList(shape, r, outerList);
271 for(j=0; j<shape->numlines; j++)
272 if(innerList[j] == MS_TRUE) numInnerRings++;
273
274 if(numInnerRings > 0) {
275 k = 0; /* inner ring counter */
276
277 innerRings = msSmallMalloc(numInnerRings*sizeof(GEOSGeom));
278
279 for(j=0; j<shape->numlines; j++) {
280 if(innerList[j] == MS_FALSE) continue;
281
282 coords = GEOSCoordSeq_create_r(handle,shape->line[j].numpoints, 2); /* todo handle z's */
283 if(!coords) {
284 free(innerRings);
285 free(innerList);
286 return NULL; /* todo, this will leak memory (shell + allocated holes) */
287 }
288
289 for(i=0; i<shape->line[j].numpoints; i++) {
290 GEOSCoordSeq_setX_r(handle,coords, i, shape->line[j].point[i].x);
291 GEOSCoordSeq_setY_r(handle,coords, i, shape->line[j].point[i].y);
292 /* GEOSCoordSeq_setZ(coords, i, shape->line[j].point[i].z); */
293 }
294
295 innerRings[k] = GEOSGeom_createLinearRing_r(handle,coords); /* innerRings[k] owns the coordinates in coords */
296 k++;
297 }
298 }
299
300 g = GEOSGeom_createPolygon_r(handle,outerRing, innerRings, numInnerRings);
301
302 free(innerList); /* clean up */
303 free(innerRings); /* clean up */
304
305 return g;
306 }
307
msGEOSShape2Geometry_polygon(shapeObj * shape)308 static GEOSGeom msGEOSShape2Geometry_polygon(shapeObj *shape)
309 {
310 int i, j;
311 GEOSGeom *polygons;
312 int *outerList, numOuterRings=0, lastOuterRing=0;
313 GEOSGeom g;
314 GEOSContextHandle_t handle = msGetGeosContextHandle();
315
316 outerList = msGetOuterList(shape);
317 for(i=0; i<shape->numlines; i++) {
318 if(outerList[i] == MS_TRUE) {
319 numOuterRings++;
320 lastOuterRing = i; /* save for the simple case */
321 }
322 }
323
324 if(numOuterRings == 1) {
325 g = msGEOSShape2Geometry_simplepolygon(shape, lastOuterRing, outerList);
326 } else { /* a true multipolygon */
327 polygons = msSmallMalloc(numOuterRings*sizeof(GEOSGeom));
328
329 j = 0; /* part counter */
330 for(i=0; i<shape->numlines; i++) {
331 if(outerList[i] == MS_FALSE) continue;
332 polygons[j] = msGEOSShape2Geometry_simplepolygon(shape, i, outerList); /* TODO: account for NULL return values */
333 j++;
334 }
335
336 g = GEOSGeom_createCollection_r(handle,GEOS_MULTIPOLYGON, polygons, numOuterRings);
337 free(polygons);
338 }
339
340 free(outerList);
341 return g;
342 }
343
msGEOSShape2Geometry(shapeObj * shape)344 GEOSGeom msGEOSShape2Geometry(shapeObj *shape)
345 {
346 if(!shape)
347 return NULL; /* a NULL shape generates a NULL geometry */
348
349 switch(shape->type) {
350 case MS_SHAPE_POINT:
351 if(shape->numlines == 0 || shape->line[0].numpoints == 0) /* not enough info for a point */
352 return NULL;
353
354 if(shape->line[0].numpoints == 1) /* simple point */
355 return msGEOSShape2Geometry_point(&(shape->line[0].point[0]));
356 else /* multi-point */
357 return msGEOSShape2Geometry_multipoint(&(shape->line[0]));
358 break;
359 case MS_SHAPE_LINE:
360 if(shape->numlines == 0 || shape->line[0].numpoints < 2) /* not enough info for a line */
361 return NULL;
362
363 if(shape->numlines == 1) /* simple line */
364 return msGEOSShape2Geometry_line(&(shape->line[0]));
365 else /* multi-line */
366 return msGEOSShape2Geometry_multiline(shape);
367 break;
368 case MS_SHAPE_POLYGON:
369 if(shape->numlines == 0 || shape->line[0].numpoints < 4) /* not enough info for a polygon (first=last) */
370 return NULL;
371
372 return msGEOSShape2Geometry_polygon(shape); /* simple and multipolygon cases are addressed */
373 break;
374 default:
375 break;
376 }
377
378 return NULL; /* should not get here */
379 }
380
msGEOSGeometry2Shape_point(GEOSGeom g)381 static shapeObj *msGEOSGeometry2Shape_point(GEOSGeom g)
382 {
383 GEOSCoordSeq coords;
384 shapeObj *shape=NULL;
385 GEOSContextHandle_t handle = msGetGeosContextHandle();
386
387 if(!g) return NULL;
388
389 shape = (shapeObj *) malloc(sizeof(shapeObj));
390 msInitShape(shape);
391
392 shape->type = MS_SHAPE_POINT;
393 shape->line = (lineObj *) malloc(sizeof(lineObj));
394 shape->numlines = 1;
395 shape->line[0].point = (pointObj *) malloc(sizeof(pointObj));
396 shape->line[0].numpoints = 1;
397 shape->geometry = (GEOSGeom) g;
398
399 coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,g);
400
401 GEOSCoordSeq_getX_r(handle,coords, 0, &(shape->line[0].point[0].x));
402 GEOSCoordSeq_getY_r(handle,coords, 0, &(shape->line[0].point[0].y));
403 /* GEOSCoordSeq_getZ(coords, 0, &(shape->line[0].point[0].z)); */
404
405 shape->bounds.minx = shape->bounds.maxx = shape->line[0].point[0].x;
406 shape->bounds.miny = shape->bounds.maxy = shape->line[0].point[0].y;
407
408 return shape;
409 }
410
msGEOSGeometry2Shape_multipoint(GEOSGeom g)411 static shapeObj *msGEOSGeometry2Shape_multipoint(GEOSGeom g)
412 {
413 int i;
414 int numPoints;
415 GEOSCoordSeq coords;
416 GEOSGeom point;
417
418 shapeObj *shape=NULL;
419 GEOSContextHandle_t handle = msGetGeosContextHandle();
420
421 if(!g) return NULL;
422 numPoints = GEOSGetNumGeometries_r(handle,g); /* each geometry has 1 point */
423
424 shape = (shapeObj *) malloc(sizeof(shapeObj));
425 msInitShape(shape);
426
427 shape->type = MS_SHAPE_POINT;
428 shape->line = (lineObj *) malloc(sizeof(lineObj));
429 shape->numlines = 1;
430 shape->line[0].point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
431 shape->line[0].numpoints = numPoints;
432 shape->geometry = (GEOSGeom) g;
433
434 for(i=0; i<numPoints; i++) {
435 point = (GEOSGeom) GEOSGetGeometryN_r(handle,g, i);
436 coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,point);
437
438 GEOSCoordSeq_getX_r(handle,coords, 0, &(shape->line[0].point[i].x));
439 GEOSCoordSeq_getY_r(handle,coords, 0, &(shape->line[0].point[i].y));
440 /* GEOSCoordSeq_getZ(coords, 0, &(shape->line[0].point[i].z)); */
441 }
442
443 msComputeBounds(shape);
444
445 return shape;
446 }
447
msGEOSGeometry2Shape_line(GEOSGeom g)448 static shapeObj *msGEOSGeometry2Shape_line(GEOSGeom g)
449 {
450 shapeObj *shape=NULL;
451 GEOSContextHandle_t handle = msGetGeosContextHandle();
452
453 int i;
454 int numPoints;
455 GEOSCoordSeq coords;
456
457 if(!g) return NULL;
458 numPoints = GEOSGetNumCoordinates_r(handle,g);
459 coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,g);
460
461 shape = (shapeObj *) malloc(sizeof(shapeObj));
462 msInitShape(shape);
463
464 shape->type = MS_SHAPE_LINE;
465 shape->line = (lineObj *) malloc(sizeof(lineObj));
466 shape->numlines = 1;
467 shape->line[0].point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
468 shape->line[0].numpoints = numPoints;
469 shape->geometry = (GEOSGeom) g;
470
471 for(i=0; i<numPoints; i++) {
472 GEOSCoordSeq_getX_r(handle,coords, i, &(shape->line[0].point[i].x));
473 GEOSCoordSeq_getY_r(handle,coords, i, &(shape->line[0].point[i].y));
474 /* GEOSCoordSeq_getZ(coords, i, &(shape->line[0].point[i].z)); */
475 }
476
477 msComputeBounds(shape);
478
479 return shape;
480 }
481
msGEOSGeometry2Shape_multiline(GEOSGeom g)482 static shapeObj *msGEOSGeometry2Shape_multiline(GEOSGeom g)
483 {
484 int i, j;
485 int numPoints, numLines;
486 GEOSCoordSeq coords;
487 GEOSGeom lineString;
488
489 shapeObj *shape=NULL;
490 lineObj line;
491 GEOSContextHandle_t handle = msGetGeosContextHandle();
492
493 if(!g) return NULL;
494 numLines = GEOSGetNumGeometries_r(handle,g);
495
496 shape = (shapeObj *) malloc(sizeof(shapeObj));
497 msInitShape(shape);
498
499 shape->type = MS_SHAPE_LINE;
500 shape->geometry = (GEOSGeom) g;
501
502 for(j=0; j<numLines; j++) {
503 lineString = (GEOSGeom) GEOSGetGeometryN_r(handle,g, j);
504 numPoints = GEOSGetNumCoordinates_r(handle,lineString);
505 coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,lineString);
506
507 line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
508 line.numpoints = numPoints;
509
510 for(i=0; i<numPoints; i++) {
511 GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
512 GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
513 /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
514 }
515
516 msAddLineDirectly(shape, &line);
517 }
518
519 msComputeBounds(shape);
520
521 return shape;
522 }
523
msGEOSGeometry2Shape_polygon(GEOSGeom g)524 static shapeObj *msGEOSGeometry2Shape_polygon(GEOSGeom g)
525 {
526 shapeObj *shape=NULL;
527 lineObj line;
528 int numPoints, numRings;
529 int i, j;
530
531 GEOSCoordSeq coords;
532 GEOSGeom ring;
533 GEOSContextHandle_t handle = msGetGeosContextHandle();
534
535 if(!g) return NULL;
536
537 shape = (shapeObj *) malloc(sizeof(shapeObj));
538 msInitShape(shape);
539 shape->type = MS_SHAPE_POLYGON;
540 shape->geometry = (GEOSGeom) g;
541
542 /* exterior ring */
543 ring = (GEOSGeom) GEOSGetExteriorRing_r(handle,g);
544 numPoints = GEOSGetNumCoordinates_r(handle,ring);
545 coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,ring);
546
547 line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
548 line.numpoints = numPoints;
549
550 for(i=0; i<numPoints; i++) {
551 GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
552 GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
553 /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
554 }
555 msAddLineDirectly(shape, &line);
556
557 /* interior rings */
558 numRings = GEOSGetNumInteriorRings_r(handle,g);
559 for(j=0; j<numRings; j++) {
560 ring = (GEOSGeom) GEOSGetInteriorRingN_r(handle,g, j);
561 if(GEOSisRing_r(handle,ring) != 1) continue; /* skip it */
562
563 numPoints = GEOSGetNumCoordinates_r(handle,ring);
564 coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,ring);
565
566 line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
567 line.numpoints = numPoints;
568
569 for(i=0; i<numPoints; i++) {
570 GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
571 GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
572 /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
573 }
574 msAddLineDirectly(shape, &line);
575 }
576
577 msComputeBounds(shape);
578
579 return shape;
580 }
581
msGEOSGeometry2Shape_multipolygon(GEOSGeom g)582 static shapeObj *msGEOSGeometry2Shape_multipolygon(GEOSGeom g)
583 {
584 int i, j, k;
585 shapeObj *shape=NULL;
586 lineObj line;
587 int numPoints, numRings, numPolygons;
588
589 GEOSCoordSeq coords;
590 GEOSGeom polygon, ring;
591 GEOSContextHandle_t handle = msGetGeosContextHandle();
592
593 if(!g) return NULL;
594 numPolygons = GEOSGetNumGeometries_r(handle,g);
595
596 shape = (shapeObj *) malloc(sizeof(shapeObj));
597 msInitShape(shape);
598 shape->type = MS_SHAPE_POLYGON;
599 shape->geometry = (GEOSGeom) g;
600
601 for(k=0; k<numPolygons; k++) { /* for each polygon */
602 polygon = (GEOSGeom) GEOSGetGeometryN_r(handle,g, k);
603
604 /* exterior ring */
605 ring = (GEOSGeom) GEOSGetExteriorRing_r(handle,polygon);
606 numPoints = GEOSGetNumCoordinates_r(handle,ring);
607 coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,ring);
608
609 line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
610 line.numpoints = numPoints;
611
612 for(i=0; i<numPoints; i++) {
613 GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
614 GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
615 /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
616 }
617 msAddLineDirectly(shape, &line);
618
619 /* interior rings */
620 numRings = GEOSGetNumInteriorRings_r(handle,polygon);
621
622 for(j=0; j<numRings; j++) {
623 ring = (GEOSGeom) GEOSGetInteriorRingN_r(handle,polygon, j);
624 if(GEOSisRing_r(handle,ring) != 1) continue; /* skip it */
625
626 numPoints = GEOSGetNumCoordinates_r(handle,ring);
627 coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,ring);
628
629 line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
630 line.numpoints = numPoints;
631
632 for(i=0; i<numPoints; i++) {
633 GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
634 GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
635 /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
636 }
637 msAddLineDirectly(shape, &line);
638 }
639 } /* next polygon */
640
641 msComputeBounds(shape);
642
643 return shape;
644 }
645
msGEOSGeometry2Shape(GEOSGeom g)646 shapeObj *msGEOSGeometry2Shape(GEOSGeom g)
647 {
648 int type;
649 GEOSContextHandle_t handle = msGetGeosContextHandle();
650
651 if(!g)
652 return NULL; /* a NULL geometry generates a NULL shape */
653
654 type = GEOSGeomTypeId_r(handle,g);
655 switch(type) {
656 case GEOS_POINT:
657 return msGEOSGeometry2Shape_point(g);
658 break;
659 case GEOS_MULTIPOINT:
660 return msGEOSGeometry2Shape_multipoint(g);
661 break;
662 case GEOS_LINESTRING:
663 return msGEOSGeometry2Shape_line(g);
664 break;
665 case GEOS_MULTILINESTRING:
666 return msGEOSGeometry2Shape_multiline(g);
667 break;
668 case GEOS_POLYGON:
669 return msGEOSGeometry2Shape_polygon(g);
670 break;
671 case GEOS_MULTIPOLYGON:
672 return msGEOSGeometry2Shape_multipolygon(g);
673 break;
674 case GEOS_GEOMETRYCOLLECTION:
675 if (!GEOSisEmpty_r(handle,g))
676 {
677 int i, j, numGeoms;
678 shapeObj* shape;
679
680 numGeoms = GEOSGetNumGeometries_r(handle,g);
681
682 shape = (shapeObj *) malloc(sizeof(shapeObj));
683 msInitShape(shape);
684 shape->type = MS_SHAPE_LINE;
685 shape->geometry = (GEOSGeom) g;
686
687 numGeoms = GEOSGetNumGeometries_r(handle,g);
688 for(i = 0; i < numGeoms; i++) { /* for each geometry */
689 shapeObj* shape2 = msGEOSGeometry2Shape((GEOSGeom)GEOSGetGeometryN_r(handle,g, i));
690 if (shape2) {
691 for (j = 0; j < shape2->numlines; j++)
692 msAddLineDirectly(shape, &shape2->line[j]);
693 shape2->numlines = 0;
694 shape2->geometry = NULL; /* not owned */
695 msFreeShape(shape2);
696 }
697 }
698 msComputeBounds(shape);
699 return shape;
700 }
701 break;
702 default:
703 msSetError(MS_GEOSERR, "Unsupported GEOS geometry type (%d).", "msGEOSGeometry2Shape()", type);
704 }
705 return NULL;
706 }
707 #endif
708
709 /*
710 ** Maintenence functions exposed to MapServer/MapScript.
711 */
712
msGEOSFreeGeometry(shapeObj * shape)713 void msGEOSFreeGeometry(shapeObj *shape)
714 {
715 #ifdef USE_GEOS
716 GEOSGeom g=NULL;
717 GEOSContextHandle_t handle = msGetGeosContextHandle();
718
719 if(!shape || !shape->geometry)
720 return;
721
722 g = (GEOSGeom) shape->geometry;
723 GEOSGeom_destroy_r(handle,g);
724 shape->geometry = NULL;
725 #else
726 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSFreeGEOSGeom()");
727 return;
728 #endif
729 }
730
731 /*
732 ** WKT input and output functions
733 */
msGEOSShapeFromWKT(const char * wkt)734 shapeObj *msGEOSShapeFromWKT(const char *wkt)
735 {
736 #ifdef USE_GEOS
737 GEOSGeom g;
738 GEOSContextHandle_t handle = msGetGeosContextHandle();
739
740 if(!wkt)
741 return NULL;
742
743 g = GEOSGeomFromWKT_r(handle,wkt);
744 if(!g) {
745 msSetError(MS_GEOSERR, "Error reading WKT geometry \"%s\".", "msGEOSShapeFromWKT()", wkt);
746 return NULL;
747 } else {
748 return msGEOSGeometry2Shape(g);
749 }
750
751 #else
752 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSShapeFromWKT()");
753 return NULL;
754 #endif
755 }
756
757 /* Return should be freed with msGEOSFreeWKT */
msGEOSShapeToWKT(shapeObj * shape)758 char *msGEOSShapeToWKT(shapeObj *shape)
759 {
760 #ifdef USE_GEOS
761 GEOSGeom g;
762 GEOSContextHandle_t handle = msGetGeosContextHandle();
763
764 if(!shape)
765 return NULL;
766
767 /* if we have a geometry, we should update it*/
768 msGEOSFreeGeometry(shape);
769
770 shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
771 g = (GEOSGeom) shape->geometry;
772 if(!g) return NULL;
773
774 return GEOSGeomToWKT_r(handle,g);
775 #else
776 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSShapeToWKT()");
777 return NULL;
778 #endif
779 }
780
msGEOSFreeWKT(char * pszGEOSWKT)781 void msGEOSFreeWKT(char* pszGEOSWKT)
782 {
783 #ifdef USE_GEOS
784 GEOSContextHandle_t handle = msGetGeosContextHandle();
785 #if GEOS_VERSION_MAJOR > 3 || (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 2)
786 GEOSFree_r(handle,pszGEOSWKT);
787 #endif
788 #else
789 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSFreeWKT()");
790 #endif
791 }
792
msGEOSOffsetCurve(shapeObj * p,double offset)793 shapeObj *msGEOSOffsetCurve(shapeObj *p, double offset) {
794 #if defined USE_GEOS && (GEOS_VERSION_MAJOR > 3 || (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 3))
795 GEOSGeom g1, g2;
796 GEOSContextHandle_t handle = msGetGeosContextHandle();
797
798 if(!p)
799 return NULL;
800
801 if(!p->geometry) /* if no geometry for the shape then build one */
802 p->geometry = (GEOSGeom) msGEOSShape2Geometry(p);
803
804 g1 = (GEOSGeom) p->geometry;
805 if(!g1) return NULL;
806
807 if (GEOSGeomTypeId_r(handle,g1) == GEOS_MULTILINESTRING)
808 {
809 GEOSGeom *lines = malloc(p->numlines*sizeof(GEOSGeom));
810 if (!lines) return NULL;
811 for(int i=0; i<p->numlines; i++)
812 {
813 lines[i] = GEOSOffsetCurve_r(handle, GEOSGetGeometryN_r(handle,g1,i),
814 offset, 4, GEOSBUF_JOIN_MITRE, fabs(offset*1.5));
815 }
816 g2 = GEOSGeom_createCollection_r(handle,GEOS_MULTILINESTRING, lines, p->numlines);
817 free(lines);
818 }
819 else
820 {
821 g2 = GEOSOffsetCurve_r(handle,g1, offset, 4, GEOSBUF_JOIN_MITRE, fabs(offset*1.5));
822 }
823 return msGEOSGeometry2Shape(g2);
824 #else
825 msSetError(MS_GEOSERR, "GEOS Offset Curve support is not available.", "msGEOSingleSidedBuffer()");
826 return NULL;
827 #endif
828 }
829
830 /*
831 ** Analytical functions exposed to MapServer/MapScript.
832 */
833
msGEOSBuffer(shapeObj * shape,double width)834 shapeObj *msGEOSBuffer(shapeObj *shape, double width)
835 {
836 #ifdef USE_GEOS
837 GEOSGeom g1, g2;
838 GEOSContextHandle_t handle = msGetGeosContextHandle();
839
840 if(!shape)
841 return NULL;
842
843 if(!shape->geometry) /* if no geometry for the shape then build one */
844 shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
845
846 g1 = (GEOSGeom) shape->geometry;
847 if(!g1) return NULL;
848
849 g2 = GEOSBuffer_r(handle,g1, width, 30);
850 return msGEOSGeometry2Shape(g2);
851 #else
852 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSBuffer()");
853 return NULL;
854 #endif
855 }
856
msGEOSSimplify(shapeObj * shape,double tolerance)857 shapeObj *msGEOSSimplify(shapeObj *shape, double tolerance)
858 {
859 #ifdef USE_GEOS
860 GEOSGeom g1, g2;
861 GEOSContextHandle_t handle = msGetGeosContextHandle();
862
863 if(!shape)
864 return NULL;
865
866 if(!shape->geometry) /* if no geometry for the shape then build one */
867 shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
868
869 g1 = (GEOSGeom) shape->geometry;
870 if(!g1) return NULL;
871
872 g2 = GEOSSimplify_r(handle,g1, tolerance);
873 return msGEOSGeometry2Shape(g2);
874 #else
875 msSetError(MS_GEOSERR, "GEOS Simplify support is not available.", "msGEOSSimplify()");
876 return NULL;
877 #endif
878 }
879
msGEOSTopologyPreservingSimplify(shapeObj * shape,double tolerance)880 shapeObj *msGEOSTopologyPreservingSimplify(shapeObj *shape, double tolerance)
881 {
882 #ifdef USE_GEOS
883 GEOSGeom g1, g2;
884 GEOSContextHandle_t handle = msGetGeosContextHandle();
885
886 if(!shape)
887 return NULL;
888
889 if(!shape->geometry) /* if no geometry for the shape then build one */
890 shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
891
892 g1 = (GEOSGeom) shape->geometry;
893 if(!g1) return NULL;
894
895 g2 = GEOSTopologyPreserveSimplify_r(handle,g1, tolerance);
896 return msGEOSGeometry2Shape(g2);
897 #else
898 msSetError(MS_GEOSERR, "GEOS Simplify support is not available.", "msGEOSTopologyPreservingSimplify()");
899 return NULL;
900 #endif
901 }
902
msGEOSConvexHull(shapeObj * shape)903 shapeObj *msGEOSConvexHull(shapeObj *shape)
904 {
905 #ifdef USE_GEOS
906 GEOSGeom g1, g2;
907 GEOSContextHandle_t handle = msGetGeosContextHandle();
908
909 if(!shape) return NULL;
910
911 if(!shape->geometry) /* if no geometry for the shape then build one */
912 shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
913 g1 = (GEOSGeom) shape->geometry;
914 if(!g1) return NULL;
915
916 g2 = GEOSConvexHull_r(handle,g1);
917 return msGEOSGeometry2Shape(g2);
918 #else
919 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSConvexHull()");
920 return NULL;
921 #endif
922 }
923
msGEOSBoundary(shapeObj * shape)924 shapeObj *msGEOSBoundary(shapeObj *shape)
925 {
926 #ifdef USE_GEOS
927 GEOSGeom g1, g2;
928 GEOSContextHandle_t handle = msGetGeosContextHandle();
929
930 if(!shape) return NULL;
931
932 if(!shape->geometry) /* if no geometry for the shape then build one */
933 shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
934 g1 = (GEOSGeom) shape->geometry;
935 if(!g1) return NULL;
936
937 g2 = GEOSBoundary_r(handle,g1);
938 return msGEOSGeometry2Shape(g2);
939 #else
940 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSBoundary()");
941 return NULL;
942 #endif
943 }
944
msGEOSGetCentroid(shapeObj * shape)945 pointObj *msGEOSGetCentroid(shapeObj *shape)
946 {
947 #ifdef USE_GEOS
948 GEOSGeom g1, g2;
949 GEOSCoordSeq coords;
950 pointObj *point;
951 GEOSContextHandle_t handle = msGetGeosContextHandle();
952
953 if(!shape) return NULL;
954
955 if(!shape->geometry) /* if no geometry for the shape then build one */
956 shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
957 g1 = (GEOSGeom) shape->geometry;
958 if(!g1) return NULL;
959
960 g2 = GEOSGetCentroid_r(handle,g1);
961 if (!g2) return NULL;
962
963 point = (pointObj *) malloc(sizeof(pointObj));
964
965 coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,g2);
966
967 GEOSCoordSeq_getX_r(handle,coords, 0, &(point->x));
968 GEOSCoordSeq_getY_r(handle,coords, 0, &(point->y));
969 /* GEOSCoordSeq_getZ(coords, 0, &(point->z)); */
970
971 GEOSGeom_destroy_r(handle, g2);
972
973 return point;
974 #else
975 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSGetCentroid()");
976 return NULL;
977 #endif
978 }
979
msGEOSUnion(shapeObj * shape1,shapeObj * shape2)980 shapeObj *msGEOSUnion(shapeObj *shape1, shapeObj *shape2)
981 {
982 #ifdef USE_GEOS
983 GEOSGeom g1, g2, g3;
984 GEOSContextHandle_t handle = msGetGeosContextHandle();
985
986 if(!shape1 || !shape2)
987 return NULL;
988
989 if(!shape1->geometry) /* if no geometry for the shape then build one */
990 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
991 g1 = (GEOSGeom) shape1->geometry;
992 if(!g1) return NULL;
993
994 if(!shape2->geometry) /* if no geometry for the shape then build one */
995 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
996 g2 = (GEOSGeom) shape2->geometry;
997 if(!g2) return NULL;
998
999 g3 = GEOSUnion_r(handle,g1, g2);
1000 return msGEOSGeometry2Shape(g3);
1001 #else
1002 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSUnion()");
1003 return NULL;
1004 #endif
1005 }
1006
msGEOSIntersection(shapeObj * shape1,shapeObj * shape2)1007 shapeObj *msGEOSIntersection(shapeObj *shape1, shapeObj *shape2)
1008 {
1009 #ifdef USE_GEOS
1010 GEOSGeom g1, g2, g3;
1011 GEOSContextHandle_t handle = msGetGeosContextHandle();
1012
1013 if(!shape1 || !shape2)
1014 return NULL;
1015
1016 if(!shape1->geometry) /* if no geometry for the shape then build one */
1017 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1018 g1 = (GEOSGeom) shape1->geometry;
1019 if(!g1) return NULL;
1020
1021 if(!shape2->geometry) /* if no geometry for the shape then build one */
1022 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1023 g2 = (GEOSGeom) shape2->geometry;
1024 if(!g2) return NULL;
1025
1026 g3 = GEOSIntersection_r(handle,g1, g2);
1027 return msGEOSGeometry2Shape(g3);
1028 #else
1029 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSIntersection()");
1030 return NULL;
1031 #endif
1032 }
1033
msGEOSDifference(shapeObj * shape1,shapeObj * shape2)1034 shapeObj *msGEOSDifference(shapeObj *shape1, shapeObj *shape2)
1035 {
1036 #ifdef USE_GEOS
1037 GEOSGeom g1, g2, g3;
1038 GEOSContextHandle_t handle = msGetGeosContextHandle();
1039
1040 if(!shape1 || !shape2)
1041 return NULL;
1042
1043 if(!shape1->geometry) /* if no geometry for the shape then build one */
1044 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1045 g1 = (GEOSGeom) shape1->geometry;
1046 if(!g1) return NULL;
1047
1048 if(!shape2->geometry) /* if no geometry for the shape then build one */
1049 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1050 g2 = (GEOSGeom) shape2->geometry;
1051 if(!g2) return NULL;
1052
1053 g3 = GEOSDifference_r(handle,g1, g2);
1054 return msGEOSGeometry2Shape(g3);
1055 #else
1056 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSDifference()");
1057 return NULL;
1058 #endif
1059 }
1060
msGEOSSymDifference(shapeObj * shape1,shapeObj * shape2)1061 shapeObj *msGEOSSymDifference(shapeObj *shape1, shapeObj *shape2)
1062 {
1063 #ifdef USE_GEOS
1064 GEOSGeom g1, g2, g3;
1065 GEOSContextHandle_t handle = msGetGeosContextHandle();
1066
1067 if(!shape1 || !shape2)
1068 return NULL;
1069
1070 if(!shape1->geometry) /* if no geometry for the shape then build one */
1071 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1072 g1 = (GEOSGeom) shape1->geometry;
1073 if(!g1) return NULL;
1074
1075 if(!shape2->geometry) /* if no geometry for the shape then build one */
1076 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1077 g2 = (GEOSGeom) shape2->geometry;
1078 if(!g2) return NULL;
1079
1080 g3 = GEOSSymDifference_r(handle,g1, g2);
1081 return msGEOSGeometry2Shape(g3);
1082 #else
1083 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSSymDifference()");
1084 return NULL;
1085 #endif
1086 }
1087
1088 /*
1089 ** Binary predicates exposed to MapServer/MapScript
1090 */
1091
1092 /*
1093 ** Does shape1 contain shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1094 */
msGEOSContains(shapeObj * shape1,shapeObj * shape2)1095 int msGEOSContains(shapeObj *shape1, shapeObj *shape2)
1096 {
1097 #ifdef USE_GEOS
1098 GEOSGeom g1, g2;
1099 int result;
1100 GEOSContextHandle_t handle = msGetGeosContextHandle();
1101
1102 if(!shape1 || !shape2)
1103 return -1;
1104
1105 if(!shape1->geometry) /* if no geometry for shape1 then build one */
1106 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1107 g1 = shape1->geometry;
1108 if(!g1) return -1;
1109
1110 if(!shape2->geometry) /* if no geometry for shape2 then build one */
1111 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1112 g2 = shape2->geometry;
1113 if(!g2) return -1;
1114
1115 result = GEOSContains_r(handle,g1, g2);
1116 return ((result==2) ? -1 : result);
1117 #else
1118 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSContains()");
1119 return -1;
1120 #endif
1121 }
1122
1123 /*
1124 ** Does shape1 overlap shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1125 */
msGEOSOverlaps(shapeObj * shape1,shapeObj * shape2)1126 int msGEOSOverlaps(shapeObj *shape1, shapeObj *shape2)
1127 {
1128 #ifdef USE_GEOS
1129 GEOSGeom g1, g2;
1130 int result;
1131 GEOSContextHandle_t handle = msGetGeosContextHandle();
1132
1133 if(!shape1 || !shape2)
1134 return -1;
1135
1136 if(!shape1->geometry) /* if no geometry for shape1 then build one */
1137 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1138 g1 = shape1->geometry;
1139 if(!g1) return -1;
1140
1141 if(!shape2->geometry) /* if no geometry for shape2 then build one */
1142 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1143 g2 = shape2->geometry;
1144 if(!g2) return -1;
1145
1146 result = GEOSOverlaps_r(handle,g1, g2);
1147 return ((result==2) ? -1 : result);
1148 #else
1149 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSOverlaps()");
1150 return -1;
1151 #endif
1152 }
1153
1154 /*
1155 ** Is shape1 within shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1156 */
msGEOSWithin(shapeObj * shape1,shapeObj * shape2)1157 int msGEOSWithin(shapeObj *shape1, shapeObj *shape2)
1158 {
1159 #ifdef USE_GEOS
1160 GEOSGeom g1, g2;
1161 int result;
1162 GEOSContextHandle_t handle = msGetGeosContextHandle();
1163
1164 if(!shape1 || !shape2)
1165 return -1;
1166
1167 if(!shape1->geometry) /* if no geometry for shape1 then build one */
1168 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1169 g1 = shape1->geometry;
1170 if(!g1) return -1;
1171
1172 if(!shape2->geometry) /* if no geometry for shape2 then build one */
1173 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1174 g2 = shape2->geometry;
1175 if(!g2) return -1;
1176
1177 result = GEOSWithin_r(handle,g1, g2);
1178 return ((result==2) ? -1 : result);
1179 #else
1180 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSWithin()");
1181 return -1;
1182 #endif
1183 }
1184
1185 /*
1186 ** Does shape1 cross shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1187 */
msGEOSCrosses(shapeObj * shape1,shapeObj * shape2)1188 int msGEOSCrosses(shapeObj *shape1, shapeObj *shape2)
1189 {
1190 #ifdef USE_GEOS
1191 GEOSGeom g1, g2;
1192 int result;
1193 GEOSContextHandle_t handle = msGetGeosContextHandle();
1194
1195 if(!shape1 || !shape2)
1196 return -1;
1197
1198 if(!shape1->geometry) /* if no geometry for shape1 then build one */
1199 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1200 g1 = shape1->geometry;
1201 if(!g1) return -1;
1202
1203 if(!shape2->geometry) /* if no geometry for shape2 then build one */
1204 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1205 g2 = shape2->geometry;
1206 if(!g2) return -1;
1207
1208 result = GEOSCrosses_r(handle,g1, g2);
1209 return ((result==2) ? -1 : result);
1210 #else
1211 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSCrosses()");
1212 return -1;
1213 #endif
1214 }
1215
1216 /*
1217 ** Does shape1 intersect shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1218 */
msGEOSIntersects(shapeObj * shape1,shapeObj * shape2)1219 int msGEOSIntersects(shapeObj *shape1, shapeObj *shape2)
1220 {
1221 #ifdef USE_GEOS
1222 GEOSGeom g1, g2;
1223 int result;
1224 GEOSContextHandle_t handle = msGetGeosContextHandle();
1225
1226 if(!shape1 || !shape2)
1227 return -1;
1228
1229 if(!shape1->geometry) /* if no geometry for shape1 then build one */
1230 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1231 g1 = (GEOSGeom) shape1->geometry;
1232 if(!g1) return -1;
1233
1234 if(!shape2->geometry) /* if no geometry for shape2 then build one */
1235 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1236 g2 = (GEOSGeom) shape2->geometry;
1237 if(!g2) return -1;
1238
1239 result = GEOSIntersects_r(handle,g1, g2);
1240 return ((result==2) ? -1 : result);
1241 #else
1242 if(!shape1 || !shape2)
1243 return -1;
1244
1245 switch(shape1->type) { /* todo: deal with point shapes */
1246 case(MS_SHAPE_LINE):
1247 switch(shape2->type) {
1248 case(MS_SHAPE_LINE):
1249 return msIntersectPolylines(shape1, shape2);
1250 case(MS_SHAPE_POLYGON):
1251 return msIntersectPolylinePolygon(shape1, shape2);
1252 }
1253 break;
1254 case(MS_SHAPE_POLYGON):
1255 switch(shape2->type) {
1256 case(MS_SHAPE_LINE):
1257 return msIntersectPolylinePolygon(shape2, shape1);
1258 case(MS_SHAPE_POLYGON):
1259 return msIntersectPolygons(shape1, shape2);
1260 }
1261 break;
1262 }
1263
1264 return -1;
1265 #endif
1266 }
1267
1268 /*
1269 ** Does shape1 touch shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1270 */
msGEOSTouches(shapeObj * shape1,shapeObj * shape2)1271 int msGEOSTouches(shapeObj *shape1, shapeObj *shape2)
1272 {
1273 #ifdef USE_GEOS
1274 GEOSGeom g1, g2;
1275 int result;
1276 GEOSContextHandle_t handle = msGetGeosContextHandle();
1277
1278 if(!shape1 || !shape2)
1279 return -1;
1280
1281 if(!shape1->geometry) /* if no geometry for shape1 then build one */
1282 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1283 g1 = (GEOSGeom) shape1->geometry;
1284 if(!g1) return -1;
1285
1286 if(!shape2->geometry) /* if no geometry for shape2 then build one */
1287 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1288 g2 = (GEOSGeom) shape2->geometry;
1289 if(!g2) return -1;
1290
1291 result = GEOSTouches_r(handle,g1, g2);
1292 return ((result==2) ? -1 : result);
1293 #else
1294 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSTouches()");
1295 return -1;
1296 #endif
1297 }
1298
1299 /*
1300 ** Does shape1 equal shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1301 */
msGEOSEquals(shapeObj * shape1,shapeObj * shape2)1302 int msGEOSEquals(shapeObj *shape1, shapeObj *shape2)
1303 {
1304 #ifdef USE_GEOS
1305 GEOSGeom g1, g2;
1306 int result;
1307 GEOSContextHandle_t handle = msGetGeosContextHandle();
1308
1309 if(!shape1 || !shape2)
1310 return -1;
1311
1312 if(!shape1->geometry) /* if no geometry for shape1 then build one */
1313 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1314 g1 = (GEOSGeom) shape1->geometry;
1315 if(!g1) return -1;
1316
1317 if(!shape2->geometry) /* if no geometry for shape2 then build one */
1318 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1319 g2 = (GEOSGeom) shape2->geometry;
1320 if(!g2) return -1;
1321
1322 result = GEOSEquals_r(handle,g1, g2);
1323 return ((result==2) ? -1 : result);
1324 #else
1325 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSEquals()");
1326 return -1;
1327 #endif
1328 }
1329
1330 /*
1331 ** Are shape1 and shape2 disjoint, returns MS_TRUE/MS_FALSE or -1 for an error.
1332 */
msGEOSDisjoint(shapeObj * shape1,shapeObj * shape2)1333 int msGEOSDisjoint(shapeObj *shape1, shapeObj *shape2)
1334 {
1335 #ifdef USE_GEOS
1336 GEOSGeom g1, g2;
1337 int result;
1338 GEOSContextHandle_t handle = msGetGeosContextHandle();
1339
1340 if(!shape1 || !shape2)
1341 return -1;
1342
1343 if(!shape1->geometry) /* if no geometry for shape1 then build one */
1344 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1345 g1 = (GEOSGeom) shape1->geometry;
1346 if(!g1) return -1;
1347
1348 if(!shape2->geometry) /* if no geometry for shape2 then build one */
1349 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1350 g2 = (GEOSGeom) shape2->geometry;
1351 if(!g2) return -1;
1352
1353 result = GEOSDisjoint_r(handle,g1, g2);
1354 return ((result==2) ? -1 : result);
1355 #else
1356 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSDisjoint()");
1357 return -1;
1358 #endif
1359 }
1360
1361 /*
1362 ** Useful misc. functions that return -1 on failure.
1363 */
msGEOSArea(shapeObj * shape)1364 double msGEOSArea(shapeObj *shape)
1365 {
1366 #if defined(USE_GEOS) && defined(GEOS_CAPI_VERSION_MAJOR) && defined(GEOS_CAPI_VERSION_MINOR) && (GEOS_CAPI_VERSION_MAJOR > 1 || GEOS_CAPI_VERSION_MINOR >= 1)
1367 GEOSGeom g;
1368 double area;
1369 int result;
1370 GEOSContextHandle_t handle = msGetGeosContextHandle();
1371
1372 if(!shape) return -1;
1373
1374 if(!shape->geometry) /* if no geometry for the shape then build one */
1375 shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
1376 g = (GEOSGeom) shape->geometry;
1377 if(!g) return -1;
1378
1379 result = GEOSArea_r(handle,g, &area);
1380 return ((result==0) ? -1 : area);
1381 #elif defined(USE_GEOS)
1382 msSetError(MS_GEOSERR, "GEOS support enabled, but old version lacks GEOSArea().", "msGEOSArea()");
1383 return -1;
1384 #else
1385 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSArea()");
1386 return -1;
1387 #endif
1388 }
1389
msGEOSLength(shapeObj * shape)1390 double msGEOSLength(shapeObj *shape)
1391 {
1392 #if defined(USE_GEOS) && defined(GEOS_CAPI_VERSION_MAJOR) && defined(GEOS_CAPI_VERSION_MINOR) && (GEOS_CAPI_VERSION_MAJOR > 1 || GEOS_CAPI_VERSION_MINOR >= 1)
1393
1394 GEOSGeom g;
1395 double length;
1396 int result;
1397 GEOSContextHandle_t handle = msGetGeosContextHandle();
1398
1399 if(!shape) return -1;
1400
1401 if(!shape->geometry) /* if no geometry for the shape then build one */
1402 shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
1403 g = (GEOSGeom) shape->geometry;
1404 if(!g) return -1;
1405
1406 result = GEOSLength_r(handle,g, &length);
1407 return ((result==0) ? -1 : length);
1408 #elif defined(USE_GEOS)
1409 msSetError(MS_GEOSERR, "GEOS support enabled, but old version lacks GEOSLength().", "msGEOSLength()");
1410 return -1;
1411 #else
1412 msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSLength()");
1413 return -1;
1414 #endif
1415 }
1416
msGEOSDistance(shapeObj * shape1,shapeObj * shape2)1417 double msGEOSDistance(shapeObj *shape1, shapeObj *shape2)
1418 {
1419 #ifdef USE_GEOS
1420 GEOSGeom g1, g2;
1421 double distance;
1422 int result;
1423 GEOSContextHandle_t handle = msGetGeosContextHandle();
1424
1425 if(!shape1 || !shape2)
1426 return -1;
1427
1428 if(!shape1->geometry) /* if no geometry for shape1 then build one */
1429 shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1430 g1 = (GEOSGeom) shape1->geometry;
1431 if(!g1) return -1;
1432
1433 if(!shape2->geometry) /* if no geometry for shape2 then build one */
1434 shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1435 g2 = (GEOSGeom) shape2->geometry;
1436 if(!g2) return -1;
1437
1438 result = GEOSDistance_r(handle,g1, g2, &distance);
1439 return ((result==0) ? -1 : distance);
1440 #else
1441 return msDistanceShapeToShape(shape1, shape2); /* fall back on brute force method (for MapScript) */
1442 #endif
1443 }
1444