Je suis tombé par hasard sur un message de kernadec sur ce forum dans lequel il avait la gentillesse de citer un travail réalisé sur le forum OpenOffice.
J'ai réalisé un petit utilitaire en PureBasic qui effectue les conversions entre systèmes géographiques ou projections cartographiques et permet de visualiser sur le web les points dont on fournit les coordonnées.
Peut-être cela intéressera-t-il quelqu'un ...
PK1157
Voici le source (GeKo.pb), je posterai le fichier d'aide GeKo_Aide.htm et le fichier de test GekoTest.csv dans un message suivant (taille !)
Edition le 10 juin 2015 18:53 : correction caractère accentué<=>entité HTML et lancement direct du navigateur pour les fichiers HTML produits.
Code : Tout sélectionner
Global NbParam.l=CountProgramParameters()
Global Dim LCmd$(7)
Global Dim GeoSys$(14)
Global Dim GeoPrm.l(14)
Global InSys$=""
Global OutSys$=""
Global NomPoint$=""
Global Result$=""
Global ResWeb$=""
Global$ RC$=Chr(13)+Chr(10)
Global$ TAB$=Chr(9)
Global.d LatWgs.d=0
Global.d LonWgs.d=0
Global.d LatNtf.d=0
Global.d LonNtf.d=0
Global.d XProj.d=0
Global.d YProj.d=0
Global NumFus=0
Global$ LettreZon$=""
Global$ WEB1$=""
Global$ WEB2$=""
Global$ NomWeb$=""
Global OS=0
Global$ Sys$=""
;Conversions Degrés <---> Radians
Global.d Rad2Deg.d=180/#PI
Global.d Deg2Rad.d=#PI/180
;système WGS84 (Ellipsoïde WGS84)
Global.d aWGS.d=6378137:;Demi Grand-Axe de l'ellipsoïde WGS84 (m)
Global.d bWGS.d=6356752.314245179:;Demi Petit-Axe de l'ellipsoïde WGS84 (m)
Global.d eWGS.d=Sqr(1-((bWGS/aWGS)*(bWGS/aWGS))):;Première Excentricité de l'ellipsoïde WGS84
;système RGF93 (Ellipsoïde IAG-GRS80)
Global.d aRGF.d=6378137:;Demi Grand-Axe de l'ellipsoïde IAG-GRS80 (m)
Global.d bRGF.d=6356752.314140356:;Demi Petit-Axe de l'ellipsoïde IAG-GRS80 (m)
Global.d eRGF.d=Sqr(1-((bRGF/aRGF)*(bRGF/aRGF))):;Première Excentricité de l'ellipsoïde IAG-GRS80
;système NTF (Ellipsoïde Clarke 1880)
Global.d aNTF.d=6378249.2:;Demi Grand-Axe de l'ellipsoïde Clarke 1880 (m)
Global.d bNTF.d=6356515:;Demi Petit-Axe de l'ellipsoïde Clarke 1880 (m)
Global.d eNTF.d=Sqr(1-((bNTF/aNTF)*(bNTF/aNTF))):;Première Excentricité de l'ellipsoïde Clarke 1880
Global.d MeridParis.d=2.33722916666667:;2°20'14,025" E de Greenwich
Procedure InitGeko()
OS=OSVersion()
Select OS
Case #PB_OS_Windows_NT3_51 To #PB_OS_Windows_Future
Sys$="W"
Case #PB_OS_Linux_2_2 To #PB_OS_Linux_Future
Sys$="X"
Default
Sys$="N"
EndSelect
;Tableau des systèmes Géo/Carto utilisables
;RGC et HTML ne'sont utilisables qu'en sortie !
;CSV n'est utilisable qu'en entrée !
For i=0 To 7
LCmd$(i)=ProgramParameter(i)
If i<>1:LCmd$(i)=ReplaceString(LCmd$(i),",","."):LCmd$(i)=UCase(LCmd$(i)):EndIf
Next i
GeoSys$(0)="CSV":GeoPrm(0)=1:;Utilisable en entrée (SOURCE) seulement
GeoSys$(1)="WGS":GeoPrm(1)=2:;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(2)="WGSD":GeoPrm(2)=2:;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(3)="NTF":GeoPrm(3)=2: ;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(4)="NTFD":GeoPrm(4)=2:;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(5)="L1":GeoPrm(5)=2: ;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(6)="L2":GeoPrm(6)=2: ;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(7)="L2E":GeoPrm(7)=2: ;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(8)="L3":GeoPrm(8)=2:;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(9)="L4":GeoPrm(9)=2:;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(10)="L93":GeoPrm(10)=2:;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(11)="DFCI":GeoPrm(11)=1:;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(12)="UTM":GeoPrm(12)=4: ;Utilisable en entrée (SOURCE) comme en sortie (CIBLE)
GeoSys$(13)="RGC":GeoPrm(13)=0;Utilisable en sortie (CIBLE) seulement
GeoSys$(14)="WEB":GeoPrm(14)=0;Utilisable en sortie (CIBLE) seulement
WEB1$="<!DOCTYPE html PUBLIC "+Chr(34)+"-//W3C//DTD XHTML 1.1//FR"+Chr(34)+" "+Chr(34)+"http://www.w3.org/TR/xhtml11/DTD/xhtml11.dtd"+Chr(34)+">"
WEB1$=WEB1$+RC$+"<html xmlns="+Chr(34)+"http://www.w3.org/1999/xhtml"+Chr(34)+" xml:lang="+Chr(34)+"fr"+Chr(34)+" >"
WEB1$=WEB1$+RC$+" <head>"
WEB1$=WEB1$+RC$+" <title>GeKo</title>"
WEB1$=WEB1$+RC$+" <meta http-equiv="+Chr(34)+"Content-Type"+Chr(34)+" content="+Chr(34)+"text/html; charset=UTF-8"+Chr(34)+" />"
WEB1$=WEB1$+RC$+" <meta name="+Chr(34)+"author"+Chr(34)+" content="+Chr(34)+"PK1157"+Chr(34)+">"
WEB1$=WEB1$+RC$+" <meta name="+Chr(34)+"copyright"+Chr(34)+" content="+Chr(34)+"PK1157"+Chr(34)+">"
WEB1$=WEB1$+RC$+" <style type="+Chr(34)+"text/css"+Chr(34)+">"
WEB1$=WEB1$+RC$+" table {margin: 0 auto; font-size: normal; text-align:center;border:1px solid black;}"
WEB1$=WEB1$+RC$+" caption {color:#801020; font-weight:bold;}"
WEB1$=WEB1$+RC$+" h1 {text-align:center;border:1px solid black;}"
WEB1$=WEB1$+RC$+" h2 {text-align:center;}"
WEB1$=WEB1$+RC$+" h3 {text-align:center;}"
WEB1$=WEB1$+RC$+" h4 {text-decoration:underline; color:teal;}"
WEB1$=WEB1$+RC$+" h5 {background: #C0C0F0; margin-left:40px; font-size: small; font-style : italic; color:teal;}"
WEB1$=WEB1$+RC$+" </style>"
WEB1$=WEB1$+RC$+" </head>"
WEB1$=WEB1$+RC$+" <body>"
WEB1$=WEB1$+RC$+" <h1>GeKo</h1>"
WEB1$=WEB1$+RC$+" <h2>Conversion Géo/Cartographique<br>(<i>Geodesic Konversion ?</i>)<br>Utilitaire en ligne de commande</h2>"
WEB1$=WEB1$+RC$+" <div>"
WEB1$=WEB1$+RC$+" <table width="+Chr(34)+"1000"+Chr(34)+" cols="+Chr(34)+"6"+Chr(34)+" border="+Chr(34)+"1"+Chr(34)+">"
WEB1$=WEB1$+RC$+" <caption>Conversion GeKo</caption>"
WEB1$=WEB1$+RC$+" <tr bgcolor="+Chr(34)+"#80DF80"+Chr(34)+"><th>Point</th><th>Lat.WGS84</th><th>lon.WGS84</th><th>OpenStreetMap<SUP>®</SUP></th><th>GeoPortail I.G.N.<SUP>®</SUP></th><th>GoogleMaps<SUP>®</SUP></th></tr>"
WEB2$=" </table>"
WEB2$=WEB2$+RC$+" </div>"
WEB2$=WEB2$+RC$+" <h2><a href="+Chr(34)+"mailto:pierrepetit57@laposte.net"+Chr(34)+" title="+Chr(34)+"Heureux d'Être Utile
(et un mèl de remerciement fait toujours plaisir !)"+Chr(34)+">HEU ...</a></h2>"
WEB2$=WEB2$+RC$+" </body>"
WEB2$=WEB2$+RC$+"</html>"
EndProcedure
Procedure Aide()
AppDir$=GetPathPart(ProgramFilename())
CurDir$=GetCurrentDirectory()
PrintN("Paramètre(s) absent(s) ou incorrects")
If FileSize(CurDir$+"GeKo_Aide.htm")=-1
If FileSize(AppDir$+"GeKo_Aide.htm")=-1
CopyFile(AppDir$+"GeKo_Aide.htm",CurDir$+"GeKo_Aide.htm")
EndIf
EndIf
If FileSize(CurDir$+"GeKo_Aide.htm")<>-1
Select Sys$
Case "W" : RunProgram(CurDir$+"Geko_Aide.htm")
Case "X" : RunProgram("xdg-open",CurDir$+"Geko_Aide.htm","")
Default :PrintN("Consultez le fichier "+CurDir$+"GeKo_Aide.htm avec votre navigateur.")
EndSelect
Else
PrintN("Impossible d'ouvrir le fichier d'aide !")
EndIf
EndProcedure
Procedure.d LatIsom(LatitDec.d, PremExcEllips.d)
;LatitDec:Latitude Géographique en Degrés Décimaux
;PremExEllips:Première excentricité de l'ellipsoïde
LatitRad.d=LatitDec*Deg2Rad
s.d=PremExcEllips*Sin(LatitRad)
s=(1-s)/(1+s)
s=Log(s)
s=s*(PremExcEllips/2)
s=Exp(s)
s2.d=Tan(#PI/4+LatitRad/2)
s=s2*s
L.d=Log(s)
ProcedureReturn(L)
EndProcedure
Procedure.d IsomLat(LatitIsom.d, PremExcEllips.d)
;PremExEllips:Première excentricité de l'ellipsoïde (par Défaut, IAG/GRS 1980 pour WGS84)
; ;----------------------------------
EL.d=Exp(LatitIsom)
rTolConv.d=0.00000000001
s0.d=0:s1.d=2*ATan(EL)-#PI/2
While (Abs(s0-s1)) > TolConverg
s0=s1
d.d=Log(((1+PremExcEllips*Sin(s0))/(1-PremExcEllips*Sin(s0))))
d=Exp(d*(PremExcEllips/2))
s1=2*ATan(d*el)-#PI/2
Wend
s1=s1*Rad2Deg
ProcedureReturn(s1)
EndProcedure
Procedure.d T2N(Coord$)
;Conversion coordonnée Texte ([[D]D]D.[M]M.[S]S[.sss]) en Décimal(999.999999)
;4 groupes au plus, séparés par un point
Res.d=0:Deg.d=0:Min.d=0:Sec.d=0:Mil.d=0
Deg=Val(StringField(Coord$,1,"."))
Min=ValD(StringField(Coord$,2,"."))/60
Sec=ValD(StringField(Coord$,3,"."))/3600
Mil=ValD(Left(StringField(Coord$,4,".")+"000",3))/3600000
Res=Deg+Min+Sec+Mil
ProcedureReturn(Res)
EndProcedure
Procedure$ T2D(Coord$)
;Conversion coordonnée Texte ([[D]D]D.[M]M.[S]S[.sss]) en DMS (DDD° MM' SS.sss")
Res$=Coord$
Res$=ReplaceString(Res$,".","° ",#PB_String_NoCase,1,1)
Res$=ReplaceString(Res$,".","' ",#PB_String_NoCase,1,1)
Res$=Res$+"""
ProcedureReturn(Res$)
EndProcedure
Procedure$ N2T(Coord.d)
;Conversion coordonnée Décimales(999.999999) en Texte ([[D]D]D.[M]M.[S]S[.sss])
;4 groupes au plus, séparés par un point
ValIn.d=Coord
Res$="":Deg=0:Min=0:Sec=0:Mil=0
Deg=Int(ValIn):ValIn=(Valin-Deg)*60
Min=Int(ValIn):ValIn=(Valin-min)*60
Sec=Int(ValIn):ValIn=(Valin-Sec)*1000
Mil=Round(ValIn,#PB_Round_Nearest)
Res$=StrD(Deg)+"."+Right("00"+StrD(Min),2)+"."+Right("00"+StrD(Sec),2)+"."+Right("000"+StrD(Mil),3)
ProcedureReturn(Res$)
EndProcedure
Procedure WGS_NTF(GLat.d,GLon.d)
;Calcul des coordonnées NTF à partir des coordonnées géographiques WGS84 en Degrés décimaux
;GLat, GLon en Degrés décimaux
;h=hauteur sur l'ellipsoïde-mise à zéro
Phi.d=GLat*Deg2Rad:Lambda.d=GLon*Deg2Rad:h.d=0
AA.d=aWGS:e2.d=eWGS*eWGS
v.d=AA/Sqr(1-e2*Sin(Phi)*Sin(Phi))
;Coordonnée Géocentriques WGS84:X, Y, Z
X.d=(v+h)*Cos(Phi)*Cos(Lambda)
Y.d=(v+h)*Cos(Phi)*Sin(Lambda)
Z.d=((1-e2)*v+h)*Sin(Phi)
;Changement de référentiel WGS84 vers NTF
f.d=(aNTF-bNTF)/aNTF
X=X+168:Y=Y+60:Z=Z-320
AA=aNTF:b=bNTF:e2=eNTF*eNTF
p.d=Sqr(X*X+Y*Y):r.d=p+Z*Z
u.d=ATan((Z/p)*((1-f)+(e2*AA/r)))
Phi=ATan((Z*(1-f)+e2*AA*Sin(u)*Sin(u)*Sin(u))/((1-f)*(p-e2*AA*Cos(u)*Cos(u)*Cos(u))))
Phi=Phi*Rad2Deg
Lambda=ATan(Y/X)
Lambda=Lambda-Deg2Rad*MeridParis
If X < 0:Lambda=Lambda+#PI:EndIf
Lambda=Lambda*Rad2Deg
LatNtf=Phi:LonNtf=Lambda
EndProcedure
Procedure NTF_WGS (GLat.d,GLon.d)
;Calcul des coordonées WGS84 à partir des coordonnées géographiques NTF en Degrés décimaux
;GLat, GLon en Degrés décimaux
;h=hauteur sur l'ellipsoïde-mise à zéro
Phi.d=GLat*Deg2Rad:Lambda.d=(GLon+MeridParis)*Deg2Rad:h.d=0
AA.d=aNTF:e2.d=eNTF*eNTF:f.d=(AA-bNTF)/AA
v.d=AA/Sqr(1-e2*Sin(Phi)*Sin(Phi))
;Coordonnée Géocentriques NTF:X, Y, Z
X.d=(v+h)*Cos(Phi)*Cos(Lambda)
Y.d=(v+h)*Cos(Phi)*Sin(Lambda)
Z.d=((1-e2)*v+h)*Sin(Phi)
;Changement de référentiel NGF vers WGS84
X=X-168:Y.d=Y-60:Z.d=Z+320:AA=aWGS:b=bWGS:e2=eWGS*eWGS:f=(AA-b)/AA
p.d=Sqr(X*X+Y*Y):r.d=p+Z*Z
u.d=ATan((Z/p)*((1-f)+(e2*AA/r)))
Phi=ATan((Z*(1-f)+e2*AA*Sin(u)*Sin(u)*Sin(u))/((1-f)*(p-e2*AA*Cos(u)*Cos(u)*Cos(u))))
Phi=Phi*Rad2Deg:LatWgs=Phi
Lambda.d=ATan(Y/X)
If X < 0:Lambda=Lambda+#PI:EndIf
Lambda=Lambda*Rad2Deg:LonWgs=Lambda
EndProcedure
Procedure NTF_Lamb(GLat.d,GLon.d, TypLamb$)
LIsom.d=LatIsom(GLat,eNTF)
n.d=0:c.d=0:Xs.d=0:Ys.d=0
Select TypLamb$
Case "L1"
n=0.7604059656:c=11603796.98:Xs=600000:Ys=5657616.674
Case "L2"
n=0.7289686274:c=11745793.39:Xs=600000:Ys=6199695.768
Case "L3"
n=0.6959127966:c=11947992.52:Xs=600000:Ys=6791905.085
Case "L4"
n=0.6712679322:c=12136281.99:Xs=234.358:Ys=7239161.542
Case "L2E"
n=0.7289686274:c=11745793.39:Xs=600000:Ys=8199695.768
EndSelect
X.d=Xs+c*Exp(-n*LIsom)*Sin(n*GLon*Deg2Rad)
X=Round(X,#PB_Round_Nearest):XProj=X
Y.d=Ys-c*Exp(-n*LIsom)*Cos(n*GLon*Deg2Rad)
Y=Round(Y,#PB_Round_Nearest):YProj=Y
EndProcedure
Procedure Lamb_NTF(XLambert.d,YLambert.d,TypLamb$)
n.d=0:c.d=0:Xs.d=0:Ys.d=0
Select TypLamb$
Case "L1"
n=0.7604059656:c=11603796.98:Xs=600000:Ys=5657616.674
Case "L2"
n=0.7289686274:c=11745793.39:Xs=600000:Ys=6199695.768
Case "L3"
n=0.6959127966:c=11947992.52:Xs=600000:Ys=6791905.085
Case "L4"
n=0.6712679322:c=12136281.99:Xs=234.358:Ys=7239161.542
Case "L2E"
n=0.7289686274:c=11745793.39:Xs=600000:Ys=8199695.768
EndSelect
R.d=Sqr((XLambert-Xs)*(XLambert-Xs)+(YLambert-Ys)*(YLambert-Ys))
g.d=ATan((XLambert-Xs)/(Ys-YLambert))
LX.d= g/n
LX=LX*Rad2Deg:LonNtf=LX
L.d=-1/n*Log(Abs(R/c))
LY.d=IsomLat(L, eNTF):LatNtf=LY
EndProcedure
Procedure WGS_L93(GLat.d,GLon.d)
ra.d=aRGF:;demi grand axe de l'ellipsoïde (m)
re.d=eRGF:;première excentricité de l'ellipsoïde
;paramètres de projection
rlc.d=3*deg2rad:;Méridien central:Lambda0=3° Est de Greenwich pour Lambert93
rPhi0.d=46.5*Deg2Rad:;latitude Origine pour Lambert93
rPhi1.d=44*Deg2Rad:;1er parallèle automécoïque
rPhi2.d=49*Deg2Rad:;2ème parallèle automécoïque
;coordonnées à l'origine
rx0.d=700000:rRy0.d=6600000
; coordonnées du point à traduire
rPhi.d=GLat*Deg2Rad:rl.d=GLon*Deg2Rad
; calcul des grandes normales
rgN1.d=ra/Sqr(1-re*re*Sin(rPhi1)*Sin(rPhi1))
rgN2.d=ra/Sqr(1-re*re*Sin(rPhi2)*Sin(rPhi2))
;calcul des latitudes isométriques
rRgl1.d=Log(Tan(#PI/4+rPhi1/2)*Exp(Log((1-re*Sin(rPhi1))/(1+re*Sin(rPhi1)))*re/2))
rRgl2.d=Log(Tan(#PI/4+rPhi2/2)*Exp(Log((1-re*Sin(rPhi2))/(1+re*Sin(rPhi2)))*re/2))
rRgl0.d=Log(Tan(#PI/4+rPhi0/2)*Exp(Log((1-re*Sin(rPhi0))/(1+re*Sin(rPhi0)))*re/2))
rRgl.d=Log(Tan(#PI/4+rPhi/2)*Exp(Log((1-re*Sin(rPhi))/(1+re*Sin(rPhi)))*re/2))
; calcul de l'exposant de la projection
rn.d=(Log((rgN2*Cos(rPhi2))/(rgN1*Cos(rPhi1))))/(rRgl1-rRgl2)
; calcul de la constante de projection
rcc.d=((rgN1*Cos(rPhi1))/rn)*Exp(rn*rRgl1)
; calcul des coordonnées Lambert93
X93.d=Round(rx0+rcc*Exp(-1*rn*rRgl)*Sin(rn*(rl-rlc)),#PB_Round_Nearest):XProj=X93
rys.d=rRy0+rcc*Exp(-1*rn*rRgl0)
Y93.d=Round(rys-rcc*Exp(-1*rn*rRgl)*Cos(rn*(rl-rlc)),#PB_Round_Nearest):YProj=Y93
EndProcedure
Procedure L93_WGS(XL93.d,YL93.d)
;Méridien central:Lambda0=3° Est de Greenwich pour Lambert93
n.d=0.725607765:c.d=11754255.426:Xs.d=700000:Ys.d=12655612.05
R.d=Sqr((XL93-Xs)*(XL93-Xs)+(YL93-Ys)*(YL93-Ys))
L.d=-1/n*Log(Abs(R/c)):LatWgs=IsomLat(L, eWGS)
g.d=ATan((XL93-Xs)/(Ys-YL93))
LX.d=g/n:LX=LX*Rad2Deg+3:LonWgs=LX
EndProcedure
Procedure$ WGS_DFCI(LatWGS.d, LonWGS.d)
;Calcul des coordonnées Cartographiques DFCI à partir des coordonnées Géographiques WGS84
;Lambert 2 étendu !
;1ERE LETTRE
NTF_Lamb(LatNtf,LonNtf,"L2E")
XLamb.d=XProj:YLamb.d=YProj-1500000
a=Int(XLamb/100000):XLamb=XLamb-a*100000:If a>7:a=a+2:EndIf
x$=Chr(a+65):res$=x$
;2EME LETTRE
a=Int((YLamb)/100000):YLamb=YLamb-a*100000:If a>7:a=a+2:EndIf
x$=Chr(a+65):res$=res$+x$
;3EME LETTRE (CHIFFRE)
a=Int(XLamb/20000):XLamb=XLamb-a*20000:res$=res$+Str(a*2)
;4EME LETTRE (CHIFFRE)
a=Int(YLamb/20000):YLamb=YLamb-a*20000:res$=res$+Str(a*2)
;5EME LETTRE
a=Int(XLamb/2000):XLamb=XLamb-a*2000:If a>7:a=a+2:EndIf
x$=Chr(a+65):res$=res$+x$
;6EME LETTRE (CHIFFRE)
a=Int(YLamb/2000):YLamb=YLamb-a*2000:res$=res$+a
If XLamb>500 And XLamb<1501 And YLamb>500 And YLamb<1501
res$=res$+".5"
Else
If XLamb<1000
If YLamb<1000:res$=res$+".4":Else:res$=res$+".1" :EndIf
Else
If YLamb<1000:res$=res$+".3":Else:res$=res$+".2":EndIf
EndIf
EndIf
ProcedureReturn(res$)
EndProcedure
Procedure DFCI_WGS(DFCI$)
YLamb.d=1500000 : XLamb.d=0
Carr100$=UCase(Left(DFCI$,2)): Carr20$=UCase(Mid(DFCI$,3,2)) : Carr2$=UCase(Mid(DFCI$,5,2)) : Zon5$=UCase(Mid(DFCI$,7,2))
If Len(Carr100$)=2
x$=Left(Carr100$,1) : a=FindString("ABCDEFGHKLMN",x$)-1 : XLamb=XLamb+a*100000
x$=Right(Carr100$,1) : a=FindString("ABCDEFGHKLMN",x$)-1 : YLamb=YLamb+a*100000
If Len(Carr20$)=2
x$=Left(Carr20$,1) : a=Val(x$): XLamb=XLamb+a*10000
x$=Right(Carr20$,1) : a=Val(x$) : YLamb=YLamb+a*10000
If Len(Carr2$)=2
x$=Left(Carr2$,1) : a=FindString("ABCDEFGHKLMN",x$)-1 : XLamb=XLamb+a*2000
x$=Right(Carr2$,1) : a=Val(x$) : YLamb=YLamb+a*2000
If Len(Zon5$)=2
x$=Right(Zon5$,1) : a=Val(x$)
Select a
Case 1
XLamb=XLamb+500 : YLamb=YLamb+1500
Case 2
XLamb=XLamb+1500 : YLamb=YLamb+1500
Case 3
XLamb=XLamb+500 : YLamb=YLamb+500
Case 4
XLamb=XLamb+500 : YLamb=YLamb+500
Default
XLamb=XLamb+1000 : YLamb=YLamb+1000
EndSelect
Else
XLamb=XLamb+1000 : YLamb=YLamb+1000
EndIf
Else
XLamb=XLamb+5000 : YLamb=YLamb+5000
EndIf
Else
XLamb=XLamb+50000 : YLamb=YLamb+50000
EndIf
EndIf
Lamb_NTF(XLamb,YLamb,"L2E")
NTF_WGS(LatNtf,LonNtf)
EndProcedure
Procedure WGS_UTM(LatWGS.d, LonWGS.d)
GabZon$="AACCCDDEEFFGGHHJJKKLLMMNNPPQQRRSSTTUUVVWWXXXZZ"
PrExc2.d=eWGS*eWGS:PrExc4.d=PrExc2*PrExc2:PrExc6.d=PrExc4*PrExc4
Lambda.d=LonWGS*Deg2Rad:Phi.d=LatWGS*Deg2Rad
SinPhi.d=Sin(Phi):S2.d=SinPhi*SinPhi:CosPhi.d=Cos(Phi):C2.d=CosPhi*CosPhi:TanPhi.d=Tan(Phi):T2.d=TanPhi*TanPhi
Lambda0.d=Int((LonWGS+180)/6)+1:NumFus+Lambda0:;Fuseau UTM
LettreZon$=Mid(GabZon$,Int((LatWGS+92)/4)+1,1):;Zone UTM
Lambda0=(6*Lambda0-183)*Deg2Rad:;Méridien de Référence de la Projection
A.d=(Lambda-Lambda0)*CosPhi
A2.d=A*A:A3.d=A*A2:A4.d=A*A3:A5.d=A*A4:A6.d=A*A5
k0.d=0.9996
k1.d=1-PrExc2/4-3*PrExc4/64-5*PrExc6/256
k2.d=3*PrExc2/8+3*PrExc4/32+45*PrExc6/1024
k3.d=15*PrExc4/256+45*PrExc6/1024
k4.d=35*PrExc6/3072
C.d=PrExc2/(1-PrExc2)*C2
sPhi.d=k1*Phi-k2*Sin(2*Phi)+k3*Sin(4*Phi)-k4*Sin(6*Phi)
vPhi.d=1/Sqr(1-PrExc2*S2)
tPhi.d=TanPhi*(A2/2+(5-T2+9*C+4*C*C)*A4/24+(61-58*T2+T2*T2)*A6/720)
Nz.d=0:If LatWGS<0:Nz=10000000:EndIf
N.d=Nz+k0*aWGS*(sPhi+vPhi*tPhi)
N=Round(N,#PB_Round_Nearest):YProj=N
E.d=500000+k0*aWGS*vPhi*(A+(1-T2+C)*A3/6+(5-18*T2+T2*T2)*A5/120)
E=Round(E,#PB_Round_Nearest):XProj=E
EndProcedure
Procedure UTM_WGS(UtmX.d, UtmY.d, UtmFuseau, UtmZone$)
;Calcul des coordonnées géographiques WGS84 en degrés décimaux à partir des coordonnées UTM
;d'après Algorithme de Steve Dutch-University of GreenBay, Wisconsin
;Longitudes Positives vers l'Est, Négatives vers l'Ouest
;Latitudes Positives vers le Nord, Négatives vers le Sud
GabZon$="AACCCDDEEFFGGHHJJKKLLMMNNPPQQRRSSTTUUVVWWXXXZZ"
e2.d=eWGS*eWGS:e4.d=e2*e2:e6.d=e4*e2
Hemisphere$="S":If UCase(UtmZone$)>"M":Hemisphere$="N":EndIf
East.d=UtmX:LatD.d=UtmY:If Hemisphere$="S":LatD=10000000-UtmY:EndIf
MCFuseau=6*UtmFuseau-183:;Méridien Central du Fuseau
e1sq.d=e2/(1-e2): k0.d=0.9996
arc.d=LatD/k0
mu.d=arc/(aWGS*(1-e2/4-3*e4/64-5*e6/256))
e1.d=(1-Sqr(1-e2))/(1+Sqr(1-e2))
ca.d=3*e1/2-27*Pow(e1,3)/32
cb.d=21*e1*e1/16-55*Pow(e1,4)/32
cc.d=151*Pow(e1,3)/96
cd.d=1097*Pow(e1,4)/512
phi1.d=mu+ca*Sin(2*mu)+cb*Sin(4*mu)+cc*Sin(6*mu)+cd*Sin(8*mu)
q0.d=e1sq*Cos(phi1)*Cos(phi1):t0.d=Tan(phi1)*Tan(phi1)
n0.d=aWGS/Sqr(1-Pow(eWGS*Sin(phi1),2))
r0.d=aWGS*(1-eWGS*eWGS)/Pow(1-Pow(eWGS*Sin(phi1),2),3/2)
dd0.d=(500000-East)/(n0*k0)
fact1.d=dd0
fact2.d=(1+2*t0+q0)*Pow(dd0,3)/6
fact3.d=(5-2*q0+28*t0-3*q0*q0+8*e1sq+24*t0*t0)*Pow(dd0,5)/120
LonWgs=((MCFuseau*Deg2Rad)-(fact1-fact2+fact3)/Cos(phi1))*Rad2Deg
fact1.d=n0*Tan(phi1)/r0:fact2=dd0*dd0/2
fact3=(5+3*t0+10*q0-4*q0*q0-9*e1sq)*Pow(dd0,4)/24
fact4=(61+90*t0+298*q0+45*t0*t0-252*e1sq-3*q0*q0)*Pow(dd0,6)/720
If Hemisphere$="N"
LatWgs=(phi1-fact1*(fact2+fact3+fact4))*Rad2Deg
Else
LatWgs=(-(phi1-fact1*(fact2+fact3+fact4)))*Rad2Deg
EndIf
EndProcedure
Procedure QuelleCible(CodCible$)
OutSys$=""
For i=1 To 14
If CodCible$=GeoSys$(i):OutSys$=CodCible$:EndIf
Next i
EndProcedure
Procedure QuelleSource(NPrm)
NomPoint$=LCmd$(NPrm-1):;Nom du Point
P2$=LCmd$(NPrm):;Système-Source
InSys$="":;Quel est le Système-Source ?
For i=0 To 12
If P2$=GeoSys$(i):InSys$=P2$:EndIf
Next i
If InSys$=""
Result$="Système-Source Inconnu : "+P2$
Result$=Result$+RC$+"Tapez GeKo pour l'aide"
Else
;Coordonnée(s)-Source
NPrm=2:NomWeb$=""
Select InSys$
Case "WGS"
NomWeb$="<br>Source : "+InSys$
Coord$=LCmd$(NPrm+1):LatWgs=ValD(Coord$)
NomWeb$=NomWeb$+"<br>Lat : "+Coord$
Coord$=LCmd$(NPrm+2):LonWgs=ValD(Coord$)
NomWeb$=NomWeb$+"<br>lon : "+Coord$
WGS_NTF(LatWgs,LonWgs)
Case "WGSD"
NomWeb$="<br>Source : "+InSys$
Coord$=LCmd$(NPrm+1):LatWgs=T2N(Coord$)
NomWeb$=NomWeb$+"<br>Lat : "+T2D(Coord$)
Coord$=LCmd$(NPrm+2):LonWgs=T2N(Coord$)
NomWeb$=NomWeb$+"<br>lon : "+T2D(Coord$)
WGS_NTF(LatWgs,LonWgs)
Case "NTF"
NomWeb$="<br>Source : "+InSys$
Coord$=LCmd$(NPrm+1):LatNtf=ValD(Coord$)
NomWeb$=NomWeb$+"<br>Lat : "+Coord$
Coord$=LCmd$(NPrm+2):LonNtf=ValD(Coord$)
NomWeb$=NomWeb$+"<br>lon : "+Coord$
NTF_WGS(LatNtf,LonNtf)
Case "NTFD"
NomWeb$="<br>Source : "+InSys$
Coord$=LCmd$(NPrm+1):LatNtf=T2N(Coord$)
NomWeb$=NomWeb$+"<br>Lat : "+T2D(Coord$)
Coord$=LCmd$(NPrm+2):LonNtf=T2N(Coord$)
NomWeb$=NomWeb$+"<br>lon : "+T2D(Coord$)
NTF_WGS(LatNtf,LonNtf)
Case "L1","L2","L3","L4","L2E"
NomWeb$="<br>Source : "+InSys$
Coord$=LCmd$(NPrm+1):XProj=ValD(Coord$)
NomWeb$=NomWeb$+"<br>X - Easting : "+Coord$
Coord$=LCmd$(NPrm+2):YProj=ValD(Coord$)
NomWeb$=NomWeb$+"<br>Y - Northing : "+Coord$
Lamb_NTF(XProj,YProj,InSys$)
NTF_WGS(LatNtf,LonNtf)
Case "L93"
NomWeb$="<br>Source : "+InSys$
Coord$=LCmd$(NPrm+1):XProj=ValD(Coord$)
NomWeb$=NomWeb$+"<br>X - Easting : "+Coord$
Coord$=LCmd$(NPrm+2):YProj=ValD(Coord$)
NomWeb$=NomWeb$+"<br>Y - Northing : "+Coord$
L93_WGS(XProj,YProj):WGS_NTF(LatWgs,LonWgs)
Case "DFCI"
NomWeb$="<br>Source : "+InSys$
Coord$=LCmd$(NPrm+1):DFCI_WGS(Coord$)
NomWeb$=NomWeb$+"<br>Carreau : "+Coord$
Case "UTM"
NomWeb$="<br>Source : "+InSys$
Coord$=LCmd$(NPrm+1):XProj=ValD(Coord$)
NomWeb$=NomWeb$+"<br>X - Easting : "+Coord$
Coord$=LCmd$(NPrm+2):YProj=ValD(Coord$)
NomWeb$=NomWeb$+"<br>Y - Northing : "+Coord$
Coord$=LCmd$(NPrm+3):NumFus=Val(Coord$)
NomWeb$=NomWeb$+"<br>Fuseau : "+Coord$
Coord$=LCmd$(NPrm+4):LettreZon$=Coord$
NomWeb$=NomWeb$+" - Zone : "+Coord$
UTM_WGS(XProj,YProj,NumFus,LettreZon$):WGS_NTF(LatWgs,LonWgs)
EndSelect
EndIf
EndProcedure
Procedure$ GetRGC(Lat.d,Lon.d)
InitNetwork()
x$="":NomFich$="GeKoAdr.txt"
ReqAdr$="https://maps.googleapis.com/maps/api/geocode/xml?latlng="+StrD(Lat)+","+StrD(Lon)
If ReceiveHTTPFile(ReqAdr$, NomFich$)
If ReadFile(2,NomFich$)
While Eof(2)=0
x$=ReadString(2):i=FindString(x$,"<formatted_address>",#PB_String_NoCase)
If i>0
x$=Mid(x$,3)
x$=RemoveString(x$,"<formatted_address>",#PB_String_NoCase)
x$=RemoveString(x$,"</formatted_address>",#PB_String_NoCase)
Break
EndIf
Wend
CloseFile(2)
DeleteFile(NomFich$)
EndIf
EndIf
ProcedureReturn(x$)
EndProcedure
Procedure GetResConv()
If OutSys$="WGS" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";WGS;"+StrD(LatWgs,6)+";"+StrD(LonWgs,6)
EndIf
If OutSys$="WGSD" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";WGSD;"+N2T(LatWgs)+";"+N2T(LonWgs)
EndIf
If OutSys$="NTF" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";NTF;"+StrD(LatNtf,6)+";"+StrD(LonNtf,6)
EndIf
If OutSys$="NTFD" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";NTFD;"+N2T(LatNtf)+";"+N2T(LonNtf)
EndIf
If OutSys$="L2E" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
NTF_Lamb(LatNtf,LonNtf,"L2E")
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";L2E;"+XProj+";"+YProj;
EndIf
If OutSys$="L1" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
NTF_Lamb(LatNtf,LonNtf,"L1")
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";L1;"+XProj+";"+YProj
EndIf
If OutSys$="L2" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
NTF_Lamb(LatNtf,LonNtf,"L2")
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";L2;"+XProj+";"+YProj
EndIf
If OutSys$="L3" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
NTF_Lamb(LatNtf,LonNtf,"L3")
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";L3;"+XProj+";"+YProj
EndIf
If OutSys$="L4" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
NTF_Lamb(LatNtf,LonNtf,"L4")
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";L4;"+XProj+";"+YProj
EndIf
If OutSys$="L93" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
WGS_L93(LatWgs,LonWgs)
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";L93;"+XProj+";"+YProj
EndIf
If OutSys$="DFCI" Or OutSys$="*"
If Result$<>"":Result$=Result$+RC$:EndIf
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";DFCI;"+WGS_DFCI(LatWgs,LonWgs)
EndIf
If OutSys$="UTM" Or OutSys$="*"
WGS_UTM(LatWgs,LonWgs)
If Result$<>"":Result$=Result$+RC$:EndIf
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";UTM;"+NumFus+";"+LettreZon$+";"+XProj+";"+YProj
EndIf
If OutSys$="RGC"
If Result$<>"":Result$=Result$+RC$:EndIf
Print(NomPoint$)
x$=GetRGC(LatWgs,LonWgs)
If x$<>""
PrintN(TAB$+"OK")
Result$=Result$+"GeKo"+";"+Chr(34)+NomPoint$+Chr(34)+";"+Chr(34)+x$+Chr(34)
Else
PrintN(TAB$+"KO")
Result$=Result$+"<RGC>"+TAB$+"<<< geocodage inverse impossible avec GM >>>"
EndIf
; Attente 1 seconde et demie avant une nouvelle requête (voir Aide - serveur GoogleMaps)
Delay(1500)
EndIf
If OutSys$="WEB"
If ResWeb$<>"":ResWeb$=ResWeb$+RC$:EndIf
ResWeb$=ResWeb$+TAB$+TAB$+TAB$+TAB$+"<tr>"
ResWeb$=ResWeb$+RC$+TAB$+TAB$+TAB$+TAB$+TAB$+"<td bgcolor="+Chr(34)+"#FFFFB0"+Chr(34)+"><b>"+NomPoint$+"</b>"+NomWeb$+"</td><td>"+StrD(LatWgs)+"</td><td>"+StrD(LonWgs)+"</td>"
ResWeb$=ResWeb$+RC$+TAB$+TAB$+TAB$+TAB$+TAB$+"<td><a target="+Chr(34)+"blank"+Chr(34)+" href="+Chr(34)+"http://www.openstreetmap.org/?lat="+StrD(LatWgs)+"&lon="+StrD(LonWgs)+"&zoom=17"+Chr(34)+">Carte OSM<SUP>®</SUP></a></td>"
ResWeb$=ResWeb$+RC$+TAB$+TAB$+TAB$+TAB$+TAB$+"<td><a target="+Chr(34)+"blank"+Chr(34)+" href="+Chr(34)+"http://www.geoportail.gouv.fr/embed/visu.html?c="+StrD(LonWgs)+","+StrD(LatWgs)+"&z=0.00002&l=GEOGRAPHICALGRIDSYSTEMS.MAPS.3D::GEOPORTAIL:OGC:WMTS==aggregate%281%29"+Chr(34)+">Carte IGN<SUP>®</SUP></a></td>"
ResWeb$=ResWeb$+RC$+TAB$+TAB$+TAB$+TAB$+TAB$+"<td><a target="+Chr(34)+"blank"+Chr(34)+" href="+Chr(34)+"http://maps.google.fr/maps?q="+StrD(LatWgs)+","+StrD(LonWgs)+Chr(34)+">Carte GM<SUP>®</SUP></a></td>"
ResWeb$=ResWeb$+RC$+TAB$+TAB$+TAB$+TAB$+"</tr>"
EndIf
EndProcedure
OpenConsole()
InitGeko()
NomPoint$="":OptHlp=0:P1$="":P2$=""
If NbParam=0
;Pas de paramètre - appel à l'aide
OptHlp=1
Else
P1$=LCmd$(0)
If P1$="/H" Or P1$="-H" Or P1$="/?" Or P1$="-?"
;1er paramètre - appel à l'aide
OptHlp=1
Else
;1er paramètre - CSV ?
If P1$="CSV"
;Mode Conversion par lot (d'un fichier CSV)
P2$=LCmd$(1)
;2ème paramètre - Nom et chemin du Fichier-Source
If FileSize(P2$)>0
;Le fichier existe, on le traite ...
l$=""
If ReadFile(0,P2$)
While Eof(0)=0
l$=ReadString(0)
If l$<>""
For i=0 To 7
LCmd$(i)=RemoveString(StringField(l$,i+1,";"),Chr(34))
If i<>1:LCmd$(i)=ReplaceString(LCmd$(i),",","."):LCmd$(i)=UCase(LCmd$(i)):EndIf
Next i
P1$=LCmd$(0)
QuelleCible(P1$)
QuelleSource(2)
GetResConv()
EndIf
Wend
CloseFile(0)
EndIf
Else
;Le fichier n'existe pas ou est vide, on conseille de lire la doc !
PrintN("Fichier CSV vide ou inexistant : "+P2$)
OptHlp=1
EndIf
Else
;Mode Immédiat Console
;1er paramètre - Système-Cible - Quel est le Système-Cible ?
If P1$="*":OutSys$="*":Else:QuelleCible(P1$):EndIf
If OutSys$=""
Result$="Système-Cible Inconnu : "+P1$
Result$=Result$+RC$+"Tapez GeKo pour l'aide"
Else
QuelleSource(2)
GetResConv()
EndIf
EndIf
EndIf
EndIf
If OptHlp=1
Aide()
Else
If OutSys$="WEB"
CurDir$=GetCurrentDirectory():FilRes$=CurDir$+"GeKoRes.htm"
CreateFile(1,FilRes$)
WriteString(1,WEB1$+RC$+ResWeb$+RC$+WEB2$)
CloseFile(1)
If Result$<>"":Result$=Result$+RC$:EndIf
Result$=Result$+RC$+"Ouvrez le fichier "+FilRes$+" avec votre navigateur WEB"+RC$
Select Sys$
Case "W" : RunProgram(FilRes$)
Case "X" : RunProgram("xdg-open",FilRes$,""):PrintN(Result$):PrintN("Appuyez sur [Entree] pour quitter"):Input()
Default :PrintN(Result$)
EndSelect
Else
PrintN(Result$)
PrintN("Appuyez sur [Entree] pour quitter"):Input()
EndIf
EndIf
CloseConsole()
