What algorithms are internally calculated in PostGIS? Can you explain it with simple code? (I want to know the standard of true north.)
There are two different
ST_Azimuth functions for
geography types, which behave differently.
Here are their similarities:
- Determines a bearing angle between two points, such that if you wanted to travel from point A to point B, the angle you'd set your compass to is the azimuth
- North is at zero (like a compass)
- Positive angles are clockwise (also like a compass)
- Angle units are in radians (not degrees; use
Now to their differences for the two functions.
Simple. It calculates the angle between two points in Cartesian space. You can calculate this using atan2, just be sure to account for the differences between the navigational vs mathematical framing of azimuth (i.e. positive CW vs CCW, zero north vs zero east). An angle of zero is grid north—not true north. (Determining true north from grid north at a location is a different question.)
ST_Azimuth(geometry, geometry) is best used for flat projected maps (e.g. UTM), but will give the wrong answer for global-scale positions. Here's a quick demo of some obvious angles:
SELECT degrees(ST_Azimuth(origin, north)) AS a_north, degrees(ST_Azimuth(origin, northeast)) AS a_northeast, degrees(ST_Azimuth(origin, east)) AS a_east FROM ( SELECT 'POINT(0 0)'::geometry AS origin, 'POINT(0 1)'::geometry AS north, 'POINT(1 1)'::geometry AS northeast, 'POINT(1 0)'::geometry AS east ) d; a_north | a_northeast | a_east ---------+-------------+-------- 0 | 45 | 90 (1 row)
Azimuth angles are determined from the first point to the second point along a geodesic on a ellipsoid of revolution (e.g. WGS 84). True north is always to the top rotational axis. This is a much more complicated calculation, as the azimuth angle normally changes while moving along the path between the points (unlike rhumb lines, which have a constant azimuth). Think of an airplane traversing the Atlantic—it arcs north—then south. Here are the azimuths to two airports on each side of the Atlantic Ocean:
SELECT degrees(ST_Azimuth(JFK, LHR)) AS jfk_to_lhr, degrees(ST_Azimuth(LHR, JFK)) AS lhr_to_jfk FROM ( SELECT 'POINT(-73.778889 40.639722)'::geography AS JFK, 'POINT(-0.461389 51.4775)'::geography AS LHR ) d; jfk_to_lhr | lhr_to_jfk ------------------+------------------ 51.3728787133603 | 287.971307265745 (1 row)
You may note that they are not symmetrical if you assume they are on a flat piece of paper (i.e. they don't have opposite angles to each other in Cartesian space).
Internally, PostGIS uses a few algorithms from GeographicLib for this calculation.