Calculating Azimuth between two points in PostGIS?

Question

I have a GIS API based on OpenLayers. I have been trying to implement Azimuth calculation in JavaScript and I needed some way to calculate the Azimuth in order to perform tests.

I started using PostGIS, but there seem to be many ways to calculate the Azimuth between two points. I show you three of them and, some return different results.

-- Example 1 - Result 90

select st_azimuth(
st_transform(st_geomfromtext('POINT(-81328.998084106 7474929.8690234)', 900913), 4326),
st_transform(st_geomfromtext('POINT(4125765.0381464 7474929.8690234)', 900913), 4326)
)/pi()*180

-- Example 2 - Result 155.692090425822
SELECT degrees(
ST_Azimuth(
ST_MakePoint(-81328.998084106, 7474929.8690234)::geography,
ST_MakePoint(4125765.0381464, 7474929.8690234)::geography)
)

-- Example 3 - Result 90

SELECT degrees(
ST_Azimuth(
ST_MakePoint(-81328.998084106, 7474929.8690234),
ST_MakePoint(4125765.0381464, 7474929.8690234))
)

What is the correct way to calculate the Azimuth in PostGIS?

Of course, I want to consider geodesic coordinates.

Is there a way to know how PostGIS realizes the calculations?

I mean, is it possible to see the way the "ST_Azimuth" function is implemented?

Answer 1

For an answer on the geodesic, your example 2 is the closest, but when you ran it you should have seen a couple warnings,

NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY
NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY

because your input coordinates are not in lon/lat, they are in web mercator. So, hybridizing your example 1 and 2, we get,

with pts as (
  select 
    st_transform(st_setsrid(st_makepoint(-81328.998084106,7474929.8690234), 900913), 4326)::geography as a,
    st_transform(st_setsrid(st_makepoint(4125765.0381464,7474929.8690234), 900913), 4326)::geography as b
)
select st_astext(a) as a, st_astext(b) as b, degrees(st_azimuth(a,b)) from pts;

which returns (less some extraneous digits)

        a             |           b          | degrees     
----------------------+----------------------+--------
 POINT(-0.730 55.578) | POINT(37.062 55.578) | 74.231

and perhaps makes you wonder why, if a and b are at the same latitude the azimuth is not 90? The answer is that on the geodesic shortest paths are great circles, so really it's the same reason that flights from Seattle to London go over Greenland, but not over the Great Lakes. (Any west/east trip in the northern hemisphere will bend to the north, and the bend will be bigger the further north you go. Vice versa in the southern hemisphere.)