坐标转换一直是空间数据处理里面一个非常重要的内容,特别是目前我国已经全面启用了CGCS2000坐标系统,以往那些54和80的坐标,未来都要统一转换到2000上面,所以很多数据处理的单位和同学,都非常关心坐标转换的问题。
虾神曾经听说地理所的一个大牛有过这样的论点——GIS大部分东西,都能在计算机专业里面找到影子,只有空间参考和投影是属于GIS自己所特有的东西。所以这个东西从来就是非地理专业与地理专业在学习和使用GIS中的一个分水岭(话说虾神作为一个纯粹的计算机专业出身的码农,当年学的时候也很痛苦……地图学原理看了好多遍,才明白了个大概)。
ArcGIS作为世界上应用最广的GIS软件,在投影转换方面的技术已经非常成熟了,但是因为中国特有的国情,导致很多国内特有的东西,他不具备——比如没有内置各种坐标系转换到CGCS2000的转换参数(一些国际特别是北美通用的转换参数,是内置的了),当然,还有国内特有的标准图幅号这种东西……
下面我们来看看,如何进行转换。
首先,转换的原理就不在这里掉书袋了,网络上很多,贴一张图意思一下:
实际上两个不同坐标系之间的转换,就是平移、旋转和比例尺度的的变化。
那么转换的方法,通常在大范围下,都是通过布尔沙沃尔夫七参数来进行转换的,数学原理(此处省略一万字和若干数学公式)……
理论研究的同学请去查阅《地图学原理》一书相关章节,下面进入工程实践操作:
ArcGIS里面,对于同椭球体下面的转换,是不需要任何参数的,比如我用WGS84(wkid:4326)转WGS84 Web Mercator(wkid: 3857),是不需要任何参数的:
但是要是换一个椭球体的话,比如换成cgcs2000,那么就需要定义地理转换参数了,如下:
当然,在新版本(10.4之后)的ArcGIS中,如果你不设定转换参数,也可以强转,只是转完之后,不保证精确度而已,而在比较老的版本里面,不设置转换参数,就直接不允许执行的。
所以,在转换的时候,定义转换参数是非常重要的,下面我们来看看怎么计算并且设定转换参数。
首先贴出公式,数学恐惧症的同学略过:
实际上就是一个七元方程组……好吧,算我没说,我们直接用软件来算。
要计算七元方程式组,就最少需要3个点来进行解算,限于写代码很痛苦,所以我们目前直接采用网上最流行的工具COORD GM来进行计算(有关心算法的,以后有机会再给大家讨论)。
这个神奇的软件如下:(PS:我很喜欢作者设置的这个图标)。
执行流程如下:
首先,准备三组坐标……准备的方法就大家各显身手了,比如用差分GPS去测量……,下面给出的三组模拟坐标,肯定是假的,仅供测试。
其中一组是WGS84的经纬度,一组是CGCS2000 3度分带 117E的投影坐标,我们就用这两组坐标来看看如何反算七参数。
首先,用ARCGIS的同学,要先明确一个问题,就是ArcGIS里面的XY坐标,与传统测量里面的XY坐标,是正好相反的,如下所示:
为什么会这样?因为传统测绘领域(包括OGC标准的GeoJSON里面,都是定义纬度在前,经度在后,原因是纬度的英文单词是latitude,简写为LAT;而经度是longitude,简写为LON……字母排序,你懂的)
明确这个问题之后,进入COORD GM操作时间:
步骤一:在COORD GM中设定要CGCS2000的椭球体参数:
CGCS2000椭球体的参数是公开的,在ArcGIS里面就有,这里只要需要输入长半轴和扁率就可以了
然后开始计算七参数了:
选定好各种参数,注意输入顺序
三组坐标输入完成之后,点击计算,得到结果:
这个就是计算出来的七参数,但是目前在ArcGIS里面还不能直接用,因为ArcGIS里面的单位与COORD里面是不一样的:
其中,ArcGIS里面的选择角度的单位是秒,而COORD里面是弧度;而ArcGIS里面的缩放比例尺用的是标准单位,而COORD里面用的百万分之,所以要进行如下转换:
后面那蓝色的就是可以在直接ArcGIS里面使用的参数了,下面把这些参数设置到ArcGIS中去:
Toolbox里面的投影转换找到创建地理变化:
按照如下信息,输入参数:
点击OK,完成创建
创建完成之后,在你的电脑的如下位置,会生成一个转换的配置文件:
C:用户你的用户名(比如Admin)AppDataRoamingEsriDesktop10.x(ArcGIS的版本)ArcToolboxCustomTransformations
比如我的机器上面,就在这个地方,配置文件的名称就是你设置的转换方法的名称:
下面就可以进行投影转换了:
然后你设置了数据和输入输出坐标系之后,系统会自动找到你设置的地理变换方法:
下面点击OK,完成坐标转换:
下面给出我这里的三组测试坐标系,大家可以用来练习一下:
No.1:
40.420897457099997,116.536157668000001
4476368.246360000222921,460634.908516999974381
No.2:
39.954091684500000,116.018334870999993,
4424893.136629999615252,416114.988767999981064
No.3:
39.834201677400003,116.712166401999994,
4411159.581939999945462,475361.417197000002488
输出结果:
DX=-0.016777
DY=0.033226
DZ=0.031063
WX=-0.0000000034
WY=0.0000000069
WZ=0.0000000064
K=-0.000000007618
ArcGIS参数
XT = -0.016777
YT = 0.033226
ZT = 0.031063
XR = -0.07015968
YR = 0.14238288
ZR = 0.13206528
SD = -0.007618
计算结果如上,如果计算结果一致,就表示你的操作步骤正确了,当然,输入的数值精度小数位不同,结果也会出现细微偏差,这个可以忽略。