1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
| public final static double[] WGS84toLambert2e(double d_long, double m_long, double s_long, char orientation_long, double d_lat, double m_lat, double s_lat , char orientation_lat)
{
double lambda_w, phi_w;
/**************************************************************************************************************/
/** 0) degres-minutes-secondes + orientation (d,m,s,o) -> radian **/
/**************************************************************************************************************/
// Pour la longitude
lambda_w = d_long + m_long/60 + s_long/3600;
if(orientation_long == 'W') lambda_w = -1 * lambda_w; // Le système de coordonnées géographiques utilisé est postif vers le Nord et vers l'Est
lambda_w = lambda_w*Math.PI/180 ;
// Pour la latitude
phi_w = d_lat + m_lat /60 + s_lat /3600;
if(orientation_lat == 'S') phi_w = -1 * phi_w; // Le système de coordonnées géographiques utilisé est postif vers le Nord et vers l'Est
phi_w = phi_w*Math.PI/180 ;
/**************************************************************************************************************/
/** 1) coordonnées géographiques WGS84 (phi_w,lambda_w) -> coordonnées cartésiennes WGS84 (X_w,Y_w,Z_w) **/
/**************************************************************************************************************/
// Nous avons utilisés 2 formules données par l'IGN dans un de leur document ainsi que deux constantes de
// l'ellipsoide de référence associé au système WGS84 (les deux demi-axes) :
double a_w = 6378137.0;// demi grand axe
double b_w = 6356752.314;// demi petit axe
// d'où
double e2_w = (a_w*a_w-b_w*b_w)/(a_w*a_w);
// et on a la grande normale de l'ellipsoide WGS84
double N = a_w/Math.sqrt(1 - e2_w * Math.pow(Math.sin(phi_w),2));
// ainsi on obtient
double X_w = N * Math.cos(phi_w) * Math.cos(lambda_w);
double Y_w = N * Math.cos(phi_w) * Math.sin(lambda_w);
double Z_w = N * (1-e2_w) * Math.sin(phi_w);
/**************************************************************************************************************/
/** 2) coordonnées cartésiennes WGS84 (X_w,Y_w,Z_w) -> coordonnées cartésiennes NTF (X_n,Y_n,Z_n) **/
/**************************************************************************************************************/
// Là il n'y a qu'une translation à effectuer :
// on a donc
double dX = 168.0;
double dY = 60.0;
double dZ = -320.0;
double X_n = X_w + dX;
double Y_n = Y_w + dY;
double Z_n = Z_w + dZ;
/**************************************************************************************************************/
/** 3) coordonnées cartésiennes NTF (X_n,Y_n,Z_n) -> coordonnées géographiques NTF (phi_n,lambda_n) **/
/**************************************************************************************************************/
// J'ai utilisé 1 formule donnée par l'IGN toujours dans le même document ainsi que deux constantes de l'ellipsoide
// de référence du système NTF, Clarke 1880 :
double a_n = 6378249.2;
double b_n = 6356515.0;
// d'où
double e2_n = (a_n*a_n-b_n*b_n)/(a_n*a_n);
// on définit une tolérance de convergence
double epsilon = Math.pow(10,-10);
// puis on amorce une boucle de calcul
double phi0 = Math.atan(Z_n/Math.sqrt(X_n*X_n+Y_n*Y_n)*(1-(a_n*e2_n)/(Math.sqrt(X_n*X_n+Y_n*Y_n+Z_n*Z_n))));
double phi1 = Math.atan((Z_n/Math.sqrt(X_n*X_n+Y_n*Y_n))/(1-(a_n*e2_n*Math.cos(phi0))/(Math.sqrt((X_n*X_n+Y_n*Y_n)*(1-e2_n*Math.pow(Math.sin(phi0),2))))));
while(!(Math.abs(phi1-phi0)<epsilon))
{
phi0 = phi1; phi1 = Math.atan((Z_n/Math.sqrt(X_n*X_n+Y_n*Y_n))/(1-(a_n*e2_n*Math.cos(phi0))/(Math.sqrt((X_n*X_n+Y_n*Y_n)*(1-e2_n*Math.pow(Math.sin(phi0),2))))));
}
double phi_n = phi1;
double lambda_n = Math.atan(Y_n/X_n);
/**********************************************************************************************************************/
/** 4) coordonnées géographiques NTF (phi_n,lambda_n) coordonnées projetées en Lambert II étendu (X_l2e, Y_l2e) **/
/**********************************************************************************************************************/
// Nous utilisons les formules de projection et les constantes fournies par l'IGN dans un autre document :
// avant tout on définit quelques constantes
double n = 0.7289686274;
double c = 11745793.39;
double Xs = 600000.0;
double Ys = 8199695.768;
double e_n = Math.sqrt(e2_n);
double lambda0 = 0.04079234433198; //correspond à la longitude en radian de Paris (2°20'14.025" E) par rapport à Greenwich
// puis on calcule la latitude isométrique
double L = Math.log(Math.tan(Math.PI/4 + phi_n/2) * Math.pow(((1-e_n*Math.sin(phi_n))/(1+e_n*Math.sin(phi_n))),(e_n/2)));
// enfin on projette
double X_l2e = Xs + c*Math.exp((-n*L))*Math.sin(n*(lambda_n-lambda0));
double Y_l2e = Ys - c*Math.exp((-n*L))*Math.cos(n*(lambda_n-lambda0));
double[]tabXY = new double[2];
tabXY[0] = X_l2e;
tabXY[1] = Y_l2e;
return tabXY;
}
} |
Partager