Как рассчитать азимут (угол на север) между двумя координатами WGS84

Question

У меня есть две координаты WGS84, широта и долгота в градусах. Эти точки довольно близко друг к другу, например, на расстоянии всего одного метра.

Есть ли простой способ рассчитать азимут линии между этими точками, то есть угол к северу?

Наивным подходом было бы предположить, что декартова система координат (потому что эти точки расположены так близко друг к другу) и просто использовать

sin (a) = abs (L2-L1) / sqrt (sqr (L2-L1) + sqr (B2-B1))

а = азимут L1, L2 = долгота B1, B2 = широта

Ошибка будет больше по мере удаления координат от экватора, потому что там расстояние между двумя продольными градусами становится все меньше, чем расстояние между двумя широтными градусами (которое остается постоянным).

Я нашел несколько довольно сложных формул, которые на самом деле не хочу реализовывать, потому что они кажутся излишними для точек, которые находятся так близко друг к другу, и мне не нужна очень высокая точность (достаточно двух десятичных знаков, одна, вероятно, подойдет либо потому, что есть другие факторы, которые в любом случае снижают точность, например, тот, который возвращает GPS).

Может быть, я мог бы просто определить приблизительный коэффициент продольной коррекции в зависимости от широты и использовать что-то вроде этого:

sin (a) = abs (L2 * f-L1 * f) / sqrt (sqr (L2 * f-L1 * f) + sqr (B2-B1))

где f - поправочный коэффициент

Есть подсказки?

(Я не хочу использовать для этого какие-либо библиотеки, особенно те, которые не требуют лицензий времени выполнения. Любой исходный код Delphi MPLed подойдет.)

Answer 1

Формулы, на которые вы ссылаетесь в тексте, предназначены для расчета расстояния большого круга между двумя точками. Вот как я вычисляю угол между точками:

uses Math, ...;
...

const
  cNO_ANGLE=-999;

...

function getAngleBetweenPoints(X1,Y1,X2,Y2:double):double;
var
  dx,dy:double;
begin
  dx := X2 - X1;
  dy := Y2 - Y1;

  if (dx > 0) then  result := (Pi*0.5) - ArcTan(dy/dx)   else
  if (dx < 0) then  result := (Pi*1.5) - ArcTan(dy/dx)   else
  if (dy > 0) then  result := 0                          else
  if (dy < 0) then  result := Pi                         else
                    result := cNO_ANGLE; // the 2 points are equal

  result := RadToDeg(result);
end;

• Не забудьте обработать ситуацию, когда 2 точки равны (проверьте, равен ли результат cNO_ANGLE, или измените функцию, чтобы вызвать исключение);

• Эта функция предполагает, что вы находитесь на плоской поверхности. С небольшими расстояниями, о которых вы упомянули, это все хорошо, но если вы собираетесь рассчитывать курс между городами по всему миру, вы можете захотеть взглянуть на что-то, что примет форму земли;

• Лучше предоставить этой функции координаты, которые уже сопоставлены с плоской поверхностью. Вы можете указать широту WGS84 непосредственно в Y (и продлить в X), чтобы получить приблизительное приближение.

Answer 2

Вот решение C #. Проверено на углы 0, 45, 90, 135, 180, 225, 270 и 315.

Редактировать Я заменил свое предыдущее уродливое решение переводом Wouter на C #:

public double GetAzimuth(LatLng destination)
{
    var longitudinalDifference = destination.Lng - this.Lng;
    var latitudinalDifference = destination.Lat - this.Lat;
    var azimuth = (Math.PI * .5d) - Math.Atan(latitudinalDifference / longitudinalDifference);
    if (longitudinalDifference > 0) return azimuth;
    else if (longitudinalDifference < 0) return azimuth + Math.PI;
    else if (latitudinalDifference < 0) return Math.PI;
    return 0d;
}

public double GetDegreesAzimuth(LatLng destination)
{
    return RadiansToDegreesConversionFactor * GetAzimuth(destination);
}

Answer 3

Я нашел эту ссылку

http://williams.best.vwh.net/avform.htm

дано в ответе

широта / долгота + расстояние + курс -> широта / долгота

Это выглядит многообещающе, особенно в приближении плоской земли, приведенном ближе к концу.

Answer 4

Это будет работать только для небольших различий. В противном случае вы не можете просто "latitudeDifference / long ПродольныйDifference".

Answer 5

Кто-нибудь проверял это? Не возвращает правильных ответов

Эта функция предполагает, что вы находитесь на плоской поверхности. С небольшими расстояниями, о которых вы упомянули, это все хорошо, но если вы собираетесь рассчитывать курс между городами по всему миру, вы можете захотеть взглянуть на что-то, что примет форму земли; *

Твоя плоская земля не имеет к этому никакого отношения. Ошибка, как вы это называете, заключается в том, что вы вычисляете начальный азимут из точки. Если вы не направляетесь прямо к полюсу, ваше отношение к полюсу будет меняться с расстоянием. Несмотря на это, вышеприведенная программа не возвращает правильные результаты.

Answer 6

Я бы рекомендовал ввести поправочный коэффициент на основе долготы. Я однажды реализовал симуляционную процедуру, чтобы вернуть все геокодированные записи в пределах x миль от конкретного места, и столкнулся с симуляцией. К сожалению, у меня больше нет кода, и я не могу вспомнить, как я добрался до номера коррекции, но вы на правильном пути.