高斯坐标转经纬度的c++实现,不使用迭代法
时间: 2024-03-25 13:35:34 浏览: 159
高斯投影坐标转换成经纬度的公式是复杂的,但是可以使用以下的公式进行计算:
```c++
double PI = 3.14159265358979323846;
double a = 6378137.0; // 长半轴
double b = 6356752.3142; // 短半轴
double e = sqrt(1 - pow(b, 2) / pow(a, 2)); // 第一偏心率
double e2 = sqrt(pow(a, 2) / pow(b, 2) - 1); // 第二偏心率
double n = (a - b) / (a + b);
double n2 = pow(n, 2);
double n3 = pow(n, 3);
double n4 = pow(n, 4);
double n5 = pow(n, 5);
double M = 0;
if (phi == 0) {
M = 0;
}
else {
M = a * ((1 - pow(e, 2) / 4 - 3 * pow(e, 4) / 64 - 5 * pow(e, 6) / 256) * phi
- (3 * pow(e, 2) / 8 + 3 * pow(e, 4) / 32 + 45 * pow(e, 6) / 1024) * sin(2 * phi)
+ (15 * pow(e, 4) / 256 + 45 * pow(e, 6) / 1024) * sin(4 * phi)
- (35 * pow(e, 6) / 3072) * sin(6 * phi));
}
double n2sin2phi = n2 * pow(sin(phi), 2);
double nu = a / sqrt(1 - n2sin2phi);
double p = lon - lon0;
double sinphi = sin(phi);
double cosphi = cos(phi);
double A = a * (1 - n + 5 * (n2 - n3) / 4 + 81 * (n4 - n5) / 64) * phi
- a * (3 * n - 3 * (n2 + n3) / 2 + 27 * (n4 - n5) / 32) * sin(2 * phi)
+ a * (15 * (n2 - n3) / 16 - 15 * (n4 - n5) / 32) * sin(4 * phi)
- a * 35 * (n4 - n5) / 48 * sin(6 * phi);
double T = pow(tan(phi), 2);
double C = e2 * pow(cosphi, 2) / (1 - pow(e, 2));
double A1 = (lon - lon0) * cosphi;
double A2 = A1 * pow(cosphi, 2) * (1 - T + C);
double A3 = A1 * pow(cosphi, 4) * (5 - 18 * T + pow(T, 2) + 72 * C - 58 * n2sin2phi);
double A4 = A1 * pow(cosphi, 6) * (61 - 58 * T + pow(T, 2) + 600 * C - 330 * n2sin2phi);
double x = A + nu / pow(2, 1) * sinphi * cosphi * A2
+ nu / pow(4, 1) * sinphi * pow(cosphi, 3) * A3
+ nu / pow(6, 1) * sinphi * pow(cosphi, 5) * A4;
double sinphi1 = x / (nu * (1 - n2sin2phi / (a * pow(cosphi, 2))));
double phi1 = asin(sinphi1);
double sinphi2 = sinphi1;
double phi2 = phi1 + 1;
while (abs(phi2 - phi1) > pow(10, -15)) {
phi1 = phi2;
nu = a / sqrt(1 - n2 * pow(sinphi1, 2));
sinphi1 = x / (nu * (1 - n2 * pow(cosphi, 2)));
phi2 = asin(sinphi1);
}
double lambda = lambda0 + atan2(A2 - pow(A1, 2) * sinphi * cosphi * (1 + C) / nu, pow(cosphi, 3) * (1 + C) * (A3 / nu - pow(A1, 2) * sinphi * pow(cosphi, 2) * (2 + C) / (nu * pow(2, 1)))
+ pow(cosphi, 5) * (A4 / nu - pow(A1, 2) * sinphi * pow(cosphi, 4) * (9 + 4 * C) / (nu * pow(4, 1))));
double lat = phi2 / PI * 180;
double lng = lambda / PI * 180;
```
其中,输入的参数为高斯投影坐标的x、y值,以及中央子午线的经度lon0。输出的结果为经度lng和纬度lat。
阅读全文