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?

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

<br><a href="http://www.thescripts.com/forum/thread452165.html">http://www.thescripts.com/forum/thread452165.html</a><br><br>Here's the code:<br><br>CREATE or replace FUNCTION get_bearing(p1 geometry, p2 geometry) RETURNS double precision AS $$

<br>DECLARE<br><br>a_x float8;<br>a_y float8;<br>b_x float8;<br>b_y float8;<br>res float8;<br><br>BEGIN<br><br>--special case - identical geometries, return 0.<br>if Equals(p1, p2) then<br>    RETURN 0;<br>end if; <br><br>

--reference vector representing north from origin <br>b_x := 0;<br>b_y := 1;<br><br>--second vector, represents the translation of p1 to p2<br>a_x :=  x(p2) - x(p1);<br>a_y :=  y(p2) - y(p1);<br><br>--special case: x is less than 0. 

<br>--the smallest angle between the vectors is calculated,<br>--therefore need to subtract this angle from 360<br>if a_x < 0 then <br>    return 360-((acos((a_x*b_x+a_y*b_y)/sqrt(a_x*a_x+a_y*a_y))/pi())*180);<br>else<br>

    return (acos((a_x*b_x+a_y*b_y)/sqrt(a_x*a_x+a_y*a_y))/pi())*180;<br>end if;<br><br>END;<br>$$ LANGUAGE plpgsql;<br><br>Cheers,<br><br>Will Temperley