高斯正反算—投影坐标转大地坐标、大地坐标转投影坐标(附有完整代码及测试结果)

大家好，欢迎来到IT知识分享网。

        本文意在介绍高斯正反算的基本原理和代码实现！针对高斯正反算网上给出的方法很多，但是我试了之后发现多少都有些问题：或公式原理问题、或精度问题！

        通过查找资料对其进行了总结与测试！原理！代码！测试结果！本文都一一给出，此外本文还对常见坐标系：北京54、西安80、WGS84、CGCS2000等坐标系的高斯转换都做出了实现并使用 QT 进行封装可视化！

        对比公式和代码实现可快速理解高斯正反算！

        精度可达 0.0001!堪称最全公式原理总结！

        坐标转换整套流程包括：像素坐标转投影坐标、投影坐标转大地坐标、大地坐标转空间直角坐标、七参数转换、空间直角坐标转大地坐标、大地坐标转投影坐标、投影坐标转像素坐标； 本人均已实现，且每一个环节都已经过测试、如有需要欢迎在下方留言评论！！！

一、常用椭球参数

二、高斯正反算原理 

1、高斯正算(大地坐标转投影坐标) 

 2、高斯反算(投影坐标转大地坐标)

三、高斯正反算代码实现 

#include<iostream> #include<cmath> #include "stdio.h" #define pi 3. #define p0 .51564 //wgs84参考椭球 const double e = 0.000; const double e1 = 0.008; const double b = .3141; const double a = .0; using namespace std; //大地坐标转投影坐标 void DadiPoint2ProjectPoint(double B, double L) { //把度转化为弧度 B = B * pi / 180; L = L * pi / 180; double N, t, n, c, V, Xz, m1, m2, m3, m4, m5, m6, a0, a2, a4, a6, a8, M0, M2, M4, M6, M8, x0, y0, l; int L_num; double L_center; //中央子午线经度，6°带 L_num = (int)(L * 180 / pi / 6.0) + 1; L_center = 6 * L_num - 3; //中央子午线经度，3°带 //L_num = (int)(L * 180 / pi / 3.0 + 0.5); //L_center = 3 * L_num; l = (L / pi * 180 - L_center) * 3600; //求带号、中央经线、经差 M0 = a * (1 - e); M2 = 3.0 / 2.0 * e * M0; M4 = 5.0 / 4.0 * e * M2; M6 = 7.0 / 6.0 * e * M4; M8 = 9.0 / 8.0 * e * M6; a0 = M0 + M2 / 2.0 + 3.0 / 8.0 * M4 + 5.0 / 16.0 * M6 + 35.0 / 128.0 * M8; a2 = M2 / 2.0 + M4 / 2 + 15.0 / 32.0 * M6 + 7.0 / 16.0 * M8; a4 = M4 / 8.0 + 3.0 / 16.0 * M6 + 7.0 / 32.0 * M8; a6 = M6 / 32.0 + M8 / 16.0; a8 = M8 / 128.0; Xz = a0 * B - a2 / 2.0 * sin(2 * B) + a4 / 4.0 * sin(4 * B) - a6 / 6.0 * sin(6 * B) + a8 / 8.0 * sin(8 * B); //计算子午线弧长 c = a * a / b; V = sqrt(1 + e1 * cos(B) * cos(B)); N = c / V; t = tan(B); n = e1 * cos(B) * cos(B); m1 = N * cos(B); m2 = N / 2.0 * sin(B) * cos(B); m3 = N / 6.0 * pow(cos(B), 3) * (1 - t * t + n); m4 = N / 24.0 * sin(B) * pow(cos(B), 3) * (5 - t * t + 9 * n); m5 = N / 120.0 * pow(cos(B), 5) * (5 - 18 * t * t + pow(t, 4) + 14 * n - 58 * n * t * t); m6 = N / 720.0 * sin(B) * pow(cos(B), 5) * (61 - 58 * t * t + pow(t, 4)); x0 = Xz + m2 * l * l / pow(p0, 2) + m4 * pow(l, 4) / pow(p0, 4) + m6 * pow(l, 6) / pow(p0, 6); y0 = m1 * l / p0 + m3 * pow(l, 3) / pow(p0, 3) + m5 * pow(l, 5) / pow(p0, 5); //计算x y坐标 double x = x0; //double y = y0 +  +  * L_num; //化为国家统一坐标 double y = y0 + ; //化为国家统一坐标 cout << "方法一 x=" << x << endl; cout << "方法一 y=" << y << endl; } //投影坐标转大地坐标 void ProjectPoint2DadiPoint(double x, double y, double l0) { //l0为中央经度 double Bf, B0, FBf, M, N, V, t, n, c, y1, n1, n2, n3, n4, n5, n6, a0, a2, a4, a6, M0, M2, M4, M6, M8, l; int L_num, L_center; L_num = (int)(x / .0); y1 = y - ; //y1 = y -  - L_num * ; //L_center = ((L_num + 1) * 6 - 3)*pi*180; //中央子午线经度，6°带 //cout<<"L_center="<<L_center<<endl; //L_center = L_num * 3; //中央子午线经度，3°带 M0 = a * (1 - e); M2 = 3.0 / 2.0 * e * M0; M4 = 5.0 / 4.0 * e * M2; M6 = 7.0 / 6.0 * e * M4; M8 = 9.0 / 8.0 * e * M6; a0 = M0 + M2 / 2.0 + 3.0 / 8.0 * M4 + 5.0 / 16.0 * M6 + 35.0 / 128.0 * M8; a2 = M2 / 2.0 + M4 / 2 + 15.0 / 32.0 * M6 + 7.0 / 16.0 * M8; a4 = M4 / 8.0 + 3.0 / 16.0 * M6 + 7.0 / 32.0 * M8; a6 = M6 / 32.0 + M8 / 16.0; cout << "a0=" << a0 << endl; cout << "a2=" << a2 << endl; cout << "a4=" << a4 << endl; cout << "a6=" << a6 << endl; Bf = x / a0; B0 = Bf; cout<<"B0="<<B0<<endl; cout<<"sin(2 * B0)="<<sin(2 * B0)/2<<endl; while ((fabs(Bf - B0) > 0.0000001) || (B0 == Bf)) { B0 = Bf; FBf = -a2 / 2.0 * sin(2 * B0) + a4 / 4.0 * sin(4 * B0) - a6 / 6.0 * sin(6 * B0); Bf = (x - FBf) / a0; } //迭代求数值为x坐标的子午线弧长对应的底点纬度 cout<<"Bf="<<Bf<<endl; t = tan(Bf); //一样 c = a * a / b; V = sqrt(1 + e1 * cos(Bf) * cos(Bf)); //一样 N = c / V; //一样 M = c / pow(V, 3); //一样 n = e1 * cos(Bf) * cos(Bf); //一样(为n的平方) n1 = 1 / (N * cos(Bf)); n2 = -t / (2.0 * M * N); n3 = -(1 + 2 * t * t + n) / (6.0 * pow(N, 3) * cos(Bf)); n4 = t * (5 + 3 * t * t + n - 9 * n * t * t) / (24.0 * M * pow(N, 3)); n5 = (5 + 28 * t * t + 24 * pow(t, 4) + 6 * n + 8 * n * t * t) / (120.0 * pow(N, 5) * cos(Bf)); n6 = -t * (61 + 90 * t * t + 45 * pow(t, 4)) / (720.0 * M * pow(N, 5)); //秒 double B = (Bf + n2 * y1 * y1 + n4 * pow(y1, 4) + n6 * pow(y1, 6)) / pi * 180; double L0=l0; l = n1 * y1 + n3 * pow(y1, 3) + n5 * pow(y1, 5); //double L = L_center + l / pi * 180; //反算得大地经纬度 double L = L0 + l / pi * 180; //反算得大地经纬度 cout << "方法一 B=" << B << endl; cout << "方法一 L=" << L << endl; } 

说明：高斯正算中的输入为度； 

 四、高斯正反算结果

1、已知旧坐标系投影坐标数据

 2、高斯反算(投影坐标转大地坐标)

 3、高斯正算(大地坐标转投影坐标)

 说明：高斯正反算经度由上述转换结果对比可知；

常见坐标系：北京54、西安80、WGS84、CGCS2000等坐标系的高斯转换都做出了实现并使用 QT 进行封装可视化：高斯正反算设计实现!!!-C++文档类资源-CSDN下载

坐标转换整套流程包括：像素坐标转投影坐标、投影坐标转大地坐标、大地坐标转空间直角坐标、七参数转换、空间直角坐标转大地坐标、大地坐标转投影坐标、投影坐标转像素坐标； 本人均已实现，且每一个环节都已经过测试、如有需要欢迎在下方留言评论！！！

如果需要代码包，请在评论区留言！！！ 

如果需要代码包，请在评论区留言！！！ 

如果需要代码包，请在评论区留言！！！ ​​​​​​​