Cleaned up distance_sphere and updated documentation.

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

diff --git a/CHANGES b/CHANGES
index 1acbe9d3fbf96d1eb101263efbd35bdbd38a08ee..74ad6ac3e258921d911b4adbf6bf082d1d9b3882 100644 (file)
--- a/CHANGES
+++ b/CHANGES
@@ -11,6 +11,7 @@ PostGIS 1.0.0
        - build scripts refinements
        - MultiLine handling BUG fix in loader and dumper
        - improved version handling, central Version.config
+        - added distance_sphere function 
 
 PostGIS 1.0.0RC3
 2005/02/24

diff --git a/doc/postgis.xml b/doc/postgis.xml
index eecaa67112b714a47d55f04d17b391a8ef337344..bec1db92e2d4ffe6e7ac5a93b58023596f598efa 100644 (file)
--- a/doc/postgis.xml
+++ b/doc/postgis.xml
@@ -3834,7 +3834,8 @@ dimension
             <para>Returns linear distance in meters between two lat/lon
             points. Uses a spherical earth and radius of 6370986 meters.
             Faster than <link linkend="distance_spheroid">distance_spheroid()</link>, but
-            less accurate.</para>
+            less accurate.
+                       Only implemented for points.</para>                     
           </listitem>
         </varlistentry>
 

diff --git a/lwgeom/lwgeom_spheroid.c b/lwgeom/lwgeom_spheroid.c
index 8a302917fdc8bbaa42e4d09a86802269edda7d3c..6b9f85f95ae05ffe39a1a34ae5e004678f166a23 100644 (file)
--- a/lwgeom/lwgeom_spheroid.c
+++ b/lwgeom/lwgeom_spheroid.c
@@ -547,16 +547,12 @@ Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS)
        double sino;
         double result;
 
-        //elog(NOTICE, "LWGEOM_distance_sphere called\n");
-
        if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2))
        {
                elog(ERROR, "LWGEOM_distance_sphere Operation on two GEOMETRIES with differenc SRIDs\n");
                PG_RETURN_NULL();
        }
 
-        //elog(NOTICE, "LWGEOM_distance_sphere passed srid check: %i, %i\n", pglwgeom_getSRID(geom1), pglwgeom_getSRID(geom2));
-
         lwpt1 = lwpoint_deserialize(SERIALIZED_FORM(geom1));
        if (lwpt1 == NULL)
        {
@@ -571,20 +567,12 @@ Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS)
                PG_RETURN_NULL();
        }
 
-        //elog(NOTICE, "LWGEOM_distance_sphere passed point check\n");
-
         pt1 = palloc(sizeof(POINT2D));
         pt2 = palloc(sizeof(POINT2D));
 
         lwpoint_getPoint2d_p(lwpt1, pt1);
         lwpoint_getPoint2d_p(lwpt2, pt2);
 
-        //elog(NOTICE, "LWGEOM_distance_sphere created point2d objects: %f, %f %f, %f\n", pt1->x, pt1->y, pt2->x, pt2->y);
-
-        //elog(NOTICE, "LWGEOM_distance_sphere beginning calculation\n");
-        //elog(NOTICE, "LWGEOM_distance_sphere TWO_PI = %f\n", TWO_PI);
-        //elog(NOTICE, "LWGEOM_distance_sphere M_PI = %f\n", M_PI);
-
        /*
         * Start geo_distance code.  Longitude is degrees west of
         * Greenwich, and thus is negative from what normal things
@@ -596,8 +584,6 @@ Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS)
        long2 = -2 * (pt2->x / 360.0) * M_PI;
        lat2 = 2 * (pt2->y / 360.0) * M_PI;
 
-        //elog(NOTICE, "LWGEOM_distance_sphere normalized: %f, %f %f, %f", long1, lat1, long2, lat2);
-
         /* compute difference in longitudes - want < 180 degrees */
         longdiff = fabs(long1 - long2);
         if (longdiff > M_PI)
@@ -605,8 +591,6 @@ Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS)
                 longdiff = (2 * M_PI) - longdiff;
         }
 
-        //elog(NOTICE, "LWGEOM_distance_sphere longdiff = %f\n", 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.)
@@ -614,11 +598,8 @@ Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS)
                 sino = 1.;
         }
 
-        //elog(NOTICE, "LWGEOM_distance_sphere sino = %f\n", sino);          
-
         result = 2. * EARTH_RADIUS * asin(sino);
         
-        //elog(NOTICE, "LWGEOM_distance_sphere finished calculation: %f\n", result);
         pfree(pt1);
         pfree(pt2);
 

diff --git a/lwgeom/lwpostgis.sql.in b/lwgeom/lwpostgis.sql.in
index e2fb60a75f27335c5a35db636e18f7de1596de9d..5fa00b0c566796ab0e4e2d4dae99e61c6f53a691 100644 (file)
--- a/lwgeom/lwpostgis.sql.in
+++ b/lwgeom/lwpostgis.sql.in
@@ -1558,41 +1558,11 @@ CREATEFUNCTION distance_spheroid(geometry,geometry,spheroid)
        AS '@MODULE_FILENAME@','LWGEOM_distance_ellipsoid_point'
        LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict);
 
-CREATEFUNCTION distance_sphere_latlon(geometry,geometry)
+CREATEFUNCTION distance_sphere(geometry,geometry)
        RETURNS FLOAT8
        AS '@MODULE_FILENAME@','LWGEOM_distance_sphere'
        LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict);
 
-CREATEFUNCTION distance_sphere(geometry,geometry)
-        RETURNS FLOAT8 AS
-'
-DECLARE
-  proj1 integer;
-  proj2 integer;
-  projcount integer;
-BEGIN
-  SELECT INTO proj1 SRID($1);
-  SELECT INTO proj2 SRID($2);
-  IF proj1 != proj2 THEN
-    RAISE EXCEPTION ''Geometry projections must match.'';
-  END IF;
-  IF proj1 != -1 AND proj2 != -1 THEN
-    SELECT INTO projcount count(*) FROM spatial_ref_sys 
-    WHERE srid = proj1 AND proj4text like ''%longlat%'';
-    IF projcount = 0 THEN
-      RAISE EXCEPTION ''Geometries must be in a long/lat format.'';
-    END IF;
-    SELECT INTO projcount count(*) FROM spatial_ref_sys
-    WHERE srid = proj2 AND proj4text like ''%longlat%'';
-    IF projcount = 0 THEN
-      RAISE EXCEPTION ''Geometries must be in a long/lat format.'';
-    END IF;
-  END IF;
-  RETURN distance_sphere_latlon($1, $2);
-END;
-'
-LANGUAGE 'plpgsql' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable);
-
 -- Minimum distance. 2d only.
 CREATEFUNCTION distance(geometry,geometry)
        RETURNS float8