How to calculate the Azimuth between two points in

问题:

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?

回答1:

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.