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;
}
} |
Partager