Change units to radians and add some tests (#657 and #1305)

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

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 626f099ffbbf74bf289da5abdbd166d785429da5..cd7918b0b01e9e63678acf5f8bfe521f86059db0 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1793,7 +1793,7 @@ double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
 * @param r - location of first point.
 * @param spheroid - spheroid definition.
 * @param distance - distance, in units of the spheroid def'n.
-* @param azimuth - azimuth in degrees.
+* @param azimuth - azimuth in radians.
 * @return s - location of projected point.
 * 
 */
@@ -1801,18 +1801,16 @@ LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, dou
 {
        GEOGRAPHIC_POINT geo_source, geo_dest;
        POINT4D pt_dest;
-       double azimuth_radians;
        double x, y;
        POINTARRAY *pa;
        LWPOINT *lwp;
 
        /* Check the azimuth validity, convert to radians */
-       if ( azimuth < -360.0 || azimuth > 360.0 ) 
+       if ( azimuth < -2.0 * M_PI || azimuth > 2.0 * M_PI ) 
        {
-               lwerror("Azimuth must be between -360 and 360");
+               lwerror("Azimuth must be between -2PI and 2PI");
                return NULL;
        }
-       azimuth_radians = deg2rad(azimuth);
 
        /* Check the distance validity */
        if ( distance < 0.0 || distance > (M_PI * spheroid->radius) )
@@ -1827,7 +1825,7 @@ LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, dou
        geographic_point_init(x, y, &geo_source);
        
        /* Try the projection */
-       if( spheroid_project(&geo_source, spheroid, distance, azimuth_radians, &geo_dest) == LW_FAILURE ) 
+       if( spheroid_project(&geo_source, spheroid, distance, azimuth, &geo_dest) == LW_FAILURE ) 
        {
                LWDEBUGF(3, "Unable to project from (%g %g) with azimuth %g and distance %g", x, y, azimuth, distance);
                lwerror("Unable to project from (%g %g) with azimuth %g and distance %g", x, y, azimuth, distance);
@@ -1836,22 +1834,23 @@ LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, dou
        
        /* Build the output LWPOINT */
        pa = ptarray_construct(0, 0, 1);
-       pt_dest.x = rad2deg(geo_dest.lon);
-       pt_dest.y = rad2deg(geo_dest.lat);
+       pt_dest.x = rad2deg(longitude_radians_normalize(geo_dest.lon));
+       pt_dest.y = rad2deg(latitude_radians_normalize(geo_dest.lat));
        pt_dest.z = pt_dest.m = 0.0;
        ptarray_set_point4d(pa, 0, &pt_dest);
        lwp = lwpoint_construct(r->srid, NULL, pa);
+       lwgeom_set_geodetic(lwpoint_as_lwgeom(lwp), LW_TRUE);
        return lwp;
 }
 
 
 /**
-* Calculate a projected point given a source point, a distance and a bearing.
+* Calculate a bearing (azimuth) given a source and destination point.
 * @param r - location of first point.
+* @param s - location of second point.
 * @param spheroid - spheroid definition.
-* @param distance - distance, in units of the spheroid def'n.
-* @param azimuth - azimuth in degrees.
-* @return s - location of projected point.
+* @param azimuth - azimuth in radians.
+* @return 
 * 
 */
 double lwgeom_azumith_spheroid(const LWPOINT *r, const LWPOINT *s, const SPHEROID *spheroid)
@@ -1869,8 +1868,14 @@ double lwgeom_azumith_spheroid(const LWPOINT *r, const LWPOINT *s, const SPHEROI
        y2 = lwpoint_get_y(s);
        geographic_point_init(x2, y2, &g2);
        
+       /* Same point, return a standard azimuth instead of NULL or NaN */
+       if ( FP_EQUALS(x1, x2) && FP_EQUALS(y1, y2) )
+       {
+               return 0.0;
+       }
+       
        /* Do the direction calculation */
-       return rad2deg(spheroid_direction(&g1, &g2, spheroid));
+       return spheroid_direction(&g1, &g2, spheroid);
 }
 
 /**

diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index 77467c3abf2c969a9908d9df9d73cfa64698fbaf..a7a463c03232a01bf1294ce81ce45238d9babce6 100644 (file)
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -601,7 +601,7 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS)
 /*
 ** geography_project(GSERIALIZED *g, distance, azimuth)
 ** returns point of projection given start point, 
-** azimuth (bearing) and distance
+** azimuth in radians (bearing) and distance in meters
 */
 PG_FUNCTION_INFO_V1(geography_project);
 Datum geography_project(PG_FUNCTION_ARGS)
@@ -636,8 +636,8 @@ Datum geography_project(PG_FUNCTION_ARGS)
        }
        
        /* Read the other parameters */
-       distance = PG_GETARG_FLOAT8(1);
-       azimuth = PG_GETARG_FLOAT8(2);
+       distance = PG_GETARG_FLOAT8(1); /* Meters */
+       azimuth = PG_GETARG_FLOAT8(2); /* Radians */
 
        /* Initialize spheroid */
        spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
@@ -709,7 +709,7 @@ Datum geography_azimuth(PG_FUNCTION_ARGS)
        /* Calculate the direction */
        azimuth = lwgeom_azumith_spheroid(lwgeom_as_lwpoint(lwgeom1), lwgeom_as_lwpoint(lwgeom2), &s);
 
-       /* Clean up, but not all the way to the point arrays */
+       /* Clean up */
        lwgeom_free(lwgeom1);
        lwgeom_free(lwgeom2);
 

diff --git a/postgis/lwgeom_functions_basic.c b/postgis/lwgeom_functions_basic.c
index 68ae03e7ce17c8d163e3865fce28ef9b8758f25a..62e971d5311f92f76154a98163d5db49cae5267f 100644 (file)
--- a/postgis/lwgeom_functions_basic.c
+++ b/postgis/lwgeom_functions_basic.c
@@ -2325,7 +2325,7 @@ Datum LWGEOM_azimuth(PG_FUNCTION_ARGS)
                lwerror("Error extracting point");
                PG_RETURN_NULL();
        }
