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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
|
class ProjectionHelper
{
/// <summary>
///
/// </summary>
/// <param name="d_long"></param>
/// <param name="m_long"></param>
/// <param name="s_long"></param>
/// <param name="orientation_long"></param>
/// <param name="d_lat"></param>
/// <param name="m_lat"></param>
/// <param name="s_lat"></param>
/// <param name="orientation_lat"></param>
/// <returns></returns>
public 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) **/
/**************************************************************************************************************/
// J'ai utilisé 2 formules données par l'IGN dans un de leur document ainsi que deux constantes de
// l'ellipsoide de référence du système WGS84 (les deux demi-axes) :
double a_w = 6378137.0;
double b_w = 6356752.314;
// 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'un translation à effectuer :
// on a donc
double dX = 168.0;
double dY = 60.0;
double dZ = -320.0;
// et...
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 p0 = 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 p1 = Math.Atan((Z_n / Math.Sqrt(X_n * X_n + Y_n * Y_n)) / (1 - (a_n * e2_n * Math.Cos(p0)) / (Math.Sqrt((X_n * X_n + Y_n * Y_n) * (1 - e2_n * Math.Pow(Math.Sin(p0), 2))))));
while (!(Math.Abs(p1 - p0) < epsilon))
{
p0 = p1; p1 = Math.Atan((Z_n / Math.Sqrt(X_n * X_n + Y_n * Y_n)) / (1 - (a_n * e2_n * Math.Cos(p0)) / (Math.Sqrt((X_n * X_n + Y_n * Y_n) * (1 - e2_n * Math.Pow(Math.Sin(p0), 2))))));
}
double phi_n = p1;
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) **/
/**********************************************************************************************************************/
// J'utilise 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;
}
/// <summary>
///
/// </summary>
/// <param name="lambda_w"></param>
/// <param name="phi_w"></param>
/// <returns></returns>
public static double[] WGS84ToLambert2e(double longitude, double latitude)
{
// Le système de coordonnées géographiques utilisé est postif vers le Nord et vers l'Est
var orientation_lat = 'N';
if (latitude < 0)
{
orientation_lat = 'S';
latitude = -1 * latitude;
}
var degree_lat = Math.Truncate(latitude);
var decPart = latitude - degree_lat;
var multipliedBy3600 = decPart * 3600;
var dividedBy60 = multipliedBy3600 / 60.0;
var minute_lat = Math.Truncate(dividedBy60);
decPart = dividedBy60 - minute_lat;
var seconde_lat = decPart * 60;
// Le système de coordonnées géographiques utilisé est postif vers le Nord et vers l'Est
var orientation_long = 'E';
if (longitude < 0)
{
orientation_long = 'W';
longitude = -1 * longitude;
}
var degree_long = Math.Truncate(longitude);
decPart = longitude - degree_long;
multipliedBy3600 = decPart * 3600;
dividedBy60 = multipliedBy3600 / 60.0;
var minute_long = Math.Truncate(dividedBy60);
decPart = dividedBy60 - minute_long;
var seconde_long = decPart * 60;
return WGS84toLambert2e(degree_long, minute_long, seconde_long, orientation_long, degree_lat, minute_lat, seconde_lat, orientation_lat);
}
/// <summary>
///
/// </summary>
/// <param name="latitude"></param>
/// <param name="longitude"></param>
/// <returns></returns>
public static double[] WGS84ToLambert93(double latitude, double longitude)
{
/**** Conversion latitude,longitude en coordonée lambert 93 ****/
// Projection conforme sécante, algo détailler dans NTG_71.pdf : http://www.ign.fr/affiche_rubrique.asp?rbr_id=1700&lng_id=FR
// > ACCUEIL > L'offre IGN Pro > Géodésie > RGF93 > Outils
//variables:
var a = 6378137; //demi grand axe de l'ellipsoide (m)
var e = 0.08181919106; //première excentricité de l'ellipsoide
var l0 = (Math.PI / 180) * 3;
var lc = l0;
var phi0 = (Math.PI / 180) * 46.5; //latitude d'origine en radian
var phi1 = (Math.PI / 180) * 44; //1er parallele automécoïque
var phi2 = (Math.PI / 180) * 49; //2eme parallele automécoïque
var x0 = 700000; //coordonnées à l'origine
var y0 = 6600000; //coordonnées à l'origine
var phi = (Math.PI / 180) * latitude;
var l = (Math.PI / 180) * longitude;
//calcul des grandes normales
var gN1 = a / Math.Sqrt(1 - e * e * Math.Sin(phi1) * Math.Sin(phi1));
var gN2 = a / Math.Sqrt(1 - e * e * Math.Sin(phi2) * Math.Sin(phi2));
//calculs des latitudes isométriques
var gl1 = Math.Log(Math.Tan(Math.PI / 4 + phi1 / 2) * Math.Pow((1 - e * Math.Sin(phi1)) / (1 + e * Math.Sin(phi1)), e / 2));
var gl2 = Math.Log(Math.Tan(Math.PI / 4 + phi2 / 2) * Math.Pow((1 - e * Math.Sin(phi2)) / (1 + e * Math.Sin(phi2)), e / 2));
var gl0 = Math.Log(Math.Tan(Math.PI / 4 + phi0 / 2) * Math.Pow((1 - e * Math.Sin(phi0)) / (1 + e * Math.Sin(phi0)), e / 2));
var gl = Math.Log(Math.Tan(Math.PI / 4 + phi / 2) * Math.Pow((1 - e * Math.Sin(phi)) / (1 + e * Math.Sin(phi)), e / 2));
//calcul de l'exposant de la projection
var n = (Math.Log((gN2 * Math.Cos(phi2)) / (gN1 * Math.Cos(phi1)))) / (gl1 - gl2);//ok
//calcul de la constante de projection
var c = ((gN1 * Math.Cos(phi1)) / n) * Math.Exp(n * gl1);//ok
//calcul des coordonnées
var ys = y0 + c * Math.Exp(-1 * n * gl0);
var x93 = x0 + c * Math.Exp(-1 * n * gl) * Math.Sin(n * (l - lc));
var y93 = ys - c * Math.Exp(-1 * n * gl) * Math.Cos(n * (l - lc));
double[] tabXY = new double[2];
tabXY[0] = x93;
tabXY[1] = y93;
return tabXY;
}
} |
Conversion de projections WGS84 à Lambert II étendu et Lambert 93
Répondre avec citation
Partager