Jak obliczyć Azymut (kąt na północ) między dwiema współrzędnymi WGS84

Question

Mam dwie współrzędne WGS84, szerokość i długość geograficzną w stopniach. Punkty te są raczej blisko siebie, np. tylko jeden metr od siebie.

Czy istnieje prosty sposób obliczenia azymutu linii między tymi punktami, to znaczy kąta na północ?

Naiwny podejście byłoby zakładać kartezjańskim układzie współrzędnych (ponieważ te punkty są tak blisko siebie) i po prostu użyć

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

a = azymut L1, L2 = długość B1, B2 = szerokość

błąd będzie większy jako współrzędne odchodzenia od równika, ponieważ odległość pomiędzy dwa stopnie podłużne stają się coraz mniejsze niż między dwoma równaniami równoleżnikowymi (które pozostają w c onstant).

Znalazłem kilka dość skomplikowanych formuł, których naprawdę nie chcę zaimplementować, ponieważ wydają się być przesadne dla punktów, które są tak blisko siebie i nie potrzebuję bardzo dużej precyzji (dwa miejsca po przecinku wystarczą, jeden prawdopodobnie dobrze, ponieważ istnieją inne czynniki, które i tak zmniejszają precyzję, jak ta, którą zwraca GPS).

Może po prostu określić przybliżoną współczynnik korekcyjny podłużny zależnie od szerokości geograficznej i użyć somthing tak:

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

gdzie f jest współczynnikiem korekcji

Wszelkie wskazówki?

(nie chcę używać żadnych bibliotek do tego, zwłaszcza te, które nie wymagają licencji uruchomieniowe. Wszelkie MPLed Delphi Źródło byłoby świetnie.)

Answer 1

Formuły, do których odwołujesz się w tekście, służą do obliczenia odległości między dwoma punktami.Oto w jaki sposób obliczyć kąt między punktami:

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; 

• Pamiętaj, aby poradzić sobie z sytuacją, gdzie 2 punkty są równe (sprawdzić, czy wynik jest równy cNO_ANGLE lub zmodyfikować funkcję wyjątek);

• Ta funkcja zakłada, że ​​użytkownik znajduje się na płaskiej powierzchni. Przy niewielkich odległościach, o których wspomniałeś, wszystko jest w porządku, ale jeśli zamierzasz obliczyć kurs między miastami na całym świecie, możesz chcieć spojrzeć na coś, co przybiera kształt ziemi;

• Najlepiej podać tę funkcję z współrzędnymi, które są już zmapowane na płaskiej powierzchni. Możesz podać WGS84 Latitude bezpośrednio na Y (i na X), aby uzyskać przybliżone przybliżenie.

Answer 2

Polecam realizacji współczynnika korekcyjnego na podstawie długości geograficznej. Zaimplementowałem jednorazową rutynę, aby zwrócić wszystkie geokodowane rekordy w promieniu x mil od określonego miejsca i pojawiły się w nim różne problemy. Niestety nie mam już kodu i nie mogę sobie przypomnieć, jak dotarłem do numeru korekty, ale jesteś na dobrej drodze.

Answer 3

znalazłem ten link

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

podane w odpowiedzi na

Lat/Lon + Distance + Heading --> Lat/Lon

To wygląda obiecująco, zwłaszcza Flat Earth przybliżenie danego blisko końca.

Answer 4

Oto rozwiązanie C#. Testowany pod kątem 0, 45, 90, 135, 180, 225, 270 i 315.

Edit Wymieniłem mój poprzedni brzydki rozwiązanie, przez # tłumaczeniu C roztworu Wouter za:

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 5

ta będzie działać tylko do niewielkich różnic. W przeciwnym razie nie można po prostu "równania południkowego/podłużnej różnicy".