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

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. 编写转换函数

1int 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}