How to calculate the 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?

You've used different coordinates for the second and third examples - looks like you repeated the Y-coordinates. That makes '90' the right answer for Example 3 since the points make a horizontal line. Example 2 I don't understand, there's no coordinate reference so I don't know what ::geography will make of it. Example 1 looks okay. — Spacedman, Aug 27 '14 at 12:14
@Spacedman thank you. I just corrected example 1. nevertheless, I still wanted to know how this is calculated internally in PostGIS, and why example 2 returns a different result. I found some nformation about this [here](http://postgis.net/docs/manual-2.1/using_postgis_dbmanagement.html#PostGIS_Geography). But I think the three should retrn the same result. — joaorodr84, Aug 27 '14 at 12:47
The PostGIS source code is mirrored on github, so a handy place to look is here: https://github.com/postgis/postgis/search?utf8=%E2%9C%93&q=ST_Azimuth but its a complex beast so hard to track down exactly where the numbers are crunched. — Spacedman, Aug 27 '14 at 12:52

John Powell · Accepted Answer · 2014-08-27T13:46:47.500

To answer your second question first, there are basically two types of azimuth calculation, those done on planar coordinates, and those done on a sphere or spheroid, as is the case with the geography datatype. The wikipedia azimuth article explains this well.

In the first case the calculation is extremely simple, being essentially just the atan between the two pairs of points. You can see the source code on github in the method azimuth_pt_pt.

For the 2nd case, the calculation is a bit more complex, being on a sphere, and you can find the actual calculation again on github in the method spheroid_direction.

The reason 2 is different is because you are using the geography datatype which should be in the range [-180,180],[-90,90]. Actually, I am surprised it doesn't give an error.

Examples 1 and 3 are essentially the same, as you have simply switched from one coordinate reference system to another, so the relative direction of the points is unchanged.

There is no correct way to do it, but if you want to use geodetic coordinates, then use the geography datatype.

Note there is a difference between:

select degrees(
   st_azimuth(
       st_makepoint(0, 0)::geography, 
       st_makepoint(45, 45)::geography)
);

which yields 35.4100589051161

while,

select degrees(
    st_azimuth(
       st_setsrid(st_makepoint(0, 0),4326),
       st_setsrid(st_makepoint(45, 45),4326))
);

yields 45. One is doing point to point azimuth on the plane, while the other is doing the spheroid azimuth, as suggested by their function names.

Thank you very much for your detailed answer. I'll analyse all the information you provided. — joaorodr84, Aug 27 '14 at 14:06
Would you please take a look at this other [post](http://stackoverflow.com/questions/25529222/is-the-azimuth-on-the-equator-equal-to-the-one-not-on-the-equator) of mine? Thanks :) — joaorodr84, Aug 27 '14 at 14:34

How to calculate the Azimuth between two points in PostGIS?

1 Answers1

Linked