Quellcode

Option Explicit

Global System_Coord As Integer

'---------------------------------------------------------------------------
' Point en 2D
'---------------------------------------------------------------------------
Public Type Point2
    X As Double
    Y As Double
End Type


' Convertion de données en lambert II vers position GPS
' X et Y en mètres
' Option Système de coordonnées
' |---------------------------------------------------------------------------------------------------------------|
' | Const | 1 'Lambert I | 2 'Lambert II | 3 'Lambert III | 4 'Lambert IV | 5 'Lambert II Etendue | 6 'Lambert 93 |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    n  | 0.7604059656 |  0.7289686274 |   0.6959127966 | 0.6712679322  |    0.7289686274       |  0.7256077650 |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    c  | 11603796.98  |  11745793.39  |   11947992.52  | 12136281.99   |    11745793.39        |  11754255.426 |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    Xs |   600000.0   |    600000.0   |   600000.0     |      234.358  |    600000.0           |     700000.0  |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    Ys | 5657616.674  |  6199695.768  |   6791905.085  |  7239161.542  |    8199695.768        | 12655612.050  |
' |---------------------------------------------------------------------------------------------------------------|
'
'Function LTversWGS84(ByVal X As Double, ByVal Y As Double, Optional SystemCoord As Integer = 5) As Point2
Function LTversWGS84(ByVal X As Double, ByVal Y As Double, Optional SystemCoord As Integer = 5) As String

Dim longitude As Double
Dim latitude As Double
Dim Longi As Double
Dim L As Double
Dim phi As Double
Dim phi0 As Double
Dim phii As Double
Dim phiprec As Double
Dim R As Double
Dim g As Double
Dim VarN As Double
Dim X_cart As Double
Dim Y_cart As Double
Dim Z_cart As Double
Dim XWGS84 As Double
Dim YWGS84 As Double
Dim ZWGS84 As Double
Dim p As Double
Dim phi840 As Double
Dim phi84prec As Double
Dim phi84i As Double
Dim phi84 As Double
Dim l84 As Double
Dim l840 As Double
Dim Const_a As Double
Dim Const_e As Double

Dim Decalage_X_cart As Double
Dim Decalage_Y_cart As Double
Dim Decalage_Z_cart As Double

Dim l0 As Double

'Debug.Print "Conversion Lambert 2 - NTF géographique"

'--  Quelques constantes ...  -->
Dim n As Double
Dim C As Double           '--   En mètres   --
Dim Xs As Double          '--   En mètres   --
Dim Ys As Double          '--   En mètres   --
    
Dim Pi As Double

Pi = 3.14159265

Const eps = 0.000000000001        '--   précision   --
Const h = 100                   '--   En mètres   --

l0 = 0                          '--  correspond à la longitude en radian de Paris (2°20'14.025" E) par rapport à Greenwich   --
l840 = 2.337229167 / 180 * Pi '0.04079234433            '--   0.04079234433 pour passer dans un référentiel par rapport au méridien   --
                                '--   de Greenwich, sinon mettre 0   --

Const_a = 6378249.2
Const_e = 0.08248325676         '--  Const_e du NTF (on le change après pour passer en WGS)   --
                                '-- (première excentricité de l’ellipsoïde Clarke 1880 français). ---
                                
 Decalage_X_cart = -168         '--   En mètres   --
 Decalage_Y_cart = -60          '--   En mètres   --
 Decalage_Z_cart = 320          '--   En mètres   --


'Definition selon systèmes de coordonnées Lambert
Select Case SystemCoord

Case 1 'Lambert I

     n = 0.7604059656
     C = 11603796.98
     Xs = 600000
     Ys = 5657616.674
     
Case 2 'Lambert II

     n = 0.7289686274
     C = 11745793.39
     Xs = 600000
     Ys = 6199695.768

Case 3 'Lambert III
     n = 0.6959127966
     C = 11947992.52
     Xs = 600000
     Ys = 6791905.085
     
Case 4 'Lambert IV
     n = 0.6712679322
     C = 12136281.99
     Xs = 234.358
     Ys = 7239161.542
     
Case 5 'Lambert II Etendue
     n = 0.7289686274
     C = 11745793.39
     Xs = 600000
     Ys = 8199695.768
     
Case 6 'Lambert 93
     
     
        n = 0.725607765
        C = 11754255.426
        Xs = 700000
        Ys = 12655612.05


        'Pour lambert 93 Longitude origine 3° Est Greenwich
        l0 = (3 / 180 * Pi) ' - 0.04079234433198
        l840 = 0
        
        
        Decalage_X_cart = 0
        Decalage_Y_cart = 0
        Decalage_Z_cart = 0

End Select

R = Sqr(((X - Xs) * (X - Xs)) + ((Y - Ys) * (Y - Ys)))
g = Atn((X - Xs) / (Ys - Y))

Longi = l0 + (g / n)
L = -(1 / n) * Log(Abs(R / C))


