Changeset: d7a76882c516 for MonetDB
URL: http://dev.monetdb.org/hg/MonetDB?cmd=changeset;node=d7a76882c516
Added Files:
geom/sql/Tests/transform.sql
geom/sql/Tests/transform.stable.err
geom/sql/Tests/transform.stable.out
Modified Files:
geom/lib/libgeom.h
geom/monetdb5/geom.c
geom/sql/Tests/All
Branch: geo
Log Message:
transform works also for polygon
diffs (truncated from 714 to 300 lines):
diff --git a/geom/lib/libgeom.h b/geom/lib/libgeom.h
--- a/geom/lib/libgeom.h
+++ b/geom/lib/libgeom.h
@@ -101,6 +101,7 @@ typedef enum wkb_type {
//wkbGeometry = 0,
wkbPoint = 1,
wkbLineString = 2,
+ wkbLinearRing = 3,
wkbPolygon = 4,
wkbMultiPoint = 5,
wkbMultiLineString = 6,
diff --git a/geom/monetdb5/geom.c b/geom/monetdb5/geom.c
--- a/geom/monetdb5/geom.c
+++ b/geom/monetdb5/geom.c
@@ -53,6 +53,7 @@
#define GEOMETRY_HAS_Z(info)(info & 0x02)
#define GEOMETRY_HAS_M(info)(info & 0x01)
+#define PI 3.14159265358979323846
/* the first argument in the functions is the return variable */
@@ -152,95 +153,286 @@ geom_export str wkbNumGeometries(int* ou
geom_export str wkbTransform(wkb**, wkb**, int*, char**, char**);
-/* It gets a geometry and transforms its coordinates to the provided srid */
-str wkbTransform(wkb** transformedWKB, wkb** geomWKB, int* srid, char**
proj4_src_str, char** proj4_dst_str) {
- projPJ proj4_src = pj_init_plus(*proj4_src_str);
- projPJ proj4_dst = pj_init_plus(*proj4_dst_str);
+/** convert degrees to radians */
+static void degrees2radians(double *x, double *y, double *z) {
+ (*x) *= PI/180.0;
+ (*y) *= PI/180.0;
+ (*z) *= PI/180.0;
+}
+
+/** convert radians to degrees */
+static void radians2degrees(double *x, double *y, double *z) {
+ (*x) *= 180.0/PI;
+ (*y) *= 180.0/PI;
+ (*z) *= 180.0/PI;
+}
+
+
+static str transformCoordSeq(int idx, int coordinatesNum, projPJ proj4_src,
projPJ proj4_dst, const GEOSCoordSequence* gcs_old, GEOSCoordSequence**
gcs_new){
+ double x=0, y=0, z=0;
+
+ GEOSCoordSeq_getX(gcs_old, idx, &x);
+ GEOSCoordSeq_getY(gcs_old, idx, &y);
+
+ //fprintf(stderr, "transforming POINT(%f %f) from '%s' to '%s'\n", x,
y, pj_get_def(proj4_src,0), pj_get_def(proj4_dst,0));
+
+ if(coordinatesNum > 2)
+ GEOSCoordSeq_getZ(gcs_old, idx, &z);
+
+ /* check if the passed reference system is geographic (proj=latlong)
+ * and change the degrees to radians because pj_transform works with
radians*/
+ if (pj_is_latlong(proj4_src)) degrees2radians(&x, &y, &z) ;
+
+
+ pj_transform(proj4_src, proj4_dst, 1, 0, &x, &y, &z);
+
+ /* check if the destination reference system is geographic and change
+ * the destination coordinates from radians to degrees */
+ if (pj_is_latlong(proj4_dst)) radians2degrees(&x, &y, &z);
+
+
+ //fprintf(stderr, "transformed to POINT(%f %f)\n", x, y);
+
+ GEOSCoordSeq_setX(*gcs_new, idx, x);
+ GEOSCoordSeq_setY(*gcs_new, idx, y);
- GEOSGeom geosGeometry, transformedGeosGeometry;
-// int geometriesNum = 0;
+ if(coordinatesNum > 2)
+ GEOSCoordSeq_setZ(*gcs_new, idx, z);
+
+
+
+ return MAL_SUCCEED;
+}
+
+static str transformPoint(GEOSGeometry** transformedGeometry, const
GEOSGeometry* geosGeometry, projPJ proj4_src, projPJ proj4_dst) {
int coordinatesNum = 0;
const GEOSCoordSequence* gcs_old;
GEOSCoordSeq gcs_new;
+
+ /* get the coordinates of the points comprising the geometry */
+ gcs_old = GEOSGeom_getCoordSeq(geosGeometry);
+
+ if(gcs_old == NULL) {
+ *transformedGeometry = NULL;
+ return createException(MAL, "geom.wkbTransform",
"GEOSGeom_getCoordSeq failed");
+ }
+
+ /* create the coordinates sequence for the transformed geometry */
+ gcs_new = GEOSCoordSeq_create(1, coordinatesNum);
+
+ /* create the transformed coordinates */
+ transformCoordSeq(0, coordinatesNum, proj4_src, proj4_dst, gcs_old,
&gcs_new);
+
+ /* create the geometry from the coordinates seqience */
+ *transformedGeometry = GEOSGeom_createPoint(gcs_new);
+
+ return MAL_SUCCEED;
+}
+
+static str transformLine(GEOSCoordSeq *gcs_new, const GEOSGeometry*
geosGeometry, projPJ proj4_src, projPJ proj4_dst) {
+ int coordinatesNum = 0;
+ const GEOSCoordSequence* gcs_old;
unsigned int pointsNum =0, i=0;
- double x=0, y=0, z=0;
-
-
- if(*geomWKB == NULL) {
- *transformedWKB = wkb_nil;
- throw(MAL, "geom.Transform", "wkb is null");
- }
-
- /* get the geosGeometry from the wkb */
- geosGeometry = wkb2geos(*geomWKB);
-// /*get the number of geometries */
-// geometriesNum = GEOSGetNumGeometries(geosGeometry);
+
/* get the number of coordinates the geometry has */
coordinatesNum = GEOSGeom_getCoordinateDimension(geosGeometry);
/* get the coordinates of the points comprising the geometry */
gcs_old = GEOSGeom_getCoordSeq(geosGeometry);
- if(gcs_old == NULL) {
+
+ if(gcs_old == NULL)
+ return createException(MAL, "geom.wkbTransform",
"GEOSGeom_getCoordSeq failed");
+
+ /* get the number of points in the geometry */
+ GEOSCoordSeq_getSize(gcs_old, &pointsNum);
+
+ /* create the coordinates sequence for the transformed geometry */
+ *gcs_new = GEOSCoordSeq_create(pointsNum, coordinatesNum);
+
+ /* create the transformed coordinates */
+ for(i=0; i<pointsNum; i++)
+ transformCoordSeq(i, coordinatesNum, proj4_src, proj4_dst,
gcs_old, gcs_new);
+
+ return MAL_SUCCEED;
+}
+
+static str transformLineString(GEOSGeometry** transformedGeometry, const
GEOSGeometry* geosGeometry, projPJ proj4_src, projPJ proj4_dst) {
+ GEOSCoordSeq coordSeq;
+ str ret = MAL_SUCCEED;
+
+ ret = transformLine(&coordSeq, geosGeometry, proj4_src, proj4_dst);
+
+ if(ret != MAL_SUCCEED) {
+ *transformedGeometry = NULL;
+ return ret;
+ }
+
+ /* create the geometry from the coordinates sequence */
+ *transformedGeometry = GEOSGeom_createLineString(coordSeq);
+
+ return ret;
+}
+
+static str transformLinearRing(GEOSGeometry** transformedGeometry, const
GEOSGeometry* geosGeometry, projPJ proj4_src, projPJ proj4_dst) {
+ GEOSCoordSeq coordSeq;
+ str ret = MAL_SUCCEED;
+
+ ret = transformLine(&coordSeq, geosGeometry, proj4_src, proj4_dst);
+
+ if(ret != MAL_SUCCEED) {
+ *transformedGeometry = NULL;
+ return ret;
+ }
+
+ /* create the geometry from the coordinates sequence */
+ *transformedGeometry = GEOSGeom_createLinearRing(coordSeq);
+
+ return ret;
+}
+
+static str transformPolygon(GEOSGeometry** transformedGeometry, const
GEOSGeometry* geosGeometry, projPJ proj4_src, projPJ proj4_dst, int srid) {
+ const GEOSGeometry* exteriorRingGeometry;
+ GEOSGeometry* transformedExteriorRingGeometry = NULL;
+ GEOSGeometry** transformedInteriorRingGeometries = NULL;
+ int numInteriorRings=0, i=0;
+
+ /* get the exterior ring of the polygon */
+ exteriorRingGeometry = GEOSGetExteriorRing(geosGeometry);
+ if(!exteriorRingGeometry) {
+ *transformedGeometry = NULL;
+ return createException(MAL,
"geom.wkbTransform","GEOSGetExteriorRing failed");
+ }
+
+ transformLinearRing(&transformedExteriorRingGeometry,
exteriorRingGeometry, proj4_src, proj4_dst);
+ GEOSSetSRID(transformedExteriorRingGeometry, srid);
+
+ numInteriorRings = GEOSGetNumInteriorRings(geosGeometry);
+ if (numInteriorRings == -1 ) {
+ *transformedGeometry = NULL;
+ GEOSGeom_destroy(transformedExteriorRingGeometry);
+ return createException(MAL, "geom.wkbTransform",
"GEOSGetInteriorRingN failed.");
+ }
+
+ /* iterate over the interiorRing and transform each one of them */
+ transformedInteriorRingGeometries =
GDKmalloc(numInteriorRings*sizeof(GEOSGeometry*));
+ for(i=0; i<numInteriorRings; i++) {
+ transformLinearRing(&(transformedInteriorRingGeometries[i]),
GEOSGetInteriorRingN(geosGeometry, i), proj4_src, proj4_dst);
+ GEOSSetSRID(transformedInteriorRingGeometries[i], srid);
+ }
+
+ *transformedGeometry =
GEOSGeom_createPolygon(transformedExteriorRingGeometry,
transformedInteriorRingGeometries, numInteriorRings);
+ return MAL_SUCCEED;
+}
+
+/* the following function is used in postgis to get projPJ from str.
+ * it is necessary to do it ina detailed way like that because pj_init_plus
+ * does not set all parameters correctly and I cannot test whether the
+ * coordinate reference systems are geographic or not */
+ static projPJ projFromStr(char* projStr) {
+ int t;
+ char *params[1024]; // one for each parameter
+ char *loc;
+ char *str;
+ size_t slen;
+ projPJ result;
+
+
+ if (projStr == NULL) return NULL;
+
+ slen = strlen(projStr);
+
+ if (slen == 0) return NULL;
+
+ str = GDKmalloc(slen+1);
+ strcpy(str, projStr);
+
+ // first we split the string into a bunch of smaller strings,
+ // based on the " " separator
+
+ params[0] = str; // 1st param, we'll null terminate at the " " soon
+
+ loc = str;
+ t = 1;
+ while ((loc != NULL) && (*loc != 0) )
+ {
+ loc = strchr(loc, ' ');
+ if (loc != NULL)
+ {
+ *loc = 0; // null terminate
+ params[t] = loc+1;
+ loc++; // next char
+ t++; //next param
+ }
+ }
+
+ if (!(result=pj_init(t, params)))
+ {
+ GDKfree(str);
+ return NULL;
+ }
+ GDKfree(str);
+ return result;
+
+}
+
+/* It gets a geometry and transforms its coordinates to the provided srid */
+str wkbTransform(wkb** transformedWKB, wkb** geomWKB, int* srid, char**
proj4_src_str, char** proj4_dst_str) {
+ projPJ proj4_src = /*pj_init_plus*/projFromStr(*proj4_src_str);
+ projPJ proj4_dst = /*pj_init_plus*/projFromStr(*proj4_dst_str);
+
+ GEOSGeom geosGeometry, transformedGeosGeometry;
+ int geometryType = -1;
+
+ str ret = MAL_SUCCEED;
+
+//fprintf(stderr, "source:\t%s\n", *proj4_src_str);
+//fprintf(stderr, "destination:\t%s\n", *proj4_dst_str);
+
+ if(*geomWKB == NULL) {
*transformedWKB = wkb_nil;
- throw(MAL, "geom.wkbTransform", "GEOSGeom_getCoordSeq failed");
+ pj_free(proj4_src);
+ pj_free(proj4_dst);
+ throw(MAL, "geom.Transform", "wkb is null");
}
- /* get the number of points in the geomtry */
- GEOSCoordSeq_getSize(gcs_old, &pointsNum);
-
- /* create the coordinates sequence for the transformed geometry */
- gcs_new = GEOSCoordSeq_create(pointsNum, coordinatesNum);
+
+ /* get the geosGeometry from the wkb */
+ geosGeometry = wkb2geos(*geomWKB);
+ /* get the type of the geometry */
+ geometryType = GEOSGeomTypeId(geosGeometry)+1;
+
+ switch(geometryType) {
+ case wkbPoint:
+ ret = transformPoint(&transformedGeosGeometry,
geosGeometry, proj4_src, proj4_dst);
_______________________________________________
checkin-list mailing list
checkin-list@monetdb.org
https://www.monetdb.org/mailman/listinfo/checkin-list
