Changeset: db68ff963637 for MonetDB
URL: https://dev.monetdb.org/hg/MonetDB/rev/db68ff963637
Modified Files:
geom/monetdb5/geom.c
geom/monetdb5/geom.h
geom/monetdb5/geom.mal
geom/sql/40_geom.sql
Branch: geo-update
Log Message:
First version of distance calculation for geographic points. Geographic
distance implemented for all (non-multi) geometric data types. Missing
perpendicular calculation on Line/Ring distance and Polygon/Polygon
intersection check.
diffs (truncated from 7219 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
@@ -15,6 +15,415 @@
#include "gdk_logger.h"
#include "mal_exception.h"
+/**
+ * Convertions
+ *
+ **/
+
+/* Converts a latitude value in degrees to radians */
+static double longitudeDegreesToRadians(double lon_degrees)
+{
+ double lon = M_PI * lon_degrees / 180.0;
+ if (lon == -1.0 * M_PI)
+ return M_PI;
+ if (lon == -2.0 * M_PI)
+ return 0.0;
+
+ if (lon > 2.0 * M_PI)
+ lon = remainder(lon, 2.0 * M_PI);
+
+ if (lon < -2.0 * M_PI)
+ lon = remainder(lon, -2.0 * M_PI);
+
+ if (lon > M_PI)
+ lon = -2.0 * M_PI + lon;
+
+ if (lon < -1.0 * M_PI)
+ lon = 2.0 * M_PI + lon;
+
+ if (lon == -2.0 * M_PI)
+ lon *= -1.0;
+
+ return lon;
+}
+
+/* Converts a latitude value in degrees to radians */
+static double latitudeDegreesToRadians(double lat_degrees)
+{
+ double lat = M_PI * lat_degrees / 180.0;
+ if (lat > 2.0 * M_PI)
+ lat = remainder(lat, 2.0 * M_PI);
+
+ if (lat < -2.0 * M_PI)
+ lat = remainder(lat, -2.0 * M_PI);
+
+ if (lat > M_PI)
+ lat = M_PI - lat;
+
+ if (lat < -1.0 * M_PI)
+ lat = -1.0 * M_PI - lat;
+
+ if (lat > M_PI_2)
+ lat = M_PI - lat;
+
+ if (lat < -1.0 * M_PI_2)
+ lat = -1.0 * M_PI - lat;
+
+ return lat;
+}
+
+/**
+* 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);
+}
+
+/**
+* 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);
+}
+
+/* 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);
+}
+
+/**
+ * Distance functions
+ *
+ **/
+
+/* The haversine formula calculate the distance between two lat/lon points in
meters.
+This formula assumes a spherical model of the earth, which can lead to an
error of about 0.3% compared to a ellipsoidal model.*/
+static double haversine(double lat1, double lon1, double lat2, double lon2)
+{
+ double d_lon = lon2 - lon1;
+ double d_lat = lat2 - lat1;
+ double a = sin(d_lat / 2) * sin(d_lat / 2) + sin(d_lon / 2) * sin(d_lon
/ 2) * cos(lat2) * cos(lat1);
+ double c = 2 * atan2(sqrt(a), sqrt(1 - a));
+ //TODO: Same as the previous line (which one is best?)
+ //double c = 2 * asin(sqrt(a));
+ double r = 6371009;
+ return r * c;
+}
+
+static double geoDistancePointPoint(GEOSGeom a, GEOSGeom b)
+{
+ double lat1_r = 0, lon1_r = 0, lat2_r = 0, lon2_r = 0;
+ pointToRadian(a, &lat1_r, &lon1_r);
+ pointToRadian(b, &lat2_r, &lon2_r);
+ return haversine(lat1_r, lon1_r, lat2_r, lon2_r);
+}
+
+//TODO -> First try in geom_todo.c
+static double calculatePerpendicular(GEOSGeom a, GEOSGeom b) { return INT_MAX;
}
+
+/* Distance between Point and a simple Line (only one Line segment) */
+static double geoDistancePointLineSingle(GEOSGeom point, GEOSGeom line)
+{
+ double distancePerpendicular, distanceStart, distanceEnd;
+
+ /* Calculate perpendicular of point in Line */
+ //TODO Implement this correctly
+ distancePerpendicular = calculatePerpendicular(a, line);
+
+ /* Calculate distance of point to start and end points of line */
+ distanceStart = geoDistancePointPoint(point,
GEOSGeomGetStartPoint(line));
+ distanceEnd = geoDistancePointPoint(point, GEOSGeomGetEndPoint(line));
+
+ /* Determine which of the distances is smaller */
+ if (distanceStart < distancePerpendicular && distanceStart <
distanceEnd)
+ {
+ return distanceStart;
+ }
+ else if (distanceEnd < distancePerpendicular && distanceEnd <
distanceStart)
+ {
+ return distanceEnd;
+ }
+ else
+ {
+ return distancePerpendicular;
+ }
+}
+
+/* Given a Line/LinearRing geometry with multiple segments and the index to
fetch, returns a single Line segment */
+static GEOSGeometry *getLineSegment(GEOSCoordSequence *multiLineCoords, int
segmentIndex)
+{
+ double x1, y1, x2, y2;
+ GEOSCoordSequence *segment;
+ GEOSCoordSeq_getXY(multiLineCoords, segmentIndex, &x1, &y1);
+ GEOSCoordSeq_getXY(multiLineCoords, segmentIndex + 1, &x2, &y2);
+ segment = GEOSCoordSeq_create(2, 2);
+ GEOSCoordSeq_setXY(segment, 0, x1, y1);
+ GEOSCoordSeq_setXY(segment, 1, x2, y2);
+ return GEOSGeom_createLineString(segment);
+}
+
+/* Distance between Point and Line (with multiple line segments) */
+static double geoDistancePointLineMulti(GEOSGeom point, GEOSGeom line, int
lineSegments)
+{
+ const GEOSCoordSequence *gcs = GEOSGeom_getCoordSeq(line);
+ GEOSGeometry *geo_segment;
+ double distance, min_distance = INT_MAX;
+ for (int i = 0; i < lineSegments; i++)
+ {
+ geo_segment = getLineSegment((GEOSCoordSequence *)gcs, i);
+ distance = geoDistancePointLineSingle(point, geo_segment);
+ if (distance < min_distance)
+ min_distance = distance;
+ if (geo_segment != NULL)
+ GEOSGeom_destroy(geo_segment);
+ }
+ return min_distance;
+}
+
+/* Distance between Point and Line */
+static double geoDistancePointLine(GEOSGeom point, GEOSGeom line)
+{
+ int numPoints = GEOSGeomGetNumPoints(line);
+ if (numPoints > 2)
+ {
+ /* For lines with multiple segments, we use a different
function */
+ return geoDistancePointLineMulti(point, line, numPoints - 1);
+ }
+ else
+ {
+ /* For lines with only one segment, we calculate the distance
directly */
+ return geoDistancePointLineSingle(point, line);
+ }
+}
+
+/* Distance between two Lines. */
+static double geoDistanceLineLine(GEOSGeom line1, GEOSGeom line2)
+{
+ int numPoints = GEOSGeomGetNumPoints(line1);
+ GEOSGeometry *linePoint;
+ double distance, min_distance = INT_MAX;
+ for (int i = 0; i < numPoints; i++)
+ {
+ linePoint = GEOSGeomGetPointN(line1, i);
+ distance = geoDistancePointLine(linePoint, line2);
+ if (distance < min_distance)
+ min_distance = distance;
+ if (linePoint != NULL)
+ GEOSGeom_destroy(linePoint);
+ }
+ return min_distance;
+}
+
+//TODO Check if this works always
+//For fast testing, we could use the polygon's minimum bounding box
+static bool pointWithinPolygonRing(GEOSGeom point, GEOSGeom polygon)
+{
+ int intersectionNum = 0, pointsPolygon;
+ double lat_p, lon_p, lat_poly1, lon_poly1, lat_poly2, lon_poly2;
+ GEOSGeometry *segmentPolygon, *intersectionPoints;
+ const GEOSCoordSequence *polygonSeq;
+
+ //Get an point that's outside the polygon
+ //TODO Get the outside point using the polygon's bounding box, instead
of being static
+ double lat_o = 48.193, lon_o = -4.551;
+
+ //Construct a line between the outside point and the input point
+ GEOSGeomGetX(point, &lon_p);
+ GEOSGeomGetY(point, &lat_p);
+ GEOSGeometry *outInLine = geographicPointsToCartesianLine(lat_p, lon_p,
lat_o, lon_o);
+
+ //Count the number of intersections between the polygon and the
constructed line
+ pointsPolygon = GEOSGeomGetNumPoints(polygon);
+ polygonSeq = GEOSGeom_getCoordSeq(polygon);
+ for (int i = 0; i < pointsPolygon; i++)
+ {
+ GEOSCoordSeq_getXY(polygonSeq, i, &lon_poly1, &lat_poly1);
+ GEOSCoordSeq_getXY(polygonSeq, (i + 1) % pointsPolygon,
&lon_poly2, &lat_poly2);
+ segmentPolygon = geographicPointsToCartesianLine(lat_poly1,
lon_poly1, lat_poly2, lon_poly2);
+
+ intersectionPoints = GEOSIntersection(segmentPolygon,
outInLine);
+
+ //If there is an intersection, a point will be returned (line
when there is none)
+ if (GEOSGeomTypeId(intersectionPoints) == GEOS_POINT)
+ {
+ double x_i, y_i, z_i, lat_i, lon_i;
+ GEOSGeomGetX(intersectionPoints, &x_i);
+ GEOSGeomGetY(intersectionPoints, &y_i);
+ GEOSGeomGetZ(intersectionPoints, &z_i);
+ cart2geog(x_i, y_i, z_i, &lon_i, &lat_i);
+ printf("\nPoint Intersection: (%f %f) (%f %f %f)\n\n",
lon_i, lat_i, x_i, y_i, z_i);
+ fflush(stdout);
+ intersectionNum++;
+ }
+
+ if (intersectionPoints != NULL)
+ GEOSGeom_destroy(intersectionPoints);
+ if (segmentPolygon != NULL)
+ GEOSGeom_destroy(segmentPolygon);
+ }
+
+ //If even, the point is not within the polygon. If odd, it is within
+ return intersectionNum % 2 == 1;
+}
+
+/* Distance between Point and Polygon.*/
+static double geoDistancePointPolygon(GEOSGeom point, GEOSGeom polygon)
+{
+ const GEOSGeometry *polygon_ring;
+ polygon_ring = GEOSGetExteriorRing(polygon);
+
+ //Check if point is in polygon
+ if (pointWithinPolygonRing(point, (GEOSGeometry *)polygon_ring))
+ return 0;
+
+ //Compare Point to the various polygon segments
+ return geoDistancePointLine(point, (GEOSGeometry *)polygon_ring);
+}
+
+/* Distance between Line and Polygon. */
+static double geoDistanceLinePolygon(GEOSGeom line, GEOSGeom polygon)
+{
_______________________________________________
checkin-list mailing list
[email protected]
https://www.monetdb.org/mailman/listinfo/checkin-list