phi0 = 2 * Atn(Exp(L)) - (Pi / 2#)
phiprec = phi0
phii = 2 * Atn(((((1 + Const_e * Sin(phiprec)) / (1 - Const_e * Sin(phiprec))) ^ (Const_e / 2#)) * Exp(L))) - (Pi / 2#)

While Not (Abs(phii - phiprec) < eps)
        phiprec = phii
        phii = 2 * Atn((((((1 + Const_e * Sin(phiprec)) / (1 - Const_e * Sin(phiprec))) ^ (Const_e / 2#)) * Exp(L)))) - (Pi / 2#)

Wend


phi = phii

'Debug.Print "Lambda = " & Longi & " rad = " & Longi * 200 / Pi & "gr"
'Debug.Print "Phi    = " & phi & " rad = " & phi * 200 / Pi & "gr"



'-- Conversion NTF géographique - NTF cartésien : ALG0009   --
'Debug.Print "Conversion NTF géographique - NTF cartésien"


VarN = Const_a / ((1 - (Const_e * Const_e) * (Sin(phi) * Sin(phi))) ^ 0.5)
X_cart = (VarN + h) * Cos(phi) * Cos(Longi)
Y_cart = (VarN + h) * Cos(phi) * Sin(Longi)
Z_cart = ((VarN * (1 - (Const_e * Const_e))) + h) * Sin(phi)

'Debug.Print "X cartésien NTF = " & X_cart
'Debug.Print "Y cartésien NTF = " & Y_cart
'Debug.Print "Z cartésien NTF = " & Z_cart

'-- Conversion NTF cartésien - WGS84 cartésien : ALG0013   --
'Debug.Print "Conversion NTF cartésien - WGS84 cartésien"

'--   Il s'agit d'une simple translation   --
XWGS84 = X_cart + Decalage_X_cart
YWGS84 = Y_cart + Decalage_Y_cart
ZWGS84 = Z_cart + Decalage_Z_cart


'Debug.Print "X cartésien WGS84 = " & XWGS84
'Debug.Print "Y cartésien WGS84 = " & YWGS84
'Debug.Print "Z cartésien WGS84 = " & ZWGS84


'-- Conversion WGS84 cartésien - WGS84 géographique : ALG0012   --
'Debug.Print "Conversion WGS84 cartésien - WGS84 géographique"

Const_e = 0.08181919106              '--   On change Const_e pour le mettre dans le système WGS84 au lieu de NTF   --
Const_a = 6378137


p = Sqr((XWGS84 * XWGS84) + (YWGS84 * YWGS84))

l84 = l840 + Atn(YWGS84 / XWGS84)

phi840 = Atn(ZWGS84 / (p * (1 - ((Const_a * Const_e * Const_e)) / Sqr((XWGS84 * XWGS84) + (YWGS84 * YWGS84) + (ZWGS84 * ZWGS84)))))

phi84prec = phi840


phi84i = Atn((ZWGS84 / p) / (1 - ((Const_a * Const_e * Const_e * Cos(phi84prec)) / (p * Sqr(1 - Const_e * Const_e * (Sin(phi84prec) * Sin(phi84prec)))))))

While Not (Abs(phi84i - phi84prec) < eps)
        phi84prec = phi84i
        phi84i = Atn((ZWGS84 / p) / (1 - ((Const_a * Const_e * Const_e * Cos(phi84prec)) / (p * Sqr(1 - ((Const_e * Const_e) * (Sin(phi84prec) * Sin(phi84prec))))))))
Wend

phi84 = phi84i

' LTversWGS84.X = l84 * 180 / Pi
' LTversWGS84.Y = phi84 * 180 / Pi
LTversWGS84 = "N" & Str(phi84 * 180 / Pi) & " E" & Str(l84 * 180 / Pi)

'Debug.Print "latitude WGS84  = " & l84 & " rad = " & l84 * 180 / Pi & " deg"
'Debug.Print "longitude WGS84 = " & phi84 & " rad = " & phi84 * 180 / Pi & " deg"
End Function

'conversion longitude latitude en texte
Function Texte_Position(longitude As Double, latitude As Double) As String
Dim TexteTempo As String
Dim Temp As Double
Dim degres As Integer
Dim minutes As Integer
Dim secondes As Integer
Dim Signe As String

             Temp = Abs(longitude)
             degres = Int(Temp)
             minutes = Int((Temp - degres) * 60)
             secondes = Round((Temp - degres - minutes / 60) * 3600)
                
            If (Temp > 0) Then
                Signe = " E "
            Else
                Signe = " O "
            End If
                
             TexteTempo = degres & "° " & minutes & "' " & secondes & "'' " & Signe
                       
            Temp = Abs(latitude)
            degres = Int(Temp)
            minutes = Int((Temp - degres) * 60)
            secondes = Round((Temp - degres - minutes / 60) * 3600)
            
            If (Temp > 0) Then
                Signe = " N "
            Else
                Signe = " S "
            End If
                
            TexteTempo = TexteTempo & vbCrLf & degres & "° " & minutes & "' " & secondes & "'' " & Signe
            
            Texte_Position = TexteTempo

End Function