Lambert IIe -> Latitude, Longitude Conversion

Just starting out? Need help? Post your questions and find answers here.
Post Reply
User avatar
GG
Enthusiast
Enthusiast
Posts: 266
Joined: Tue Jul 26, 2005 12:02 pm
Location: Lieusaint (77), France

Lambert IIe -> Latitude, Longitude Conversion

Post by GG »

Hi All,

Has anyone ever coded a function which converts geographic coordinates from Lambert II extended (for example x=1910245,23354 y=2410356,98523) to Latitude and Longitude coordinates (Lat=N 45.95255, Long= E 2.17552) ?
I've heard about WDGPS.dll, but i don't know the name of the function to call (only know about the reverse function LatLon_WGS84_2_LT2e).

Thanks in advance !
Purebasic 6.12 64 bits - Windows 11 Pro 64 bits 23H2
Top
User avatar
idle
Always Here
Always Here
Posts: 6029
Joined: Fri Sep 21, 2007 5:52 am
Location: New Zealand

Post by idle »

this might be of help http://forum.worldwindcentral.com/archi ... -2743.html
Top
User avatar
GG
Enthusiast
Enthusiast
Posts: 266
Joined: Tue Jul 26, 2005 12:02 pm
Location: Lieusaint (77), France

Post by GG »

Thanks Idle ,

But this example uses some weird external functions I won't be able to use. :oops:
Purebasic 6.12 64 bits - Windows 11 Pro 64 bits 23H2
Top
pif
New User
New User
Posts: 6
Joined: Sat Jun 21, 2008 1:06 pm

Post by pif »

Hello

there is a VB6 code on this page:

http://www.forumsig.org/showthread.php?t=7418&page=3

Cuq's post 21/03/2008, 16h23

@+
Top
User avatar
GG
Enthusiast
Enthusiast
Posts: 266
Joined: Tue Jul 26, 2005 12:02 pm
Location: Lieusaint (77), France

Post by GG »

Absolutely, I'm working on this code for a few weeks...

Code: Select all

; Lambert2WGS84
;
;
; |---------------------------------------------------------------------------------------------------------------|
; | 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  |
; |---------------------------------------------------------------------------------------------------------------|
Procedure Lambert2WGS84(X.d,Y.d, Type_Lambert)
; longitude.d
; latitude.d
 Longi.d
 L.d
 phi.d
 phi0.d
 phii.d
 phiprec.d
 R.d
 g.d
 VarN.d
 X_cart.d
 Y_cart.d
 Z_cart.d
 XWGS84.d
 YWGS84.d
 ZWGS84.d
 p.d
 phi840.d
 phi84prec.d
 phi84i.d
 phi84.d
 l84.d
 l840.d
 Const_a.d
 Const_e.d

 Decalage_X_cart.d
 Decalage_Y_cart.d
 Decalage_Z_cart.d

 l0.d


;--  Quelques constantes ...  -->
 n.d
 C.d           ;--   En mètres   --
 Xs.d          ;--   En mètres   --
 Ys.d          ;--   En mètres   --
   
   

eps.d = 0.000000000001        ;--   précision   --
h = 100                   ;--   En mètres   --

l0 = 0                          ;--  correspond à la longitude en radian de Paris (2°20'14.025" E) par rapport à Greenwich   --
pi.d=3.1415926535897
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   --
Select Type_Lambert

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
EndSelect

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

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

;phi0 = 2 * ATan((Pow(e,L)) - (Pi / 2))
phi0 = 2 * ATan((Pow(e,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#)
phii = 2 * ATan(((Pow(((1 + Const_e * Sin(phiprec)) / (1 - Const_e * Sin(phiprec))) , (Const_e / 2)) * Pow(e,L)))) - (Pi / 2)
;phii = 2 * ATan(((((1 + Const_e * Sin(phiprec)) / Pow((1 - Const_e * Sin(phiprec)) , (Const_e / 2)) * Pow(e,L))) - (Pi / 2)

While Not (Abs(phii - phiprec) < eps)
        phiprec = phii
        ;phii = 2 * ATan((((((1 + Const_e * Sin(phiprec)) / Pow(((1 - Const_e * Sin(phiprec))) , (Const_e / 2)))) * Pow(e,L)))) - (Pi / 2)
				phii = 2 * ATan(((Pow(((1 + Const_e * Sin(phiprec)) / (1 - Const_e * Sin(phiprec))) , (Const_e / 2)) * Pow(e,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)
VarN = Const_a / Pow(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 + ATan(YWGS84 / XWGS84)

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

phi84prec = phi840


phi84i = ATan((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 = ATan((ZWGS84 / p) / (1 - ((Const_a * Const_e * Const_e * Cos(phi84prec)) / (p * Sqr(1 - ((Const_e * Const_e) * (Sin(phi84prec) * Sin(phi84prec))))))))
        phi84i = ATan((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

LTversWGS84X.d = l84 * 180 / Pi
LTversWGS84Y.d = phi84 * 180 / Pi
OpenConsole()
PrintN("X="+StrD(LTversWGS84X)+" Y="+StrD(LTversWGS84Y))

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


EndProcedure

Lambert2WGS84(616388,2404200,5)
Maybe I'm wrong in converting with parenthesis...
EDIT :

Error found, it's OK now.

Thanks ! :)
Purebasic 6.12 64 bits - Windows 11 Pro 64 bits 23H2
Top
Post Reply

Return to “Coding Questions”