本程序是海洋地质制图常用地图投影系列小程序之一,程序能用于不同基准面、3°或6°分带、单点及批量数据的高斯-克吕格投影正、反转换,正投影时的输入经纬度数据可以是度、度分及度分秒格式。
高斯-克吕格投影正、反转换动态连接库接口:
正投影接口:GS(latitude,longitude:double; var xx,yy:Double; centralat,centralon:Double; Datum,Deg:integer)
反投影接口:GSinv(xx,yy:double;var latitude,longitude:Double;centralat,centralon:Double; Datum,Deg:integer)
参数说明:latitude,longitude--纬度与经度,单位度;xx,yy--横向(经度方向)及纵向(纬度方向)坐标,单位米;centralat,centralon--原点纬度(一般为0)及中央经度,单位度;Datum--椭球体,1-北京54,2-西安80,3-WGS84;Deg:3-三度带,6-六度带。
本套系列程序最初在2004年4月完成,包括“3°、6°带高斯-克吕格投影正反转换程序”、“墨卡托投影正反转换程序”及“兰勃特等角投影正反转换程序”。2005年3月根据用户反馈作了更新,更新后增加了“UTM投影正反转换程序”,其中包括UTM直接转高斯-克吕格投影的功能,更新后反投影精度提高,反算能完全精确到小数后六位的度。2005年4月根据北京航遥一用户需求增加了“Albers等面积投影正、反转换程序”。2006年11月根据大连环境检测中心一用户要求增加了“UTM与墨卡托投影正反转换程序”,用于多波束数据的坐标转换。
此次更新在帮助文件中增加了动态连接库接口说明,供开发调用。另外在打包文件中增加了“常用地图投影转换公式”及翻译的“坐标系转换公式”,供参考。需要说明的是本套程序没有涉及坐标系的转换,只是同坐标系的投影转换。
编制这套程序原因有三:之一,本人工作中常需要投影计算,现有软件使用不太方便;之二,常发现用十进制度坐标数据作正式成果图的现象,可能是手头没有合适的投影软件所至;之三,常发现WGS84定位数据被当作北京54(克拉索夫斯基椭球体)坐标数据投影,可能是沿用早年的投影转换程序所至。这些原因促成了我编制一套简单实用、在Windows环境下的常用地图投影小程序的想法。编制完成后放在兔八哥的GIS空间站上( http://www.gissky.net),收到了不少用户反馈,我写“常用地图投影转换公式”及翻译“坐标系转换公式”应该说也是在用户建议下完成的,有个陕西的用户还给我发来了最新的“Coordinate Conversions and Transformation including Formulas”(POSC,http://www.posc.org ,国际石油技术软件开放公司),在此一并表示感谢。
有用户建议将这些小程序合并成一个总的投影转换软件,我考虑了一下最后还是保留了原来状态,一方面实际中同时用多种投影的情况不多,另一方面我的初衷是做小程序,简单而实用。
2006年12月,青岛海洋地质研究所 戴勤奋,Email:qddqinfen@cgs.gov.cn
选择投影的目的在于使所选投影的性质、特点适合于地图的用途,同时考虑地图在图廓范围内变形较小而且变形分布均匀。
我国的基本比例尺地形图(1:5千,1:1万,1:2.5万,1:5万,1:10万,1:25万,1:50万,1:100万)中,大于等于50万的多采用高斯-克吕格投影(Gauss-Kruger),这是一个等角横切椭圆柱投影,是横轴墨卡托投影(Transverse Mercator)的一个变种;小于50万的地形图采用等角正轴割园锥投影,又叫兰勃特等角投影(Lambert Conformal Conic);海上小于50万的地形图多用等角正轴圆柱投影,又叫墨卡托投影(Mercator)。一般应该采用与我国基本比例尺地形图系列一致的地图投影系统。
地图坐标系由大地基准面和地图投影确定,大地基准面是利用特定椭球体对特定地区地球表面的逼近,因此每个国家或地区均有各自的大地基准面,我们通常称谓的北京54坐标系、西安80坐标系实际上指的是我国的两个大地基准面。我国参照前苏联从1953年起采用克拉索夫斯基(Krassovsky)椭球体建立了我国的北京54坐标系,1978年采用国际大地测量协会推荐的IAG
75地球椭球体建立了我国新的大地坐标系--西安80坐标系,
目前GPS定位所得出的结果都属于WGS84坐标系统,WGS84基准面采用WGS84椭球体,它是一地心坐标系,即以地心作为椭球体中心的坐标系。相对同一地理位置,不同的大地基准面,它们的经纬度坐标是有差异的。需要说明的是,在“常用地图投影系列小程序”中,程序界面上的所谓“北京1954“西安1980”及“WGS
84”在实际计算中只涉及了相应的椭球体参数,与实际基准面无关。
本程序中采用的3个椭球体参数如下(源自“全球定位系统测量规范 GB/T 18314-2001”):
| 椭球体 | 长半轴 | 短半轴 |
| Krassovsky | 6378245 | 6356863.0188 |
| IAG 75 | 6378140 | 6356755.2882 |
| WGS 84 | 6378137 | 6356752.3142 |
椭球体与大地基准面之间的关系是一对多的关系,也就是基准面是在椭球体基础上建立的,但椭球体不能代表基准面,同样的椭球体能定义不同的基准面,如前苏联的Pulkovo 1942、索马里的Afgooye基准面都采用了Krassovsky椭球体,但它们的大地基准面显然是不同的。在目前的GIS商用软件中,大地基准面都通过当地基准面向WGS84的转换7参数来定义,即三个平移参数ΔX、ΔY、ΔZ表示两坐标原点的平移值;三个旋转参数εx、εy、εz表示当地坐标系旋转至与地心坐标系平行时,分别绕Xt、Yt、Zt的旋转角;最后是比例校正因子,用于调整椭球大小。我国的北京54、西安80相对WGS84的转换参数保密,至今没有公开。实际工作中可通过当地测绘部门获取相应参数,或利用工作区内已知的北京54或西安80坐标控制点计算,计算方法有“三参数莫洛金斯基(Molodenski)”法,“七参数布尔莎(Bursa-Wolf)”法,及多项式或线性转换等。
以中央经度为123°的(32°,121°)高斯-克吕格投影结果为例,北京54及WGS84基准面,两者投影结果在南北方向差距约63米(见下表),东西方向基本一致;西安80与WGS84投影结果基本一致。
| 输入坐标(度) |
北京54 高斯(米) |
西安80 高斯(米) | WGS84 高斯(米) | |
| 纬度值(X) | 32 | 3543664 | 3543603 | 3543601 |
| 经度值(Y) | 121 | 21310994 | 21310997 | 21310997 |
(1)高斯-克吕格投影性质
高斯-克吕格(Gauss-Kruger)投影简称“高斯投影”,又名"等角横切椭圆柱投影”,地球椭球面和平面间正形投影的一种。德国数学家、物理学家、天文学家高斯(Carl FriedrichGauss,1777一 1855)于十九世纪二十年代拟定,后经德国大地测量学家克吕格(Johannes Kruger,1857~1928)于 1912年对投影公式加以补充,故名。该投影按照投影带中央子午线投影为直线且长度不变和赤道投影为直线的条件,确定函数的形式,从而得到高斯一克吕格投影公式。投影后,除中央子午线和赤道为直线外, 其他子午线均为对称于中央子午线的曲线。设想用一个椭圆柱横切于椭球面上投影带的中央子午线,按上述投影条件,将中央子午线两侧一定经差范围内的椭球面正形投影于椭圆柱面。将椭圆柱面沿过南北极的母线剪开展平,即为高斯投影平面。取中央子午线与赤道交点的投影为原点,中央子午线的投影为纵坐标x轴,赤道的投影为横坐标y轴,构成高斯克吕格平面直角坐标系。
高斯-克吕格投影在长度和面积上变形很小,中央经线无变形,自中央经线向投影带边缘,变形逐渐增加,变形最大之处在投影带内赤道的两端。由于其投影精度高,变形小,而且计算简便(各投影带坐标一致,只要算出一个带的数据,其他各带都能应用),因此在大比例尺地形图中应用,可以满足军事上各种需要,能在图上进行精确的量测计算。
(2)高斯-克吕格投影分带
按一定经差将地球椭球面划分成若干投影带,这是高斯投影中限制长度变形的最有效方法。分带时既要控制长度变形使其不大于测图误差,又要使带数不致过多以减少换带计算工作,据此原则将地球椭球面沿子午线划分成经差相等的瓜瓣形地带,以便分带投影。通常按经差6度或3度分为六度带或三度带。六度带自0度子午线起每隔经差6度自西向东分带,带号依次编为第 1、2…60带。三度带是在六度带的基础上分成的,它的中央子午线与六度带的中央子午线和分带子午线重合,即自 1.5度子午线起每隔经差3度自西向东分带,带号依次编为三度带第 1、2…120带。我国的经度范围西起 73°东至135°,可分成六度带十一个,各带中央经线依次为75°、81°、87°、……、117°、123°、129°、135°,或三度带二十二个。六度带可用于中小比例尺(如 1:250000)测图,三度带可用于大比例尺(如 1:10000)测图,城建坐标多采用三度带的高斯投影。
(3)高斯-克吕格投影坐标系
高斯- 克吕格投影是按分带方法各自进行投影,故各带坐标成独立系统。以中央经线投影为纵轴X, 赤道投影为横轴Y,两轴交点即为各带的坐标原点。纵坐标以赤道为零起算,赤道以北为正,以南为负。我国位于北半球,纵坐标均为正值。横坐标如以中央经线为零起算,中央经线以东为正,以西为负,横坐标出现负值,使用不便,故规定将坐标纵轴西移500公里当作起始轴,凡是带内的横坐标值均加 500公里。由于高斯-克吕格投影每一个投影带的坐标都是对本带坐标原点的相对值,所以各带的坐标完全相同,为了区别某一坐标系统属于哪一带,在横轴坐标前加上带号,如(4231898m,21655933m),其中21即为带号。
(4)高斯-克吕格投影与UTM投影
高斯-克吕格(Gauss-Kruger)投影与UTM投影(Universal Transverse Mercator,通用横轴墨卡托投影)都是横轴墨卡托投影的变种,目前一些国外的软件或国外进口仪器的配套软件往往不支持高斯-克吕格投影,但支持UTM投影,因此常有把UTM投影当作高斯-克吕格投影的现象。
从投影几何方式看,高斯-克吕格投影是“等角横切圆柱投影”,投影后中央经线保持长度不变,即比例系数为1;UTM投影是“等角横轴割圆柱投影”,圆柱割地球于南纬80度、北纬84度两条等高圈,投影后两条割线上没有变形,中央经线上长度比0.9996。从计算结果看,两者主要差别在比例因子上,高斯-克吕格投影中央经线上的比例系数为1, UTM投影为0.9996,高斯-克吕格投影与UTM投影可近似采用 X[UTM]=0.9996 * X[高斯],Y[UTM]=0.9996 * Y[高斯],进行坐标转换(见下表)。从分带方式看,两者的分带起点不同,高斯-克吕格投影自0度子午线起每隔经差6度自西向东分带,第1带的中央经度为3°;UTM投影自西经180°起每隔经差6度自西向东分带,第1带的中央经度为-177°,因此高斯-克吕格投影的第1带是UTM的第31带。此外,两投影的东伪偏移都是500公里,高斯-克吕格投影北伪偏移为零,UTM北半球投影北伪偏移为零,南半球则为10000公里。
高斯-克吕格投影与UTM投影可近似采用 Xutm=0.9996 * X高斯,Yutm=0.9996 * Y高斯进行坐标转换。以下举例说明(基准面为WGS84,中央经度为123°):
| 输入坐标(度) | 高斯投影(米) | UTM投影(米) |
Xutm=0.9996 * X高斯, Yutm=0.9996 * Y高斯 |
|
| 纬度值(X) | 32 | 3543600.9 | 3542183.5 | 3543600.9*0.9996 ≈ 3542183.5 |
| 经度值(Y) | 121 | 21310996.8 | 311072.4 | (310996.8-500000)*0.9996+500000 ≈ 311072.4 |
注:坐标点(32,121)位于高斯投影的21带,高斯投影Y值21310996.8中前两位“21”为带号;坐标点(32,121)位于UTM投影的51带,上表中UTM投影的Y值没加带号。因坐标纵轴西移了500000米,转换时必须将Y值减去500000乘上比例因子后再加500000。
单点转换步骤如下:
(1)选择是高斯正转换还是反转换,缺省为经纬度转换到高斯投影坐标,投影坐标单位为米。
(2)选择大地基准面,缺省北京54,如果是GPS定位数据别忘了切换为WGS84。
(3)选择分带,3度或6度, 缺省为6度。
(4)输入中央经度,20带(114°E~120°E)中央经度为117度,21带(120°E~126°E)中央经度为123度。
(5)如正向投影,选择经纬度输入数据格式,有三个选项,缺省为十进制度格式。具体输入方式如下例:
| 格 式 | 原始纬度值 | 原始经度值 | 输入纬度值 | 输入经度值 |
| 十进制度 | 35.445901° | 122.997344° | 35.445901 | 122.997344 |
| 度分 | 35°26.7541′ | 122°59.8406′ | 3526.7541 | 12259.8406 |
| 度分秒 | 35°26′45.245″ |
122°59′50.438″ |
352645.245 |
1225950.438 |
(6)正投影按选定格式在“输入”栏输入经纬度值,反投影输入以米为单位的X、Y坐标值。
(7)单击“单点转换”按钮。
(8)在“输出”栏查看计算结果。
批量转换步骤如下:
(1)准备好需要转换的输入数据文件,要求是文本文件,分两列,第一列纬度值或纵向坐标值,第二列经度值或横向坐标值,两列之间用分隔符分开。正向投影时,纬度值及经度值格式可以有三种选择(见表),缺省当作十进制度处理;反向投影时,纵向及横向坐标值必须以米为单位。
下例为度分秒格式(WGS84)的6°带正投影输入数据文件 testdata.txt
352645.245 1225950.438
353800.402 1230000.378
351600.519 1225959.506
345800.101 1225959.8
343600.336 1230000.26
341400.018 1225959.897
335159.17 1225959.46
333000.08 1230000.28
(2)选择是高斯正转换还是反转换,缺省为经纬度转换到高斯投影坐标,投影坐标单位为米。
(3)选择大地基准面,缺省为WGS84。
(4)选择分带,3度或6度, 缺省为6度。
(5)输入中央经度,20带(114°E~120°E)中央经度为117度,21带(120°E~126°E)中央经度为123度。
(6)如正向投影,选择输入数据文件中的经纬度输入数据格式,有三个选项,缺省为十进制度格式。
(7)单击“批量转换”按钮。弹出打开文件对话框,输入你的数据文件名。
(8)输入转换结果文件名,单击“保存”后,程序开始进行计算。
(9)打开输出文件查看计算结果,结果分五列,第一序号,第二列输入纬度值或纵向坐标值,第三列输入经度值或横向坐标值,第四列转换后纬度值或纵向坐标值,第五列转换后经度值或横向坐标值。
下例为度分秒格式(WGS84)的6°带,中央经度123°正投影转换结果数据文件 result.txt
1
352645.245 1225950.438 3924063.29 21499758.85
2 353800.402 1230000.378 3944871.34
21500009.51
3 351600.519 1225959.506 3904193.7
21499987.51
4 345800.101 1225959.8
3870898.01 21499994.93
5 343600.336 1230000.26 3830228.49
21500006.62
6 341400.018 1225959.897 3789544.38
21499997.36
7 335159.17 1225959.46
3748846.38 21499986.12
8 333000.08 1230000.28
3708204.97 21500007.23