Added the distance_sphere function to calculate the distance between two points

on an earth-sized sphere using an algorithm implemented by Bruno Wolff III.
Added the postgresql loader function.

git-svn-id: http://svn.osgeo.org/postgis/trunk@828 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/CREDITS b/CREDITS
index ee2d8c949e777e906d146c9ce764a98890c4ab77..a8ff797cb5540edb61b1acdb40c56bf133311378 100644 (file)
--- a/CREDITS
+++ b/CREDITS
@@ -40,3 +40,4 @@ Version 0.8 of PostGIS includes patches from:
 Version 0.9 of PostGIS includes patches from:
   Klaus Foerster <klaus@svg.cc> on AsSVG()
   Olivier Courtin <pnine@free.fr> on AsSVG()
+  Bruno Wolff III on the distance_sphere algorithm

diff --git a/postgis.h b/postgis.h
index efe2962994f42993232cc615afa177bf552b591b..dd4055588815a7c2ae1e9f00d2570295df75cc89 100644 (file)
--- a/postgis.h
+++ b/postgis.h
@@ -11,6 +11,11 @@
  *
  **********************************************************************
  * $Log$
+ * Revision 1.51  2004/09/16 15:50:59  mleslie
+ * Added the distance_sphere function to calculate the distance between two points
+ * on an earth-sized sphere using an algorithm implemented by Bruno Wolff III.
+ * Added the postgresql loader function.
+ *
  * Revision 1.50  2004/08/19 06:15:58  strk
  * USE_VERSION gets 80 where it got 75
  *
@@ -585,6 +590,7 @@ Datum ellipsoid_in(PG_FUNCTION_ARGS);
 Datum length_ellipsoid(PG_FUNCTION_ARGS);
 Datum length3d_ellipsoid(PG_FUNCTION_ARGS);
 Datum distance_ellipsoid(PG_FUNCTION_ARGS);
+Datum distance_sphere(PG_FUNCTION_ARGS);
 
 Datum point_inside_circle(PG_FUNCTION_ARGS);
 Datum distance(PG_FUNCTION_ARGS);

diff --git a/postgis.sql.in b/postgis.sql.in
index 8c0930897f71663784d65be1e18418458d7ee982..c6fb4310c1277fed16590e15d11948e7fa82e912 100644 (file)
--- a/postgis.sql.in
+++ b/postgis.sql.in
@@ -1301,6 +1301,11 @@ CREATEFUNCTION distance_spheroid(geometry,geometry,spheroid)
        AS '@MODULE_FILENAME@','distance_ellipsoid'
        LANGUAGE 'C' WITH (isstrict);
 
+CREATEFUNCTION distance_sphere(geometry,geometry)
+       RETURNS FLOAT8
+       AS '@MODULE_FILENAME@', 'distance_sphere'
+       LANGUAGE 'C' WITH (isstrict);
+
 --
 -- Generic operations
 --

diff --git a/postgis_proj.c b/postgis_proj.c
index f4560444d51375ff40b6b7eb1bd324ad4d2aa1c4..058a8e789844f408acb27a575537b6ac06e78c64 100644 (file)
--- a/postgis_proj.c
+++ b/postgis_proj.c
@@ -11,6 +11,11 @@
  *
  **********************************************************************
  * $Log$
+ * Revision 1.14  2004/09/16 15:50:59  mleslie
+ * Added the distance_sphere function to calculate the distance between two points
+ * on an earth-sized sphere using an algorithm implemented by Bruno Wolff III.
+ * Added the postgresql loader function.
+ *
  * Revision 1.13  2004/04/28 22:26:02  pramsey
  * Fixed spelling mistake in header text.
  *
@@ -418,6 +423,86 @@ Datum length3d_ellipsoid(PG_FUNCTION_ARGS)
        PG_RETURN_FLOAT8(dist);
 }
 
+/*
+ * This algorithm was taken from the geo_distance function of the 
+ * earthdistance package contributed by Bruno Wolff III.
+ * It was altered to accept GEOMETRY objects and return results in
+ * meters.
+ */
+PG_FUNCTION_INFO_V1(distance_sphere);
+Datum distance_sphere(PG_FUNCTION_ARGS)
+{
+               const double EARTH_RADIUS = 6370986.884258304;
+               const double TWO_PI = 2.0 * M_PI;
+               GEOMETRY                *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+               GEOMETRY                *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+               POINT3D                 *pt1, *pt2;
+               int32                   *offsets1;
+               int32                   *offsets2;
+               char                    *o;
+
+               double                  long1, lat1, long2, lat2;
+               double                  longdiff;
+               double                  sino;
+
+       if (geom1->SRID != geom2->SRID)
+       {
+               elog(ERROR, "optimistic_overlap:Operation on two GEOMETRIES with differenc SRIDs\n");
+               PG_RETURN_NULL();
+       }
+
+       if (geom1->type != POINTTYPE)
+       {
+               elog(ERROR, "optimistic_overlap: first arg isnt a point\n");
+               PG_RETURN_NULL();
+       }
+
+       if (geom2->type != POINTTYPE)
+       {
+               elog(ERROR, "optimistic_overlap: second arg isnt a point\n");
+               PG_RETURN_NULL();
+       }
+
+       offsets1 = (int32 *) ( ((char *) &(geom1->objType[0] ))+ sizeof(int32) * geom1->nobjs );
+       offsets2 = (int32 *) ( ((char *) &(geom2->objType[0] ))+ sizeof(int32) * geom2->nobjs );
+       o = (char *) geom1 + offsets1[0];
+       pt1 = (POINT3D *) o;
+
+       o = (char *) geom2 + offsets2[0];
+       pt2 = (POINT3D *) o;
+
+       /*
+        * Start geo_distance code.  Longitude is degrees west of
+        * Greenwich, and thus is negative from what normal things
+        * will supply the function.
+        */
+        long1 = -1 * (pt1->x / 360.0) * TWO_PI;
+        lat1 = (pt1->y / 360.0) * TWO_PI;
+
+       long2 = -1 * (pt2->x / 360.0) * TWO_PI;
+       lat2 = (pt2->y / 360.0) * TWO_PI;
+
+        /* compute difference in longitudes - want < 180 degrees */
+        longdiff = fabs(long1 - long2);
+        if (longdiff > M_PI)
+                longdiff = TWO_PI - longdiff;
+
+        sino = sqrt(sin(fabs(lat1 - lat2) / 2.) * sin(fabs(lat1 - lat2) / 2.) +
+                cos(lat1) * cos(lat2) * sin(longdiff / 2.) * sin(longdiff / 2.));
+        if (sino > 1.)
+                sino = 1.;
+        PG_RETURN_FLOAT8(2. * EARTH_RADIUS * asin(sino));
+
+/*     PG_RETURN_FLOAT8(distance_sphere_method(pt1->y*M_PI/180.0 ,
+                                               pt1->x*M_PI/180.0 ,
+                                               pt2->y*M_PI/180.0 ,
+                                               pt2->x*M_PI/180.0 ,
+                                               sphere));
+*/
+
+}
+
+
 //distance (geometry,geometry, sphere)
 // -geometrys MUST be points
 PG_FUNCTION_INFO_V1(distance_ellipsoid);