1
/******************************************************************************
4
* Purpose: MapServer-GEOS integration.
5
* Author: Steve Lime and the MapServer team.
7
******************************************************************************
8
* Copyright (c) 1996-2005 Regents of the University of Minnesota.
10
* Permission is hereby granted, free of charge, to any person obtaining a
11
* copy of this software and associated documentation files (the "Software"),
12
* to deal in the Software without restriction, including without limitation
13
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
14
* and/or sell copies of the Software, and to permit persons to whom the
15
* Software is furnished to do so, subject to the following conditions:
17
* The above copyright notice and this permission notice shall be included in
18
* all copies of this Software or works derived from this Software.
20
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26
* DEALINGS IN THE SOFTWARE.
27
******************************************************************************
29
* $Log: mapgeos.cpp,v $
30
* Revision 1.14.2.1 2006/01/18 07:14:38 sdlime
31
* Update GEOS toWKT function to produce a malloc'd pointer. Updated shape.i to produce a new object. (bug 1466)
33
* Revision 1.14 2005/10/31 04:59:43 sdlime
34
* Now catching exceptions when reading WKT strings.
36
* Revision 1.13 2005/10/30 05:05:07 sdlime
37
* Initial support for WKT via GEOS. The reader is only integrated via the map file reader, with MapScript, CGI and URL support following ASAP. (bug 1466)
39
* Revision 1.12 2005/06/14 16:03:33 dan
40
* Updated copyright date to 2005
42
* Revision 1.11 2005/06/14 04:28:44 sdlime
43
* Fixed GEOS to shapeObj for multipolgon geometries.
45
* Revision 1.10 2005/05/19 05:57:08 sdlime
46
* More interface clean up, added msGEOSFreeGeometry and updated geometry2shape code to so the new shapeObj contains a reference to the geometry used to create it.
48
* Revision 1.9 2005/05/17 13:59:57 sdlime
49
* More GEOS interface improvements (global GeometryFactory).
51
* Revision 1.8 2005/05/11 07:10:21 sdlime
52
* Finished the rest of the wrapper funtions for initial GEOS support which basically ammounts to buffer and convex hull creation from MapScript at the moment. There are likely a number of memory leaks associated with the implementation at the moment.
54
* Revision 1.7 2005/02/23 05:16:50 sdlime
55
* Added GEOS=>shape conversion for GEOS_MULTILINE geometries.
57
* Revision 1.6 2005/02/23 04:52:33 sdlime
58
* Added GEOS=>shape conversion for GEOS_MULTIPOINT geometries.
60
* Revision 1.5 2005/02/23 04:40:17 sdlime
61
* Added wrapper for creating convex hulls to GEOS support. Added to MapScript as well.
63
* Revision 1.4 2005/02/22 18:33:37 dan
64
* Fixes to build without USE_GEOS
66
* Revision 1.3 2005/02/22 07:40:27 sdlime
67
* A bunch of updates to GEOS integration. Can move many primatives between
68
* MapServer and GEOS, still need to do collections (e.g. multi-point/line/polygon).
69
* Added buffer method to mapscript (mapscript/shape.i).
71
* Revision 1.2 2004/10/28 02:23:35 frank
72
* added standard header
85
MS_CVSID("$Id: mapgeos.cpp,v 1.14.2.1 2006/01/18 07:14:38 sdlime Exp $")
90
** Global variables: only the read-only GeometryFactory.
92
GeometryFactory *gf=NULL;
94
static void msGEOSCreateGeometryFactory()
97
gf = new GeometryFactory();
101
** Translation functions
103
static Geometry *msGEOSShape2Geometry_point(pointObj *point)
106
Coordinate *c = new Coordinate(point->x, point->y);
107
Geometry *g = gf->createPoint(*c);
110
} catch (GEOSException *ge) {
111
msSetError(MS_GEOSERR, "%s", "msGEOSShape2Geometry_point()", (char *) ge->toString().c_str());
119
static Geometry *msGEOSShape2Geometry_multipoint(lineObj *multipoint)
123
vector<Geometry *> *parts=NULL;
125
parts = new vector<Geometry *>(multipoint->numpoints);
126
for(i=0; i<multipoint->numpoints; i++)
127
(*parts)[i] = msGEOSShape2Geometry_point(&(multipoint->point[i]));
129
Geometry *g = gf->createMultiPoint(parts);
132
} catch (GEOSException *ge) {
133
msSetError(MS_GEOSERR, "%s", "msGEOSShape2Geometry_multiline()", (char *) ge->toString().c_str());
141
static Geometry *msGEOSShape2Geometry_line(lineObj *line)
147
DefaultCoordinateSequence *coords = new DefaultCoordinateSequence(line->numpoints);
149
c.z = DoubleNotANumber; /* same for all points in the line (TODO: support z) */
150
for(i=0; i<line->numpoints; i++) {
151
c.x = line->point[i].x;
152
c.y = line->point[i].y;
156
Geometry *g = gf->createLineString(coords);
159
} catch (GEOSException *ge) {
160
msSetError(MS_GEOSERR, "%s", "msGEOSShape2Geometry_line()", (char *) ge->toString().c_str());
168
static Geometry *msGEOSShape2Geometry_multiline(shapeObj *multiline)
172
vector<Geometry *> *parts=NULL;
174
parts = new vector<Geometry *>(multiline->numlines);
175
for(i=0; i<multiline->numlines; i++)
176
(*parts)[i] = msGEOSShape2Geometry_line(&(multiline->line[i]));
178
Geometry *g = gf->createMultiLineString(parts);
181
} catch (GEOSException *ge) {
182
msSetError(MS_GEOSERR, "%s", "msGEOSShape2Geometry_multiline()", (char *) ge->toString().c_str());
190
static Geometry *msGEOSShape2Geometry_simplepolygon(shapeObj *shape, int r, int *outerList)
195
LinearRing *outerRing, *innerRing;
196
vector<Geometry *> *innerRings=NULL;
197
int numInnerRings=0, *innerList;
199
/* build the outer ring */
200
DefaultCoordinateSequence *coords = new DefaultCoordinateSequence(shape->line[r].numpoints);
202
c.z = DoubleNotANumber; /* same for all points in the line (TODO: support z) */
203
for(i=0; i<shape->line[r].numpoints; i++) {
204
c.x = shape->line[r].point[i].x;
205
c.y = shape->line[r].point[i].y;
209
outerRing = (LinearRing *) gf->createLinearRing(coords);
211
/* build the inner rings */
212
innerList = msGetInnerList(shape, r, outerList);
213
for(j=0; j<shape->numlines; j++)
214
if(innerList[j] == MS_TRUE) numInnerRings++;
216
/* printf("\tnumber of inner rings=%d\n", numInnerRings); */
218
if(numInnerRings > 0) {
219
k = 0; /* inner ring counter */
220
innerRings = new vector<Geometry *>(numInnerRings);
221
for(j=0; j<shape->numlines; j++) {
222
if(innerList[j] == MS_FALSE) continue;
224
DefaultCoordinateSequence *coords = new DefaultCoordinateSequence(shape->line[j].numpoints);
226
c.z = DoubleNotANumber; /* same for all points in the line (TODO: support z) */
227
for(i=0; i<shape->line[j].numpoints; i++) {
228
c.x = shape->line[j].point[i].x;
229
c.y = shape->line[j].point[i].y;
233
innerRing = (LinearRing *) gf->createLinearRing(coords);
234
(*innerRings)[k] = innerRing;
239
Geometry *g = gf->createPolygon(outerRing, innerRings);
241
free(innerList); /* clean up */
244
} catch (GEOSException *ge) {
245
msSetError(MS_GEOSERR, "%s", "msGEOSShape2Geometry_simplepolygon()", (char *) ge->toString().c_str());
253
static Geometry *msGEOSShape2Geometry_polygon(shapeObj *shape)
257
vector<Geometry *> *parts=NULL;
258
int *outerList, numOuterRings=0, lastOuterRing=0;
261
outerList = msGetOuterList(shape);
262
for(i=0; i<shape->numlines; i++) {
263
if(outerList[i] == MS_TRUE) {
265
lastOuterRing = i; /* save for the simple case */
269
/* printf("number of outer rings=%d\n", numOuterRings); */
271
if(numOuterRings == 1) {
272
g = msGEOSShape2Geometry_simplepolygon(shape, lastOuterRing, outerList);
273
} else { /* a true multipolygon */
274
parts = new vector<Geometry *>(numOuterRings);
276
j = 0; /* part counter */
277
for(i=0; i<shape->numlines; i++) {
278
if(outerList[i] == MS_FALSE) continue;
279
/* printf("working on ring, id=%d (part %d)\n", i, j); */
280
(*parts)[j] = msGEOSShape2Geometry_simplepolygon(shape, i, outerList); /* TODO: account for NULL return values */
284
g = gf->createMultiPolygon(parts);
289
} catch (GEOSException *ge) {
290
msSetError(MS_GEOSERR, "%s", "msGEOSShape2Geometry_line()", (char *) ge->toString().c_str());
298
Geometry *msGEOSShape2Geometry(shapeObj *shape)
301
return NULL; /* a NULL shape generates a NULL geometry */
303
/* if there is not an instance of a GeometeryFactory, create one */
305
msGEOSCreateGeometryFactory();
307
switch(shape->type) {
309
if(shape->numlines == 0 || shape->line[0].numpoints == 0) /* not enough info for a point */
312
if(shape->line[0].numpoints == 1) /* simple point */
313
return msGEOSShape2Geometry_point(&(shape->line[0].point[0]));
314
else /* multi-point */
315
return msGEOSShape2Geometry_multipoint(&(shape->line[0]));
318
if(shape->numlines == 0 || shape->line[0].numpoints < 2) /* not enough info for a line */
321
if(shape->numlines == 1) /* simple line */
322
return msGEOSShape2Geometry_line(&(shape->line[0]));
323
else /* multi-line */
324
return msGEOSShape2Geometry_multiline(shape);
326
case MS_SHAPE_POLYGON:
327
if(shape->numlines == 0 || shape->line[0].numpoints < 4) /* not enough info for a polygon (first=last) */
330
return msGEOSShape2Geometry_polygon(shape); /* simple and multipolygon cases are addressed */
336
return NULL; /* should not get here */
339
static shapeObj *msGEOSGeometry2Shape_point(Geometry *g)
342
const Coordinate *c = g->getCoordinate();
343
shapeObj *shape=NULL;
345
shape = (shapeObj *) malloc(sizeof(shapeObj));
348
shape->type = MS_SHAPE_POINT;
349
shape->line = (lineObj *) malloc(sizeof(lineObj));
351
shape->line[0].point = (pointObj *) malloc(sizeof(pointObj));
352
shape->line[0].numpoints = 1;
355
shape->line[0].point[0].x = c->x;
356
shape->line[0].point[0].y = c->y;
357
/* shape->line[0].point[0].z = c->z; */
360
} catch (GEOSException *ge) {
361
msSetError(MS_GEOSERR, "%s", "msGEOSGeometry2Shape_point()", (char *) ge->toString().c_str());
369
static shapeObj *msGEOSGeometry2Shape_multipoint(Geometry *g)
373
int numPoints = g->getNumPoints();
374
CoordinateSequence *coords = g->getCoordinates(); /* would be nice to have read-only access */
376
shapeObj *shape=NULL;
378
shape = (shapeObj *) malloc(sizeof(shapeObj));
381
shape->type = MS_SHAPE_POINT;
382
shape->line = (lineObj *) malloc(sizeof(lineObj));
384
shape->line[0].point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
385
shape->line[0].numpoints = numPoints;
388
for(i=0; i<numPoints; i++) {
389
c = coords->getAt(i);
391
shape->line[0].point[i].x = c.x;
392
shape->line[0].point[i].y = c.y;
393
/* shape->line[0].point[i].z = c.z; */
398
} catch (GEOSException *ge) {
399
msSetError(MS_GEOSERR, "%s", "msGEOSGeometry2Shape_multipoint()", (char *) ge->toString().c_str());
407
static shapeObj *msGEOSGeometry2Shape_line(LineString *g)
409
shapeObj *shape=NULL;
413
int numPoints = g->getNumPoints();
414
const CoordinateSequence *coords = g->getCoordinatesRO(); /* a pointer to coordinates, do not delete */
417
shape = (shapeObj *) malloc(sizeof(shapeObj));
420
shape->type = MS_SHAPE_LINE;
421
shape->line = (lineObj *) malloc(sizeof(lineObj));
423
shape->line[0].point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
424
shape->line[0].numpoints = numPoints;
427
for(i=0; i<numPoints; i++) {
428
c = coords->getAt(i);
430
shape->line[0].point[i].x = c.x;
431
shape->line[0].point[i].y = c.y;
432
/* shape->line[0].point[i].z = c.z; */
436
} catch (GEOSException *ge) {
437
msSetError(MS_GEOSERR, "%s", "msGEOSGeometry2Shape_line()", (char *) ge->toString().c_str());
445
static shapeObj *msGEOSGeometry2Shape_multiline(MultiLineString *g)
449
int numPoints, numLines = g->getNumGeometries();
450
const CoordinateSequence *coords;
451
const LineString *lineString;
454
shapeObj *shape=NULL;
457
shape = (shapeObj *) malloc(sizeof(shapeObj));
460
shape->type = MS_SHAPE_LINE;
461
shape->line = (lineObj *) malloc(sizeof(lineObj)*numLines);
462
shape->numlines = numLines;
465
for(j=0; j<numLines; j++) {
466
lineString = (LineString *) g->getGeometryN(j);
467
coords = lineString->getCoordinatesRO();
468
numPoints = lineString->getNumPoints();
470
line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
471
line.numpoints = numPoints;
473
for(i=0; i<numPoints; i++) {
474
c = coords->getAt(i);
476
line.point[i].x = c.x;
477
line.point[i].y = c.y;
478
/* line.point[i].z = c.z; */
480
msAddLine(shape, &line);
485
} catch (GEOSException *ge) {
486
msSetError(MS_GEOSERR, "%s", "msGEOSGeometry2Shape_multiline()", (char *) ge->toString().c_str());
494
static shapeObj *msGEOSGeometry2Shape_polygon(Polygon *g)
497
shapeObj *shape=NULL;
499
int numPoints, numRings;
503
const CoordinateSequence *coords;
504
const LineString *ring;
506
shape = (shapeObj *) malloc(sizeof(shapeObj));
508
shape->type = MS_SHAPE_POLYGON;
512
ring = g->getExteriorRing();
513
coords = ring->getCoordinatesRO();
514
numPoints = ring->getNumPoints();
516
line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
517
line.numpoints = numPoints;
519
for(i=0; i<numPoints; i++) {
520
c = coords->getAt(i);
522
line.point[i].x = c.x;
523
line.point[i].y = c.y;
524
/* line.point[i].z = c.z; */
526
msAddLine(shape, &line);
530
numRings = g->getNumInteriorRing();
531
for(j=0; j<numRings; j++) {
532
ring = g->getInteriorRingN(j);
533
coords = ring->getCoordinatesRO();
534
numPoints = ring->getNumPoints();
536
line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
537
line.numpoints = numPoints;
539
for(i=0; i<numPoints; i++) {
540
c = coords->getAt(i);
542
line.point[i].x = c.x;
543
line.point[i].y = c.y;
544
/* line.point[i].z = 0;
545
line.point[i].m = 0; */
547
msAddLine(shape, &line);
552
} catch (GEOSException *ge) {
553
msSetError(MS_GEOSERR, "%s", "msGEOSGeometry2Shape_polygon()", (char *) ge->toString().c_str());
561
static shapeObj *msGEOSGeometry2Shape_multipolygon(MultiPolygon *g)
565
shapeObj *shape=NULL;
567
int numPoints, numRings, numPolygons=g->getNumGeometries();
570
const CoordinateSequence *coords;
571
const LineString *ring;
572
const Polygon *polygon;
574
shape = (shapeObj *) malloc(sizeof(shapeObj));
577
shape->type = MS_SHAPE_POLYGON;
580
for(k=0; k<numPolygons; k++) { /* for each polygon */
581
polygon = (Polygon *) g->getGeometryN(k);
584
ring = polygon->getExteriorRing();
585
coords = ring->getCoordinatesRO();
586
numPoints = ring->getNumPoints();
588
line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
589
line.numpoints = numPoints;
591
for(i=0; i<numPoints; i++) {
592
c = coords->getAt(i);
594
line.point[i].x = c.x;
595
line.point[i].y = c.y;
596
/* line.point[i].z = c.z;
597
line.point[i].m = 0; */
599
msAddLine(shape, &line);
603
numRings = polygon->getNumInteriorRing();
604
for(j=0; j<numRings; j++) {
605
ring = polygon->getInteriorRingN(j);
606
coords = ring->getCoordinatesRO();
607
numPoints = ring->getNumPoints();
609
line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
610
line.numpoints = numPoints;
612
for(i=0; i<numPoints; i++) {
613
c = coords->getAt(i);
615
line.point[i].x = c.x;
616
line.point[i].y = c.y;
617
/* line.point[i].z = 0;
618
line.point[i].m = 0; */
620
msAddLine(shape, &line);
626
} catch (GEOSException *ge) {
627
msSetError(MS_GEOSERR, "%s", "msGEOSGeometry2Shape_multipolygon()", (char *) ge->toString().c_str());
635
shapeObj *msGEOSGeometry2Shape(Geometry *g)
638
return NULL; /* a NULL geometry generates a NULL shape */
641
switch(g->getGeometryTypeId()) {
643
return msGEOSGeometry2Shape_point(g);
645
case GEOS_MULTIPOINT:
646
return msGEOSGeometry2Shape_multipoint(g);
648
case GEOS_LINESTRING:
649
return msGEOSGeometry2Shape_line((LineString *)g);
651
case GEOS_MULTILINESTRING:
652
return msGEOSGeometry2Shape_multiline((MultiLineString *)g);
655
return msGEOSGeometry2Shape_polygon((Polygon *)g);
657
case GEOS_MULTIPOLYGON:
658
return msGEOSGeometry2Shape_multipolygon((MultiPolygon *)g);
661
msSetError(MS_GEOSERR, "Unsupported GEOS geometry type (%d).", "msGEOSGeometry2Shape()", g->getGeometryTypeId());
664
} catch (GEOSException *ge) {
665
msSetError(MS_GEOSERR, "%s", "msGEOSGeometry2Shape()", (char *) ge->toString().c_str());
675
** Maintenence functions exposed to MapServer/MapScript.
678
void msGEOSFreeGeometry(shapeObj *shape)
681
if(!shape || !shape->geometry || !gf)
685
Geometry *g = (Geometry *) shape->geometry;
686
gf->destroyGeometry(g);
687
} catch (GEOSException *ge) {
688
msSetError(MS_GEOSERR, "%s", "msGEOSDeleteGeometry()", (char *) ge->toString().c_str());
695
msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSFreeGeometry()");
701
** WKT input and output functions
703
shapeObj *msGEOSShapeFromWKT(const char *string)
709
/* if there is not an instance of a GeometeryFactory, create one */
711
msGEOSCreateGeometryFactory();
714
WKTReader *r = new WKTReader(gf);
715
Geometry *g = r->read(string);
716
return msGEOSGeometry2Shape(g);
717
} catch (GEOSException *ge) {
718
msSetError(MS_GEOSERR, "%s", "msGEOSShapeFromWKT()", (char *) ge->toString().c_str());
726
msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSShapeFromWKT()");
731
char *msGEOSShapeToWKT(shapeObj *shape)
739
if(!shape->geometry) /* if no geometry for the shape then build one */
740
shape->geometry = msGEOSShape2Geometry(shape);
741
Geometry *g = (Geometry *) shape->geometry;
745
wkt = strdup((char *) g->toString().c_str());
749
msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSShapeToWKT()");
755
** Analytical functions exposed to MapServer/MapScript.
758
shapeObj *msGEOSBuffer(shapeObj *shape, double width)
764
if(!shape->geometry) /* if no geometry for the shape then build one */
765
shape->geometry = msGEOSShape2Geometry(shape);
766
Geometry *g = (Geometry *) shape->geometry;
771
Geometry *bg = g->buffer(width);
772
return msGEOSGeometry2Shape(bg);
773
} catch (GEOSException *ge) {
774
msSetError(MS_GEOSERR, "%s", "msGEOSBuffer()", (char *) ge->toString().c_str());
781
msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSBuffer()");
786
shapeObj *msGEOSConvexHull(shapeObj *shape)
792
if(!shape->geometry) /* if no geometry for the shape then build one */
793
shape->geometry = msGEOSShape2Geometry(shape);
794
Geometry *g = (Geometry *) shape->geometry;
799
Geometry *bg = g->convexHull();
800
return msGEOSGeometry2Shape(bg);
801
} catch (GEOSException *ge) {
802
msSetError(MS_GEOSERR, "%s", "msGEOSConvexHull()", (char *) ge->toString().c_str());
809
msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSConvexHull()");
added
removed
mapgeos.cpp
