Changeset: 8a27bd9efb0f for MonetDB
URL: https://dev.monetdb.org/hg/MonetDB/rev/8a27bd9efb0f
Modified Files:
geom/monetdb5/geom.c
geom/monetdb5/geom.h
Branch: geo-update
Log Message:
Added simpler data types to represent geometries (instead of using GEOSGeom)
and refactored the code to use them.
diffs (truncated from 913 to 300 lines):
diff --git a/geom/monetdb5/geom.c b/geom/monetdb5/geom.c
--- a/geom/monetdb5/geom.c
+++ b/geom/monetdb5/geom.c
@@ -16,54 +16,22 @@
#include "mal_exception.h"
/**
- * Auxilary Functions
- *
- **/
-/* Destroy a Geos geometry */
-static void destroyGeom(GEOSGeom *g)
-{
- if (*g)
- GEOSGeom_destroy((*g));
- (*g) = NULL;
-}
-
-/* Converts wkb into Geos Geometries, if they are not NULL and have the same
SRID (used for geographic functions) */
-static str wkbGetComplatibleGeometries(wkb **a, wkb **b, GEOSGeom *ga,
GEOSGeom *gb)
-{
- str err = MAL_SUCCEED;
-
- if (is_wkb_nil(*a) || is_wkb_nil(*b))
- {
- destroyGeom(ga);
- destroyGeom(gb);
- return MAL_SUCCEED;
- }
- (*ga) = wkb2geos(*a);
- (*gb) = wkb2geos(*b);
- if ((*ga) == NULL || (*gb) == NULL)
- {
- destroyGeom(ga);
- destroyGeom(gb);
- err = createException(MAL, "geom.GetComplatibleGeometries",
SQLSTATE(38000) "Geos operation wkb2geos failed");
- }
- else if (GEOSGetSRID((*ga)) != GEOSGetSRID((*gb)))
- {
- destroyGeom(ga);
- destroyGeom(gb);
- err = createException(MAL, "geom.GetComplatibleGeometries",
SQLSTATE(38000) "Geometries of different SRID");
- }
- return err;
-}
-
-/**
* Convertions
*
**/
-/* Converts a latitude value in degrees to radians */
-static double longitudeDegreesToRadians(double lon_degrees)
-{
+const double earth_radius = 6371.009;
+const double earth_radius_meters = 6371009;
+
+/**
+ * Converts a longitude value in degrees to radians
+ * The normalization part was taken from PostGIS
+ */
+static double deg2RadLongitude(double lon_degrees)
+{
+ //Convert
double lon = M_PI * lon_degrees / 180.0;
+ //Normalize
if (lon == -1.0 * M_PI)
return M_PI;
if (lon == -2.0 * M_PI)
@@ -87,10 +55,15 @@ static double longitudeDegreesToRadians(
return lon;
}
-/* Converts a latitude value in degrees to radians */
-static double latitudeDegreesToRadians(double lat_degrees)
-{
+/**
+ * Converts a latitude value in degrees to radians
+ * The normalization part was taken from PostGIS
+ */
+static double deg2RadLatitude(double lat_degrees)
+{
+ //Convert
double lat = M_PI * lat_degrees / 180.0;
+ //Normalize
if (lat > 2.0 * M_PI)
lat = remainder(lat, 2.0 * M_PI);
@@ -112,135 +85,264 @@ static double latitudeDegreesToRadians(d
return lat;
}
+/* Converts the GeoPoint from degrees to radians latitude and longitude*/
+static GeoPoint deg2RadPoint(GeoPoint geo)
+{
+ geo.lon = deg2RadLongitude(geo.lon);
+ geo.lat = deg2RadLatitude(geo.lat);
+ return geo;
+}
+
+#if 0
/**
-* Convert spherical coordinates to cartesian coordinates on unit sphere
-* From PostGIS
-*/
-static void geog2cart(double lon, double lat, double *x, double *y, double *z)
-{
- (*x) = cos(lat) * cos(lon);
- (*y) = cos(lat) * sin(lon);
- (*z) = sin(lat);
-}
-
-static void geog2cart_r(double lon, double lat, double *x, double *y, double
*z)
-{
- double r = 6371.009;
- (*x) = r * cos(lat) * cos(lon);
- (*y) = r * cos(lat) * sin(lon);
- (*z) = r * sin(lat);
+ * Converts a longitude value in radians to degrees
+ * The normalization part was taken from PostGIS
+ */
+static double rad2DegLongitude(double lon_radians)
+{
+ //Convert
+ double lon = lon_radians * 180.0 / M_PI;
+ //Normalize
+ if (lon > 360.0)
+ lon = remainder(lon, 360.0);
+
+ if (lon < -360.0)
+ lon = remainder(lon, -360.0);
+
+ if (lon > 180.0)
+ lon = -360.0 + lon;
+
+ if (lon < -180.0)
+ lon = 360 + lon;
+
+ if (lon == -180.0)
+ return 180.0;
+
+ if (lon == -360.0)
+ return 0.0;
+
+ return lon;
}
/**
-* Convert cartesian coordinates on unit sphere to spherical coordinates
-* From PostGIS
-*/
-/*static void cart2geog(double x, double y, double z, double *lon, double *lat)
-{
- (*lon) = atan2(y, x);
- (*lat) = asin(z);
+ * Converts a latitude value in radians to degrees
+ * The normalization part was taken from PostGIS
+ */
+static double rad2DegLatitude(double lat_radians)
+{
+ //Convert
+ double lat = lat_radians * 180.0 / M_PI;
+ //Normalize
+ if (lat > 360.0)
+ lat = remainder(lat, 360.0);
+
+ if (lat < -360.0)
+ lat = remainder(lat, -360.0);
+
+ if (lat > 180.0)
+ lat = 180.0 - lat;
+
+ if (lat < -180.0)
+ lat = -180.0 - lat;
+
+ if (lat > 90.0)
+ lat = 180.0 - lat;
+
+ if (lat < -90.0)
+ lat = -180.0 - lat;
+
+ return lat;
+}
+
+/* Converts the GeoPoint from degrees to radians latitude and longitude*/
+static void rad2DegPoint(GeoPoint *geo)
+{
+ geo->lon = rad2DegLongitude(geo->lon);
+ geo->lat = rad2DegLatitude(geo->lat);
+}
+#endif
+
+static GeoPoint geoPointFromGeom(GEOSGeom geom)
+{
+ GeoPoint geo;
+ GEOSGeomGetX(geom, &(geo.lon));
+ GEOSGeomGetY(geom, &(geo.lat));
+ return geo;
+}
+
+/*static GeoLine geoLineFromGeom(GEOSGeom geom)
+{
+ GeoLine geo;
+ GEOSGeomGetX(GEOSGeomGetStartPoint(geom), &(geo.start.lon));
+ GEOSGeomGetY(GEOSGeomGetStartPoint(geom), &(geo.start.lat));
+ GEOSGeomGetX(GEOSGeomGetEndPoint(geom), &(geo.end.lon));
+ GEOSGeomGetY(GEOSGeomGetEndPoint(geom), &(geo.end.lat));
+ return geo;
}*/
-static void cart2geog_r(double x, double y, double z, double *lon, double *lat)
-{
- double r = 6371.009;
- (*lon) = atan2(y, x);
- (*lat) = asin(z / r);
-}
-
-/* Returns the latitude and longitude (in radians) from a geographic Geom
point */
-static void pointToRadian(GEOSGeom geom, double *lat_r, double *lon_r)
-{
- double lat_d, lon_d;
- GEOSGeomGetX(geom, &lon_d);
- GEOSGeomGetY(geom, &lat_d);
-
- (*lat_r) = latitudeDegreesToRadians(lat_d);
- (*lon_r) = longitudeDegreesToRadians(lon_d);
-}
-
-/* Converts two lat/lon points into cartesian coordinates and creates a Line
geometry */
-static GEOSGeom geographicPointsToCartesianLine(double lat1, double lon1,
double lat2, double lon2)
-{
- double x1, y1, z1, x2, y2, z2;
- geog2cart(lon1, lat1, &x1, &y1, &z1);
- geog2cart(lon2, lat2, &x2, &y2, &z2);
- GEOSCoordSequence *lineSeq = GEOSCoordSeq_create(2, 3);
- GEOSCoordSeq_setXYZ(lineSeq, 0, x1, y1, z1);
- GEOSCoordSeq_setXYZ(lineSeq, 1, x2, y2, z2);
- return GEOSGeom_createLineString(lineSeq);
-}
+static GeoLines geoLinesFromGeom(GEOSGeom geom)
+{
+ const GEOSCoordSequence *gcs = GEOSGeom_getCoordSeq(geom);
+ int segmentCount = GEOSGeomGetNumPoints(geom) - 1;
+ GeoLines geo;
+ geo.segmentCount = segmentCount;
+ //TODO Malloc fail exception?
+ geo.segments = GDKmalloc(sizeof(GeoLine) * segmentCount);
+ for (int i = 0; i < segmentCount; i++)
+ {
+ GEOSCoordSeq_getXY(gcs, i, &geo.segments[i].start.lon,
&geo.segments[i].start.lat);
+ GEOSCoordSeq_getXY(gcs, i + 1, &geo.segments[i].end.lon,
&geo.segments[i].end.lat);
+ }
+ return geo;
+}
+
+static GeoPoint geoPointFromLatLon(double lon, double lat)
+{
+ GeoPoint geo;
+ geo.lon = lon;
+ geo.lat = lat;
+ return geo;
+}
+
+static CartPoint cartPointFromXYZ(double x, double y, double z)
+{
+ CartPoint cart;
+ cart.x = x;
+ cart.y = y;
+ cart.z = z;
+ return cart;
+}
+
+/* Converts Well-Known Bytes into Geos Geometries, if they are not NULL and
have the same SRID (used for geographic functions) */
+static str wkbGetComplatibleGeometries(wkb **a, wkb **b, GEOSGeom *ga,
GEOSGeom *gb)
+{
+ str err = MAL_SUCCEED;
+
+ if (is_wkb_nil(*a) || is_wkb_nil(*b))
+ {
+ (*ga) = NULL;
+ (*gb) = NULL;
+ return MAL_SUCCEED;
+ }
+ (*ga) = wkb2geos(*a);
+ (*gb) = wkb2geos(*b);
+ if ((*ga) == NULL || (*gb) == NULL)
+ {
+ err = createException(MAL, "geom.wkbGetComplatibleGeometries",
SQLSTATE(38000) "Geos operation wkb2geos failed");
+ }
+ else if (GEOSGetSRID((*ga)) != GEOSGetSRID((*gb)))
+ {
+ GEOSGeom_destroy(*ga);
+ GEOSGeom_destroy(*gb);
+ err = createException(MAL, "geom.wkbGetComplatibleGeometries",
SQLSTATE(38000) "Geometries of different SRID");
+ }
+ return err;
+}
+
+/**
+* Convert spherical coordinates to cartesian coordinates on unit sphere
_______________________________________________
checkin-list mailing list
[email protected]
https://www.monetdb.org/mailman/listinfo/checkin-list
