-
/**********************************************************************
* $Id$
*
*
**********************************************************************
* $Log$
+ * Revision 1.5 2003/08/01 23:58:08 dblasby
+ * added the functionality to convert GEOS->PostGIS geometries. Added those geos
+ * functions to postgis.
+ *
* Revision 1.4 2003/07/01 18:30:55 pramsey
* Added CVS revision headers.
*
*
**********************************************************************/
+
+//--------------------------------------------------------------------------
+//
#ifdef USE_GEOS
#include "postgres.h"
extern char GEOSrelateOverlaps(Geometry *g1, Geometry*g2);
extern char *GEOSasText(Geometry *g1);
+extern char GEOSisEmpty(Geometry *g1);
+extern char *GEOSGeometryType(Geometry *g1);
extern void GEOSdeleteChar(char *a);
extern char *throw_exception(Geometry *g);
+extern Geometry *GEOSIntersection(Geometry *g1, Geometry *g2);
+
+
+extern POINT3D *GEOSGetCoordinate(Geometry *g1);
+extern POINT3D *GEOSGetCoordinates(Geometry *g1);
+extern int GEOSGetNumCoordinate(Geometry *g1);
+extern Geometry *GEOSGetGeometryN(Geometry *g1, int n);
+extern Geometry *GEOSGetExteriorRing(Geometry *g1);
+extern Geometry *GEOSGetInteriorRingN(Geometry *g1,int n);
+extern int GEOSGetNumInteriorRings(Geometry *g1);
+extern int GEOSGetSRID(Geometry *g1);
+extern int GEOSGetNumGeometries(Geometry *g1);
+
+extern Geometry *GEOSBuffer(Geometry *g1,double width);
+extern Geometry *GEOSConvexHull(Geometry *g1);
+extern Geometry *GEOSDifference(Geometry *g1,Geometry *g2);
+extern Geometry *GEOSBoundary(Geometry *g1);
+extern Geometry *GEOSSymDifference(Geometry *g1,Geometry *g2);
+extern Geometry *GEOSUnion(Geometry *g1,Geometry *g2);
+
Datum isvalid(PG_FUNCTION_ARGS);
+Datum buffer(PG_FUNCTION_ARGS);
+Datum intersection(PG_FUNCTION_ARGS);
+Datum convexhull(PG_FUNCTION_ARGS);
+Datum difference(PG_FUNCTION_ARGS);
+Datum boundary(PG_FUNCTION_ARGS);
+Datum symdifference(PG_FUNCTION_ARGS);
+Datum geomunion(PG_FUNCTION_ARGS);
+
Geometry *POSTGIS2GEOS(GEOMETRY *g);
void errorIfGeometryCollection(GEOMETRY *g1, GEOMETRY *g2);
+GEOMETRY *GEOS2POSTGIS(Geometry *g, char want3d );
+
+POLYGON3D *PolyFromGeometry(Geometry *g, int *size);
+LINE3D *LineFromGeometry(Geometry *g, int *size);
+
+//-----------------------------------------------
// return a GEOS Geometry from a POSTGIS GEOMETRY
+//----------------------------------------------
+
+
+//select geomunion('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))');
+PG_FUNCTION_INFO_V1(geomunion);
+Datum geomunion(PG_FUNCTION_ARGS)
+{
+ GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+ Geometry *g1,*g2,*g3;
+ GEOMETRY *result;
+ char empty;
+
+ initGEOS(MAXIMUM_ALIGNOF);
+
+ g1 = POSTGIS2GEOS(geom1 );
+ g2 = POSTGIS2GEOS(geom2 );
+ g3 = GEOSUnion(g1,g2);
+
+ if (g3 == NULL)
+ {
+ elog(ERROR,"GEOS union() threw an error!");
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+ empty = GEOSisEmpty(g3);
+ if (empty ==2)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS union() threw an error (couldnt test empty on result)!");
+ PG_RETURN_NULL(); //never get here
+ }
+ if (empty)
+ {
+ PG_RETURN_NULL();
+ }
+
+// elog(NOTICE,"result: %s", GEOSasText(g3) ) ;
+
+ result = GEOS2POSTGIS(g3, geom1->is3d || geom2->is3d);
+ if (result == NULL)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS union() threw an error (result postgis geometry formation)!");
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+
+ PG_RETURN_POINTER(result);
+}
+
+
+// select symdifference('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))');
+PG_FUNCTION_INFO_V1(symdifference);
+Datum symdifference(PG_FUNCTION_ARGS)
+{
+ GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+ Geometry *g1,*g2,*g3;
+ GEOMETRY *result;
+ char empty;
+
+ initGEOS(MAXIMUM_ALIGNOF);
+
+ g1 = POSTGIS2GEOS(geom1 );
+ g2 = POSTGIS2GEOS(geom2 );
+ g3 = GEOSSymDifference(g1,g2);
+
+ if (g3 == NULL)
+ {
+ elog(ERROR,"GEOS symdifference() threw an error!");
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+ empty = GEOSisEmpty(g3);
+ if (empty ==2)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS symdifference() threw an error (couldnt test empty on result)!");
+ PG_RETURN_NULL(); //never get here
+ }
+ if (empty)
+ {
+ PG_RETURN_NULL();
+ }
+
+// elog(NOTICE,"result: %s", GEOSasText(g3) ) ;
+
+ result = GEOS2POSTGIS(g3, geom1->is3d || geom2->is3d);
+ if (result == NULL)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS symdifference() threw an error (result postgis geometry formation)!");
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+
+ PG_RETURN_POINTER(result);
+}
+
+
+PG_FUNCTION_INFO_V1(boundary);
+Datum boundary(PG_FUNCTION_ARGS)
+{
+ GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+
+ Geometry *g1,*g3;
+ GEOMETRY *result;
+ char empty;
+
+ initGEOS(MAXIMUM_ALIGNOF);
+
+ g1 = POSTGIS2GEOS(geom1 );
+ g3 = GEOSBoundary(g1);
+
+ if (g3 == NULL)
+ {
+ elog(ERROR,"GEOS bounary() threw an error!");
+ GEOSdeleteGeometry(g1);
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+ empty = GEOSisEmpty(g3);
+ if (empty ==2)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS bounary() threw an error (couldnt test empty on result)!");
+ PG_RETURN_NULL(); //never get here
+ }
+ if (empty)
+ {
+ PG_RETURN_NULL();
+ }
+
+// elog(NOTICE,"result: %s", GEOSasText(g3) ) ;
+
+ result = GEOS2POSTGIS(g3, geom1->is3d);
+ if (result == NULL)
+ {
+ GEOSdeleteGeometry(g1);
+
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS bounary() threw an error (result postgis geometry formation)!");
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g3);
+
+ PG_RETURN_POINTER(result);
+}
+
+PG_FUNCTION_INFO_V1(convexhull);
+Datum convexhull(PG_FUNCTION_ARGS)
+{
+ GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ char empty;
+ Geometry *g1,*g3;
+ GEOMETRY *result;
+
+ initGEOS(MAXIMUM_ALIGNOF);
+
+ g1 = POSTGIS2GEOS(geom1 );
+ g3 = GEOSConvexHull(g1);
+
+
+ if (g3 == NULL)
+ {
+ elog(ERROR,"GEOS convexhull() threw an error!");
+ GEOSdeleteGeometry(g1);
+ PG_RETURN_NULL(); //never get here
+ }
+
+ empty = GEOSisEmpty(g3);
+ if (empty ==2)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS convexhull() threw an error (couldnt test empty on result)!");
+ PG_RETURN_NULL(); //never get here
+ }
+ if (empty)
+ {
+ PG_RETURN_NULL();
+ }
+
+// elog(NOTICE,"result: %s", GEOSasText(g3) ) ;
+
+ result = GEOS2POSTGIS(g3, geom1->is3d);
+ if (result == NULL)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS convexhull() threw an error (result postgis geometry formation)!");
+ PG_RETURN_NULL(); //never get here
+ }
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g3);
+
+ PG_RETURN_POINTER(result);
+
+}
+
+PG_FUNCTION_INFO_V1(buffer);
+Datum buffer(PG_FUNCTION_ARGS)
+{
+ GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ double size = PG_GETARG_FLOAT8(1);
+ char empty;
+ Geometry *g1,*g3;
+ GEOMETRY *result;
+
+ initGEOS(MAXIMUM_ALIGNOF);
+
+ g1 = POSTGIS2GEOS(geom1 );
+ g3 = GEOSBuffer(g1,size);
+
+
+ if (g3 == NULL)
+ {
+ elog(ERROR,"GEOS buffer() threw an error!");
+ GEOSdeleteGeometry(g1);
+ PG_RETURN_NULL(); //never get here
+ }
+
+ empty = GEOSisEmpty(g3);
+ if (empty ==2)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS buffer() threw an error (couldnt test empty on result)!");
+ PG_RETURN_NULL(); //never get here
+ }
+ if (empty)
+ {
+ PG_RETURN_NULL();
+ }
+
+// elog(NOTICE,"result: %s", GEOSasText(g3) ) ;
+
+ result = GEOS2POSTGIS(g3, geom1->is3d);
+ if (result == NULL)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS buffer() threw an error (result postgis geometry formation)!");
+ PG_RETURN_NULL(); //never get here
+ }
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g3);
+
+ PG_RETURN_POINTER(result);
+
+}
+
+
+//select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))');
+//select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))');
+//select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','POLYGON((25 5, 35 5, 35 7, 25 7, 25 5))');
+//select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','POINT(5 5)');
+//select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','LINESTRING(5 5, 10 10)');
+//select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','MULTIPOINT(5 5, 7 7, 9 9, 10 10, 11 11)');
+//select intersection('POLYGON(( 0 0, 10 0, 10 10, 0 10, 0 0))','POLYGON((5 5, 15 5, 15 7, 5 7, 5 5 ),(6 6,6.5 6, 6.5 6.5,6 6.5,6 6))');
+
+
+////select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','MULTIPOINT(5 5, 7 7, 9 9, 10 10, 11 11)');
+// select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','MULTILINESTRING((5 5, 10 10),(1 1, 2 2) )');
+//select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','MULTILINESTRING((5 5, 10 10),(1 1, 2 2) )');
+//select intersection('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','MULTIPOLYGON(((5 5, 15 5, 15 7, 5 7, 5 5)),((1 1,1 2,2 2,1 2, 1 1)))');
+PG_FUNCTION_INFO_V1(intersection);
+Datum intersection(PG_FUNCTION_ARGS)
+{
+ GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+ Geometry *g1,*g2,*g3;
+ GEOMETRY *result;
+ char empty;
+
+ initGEOS(MAXIMUM_ALIGNOF);
+
+ g1 = POSTGIS2GEOS(geom1 );
+ g2 = POSTGIS2GEOS(geom2 );
+ g3 = GEOSIntersection(g1,g2);
+
+ if (g3 == NULL)
+ {
+ elog(ERROR,"GEOS Intersection() threw an error!");
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+ empty = GEOSisEmpty(g3);
+ if (empty ==2)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS Intersection() threw an error (couldnt test empty on result)!");
+ PG_RETURN_NULL(); //never get here
+ }
+ if (empty)
+ {
+ PG_RETURN_NULL();
+ }
+
+// elog(NOTICE,"result: %s", GEOSasText(g3) ) ;
+
+ result = GEOS2POSTGIS(g3, geom1->is3d || geom2->is3d);
+ if (result == NULL)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS Intersection() threw an error (result postgis geometry formation)!");
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+
+ PG_RETURN_POINTER(result);
+}
+
+//select difference('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))','POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))');
+PG_FUNCTION_INFO_V1(difference);
+Datum difference(PG_FUNCTION_ARGS)
+{
+ GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+ Geometry *g1,*g2,*g3;
+ GEOMETRY *result;
+ char empty;
+
+ initGEOS(MAXIMUM_ALIGNOF);
+
+ g1 = POSTGIS2GEOS(geom1 );
+ g2 = POSTGIS2GEOS(geom2 );
+ g3 = GEOSDifference(g1,g2);
+
+ if (g3 == NULL)
+ {
+ elog(ERROR,"GEOS difference() threw an error!");
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+ empty = GEOSisEmpty(g3);
+ if (empty ==2)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS difference() threw an error (couldnt test empty on result)!");
+ PG_RETURN_NULL(); //never get here
+ }
+ if (empty)
+ {
+ PG_RETURN_NULL();
+ }
+
+// elog(NOTICE,"result: %s", GEOSasText(g3) ) ;
+
+ result = GEOS2POSTGIS(g3, geom1->is3d || geom2->is3d);
+ if (result == NULL)
+ {
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+ elog(ERROR,"GEOS difference() threw an error (result postgis geometry formation)!");
+ PG_RETURN_NULL(); //never get here
+ }
+
+
+
+ GEOSdeleteGeometry(g1);
+ GEOSdeleteGeometry(g2);
+ GEOSdeleteGeometry(g3);
+
+ PG_RETURN_POINTER(result);
+}
+
+
+
+//----------------------------------------------
+
void errorIfGeometryCollection(GEOMETRY *g1, GEOMETRY *g2)
}
+POLYGON3D *PolyFromGeometry(Geometry *g, int *size)
+{
+
+ int ninteriorrings =GEOSGetNumInteriorRings(g);
+ POINT3D *pts;
+ int *pts_per_ring;
+ int t;
+ POLYGON3D *poly;
+ int npoints;
+
+
+ npoints = GEOSGetNumCoordinate(g);
+ pts = GEOSGetCoordinates(g);
+ if (npoints <3)
+ {
+ GEOSdeleteChar( (char*) pts);
+ return NULL;
+ }
+
+ pts_per_ring = palloc(sizeof(int) * (ninteriorrings+1));
+ pts_per_ring[0] = GEOSGetNumCoordinate(GEOSGetExteriorRing(g));
+
+ for (t=0;t<ninteriorrings;t++)
+ {
+ pts_per_ring[t+1] = GEOSGetNumCoordinate(GEOSGetInteriorRingN(g,t));
+ }
+
+
+
+ poly = make_polygon( ninteriorrings+1, pts_per_ring,
+ pts, GEOSGetNumCoordinate(g), size);
+
+
+
+ GEOSdeleteChar( (char*) pts);
+ return poly;
+
+}
+
+
+
+LINE3D *LineFromGeometry(Geometry *g,int *size)
+{
+ POINT3D *pts = GEOSGetCoordinates(g);
+
+ LINE3D *line;
+ int npoints = GEOSGetNumCoordinate(g);
+
+ if (npoints <2)
+ {
+ GEOSdeleteChar( (char*) pts);
+ return NULL;
+ }
+
+ line = make_line(npoints, pts, size);
+ GEOSdeleteChar( (char*) pts);
+ return line;
+}
+
+
+
+
+GEOMETRY *GEOS2POSTGIS(Geometry *g,char want3d)
+{
+ char *type = GEOSGeometryType(g) ;
+ GEOMETRY *result = NULL;
+
+ if (strcmp(type,"Point") ==0 )
+ {
+ POINT3D *pt = GEOSGetCoordinate(g);
+
+ result = make_oneobj_geometry(sizeof(POINT3D),
+ (char *) pt,
+ POINTTYPE, want3d, GEOSGetSRID(g),1.0, 0.0, 0.0
+ );
+
+
+ GEOSdeleteChar( (char*) pt);
+ return result;
+ }
+ else if (strcmp(type,"LineString")==0)
+ {
+ LINE3D *line;
+ int size;
+
+ line = LineFromGeometry(g,&size);
+ if (line == NULL)
+ return NULL;
+
+
+ result = make_oneobj_geometry(size,
+ (char *) line,
+ LINETYPE, want3d, GEOSGetSRID(g),1.0, 0.0, 0.0
+ );
+ return result;
+ }
+ else if (strcmp(type,"Polygon")==0)
+ {
+
+ int size;
+ POLYGON3D *poly;
+
+
+
+ poly = PolyFromGeometry(g,&size);
+ if (poly == NULL)
+ return NULL;
+
+ result = make_oneobj_geometry(size,
+ (char *) poly,
+ POLYGONTYPE, want3d, GEOSGetSRID(g),1.0, 0.0, 0.0
+ );
+ return result;
+ }
+ else if (strcmp(type,"MultiPoint")==0)
+ {
+ int ngeoms,t;
+ POINT3D *pts;
+ GEOMETRY *g_new,*g_old;
+
+ g_new = NULL;
+
+ ngeoms = GEOSGetNumGeometries(g);
+ if (ngeoms ==0)
+ {
+ return NULL;
+ }
+ pts = GEOSGetCoordinates(g);
+
+ g_old = make_oneobj_geometry(sizeof(POINT3D),
+ (char *) pts,
+ MULTIPOINTTYPE, want3d, GEOSGetSRID(g),1.0, 0.0, 0.0
+ );
+
+ for (t=1;t<ngeoms;t++)
+ {
+ g_new= add_to_geometry(g_old,sizeof(POINT3D), (char*)&pts[t], POINTTYPE);
+ pfree(g_old);
+ g_old =g_new;
+ }
+
+ GEOSdeleteChar( (char*) pts);
+ return g_new;
+ }
+ else if (strcmp(type,"MultiLineString")==0)
+ {
+ int ngeoms,t,size;
+ LINE3D *line;
+ GEOMETRY *g_old;
+
+
+ ngeoms = GEOSGetNumGeometries(g);
+ if (ngeoms ==0)
+ {
+ return NULL;
+ }
+ for (t=0;t<ngeoms;t++)
+ {
+ if (t==0)
+ {
+ line = LineFromGeometry(GEOSGetGeometryN(g,0) ,&size);
+ result = make_oneobj_geometry(size,
+ (char *) line,
+ MULTILINETYPE, want3d, GEOSGetSRID(g),1.0, 0.0, 0.0
+ );
+ }
+ else
+ {
+ line = LineFromGeometry(GEOSGetGeometryN(g,t) ,&size);
+ g_old = result;
+ result = add_to_geometry(g_old,size, (char*) line, LINETYPE);
+ pfree(g_old);
+ }
+ }
+ return result;
+
+ }
+ else if (strcmp(type,"MultiPolygon")==0)
+ {
+ int ngeoms,t,size;
+ POLYGON3D *poly;
+ GEOMETRY *g_old;
+
+
+ ngeoms = GEOSGetNumGeometries(g);
+ if (ngeoms ==0)
+ {
+ return NULL;
+ }
+ for (t=0;t<ngeoms;t++)
+ {
+ if (t==0)
+ {
+ poly = PolyFromGeometry(GEOSGetGeometryN(g,0) ,&size);
+ result = make_oneobj_geometry(size,
+ (char *) poly,
+ MULTIPOLYGONTYPE, want3d, GEOSGetSRID(g),1.0, 0.0, 0.0
+ );
+ }
+ else
+ {
+ poly = PolyFromGeometry(GEOSGetGeometryN(g,t) ,&size);
+ g_old = result;
+ result = add_to_geometry(g_old,size, (char*) poly, POLYGONTYPE);
+ pfree(g_old);
+ }
+ }
+ return result;
+ }
+ else if (strcmp(type,"GeometryCollection")==0)
+ {
+ //this is more difficult because GEOS allows GCs of GCs
+ int ngeoms = GEOSGetNumGeometries(g);
+ GEOMETRY *geom, *g2, *r;
+ int t;
+
+ if (ngeoms ==0)
+ {
+ return NULL;
+ }
+ if (ngeoms == 1)
+ {
+ return GEOS2POSTGIS(GEOSGetGeometryN(g,0) , want3d); // short cut!
+ }
+ geom = GEOS2POSTGIS(GEOSGetGeometryN(g,0) , want3d);
+ for (t=1;t<ngeoms;t++)
+ {
+ g2 = GEOS2POSTGIS(GEOSGetGeometryN(g,t) , want3d);
+ r = geom;
+ geom = (GEOMETRY *)
+ DatumGetPointer(
+ DirectFunctionCall2(collector,
+ PointerGetDatum(geom),PointerGetDatum(g2)
+ )
+ );
+ pfree(r);
+ pfree(g2);
+ }
+ return geom;
+ }
+ return NULL;
+}
+
+
+
//BBOXONLYTYPE -> returns as a 2d polygon
Geometry *POSTGIS2GEOS(GEOMETRY *g)
{
#include "io.h"
//#include "opRelate.h"
+using namespace geos;
+
//WARNING THIS *MUST* BE SET CORRECTLY.
extern "C" char *GEOSasText(Geometry *g1);
+extern "C" char GEOSisEmpty(Geometry *g1);
+extern "C" char *GEOSGeometryType(Geometry *g1);
+
extern "C" char *throw_exception(Geometry *g);
+extern "C" Geometry *GEOSIntersection(Geometry *g1,Geometry *g1);
+extern "C" Geometry *GEOSBuffer(Geometry *g1,double width);
+extern "C" Geometry *GEOSConvexHull(Geometry *g1);
+extern "C" Geometry *GEOSDifference(Geometry *g1,Geometry *g2);
+extern "C" Geometry *GEOSBoundary(Geometry *g1);
+extern "C" Geometry *GEOSSymDifference(Geometry *g1,Geometry *g2);
+extern "C" Geometry *GEOSUnion(Geometry *g1,Geometry *g2);
+
+
+extern "C" POINT3D *GEOSGetCoordinate(Geometry *g1);
+extern "C" POINT3D *GEOSGetCoordinates(Geometry *g1);
+extern "C" int GEOSGetNumCoordinate(Geometry *g1);
+extern "C" Geometry *GEOSGetGeometryN(Geometry *g1, int n);
+extern "C" Geometry *GEOSGetExteriorRing(Geometry *g1);
+extern "C" Geometry *GEOSGetInteriorRingN(Geometry *g1, int n);
+extern "C" int GEOSGetNumInteriorRings(Geometry *g1);
+extern "C" int GEOSGetSRID(Geometry *g1);
+extern "C" int GEOSGetNumGeometries(Geometry *g1);
+
//###########################################################
c.y = line->points[t].y;
c.z = DoubleNotANumber;
coords->setAt( c ,t);
+
}
}
}
}
+char GEOSisEmpty(Geometry *g1)
+{
+ try
+ {
+ return g1->isEmpty();
+ }
+ catch (...)
+ {
+ return 2;
+ }
+}
+
+char *GEOSGeometryType(Geometry *g1)
+{
+ try
+ {
+ string s = g1->getGeometryType();
+
+
+ char *result;
+ result = (char*) malloc( s.length() + 1);
+ strcpy(result, s.c_str() );
+ return result;
+ }
+ catch (...)
+ {
+ return NULL;
+ }
+}
+
+
+
+
+//-------------------------------------------------------------------
+// GEOS functions that return geometries
+//-------------------------------------------------------------------
+
+Geometry *GEOSIntersection(Geometry *g1,Geometry *g2)
+{
+ try
+ {
+ Geometry *g3 = g1->intersection(g2);
+ return g3;
+ }
+ catch (...)
+ {
+ return NULL;
+ }
+}
+
+Geometry *GEOSBuffer(Geometry *g1,double width)
+{
+ try
+ {
+ Geometry *g3 = g1->buffer(width);
+ return g3;
+ }
+ catch (...)
+ {
+ return NULL;
+ }
+}
+
+Geometry *GEOSConvexHull(Geometry *g1)
+{
+ try
+ {
+ Geometry *g3 = g1->convexHull();
+ return g3;
+ }
+ catch (...)
+ {
+ return NULL;
+ }
+}
+
+Geometry *GEOSDifference(Geometry *g1,Geometry *g2)
+{
+ try
+ {
+ Geometry *g3 = g1->difference(g2);
+ return g3;
+ }
+ catch (...)
+ {
+ return NULL;
+ }
+}
+
+Geometry *GEOSBoundary(Geometry *g1)
+{
+ try
+ {
+ Geometry *g3 = g1->getBoundary();
+ return g3;
+ }
+ catch (...)
+ {
+ return NULL;
+ }
+}
+
+Geometry *GEOSSymDifference(Geometry *g1,Geometry *g2)
+{
+ try
+ {
+ Geometry *g3 = g1->symDifference(g2);
+ return g3;
+ }
+ catch (...)
+ {
+ return NULL;
+ }
+}
+
+Geometry *GEOSUnion(Geometry *g1,Geometry *g2)
+{
+ try
+ {
+ Geometry *g3 = g1->Union(g2);
+ return g3;
+ }
+ catch (...)
+ {
+ return NULL;
+ }
+}
+
+
+
+
+
+
+
+
//-------------------------------------------------------------------
// memory management functions
//------------------------------------------------------------------
void GEOSdeleteChar(char *a)
{
try{
- delete a;
+ free(a);
}
catch(...)
{
// do nothing!
}
+}
+
+
+//-------------------------------------------------------------------
+//GEOS => POSTGIS conversions
+//-------------------------------------------------------------------
+
+// free the result when done!
+// g1 must be a point
+POINT3D *GEOSGetCoordinate(Geometry *g1)
+{
+ try{
+ POINT3D *result = (POINT3D*) malloc (sizeof(POINT3D));
+ Coordinate *c =g1->getCoordinate();
+
+ result->x = c->x;
+ result->y = c->y;
+ result->z = c->z;
+ return result;
+ }
+ catch(...)
+ {
+ return NULL;
+ }
+
}
+
+//must free the result when done
+// result is an array length g1->getNumCoordinates()
+// only call on linestrings or linearrings
+POINT3D *GEOSGetCoordinates(Geometry *g1)
+{
+ try {
+ int numPoints = g1->getNumPoints();
+ POINT3D *result = (POINT3D*) malloc (sizeof(POINT3D) * numPoints );
+ int t;
+ CoordinateList *cl = g1->getCoordinates();
+ Coordinate c;
+
+ for (t=0;t<numPoints;t++)
+ {
+ c =cl->getAt(t);
+
+ result[t].x = c.x;
+ result[t].y = c.y;
+ result[t].z = c.z;
+ }
+
+ return result;
+ }
+ catch(...)
+ {
+ return NULL;
+ }
+
+}
+
+
+
+int GEOSGetNumCoordinate(Geometry *g1)
+{
+ try{
+ return g1->getNumPoints();
+ }
+ catch(...)
+ {
+ return 0;
+ }
+}
+
+int GEOSGetNumInteriorRings(Geometry *g1)
+{
+ try{
+ Polygon *p = (Polygon *) g1;
+ return p->getNumInteriorRing();
+ }
+ catch(...)
+ {
+ return 0;
+ }
+}
+
+
+//only call on GCs
+int GEOSGetNumGeometries(Geometry *g1)
+{
+ try{
+ GeometryCollection *gc = (GeometryCollection *) g1;
+ return gc->getNumGeometries();
+ }
+ catch(...)
+ {
+ return 0;
+ }
+}
+
+
+//call only on GEOMETRYCOLLECTION or MULTI*
+Geometry *GEOSGetGeometryN(Geometry *g1, int n)
+{
+ try{
+ GeometryCollection *gc = (GeometryCollection *) g1;
+ return gc->getGeometryN(n);
+ }
+ catch(...)
+ {
+ return NULL;
+ }
+}
+
+
+//call only on polygon
+Geometry *GEOSGetExteriorRing(Geometry *g1)
+{
+ try{
+ Polygon *p = (Polygon *) g1;
+ return p->getExteriorRing();
+ }
+ catch(...)
+ {
+ return 0;
+ }
+}
+
+//call only on polygon
+Geometry *GEOSGetInteriorRingN(Geometry *g1,int n)
+{
+ try{
+ Polygon *p = (Polygon *) g1;
+ return p->getInteriorRingN(n);
+ }
+ catch(...)
+ {
+ return NULL;
+ }
+}
+
+
+int GEOSGetSRID(Geometry *g1)
+{
+ try{
+ return g1->getSRID();
+ }
+ catch(...)
+ {
+ return 0;
+ }
+}
+
+
projects / postgis / commitdiff