-       lwgeom_release((LWGEOM *)lwpoint);
+       lwpoint_free(lwpoint);
        PG_FREE_IF_COPY(geom, 0);
 
        /* Extract second point */
@@ -2349,9 +2349,15 @@ Datum LWGEOM_azimuth(PG_FUNCTION_ARGS)
                lwerror("Error extracting point");
                PG_RETURN_NULL();
        }
-       lwgeom_release((LWGEOM *)lwpoint);
+       lwpoint_free(lwpoint);
        PG_FREE_IF_COPY(geom, 1);
 
+       /* Standard return value for equality case */
+       if ( FP_EQUALS(p1.x, p2.x) && FP_EQUALS(p1.y, p2.y) )
+       {
+               PG_RETURN_FLOAT8(0.0);
+       }
+
        /* Compute azimuth */
        if ( ! azimuth_pt_pt(&p1, &p2, &result) )
        {

diff --git a/regress/tickets.sql b/regress/tickets.sql
index 9dbb675aab05f28d8377b9fe0f5984cf87bc9297..9cbbb3abc3e8662bca87dac8577037806d8ba405 100644 (file)
--- a/regress/tickets.sql
+++ b/regress/tickets.sql
@@ -469,6 +469,16 @@ select '#1344', octet_length(ST_AsEWKB(st_makeline(g))) FROM ( values ('POINT(0
 -- #1385
 SELECT '#1385', ST_Extent(g) FROM ( select null::geometry as g ) as foo; 
 
+-- #657
+SELECT '#657.1',Round(ST_X(ST_Project('POINT(175 10)'::geography, 2000000, 3.1415/2)::GEOMETRY)::numeric,2);
+SELECT '#657.2',Round(ST_Distance(ST_Project('POINT(10 10)'::geography, 10, 0), 'POINT(10 10)'::geography)::numeric,2);
+SELECT '#657.3',ST_DWithin(ST_Project('POINT(10 10)'::geography, 2000, pi()/2), 'POINT(10 10)'::geography, 2000);
+
+-- #1305
+SELECT '#1305.1',ST_AsText(ST_Project('POINT(10 10)'::geography, 0, 0));
+WITH pts AS ( SELECT 'POINT(0 45)'::geography AS s, 'POINT(45 45)'::geography AS e ) 
+SELECT '#1305.2',abs(ST_Distance(e, ST_Project(s, ST_Distance(s, e), ST_Azimuth(s, e)))) < 0.001 FROM pts;
+
 -- Clean up
 DELETE FROM spatial_ref_sys;
 

diff --git a/regress/tickets_expected b/regress/tickets_expected
index 7ac5a649cec610418be112ca349ab2e1eba8b513..5938a0e344bd11d6207bc294e1355c7a9cfdc051 100644 (file)
--- a/regress/tickets_expected
+++ b/regress/tickets_expected
@@ -143,3 +143,8 @@ ERROR:  Geometry type (Polygon) does not match column type (MultiPolygon)
 </#1320>
 #1344|25
 #1385|
+#657.1|-166.78
+#657.2|10.00
+#657.3|t
+#1305.1|POINT(10 10)
+#1305.2|t