2009-03-13 24 views
7

J'ai deux coordonnées WGS84, latitude et longitude en degrés. Ces points sont plutôt proches les uns des autres, par ex. seulement un mètre d'intervalle.Comment calculer l'azimut (angle au nord) entre deux coordonnées WGS84

Existe-t-il un moyen facile de calculer l'azimut de la ligne entre ces points, c'est-à-dire l'angle vers le nord?

L'approche naïve serait de supposer un système de coordonnées cartésiennes (parce que ces points sont tellement rapprochés) et il suffit d'utiliser

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

a = azimut L1, L2 = longitude B1, B2 = latitude

l'erreur est plus grande que les coordonnées se déplacer loin de l'équateur, car il la distance entre deux degrés longitudinaux deviennent de plus en plus petits que celui entre deux degrés latitudinaux (qui reste onstant).

J'ai trouvé des formules assez complexes que je ne veux pas vraiment implémenter car elles semblent trop compliquées pour des points aussi proches les uns des autres et je n'ai pas besoin de très haute précision (deux décimales suffisent, probablement une bien soit parce qu'il y a d'autres facteurs qui réduisent la précision de toute façon, comme celui que le GPS retourne).

Peut-être que pourrait simplement déterminer un facteur de correction longitudinale approximative en fonction de la latitude et utiliser somthing comme ceci:

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

où f est le facteur de correction

Tout conseils?

(Je ne veux pas utiliser des bibliothèques pour cela, surtout pas ceux qui nécessitent pas MPLed Delphi licences d'exécution. Source serait grande.)

+1

Pour ce que ça vaut, le terme que vous recherchez est "heading". – hobbs

Répondre

10

Les formules auxquelles vous faites référence dans le texte doivent calculer la distance du grand cercle entre 2 points.Voilà comment je calcule l'angle entre les points:

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; 
  • Souvenez-vous de gérer la situation où 2 points sont égaux (vérifier si le résultat est égal à cNO_ANGLE, ou modifier la fonction de lancer une exception);

  • Cette fonction suppose que vous êtes sur une surface plane. Avec les petites distances que vous avez mentionnées, tout va bien, mais si vous calculez le cap entre les villes du monde entier, vous pourriez vouloir regarder quelque chose qui prend la forme de la Terre en nombre;

  • Il est préférable de fournir cette fonction avec des coordonnées qui sont déjà mappées sur une surface plane. Vous pouvez utiliser WGS84 Latitude directement dans Y (et lon dans X) pour obtenir une approximation approximative.

+0

Excellente solution! Je l'ai traduit en C# pour remplacer ma réponse précédente. –

+1

N'est-ce pas juste atan2()? –

+1

@Martin Beckett: Oui, atan2 le fait aussi avec 'Atan2 (dy; dx)', mais bien sûr, un dx négatif renverra une valeur négative qui devra être ajoutée à 360. Et les 0 doivent être traités séparément. – MPelletier

0

je recommande la mise en œuvre d'un facteur de correction en fonction de la longitude. J'ai implémenté une routine simulaire une fois pour renvoyer tous les enregistrements géocodés dans des x miles d'un endroit spécifique et rencontré des problèmes similaires. Malheureusement, je n'ai plus le code et je n'arrive pas à me souvenir comment je suis arrivé au numéro de correction, mais vous êtes sur la bonne voie.

5

Voici la solution C#. Testé pour les angles 0, 45, 90, 135, 180, 225, 270 et 315.

Modifier J'ai remplacé ma précédente solution laide, par la traduction C# de la solution de 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); 
} 
+1

N'est-il pas préférable d'utiliser Math.Atan2 plutôt que Math.Atan pour éviter la division, et en particulier éviter la division par zéro. L'avantage supplémentaire est qu'il donne une réponse avec une gamme de -π .. + π, évitant ainsi les corrections que vous faites à la fin. –

2

Cela ne peut fonctionner que pour les petites différences. Sinon vous ne pouvez pas simplement "latitudinalDifference/longitudinalDifference".

Questions connexes