ST_Azimuth on the spheroid (#1305)

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

diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index 5a5642a7105675f7cc19ce7be90b97f15a33f9cd..eaf9ba5ce259127de0ec4e840cf6399d0f081aef 100644 (file)
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -1322,10 +1322,15 @@ extern void spheroid_init(SPHEROID *s, double a, double b);
 extern double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance);
 
 /**
-* Calculate the location of a point, give a start point, bearing and distance.
+* Calculate the location of a point on a spheroid, give a start point, bearing and distance.
 */
 extern LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, double distance, double azimuth);
 
+/**
+* Calculate the bearing between two points on a spheroid.
+*/
+extern double lwgeom_azumith_spheroid(const LWPOINT *r, const LWPOINT *s, const SPHEROID *spheroid);
+
 /**
 * Calculate the geodetic area of a lwgeom on the sphere. The result
 * will be multiplied by the average radius of the supplied spheroid.

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 2935c2e2c232fea9def250d2737e3486c9525d4e..626f099ffbbf74bf289da5abdbd166d785429da5 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1844,6 +1844,35 @@ LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, dou
        return lwp;
 }
 
+
+/**
+* Calculate a projected point given a source point, a distance and a bearing.
+* @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.
+* @return s - location of projected point.
+* 
+*/
+double lwgeom_azumith_spheroid(const LWPOINT *r, const LWPOINT *s, const SPHEROID *spheroid)
+{
+       GEOGRAPHIC_POINT g1, g2;
+       double x1, y1, x2, y2;
+
+       /* Convert r to a geodetic point */
+       x1 = lwpoint_get_x(r);
+       y1 = lwpoint_get_y(r);
+       geographic_point_init(x1, y1, &g1);
+
+       /* Convert s to a geodetic point */
+       x2 = lwpoint_get_x(s);
+       y2 = lwpoint_get_y(s);
+       geographic_point_init(x2, y2, &g2);
+       
+       /* Do the direction calculation */
+       return rad2deg(spheroid_direction(&g1, &g2, spheroid));
+}
+
 /**
 * Calculate the distance between two LWGEOMs, using the coordinates are
 * longitude and latitude. Return immediately when the calulated distance drops

diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c
index c29b2f506884ff15f48301efb0e003150443a9df..5df8ab475345e7f6b1ef2ddbc99d54a1b2534c43 100644 (file)
--- a/postgis/geography.sql.in.c
+++ b/postgis/geography.sql.in.c
@@ -677,6 +677,13 @@ CREATE OR REPLACE FUNCTION ST_Project(geog geography, distance float8, azimuth f
        LANGUAGE 'C' IMMUTABLE STRICT
        COST 100;
 
+-- Availability: 2.0.0
+CREATE OR REPLACE FUNCTION ST_Azimuth(geog1 geography, geog2 geography)
+       RETURNS float8
+       AS 'MODULE_PATHNAME','geography_azimuth'
+       LANGUAGE 'C' IMMUTABLE STRICT
+       COST 100;
+
 -- Availability: 2.0.0
 CREATE OR REPLACE FUNCTION ST_Perimeter(geog geography, use_spheroid boolean DEFAULT true)
        RETURNS float8

diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index 9825635f9b05d8d7ae3ccfeccb6fa2b580043a07..77467c3abf2c969a9908d9df9d73cfa64698fbaf 100644 (file)
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -36,6 +36,7 @@ Datum geography_covers(PG_FUNCTION_ARGS);
 Datum geography_bestsrid(PG_FUNCTION_ARGS);
 Datum geography_perimeter(PG_FUNCTION_ARGS);
 Datum geography_project(PG_FUNCTION_ARGS);
+Datum geography_azimuth(PG_FUNCTION_ARGS);
 
 /*
 ** geography_distance(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
@@ -660,3 +661,61 @@ Datum geography_project(PG_FUNCTION_ARGS)
        PG_RETURN_POINTER(g_out);
 }
 
+
+/*
+** geography_azimuth(GSERIALIZED *g1, GSERIALIZED *g2)
+** returns direction between points (north = 0)
+** azimuth (bearing) and distance
+*/
+PG_FUNCTION_INFO_V1(geography_azimuth);
+Datum geography_azimuth(PG_FUNCTION_ARGS)
+{
+       LWGEOM *lwgeom1 = NULL;
+       LWGEOM *lwgeom2 = NULL;
+       GSERIALIZED *g1 = NULL;
+       GSERIALIZED *g2 = NULL;
+       double azimuth;
+       SPHEROID s;
+       uint32_t type1, type2;
+
+       /* Get our geometry object loaded into memory. */
+       g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+       
+       /* Only return for points. */
+    type1 = gserialized_get_type(g1);
+    type2 = gserialized_get_type(g2);
+    if ( type1 != POINTTYPE || type2 != POINTTYPE )
+    {
+               elog(ERROR, "ST_Azimuth(geography, geography) is only valid for point inputs");
+               PG_RETURN_NULL();
+    }
+       
+       lwgeom1 = lwgeom_from_gserialized(g1);
+       lwgeom2 = lwgeom_from_gserialized(g2);
+
+       /* EMPTY things cannot be used */
+       if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
+       {
+               lwgeom_free(lwgeom1);
+               lwgeom_free(lwgeom2);
+               elog(ERROR, "ST_Azimuth(geography, geography) cannot work with empty points");
+               PG_RETURN_NULL();
+       }
+       
+       /* Initialize spheroid */
+       spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
+
+       /* 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 */
+       lwgeom_free(lwgeom1);
+       lwgeom_free(lwgeom2);
+
+       PG_FREE_IF_COPY(g1, 0);
+       PG_FREE_IF_COPY(g2, 1);
+       PG_RETURN_FLOAT8(azimuth);
+}
+
+