Как рассчитать азимут между двумя точками в PostGIS?

Question

У меня есть ГИС API на основе OpenLayers. Я пытался реализовать вычисление азимута в JavaScript, и мне нужно было каким-то образом рассчитать азимут для выполнения тестов.

Я начал использовать PostGIS, но, похоже, есть много способов вычислить азимут между двумя точками. Я покажу вам три из них, и некоторые из них дают разные результаты.

-- 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))
)

Как правильно рассчитать азимут в PostGIS? Конечно, я хочу рассмотреть геодезические координаты.

Есть ли способ узнать, как PostGIS выполняет расчеты? Я имею в виду, можно ли увидеть способ реализации функции "ST_Azimuth"?

Answer 1

Чтобы ответить на ваш второй вопрос в первую очередь, есть в основном два типа вычисления азимута: те, которые выполняются на плоских координатах, и те, которые выполняются на сфере или сфероиде, как в случае с типом данных географии. В азимутальной статье Википедии это хорошо объясняется.

В первом случае расчет предельно прост, будучи по сути просто atan между двумя парами точек. Вы можете увидеть исходный код на github в методе azimuth_pt_pt,

Во втором случае вычисления немного сложнее, они находятся на сфере, и вы можете снова найти фактический расчет на github в методе. spheroid_direction,

Причина 2 в том, что вы используете тип данных geography, который должен находиться в диапазоне [-180,180],[-90,90]. На самом деле, я удивлен, что это не дает ошибку.

Примеры 1 и 3 практически одинаковы, так как вы просто переключились с одной системы координат на другую, поэтому относительное направление точек не изменилось.

Нет правильного способа сделать это, но если вы хотите использовать geodetic coordinates затем используйте geography datatype,

Обратите внимание, что есть разница между:

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

что дает 35,4100589051161

в то время как,

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

дает 45. Один выполняет указание азимута на плоскости, в то время как другой выполняет азимут сфероида, как следует из названий их функций.