高斯克吕格投影，将经纬度转换为投影坐标

Posted 2021-01-17 xingzhensun

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     }