Come calcolare Azimuth (angolo a nord) tra due coordinate WGS84

Question

Ho due coordinate WGS84, latitudine e longitudine in gradi. Questi punti sono piuttosto vicini tra loro, ad es. solo un metro di distanza.

C'è un modo semplice per calcolare l'azimut della linea tra questi punti, cioè l'angolo verso nord?

L'approccio più semplice sarebbe quella di assumere un sistema di coordinate cartesiane (perché questi punti sono così vicini) e solo uso

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

a = azimut L1, L2 = longitudine B1, B2 = latitudine

l'errore sarà più grande come le coordinate allontanano dall'equatore perché la distanza tra due gradi longitudinali diventano sempre più piccoli di quello tra due gradi latitudinali (che rimane c onstant).

Ho trovato alcune formule abbastanza complesse che non voglio veramente implementare perché sembrano essere eccessive per punti così vicini e non ho bisogno di una precisione molto alta (sono sufficienti due decimali, uno è probabilmente bene perché ci sono altri fattori che riducono comunque la precisione, come quello che il GPS restituisce).

Forse ho potuto solo determinare un fattore di correzione longitudinale approssimativa seconda della latitudine e utilizzando somthing simili:

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

dove f è il fattore di correzione

Eventuali suggerimenti?

(io non voglio usare le librerie per questo, soprattutto non quelli che richiedono le licenze di runtime. Qualsiasi MPLed Delphi Fonte sarebbe grande.)

Answer 1

Le formule a cui si fa riferimento nel testo sono per calcolare la distanza del cerchio grande tra 2 punti.Ecco come a calcolare l'angolo tra i punti:

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; 

• Ricordati di gestire la situazione in cui 2 punti sono uguali (controllare se il risultato è uguale a cNO_ANGLE, o modificare la funzione per generare un'eccezione);

• Questa funzione presuppone che tu sia su una superficie piana. Con le piccole distanze che hai menzionato, tutto va bene, ma se stai calcolando la rotta tra le città di tutto il mondo potresti voler esaminare qualcosa che prende in considerazione la forma della terra;

• È meglio fornire questa funzione con coordinate già mappate su una superficie piana. Potresti alimentare WGS84 Latitude direttamente in Y (e lon in X) per ottenere un'approssimazione approssimativa.

Answer 2

mi sento di raccomandare l'attuazione di un fattore di correzione in base alla longitudine. Ho implementato una routine simulata una volta per restituire tutti i record geocodificati entro x miglia da un punto specifico e ho riscontrato problemi similari. Sfortunatamente non ho più il codice e non riesco a ricordare come ho raggiunto il numero di correzione ma sei sulla strada giusta.

Answer 3

Ho trovato questo link

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

dato nella risposta alla

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

Questo sembra essere molto promettente, in particolare la terra piatta approssimazione data verso la fine.

Answer 4

Ecco la soluzione C#. Testato per 0, 45, 90, 135, 180, 225, 270 e 315 angoli.

Edit ho sostituito il mio precedente soluzione brutto, dalla C# traduzione della soluzione di Wouter:

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

questo dovrebbe funzionare solo per piccole differenze. Altrimenti non puoi semplicemente "latitudinalDifference/longitudinalDifference".