Changeset: 26c8480c79a9 for MonetDB
URL: https://dev.monetdb.org/hg/MonetDB/rev/26c8480c79a9
Added Files:
geom/sql/functions/Tests/ST_Transform.reqtests
geom/sql/functions/Tests/ST_Transform.test
Modified Files:
geom/lib/libgeom.h
geom/monetdb5/geom.c
geom/sql/functions/Tests/All
Branch: geo-update
Log Message:
Updated PROJ library calls to latest version (we now only support versions from
PROJ 6) in st_transform functions (+ the library include); added tests for
ST_Transform
diffs (truncated from 436 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
@@ -27,10 +27,8 @@
#endif
#include <geos_c.h>
-
#ifdef HAVE_PROJ
-#define ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
-#include <proj_api.h> //it is needed to transform from one srid to another
+#include <proj.h>
#endif
/* geos does not support 3d envelope */
diff --git a/geom/monetdb5/geom.c b/geom/monetdb5/geom.c
--- a/geom/monetdb5/geom.c
+++ b/geom/monetdb5/geom.c
@@ -178,7 +178,6 @@ wkbCollectAggrSubGroupedCand(bat *outid,
geomCollectionType =
GEOSGeom_getCollectionType(GEOSGeomTypeId(inGEOM));
//srid = GEOSGetSRID(inGEOM);
}
-
unionGroup[geomCount] = inGEOM;
geomCount += 1;
if (geomCollectionType != GEOS_GEOMETRYCOLLECTION &&
GEOSGeom_getCollectionType(GEOSGeomTypeId(inGEOM)) != geomCollectionType)
@@ -311,59 +310,57 @@ geometryHasM(int info)
#define M_PI ((double) 3.14159265358979323846) /* pi */
#endif
+//TODO Remove?
/** convert degrees to radians */
-static inline void
+/*static inline void
degrees2radians(double *x, double *y, double *z)
{
double val = M_PI / 180.0;
*x *= val;
*y *= val;
*z *= val;
-}
-
+}*/
+
+//TODO Remove?
/** convert radians to degrees */
-static inline void
+/*static inline void
radians2degrees(double *x, double *y, double *z)
{
double val = 180.0 / M_PI;
*x *= val;
*y *= val;
*z *= val;
-}
+}*/
static str
-transformCoordSeq(int idx, int coordinatesNum, projPJ proj4_src, projPJ
proj4_dst, const GEOSCoordSequence *gcs_old, GEOSCoordSeq gcs_new)
+transformCoordSeq(int idx, int coordinatesNum, PJ *P, const GEOSCoordSequence
*gcs_old, GEOSCoordSeq gcs_new)
{
double x = 0, y = 0, z = 0;
- int *errorNum = 0;
+ int errorNum = 0;
+ PJ_COORD c, c_out;
if (!GEOSCoordSeq_getX(gcs_old, idx, &x) ||
!GEOSCoordSeq_getY(gcs_old, idx, &y) ||
(coordinatesNum > 2 && !GEOSCoordSeq_getZ(gcs_old, idx, &z)))
throw(MAL, "geom.Transform", SQLSTATE(38000) "Geos cannot get
coordinates");
- /* 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);
-
- errorNum = pj_get_errno_ref();
- if (*errorNum != 0) {
+ c.lpzt.lam=x;
+ c.lpzt.phi=y;
+ c.lpzt.z=z;
+ c.lpzt.t = HUGE_VAL;
+
+ c_out = proj_trans(P, PJ_FWD, c);
+
+ errorNum = proj_errno(P);
+ if (errorNum != 0) {
if (coordinatesNum > 2)
- throw(MAL, "geom.Transform", SQLSTATE(38000) "Geos
cannot transform point (%f %f %f): %s\n", x, y, z, pj_strerrno(*errorNum));
+ throw(MAL, "geom.Transform", SQLSTATE(38000) "Geos
cannot transform point (%f %f %f): %s\n", x, y, z, proj_errno_string(errorNum));
else
- throw(MAL, "geom.Transform", SQLSTATE(38000) "Geos
cannot transform point (%f %f): %s\n", x, y, pj_strerrno(*errorNum));
- }
-
- /* 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);
-
- if (!GEOSCoordSeq_setX(gcs_new, idx, x) ||
- !GEOSCoordSeq_setY(gcs_new, idx, y) ||
+ throw(MAL, "geom.Transform", SQLSTATE(38000) "Geos
cannot transform point (%f %f): %s\n", x, y, proj_errno_string(errorNum));
+ }
+
+ if (!GEOSCoordSeq_setX(gcs_new, idx,c_out.xy.x) ||
+ !GEOSCoordSeq_setY(gcs_new, idx,c_out.xy.y) ||
(coordinatesNum > 2 && !GEOSCoordSeq_setZ(gcs_new, idx, z)))
throw(MAL, "geom.Transform", SQLSTATE(38000) "Geos cannot set
coordinates");
@@ -371,7 +368,7 @@ transformCoordSeq(int idx, int coordinat
}
static str
-transformPoint(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, projPJ proj4_src, projPJ proj4_dst)
+transformPoint(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P)
{
int coordinatesNum = 0;
const GEOSCoordSequence *gcs_old;
@@ -394,7 +391,7 @@ transformPoint(GEOSGeometry **transforme
throw(MAL, "geom.Transform", SQLSTATE(38000) "Geos operation
GEOSGeom_getCoordSeq failed");
/* create the transformed coordinates */
- ret = transformCoordSeq(0, coordinatesNum, proj4_src, proj4_dst,
gcs_old, gcs_new);
+ ret = transformCoordSeq(0, coordinatesNum, P, gcs_old, gcs_new);
if (ret != MAL_SUCCEED) {
GEOSCoordSeq_destroy(gcs_new);
return ret;
@@ -411,7 +408,7 @@ transformPoint(GEOSGeometry **transforme
}
static str
-transformLine(GEOSCoordSeq *gcs_new, const GEOSGeometry *geosGeometry, projPJ
proj4_src, projPJ proj4_dst)
+transformLine(GEOSCoordSeq *gcs_new, const GEOSGeometry *geosGeometry, PJ *P)
{
int coordinatesNum = 0;
const GEOSCoordSequence *gcs_old;
@@ -437,7 +434,7 @@ transformLine(GEOSCoordSeq *gcs_new, con
/* create the transformed coordinates */
for (i = 0; i < pointsNum; i++) {
- ret = transformCoordSeq(i, coordinatesNum, proj4_src,
proj4_dst, gcs_old, *gcs_new);
+ ret = transformCoordSeq(i, coordinatesNum, P, gcs_old,
*gcs_new);
if (ret != MAL_SUCCEED) {
GEOSCoordSeq_destroy(*gcs_new);
*gcs_new = NULL;
@@ -449,12 +446,12 @@ transformLine(GEOSCoordSeq *gcs_new, con
}
static str
-transformLineString(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, projPJ proj4_src, projPJ proj4_dst)
+transformLineString(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P)
{
GEOSCoordSeq coordSeq;
str ret = MAL_SUCCEED;
- ret = transformLine(&coordSeq, geosGeometry, proj4_src, proj4_dst);
+ ret = transformLine(&coordSeq, geosGeometry, P);
if (ret != MAL_SUCCEED) {
*transformedGeometry = NULL;
return ret;
@@ -471,12 +468,12 @@ transformLineString(GEOSGeometry **trans
}
static str
-transformLinearRing(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, projPJ proj4_src, projPJ proj4_dst)
+transformLinearRing(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P)
{
GEOSCoordSeq coordSeq = NULL;
str ret = MAL_SUCCEED;
- ret = transformLine(&coordSeq, geosGeometry, proj4_src, proj4_dst);
+ ret = transformLine(&coordSeq, geosGeometry, P);
if (ret != MAL_SUCCEED) {
*transformedGeometry = NULL;
return ret;
@@ -493,7 +490,7 @@ transformLinearRing(GEOSGeometry **trans
}
static str
-transformPolygon(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, projPJ proj4_src, projPJ proj4_dst, int srid)
+transformPolygon(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P, int srid)
{
const GEOSGeometry *exteriorRingGeometry;
GEOSGeometry *transformedExteriorRingGeometry = NULL;
@@ -508,7 +505,7 @@ transformPolygon(GEOSGeometry **transfor
throw(MAL, "geom.Transform", SQLSTATE(38000) "Geos operation
GEOSGetExteriorRing failed");
}
- ret = transformLinearRing(&transformedExteriorRingGeometry,
exteriorRingGeometry, proj4_src, proj4_dst);
+ ret = transformLinearRing(&transformedExteriorRingGeometry,
exteriorRingGeometry, P);
if (ret != MAL_SUCCEED) {
*transformedGeometry = NULL;
return ret;
@@ -532,7 +529,7 @@ transformPolygon(GEOSGeometry **transfor
throw(MAL, "geom.Transform", SQLSTATE(HY013)
MAL_MALLOC_FAIL);
}
for (i = 0; i < numInteriorRings; i++) {
- ret =
transformLinearRing(&transformedInteriorRingGeometries[i],
GEOSGetInteriorRingN(geosGeometry, i), proj4_src, proj4_dst);
+ ret =
transformLinearRing(&transformedInteriorRingGeometries[i],
GEOSGetInteriorRingN(geosGeometry, i), P);
if (ret != MAL_SUCCEED) {
while (--i >= 0)
GEOSGeom_destroy(transformedInteriorRingGeometries[i]);
@@ -546,19 +543,19 @@ transformPolygon(GEOSGeometry **transfor
}
*transformedGeometry =
GEOSGeom_createPolygon(transformedExteriorRingGeometry,
transformedInteriorRingGeometries, numInteriorRings);
+
if (*transformedGeometry == NULL) {
for (i = 0; i < numInteriorRings; i++)
GEOSGeom_destroy(transformedInteriorRingGeometries[i]);
ret = createException(MAL, "geom.Transform", SQLSTATE(38000)
"Geos operation GEOSGeom_createPolygon failed");
}
GDKfree(transformedInteriorRingGeometries);
- GEOSGeom_destroy(transformedExteriorRingGeometry);
-
+ //GEOSGeom_destroy(transformedExteriorRingGeometry);
return ret;
}
static str
-transformMultiGeometry(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, projPJ proj4_src, projPJ proj4_dst, int srid, int geometryType)
+transformMultiGeometry(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P, int srid, int geometryType)
{
int geometriesNum, subGeometryType, i;
GEOSGeometry **transformedMultiGeometries = NULL;
@@ -580,16 +577,16 @@ transformMultiGeometry(GEOSGeometry **tr
else {
switch (subGeometryType) {
case wkbPoint_mdb:
- ret =
transformPoint(&transformedMultiGeometries[i], multiGeometry, proj4_src,
proj4_dst);
+ ret =
transformPoint(&transformedMultiGeometries[i], multiGeometry, P);
break;
case wkbLineString_mdb:
- ret =
transformLineString(&transformedMultiGeometries[i], multiGeometry, proj4_src,
proj4_dst);
+ ret =
transformLineString(&transformedMultiGeometries[i], multiGeometry, P);
break;
case wkbLinearRing_mdb:
- ret =
transformLinearRing(&transformedMultiGeometries[i], multiGeometry, proj4_src,
proj4_dst);
+ ret =
transformLinearRing(&transformedMultiGeometries[i], multiGeometry, P);
break;
case wkbPolygon_mdb:
- ret =
transformPolygon(&transformedMultiGeometries[i], multiGeometry, proj4_src,
proj4_dst, srid);
+ ret =
transformPolygon(&transformedMultiGeometries[i], multiGeometry, P, srid);
break;
default:
ret = createException(MAL, "geom.Transform",
SQLSTATE(38000) "Geos unknown geometry type");
@@ -618,48 +615,6 @@ transformMultiGeometry(GEOSGeometry **tr
return ret;
}
-
-/* the following function is used in postgis to get projPJ from str.
- * it is necessary to do it in a 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(const char *projStr)
-{
- int t;
- char *params[1024]; // one for each parameter
- char *loc;
- char *str;
- projPJ result;
-
- if (projStr == NULL)
- return NULL;
-
- str = GDKstrdup(projStr);
- if (str == NULL)
- return NULL;
-
- // 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
-
- t = 1;
- for (loc = strchr(str, ' '); loc != NULL; loc = strchr(loc, ' ')) {
- if (t == (int) (sizeof(params) / sizeof(params[0]))) {
- /* too many parameters */
- GDKfree(str);
- return NULL;
- }
- *loc++ = 0; // null terminate and advance
- params[t++] = loc;
- }
-
- result = pj_init(t, params);
- GDKfree(str);
-
- return result;
-}
#endif
/* It gets a geometry and transforms its coordinates to the provided srid */
@@ -673,9 +628,9 @@ wkbTransform(wkb **transformedWKB, wkb *
(void) srid_dst;
(void) proj4_src_str;
_______________________________________________
checkin-list mailing list -- checkin-list@monetdb.org
To unsubscribe send an email to checkin-list-le...@monetdb.org
