Donne la distance (en m) entre deux points de coordonnées géographiques connues

Soyez le premier à donner votre avis sur cette source.

Snippet vu 34 950 fois - Téléchargée 26 fois

Contenu du snippet

en faisant usage de l'ellispoide WGS-84

inspiré d'un script javascript
http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
http://www.movable-type.co.uk/

Source / Exemple :


Attribute VB_Name = "modDistance_2006_05_29_B"
Option Explicit

Public Const PI      As Double = 3.14159265358979     
Public Const PI_DEMI As Double = 1.5707963267949
Public Const EPSILON As Double = 0.00001 ' 1E-05       

Type Pt_GEOGRAPHIQUE
    Lon      As Double ' radian
    Lat      As Double ' radian
    LonDegDD As Double ' decimal degree
    LatDegDD As Double ' decimal degree
End Type

Function dtR(deg As Double) As Double
    dtR = deg * (PI / 180)
End Function

Function arctan2(y As Double, x As Double) As Double
    If Abs(x) < 0.000000000000001 Then   '1E-15
        ' donc x = 0
        arctan2 = IIf(y < 0, -PI_DEMI, PI_DEMI)
    Else
        arctan2 = Arctan(y / x) + IIf(x > 0#, 0#, IIf(y < 0#, -PI, PI))
    End If
End Function
'---------------------------------------------------------------------------------------
' Procedure : distVincenty
' DateTime  : 29-05-2006  20:26
' Author    : CADRATURE
' Purpose   : Donne la distance (en m) entre deux points de coordonnées géographiques connues
'           :        en faisant usage de l'ellispoide WGS-84
'           :        inspiré d'un script javascript
'           :        http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
'           :        http://www.movable-type.co.uk/
'---------------------------------------------------------------------------------------
Public Function distVincenty(p1 As Pt_GEOGRAPHIQUE, p2 As Pt_GEOGRAPHIQUE) As Double

    On Error GoTo Error_distVincenty

    Dim sMot As String
    Dim sTemp As String
    Dim sFile As String
    Dim cP As Integer
    Dim i As Integer
    Dim k As Integer
    Dim fA As Byte
    Dim fB As Byte
    Dim aT As Double
    Dim bT As Double
    Dim f As Double
    Dim A As Double
    Dim B As Double
    Dim C As Double
    Dim cos2SigmaM As Double
    Dim cosSigma As Double
    Dim cosSqAlpha As Double
    Dim deg As Double
    Dim deltaSigma As Double
    Dim iterLimit As Double
    Dim L As Double
    Dim lambda As Double
    Dim s As Double
    Dim sigma As Double
    Dim sinAlpha As Double
    Dim sinLambda As Double
    Dim sinSigma As Double
    Dim sinU2 As Double
    Dim U1 As Double
    Dim U2 As Double
    Dim uSq As Double
    Dim sinU1 As Double
    Dim cosU1 As Double
    Dim cosU2 As Double
    Dim lambdaP As Double
    Dim cosLambda As Double
    Dim rA As Double
    Dim rB As Double
    Dim rC As Double
    Dim rD As Double
    Dim rE As Double
    Dim dTempA As Double
    Dim dTempB As Double
    Dim dTempC As Double
    Dim oldEcart As Double

    aT = 6378137
    bT = 6356752.3142
    f = 1 / 298.257223563 ' WGS-84 ellipsoid
    L = (1000000 * p2.Lon - 1000000 * p1.Lon) / 1000000
    U1 = Arctan((1 - f) * Tan(p1.Lat))
    U2 = Arctan((1 - f) * Tan(p2.Lat))
    sinU1 = sin(U1)
    cosU1 = cos(U1)
    sinU2 = sin(U2)
    cosU2 = cos(U2)

    lambda = L
    lambdaP = 2 * PI
    iterLimit = 20

    While (Abs(lambda - lambdaP) > EPSILON And iterLimit > 0)
        oldEcart = Abs(lambda - lambdaP)

        DoEvents
        iterLimit = iterLimit - 1
        sinLambda = sin(lambda)
        cosLambda = cos(lambda)
        sinSigma = Sqr((cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda))
        If (sinSigma = 0) Then
            distVincenty = 0 ' co-incident points
            Exit Function
        End If
        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
        sigma = arctan2(sinSigma, cosSigma)
        DoEvents
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha = 1 - sinAlpha * sinAlpha

        DoEvents

        If (cosSqAlpha = 0) Then
            distVincenty = Abs(aT * L) ' two points on equator
        End If

        cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha

        DoEvents
        C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
        lambdaP = lambda

        dTempA = (-1 + 2 * cos2SigmaM * cos2SigmaM)
        dTempB = (cos2SigmaM + (C * cosSigma) * dTempA)
        dTempC = (1 - C) * (f * sinAlpha)
        lambda = L + dTempC * (sigma + C * sinSigma * dTempB)

        DoEvents

        rA = (1 - C) * (f * sinAlpha)
        rB = (sigma + C * sinSigma)
        rC = (cos2SigmaM + (C * cosSigma) * (-1 + 2 * cos2SigmaM * cos2SigmaM))
        lambdaP = L + rA * rB * rC

        Debug.Print iterLimit & " " & Abs(lambda - lambdaP) * 100000#
        DoEvents
    Wend
    If (iterLimit = 0) And ((Abs(lambda - lambdaP) - oldEcart) / oldEcart) > 0.01 Then
        distVincenty = -1  ' formula failed to converge
        Debug.Assert False
        Exit Function
    End If
    uSq = cosSqAlpha * (aT * aT - bT * bT) / (bT * bT)

    'A-------------
    rA = (-768 + uSq * (320 - 175 * uSq))
    rB = (4096 + uSq * rA)
    A = 1 + uSq / 16384 * rB

    'B-------------
    B = 0
    rA = (-128 + uSq * (74 - 47 * uSq))
    rB = (256 + uSq * rA)
    B = uSq / 1024 * rB

    'deltaSigma
    rA = (-3 + 4 * cos2SigmaM * cos2SigmaM)
    rB = (-3 + 4 * sinSigma * sinSigma)
    rC = ((6 * cos2SigmaM) * rA)
    rD = (-1 + 2 * cos2SigmaM * cos2SigmaM)
    rE = (4 * (cosSigma * rD - B) / rC)
    deltaSigma = B * sinSigma * (cos2SigmaM + B / rE)

    '----------
    s = bT * A * (sigma - deltaSigma)

    distVincenty = dblRoundOff(s, 3) ' round to 1mm precision

Sortir_distVincenty:
    On Error GoTo 0
    Exit Function

Error_distVincenty:

End Function

Sub Main()
    Dim dDistance As Double
    Dim p2 As Pt_GEOGRAPHIQUE
    Dim p1 As Pt_GEOGRAPHIQUE
    
    p1.LatDegDD = 52.874
    p1.LonDegDD = 4.389
    
    p2.LatDegDD = 45.001
    p2.LonDegDD = 15.716
    
    
    p1.Lat = dtR(p1.LatDegDD)
    p1.Lon = dtR(p1.LonDegDD)
    
    p2.Lat = dtR(p2.LatDegDD)
    p2.Lon = dtR(p2.LonDegDD)
    
    
    dDistance = distVincenty(p1, p2)
    Debug.Assert False
End Sub

Conclusion :


Bon amusement.

A voir également

Ajouter un commentaire

Commentaires

Messages postés
17
Date d'inscription
jeudi 21 août 2008
Statut
Membre
Dernière intervention
10 septembre 2011

STP, j'ai un peu de mal pour réaliser la forme, que dois je mettre comme Textbox et script dans la fonction clic ?
Messages postés
15
Date d'inscription
mardi 16 mars 2004
Statut
Membre
Dernière intervention
16 avril 2008

On vient de me demander par mail comment calculer la distance loxo entre deux points sur une sphère. C'est assez simple :

A(xa,ya), B(xb,yb)

CapLoxo = Atn((xa - xb) / Log(Tan(ya / 2 + pi / 4) / Tan(yb / 2 + pi / 4)))

La distance Loxo ainsi :
DstLoxo = RayonSphere * Abs((ya - yb) / Cos(CapLoxo))
Messages postés
13
Date d'inscription
mardi 14 novembre 2000
Statut
Membre
Dernière intervention
18 janvier 2008

Merci pour cette précision !
Messages postés
15
Date d'inscription
mardi 16 mars 2004
Statut
Membre
Dernière intervention
16 avril 2008

Je vais me permettre de répondre à la place de vbsorcier.
1.609km est la valeur d'un "miles" anglais. En revanche en navigation on utilise le "mile nautique" qui vaut 1.852 km.
Messages postés
13
Date d'inscription
mardi 14 novembre 2000
Statut
Membre
Dernière intervention
18 janvier 2008

@vbsorcier : hum je me pose une question : un mile = 1.609km... pourquoi mets tu ceci " 'On multiplie par 1.853 pour convertir en Km car ce sont des miles" ???!!! o_O
Afficher les 13 commentaires

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.