GPS坐标转换

前段时间，工作需要学习了一下GPS的坐标转换，发现这个转换过程不算难，但就是概念很容易混淆，结果无论是上网查资料，还是问一些资深人士，都能碰到或多或少的一些概念上的错误。

    经过一段时间的学习，综合了不少资料，我终于搞明白了其中的道理，不敢独享，决定写出来让大家批判一下，希望对大家有所帮助，同时也希望如果各位发现我的错误，也及时通知我。

    说到转换，我不得不先提几个概念：

1.     84坐标系――美国的GPS就用这个坐标系，也就是说我们从GPS接收过来的地心大地坐标就是在这个坐标系下的；

2.     54坐标系――中国常用的坐标系，一般来说，我们就是要得到该坐标系下的平面直角坐标；

3.     地心大地坐标――就是经纬度表示的大地位置数据（B、L、H）；

4.     地心直角坐标－－是以地心为原点，但表示方式采样直角坐标（X、Y、Z）；

5.     平面直角坐标――不用多解释了，就是我们平常用来表示位置的坐标（x、y、z）；

6.     椭球参数――地心大地坐标转换到地心直角坐标需要使用的参数；

7.     七参数――不同椭球系下的地心直角坐标互换所使用到的参数；

8.     投影参数――同一椭球系下，地心大地坐标转换为平面直角坐标所使用的参数；

9.     四参数――同一椭球系下，不同平面坐标系之间的转换所使用的参数。

    有了以上概念，只要知道转换过程就OK了。下面我就不罗嗦了，直接用方向线来说明过程，其中我以常用的北京54坐标系作为目标坐标系：

从GPS中接收到84坐标系下的地心大地坐标

↓

使用84坐标系的椭球参数转换为84坐标系下的地心直角坐标

↓

使用七参数转换为54坐标系下的地心直角坐标

↓

使用54坐标系的椭球参数转换为54坐标系下的地心大地坐标

↓

使用投影参数转换为54坐标系下的平面直角坐标

↓

使用四参数转换为其它区域坐标系下的平面直角坐标

上一章已经简单说明了各坐标系的基本概念和相互的转换关系（其实也不叫学习笔记了，所有转换都做完了，只是一直在忙，没有时间一次写完，只能有时间就写一章，一点一点来，实在愧对关注的朋友）。

    这一章就说一下84坐标系下的地心大地坐标如何使用84坐标系的椭球参数转换为84坐标系下的地心直角坐标（好拗口……）。

    先来看看84坐标系的椭球参数：

        椭圆长半径A:6378137.000000

        椭圆短半径B:6356752.314245179497

    根据椭圆长短半径，我们很容易就能算出椭圆第一偏心率e1，计算公式如下：

       e1 ＝ sqrt(A×A－B×B)/ A；

    注意这里sqrt是开方的意思。

    OK，上面我们通过一个简单的公式就计算出了椭圆的第一偏心率e1，下面，我们就利用这个参数把地心大地坐标转换为地心直角坐标。

    已知：B、L、H地心大地坐标

    求：X、Y、Z地心直角坐标

    计算步骤：

1．参数N：N＝A/sqrt（1－e1×e1×sin（B）×sin（B））；

2．X：X＝（N＋H）×cos（B）×cos（L）；

3．Y：Y＝（N＋H）×cos（B）×sin（L）；

4．Z：Z＝（N×（1－e1×e1）＋H）×sin（B）；

    好了，就是这么简单，BLH方式表示的大地坐标一下就转换为XYZ表示的直角坐标了。不过大家要注意的是，就是BL这两个数值是角度值，在C编程中，cos和sin函数要用弧度值，这个角度转换问题要注意。 

上次说到在84坐标系下如何利用84坐标系的椭球参数把地心大地坐标转换为地心直角坐标。有了84坐标系下的地心直角坐标就可以考虑如何得到54坐标系下的地心直角坐标了（不一定是54坐标系，不过为了描述方便，我选了54坐标系而已）。

    所以这次就讲解一下如何把84坐标系下的地心直角坐标转换成54坐标系下的地心直角坐标。

    转换的方法有很多，但常用的是七参数法，这个方法还有一个外文翻译过来的名称，但我现在忘记了，也懒得追查，反正是一种空间转换方法，顾名思义，需要七个参数才能进行转换，七个参数分别是：

△X △Y △Z——三个坐标方向的平移参数；

αβγ——三个方向的旋转角参数

k——尺度参数

这七个参数不是固定的，而是不同的地区有不同数值，知道的当然可以直接填写，不知道的一般就会找几个已知地方坐标的位置获取GPS采集的坐标再反算这七个参数，如何反算以后我再写。

闲话少说，有了七参数后，把84坐标系下的地心直角坐标转换成54坐标系下的地心直角坐标的计算方式是：

    已知：84坐标系下的地心直角坐标x，y，z

    求：54坐标系下的地心直角坐标X、Y、Z

       X = △X + k*x - β*z + γ*y + x

    Y = △Y + k*y + α*z - γ*x + y

    Z = △Z + k*z - α*y + β*x + z

    公式是推导的，倒没有什么特别，主要是里面的计算单位不要弄错，角度就肯定是弧度的，你别用了角度之类的单位，而那个k的尺度参数单位是ppm，一般坐标都是m为单位，所以k还有除以1000000才能参与计算。