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.
- DFCI_LatWGS et DFCI_LonWGS transforment les coordonnées du carré DFCI en coordonées WGS84 quand WGS_DFCI effectue l'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.
--- 3 avril 2017 ----
J'ai ajouté la conversion UTM vers MGRS et réciproque à la demande d'un intervenant sur le forum. J'ai été un peu présomptueux dans ma réponse précédente, ça m'a pris plus de temps que prévu ...
J'en ai profité pour corriger quelques erreurs mineures (redondances inutiles) et j'ai pris en compte les particularités de la Norvège et du Svalbard dans le carroyage UTM (et donc MGRS).
La prochaine étape, à la demande de Papyzède, sera de développer les conversions avec les coniques conformes zonées utilisées en France métropolitaine.
à bientôt,
PK1157
PS : le fichier d'exemple et le code ont été mis à jour.
Code : Tout sélectionner
REM*****BASIC*****
FUNCTION DMSDec(Param AS STRING) AS DOUBLE
REM Conversion Degrés, Minutes, Secondes en Degrés Décimaux
REM Vérifié 23 mars 2017
DIM ok AS BOOLEAN
DIM x AS STRING
DIM neg AS BOOLEAN
DIM d,m,s,z AS DOUBLE
DIM DMS(2)
DIM i as INTEGER
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,CHR(34))>0 THEN
ok=1: x=REPLACE(x,CHR(34),"\")
ENDIF
IF INSTR(x,CHR(39))>0 THEN
ok=1: x=REPLACE(x,CHR(39),"\")
ENDIF
IF INSTR(x,CHR(8217))>0 THEN
ok=1: x=REPLACE(x,CHR(8217),"\")
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
IF ABS(d)>180 OR ABS(m)>59 OR INT(ABS(s))>59 THEN z=-9999
DMSDec=z
END FUNCTION
FUNCTION DecDMS(P1 AS DOUBLE, OPTIONAL P2 AS INTEGER) AS STRING
REM Conversion Degrés Décimaux en Degrés, Minutes, Secondes
REM P1: réel à convertir en chaîne
REM P2: 0 pour Latitude, valeur non nulle pour Longitude
REM Vérifié 23 mars 2017
DIM x AS STRING
DIM neg AS BOOLEAN
DIM d,m,s,z AS DOUBLE
DIM MaxVal as DOUBLE
z=P1 : MaxVal=180
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((INT(z*1000+0.5)/1000)*10)/10
x=d & "°" & m & "'" & s & CHR(34)
IF ISMISSING(P2) THEN
IF neg THEN x="-" & x
ELSE
IF P2=0 THEN
MaxVal=90
IF neg THEN x=x & " S" ELSE x=x &" N"
ELSE
IF neg THEN x=x & " W" ELSE x=x &" E"
ENDIF
ENDIF
IF ABS(z)>MaxVal THEN DecDMS="Valeur hors bornes !" ELSE DecDMS=x
END FUNCTION
FUNCTION WGSLatNTF(GLat AS DOUBLE, GLon AS DOUBLE) AS DOUBLE
REM Calcul de la Latitude NTF à partir des coordonnées géographiques WGS84 en Degrés décimaux
REM GLat, GLon en Degrés décimaux
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 Degrés <---> Radians
DIM Rad2Deg, Deg2Rad AS DOUBLE
Rad2Deg=180/PI(): Deg2Rad=PI()/180
REM ----------------------------------
REM h = hauteur sur l'ellipsoïde - mis à zéro
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 = Première Excentricité 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 = Première Excentricité de l'ellipsoïde WGS84
eWGS = SQR(1 - ((bWGS / aWGS)*(bWGS / aWGS)))
AA = aWGS: b=bWGS : e2 = eWGS * eWGS: f = (AA - b) / AA
v = AA / SQR(1 - e2 * SIN(Phi) * SIN(Phi))
REM Coordonnée Géocentriques 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 référentiel 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 coordonnées géographiques WGS84 en Degrés décimaux
REM GLat, GLon en Degrés décimaux
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 Degrés <---> Radians
DIM Rad2Deg, Deg2Rad AS DOUBLE
Rad2Deg=180/PI(): Deg2Rad=PI()/180
REM ----------------------------------
REM h = hauteur sur l'ellipsoïde - mis à zéro
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 = Première Excentricité 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 Coordonnée Géocentriques 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 référentiel 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 coordonnées géographiques NTF en Degrés décimaux
REM GLat, GLon en Degrés décimaux
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 à zéro
DIM MeridParis AS DOUBLE
MeridParis=2.33722916666667
REM 2°20'14,025" E de Greenwich
REM ----------------------------------
REM Conversions Degrés <---> 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 = Première Excentricité 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 = Première Excentricité 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 Coordonnée Géocentriques 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 référentiel 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 coordonnées géographiques NTF en Degrés décimaux
REM GLat, GLon en Degrés décimaux
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 Degrés <---> Radians
DIM Rad2Deg, Deg2Rad AS DOUBLE
Rad2Deg=180/PI(): Deg2Rad=PI()/180
REM ----------------------------------
REM h = hauteur sur l'ellipsoïde - mis à zéro
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 = Première Excentricité 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 Coordonnée Géocentriques 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 référentiel 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 Géographique en Degrés Décimaux
REM PremExEllips : Première excentricité de l'ellipsoïde (par Défaut, 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.081819191042815
REM ----------------------------------
REM Conversions Degrés <---> 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 : Première excentricité de l'ellipsoïde (par Défaut, IAG/GRS 1980 pour WGS84)
DIM Clarke1880 AS DOUBLE
DIM IagGrs80 AS DOUBLE
Clarke1880=0.08248325676342
IagGrs80=0.081819191042815
DIM s0, s1, d, EL, rPremEx, rTolConv AS DOUBLE
REM ----------------------------------
REM Conversions Degrés <---> 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 Degrés <---> 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 étendu par défaut
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 Degrés <---> 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 étendu par défaut
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 Degrés <---> 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 étendu par défaut
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 Degrés <---> 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 étendu par défaut
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 Degrés <---> Radians
DIM Rad2Deg, Deg2Rad AS DOUBLE
Rad2Deg=180/PI(): Deg2Rad=PI()/180
REM ----------------------------------
REM système WGS84
IagGrs80=0.081819191042815
rAWGS=6378137
ra=rAWGS : REM demi grand axe de l'ellipsoïde (m)
re=IagGrs80: REM première excentricité de l'ellipsoïde
REM paramètres de projection
rlc=3*deg2rad: REM Méridien central : Lambda0 = 3° Est de Greenwich pour Lambert93
rPhi0=46.5*Deg2Rad: REM latitude Origine pour Lambert93
rPhi1=44*Deg2Rad: REM 1er parallèle automécoïque
rPhi2=49*Deg2Rad: REM 2ème parallèle automécoïque
REM coordonnées à l'origine
rx0=700000: rRy0=6600000
REM coordonnées 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 isométriques
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 Degrés <---> Radians
DIM Rad2Deg, Deg2Rad AS DOUBLE
Rad2Deg=180/PI(): Deg2Rad=PI()/180
REM ----------------------------------
REM système WGS84
IagGrs80=0.081819191042815
rAWGS=6378137
ra=rAWGS : REM demi grand axe de l'ellipsoïde (m)
re=IagGrs80: REM première excentricité de l'ellipsoïde
REM paramètres de projection
rlc=3*deg2rad: REM Méridien central : Lambda0 = 3° Est de Greenwich pour Lambert93
rPhi0=46.5*Deg2Rad: REM latitude Origine pour Lambert93
rPhi1=44*Deg2Rad: REM 1er parallèle automécoïque
rPhi2=49*Deg2Rad: REM 2ème parallèle automécoïque
REM coordonnées à l'origine
rx0=700000: rRy0=6600000
REM coordonnées 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 isométriques
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 coordonnées 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 Degrés <---> Radians
DIM Rad2Deg, Deg2Rad AS DOUBLE
Rad2Deg=180/PI(): Deg2Rad=PI()/180
REM ----------------------------------
IagGrs80=0.081819191042815
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 Degrés <---> 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 Méridien 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="CDEFGHJKLMNPQRSTUVWXX"
NumFus=INT((LonWGS+180)/6)+1
IF LatWGS>=84 THEN
REM Arctique
IF NumFus<31 THEN LettreZon="Y" ELSE LettreZon="Z"
ELSE
IF LatWGS<=-80 THEN
REM Antarctique
IF NumFus<31 THEN LettreZon="A" ELSE LettreZon="B"
ELSE
REM Cas General
LettreZon=MID(GabZon,INT((LatWGS+80)/8)+1,1)
IF LatWGS>=56 AND LatWGS<64 AND LonWGS>=0 AND LonWGS<12 THEN
REM Exception Norvege
IF LonWGS>=3 THEN NumFus=32 ELSE NumFus=31
ENDIF
IF LatWGS>=72 AND LatWGS<84 AND LonWGS>=0 AND LonWGS<42 THEN
REM Exception Svalbard
NumFus=31+2*INT((LonWGS+3)/12)
ENDIF
ENDIF
ENDIF
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 Degrés <---> 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 Première Excentricité 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 Méridien de Référence de la Projection
REM ------------ Prise en compte particularites grille UTM 03/04/2017 --
IF LatWGS>=56 AND LatWGS<64 AND LonWGS>=0 AND LonWGS<12 THEN
REM Exception Norvege
IF LonWGS>=3 THEN Lambda0=7.5 ELSE Lambda0=1.5
ENDIF
IF LatWGS>=72 AND LatWGS<84 AND LonWGS>=0 AND LonWGS<42 THEN
REM Exception Svalbard
Lambda0=4.5
IF LonWGS>=9 THEN Lambda0=15
IF LonWGS>=21 THEN Lambda0=27
IF LonWGS>=33 THEN Lambda0=37.5
ENDIF
REM --------------------------------------------------------------------
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 Degrés <---> 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 Première Excentricité 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 Méridien de Référence de la Projection
REM ------------ Prise en compte particularites grille UTM 03/04/2017 --
IF LatWGS>=56 AND LatWGS<64 AND LonWGS>=0 AND LonWGS<12 THEN
REM Exception Norvege
IF LonWGS>=3 THEN Lambda0=7.5 ELSE Lambda0=1.5
ENDIF
IF LatWGS>=72 AND LatWGS<84 AND LonWGS>=0 AND LonWGS<42 THEN
REM Exception Svalbard
Lambda0=4.5
IF LonWGS>=9 THEN Lambda0=15
IF LonWGS>=21 THEN Lambda0=27
IF LonWGS>=33 THEN Lambda0=37.5
ENDIF
REM --------------------------------------------------------------------
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=10000000 ELSE Nz=0
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 géographique 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
'
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 Degrés <---> Radians
DIM Rad2Deg, Deg2Rad AS DOUBLE
Rad2Deg=180/PI(): Deg2Rad=PI()/180
REM ----------------------------------
IF UCase(UtmZone)>"M" THEN Hemisphere="N" ELSE Hemisphere="S"
REM système 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 = Première Excentricité 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 géographique 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
'
DIM LatD, East AS DOUBLE
DIM MCFuseau 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 AS DOUBLE
DIM Hemisphere AS STRING
REM ----------------------------------
REM Conversions Degrés <---> Radians
DIM Rad2Deg, Deg2Rad AS DOUBLE
Rad2Deg=180/PI(): Deg2Rad=PI()/180
REM ----------------------------------
IF UCase(UtmZone)>"M" THEN Hemisphere="N" ELSE Hemisphere="S"
REM système 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 = Première Excentricité 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
' Méridien Central du Fuseau
MCFuseau = 6 * UtmFuseau - 183
REM ------------ Prise en compte particularites grille UTM 03/04/2017 --
REM Exception Norvege
IF UtmZone="V" AND UtmFuseau=31 THEN MCFuseau=1.5
IF UtmZone="V" AND UtmFuseau=32 THEN MCFuseau=7.5
REM Exception Svalbard
IF UtmZone="X" AND UtmFuseau=31 THEN MCFuseau=4.5
IF UtmZone="X" AND UtmFuseau=33 THEN MCFuseau=15
IF UtmZone="X" AND UtmFuseau=35 THEN MCFuseau=27
IF UtmZone="X" AND UtmFuseau=37 THEN MCFuseau=37.5
REM --------------------------------------------------------------------
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
FUNCTION DFCI_LatWGS(DFCI AS STRING) AS DOUBLE
' Calcul de la Latitude geographique WGS84 en degres decimaux à partir des coordonnees DFCI
'
DIM Carr100, Carr20, Carr2, Zon5 AS STRING
DIM XLamb, YLamb, LatNTF, LonNTF AS DOUBLE
DIM x AS STRING
DIM a AS DOUBLE
REM ----------------------------------
Carr100=UCase(LEFT(DFCI,2)): Carr20=UCase(MID(DFCI,3,2)) : Carr2=UCase(MID(DFCI,5,2)) : Zon5=UCase(MID(DFCI,7,2))
YLamb=1500000 : XLamb=0
IF Len(Carr100)=2 THEN
x=left(Carr100,1) : a=InStr("ABCDEFGHKLMN",x)-1 : XLamb=XLamb+a*100000
x=right(Carr100,1) : a=InStr("ABCDEFGHKLMN",x)-1 : YLamb=YLamb+a*100000
IF Len(Carr20)=2 THEN
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 THEN
x=left(Carr2,1) : a=InStr("ABCDEFGHKLMN",x)-1 : XLamb=XLamb+a*2000
x=right(Carr2,1) : a=Val(x) : YLamb=YLamb+a*2000
IF Len(Zon5)=2 THEN
x=right(Zon5,1) : a=Val(x)
SELECT CASE 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
CASE ELSE
XLamb=XLamb+1000 : YLamb=YLamb+1000
END SELECT
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
LatNTF=LambLatNTF(XLamb,YLamb,0) : LonNTF=LambLonNTF(XLamb,YLamb,0)
DFCI_LatWGS=NTFLatWGS(LatNTF,LonNTF)
END FUNCTION
FUNCTION DFCI_LonWGS(DFCI AS STRING) AS DOUBLE
' Calcul de la Longitude geographique WGS84 en degres decimaux à partir des coordonnees DFCI
'
DIM Carr100, Carr20, Carr2, Zon5 AS STRING
DIM XLamb, YLamb, LatNTF, LonNTF AS DOUBLE
DIM x AS STRING
DIM a AS DOUBLE
REM ----------------------------------
Carr100=UCase(LEFT(DFCI,2)): Carr20=UCase(MID(DFCI,3,2)) : Carr2=UCase(MID(DFCI,5,2)) : Zon5=UCase(MID(DFCI,7,2))
YLamb=1500000 : XLamb=0
IF Len(Carr100)=2 THEN
x=left(Carr100,1) : a=InStr("ABCDEFGHKLMN",x)-1 : XLamb=XLamb+a*100000
x=right(Carr100,1) : a=InStr("ABCDEFGHKLMN",x)-1 : YLamb=YLamb+a*100000
IF Len(Carr20)=2 THEN
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 THEN
x=left(Carr2,1) : a=InStr("ABCDEFGHKLMN",x)-1 : XLamb=XLamb+a*2000
x=right(Carr2,1) : a=Val(x) : YLamb=YLamb+a*2000
IF Len(Zon5)=2 THEN
x=right(Zon5,1) : a=Val(x)
SELECT CASE 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
CASE ELSE
XLamb=XLamb+1000 : YLamb=YLamb+1000
END SELECT
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
LatNTF=LambLatNTF(XLamb,YLamb,0) : LonNTF=LambLonNTF(XLamb,YLamb,0)
DFCI_LonWGS=NTFLonWGS(LatNTF,LonNTF)
END FUNCTION
FUNCTION WGS_DFCI(LatWGS AS DOUBLE, LonWGS AS DOUBLE) AS STRING
' Calcul des coordonnées Cartographiques DFCI à partir des coordonnées Géographiques WGS84
'
DIM XLamb, YLamb, LatNTF, LonNTF AS DOUBLE
DIM x, res AS STRING
DIM a AS DOUBLE
REM ----------------------------------
LatNTF=WGSLatNTF(LatWGS,LonWGS) : LonNTF=WGSLonNTF(LatWGS,LonWGS)
REM Lambert 2 étendu !
REM 1ERE LETTRE
XLamb=NTFXLamb(LatNTF,LonNtf,0)
a=INT(XLamb/100000) : XLamb=XLamb-a*100000 : IF a >7 THEN a=a+2
x=Chr(a+65) : res=x
REM 2EME LETTRE
YLamb=NTFYLamb(LatNTF,LonNtf,0)-1500000
a=INT(YLamb/100000) : YLamb=YLamb-a*100000 : IF a >7 THEN a=a+2
x=Chr(a+65) : res=res & x
REM 3EME LETTRE (CHIFFRE)
a=INT(XLamb/20000) : XLamb=XLamb-a*20000 : res=res & a*2
REM 4EME LETTRE (CHIFFRE)
a=INT(YLamb/20000) : YLamb=YLamb-a*20000 : res=res & a*2
REM 5EME LETTRE
a=INT(XLamb/2000) : XLamb=XLamb-a*2000 : IF a >7 THEN a=a+2
x=Chr(a+65) : res=res & x
REM 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 THEN
res=res & ".5"
ELSE
IF XLamb<1000 THEN
IF YLamb<1000 THEN res=res & ".4" ELSE res=res & ".1"
ELSE
IF YLamb<1000 THEN res=res & ".3" ELSE res=res & ".2"
ENDIF
ENDIF
WGS_DFCI=res
END FUNCTION
FUNCTION GeoDeCode(GLat AS DOUBLE, GLon AS DOUBLE, OPTIONAL Serv AS INTEGER) AS STRING
REM Géocodage inverse : fournit l'adresse postale à partir des coordonnées WGS84
REM Le paramètre optionnel Serv permet de choisir entre Google(1) et OpenStreetMap(0 ou absent)
DIM T As Date
IF ISMISSING(Serv) THEN Serv=0
Svc = CreateUnoService( "com.sun.star.sheet.FunctionAccess" ) 'Crée un service pour utiliser les fonctions Calc
Y=CStr(GLat) : Virgule=InStr(Y,","): If Virgule>0 Then Mid( Y, Virgule, 1 ) = "."
X=CStr(GLon) : Virgule=InStr(X,","): If Virgule>0 Then Mid( X, Virgule, 1 ) = "."
IF Serv=0 THEN
ReqGMap="http://nominatim.openstreetmap.org/reverse?format=xml&lat=" & Y & "&lon=" & X
Balise="result"
ELSE
ReqGMap="https://maps.googleapis.com/maps/api/geocode/xml?latlng=" & Y & "," & X
Balise="formatted_address"
ENDIF
' ajoute une temporisation de 1 seconde pour empêcher que le serveur n'envoie un message de refus !
T = Timer + 1: Do Until Timer > T: DoEvents: Loop
' (les serveurs n'acceptent pas de demandes répétées trop fréquentes provenant de la même IP)
XML_String = svc.callFunction("WEBSERVICE",Array(ReqGMap))
DataDeb=instr(XML_String,"<" & Balise) : DataFin=instr(XML_String,"</" & Balise)
XML_String = Mid(XML_String,DataDeb,DataFin-DataDeb)
DataDeb=instr(XML_String,">") : XML_String = Mid(XML_String,DataDeb+1,65535)
GeoDeCode=XML_String
END FUNCTION
FUNCTION UTM_MGRS(UtmX AS DOUBLE, UtmY AS DOUBLE, UtmFuseau AS INTEGER, UtmZone AS STRING, OPTIONAL Prec AS INTEGER, OPTIONAL CarSep AS STRING) AS STRING
REM Transformation de coordonnées UTM au format militaire MGRS
REM Le parametre Prec regit la precision souhaitee (le nombre de chiffres : 1=10km 2=1km 3=100m 4=10m 5=1m) 5 par defaut
REM Le parametre CarSep correspond au separateur des groupes dans le resultat, <aucun> par defaut
DIM Sep, res AS STRING
DIM L1, L2, Offset AS LONG
DIM X, Y, Fact AS DOUBLE
IF ISMISSING(Prec) THEN
Prec=5
ELSE
IF Prec>5 THEN Prec=5
IF Prec<1 THEN Prec=1
ENDIF
Fact=0.5*10^(5-Prec)
IF ISMISSING(CarSep) THEN Sep=" " ELSE Sep=LEFT(CarSep,1)
res=FORMAT(UtmFuseau,"00") & Sep & UtmZone & Sep
L1=((UtmFuseau-1) MOD 3)*8+INT(UtmX/100000)
IF UtmFuseau MOD 2 = 0 THEN Offset=5 ELSE Offset=0
L2=((Offset+INT(UtmY/100000)) mod 20)+1
res=res & Mid("ABCDEFGHJKLMNPQRSTUVWXYZ",L1,1) & Mid("ABCDEFGHJKLMNPQRSTUV",L2,1) & Sep
X=UtmX-INT(UtmX/100000)*100000 : X=INT(X+Fact)
Y=UtmY-INT(UtmY/100000)*100000 : Y=INT(Y+Fact)
res = res & LEFT(FORMAT(X,"00000"),Prec) & Sep & LEFT(FORMAT(Y,"00000"),Prec)
UTM_MGRS=res
END FUNCTION
FUNCTION MGRS_FusUTM(Mgrs AS STRING) AS STRING
DIM l,i AS INTEGER
DIM c, VarIn AS STRING
VarIn="" : l=LEN(Mgrs)
FOR i=1 TO l
c=MID(Mgrs,i,1)
IF c>="A" AND c<="Z" THEN VarIn=VarIn & c
IF c>="0" AND c<="9" THEN VarIn=VarIn & c
NEXT i
MGRS_FusUTM=LEFT(VarIn,2)
END FUNCTION
FUNCTION MGRS_ZonUTM(Mgrs AS STRING) AS STRING
DIM l,i AS INTEGER
DIM c, VarIn AS STRING
VarIn="" : l=LEN(Mgrs)
FOR i=1 TO l
c=MID(Mgrs,i,1)
IF c>="A" AND c<="Z" THEN VarIn=VarIn & c
IF c>="0" AND c<="9" THEN VarIn=VarIn & c
NEXT i
MGRS_ZonUTM=MID(VarIn,3,1)
END FUNCTION
FUNCTION MGRS_UTMEasting(Mgrs AS STRING) AS DOUBLE
DIM l,i AS INTEGER
DIM c, VarIn, L1, XY AS STRING
DIM X AS DOUBLE
VarIn="" : l=LEN(Mgrs)
FOR i=1 TO l
c=MID(Mgrs,i,1)
IF c>="A" AND c<="Z" THEN VarIn=VarIn & c
IF c>="0" AND c<="9" THEN VarIn=VarIn & c
NEXT i
L1=MID(VarIn,4,1) : XY=MID(VarIn,6) : l=LEN(XY)/2
X=VAL(LEFT(XY,l))
OffX=Instr("ABCDEFGHJKLMNPQRSTUVWXYZ",L1)-1
OffX=(OffX MOD 8)+1
OffX=OffX*100000
MGRS_UTMEasting=OffX+X
END FUNCTION
FUNCTION MGRS_UTMNorthing(Mgrs AS STRING) AS DOUBLE
DIM l,i,j, RUtm AS INTEGER
DIM c, VarIn, L2,XY, ZonUtm AS STRING
DIM Y, FloorUtm AS DOUBLE
VarIn="" : l=LEN(Mgrs)
FOR i=1 TO l
c=MID(Mgrs,i,1)
IF c>="A" AND c<="Z" THEN VarIn=VarIn & c
IF c>="0" AND c<="9" THEN VarIn=VarIn & c
NEXT i
L2=MID(VarIn,5,1) : ZonUtm=MID(VarIn,3,1) : XY=MID(VarIn,6) : l=LEN(XY)/2
Y=VAL(RIGHT(XY,l))
OffY=Instr("ABCDEFGHJKLMNPQRSTUV",L2)-1
rem IF ZonUtm>"B" AND ZonUtm<"Y" THEN
i=VAL(LEFT(VarIn,2)) : j=INT(i/2) : IF i=2*j THEN OffY=OffY-5
rem ENDIF
OffY=OffY*100000 : OffY=OffY+Y
IF ZonUtm>="N" THEN
REM Hemisphere Nord
RUtm=Instr("ABCDEFGHJKLMNPQRSTUVWXYZ",ZonUtm)-13 : IF RUtm>10 THEN RUtm=10
FloorUtm=110500*8*RUtm
WHILE OffY<FloorUtm
OffY=OffY+2000000
WEND
ELSE
REM Hemisphere Sud
RUtm=13-Instr("ABCDEFGHJKLMNPQRSTUVWXYZ",ZonUtm)
FloorUtm=110500*(90-(RUtm*8))
WHILE OffY<FloorUtm
OffY=OffY+2000000
WEND
ENDIF
MGRS_UTMNorthing=OffY
END FUNCTION