[postgis-users] The simple task of calculating a bearing from point a to point b

[William Temperley]

Dear all,

Ok, it has taken me ages to find a solution to what should be a simple
task.  I have written a pg/plsql function to do this, and I'd like to know
where I've gone wrong, if I have. I'm not a mathematician.

I think the function will work with conformal coordinate systems, however
I'd like to know which for cases/projections/coordinate systems it wouldn't
work.

Will it work with wgs84? Or would I have to project into a planar coordinate
system before doing so?

The trig bit of this function is based
on:<http://www.thescripts.com/forum/thread452165.html>
http://en.wikipedia.org/wiki/Dot_product
And:
http://www.thescripts.com/forum/thread452165.html

Here's the code:

CREATE or replace FUNCTION get_bearing(p1 geometry, p2 geometry) RETURNS
double precision AS $$
DECLARE

a_x float8;
a_y float8;
b_x float8;
b_y float8;
res float8;

BEGIN

--special case - identical geometries, return 0.
if Equals(p1, p2) then
    RETURN 0;
end if;

--reference vector representing north from origin
b_x := 0;
b_y := 1;

--second vector, represents the translation of p1 to p2
a_x :=  x(p2) - x(p1);
a_y :=  y(p2) - y(p1);

--special case: x is less than 0.
--the smallest angle between the vectors is calculated,
--therefore need to subtract this angle from 360
if a_x < 0 then
    return 360-((acos((a_x*b_x+a_y*b_y)/sqrt(a_x*a_x+a_y*a_y))/pi())*180);
else
    return (acos((a_x*b_x+a_y*b_y)/sqrt(a_x*a_x+a_y*a_y))/pi())*180;
end if;

END;
$$ LANGUAGE plpgsql;

Cheers,

Will Temperley
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20071121/06418b2b/attachment.html>