Donne la distance (en m) entre deux points de coordonnées géographiques connues
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/
Bon amusement.
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
- Donne la distance (en m) entre deux points de coordonnées géographiques connues
- Comment calculer la distance entre deux points - Forum Visual Basic 6
- Google map distance entre deux points - Forum PHP
- Calculer un azimut entre deux points - Forum PHP
- Ephéméride entre deux dates ✓ - Forum PHP
- Calcul de distance entre 2 points géodésiques - Forum Visual Basic
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.