Clever Geek Handbook

How to calculate azimuth between two points in PostGIS?

I have a GIS API based on OpenLayers. I am trying to do an azimuth calculation in JavaScript and I needed to calculate the azimuth somehow in order to run the tests.

I started using PostGIS, but there seem to be many ways to calculate the azimuth between two points. I am showing you three of them and some of them 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 azimuth in PostGIS? Of course, I want to look at geodetic coordinates.

Is there a way to find out how PostGIS implements computation? I mean, can you see how the "ST_Azimuth" function is implemented?

0


source to share


1 answer


First of all, to answer your second question, there are basically two types of azimuthal calculations done in planar coordinates and those done on a sphere or spheroid, as in the case of a type of geography. The apimuth wikipedia explains it well.

In the first case, the calculation is extremely simple, essentially only atan

between two pairs of points. You can see the githubazimuth_pt_pt

source code in the method .

For the second case, the computation is a bit more complicated while on the sphere and you can find the actual computation again on github in the methodspheroid_direction

.

Reason 2 is different in that you are using the geography type, which must be in the range [-180,180], [- 90,90]. In fact, I'm surprised that it doesn't give an error.

Examples 1 and 3 are essentially the same since you just switched from one coordinate system to another, so the relative direction of the points did not change.

There is no correct way to do this, but if you want to use geodetic coordinates

, use geography datatype

.



Please note that there is a difference between:

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

which gives 35.4100589051161

then,

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

gives 45 . One of them points to a point azimuth on the plane, and the other performs the azimuth of a spheroid, as suggested by their function names.

+3


source






More articles:


All Articles