Convertion Longitude/Latitude en Pixel [résolu]

[MetalOS]

Salut Kernadec et merci de ton aide.

J'utilisais déjà tes formules sur mon ancien programme et celle que tu vient de poster sont identique et du coup ca ne marche pas avec le reste de mon projet car la projection n'est plus bonne.

Il faut vraiment que j'utilise la formule de mon code d'exemple ainsi que la carte sans la modifier. C'est juste la conversion de pixels en coordonnées décimal que je n'arrive pas à trouver.

[GallyHC]

Bonjour,

Vu la taille de la carte, tu n'auras jamais le résultat exact. A pars si tu peux faire un zoom pour pouvoir redéfinir l’échelle avec plus de « détail ».

Pourquoi : la raison est simple tu as des coordonnée en degré (ou même degré minute ou autres). Ta conversion réel, ne donne pas un point a X, Y pixel mais dans le calcul a X.xxx et y.yyy et c'est donc c'est .xxx et .yyy qui poseront toujours problème.

En gros il faudrait une résolution très grande de ta carte pour avoir un résultat plus ou moins acceptable.

La conversion d'un point en degré/pixel ne pose pas trop de problème car converti cela donne un résultat en pixel avec + ou – un bon pixel.

Pour faire simple voila une donné en degré minute seconde et surtout ce que cela représente

Code : Tout sélectionner

	Donc  1 degré représente à peu près 111.11 Km à l’équateur…
	Une minute ~1.85 Km
	Une seconde ~30m

Dans un calcul en degré simple minute et seconde sont dans le virgule donc 1.85 + 0.30 sois une différence envisagable de 2.15 Km avec des calculs précis de pixel.

La question comment convertir 1 en X.xxx (réel) sans aucune précision, cela est impossible (pour ma part, avec une carte comme celle-la).

Cordialement,
GallyHC

[kernadec]

bonjour MetalOS

Je suis étonné que ce code ne te convienne pas car il converti des coordonnées géographiques
en une position pixel et inversement..
je ne vois pas vraiment la raison pour laquelle il te faudrait une projection à plat pour réaliser la même chose.

cette méthode permet de gérer 5 degrés de part et d'autre d'une latitude référence idem pour la longitude.

ce code utilise des cartes Google Map directement, par exemple pour l’Allemagne on pourrait l'adapter
en prenant une latitude moyenne de 51 degrés et un méridien de 10 degrés pour la longitude
on détermine ce point de croisement sur la carte et on centre le point de référence latitude longitude dessus.

ensuite on met place quelques coordonnées villes allemandes et on ajuste la position avec l’échelle dans le programme,
et pour optimiser On peut aussi donner de l'inclinaison à ces points au cas ou.

Cordialement

[MetalOS]

Si le code me convient bien mais tu à été obliger de modifier la carte d'origine pour que sa fonctionne. Hors j'ai vraiment besoin de garder les détails d'origine de la carte pour la suite de mon projet.

Voici l'image avec ton code:



Et voici l'image d'origine:



Sur ton code l'image n'est pas net et il me manque une bordure sur le tour de la carte. C'est vrai que ta formule de départ me convenait bien mais pas pour cet carte. C'est pour ca que j'ai changer de formule pour la projection. Pense tu qu'il est possible d'adapter ta formule pour conserver ma carte d'origine ?

[kernadec]

Je connais pas la carte que tu as utilisée.
Mais il a un moyen simple de l'utiliser comme tu l'entend si c'est vraiment
une projection à plat de style Mercator.

pour infos:
je te rappel que Mercator utilise des rectangles dans ses projections, ce qui signifie que le rapport
entre la latitude et la longitude est différent.
La cause, c' est qu'un planisphère, c'est un cylindre avec un rapport (périmètre/diamètre) de la sphère

tu as probablement pensé à la solution suivante pour régler ce problème en prenant deux points connus,
un point avec Lille et pour l'autre que tu descend verticalement jusqu’à la perpendiculaire de Marseille,
on auras une distance pour la verticale en pixels et une distance en latitude,
Ce rapport permettra de passer de pixel en latitude en vice versa.

il faudra faire la même chose avec Brest et Strasbourg pour le rapport pixel en longitude.

Mais cela donnera une moins bonne précision pour les calculs de distances entre les points.


Cordialement

[MetalOS]

Non la carte n'est pas une projection de Mercator mais plane surement de type lambert mais pas sûr. J'ai juste les limites de la carte qui son dans le code que j'ai mis.

[kernadec]

bonjour MetalOS
Tu as sans doute raison,
La carte semble provenir du site GeoPortail
Mode: "Lambert II étendu" pour avoir une carte la France entière et plus...
mais en réalité ces coordonnées Lambert travaille en Zones Coniques
et ce procédé prend rapidement de la déformation au delà de zones de plus de 4 degrés.
Et c'est pour cette raison que la France est coupée en zones horizontales

Lambert I de 51,00 à 48,00 degrés
Lambert II de 48,00 à 45,30 degrés
Lambert III de 45,30 à 42,30 degrés
Lambert IV de 43,00 à 41,00 degrés pour la Corse

Alors en faisant une copie d’écran d'une carte de la même zone en copie d'écran sur le site
et que je la superpose avec La "carte marron" il y a une différence.

Pour la longitude c'est bon.. mais si l'on fait coller la partie Sud
le nord de la Bretagne jusqu’à la Hollande et l'Angleterre
prennent une différence de 0,3cm à 1cm dans la verticale avec une fenêtre 800 800.

Pour travailler avec ce type de carte,
il faudrait écrire un convertisseur Lambert en pixel et l'inverse
cela ne semble pas facile du tout avec toutes ces bandes à gérer,
je sais pas faire... mais peut être qu'il existe des codes qui le font

Peut être que ce code en java va t'aider:
http://files.codes-sources.com/fichier. ... mbert.java

Cordialement

[MetalOS]

Je vais regarder tous ca. Merci encore pour ton aide Kernadec et pour la carte que tu m'a envoyé par mail

[kernadec]

justement peux tu poster cette image merci.

comme les coordonnées Lambert ont une lecture en " mètres- kilomètres"
il te faut prendre une carte du site Geoportail et le rapport pixel coordonnées et de 1.46 environ à optimiser


j'ai oublié de te préciser pour cette carte.

Le point de référence sur la carte est en partant 0,0 de la fenêtre
X= 389 pixel et égal à 2.00 degrés de longitude EST
Y= 361 pixel et égal à 47.00 degrés de latitude NORD

toutes les mesures en pixels calculées de ce point central de 389/361 par rapport à la position de la souris dans la fenêtre
si elles sont multipliées par 1.46 donne une distance en km du point dans les X idem en Y
il faudra reprendre la formule du mile nautique mode inverse pour obtenir une coordonnée géographique.

on pourrai également recalculer la valeur en coordonnée du point 0,0 de la fenêtré
pour avoir une lecture directe depuis la position de la souris

cela dit, en prenant un point plus prés du centre on réduit la marge d'erreur avec des distances moins grandes.

[réédit] Avec des tests plus poussés je n'ai obtenu que des résultats vraiment décevants..
après réflexion cette méthode ne peut pas être juste car le Lambert est une disposition en trapèze.
sur un écran X,Y, la meilleur méthode est de considérer l'écran comme un planisphère
Mercator avait surement prévu qu'on aurai besoin un jour de sa méthode pour nos écrans

Cordialement

[MetalOS]

Voici la carte.

[kernadec]

Bonjour Metalos
comme tu utilise les carte de l'IGN, il serait raisonnable d'utiliser les coordonnées Lambert

Et je t'ai trouvé un site qui donne un code en Basic IBM
Pour convertir WGS84 en Lambert 1,2,3,4 et l' inverse, reste plus qu'a le convertir en PB.

Au passage je remercie le remarquable travail de "PK1157"

je précise que sur le site dont l'URL est en début de code
se trouve un fichier d'exemple sous OpenOffice. Calc "geo.ods"
j'ai adjoint ce texte en fin de code.

