在GIS工作中,很多非测绘学科的同学很难高清各种坐标系,下面就一些基础知识与误区进行科普。
大地坐标系
现在一提到空间坐标,大家一般想到的就是手机GPS里的经纬度,如(116.6832,39.6789),地理稍微学得好一些的同学,还可以说出坐标x正值代表东经,负值代表西经,取值范围是-180到180,y的正值代码代表北纬,负的代表南维,取值范围(-90,90),知识再深一些的还会告诉你,0度经线经过格林威治,0度纬线就是赤道。
先不论坐标正负的取法是否符合全世界人民的共识(中国在东经、北纬,当然把自己设置在第一象限),我们先存在一个问题,但我站在一个位置,我的经纬度是否唯一呢?
很多同学可能认为是唯一的,因为我就站在这里,经纬度测不准只可能是设备的问题。但其实答案是no,这里最大的原因要怪地球,因为它不是一个正球体。
假如地球是一个完美的球体,每一个地面上的位置到地心的距离都是相等的,则球面上的任意一点,可以都可以用唯一的经纬度来标识。但是因为地球赤道半径、极低半径以及地区的半径都存在差别(就是说有些地方长,有些地方扁),所以它也不是一个标准的椭球体,你能想象用经纬度来描述一个不规则土豆表皮的各个位置么?
那么问题来了,我们还是需要一种定位的方法,这就是把地球认为用一个近似的椭球体,然后在上面画出经线圈与尾线圈进行定位。但是正如前面所述,因为地球也不是一个完美的椭球体,如果都用赤道半径作长轴,极地半径做短轴,得出的椭球体在不同海拔的位置得出的经纬度误差是不同的,与这个椭球体接近的位置精度高,在与椭球面远离的位置精度就比较低。所以取决于测量精度的要求以及不同位置的高程,可以使用不同的长轴与短轴参数进行逼近,这样才能局部的经纬度误差才比较小。
我们一般手机的定位都是采用美国GPS的WGS84(ESPG:4326,世界大地坐标系统,其坐标原点在地心)坐标系,中国也制定了如北京54(大地原点在苏联,将与苏联大地网联测后我国东北边境的三个点的坐标作为我国天文大地网起算数据,然后通过天文大地网坐标计算,推算出北京一点的坐标,故命名为北京坐标系),西安80(采用1975年国际椭球,大地原点在陕西省永乐镇,椭球面与我国境内的大地水准面密合最佳),国家2000(CGCS 2000,ESPG:4052)。不过,由于北京54与西安80坐标系都存在一定的问题,不够普适,而且我国在很长一段时间内也缺乏有效测绘卫星与设备,导致WGS84成为一种事实上的标准,民间测绘与软件(如高德地图、百度地图等)限于条件一般都会采用与GPS相同的坐标系(关于中国特色的火星坐标加密我们以后再说);值得注意的是,我国北斗系统采用的CGCS2000坐标系,如果以后做大做广,很多搞GIS的同学就可能就会经常面临不同椭球面之间坐标转换的问题了。
不同的大地坐标系,对应的球心、长半轴、短半轴、大地基准面不一样,比如WGS84是使用地心作为质心,这类坐标系被称为“地心基准面”,具有较好的普适性,但在局部地区可能精度不足;而西安80等坐标系,则是采用“区域基准面”,它们可以在特定区域内与地球表面最佳吻合,大地原点是参考椭球与大地水准面相切的点。不同球面之间坐标转换一般涉及到下面7个参数:
1. 三个平移参数ΔX、ΔY、ΔZ表示两坐标原点的平移值。
2. 三个旋转参数εx、εy、εz表示当地坐标系旋转至与地心坐标系平行时,分别绕Xt、Yt、Zt的旋转角。
3. 最后是比例校正因子,用于调整椭球大小。
投影坐标系
大地坐标系的单位一般是经纬度(分、秒),它是一种球面坐标系,衡量的是与本初面的夹角。这种坐标系下我们有些基本的测量需求无法被满足,如测量地面两个位置间的距离,你肯定不希望得到的是球心夹角,而是希望获得球面距离与面积。这时就需要投影坐标系,将球面投影成平面。
地球椭球表面是一种不可能展开的曲面,要把这样一个曲面表现到平面上,就会发生裂隙或褶皱。在投影面上,可运用经纬线的“拉伸”或“压缩”(通过数学手段)来加以避免,以便形成一幅完整的地图。但不可避免会产生变形。
地图投影的变形通常有:长度变形、面积变形和角度变形。在实际应用中,根据使用地图的目的,限定某种变形。我们最常使用的投影,是横轴墨卡托投影(UTM),关于投影的介绍,可以阅读后面的参考资料,里面介绍得很详细
POSTGIS应用
程序员在使用POSTGIS时可能经常遇到计算两个Geometry之间距离的问题,我们一般有以下几种方法:
1. 使用Geography类型,将Geometry转换成Geography
select ST_Distance(Geography(ST_GeometryFromText('Point(116.733 39.567)')),Geography(ST_GeometryFromText('Point(117.733 36.567)')));
2. 转换成投影坐标系
select ST_Distance(ST_transform(ST_GeometryFromText('Point(116.733 39.567)',4326),32650),ST_Transform(ST_GeometryFromText('Point(117.733 36.567)',4326),32650));
参考资料
http://blog.csdn.net/angelazy/article/details/44085099