1. 程式人生 > >高斯克呂格投影,將經緯度轉換為投影座標

高斯克呂格投影,將經緯度轉換為投影座標

1. 預先定義橢球引數

 1         /************************************************************************/
 2         /*                                 84橢球                               */
 3         double a = 6378137.0;
 4         double b = 6356752.3142;
 5         double f = 1 / 298.257223563;
 6         double
c = a / (1 - f); 7 double e2 = 2 * f - f * f; 8 double ep2 = 1 / ((1 - f) * (1 - f)) - 1; 9 double dpi = 0.0174532925199432957692; 10 double m0 = a * (1 - e2); 11 double m2 = 1.5 * e2 * m0; 12 double m4 = 1.25 * e2 * m2; 13 double m6 = 7 * e2 * m4 / 6
; 14 double m8 = 9 * e2 * m6 / 8; 15 double a0 = m0 + m2 / 2 + 3 * m4 / 8 + 5 * m6 / 16 + 35 * m8 / 128; 16 double a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16; 17 double a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32; 18 double a6 = m6 / 32 + m8 / 16; 19 double ml = 114
; // 中央經線 根據需要進行設定 20 21 double _iproj[12] = { a, b, f, c, e2, ep2, dpi, a0, a2, a4, a6, ml };

2. 編寫轉換函式

 1     int latlon2xy(double &x, double &y)
 2     {
 3         double lat = x * _iproj[6];
 4         double dL = (y - _iproj[11]) * _iproj[6];
 5 
 6         double X = _iproj[7] * lat - _iproj[8] * sin(2 * lat) / 2 + _iproj[9] * sin(4 * lat) / 4 - _iproj[10] * sin(6 * lat) / 6;
 7         double tn = tan(lat);
 8         double tn2 = tn * tn;
 9         double tn4 = tn2 * tn2;
10 
11         double j2 = (1 / pow(1 - _iproj[2], 2) - 1) * pow(cos(lat), 2);
12         double n = _iproj[0] / sqrt(1.0 - _iproj[4] * sin(lat) * sin(lat));
13 
14         double temp[6] = { 0 };
15         temp[0] = n * sin(lat) * cos(lat) * dL * dL / 2;
16         temp[1] = n * sin(lat) * pow(cos(lat), 3) * (5 - tn2 + 9 * j2 + 4 * j2 * j2) * pow(dL, 4) / 24;
17         temp[2] = n * sin(lat) * pow(cos(lat), 5) * (61 - 58 * tn2 + tn4) * pow(dL, 6) / 720;
18         temp[3] = n * cos(lat) * dL;
19         temp[4] = n * pow(cos(lat), 3) * (1 - tn2 + j2) * pow(dL, 3) / 6;
20         temp[5] = n * pow(cos(lat), 5) * (5 - 18 * tn2 + tn4 + 14 * j2 - 58 * tn2 * j2) * pow(dL, 5) / 120;
21 
22         y = X + temp[0] + temp[1] + temp[2];
23         x = temp[3] + temp[4] + temp[5];
24     }