Cordialement

Code : Tout sélectionner

;http://forum.openoffice.org/fr/forum/viewtopic.php?f=8&t=26883

;[Calc] Conversion de coordonnées Géo ou Cartographiques

;Messagepar PK1157 » 30 Jan 2011 23:53
;J'ai développé une collection de fonctions Calc permettant les conversions de coordonnées suivantes :
;Cartographiques Lambert I, II , II étendu, III, IV vers géographiques NTF et conversion réciproque ;
;Géographiques NTF ou Cartographiques Lambert 93 vers WGS84 et conversion réciproque ;
;WGS84 vers UTM et Réciproque.
;Cela intéresse-t-il quelqu'un ?
;- LatIsom et IsomLat sont des fonctions "de cuisine interne" de conversion pour latitude isométrique
;- LAMB_LatNTF permet de trouver la latitude géographique NTF à partir d'un couple de coordonnées Lambert I, II , II étendu, III, IV ou 93
;- LAMB_LonNTF permet de trouver la longitude géographique NTF à partir d'un couple de coordonnées Lambert I, II , II étendu, III, IV ou 93
;- NTF_LatWGS permet de trouver la latitude géographique WGS84 à partir d'un couple de coordonnées géographiques NTF
;- NTF_LonWGS permet de trouver la longitude géographique WGS84 à partir d'un couple de coordonnées géographiques NTF
;- NTF_XLamb permet de trouver l'"easting" Lambert I, II , II étendu, III, IV ou 93 à partir d'un couple de coordonnées géographiques NTF
;- NTF_YLamb permet de trouver le"northing" Lambert I, II , II étendu, III, IV ou 93 à partir d'un couple de coordonnées géographiques NTF
;- WGS_LatNTF permet de trouver la latitude géographique NTF à partir d'un couple de coordonnées géographiques WGS84
;- WGS_LonNTF permet de trouver la longitude géographique NTF à partir d'un couple de coordonnées géographiques WGS84
;- WGS_XL93 et WGS_YL93 convertissent les coordonnées WGS84 en Projection Lambert 93
;- L93_LatWgs et L93_LonWgs effectuent la conversion inverse.
;- DecDMS et DMSDec traduisent des degrés décimaux en Degrés, Minutes, Secondes et réciproquement.
;-Utm_LatWGS et Utm_LonWGS permettent de convertir les coordonnées cartographiques UTM en géographiques WGS84
;- WGS_FusZonUTM, WGS_UTMEasting et WGS_UTMNorthing effectuent les conversions réciproques.

