From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Wed, 21 Dec 2011 19:07:02 +0000 (+0000)
Subject: ST_Azimuth on the spheroid (#1305)
X-Git-Tag: 2.0.0alpha1~379
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=67196edc607948c58e144f4eb3da4c13374ac968;p=postgis

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 5a5642a71..eaf9ba5ce 100644
--- 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 2935c2e2c..626f099ff 100644
--- 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 c29b2f506..5df8ab475 100644
--- 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 9825635f9..77467c3ab 100644
--- 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);
+}
+
+