From: Sandro Santilli Date: Tue, 7 Sep 2004 17:04:03 +0000 (+0000) Subject: GEOS support added. X-Git-Tag: pgis_0_9_1~19 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=6cdb7f83e2b412f41229d70eea3afd23c2604a0e;p=postgis GEOS support added. git-svn-id: http://svn.osgeo.org/postgis/trunk@775 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/lwgeom/MISSING_OBJECTS b/lwgeom/MISSING_OBJECTS index c414d4cb7..e1f670bf8 100644 --- a/lwgeom/MISSING_OBJECTS +++ b/lwgeom/MISSING_OBJECTS @@ -2,40 +2,11 @@ FUNC: KEEPING FUNCTION: [line_interpolate_point(geometry, double precision)] FUNC: KEEPING FUNCTION: [segmentize(geometry, double precision)] - FUNC: KEEPING FUNCTION: [asbinary(geometry)] FUNC: KEEPING FUNCTION: [asbinary(geometry, text)] +FUNC: KEEPING FUNCTION: [centroid(geometry)] (!GEOS implementation) --- OBSOLETED ? FUNC: KEEPING FUNCTION: [geometry_size(geometry, internal)] FUNC: KEEPING FUNCTION: [optimistic_overlap(geometry, geometry, double precision)] ---- GEOS -FUNC: KEEPING FUNCTION: [unite_garray(geometry[])] -FUNC: KEEPING FUNCTION: [centroid(geometry)] (Also !GEOS implementation) -AGGREGATE: KEEPING AGGREGATE [memgeomunion(geometry)] -AGGREGATE: KEEPING AGGREGATE [geomunion(geometry)] -FUNC: KEEPING FUNCTION: [isring(geometry)] -FUNC: KEEPING FUNCTION: [pointonsurface(geometry)] -FUNC: KEEPING FUNCTION: [buffer(geometry, double precision)] -FUNC: KEEPING FUNCTION: [contains(geometry, geometry)] -FUNC: KEEPING FUNCTION: [convexhull(geometry)] -FUNC: KEEPING FUNCTION: [crosses(geometry, geometry)] -FUNC: KEEPING FUNCTION: [difference(geometry, geometry)] -FUNC: KEEPING FUNCTION: [disjoint(geometry, geometry)] -FUNC: KEEPING FUNCTION: [within(geometry, geometry)] -FUNC: KEEPING FUNCTION: [symdifference(geometry, geometry)] -FUNC: KEEPING FUNCTION: [symmetricdifference(geometry, geometry)] -FUNC: KEEPING FUNCTION: [touches(geometry, geometry)] -FUNC: KEEPING FUNCTION: [relate(geometry, geometry)] -FUNC: KEEPING FUNCTION: [relate(geometry, geometry, text)] -FUNC: KEEPING FUNCTION: [overlaps(geometry, geometry)] -FUNC: KEEPING FUNCTION: [issimple(geometry)] -FUNC: KEEPING FUNCTION: [isvalid(geometry)] -FUNC: KEEPING FUNCTION: [geomunion(geometry, geometry)] -FUNC: KEEPING FUNCTION: [geosnoop(geometry)] -FUNC: KEEPING FUNCTION: [intersection(geometry, geometry)] -FUNC: KEEPING FUNCTION: [intersects(geometry, geometry)] -FUNC: KEEPING FUNCTION: [equals(geometry, geometry)] -FUNC: KEEPING FUNCTION: [boundary(geometry)] - diff --git a/lwgeom/Makefile b/lwgeom/Makefile index 1567c14a2..8733269f2 100644 --- a/lwgeom/Makefile +++ b/lwgeom/Makefile @@ -111,7 +111,7 @@ override CFLAGS += -DPOSTGIS_SCRIPTS_VERSION='"$(SCRIPTS_VERSION)"' ifeq ($(USE_GEOS),1) override CFLAGS += -I$(GEOS_DIR)/include -DUSE_GEOS #GEOS_RULES=detect_geos_version - GEOS_WRAPPER=../postgis_geos_wrapper.o + GEOS_WRAPPER=lwgeom_geos_wrapper.o endif ifeq ($(USE_PROJ),1) diff --git a/lwgeom/TODO b/lwgeom/TODO index 32dad4c69..a0daedf96 100644 --- a/lwgeom/TODO +++ b/lwgeom/TODO @@ -1,3 +1,6 @@ +- Implement lwexploded_serialize_buf for speedup of GEOS<=>POSTGIS conversion + and more! + - box extractor functions writing to a pointer or returning an actual object should probably find a way to return NULL (for empty geometries). diff --git a/lwgeom/lwgeom_api.c b/lwgeom/lwgeom_api.c index 1510f4b3a..45582b186 100644 --- a/lwgeom/lwgeom_api.c +++ b/lwgeom/lwgeom_api.c @@ -2691,10 +2691,21 @@ pfree_exploded(LWGEOM_EXPLODED *exploded) LWGEOM_EXPLODED * lwgeom_explode(char *serialized) { - LWGEOM_INSPECTED *inspected = lwgeom_inspect(serialized); + LWGEOM_INSPECTED *inspected; LWGEOM_EXPLODED *subexploded, *result; int i; +#ifdef DEBUG + elog(NOTICE, "lwgeom_explode called"); +#endif + + inspected = lwgeom_inspect(serialized); + + +#ifdef DEBUG +elog(NOTICE, "lwgeom_explode: serialized inspected"); +#endif + result = palloc(sizeof(LWGEOM_EXPLODED)); result->points = palloc(1); result->lines = palloc(1); @@ -2703,10 +2714,6 @@ lwgeom_explode(char *serialized) result->nlines = 0; result->npolys = 0; -#ifdef DEBUG - elog(NOTICE, "lwgeom_explode called"); -#endif - if ( ! inspected->ngeometries ) { pfree(result->points); @@ -2715,9 +2722,7 @@ lwgeom_explode(char *serialized) result->SRID = -1; result->ndims = 0; pfree_inspected(inspected); -#ifdef DEBUG elog(NOTICE, "lwgeom_explode: no geometries"); -#endif return result; } @@ -2732,8 +2737,11 @@ lwgeom_explode(char *serialized) if ( type == POINTTYPE ) { +#ifdef DEBUG +elog(NOTICE, "lwgeom_explode: it's a point"); +#endif result->points = repalloc(result->points, - (result->npoints+1)*sizeof(LWPOINT *)); + (result->npoints+1)*sizeof(char *)); result->points[result->npoints] = subgeom; result->npoints++; continue; @@ -2741,8 +2749,11 @@ lwgeom_explode(char *serialized) if ( type == LINETYPE ) { +#ifdef DEBUG +elog(NOTICE, "lwgeom_explode: it's a line"); +#endif result->lines = repalloc(result->lines, - (result->nlines+1)*sizeof(LWLINE *)); + (result->nlines+1)*sizeof(char *)); result->lines[result->nlines] = subgeom; result->nlines++; continue; @@ -2750,20 +2761,27 @@ lwgeom_explode(char *serialized) if ( type == POLYGONTYPE ) { +#ifdef DEBUG +elog(NOTICE, "lwgeom_explode: it's a polygon"); +#endif result->polys = repalloc(result->polys, - (result->npolys+1)*sizeof(LWPOLY *)); + (result->npolys+1)*sizeof(char *)); result->polys[result->npolys] = subgeom; result->npolys++; continue; } #ifdef DEBUG - elog(NOTICE, "subtype is %d, recursing", type); + elog(NOTICE, "type of subgeom %d is %d, recursing", i, type); #endif // it's a multi geometry, recurse subexploded = lwgeom_explode(subgeom); +#ifdef DEBUG + elog(NOTICE, "subgeom %d, exploded: %d point, %d lines, %d polys", i, subexploded->npoints, subexploded->nlines, subexploded->npolys); +#endif + // Re-allocate adding space for new exploded geoms // (-1 because 1 was already allocated for the collection) // Copy subgeom pointers from subexploded to current @@ -2772,21 +2790,36 @@ lwgeom_explode(char *serialized) if ( subexploded->npoints ) { result->points = repalloc(result->points, - result->npoints+subexploded->npoints-1); + sizeof(char *)*(result->npoints+subexploded->npoints-1)); + if ( ! result ) + elog(ERROR, "Out of virtual memory"); + +#ifdef DEBUG + elog(NOTICE, "repalloc'ed exploded->points"); +#endif - memcpy(result->points[result->npoints], + memcpy(&(result->points[result->npoints]), subexploded->points, subexploded->npoints*sizeof(char *)); +#ifdef DEBUG + elog(NOTICE, "memcpied exploded->points"); +#endif + result->npoints += subexploded->npoints; + +#ifdef DEBUG + elog(NOTICE, "memcopied %d points from subexploded (exploded points: %d", subexploded->npoints, result->npoints); +#endif } if ( subexploded->nlines ) { result->lines = repalloc(result->lines, - result->nlines+subexploded->nlines-1); + sizeof(char *)* + (result->nlines+subexploded->nlines-1)); - memcpy(result->lines[result->nlines], + memcpy(&(result->lines[result->nlines]), subexploded->lines, subexploded->nlines*sizeof(char *)); @@ -2796,9 +2829,10 @@ lwgeom_explode(char *serialized) if ( subexploded->npolys ) { result->polys = repalloc(result->polys, - result->npolys+subexploded->npolys-1); + sizeof(char *)* + (result->npolys+subexploded->npolys-1)); - memcpy(result->polys[result->npolys], + memcpy(&(result->polys[result->npolys]), subexploded->polys, subexploded->npolys*sizeof(char *)); @@ -2812,6 +2846,10 @@ lwgeom_explode(char *serialized) pfree_inspected(inspected); +#ifdef DEBUG +elog(NOTICE, "lwgeom_explode: returning"); +#endif + return result; } diff --git a/lwgeom/lwgeom_functions_basic.c b/lwgeom/lwgeom_functions_basic.c index 4db73c3bb..691aa23ef 100644 --- a/lwgeom/lwgeom_functions_basic.c +++ b/lwgeom/lwgeom_functions_basic.c @@ -49,7 +49,6 @@ int32 lwgeom_nrings_recursive(char *serialized); void dump_lwexploded(LWGEOM_EXPLODED *exploded); - /*------------------------------------------------------------------*/ // pt_in_ring_2d(): crossing number test for a point in a polygon diff --git a/lwgeom/lwgeom_geos.c b/lwgeom/lwgeom_geos.c index 9bc4376da..ce72dbbd1 100644 --- a/lwgeom/lwgeom_geos.c +++ b/lwgeom/lwgeom_geos.c @@ -1,13 +1,111 @@ #include "postgres.h" -#include "fmgr.h" +#include "utils/array.h" +#include "utils/builtins.h" + +#include "lwgeom.h" +#include "wktparse.h" + +Datum relate_full(PG_FUNCTION_ARGS); +Datum relate_pattern(PG_FUNCTION_ARGS); +Datum disjoint(PG_FUNCTION_ARGS); +Datum touches(PG_FUNCTION_ARGS); +Datum intersects(PG_FUNCTION_ARGS); +Datum crosses(PG_FUNCTION_ARGS); +Datum within(PG_FUNCTION_ARGS); +Datum contains(PG_FUNCTION_ARGS); +Datum overlaps(PG_FUNCTION_ARGS); +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); +Datum unite_garray(PG_FUNCTION_ARGS); +Datum issimple(PG_FUNCTION_ARGS); +Datum isring(PG_FUNCTION_ARGS); +Datum geomequals(PG_FUNCTION_ARGS); +Datum pointonsurface(PG_FUNCTION_ARGS); +Datum GEOSnoop(PG_FUNCTION_ARGS); +Datum postgis_geos_version(PG_FUNCTION_ARGS); +Datum centroid(PG_FUNCTION_ARGS); + + #if USE_GEOS + /* + * Define this to have have many notices printed + * during postgis->geos and geos->postgis conversions + */ +#undef DEBUG_CONVERTER +#undef DEBUG_POSTGIS2GEOS +#undef DEBUG_GEOS2POSTGIS + +typedef struct Geometry Geometry; + +extern const char * createGEOSPoint(POINT3D *pt); +extern void initGEOS(int maxalign); +extern char *GEOSrelate(Geometry *g1, Geometry*g2); +extern char GEOSrelatePattern(Geometry *g1, Geometry*g2,char *pat); +extern char GEOSrelateDisjoint(Geometry *g1, Geometry*g2); +extern char GEOSrelateTouches(Geometry *g1, Geometry*g2); +extern char GEOSrelateIntersects(Geometry *g1, Geometry*g2); +extern char GEOSrelateCrosses(Geometry *g1, Geometry*g2); +extern char GEOSrelateWithin(Geometry *g1, Geometry*g2); +extern char GEOSrelateContains(Geometry *g1, Geometry*g2); +extern char GEOSrelateOverlaps(Geometry *g1, Geometry*g2); +extern char *GEOSasText(Geometry *g1); +extern char GEOSisEmpty(Geometry *g1); +extern char *GEOSGeometryType(Geometry *g1); +extern int GEOSGeometryTypeId(Geometry *g1); extern char *GEOSversion(); extern char *GEOSjtsport(); +extern char GEOSisvalid(Geometry *g1); +extern Geometry *GEOSIntersection(Geometry *g1, Geometry *g2); +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); +extern char GEOSequals(Geometry *g1, Geometry*g2); +extern char GEOSisSimple(Geometry *g1); +extern char GEOSisRing(Geometry *g1); +extern Geometry *GEOSpointonSurface(Geometry *g1); +extern Geometry *GEOSGetCentroid(Geometry *g); + +extern void GEOSdeleteChar(char *a); +extern void GEOSdeleteGeometry(Geometry *a); + +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 *PostGIS2GEOS_point(const LWPOINT *point); +extern Geometry *PostGIS2GEOS_linestring(const LWLINE *line); +extern Geometry *PostGIS2GEOS_polygon(const LWPOLY *polygon); +extern Geometry *PostGIS2GEOS_multipolygon(LWPOLY *const *const geoms, uint32 ngeoms, int SRID, int is3d); +extern Geometry *PostGIS2GEOS_multilinestring(LWLINE *const *const geoms, uint32 ngeoms, int SRID, int is3d); +extern Geometry *PostGIS2GEOS_multipoint(LWPOINT *const *const geoms, uint32 ngeoms, int SRID, int is3d); +extern Geometry *PostGIS2GEOS_collection(Geometry **geoms, int ngeoms, int SRID, bool is3d); void NOTICE_MESSAGE(char *msg); +LWGEOM *GEOS2POSTGIS(Geometry *geom, char want3d); +Geometry * POSTGIS2GEOS(LWGEOM *g); +void errorIfGeometryCollection(LWGEOM *g1, LWGEOM *g2); +char *PointFromGeometry(Geometry *g, char want3d); +char *LineFromGeometry(Geometry *g, char want3d); +char *PolyFromGeometry(Geometry *g, char want3d); +void addToExploded_recursive(Geometry *geom, LWGEOM_EXPLODED *exp); void NOTICE_MESSAGE(char *msg) { @@ -26,13 +124,1603 @@ Datum postgis_geos_version(PG_FUNCTION_ARGS) PG_RETURN_POINTER(result); } -#else // ! defined USE_GEOS -PG_FUNCTION_INFO_V1(postgis_geos_version); -Datum postgis_geos_version(PG_FUNCTION_ARGS) +/* + * This is the final function for union/fastunite/geomunion + * aggregate (still discussing the name). Will have + * as input an array of Geometry pointers. + * Will iteratively call GEOSUnion on the GEOS-converted + * versions of them and return PGIS-converted version back. + * Changing combination order *might* speed up performance. + * + * Geometries in the array are pfree'd as soon as possible. + * + */ +PG_FUNCTION_INFO_V1(unite_garray); +Datum unite_garray(PG_FUNCTION_ARGS) { - //elog(ERROR,"GEOSversion:: operation not implemented - compile PostGIS with GEOS support"); - PG_RETURN_NULL(); + Datum datum; + ArrayType *array; + int is3d = 0; + int nelems, i; + LWGEOM **geoms, *result, *pgis_geom; + Geometry *g1, *g2, *geos_result; +#ifdef DEBUG + static int call=1; +#endif + +#ifdef DEBUG + call++; +#endif + + datum = PG_GETARG_DATUM(0); + + /* Null array, null geometry (should be empty?) */ + if ( (Pointer *)datum == NULL ) PG_RETURN_NULL(); + + array = (ArrayType *) PG_DETOAST_DATUM(datum); + + nelems = ArrayGetNItems(ARR_NDIM(array), ARR_DIMS(array)); + +#ifdef DEBUG + elog(NOTICE, "unite_garray: number of elements: %d", nelems); +#endif + + if ( nelems == 0 ) PG_RETURN_NULL(); + + geoms = (LWGEOM **)ARR_DATA_PTR(array); + + /* One-element union is the element itself */ + if ( nelems == 1 ) PG_RETURN_POINTER(geoms[0]); + + /* Ok, we really need geos now ;) */ + initGEOS(MAXIMUM_ALIGNOF); + + if ( lwgeom_ndims(geoms[0]->type) > 2 ) is3d = 1; + geos_result = POSTGIS2GEOS(geoms[0]); + pfree(geoms[0]); + for (i=1; itype) > 2 ) || + ( lwgeom_ndims(geom2->type) > 2 ); + + Geometry *g1,*g2,*g3; + LWGEOM *result; + + initGEOS(MAXIMUM_ALIGNOF); +//elog(NOTICE,"in geomunion"); + + g1 = POSTGIS2GEOS(geom1); + g2 = POSTGIS2GEOS(geom2); + +//elog(NOTICE,"g1=%s",GEOSasText(g1)); +//elog(NOTICE,"g2=%s",GEOSasText(g2)); + g3 = GEOSUnion(g1,g2); + +//elog(NOTICE,"g3=%s",GEOSasText(g3)); + + PG_FREE_IF_COPY(geom1, 0); + PG_FREE_IF_COPY(geom2, 1); + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + + if (g3 == NULL) + { + elog(ERROR,"GEOS union() threw an error!"); + PG_RETURN_NULL(); //never get here + } + +//elog(NOTICE,"result: %s", GEOSasText(g3) ) ; + + result = GEOS2POSTGIS(g3, is3d); + + GEOSdeleteGeometry(g3); + + if (result == NULL) + { + elog(ERROR,"GEOS union() threw an error (result postgis geometry formation)!"); + PG_RETURN_NULL(); //never get here + } + + //compressType(result); // convert multi* to single item if appropriate + 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) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + Geometry *g1,*g2,*g3; + LWGEOM *result; + int is3d = ( lwgeom_ndims(geom1->type) > 2 ) || + ( lwgeom_ndims(geom2->type) > 2 ); + + 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 + } + +// elog(NOTICE,"result: %s", GEOSasText(g3) ) ; + + result = GEOS2POSTGIS(g3, 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); + + //compressType(result); // convert multi* to single item if appropriate + + PG_RETURN_POINTER(result); +} + + +PG_FUNCTION_INFO_V1(boundary); +Datum boundary(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + Geometry *g1,*g3; + LWGEOM *result; + 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 + } + +// elog(NOTICE,"result: %s", GEOSasText(g3) ) ; + + result = GEOS2POSTGIS(g3, lwgeom_ndims(geom1->type) > 2); + 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); + + //compressType(result); // convert multi* to single item if appropriate + + PG_RETURN_POINTER(result); +} + +PG_FUNCTION_INFO_V1(convexhull); +Datum convexhull(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + Geometry *g1,*g3; + LWGEOM *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 + } + + +// elog(NOTICE,"result: %s", GEOSasText(g3) ) ; + + result = GEOS2POSTGIS(g3, lwgeom_ndims(geom1->type) > 2); + 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); + + + //compressType(result); // convert multi* to single item if appropriate + + PG_RETURN_POINTER(result); + +} + +PG_FUNCTION_INFO_V1(buffer); +Datum buffer(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + double size = PG_GETARG_FLOAT8(1); + Geometry *g1,*g3; + LWGEOM *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 + } + + +// elog(NOTICE,"result: %s", GEOSasText(g3) ) ; + + result = GEOS2POSTGIS(g3, lwgeom_ndims(geom1->type) > 2); + 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); + + + //compressType(result); // convert multi* to single item if appropriate + PG_RETURN_POINTER(result); + +} + +PG_FUNCTION_INFO_V1(intersection); +Datum intersection(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + Geometry *g1,*g2,*g3; + LWGEOM *result; + int is3d = ( lwgeom_ndims(geom1->type) > 2 ) || + ( lwgeom_ndims(geom2->type) > 2 ); + + initGEOS(MAXIMUM_ALIGNOF); + +//elog(NOTICE,"intersection() START"); + + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + +//elog(NOTICE," constructed geometrys - calling geos"); + +//elog(NOTICE,"g1 = %s",GEOSasText(g1)); +//elog(NOTICE,"g2 = %s",GEOSasText(g2)); + + +//if (g1==NULL) +// elog(NOTICE,"g1 is null"); +//if (g2==NULL) +// elog(NOTICE,"g2 is null"); +//elog(NOTICE,"g2 is valid = %i",GEOSisvalid(g2)); +//elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); + + + g3 = GEOSIntersection(g1,g2); + +//elog(NOTICE," intersection finished"); + + if (g3 == NULL) + { + elog(ERROR,"GEOS Intersection() threw an error!"); + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + PG_RETURN_NULL(); //never get here + } + + + + +// elog(NOTICE,"result: %s", GEOSasText(g3) ) ; + + result = GEOS2POSTGIS(g3, 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); + + //compressType(result); // convert multi* to single item if appropriate + + 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) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + Geometry *g1,*g2,*g3; + LWGEOM *result; + int is3d = ( lwgeom_ndims(geom1->type) > 2 ) || + ( lwgeom_ndims(geom2->type) > 2 ); + + 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 + } + +// elog(NOTICE,"result: %s", GEOSasText(g3) ) ; + + result = GEOS2POSTGIS(g3, 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); + + ////compressType(result); // convert multi* to single item if appropriate + + PG_RETURN_POINTER(result); +} + + +//select pointonsurface('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'); +PG_FUNCTION_INFO_V1(pointonsurface); +Datum pointonsurface(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + Geometry *g1,*g3; + LWGEOM *result; + + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1 ); + g3 = GEOSpointonSurface(g1); + + if (g3 == NULL) + { + elog(ERROR,"GEOS pointonsurface() threw an error!"); + GEOSdeleteGeometry(g1); + PG_RETURN_NULL(); //never get here + } + + + + +// elog(NOTICE,"result: %s", GEOSasText(g3) ) ; + + result = GEOS2POSTGIS(g3, (lwgeom_ndims(geom1->type) > 2)); + if (result == NULL) + { + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g3); + elog(ERROR,"GEOS pointonsurface() threw an error (result postgis geometry formation)!"); + PG_RETURN_NULL(); //never get here + } + + + + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g3); + + // convert multi* to single item if appropriate + //compressType(result); + + PG_RETURN_POINTER(result); +} + +PG_FUNCTION_INFO_V1(centroid); +Datum centroid(PG_FUNCTION_ARGS) +{ + LWGEOM *geom, *result; + Geometry *geosgeom, *geosresult; + + geom = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + initGEOS(MAXIMUM_ALIGNOF); + + geosgeom = POSTGIS2GEOS(geom); + + geosresult = GEOSGetCentroid(geosgeom); + if ( geosresult == NULL ) + { + GEOSdeleteGeometry(geosgeom); + elog(ERROR,"GEOS getCentroid() threw an error!"); + PG_RETURN_NULL(); + } + + result = GEOS2POSTGIS(geosresult, (lwgeom_ndims(geom->type) > 2)); + if (result == NULL) + { + GEOSdeleteGeometry(geosgeom); + GEOSdeleteGeometry(geosresult); + elog(ERROR,"Error in GEOS-PGIS conversion"); + PG_RETURN_NULL(); + } + GEOSdeleteGeometry(geosgeom); + GEOSdeleteGeometry(geosresult); + + PG_RETURN_POINTER(result); +} + + + +//---------------------------------------------- + + + +void errorIfGeometryCollection(LWGEOM *g1, LWGEOM *g2) +{ + int t1 = lwgeom_getType(g1->type); + int t2 = lwgeom_getType(g2->type); + + if ( (t1 == COLLECTIONTYPE) || (t2 == COLLECTIONTYPE) ) + elog(ERROR,"Relate Operation called with a LWGEOMCOLLECTION type. This is unsupported"); +} + +PG_FUNCTION_INFO_V1(isvalid); +Datum isvalid(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + bool result; + Geometry *g1; + + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1 ); + + result = GEOSisvalid(g1); + GEOSdeleteGeometry(g1); + if (result == 2) + { + elog(ERROR,"GEOS isvalid() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + + PG_RETURN_BOOL(result); } + +// overlaps(LWGEOM g1,LWGEOM g2) +// returns if GEOS::g1->overlaps(g2) returns true +// throws an error (elog(ERROR,...)) if GEOS throws an error +PG_FUNCTION_INFO_V1(overlaps); +Datum overlaps(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + Geometry *g1,*g2; + bool result; + + errorIfGeometryCollection(geom1,geom2); + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + + result = GEOSrelateOverlaps(g1,g2); + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + if (result == 2) + { + elog(ERROR,"GEOS overlaps() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + PG_RETURN_BOOL(result); +} + + + +PG_FUNCTION_INFO_V1(contains); +Datum contains(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + Geometry *g1,*g2; + bool result; + //POINT3D *testpoint; + //POLYGON3D *poly; + + errorIfGeometryCollection(geom1,geom2); + + /* + * short-circuit 1: if geom2 bounding box is not completely inside + * geom1 bounding box we can prematurely return FALSE + */ +// DISABLED ... TODO: use only if bbox is available +#if 0 + if ( geom2->bvol.LLB.x < geom1->bvol.LLB.x ) PG_RETURN_BOOL(FALSE); + if ( geom2->bvol.URT.x > geom1->bvol.URT.x ) PG_RETURN_BOOL(FALSE); + if ( geom2->bvol.LLB.y < geom1->bvol.LLB.y ) PG_RETURN_BOOL(FALSE); + if ( geom2->bvol.URT.y > geom1->bvol.URT.y ) PG_RETURN_BOOL(FALSE); + + /* + * short-circuit 2: if geom1 is a polygon and any corner of + * geom2 bounding box is not 'within' geom1 we can prematurely + * return FALSE + */ + //if ( geom1->type == POLYGONTYPE ) + //{ + // poly = (POLYGON3D *)geom1->objData; + // testpoint.x = geom2->bvol.LLB.x; + // testpoint.y = geom2->bvol.LLB.y; + // if ( !point_within_polygon(&testpoint, poly) ) + // PG_RETURN_BOOL(FALSE); + // testpoint.x = geom2->bvol.LLB.x; + // testpoint.y = geom2->bvol.URT.y; + // if ( !point_within_polygon(&testpoint, poly) ) + // PG_RETURN_BOOL(FALSE); + // testpoint.x = geom2->bvol.URT.x; + // testpoint.y = geom2->bvol.URT.y; + // if ( !point_within_polygon(&testpoint, poly) ) + // PG_RETURN_BOOL(FALSE); + // testpoint.x = geom2->bvol.URT.x; + // testpoint.y = geom2->bvol.LLB.y; + // if ( !point_within_polygon(&testpoint, poly) ) + // PG_RETURN_BOOL(FALSE); + //} +#endif + + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + + + + result = GEOSrelateContains(g1,g2); + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + + if (result == 2) + { + elog(ERROR,"GEOS contains() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + + + PG_RETURN_BOOL(result); +} + + +PG_FUNCTION_INFO_V1(within); +Datum within(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + Geometry *g1,*g2; + bool result; + //POINT3D testpoint; + //POLYGON3D *poly; + + errorIfGeometryCollection(geom1,geom2); + + +// DISABLED ... TODO: use only if bbox is available +#if 0 + /* + * short-circuit 1: if geom1 bounding box is not completely inside + * geom2 bounding box we can prematurely return FALSE + */ + if ( geom1->bvol.LLB.x < geom2->bvol.LLB.x ) PG_RETURN_BOOL(FALSE); + if ( geom1->bvol.URT.x > geom2->bvol.URT.x ) PG_RETURN_BOOL(FALSE); + if ( geom1->bvol.LLB.y < geom2->bvol.LLB.y ) PG_RETURN_BOOL(FALSE); + if ( geom1->bvol.URT.y > geom2->bvol.URT.y ) PG_RETURN_BOOL(FALSE); + + /* + * short-circuit 2: if geom2 is a polygon and any corner of + * geom1 bounding box is not 'within' geom2 we can prematurely + * return FALSE + */ + //if ( geom2->type == POLYGONTYPE ) + //{ + // poly = (POLYGON3D *)geom2->objData; + // testpoint.x = geom1->bvol.LLB.x; + // testpoint.y = geom1->bvol.LLB.y; + // if ( !point_within_polygon(&testpoint, poly) ) + // PG_RETURN_BOOL(FALSE); + // testpoint.x = geom1->bvol.LLB.x; + // testpoint.y = geom1->bvol.URT.y; + // if ( !point_within_polygon(&testpoint, poly) ) + // PG_RETURN_BOOL(FALSE); + // testpoint.x = geom1->bvol.URT.x; + // testpoint.y = geom1->bvol.URT.y; + // if ( !point_within_polygon(&testpoint, poly) ) + // PG_RETURN_BOOL(FALSE); + // testpoint.x = geom1->bvol.URT.x; + // testpoint.y = geom1->bvol.LLB.y; + // if ( !point_within_polygon(&testpoint, poly) ) + // PG_RETURN_BOOL(FALSE); + //} #endif + + initGEOS(MAXIMUM_ALIGNOF); + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + + result = GEOSrelateWithin(g1,g2); + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + + if (result == 2) + { + elog(ERROR,"GEOS within() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + + + PG_RETURN_BOOL(result); +} + + + +PG_FUNCTION_INFO_V1(crosses); +Datum crosses(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + + Geometry *g1,*g2; + bool result; + + errorIfGeometryCollection(geom1,geom2); + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + + + + result = GEOSrelateCrosses(g1,g2); + + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + + if (result == 2) + { + elog(ERROR,"GEOS crosses() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + + + PG_RETURN_BOOL(result); +} + + + +PG_FUNCTION_INFO_V1(intersects); +Datum intersects(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + Geometry *g1,*g2; + bool result; + + + + + errorIfGeometryCollection(geom1,geom2); + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + + + + result = GEOSrelateIntersects(g1,g2); + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + if (result == 2) + { + elog(ERROR,"GEOS intersects() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + + + PG_RETURN_BOOL(result); +} + + +PG_FUNCTION_INFO_V1(touches); +Datum touches(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + + Geometry *g1,*g2; + bool result; + + errorIfGeometryCollection(geom1,geom2); + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + + + + result = GEOSrelateTouches(g1,g2); + + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + + if (result == 2) + { + elog(ERROR,"GEOS touches() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + + + PG_RETURN_BOOL(result); +} + + +PG_FUNCTION_INFO_V1(disjoint); +Datum disjoint(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + + Geometry *g1,*g2; + bool result; + + errorIfGeometryCollection(geom1,geom2); + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + + + result = GEOSrelateDisjoint(g1,g2); + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + + if (result == 2) + { + elog(ERROR,"GEOS disjoin() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + PG_RETURN_BOOL(result); +} + + +PG_FUNCTION_INFO_V1(relate_pattern); +Datum relate_pattern(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + char *patt; + bool result; + + Geometry *g1,*g2; + + errorIfGeometryCollection(geom1,geom2); + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1); + g2 = POSTGIS2GEOS(geom2); + + patt = DatumGetCString(DirectFunctionCall1(textout, + PointerGetDatum(PG_GETARG_DATUM(2)))); + + result = GEOSrelatePattern(g1,g2,patt); + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + pfree(patt); + + if (result == 2) + { + elog(ERROR,"GEOS relate_pattern() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + PG_RETURN_BOOL(result); + + +} + + + +PG_FUNCTION_INFO_V1(relate_full); +Datum relate_full(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + Geometry *g1,*g2; + char *relate_str; + int len; + char *result; + +//elog(NOTICE,"in relate_full()"); + + errorIfGeometryCollection(geom1,geom2); + + + initGEOS(MAXIMUM_ALIGNOF); + +//elog(NOTICE,"GEOS init()"); + + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + +//elog(NOTICE,"constructed geometries "); + + + + + +if ((g1==NULL) || (g2 == NULL)) + elog(NOTICE,"g1 or g2 are null"); + +//elog(NOTICE,GEOSasText(g1)); +//elog(NOTICE,GEOSasText(g2)); + +//elog(NOTICE,"valid g1 = %i", GEOSisvalid(g1)); +//elog(NOTICE,"valid g2 = %i",GEOSisvalid(g2)); + +//elog(NOTICE,"about to relate()"); + + + relate_str = GEOSrelate(g1, g2); + +//elog(NOTICE,"finished relate()"); + + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + + + + if (relate_str == NULL) + { + //free(relate_str); + elog(ERROR,"GEOS relate() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + + len = strlen(relate_str) + 4; + + result= palloc(len); + *((int *) result) = len; + + memcpy(result +4, relate_str, len-4); + + free(relate_str); + + + PG_RETURN_POINTER(result); +} + +//============================== + +PG_FUNCTION_INFO_V1(geomequals); +Datum geomequals(PG_FUNCTION_ARGS) +{ + LWGEOM *geom1 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *geom2 = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + + Geometry *g1,*g2; + bool result; + + errorIfGeometryCollection(geom1,geom2); + initGEOS(MAXIMUM_ALIGNOF); + + g1 = POSTGIS2GEOS(geom1 ); + g2 = POSTGIS2GEOS(geom2 ); + + + result = GEOSequals(g1,g2); + GEOSdeleteGeometry(g1); + GEOSdeleteGeometry(g2); + + if (result == 2) + { + elog(ERROR,"GEOS equals() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + PG_RETURN_BOOL(result); +} + +PG_FUNCTION_INFO_V1(issimple); +Datum issimple(PG_FUNCTION_ARGS) +{ + LWGEOM *geom = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + Geometry *g1; + int result; + + if (lwgeom_getnumgeometries(SERIALIZED_FORM(geom)) == 0) + PG_RETURN_BOOL(true); + + initGEOS(MAXIMUM_ALIGNOF); + + //elog(NOTICE,"GEOS init()"); + + g1 = POSTGIS2GEOS(geom ); + + result = GEOSisSimple(g1); + GEOSdeleteGeometry(g1); + + if (result == 2) + { + elog(ERROR,"GEOS issimple() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + PG_RETURN_BOOL(result); +} + +PG_FUNCTION_INFO_V1(isring); +Datum isring(PG_FUNCTION_ARGS) +{ + LWGEOM *geom = (LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + Geometry *g1; + int result; + + if (lwgeom_getType(geom->type) != LINETYPE) + { + elog(ERROR,"isring() should only be called on a LINE"); + } + + if (lwgeom_getnumgeometries(SERIALIZED_FORM(geom)) == 0) + PG_RETURN_BOOL(false); + + initGEOS(MAXIMUM_ALIGNOF); + + //elog(NOTICE,"GEOS init()"); + + g1 = POSTGIS2GEOS(geom ); + + result = GEOSisRing(g1); + GEOSdeleteGeometry(g1); + + if (result == 2) + { + elog(ERROR,"GEOS isring() threw an error!"); + PG_RETURN_NULL(); //never get here + } + + PG_RETURN_BOOL(result); +} + + + +//= GEOS <=> POSTGIS CONVERSION ========================= + +// Return a serialized line from a GEOS geometry (no SRID, no BBOX) +char *PointFromGeometry(Geometry *g, char want3d) +{ + POINTARRAY *pa; + LWPOINT *point; + char *srl; + POINT3D *pts; + int ptsize = want3d ? sizeof(POINT3D) : sizeof(POINT2D); + + // Construct point array + pa = (POINTARRAY *)palloc(sizeof(POINTARRAY)); + pa->ndims = want3d ? 3 : 2; + pa->npoints = 1; + + // Fill point array + pa->serialized_pointlist = palloc(ptsize); + pts = GEOSGetCoordinates(g); + memcpy(pa->serialized_pointlist, pts, ptsize); + GEOSdeleteChar( (char*) pts); + + // Construct LWPOINT + point = lwpoint_construct(pa->ndims, -1, pa); + + // Serialize LWPOINT + srl = lwpoint_serialize(point); + + return srl; +} + +// Return a serialized line from a GEOS geometry (no SRID no BBOX) +char *LineFromGeometry(Geometry *g, char want3d) +{ + POINTARRAY *pa; + LWLINE *line; + int npoints; + char *srl; + POINT3D *pts, *ip, *op; + int ptsize = want3d ? sizeof(POINT3D) : sizeof(POINT2D); + int i; + + npoints = GEOSGetNumCoordinate(g); + if (npoints <2) return NULL; + + // Construct point array + pa = (POINTARRAY *)palloc(sizeof(POINTARRAY)); + pa->ndims = want3d ? 3 : 2; + pa->npoints = npoints; + + // Fill point array + pa->serialized_pointlist = palloc(ptsize*npoints); + pts = GEOSGetCoordinates(g); + for (i=0; indims, -1, pa); + + // Serialize LWPOINT + srl = lwline_serialize(line); + + return srl; +} + +// Return a serialized polygon from a GEOS geometry (no SRID no BBOX) +char *PolyFromGeometry(Geometry *g, char want3d) +{ + POINTARRAY **rings, *pa; + LWPOLY *poly; + int ndims = want3d ? 3 : 2; + int nrings; + int npoints; + int i, j; + char *srl; + POINT3D *pts, *ip, *op; + int ptoff=0; // point offset inside POINT3D * + int ptsize = 8*ndims; + + // Get number of rings, and pointlist + pts = GEOSGetCoordinates(g); + nrings = GEOSGetNumInteriorRings(g); + rings = (POINTARRAY **)palloc(sizeof(POINTARRAY *)*(nrings+1)); + + // Exterior ring + npoints = GEOSGetNumCoordinate(GEOSGetExteriorRing(g)); + rings[0] = (POINTARRAY *)palloc(sizeof(POINTARRAY)); + pa = rings[0]; + pa->ndims = ndims; + pa->npoints = npoints; + + // Fill point array + pa->serialized_pointlist = palloc(ptsize*npoints); + for (i=0; indims = ndims; + pa->npoints = npoints; + + // Fill point array + pa->serialized_pointlist = palloc(ptsize*npoints); + for (i=0; indims, -1, nrings+1, rings); + + // Serialize LWPOLY + srl = lwpoly_serialize(poly); + + // Get rid of GEOS points... + GEOSdeleteChar( (char*) pts); + + return srl; +} + +//-----=GEOS2POSTGIS= + +/* + * Recursively add a Geometry to the LWGEOM_EXPLODED structure + * The exploded struct contains note about the number of dimensions + * requested. + */ +void +addToExploded_recursive(Geometry *geom, LWGEOM_EXPLODED *exp) +{ + char *srl; + int type = GEOSGeometryTypeId(geom) ; + char want3d = exp->ndims > 2 ? 1 : 0; + int ngeoms; + int t; + + + switch (type) + { + /* From slower to faster.. compensation rule :) */ + + case COLLECTIONTYPE: + case MULTIPOLYGONTYPE: + case MULTILINETYPE: + case MULTIPOINTTYPE: + + //this is more difficult because GEOS allows GCs of GCs + ngeoms = GEOSGetNumGeometries(geom); +#ifdef DEBUG_GEOS2POSTGIS + elog(NOTICE, "GEOS2POSTGIS: It's a MULTI (of %d elems)", ngeoms); +#endif + if (ngeoms == 0) + { + return; // skip EMPTY geoms + } + if (ngeoms == 1) { + Geometry *g = GEOSGetGeometryN(geom, 0); + // short cut! + return addToExploded_recursive(g, exp); + } + for (t=0; tnpolys++; + exp->polys = repalloc(exp->polys, + sizeof(char *)*exp->npolys); + exp->polys[exp->npolys-1] = srl; + return; + + case LINETYPE: +#ifdef DEBUG_GEOS2POSTGIS + elog(NOTICE, "GEOS2POSTGIS: It's a LINE"); +#endif + srl = LineFromGeometry(geom, want3d); + if (srl == NULL) return; + exp->nlines++; + exp->lines = repalloc(exp->lines, + sizeof(char *)*exp->nlines); + exp->lines[exp->nlines-1] = srl; + return; + + case POINTTYPE: +#ifdef DEBUG_GEOS2POSTGIS + elog(NOTICE, "GEOS2POSTGIS: It's a POINT"); +#endif + srl = PointFromGeometry(geom, want3d); + if (srl == NULL) return; + exp->npoints++; + exp->points = repalloc(exp->points, + sizeof(char *)*exp->npoints); + exp->points[exp->npoints-1] = srl; + return; + + default: +#ifdef DEBUG_GEOS2POSTGIS + elog(NOTICE, "GEOS2POSTGIS: It's UNKNOWN!"); +#endif + return; + + } + +} + +LWGEOM * +GEOS2POSTGIS(Geometry *geom, char want3d) +{ + char *serialized; + LWGEOM *result; + LWGEOM_EXPLODED *oexp; + int SRID = GEOSGetSRID(geom); + int wantbbox = 0; // might as well be 1 ... + + // Initialize exploded lwgeom + oexp = (LWGEOM_EXPLODED *)palloc(sizeof(LWGEOM_EXPLODED)); + oexp->SRID = SRID; + oexp->ndims = want3d ? 3 : 2; + oexp->npoints = 0; + oexp->points = palloc(sizeof(char *)); + oexp->nlines = 0; + oexp->lines = palloc(sizeof(char *)); + oexp->npolys = 0; + oexp->polys = palloc(sizeof(char *)); + + + addToExploded_recursive(geom, oexp); + serialized = lwexploded_serialize(oexp, wantbbox); + result = LWGEOM_construct(serialized, SRID, wantbbox); + return result; +} + +Geometry * +POSTGIS2GEOS(LWGEOM *geom) +{ + LWPOINT *point; + LWLINE *line; + LWPOLY *poly; + LWGEOM_EXPLODED *exploded; + int type = lwgeom_getType(geom->type); + int i, j; + Geometry **collected; + int ncollected; + + switch (type) + { + case POINTTYPE: +#ifdef DEBUG_POSTGIS2GEOS + elog(NOTICE, "POSTGIS2GEOS: it's a point"); +#endif + point = lwpoint_deserialize(SERIALIZED_FORM(geom)); + return PostGIS2GEOS_point(point); + + case LINETYPE: +#ifdef DEBUG_POSTGIS2GEOS + elog(NOTICE, "POSTGIS2GEOS: it's a line"); +#endif + line = lwline_deserialize(SERIALIZED_FORM(geom)); + return PostGIS2GEOS_linestring(line); + + case POLYGONTYPE: +#ifdef DEBUG_POSTGIS2GEOS + elog(NOTICE, "POSTGIS2GEOS: it's a polygon"); +#endif + poly = lwpoly_deserialize(SERIALIZED_FORM(geom)); + return PostGIS2GEOS_polygon(poly); + + default: + break; + } + +#ifdef DEBUG_POSTGIS2GEOS + elog(NOTICE, "POSTGIS2GEOS: out of switch"); +#endif + + // Since it is not a base type, let's explode... + exploded = lwgeom_explode(SERIALIZED_FORM(geom)); + +#ifdef DEBUG_POSTGIS2GEOS + elog(NOTICE, "POSTGIS2GEOS: exploded produced"); +#endif + + ncollected = exploded->npoints + exploded->nlines + exploded->npolys; + +#ifdef DEBUG_POSTGIS2GEOS + elog(NOTICE, "POSTGIS2GEOS: it's a collection of %d elems", ncollected); +#endif + + collected = (Geometry **)palloc(sizeof(Geometry *)*ncollected); + + j=0; + for (i=0; inpoints; i++) + { + point = lwpoint_deserialize(exploded->points[i]); + collected[j++] = PostGIS2GEOS_point(point); + } + for (i=0; inlines; i++) + { + line = lwline_deserialize(exploded->lines[j]); + collected[j++] = PostGIS2GEOS_linestring(line); + } + for (i=0; inpolys; i++) + { + poly = lwpoly_deserialize(exploded->polys[i]); + collected[j++] = PostGIS2GEOS_polygon(poly); + } + + return PostGIS2GEOS_collection(collected, ncollected, + exploded->SRID, exploded->ndims > 2); + +} + +PG_FUNCTION_INFO_V1(GEOSnoop); +Datum GEOSnoop(PG_FUNCTION_ARGS) +{ + LWGEOM *geom; + Geometry *geosgeom; + LWGEOM *result; + + initGEOS(MAXIMUM_ALIGNOF); + + geom = (LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); +#ifdef DEBUG_CONVERTER + elog(NOTICE, "GEOSnoop: IN: %s", unparse_WKT((char *)geom, malloc, free)); +#endif + + geosgeom = POSTGIS2GEOS(geom); + if ( (Pointer *)geom != (Pointer *)PG_GETARG_DATUM(0) ) pfree(geom); + if ( ! geosgeom ) PG_RETURN_NULL(); + + result = GEOS2POSTGIS(geosgeom, lwgeom_ndims(geom->type) > 2); + GEOSdeleteGeometry(geosgeom); + +#ifdef DEBUG_CONVERTER + elog(NOTICE, "GEOSnoop: OUT: %s", unparse_WKT((char *)result, malloc, free)); +#endif + + PG_RETURN_POINTER(result); +} + +#else // ! USE_GEOS + +PG_FUNCTION_INFO_V1(postgis_geos_version); +Datum postgis_geos_version(PG_FUNCTION_ARGS) +{ + //elog(ERROR,"GEOSversion:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); +} + +PG_FUNCTION_INFO_V1(intersection); +Datum intersection(PG_FUNCTION_ARGS) +{ + elog(ERROR,"intersection:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(convexhull); +Datum convexhull(PG_FUNCTION_ARGS) +{ + elog(ERROR,"convexhull:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(difference); +Datum difference(PG_FUNCTION_ARGS) +{ + elog(ERROR,"difference:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(boundary); +Datum boundary(PG_FUNCTION_ARGS) +{ + elog(ERROR,"boundary:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(symdifference); +Datum symdifference(PG_FUNCTION_ARGS) +{ + elog(ERROR,"symdifference:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(geomunion); +Datum geomunion(PG_FUNCTION_ARGS) +{ + elog(ERROR,"geomunion:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} + +PG_FUNCTION_INFO_V1(buffer); +Datum buffer(PG_FUNCTION_ARGS) +{ + elog(ERROR,"buffer:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} + + + +PG_FUNCTION_INFO_V1(relate_full); +Datum relate_full(PG_FUNCTION_ARGS) +{ + elog(ERROR,"relate_full:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(relate_pattern); +Datum relate_pattern(PG_FUNCTION_ARGS) +{ + elog(ERROR,"relate_pattern:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(disjoint); +Datum disjoint(PG_FUNCTION_ARGS) +{ + elog(ERROR,"disjoint:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(intersects); +Datum intersects(PG_FUNCTION_ARGS) +{ + elog(ERROR,"intersects:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(touches); +Datum touches(PG_FUNCTION_ARGS) +{ + elog(ERROR,"touches:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(crosses); +Datum crosses(PG_FUNCTION_ARGS) +{ + elog(ERROR,"crosses:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(within); +Datum within(PG_FUNCTION_ARGS) +{ + elog(ERROR,"within:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(contains); +Datum contains(PG_FUNCTION_ARGS) +{ + elog(ERROR,"contains:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(overlaps); +Datum overlaps(PG_FUNCTION_ARGS) +{ + elog(ERROR,"overlaps:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(isvalid); +Datum isvalid(PG_FUNCTION_ARGS) +{ + elog(ERROR,"isvalid:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} + +PG_FUNCTION_INFO_V1(issimple); +Datum issimple(PG_FUNCTION_ARGS) +{ + elog(ERROR,"issimple:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(geomequals); +Datum geomequals(PG_FUNCTION_ARGS) +{ + elog(ERROR,"geomequals:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} +PG_FUNCTION_INFO_V1(isring); +Datum isring(PG_FUNCTION_ARGS) +{ + elog(ERROR,"isring:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} + +PG_FUNCTION_INFO_V1(pointonsurface); +Datum pointonsurface(PG_FUNCTION_ARGS) +{ + elog(ERROR,"pointonsurface:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} + +PG_FUNCTION_INFO_V1(unite_garray); +Datum unite_garray(PG_FUNCTION_ARGS) +{ + elog(ERROR,"unite_garray:: operation not implemented - compile PostGIS with GEOS support"); + PG_RETURN_NULL(); // never get here +} + +PG_FUNCTION_INFO_V1(centroid); +Datum centroid(PG_FUNCTION_ARGS) +{ + elog(ERROR, "postgis version of centroid calculation not ported yet"); + PG_RETURN_NULL(); // never get here +} + +#endif // ! USE_GEOS diff --git a/lwgeom/lwpostgis.sql.in b/lwgeom/lwpostgis.sql.in index df5670db8..f0b27c48d 100644 --- a/lwgeom/lwpostgis.sql.in +++ b/lwgeom/lwpostgis.sql.in @@ -2503,6 +2503,127 @@ CREATEFUNCTION simplify(geometry, float8) AS '@MODULE_FILENAME@', 'LWGEOM_simplify2d' LANGUAGE 'C' WITH (isstrict); +--------------------------------------------------------------- +-- GEOS +--------------------------------------------------------------- + +CREATEFUNCTION intersection(geometry,geometry) + RETURNS geometry + AS '@MODULE_FILENAME@','intersection' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION buffer(geometry,float8) + RETURNS geometry + AS '@MODULE_FILENAME@','buffer' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION convexhull(geometry) + RETURNS geometry + AS '@MODULE_FILENAME@','convexhull' + LANGUAGE 'C' WITH (isstrict); + + +CREATEFUNCTION difference(geometry,geometry) + RETURNS geometry + AS '@MODULE_FILENAME@','difference' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION boundary(geometry) + RETURNS geometry + AS '@MODULE_FILENAME@','boundary' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION symdifference(geometry,geometry) + RETURNS geometry + AS '@MODULE_FILENAME@','symdifference' + LANGUAGE 'C' WITH (isstrict); + + +CREATEFUNCTION symmetricdifference(geometry,geometry) + RETURNS geometry + AS '@MODULE_FILENAME@','symdifference' + LANGUAGE 'C' WITH (isstrict); + + +CREATEFUNCTION GeomUnion(geometry,geometry) + RETURNS geometry + AS '@MODULE_FILENAME@','geomunion' + LANGUAGE 'C' WITH (isstrict); + +CREATE AGGREGATE MemGeomUnion ( + basetype = geometry, + sfunc = geomunion, + stype = geometry + ); + +CREATEFUNCTION unite_garray (geometry[]) + RETURNS geometry + AS '@MODULE_FILENAME@' + LANGUAGE 'C'; + +CREATE AGGREGATE GeomUnion ( + sfunc = geom_accum, + basetype = geometry, + stype = geometry[], + finalfunc = unite_garray + ); + + +CREATEFUNCTION relate(geometry,geometry) + RETURNS text + AS '@MODULE_FILENAME@','relate_full' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION relate(geometry,geometry,text) + RETURNS boolean + AS '@MODULE_FILENAME@','relate_pattern' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION disjoint(geometry,geometry) + RETURNS boolean + AS '@MODULE_FILENAME@' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION touches(geometry,geometry) + RETURNS boolean + AS '@MODULE_FILENAME@' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION intersects(geometry,geometry) + RETURNS boolean + AS '@MODULE_FILENAME@' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION crosses(geometry,geometry) + RETURNS boolean + AS '@MODULE_FILENAME@' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION within(geometry,geometry) + RETURNS boolean + AS '@MODULE_FILENAME@' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION contains(geometry,geometry) + RETURNS boolean + AS '@MODULE_FILENAME@' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION overlaps(geometry,geometry) + RETURNS boolean + AS '@MODULE_FILENAME@' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION isvalid(geometry) + RETURNS boolean + AS '@MODULE_FILENAME@' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION geosnoop(geometry) + RETURNS geometry + AS '@MODULE_FILENAME@', 'GEOSnoop' + LANGUAGE 'C' WITH (isstrict); + --------------------------------------------------------------- -- END ---------------------------------------------------------------