REM*****BASIC*****
REM ==========================================
REM
REM   Remerciements et Documentation
REM
REM   Pour la plupart des algorithmes implementes, un grand merci à l'I.G.N.
REM   (Institut Geographique National - France), Etablissement Public qui honore cet adjectif depuis de longues annees
REM   en rendant accessible au plus grand nombre un savoir et une experience de si grande qualite que seuls quelques-uns
REM   ont eu la patience (passion ?), la rigueur et la tenacite de les elaborer ...
REM   L'accueil du site de l'I.G.N. : www.ign.fr
REM   Les ressources techniques (algorithmes avec jeu d'essai) : http://geodesie.ign.fr/index.php?page=algorithmes
REM   
REM   Un site de base pour ceux qui veulent s'initier à la geodesie, http://sgcaf.free.fr/pages/techniques/ign_coordonnees.htm
REM
REM   Pour la partie traitant des coordonnees UTM, une visite du site du Professeur Steve DUTCH (Universite du Wisconsin) s'impose : http://www.uwgb.edu/dutchs/index.html
REM
REM   Merci à tous les anonymes qui ont publie sur internet un "petit bout d'algorithme", souvent bogue, mais qui une fois retravaille, m'ont permis de resoudre
REM   - et surtout de comprendre - la nature du problème ...
REM
REM   Pardon, enfin, à tous ceux qui, comme moi, sont herisses par l'obligation pour raisons de compatibilite de supprimer les accents de nos voyelles francophones ! 
REM
REM ==========================================

FUNCTION DMSDec(Param As STRING) As DOUBLE
   REM Conversion Degres, Minutes, Secondes en Degres Decimaux
   Dim ok As BOOLEAN
   Dim x As STRING
   Dim   neg As BOOLEAN
   Dim d,m,s,z As DOUBLE
   Dim DMS(3)
   neg=0
   x=UCase(Trim(Param))
   If Left(x,1)="-" THEN
      neg=1: x=Mid(x,2,255)
   EndIf
   If Right(x,1)="W" THEN neg=1
   If Right(x,1)="S" THEN neg=1
   If Right(x,1)="O" THEN neg=1
   ok=1
   While ok<>0
      ok=0
      If INSTR(x,Chr(8211))>0 THEN
         ok=1: x=REPLACE(x,Chr(8211),"-"): x=Trim(x)
      EndIf
      If INSTR(x,Chr(160))>0 THEN
         ok=1: x=REPLACE(x,Chr(160)," "): x=Trim(x)
      EndIf
      If INSTR(x,"-")>0 THEN
         ok=1: x=REPLACE(x,"-"," "): x=Trim(x)
      EndIf
      If INSTR(x,"N")>0 THEN
         ok=1: x=REPLACE(x,"N","")
      EndIf
      If INSTR(x,"S")>0 THEN
         ok=1: x=REPLACE(x,"S","")
      EndIf
      If INSTR(x,"E")>0 THEN
         ok=1: x=REPLACE(x,"E","")
      EndIf
      If INSTR(x,"W")>0 THEN
         ok=1: x=REPLACE(x,"W","")
      EndIf
      If INSTR(x,"O")>0 THEN
         ok=1: x=REPLACE(x,"O","")
      EndIf
      If INSTR(x," ")>0 THEN
         ok=1: x=REPLACE(x," ","\")
      EndIf
      If INSTR(x,"°")>0 THEN
         ok=1: x=REPLACE(x,"°","\")
      EndIf
      If INSTR(x,"'")>0 THEN
         ok=1: x=REPLACE(x,"'","\")
      EndIf
      If INSTR(x,Chr(34))>0 THEN
         ok=1: x=REPLACE(x,Chr(34),"\")
      EndIf
      If INSTR(x,"-")>0 THEN
         ok=1: x=REPLACE(x,"-","\")
      EndIf
      If INSTR(x,"+")>0 THEN
         ok=1: x=REPLACE(x,"+","\")
      EndIf
      If INSTR(x,"/")>0 THEN
         ok=1: x=REPLACE(x,"/","\")
      EndIf
      If INSTR(x,",")>0 THEN
         ok=1: x=REPLACE(x,",",".")
      EndIf
      If INSTR(x,"\\")>0 THEN
         ok=1: x=REPLACE(x,"\\","\")
      EndIf
   Wend
   DMS=SPLIT(x,"\"): d=Val(DMS(0)): m=Val(DMS(1)): s=Val(DMS(2)): z=d+m/60+s/3600
   If neg THEN z=-z
   DMSDec=z
End FUNCTION

FUNCTION DecDMS(P1 As DOUBLE, OPTIONAL P2 As INTEGER) As STRING
   REM Conversion Degres Decimaux en Degres, Minutes, Secondes
   Dim x As STRING
   Dim   neg As BOOLEAN
   Dim d,m,s,z As DOUBLE
   z=P1
   If z<0 THEN
      z=Abs(z): neg=1
   EndIf
   d=Int(z): z=z-d: z=z*60
   m=Int(z): z=z-m: z=z*60
   s=Int(z*1000+0.5)/1000
   x=d & "°" & m & "'" & s & Chr(34)
   If ISMISSING(P2) THEN
      If neg THEN x="-" & x
   Else
      If P2=0 THEN
         If neg THEN x=x & " S" Else x=x &" N"
      Else
         If neg THEN x=x & " W" Else x=x &" E"
      EndIf
   EndIf
   DecDMS=x
End FUNCTION

FUNCTION WGSLatNTF(GLat As DOUBLE, GLon As DOUBLE) As DOUBLE
   REM Calcul de la Latitude NTF à partir des coordonnees geographiques WGS84 en Degres decimaux
   REM GLat, GLon en Degres decimaux
   REM
   Dim   Phi, Lambda, h, AA, aNTF, aWGS, b, bNTF, bWGS, eNTF, eWGS, e2, v, X, Y, Z, p, r, f,u As DOUBLE
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------
   REM h = hauteur sur l'ellipsoïde - mis à zero
   Phi = GLat*Deg2Rad: Lambda = GLon*Deg2Rad: h = 0
   REM aNTF = Demi Grand-Axe de l'ellipsoïde de Clarke 1880
   aNTF = 6378249.2
   REM bNTF = Demi Petit-Axe de l'ellipsoïde de Clarke 1880
   bNTF = 6356515
   REM eNTF = Premiere Excentricite de l'ellipsoïde de Clarke 1880
   eNTF = Sqr(1 - (bNTF / aNTF)*(bNTF / aNTF))
   REM aWGS = Demi Grand-Axe de l'ellipsoïde WGS84
   aWGS = 6378137
   REM bWGS = Demi Petit-Axe de l'ellipsoïde WGS84
   bWGS = 6356752.314
   REM eWGS = Premiere Excentricite de l'ellipsoïde WGS84
   eWGS = Sqr(1 - ((bWGS / aWGS)*(bWGS / aWGS)))
   AA = aWGS: e2 = eWGS * eWGS: f = (AA - b) / AA
   v = AA / Sqr(1 - e2 * Sin(Phi) * Sin(Phi))
   REM Coordonnee Geocentriques WGS84 : X, Y, Z
   X = (v + h) * Cos(Phi) * Cos(Lambda)
   Y = (v + h) * Cos(Phi) * Sin(Lambda)
   Z = ((1 - e2) * v + h) * Sin(Phi)
   REM Changement de referentiel WGS84 vers NTF
   X = X + 168: Y = Y + 60: Z = Z - 320: AA = aNTF: b = bNTF: e2 = eNTF * eNTF: f = (AA - b) / AA
   p = Sqr(X * X + Y * Y): r = p + Z * Z
   u = ATN((Z / p) * ((1 - f) + (e2 * AA / r)))
   Phi = ATN((Z * (1 - f) + e2 * AA * Sin(u) * Sin(u) * Sin(u)) / ((1 - f) * (p - e2 * AA * Cos(u) * Cos(u) * Cos(u))))
   WGSLatNTF=Phi*Rad2Deg
End FUNCTION

FUNCTION WGSLonNTF(GLat As DOUBLE,GLon As DOUBLE) As DOUBLE
   REM Calcul de la Longitude NTF à partir des coordonnees geographiques WGS84 en Degres decimaux
   REM GLat, GLon en Degres decimaux
   REM
   Dim Phi, Lambda, h, AA, aWGS, b, bWGS, eWGS, e2, v, X, Y, Z As DOUBLE
   Dim MeridParis As DOUBLE
   MeridParis=2.33722916666667*PI()/180
   REM   2°20'14,025" E de Greenwich
      
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------
   REM h = hauteur sur l'ellipsoïde - mis à zero
   Phi = GLat*Deg2Rad: Lambda = GLon*Deg2Rad: h = 0
   REM   aWGS = Demi Grand-Axe de l'ellipsoïde WGS84
   aWGS = 6378137
   REM   bWGS = Demi Petit-Axe de l'ellipsoïde WGS84
   bWGS = 6356752.314
   REM   eWGS = Premiere Excentricite de l'ellipsoïde WGS84
   eWGS = Sqr(1 - ((bWGS / aWGS)*(bWGS / aWGS)))
   AA = aWGS: e2 = eWGS * eWGS
   v = AA / Sqr(1 - e2 * Sin(Phi) * Sin(Phi))
   REM   Coordonnee Geocentriques WGS84 : X, Y, Z
   X = (v + h) * Cos(Phi) * Cos(Lambda)
   Y = (v + h) * Cos(Phi) * Sin(Lambda)
   Z = ((1 - e2) * v + h) * Sin(Phi)
   REM   Changement de referentiel WGS84 vers NTF
   X = X + 168: Y = Y + 60: Z = Z - 320
   Lambda = ATN(Y / X)
   Lambda=Lambda-MeridParis
   If X < 0 THEN Lambda = Lambda+PI()
   Lambda=Lambda*Rad2Deg
   WGSLonNTF=Lambda
End FUNCTION

FUNCTION NTFLatWGS(GLat As DOUBLE,GLon As DOUBLE) As DOUBLE
   REM Calcul de la Latitude WGS84 à partir des coordonnees geographiques NTF en Degres decimaux
   REM GLat, GLon en Degres decimaux
   REM
   Dim Phi, Lambda, h, AA, aNTF, aWGS, b, bNTF, bWGS, eNTF, eWGS, e2, v, X, Y, Z, p, r, f, u As DOUBLE
   REM h = hauteur sur l'ellipsoïde - mis à zero
   Dim MeridParis As DOUBLE
   MeridParis=2.33722916666667
   REM   2°20'14,025" E de Greenwich
      
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   Phi = GLat*Deg2Rad: Lambda = (GLon+MeridParis)*Deg2Rad: h = 0
   REM aNTF = Demi Grand-Axe de l'ellipsoïde de Clarke 1880
   aNTF = 6378249.2
   REM bNTF = Demi Petit-Axe de l'ellipsoïde de Clarke 1880
   bNTF = 6356515
   REM eNTF = Premiere Excentricite de l'ellipsoïde de Clarke 1880
   eNTF = Sqr(1 - (bNTF / aNTF)*(bNTF / aNTF))
   REM aWGS = Demi Grand-Axe de l'ellipsoïde WGS84
   aWGS = 6378137
   REM bWGS = Demi Petit-Axe de l'ellipsoïde WGS84
   bWGS = 6356752.314
   REM eWGS = Premiere Excentricite de l'ellipsoïde WGS84
   eWGS = Sqr(1 - (bWGS / aWGS)*(bWGS / aWGS))
   AA = aNTF: e2 = eNTF * eNTF: f = (AA - b) / AA
   v = AA / Sqr(1 - e2 * Sin(Phi) * Sin(Phi))
   REM Coordonnee Geocentriques NTF : X, Y, Z
   X = (v + h) * Cos(Phi) * Cos(Lambda)
   Y = (v + h) * Cos(Phi) * Sin(Lambda)
   Z = ((1 - e2) * v + h) * Sin(Phi)
   REM Changement de referentiel NGF vers WGS84
   X = X - 168: Y = Y - 60: Z = Z + 320: AA = aWGS: b = bWGS: e2 = eWGS * eWGS: f = (AA - b) / AA
   p = Sqr(X * X + Y * Y): r = p + Z * Z
   u = ATN((Z / p) * ((1 - f) + (e2 * AA / r)))
   Phi = ATN((Z * (1 - f) + e2 * AA * Sin(u) * Sin(u) * Sin(u)) / ((1 - f) * (p - e2 * AA * Cos(u) * Cos(u) * Cos(u))))
   NTFLatWGS=Phi*Rad2Deg
End FUNCTION

FUNCTION NTFLonWGS(GLat As DOUBLE,GLon As DOUBLE) As DOUBLE
   REM Calcul de la Longitude WGS84 à partir des coordonnees geographiques NTF en Degres decimaux
   REM GLat, GLon en Degres decimaux
   REM
   Dim Phi, Lambda, h, AA, aNTF, b, bNTF, eNTF, e2, v, X, Y, Z, f, u As DOUBLE
   Dim MeridParis As DOUBLE
   MeridParis=2.33722916666667
   REM   2°20'14,025" E de Greenwich
      
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
      
   REM h = hauteur sur l'ellipsoïde - mis à zero
   Phi = GLat*Deg2Rad: Lambda = (GLon+MeridParis)*Deg2Rad: h = 0
   REM aNTF = Demi Grand-Axe de l'ellipsoïde de Clarke 1880
   aNTF = 6378249.2
   REM bNTF = Demi Petit-Axe de l'ellipsoïde de Clarke 1880
   bNTF = 6356515
   REM eNTF = Premiere Excentricite de l'ellipsoïde de Clarke 1880
   eNTF= Sqr(1 - (bNTF / aNTF)*(bNTF / aNTF))
   AA = aNTF: e2 = eNTF * eNTF: f = (AA - b) / AA
   v = AA / Sqr(1 - e2 * Sin(Phi) * Sin(Phi))
   REM Coordonnee Geocentriques NTF : X, Y, Z
   X = (v + h) * Cos(Phi) * Cos(Lambda)
   Y = (v + h) * Cos(Phi) * Sin(Lambda)
   Z = ((1 - e2) * v + h) * Sin(Phi)
   REM Changement de referentiel NTF vers WGS84
   X = X - 168: Y = Y - 60: Z = Z + 320
   Lambda = ATN(Y / X)
   If X < 0 THEN Lambda = Lambda + PI()
   NTFLonWGS=Lambda*Rad2Deg
End FUNCTION

FUNCTION LatIsom(LatitDec As DOUBLE, OPTIONAL PremExcEllips As DOUBLE) As DOUBLE
   REM LatitDec : Latitude Geographique en Degres Decimaux
   REM PremExEllips : Premiere excentricite de l'ellipsoïde (par Defaut, IAG/GRS 1980 pour WGS84)
   Dim   Clarke1880 As DOUBLE
   Dim   IagGrs80 As DOUBLE
   Dim LatitRad,s,s2,L As DOUBLE
   Clarke1880=0.08248325676342
   IagGrs80=0.08182297965731
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   LatitRad=LatitDec*Deg2Rad
   If ISMISSING(PremExcEllips) THEN PremExcEllips=IagGrs80
   s = PremExcEllips * Sin(LatitRad)
   s=(1-s)/(1+s)
   s=Log(s)
   s=s*(PremExcEllips / 2)
   s=Exp(s)
   s2 = Tan(PI()/4+LatitRad / 2)
   s=s2 * s
   L=Log(s)
   LatIsom=L
End FUNCTION

FUNCTION IsomLat(LatitIsom As DOUBLE, OPTIONAL PremExcEllips As DOUBLE, OPTIONAL TolConverg As DOUBLE) As DOUBLE
   REM PremExEllips : Premiere excentricite de l'ellipsoïde (par Defaut, IAG/GRS 1980 pour WGS84)
   Dim   Clarke1880 As DOUBLE
   Dim   IagGrs80 As DOUBLE
   Clarke1880=0.08248325676342
   IagGrs80=0.08182297965731
   Dim s0, s1, d, EL, rPremEx, rTolConv As DOUBLE
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   EL=Exp(LatitIsom)
   If ISMISSING(PremExcEllips) THEN rPremEx=IagGrs80 Else rPremEx=PremExcEllips
   If ISMISSING(TolConverg) THEN rTolConv=0.00000000001 Else rTolConv=TolConverg
   s1 = 2 * ATN(EL)-PI()/2
   While (Abs(s0 - s1)) > TolConverg
      s0 = s1
rem      d = ((1 + PremExcEllips * Sin(s0)) / (1 - PremExcEllips * Sin(s0)))^(PremExcEllips / 2)
      d = Log(((1 + PremExcEllips * Sin(s0)) / (1 - PremExcEllips * Sin(s0))))
      d = Exp(d*(PremExcEllips / 2))
      s1 = 2 * ATN(d*el)-PI()/2
   Wend
   IsomLat=s1*Rad2Deg
End FUNCTION

FUNCTION NTFYLamb(GLat As DOUBLE,GLon As DOUBLE,OPTIONAL TypLamb As INTEGER) As DOUBLE
   Dim LIsom, n, c, Xs, Ys, Y As DOUBLE
   Dim Clarke1880 As DOUBLE
   Dim TL As INTEGER
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   If ISMISSING(TypLamb) THEN TL=0 Else TL=TypLamb
   Clarke1880=0.08248325676342
   LIsom=LatIsom(GLat,Clarke1880)
   REM   Lambert II etendu par defaut
   Select Case TL
      Case 1
         n = 0.7604059656: c = 11603796.98: Xs = 600000: Ys = 5657616.674
      Case 2
         n = 0.7289686274: c = 11745793.39: Xs = 600000: Ys = 6199695.768
      Case 3
         n = 0.6959127966: c = 11947992.52: Xs = 600000: Ys = 6791905.085
      Case 4
         n = 0.6712679322: c = 12136281.99: Xs = 234.358: Ys = 7239161.542
      Case Else
         n = 0.7289686274: c = 11745793.39: Xs = 600000: Ys = 8199695.768
   End Select
   Y=Ys - c * Exp(-n * LIsom) * Cos(n * GLon*Deg2Rad)
   Y=Int(Y+0.5)
   NTFYLamb=Y
End FUNCTION

FUNCTION NTFXLamb(GLat As DOUBLE,GLon As DOUBLE,OPTIONAL TypLamb As INTEGER) As DOUBLE
   Dim LIsom, n, c, Xs, Ys, X As DOUBLE
   Dim Clarke1880 As DOUBLE
   Dim TL As INTEGER
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   If ISMISSING(TypLamb) THEN TL=0 Else TL=TypLamb
   Clarke1880=0.08248325676342
   LIsom=LatIsom(GLat,Clarke1880)
   REM   Lambert II etendu par defaut
   Select Case TypLamb
      Case 1
         n = 0.7604059656: c = 11603796.98: Xs = 600000: Ys = 5657616.674
      Case 2
         n = 0.7289686274: c = 11745793.39: Xs = 600000: Ys = 6199695.768
      Case 3
         n = 0.6959127966: c = 11947992.52: Xs = 600000: Ys = 6791905.085
      Case 4
         n = 0.6712679322: c = 12136281.99: Xs = 234.358: Ys = 7239161.542
      Case Else
         n = 0.7289686274: c = 11745793.39: Xs = 600000: Ys = 8199695.768
   End Select
   X=Xs + c * Exp(-n * LIsom) * Sin(n * GLon*Deg2Rad)
   X=Int(X+0.5)
   NTFXLamb=X
End FUNCTION

FUNCTION LambLatNTF(XLambert As DOUBLE,YLambert As DOUBLE,TypLamb As INTEGER) As DOUBLE
   Dim Clarke1880 As DOUBLE
   Dim e,LY, R, L, n, c, Xs, Ys As DOUBLE
   Dim TL As INTEGER
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   If ISMISSING(TypLamb) THEN TL=0 Else TL=TypLamb
   Clarke1880=0.08248325676342
   e=Clarke1880
   REM   Lambert II etendu par defaut
   Select Case TypLamb
      Case 1
         n = 0.7604059656: c = 11603796.98: Xs = 600000: Ys = 5657616.674
      Case 2
         n = 0.7289686274: c = 11745793.39: Xs = 600000: Ys = 6199695.768
      Case 3
         n = 0.6959127966: c = 11947992.52: Xs = 600000: Ys = 6791905.085
      Case 4
         n = 0.6712679322: c = 12136281.99: Xs = 234.358: Ys = 7239161.542
      Case Else
         n = 0.7289686274: c = 11745793.39: Xs = 600000: Ys = 8199695.768
   End Select
   R = Sqr((XLambert-Xs)*(XLambert-Xs) + (YLambert-Ys)*(YLambert-Ys))
   L = -1 / n * Log(Abs(R/c))
   LY=IsomLat(L, e, 0.00000000001)
   LambLatNTF=LY
End FUNCTION

FUNCTION LambLonNTF(XLambert As DOUBLE,YLambert As DOUBLE,TypLamb As INTEGER) As DOUBLE
Dim LX, R, g, n, c, Xs, Ys As DOUBLE
   Dim TL As INTEGER
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   Dim MeridParis As DOUBLE
   MeridParis=2.33722916666667
   REM   2°20'14,025" E de Greenwich
   If ISMISSING(TypLamb) THEN TL=0 Else TL=TypLamb
   REM   Lambert II etendu par defaut
   Select Case TypLamb
      Case 1
         n = 0.7604059656: c = 11603796.98: Xs = 600000: Ys = 5657616.674
      Case 2
         n = 0.7289686274: c = 11745793.39: Xs = 600000: Ys = 6199695.768
      Case 3
         n = 0.6959127966: c = 11947992.52: Xs = 600000: Ys = 6791905.085
      Case 4
         n = 0.6712679322: c = 12136281.99: Xs = 234.358: Ys = 7239161.542
      Case Else
         n = 0.7289686274: c = 11745793.39: Xs = 600000: Ys = 8199695.768
   End Select
   R = Sqr((XLambert-Xs)*(XLambert-Xs) + (YLambert-Ys)*(YLambert-Ys))
   g = ATN((XLambert-Xs) / (Ys-YLambert))
   LX= g / n
   LX=LX*Rad2Deg
   LambLonNTF=LX
End FUNCTION

FUNCTION WgsXL93(GLat As DOUBLE,GLon As DOUBLE) As DOUBLE
   Dim X93, IagGrs80, rAWGS, ra, re, rn, rcc, rys, rRgl, rRgl0, rRgl1, rRgl2 As DOUBLE
   Dim rlc,rl,rPhi,rPhi0,rPhi1,rPhi2,rx0,rRy0,rgN1,rgN2 As DOUBLE
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   REM   systeme WGS84
   IagGrs80=0.08182297965731
   rAWGS=6378137
   
   ra=rAWGS : REM   demi grand axe de l'ellipsoïde (m)
   re=IagGrs80: REM   premiere excentricite de l'ellipsoïde
   REM   parametres de projection
   rlc=3*deg2rad:    REM   Meridien central : Lambda0 = 3° Est de Greenwich pour Lambert93
   rPhi0=46.5*Deg2Rad:    REM   latitude Origine pour Lambert93
   rPhi1=44*Deg2Rad:    REM   1er parallele automecoïque
   rPhi2=49*Deg2Rad:    REM   2eme parallele automecoïque
   
   REM   coordonnees à l'origine
   rx0=700000: rRy0=6600000
   REM   coordonnees du point à traduire
   rPhi=GLat*Deg2Rad: rl=GLon*Deg2Rad
   REM   calcul des grandes normales
   rgN1=ra/Sqr(1-re*re*Sin(rPhi1)*Sin(rPhi1))
   rgN2=ra/Sqr(1-re*re*Sin(rPhi2)*Sin(rPhi2))
   REM calcul des latitudes isometriques
   rRgl1=Log(Tan(PI()/4+rPhi1/2)*Exp(Log((1-re*Sin(rPhi1))/(1+re*Sin(rPhi1)))*re/2)
   rRgl2=Log(Tan(PI()/4+rPhi2/2)*Exp(Log((1-re*Sin(rPhi2))/(1+re*Sin(rPhi2)))*re/2)
   rRgl0=Log(Tan(PI()/4+rPhi0/2)*Exp(Log((1-re*Sin(rPhi0))/(1+re*Sin(rPhi0)))*re/2)
   rRgl=Log(Tan(PI()/4+rPhi/2)*Exp(Log((1-re*Sin(rPhi))/(1+re*Sin(rPhi)))*re/2)
   REM   calcul de l'exposant de la projection
   rn=(Log((rgN2*Cos(rPhi2))/(rgN1*Cos(rPhi1))))/(rRgl1-rRgl2)
   REM   calcul de la constante de projection
   rcc=((rgN1*Cos(rPhi1))/rn)*Exp(rn*rRgl1)
   X93=rx0+rcc*Exp(-1*rn*rRgl)*Sin(rn*(rl-rlc))
   WgsXL93=X93
End FUNCTION

FUNCTION WgsYL93(GLat As DOUBLE,GLon As DOUBLE) As DOUBLE
   Dim Y93, IagGrs80, rAWGS, ra, re, rn, rcc, rys, rRgl, rRgl0, rRgl1, rRgl2 As DOUBLE
   Dim rlc,rl,rPhi,rPhi0,rPhi1,rPhi2,rx0,rRy0,rgN1,rgN2 As DOUBLE
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   REM   systeme WGS84
   IagGrs80=0.08182297965731
   rAWGS=6378137
   
   ra=rAWGS : REM   demi grand axe de l'ellipsoïde (m)
   re=IagGrs80: REM   premiere excentricite de l'ellipsoïde
   REM   parametres de projection
   rlc=3*deg2rad:    REM   Meridien central : Lambda0 = 3° Est de Greenwich pour Lambert93
   rPhi0=46.5*Deg2Rad:    REM   latitude Origine pour Lambert93
   rPhi1=44*Deg2Rad:    REM   1er parallele automecoïque
   rPhi2=49*Deg2Rad:    REM   2eme parallele automecoïque
   
   REM   coordonnees à l'origine
   rx0=700000: rRy0=6600000
   REM   coordonnees du point à traduire
   rPhi=GLat*Deg2Rad: rl=GLon*Deg2Rad
   REM   calcul des grandes normales
   rgN1=ra/Sqr(1-re*re*Sin(rPhi1)*Sin(rPhi1))
   rgN2=ra/Sqr(1-re*re*Sin(rPhi2)*Sin(rPhi2))
   REM calcul des latitudes isometriques
   rRgl1=Log(Tan(PI()/4+rPhi1/2)*Exp(Log((1-re*Sin(rPhi1))/(1+re*Sin(rPhi1)))*re/2)
   rRgl2=Log(Tan(PI()/4+rPhi2/2)*Exp(Log((1-re*Sin(rPhi2))/(1+re*Sin(rPhi2)))*re/2)
   rRgl0=Log(Tan(PI()/4+rPhi0/2)*Exp(Log((1-re*Sin(rPhi0))/(1+re*Sin(rPhi0)))*re/2)
   rRgl=Log(Tan(PI()/4+rPhi/2)*Exp(Log((1-re*Sin(rPhi))/(1+re*Sin(rPhi)))*re/2)
   REM   calcul de l'exposant de la projection
   rn=(Log((rgN2*Cos(rPhi2))/(rgN1*Cos(rPhi1))))/(rRgl1-rRgl2)
   REM   calcul de la constante de projection
   rcc=((rgN1*Cos(rPhi1))/rn)*Exp(rn*rRgl1)
   REM   calcul des coordonnees Lambert93
   X93=rx0+rcc*Exp(-1*rn*rRgl)*Sin(rn*(rl-rlc))
   rys=rRy0+rcc*Exp(-1*rn*rRgl0)
   Y93=rys-rcc*Exp(-1*rn*rRgl)*Cos(rn*(rl-rlc))
   WgsYL93=Y93
End FUNCTION

FUNCTION L93LatWGS(XLambert As DOUBLE,YLambert As DOUBLE) As DOUBLE
   Dim IagGrs80 As DOUBLE
   Dim e,LY, R, L, n, c, Xs, Ys As DOUBLE
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   IagGrs80=0.08182297965731
   n = 0.725607765:c = 11754255.426:Xs = 700000:Ys = 12655612.05: e=IagGrs80
   R = Sqr((XLambert-Xs)*(XLambert-Xs) + (YLambert-Ys)*(YLambert-Ys))
   L = -1 / n * Log(Abs(R/c))
   LY=IsomLat(L, e, 0.00000000001)
   L93LatWGS=LY
End FUNCTION

FUNCTION L93LonWGS(XLambert As DOUBLE,YLambert As DOUBLE) As DOUBLE
   Dim LX, R, g, n, c, Xs, Ys As DOUBLE
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   n = 0.725607765: c = 11754255.426: Xs = 700000: Ys = 12655612.05
   R = Sqr((XLambert-Xs)*(XLambert-Xs) + (YLambert-Ys)*(YLambert-Ys))
   g = ATN((XLambert-Xs) / (Ys-YLambert))
   LX= g / n
   LX=LX*Rad2Deg+3:    REM   Meridien central : Lambda0 = 3° Est de Greenwich pour Lambert93
   L93LonWGS=LX
End FUNCTION

FUNCTION WGS_FusZonUTM(LatWGS As DOUBLE, LonWGS As DOUBLE) As STRING
   Dim GabZon As STRING
   Dim LettreZon As STRING
   Dim NumFus As DOUBLE

   GabZon="AACCCDDEEFFGGHHJJKKLLMMNNPPQQRRSSTTUUVVWWXXXZZ"
   NumFus=Int((LonWGS+180)/6)+1
   LettreZon=Mid(GabZon,Int((LatWGS+92)/4)+1,1)
   WGS_FusZonUTM=FORMAT(NumFus,"00") & LettreZon
End FUNCTION

FUNCTION WGS_UTMEasting(LatWGS As DOUBLE, LonWGS As DOUBLE) As DOUBLE
   Dim aWGS, PrExcWGS As DOUBLE
   Dim PrExc2, PrExc4, PrExc6 As DOUBLE
   Dim Lambda, Phi, SinPhi, S2, CosPhi, C2, TanPhi, T2, C As DOUBLE
   Dim Lambda0 As DOUBLE
   Dim A, A2,A3, A4, A5, A6 As DOUBLE
   Dim k0, k1, k2, k3, k4 As DOUBLE
   Dim sPhi, vPhi, tPhi, E As DOUBLE
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   aWGS=6378137 : REM Rayon Equatorial de l'Ellipsoïde IAG/GRS80
   PrExcWGS=0.0818191910428152 : REM Premiere Excentricite de l'Ellipsoïde IAG/GRS80
   PrExc2=PrExcWGS*PrExcWGS : PrExc4=PrExc2*PrExc4 : PrExc6=PrExc4*PrExc4
   Lambda=LonWGS*Deg2Rad : Phi=LatWGS*Deg2Rad : SinPhi=Sin(Phi) : S2=SinPhi*SinPhi : CosPhi=Cos(Phi) : C2=CosPhi*CosPhi : TanPhi=Tan(Phi) : T2=TanPhi*TanPhi
   Lambda0=Int((LonWGS+180)/6)+1 : REM Fuseau UTM
   Lambda0=6*Lambda0-183 : REM Meridien de Reference de la Projection
   Lambda0=Lambda0*Deg2Rad : REM en Radians !
   A=(Lambda-Lambda0)*CosPhi
   A2=A*A : A3=A*A2 : A4=A*A3 : A5=A*A4 : A6=A*A5
   k0=0.9996
   k1=1-PrExc2/4-3*PrExc4/64-5*PrExc6/256
   k2=3*PrExc2/8+3*PrExc4/32+45*PrExc6/1024
   k3=15*PrExc4/256+45*PrExc6/1024
   k4=35*PrExc6/3072
   C=PrExc2/(1-PrExc2)*C2
   sPhi=k1*Phi-k2*Sin(2*Phi)+k3*Sin(4*Phi)-k4*Sin(6*Phi)
   vPhi=1/Sqr(1-PrExc2*S2)
   tPhi=TanPhi*(A2/2+(5-T2+9*C+4*C*C)*A4/24+(61-58*T2+T2*T2)*A6/720)
   E=500000+k0*aWGS*vPhi*(A+(1-T2+C)*A3/6+(5-18*T2+T2*T2)*A5/120)
   WGS_UTMEasting=Int(E*10+0.5)/10
End FUNCTION

FUNCTION WGS_UTMNorthing(LatWGS As DOUBLE, LonWGS As DOUBLE) As DOUBLE
   Dim aWGS, PrExcWGS As DOUBLE
   Dim PrExc2, PrExc4, PrExc6 As DOUBLE
   Dim Lambda, Phi, SinPhi, S2, CosPhi, C2, TanPhi, T2, C As DOUBLE
   Dim Lambda0 As DOUBLE
   Dim A, A2,A3, A4, A5, A6 As DOUBLE
   Dim k0, k1, k2, k3, k4 As DOUBLE
   Dim sPhi, vPhi, tPhi, N, Nz As DOUBLE
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------   
   aWGS=6378137 : REM Rayon Equatorial de l'Ellipsoïde IAG/GRS80
   PrExcWGS=0.0818191910428152 : REM Premiere Excentricite de l'Ellipsoïde IAG/GRS80
   PrExc2=PrExcWGS*PrExcWGS : PrExc4=PrExc2*PrExc2 : PrExc6=PrExc2*PrExc4
   Lambda=LonWGS*Deg2Rad : Phi=LatWGS*Deg2Rad : SinPhi=Sin(Phi) : S2=SinPhi*SinPhi : CosPhi=Cos(Phi) : C2=CosPhi*CosPhi : TanPhi=Tan(Phi) : T2=TanPhi*TanPhi
   Lambda0=Int((LonWGS+180)/6)+1 : REM Fuseau UTM
   Lambda0=6*Lambda0-183 : REM Meridien de Reference de la Projection
   Lambda0=Lambda0*Deg2Rad : REM en Radians !
   A=(Lambda-Lambda0)*CosPhi
   A2=A*A : A3=A*A2 : A4=A*A3 : A5=A*A4 : A6=A*A5
   k0=0.9996
   k1=1-PrExc2/4-3*PrExc4/64-5*PrExc6/256
   k2=3*PrExc2/8+3*PrExc4/32+45*PrExc6/1024
   k3=15*PrExc4/256+45*PrExc6/1024
   k4=35*PrExc6/3072
   C=PrExc2/(1-PrExc2)*C2
   sPhi=k1*Phi-k2*Sin(2*Phi)+k3*Sin(4*Phi)-k4*Sin(6*Phi)
   vPhi=1/Sqr(1-PrExc2*S2)
   tPhi=TanPhi*(A2/2+(5-T2+9*C+4*C*C)*A4/24+(61-58*T2+T2*T2)*A6/720)
   If LatWGS>0 THEN Nz=0 Else Nz=10000000
   N=Nz+k0*aWGS*(sPhi+vPhi*tPhi)
   WGS_UTMNorthing=Int(N*10+0.5)/10
End FUNCTION

FUNCTION UTM_LatWGS(UtmX As DOUBLE, UtmY As DOUBLE, UtmZone As STRING) As DOUBLE
   ' Calcul de la Latitude geographique WGS84 en degres decimaux à partir des coordonnees UTM
   ' d'apres Algorithme de Steve Dutch - University of GreenBay, Wisconsin
   ' Longitudes Positives vers l'Est, Negatives vers l'Ouest
   ' Latitudes Positives vers le Nord, Negatives vers le Sud
   '
   Dim LatD, East As DOUBLE
   Dim aWGS, bWGS, eWGS, e2, e4, e6, e1sq, k0, arc, mu, e1 As DOUBLE
   Dim ca, cb, cc, cd, phi1 As DOUBLE
   Dim Sin1, q0, t0, n0, r0, dd0 As DOUBLE
   Dim fact1, fact2, fact3, fact4 As DOUBLE
   Dim Hemisphere As STRING
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------
   If UCase(UtmZone)>"M" THEN Hemisphere="N" Else Hemisphere="S"
   REM   systeme WGS84
   REM aWGS = Demi Grand-Axe de l'ellipsoïde WGS84
   aWGS = 6378137
   REM bWGS = Demi Petit-Axe de l'ellipsoïde WGS84
   bWGS = 6356752.314
   REM eWGS = Premiere Excentricite de l'ellipsoïde WGS84
   eWGS = Sqr(1 - ((bWGS / aWGS)*(bWGS / aWGS))):e2=eWGS*eWGS:e4=e2*e2:e6=e4*e2
   East = UtmX
   If Hemisphere = "S" THEN LatD = 10000000 - UtmY Else LatD = UtmY
   e1sq = e2 / (1 - e2): k0 = 0.9996
   arc = LatD / k0
   mu = arc / (aWGS * (1 - e2 / 4 - 3 * e4 / 64 - 5 * e6 / 256))
   e1 = (1 - (1 - e2) ^ 0.5) / (1 + (1 - e2) ^ 0.5)
   ca = 3 * e1 / 2 - 27 * e1 ^ 3 / 32
   cb = 21 * e1 ^ 2 / 16 - 55 * e1 ^ 4 / 32
   cc = 151 * e1 ^ 3 / 96
   cd = 1097 * e1 ^ 4 / 512
   phi1 = mu + ca * Sin(2 * mu) + cb * Sin(4 * mu) + cc * Sin(6 * mu) + cd * Sin(8 * mu)
   q0 = e1sq * Cos(phi1) ^ 2: t0 = Tan(phi1) ^ 2
   n0 = aWGS / ((1 - (eWGS * Sin(phi1)) ^ 2) ^ (0.5))
   r0 = aWGS * (1 - eWGS * eWGS) / (1 - (eWGS * Sin(phi1)) ^ 2) ^ (3 / 2)
   dd0 = (500000 - East) / (n0 * k0)
   fact1 = n0 * Tan(phi1) / r0: fact2 = dd0 * dd0 / 2
   fact3 = (5 + 3 * t0 + 10 * q0 - 4 * q0 * q0 - 9 * e1sq) * dd0 ^ 4 / 24
   fact4 = (61 + 90 * t0 + 298 * q0 + 45 * t0 * t0 - 252 * e1sq - 3 * q0 * q0) * dd0 ^ 6 / 720
   If Hemisphere = "N" THEN
      UTM_LatWGS = (phi1 - fact1 * (fact2 + fact3 + fact4))*Rad2Deg
   Else
      UTM_LatWGS = (-(phi1 - fact1 * (fact2 + fact3 + fact4)))*Rad2Deg
   EndIf
End FUNCTION

FUNCTION UTM_LonWGS(UtmX As DOUBLE, UtmY As DOUBLE, UtmFuseau As INTEGER, UtmZone As STRING) As DOUBLE
   ' Calcul de la Longitude geographique WGS84 en degres decimaux à partir des coordonnees UTM
   ' d'apres Algorithme de Steve Dutch - University of GreenBay, Wisconsin
   ' Longitudes Positives vers l'Est, Negatives vers l'Ouest
   ' Latitudes Positives vers le Nord, Negatives vers le Sud
   '
   Dim LatD, East As DOUBLE
   Dim MCFuseau As INTEGER
   Dim aWGS, bWGS, eWGS, e2, e4, e6, e1sq, k0, arc, mu, e1 As DOUBLE
   Dim ca, cb, cc, cd, phi1 As DOUBLE
   Dim Sin1, q0, t0, n0, r0, dd0 As DOUBLE
   Dim fact1, fact2, fact3 As DOUBLE
   REM   ----------------------------------
   REM   Conversions Degres <---> Radians
   Dim Rad2Deg, Deg2Rad As DOUBLE
   Rad2Deg=180/PI(): Deg2Rad=PI()/180
   REM   ----------------------------------
   REM   systeme WGS84
   REM aWGS = Demi Grand-Axe de l'ellipsoïde WGS84
   aWGS = 6378137
   REM bWGS = Demi Petit-Axe de l'ellipsoïde WGS84
   bWGS = 6356752.314
   REM eWGS = Premiere Excentricite de l'ellipsoïde WGS84
   eWGS = Sqr(1 - ((bWGS / aWGS)*(bWGS / aWGS))):e2=eWGS*eWGS:e4=e2*e2:e6=e4*e2
   East = UtmX
   If UCase(UtmZone) = "S" THEN LatD = 10000000 - UtmY Else LatD = UtmY
   ' Meridien Central du Fuseau
   MCFuseau = 6 * UtmFuseau - 183
   e1sq = e2 / (1 - e2): k0 = 0.9996
   arc = LatD / k0
   mu = arc / (aWGS * (1 - e2 / 4 - 3 * e4 / 64 - 5 * e6 / 256))
   e1 = (1 - (1 - e2) ^ 0.5) / (1 + (1 - e2) ^ 0.5)
   ca = 3 * e1 / 2 - 27 * e1 ^ 3 / 32
   cb = 21 * e1 ^ 2 / 16 - 55 * e1 ^ 4 / 32
   cc = 151 * e1 ^ 3 / 96
   cd = 1097 * e1 ^ 4 / 512
   phi1 = mu + ca * Sin(2 * mu) + cb * Sin(4 * mu) + cc * Sin(6 * mu) + cd * Sin(8 * mu)
   q0 = e1sq * Cos(phi1) ^ 2: t0 = Tan(phi1) ^ 2
   n0 = aWGS / ((1 - (eWGS * Sin(phi1)) ^ 2) ^ (0.5))
   r0 = aWGS * (1 - e2) / (1 - (eWGS * Sin(phi1)) ^ 2) ^ (3 / 2)
   dd0 = (500000 - East) / (n0 * k0)
   fact1 = dd0
   fact2 = (1 + 2 * t0 + q0) * dd0 ^ 3 / 6
   fact3 = (5 - 2 * q0 + 28 * t0 - 3 * q0 ^ 2 + 8 * e1sq + 24 * t0 ^ 2) * dd0 ^ 5 / 120
   UTM_LonWGS = ((MCFuseau * Deg2Rad) - (fact1 - fact2 + fact3) / Cos(phi1))*Rad2Deg
   End FUNCTION
   
   
 ;#######################################################################################################################

;voir ce document "Geo.ods" Calc open office sur ce site: http://forum.openoffice.org/fr/forum/viewtopic.php?f=8&t=26883

;#################### EXEMPLE D UTILISATION ############################################################################
;	Conversion Degrés, Minutes, Secondes => Degrés Décimaux :			DMS_Dec	Commentaires
;	Degrés Dddd	Degrés Décimaux	Degrés Décimaux		Param1: chaîne à convertir en réel
;	Latitude	Longitude	Latitude	Longitude	
;	42°5'4.321"	2-21 – 32,0126	#NOM ?	#NOM ?	
;	42°5'4.321" N	2  21 32,0126 W	#NOM ?	#NOM ?	
;	42°5'4,321" S	2/21/32.0126 O	#NOM ?	#NOM ?	
;					
;	Conversion Degrés Décimaux => Degrés, Minutes, Secondes : 			Dec_DMS	
;	Degrés Décimaux	Degrés Décimaux	Degrés Dddd		Param1: réel à convertir en chaîne
;	Latitude	Longitude	Latitude	Longitude	Param2: 0 pour Latitude, 1 pour Longitude
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;					
;	Déterminer les Coordonnées NTF à partir de Coordonnées WGS84 :			WGS_LatNTF
;WGS_LonNTF	
;	Degrés Dddd WGS84	Degrés Décimaux	Degrés Décimaux WGS84		Param1: Latitude WGS84 en Degrés Décimaux
;	Latitude	Longitude	Latitude	Longitude	Param2: Longitude WGS84 en Degrés Décimaux
;	43°36'40,1" N	3°52'20,6" E	#NOM ?	#NOM ?	Méridien Origine : Greenwich
;	Degrés Dddd NTF	Degrés Décimaux	Degrés Décimaux NTF		
;	Latitude	Longitude	Latitude	Longitude	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Paris (2°20'14,025" E de Greenwich)
;					
;	Déterminer les Coordonnées WGS84 à partir de Coordonnées NTF :			NTF_LatWGS
;NTF_LonWGS	
;	Degrés Dddd NTF	Degrés Décimaux	Degrés Décimaux NTF		Param1: Latitude NTF en Degrés Décimaux
;	Latitude	Longitude	Latitude	Longitude	Param2: Longitude NTF en Degrés Décimaux
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Paris (2°20'14,025" E de Greenwich)
;	Degrés Dddd WGS84	Degrés Décimaux	Degrés Décimaux WGS84		
;	Latitude	Longitude	Latitude	Longitude	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Greenwich
					
;	Déterminer les Projections Lambert (I à IV + 2E) à partir de Coordonnées NTF :			NTF_XLamb
;NTF_YLamb	Si on dispose de coordonnées WGS84, les convertir en NTF à  l'aide des fonctions WGS_LatNTF et WGS_LonNTF
;	Degrés Dddd	Degrés Décimaux	Degrés Décimaux		Param1: Latitude NTF en Degrés Décimaux
;	Latitude	Longitude	Latitude	Longitude	Param2: Longitude NTF en Degrés Décimaux
;WGS84	49°42'25.2" N	0°11'33.8" E	#NOM ?	#NOM ?	n° Zone Lambert (1,2,3 ou 4), 0 ou absent pour II étendu
;NTF	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Aiguille d'Etretat (Normandie)
;1	Lambert I (Nord)	Degrés Décimaux	Lambert II étendu (France)		
;	Y – Northing	X – Easting	Y – Northing	X – Easting	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;	Degrés Dddd	Degrés Décimaux	Degrés Décimaux		
;	Latitude	Longitude	Latitude	Longitude	
;WGS84	47°12'42.1" N	1°32'09.4" O	#NOM ?	#NOM ?	
;NTF	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Pont Willy Brandt (Nantes)
;2	Lambert II (Centre)	Degrés Décimaux	Lambert II étendu (France)		
;	Y – Northing	X – Easting	Y – Northing	X – Easting	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;	Degrés Dddd	Degrés Décimaux	Degrés Décimaux		
;	Latitude	Longitude	Latitude	Longitude	
;WGS84	42°41'46.3" N	2°52'47.2" E	#NOM ?	#NOM ?	
;NTF	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Gare de Perpignan (Centre du Monde)
;3	Lambert III (Sud)	Degrés Décimaux	Lambert II étendu (France)		
;	Y – Northing	X – Easting	Y – Northing	X – Easting	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;	Degrés Dddd	Degrés Décimaux	Degrés Décimaux		
;	Latitude	Longitude	Latitude	Longitude	
;WGS84	41°55'09.6" N	8°44'37.9" E	#NOM ?	#NOM ?	
;NTF	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Port d'Ajaccio (Corse)
;4	Lambert IV (Corse)	Degrés Décimaux	Lambert II étendu (France)		
;	Y – Northing	X – Easting	Y – Northing	X – Easting	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;					
;	Déterminer les Coordonnées NTF à partir de Projections Lambert (I à IV + 2E) :			Lamb_LatNTF
;Lamb_LonNTF	Si on désire les coordonnées WGS84, convertir le résultat NTF en WGS84 à  l'aide des fonctions NTF_LatWGS et NTF_LonWGS
;1	Lambert I (Nord)	Degrés Décimaux	Lambert II étendu (France)		Param1: X Lambert (Easting) en entier
;	Y – Northing	X – Easting	Y – Northing	X – Easting	Param2: Y Lambert (Northing) en entier SANS ZONE
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	n° Zone Lambert (1,2,3 ou 4), 0 ou absent pour II étendu
;	Degrés Dddd	Degrés Décimaux	Degrés Décimaux		
;	Latitude	Longitude	Latitude	Longitude	Aiguille d'Etretat (Normandie)
;NTF	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Paris (2°20'14,025" E de Greenwich)
;WGS84	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Greenwich
;2	Lambert II (Centre)	Degrés Décimaux	Lambert II étendu (France)		
;	Y – Northing	X – Easting	Y – Northing	X – Easting	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;	Degrés Dddd	Degrés Décimaux	Degrés Décimaux		
;	Latitude	Longitude	Latitude	Longitude	Pont Willy Brandt (Nantes)
;NTF	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Paris (2°20'14,025" E de Greenwich)
;WGS84	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Greenwich
;3	Lambert III (Sud)	Degrés Décimaux	Lambert II étendu (France)		
;	Y – Northing	X – Easting	Y – Northing	X – Easting	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;	Degrés Dddd	Degrés Décimaux	Degrés Décimaux		
;	Latitude	Longitude	Latitude	Longitude	Gare de Perpignan (Centre du Monde)
;NTF	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Paris (2°20'14,025" E de Greenwich)
;WGS84	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Greenwich
;4	Lambert IV (Corse)	Degrés Décimaux	Lambert II étendu (France)		
;	Y – Northing	X – Easting	Y – Northing	X – Easting	
;	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;	Degrés Dddd	Degrés Décimaux	Degrés Décimaux		
;	Latitude	Longitude	Latitude	Longitude	Port d'Ajaccio (Corse)
;NTF	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Paris (2°20'14,025" E de Greenwich)
;WGS84	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Méridien Origine : Greenwich
;					
;	Déterminer la projection Lambert 93 à partir de Coordonnées WGS84 :			WGS_XL93
;WGS_YL93	
;	Degrés Dddd WGS84	Degrés Décimaux	Degrés Décimaux WGS84		Param1: Latitude WGS84 en Degrés Décimaux
;	Latitude	Longitude	Latitude	Longitude	Param2: Longitude WGS84 en Degrés Décimaux
;WGS84	43°36'40,1" N	3°52'20,6" E	#NOM ?	#NOM ?	Méridien Origine : Greenwich
;	Lambert 93	Degrés Décimaux			
;	Y – Northing	X – Easting			Arc de Triomphe du Peyrou (Montpellier)
;L93	#NOM ?	#NOM ?			
;					
;	Déterminer les Coordonnées WGS84 à partir de la projection Lambert 93  :			L93_LatWGS
;L93_LonWGS	
;	Lambert 93	Degrés Décimaux			Param1: X Lambert (Easting) en entier
;	Y – Northing	X – Easting			Param2: Y Lambert (Northing) en entier
;L93	#NOM ?	#NOM ?			
;	Degrés Dddd WGS84	Degrés Décimaux	Degrés Décimaux WGS84		
;	Latitude	Longitude	Latitude	Longitude	Arc de Triomphe du Peyrou (Montpellier)
;WGS84	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;					
;	Déterminer les Coordonnées UTM à partir de Coordonnées WGS84 :			WGS_FUSZONUTM
;WGS_UTMEASTING
;WGS_UTMNORTHING	
;	Degrés Dddd WGS84	Degrés Décimaux	Degrés Décimaux WGS84		Un grand merci au Professeur Steve DUTCH pour son travail !
;	Latitude	Longitude	Latitude	Longitude	Special thanks To Prof. Steve DUTCH For his work !
;WGS84	43°36'40" N	3°52'21" E	#NOM ?	#NOM ?	lien vers le site du Prof. Steve Dutch
;	Fuseau UTM	Zone UTM	Easting UTM	Northing UTM	
;UTM	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
;					
;	Déterminer les Coordonnées WGS84 à partir de Coordonnées UTM :			UTM_LatWGS
;UTM_LonWGS	
;	Fuseau UTM	Zone UTM	Easting UTM	Northing UTM	Un grand merci au Professeur Steve DUTCH pour son travail !
;UTM	#NOM ?	#NOM ?	#NOM ?	#NOM ?	Special thanks To Prof. Steve DUTCH For his work !
;	Degrés Décimaux WGS84		Degrés Dddd WGS84		lien vers le site du Prof. Steve Dutch
;	Latitude	Longitude	Latitude	Longitude	
;WGS84	#NOM ?	#NOM ?	#NOM ?	#NOM ?	
  

[MetalOS]

Cool merci Kernadec