须知
利用GMT的mapproject
转换经纬度和UTM坐标时,首先要明确以下几点:
- UTM坐标系横轴为Y,纵轴为X,某点的坐标值为
(Y, X)
,与常规笛卡尔坐标系符号表示相反 - UTM坐标单位为
m
,Y坐标不加带号共6位,X坐标为7位
- 需要知道待转换坐标所在的UTM分带号,在东半球其计算方法如下,
n=int()+31, lon为研究区中央经线经度
经纬度 <=> UTM
本实例中,坐标点(121,32)所在UTM带号为51.
# 经纬度 => UTM
echo 121 32 | gmt mapproject -Ju51/1:1 -C -F
输出结果如下:
311072.361931 3542183.49115
# UTM => 经纬度
echo 311072.4 3542183.5 | gmt mapproject -Ju51/1:1 -C -F -I
输出结果如下:
121.000000401 32.0000000858
经纬度 <=> 高斯克吕格
- 高斯-克吕格为国内地形图常用平面投影坐标系,该投影方式和UTM非常相似,区别仅仅为中央经线的比例因子不同。
- 高斯-克吕格中央经线长度比为1,而UTM中央经线长度比为0.9996。
- GMT中使用
-Ju
投影时,默认比例因子为UTM的0.9996。如果想投影为高斯-克吕格坐标,须使用PROJ_SCALE_FACTOR
指定高斯-克吕格的比例因子1。
# 经纬度转换为Gauss-Kruger
echo 121 32 | gmt mapproject -Ju51/1:1 -C -F --PROJ_SCALE_FACTOR=1
输出结果如下:
310996.760635 3543600.93152
# 将Gauss-Kruger坐标转换为经纬度
echo 310996.760635 3543600.93152 | gmt mapproject -Ju51/1:1 -C -F -I --PROJ_SCALE_FACTOR=1
输出结果如下:
121 31.9999999997
结果对比
PS:椭球参数问题
WGS84
一般来讲,GPS设备提供的坐标(B,L,H)是1984年世界大地坐标系(WordGeodetic System 1984即WGS-84)的坐标,为地心坐标系。其中B为纬度,L为经度,H为大地高即是到WGS-84椭球面的高度(latitude_longitude)。-
北京54 & 西安80
在过去,我国地图采用的是1954北京、1980西安参心坐标系下的高斯投影坐标(x,y),也有一些采用1954北京、1980西安坐标系下的经纬度坐标(B,L),高程一般为海拔高度h。GPS的测量结果与我国的54系或80系坐标相差几十米至一百多米,随区域不同,差别也不同,经粗落统计,我国西部相差70米左右,东北部140米左右,南部75米左右,中部45米左右。
-
2000坐标系
当前我国采用的坐标系统为2000国家大地坐标系,为全球地心坐标系。以后新提交的地质资料必须采用2000坐标系。关于2000坐标系和WGS84坐标系的差别,魏子卿院士在《2000中国大地坐标系及其与WGS84的比较》一文进行了比较,认为CGCS2000与WGS84是相容的,在坐标系实现精度范围内二者坐标是一致的。
北斗坐标系
北斗导航系统独立于国家大地坐标系,采用自己的专用坐标系,但其椭球参数使用CGCS2000。详情见魏子卿院士文章《北斗坐标系》。
各坐标系椭球参数如下:
北京54 | 西安80 | WGS84 | CGCS2000 | |
---|---|---|---|---|
长半轴a | 6378245 | 6378140 | 6378137 | 6378137 |
短半轴b | 6356863.0188 | 6356755.2882 | 6356752.3142451795 | 6356752.31414 |
扁率f |
-
GMT进行坐标转换时需要注意的事情:
确定图件的椭球参数
GMT默认的椭球参数为WGS84,当使用mapproject
进行坐标转换时,如果投影平面坐标图件用的是西安80、北京54,或者国家2000椭球体等非WGS84椭球体时,需要在GMT中利用PROJ_ELLIPSOID
指定相应的椭球参数。
例如:
使用GMT将2000坐标系下的经纬度转换为UTM或高斯-克吕格坐标时,需要先使用以下代码定义椭球参数(虽然2000坐标系和GMT默认的WGS84几乎是一样的):
gmt set PROJ_ELLIPSOID=6378137.0,298.257222101
坐标转换要在同一个椭球体内进行
经纬度坐标、UTM坐标、高斯克吕格坐标,只能在同一椭球体下相互转换。跨椭球体转换需要当地控制点数据,而控制点是保密的,需要向测绘局申请。