高斯坐标反算公式

2013-05-17 13:47阅读：

http://blog.sina.cn/dpool/blog/u/1708811004

最近需要从地方平面坐标反算到WGS84的经纬度，反算的过程个人认为比正算工程要坎坷，所以特此记下来，以便以后查阅使用。
坐标反算的流程如下：地方平面坐标系四参数转换 WGS84投影下的平面坐标系高斯反算 WGS84球面坐标系。具体流程如下:
四参数转换很简单，就不做介绍，一搜网上一大堆，经过四参数转换得到WGS84投影下的平面坐标系下的X和Y。
第二步，Y值要减去500000，再进行下一步运算，这个是因为高斯投影正算的坐标Y值都加上了500000，以免Y值是负值。
第三步，高斯反算，高斯反算公式如下：

高斯坐标反算公式

在这里我们看到了底点纬度，所以第三步的第一小步应该是先求底点纬度，在网上找了很多底点纬度算法公式，大多都结果都差很多，最终找到一个比较准确的方法，链接http://hi.baidu.com/yiyiyis/item/95a590cd9959580a0bd93a2f,感谢原文主人，我自

己写的函数如下:
/// <summary>
/// 计算底点纬度公式
/// </summary>
/// <param name='m_pWGS84a'>椭球长半径</param>
/// <param name='m_pWGS84f'>椭球扁率的倒数</param>
/// <param name='x'>所计算点的横坐标</param>
/// <returns>返回底点纬度</returns>
private static double CalculateBf(double m_pWGS84a, double m_pWGS84f, double x)
{
double currf = 1 / m_pWGS84f; //扁率
double currb = m_pWGS84a * (1 - currf); //b,短半轴
double e2 = (m_pWGS84a * m_pWGS84a - currb * currb) / (m_pWGS84a * m_pWGS84a); //e的平的平方
double ee2 = (m_pWGS84a * m_pWGS84a - currb * currb) / (currb * currb); //e’的平方
double e4 = e2 * e2;
double e6 = e2 * e2 * e2;
double e8 = Math.Pow(e2, 4);
double e10 = Math.Pow(e2, 5);
double e12 = Math.Pow(e2, 6);
double e14 = Math.Pow(e2, 7);
double e16 = Math.Pow(e2, 8);
double c0 = 1 + e2 / 4 + 7 * e4 / 64 + 15 * e6 / 256 + 579 * e8 / 16384 + 1515 * e10 / 65536 + 16837 * e12 / 1048576 + 48997 * e14 / 4194304 + 9467419 * e16 / 1073741824;
c0 = m_pWGS84a / c0;
double b0 = x / c0;
double d1 = 3 * e2 / 8 + 45 * e4 / 128 + 175 * e6 / 512 + 11025 * e8 / 32768 + 43659 * e10 / 131072 + 693693 * e12 / 2097152 + 10863435 * e14 / 33554432;
double d2 = -21 * e4 / 64 - 277 * e6 / 384 - 19413 * e8 / 16384 - 56331 * e10 / 32768 - 2436477 * e12 / 1048576 - 196473 * e14 / 65536;
double d3 = 151 * e6 / 384 + 5707 * e8 / 4096 + 53189 * e10 / 163840 + 4599609 * e12 / 655360 + 15842375 * e14 / 1048576;
double d4 = -1097 * e8 / 2048 - 1687 * e10 / 640 - 3650333 * e12 / 327680 - 114459079 * e14 / 27525120;
double d5 = 8011 * e10 / 1024 + 874457 * e12 / 98304 + 216344925 * e14 / 3670016;
double d6 = -682193 * e12 / 245760 - 46492223 * e14 / 1146880;
double d7 = 36941521 * e14 / 3440640;
double bf = b0 + Math.Sin(2 * b0) * (d1 + Math.Sin(b0) * Math.Sin(b0) * (d2 + Math.Sin(b0) * Math.Sin(b0) * (d3 + Math.Sin(b0) * Math.Sin(b0) * (d4 + Math.Sin(b0) * Math.Sin(b0) * (d5 + Math.Sin(b0) * Math.Sin(b0) * (d6 + d7 * Math.Sin(b0) * Math.Sin(b0)))))));
return bf;
}
求出底点纬度后再带入到反算公式，即可得到一个B，L值。
最后一步就比较容易了，B值和L值都是弧度单位，将B值转化为我们需要的单位，L值要加上该地方的中央子午线的值，再转换为我们需要的单位即可。

举报/Report

我的更多文章

下载客户端阅读体验更佳