From f1ca2a8e8ee2c82a471c3ee7c275069a62c7e902 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Tue, 3 Jan 2012 19:21:32 +0000 Subject: [PATCH] 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 --- liblwgeom/lwgeodetic.c | 31 ++++++++++++++++++------------- postgis/geography_measurement.c | 8 ++++---- postgis/lwgeom_functions_basic.c | 10 ++++++++-- regress/tickets.sql | 10 ++++++++++ regress/tickets_expected | 5 +++++ 5 files changed, 45 insertions(+), 19 deletions(-) diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 626f099ff..cd7918b0b 100644 --- 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 77467c3ab..a7a463c03 100644 --- 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 68ae03e7c..62e971d53 100644 --- 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 9dbb675aa..9cbbb3abc 100644 --- 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 7ac5a649c..5938a0e34 100644 --- a/regress/tickets_expected +++ b/regress/tickets_expected @@ -143,3 +143,8 @@ ERROR: Geometry type (Polygon) does not match column type (MultiPolygon) #1344|25 #1385| +#657.1|-166.78 +#657.2|10.00 +#657.3|t +#1305.1|POINT(10 10) +#1305.2|t -- 2.40.0