|
本文意在介绍高斯正反算的基本原理和代码实现！针对高斯正反算网上给出的方法很多，但是我试了之后发现多少都有些问题：或公式原理问题、或精度问题！
通过查找资料对其进行了总结与测试！原理！代码！测试结果！本文都一一给出，此外本文还对常见坐标系：北京54、西安80、WGS84、CGCS2000等坐标系的高斯转换都做出了实现并使用 QT 进行封装可视化！
对比公式和代码实现可快速理解高斯正反算！
精度可达 0.0001!堪称最全公式原理总结！
坐标转换整套流程包括：像素坐标转投影坐标、投影坐标转大地坐标、大地坐标转空间直角坐标、七参数转换、空间直角坐标转大地坐标、大地坐标转投影坐标、投影坐标转像素坐标； 本人均已实现，且每一个环节都已经过测试、如有需要欢迎在下方留言评论！！！
一、常用椭球参数
北京54坐标系 西安80坐标系 WGS坐标系 CGC2000坐标系 a 6378245.0000000000 6378140.0000000000 6378137.0000000000 6378137.0000000000 b 6356863.0187730473 6356755.2881575287 6356752.3142 6356752.314 f 1/298.3 1/298.257 1/298.257223563 1/298.257222101 c 6399698.9017827110 6399596.6519880105 6399593.6258 6399593.6259 e1 0.006693421622966 0.006694384999588 0.00669437999013 0.00669438002290 e2 0.006738525414683 0.006739501819473 0.00673949674227 0.00673949677548
二、高斯正反算原理
1、高斯正算(大地坐标转投影坐标)
2、高斯反算(投影坐标转大地坐标)
三、高斯正反算代码实现
#include
#include
#include "stdio.h"
#define pi 3.141592653589793238463
#define p0 206264.8062470963551564
//wgs84参考椭球
const double e = 0.00669438002290;
const double e1 = 0.00673949677548;
const double b = 6356752.3141;
const double a = 6378137.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 + 500000 + 1000000 * L_num; //化为国家统一坐标
double y = y0 + 500000; //化为国家统一坐标
cout |