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
_______________________________________________ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users