地图投影 - Rover's Official Blog

地图投影

快两个星期没写博客了，也许有点懒，也有点小忙。前阵子帮师弟做了个毕业设计，帮朋友做了个网站，最近工作上又比较忙，而且也一直没有做什么研究，所以实在想不起博客上写点什么。
想研究投影的，但也一直没怎么搞的懂，算法就更不用讲了，一切只能按我头脑的想的，但也整理不出什么思想来。真是惭愧那，做地图的连投影都不是清楚。写这篇，最多算存点投影的资料吧。有兴趣的朋友可以研究研究。如果喜欢，当然可以把下面的评论收藏，呵呵。至于提问，就免了，我也回答不上来。愿意的话多讨论讨论。
推荐园子里的一位兄弟JETZ，写了十多篇的关于坐标投影的问题，可以看他的随笔分类：GIS-坐标
下面评论内容如下：
1.MapInfo地图投影的添加，吐鲁番地区的坐标系－李风云的毕业论文，写的不错。
2.MapInfo中如何使用地理坐标系－网络上搜索最多的文章，讲了MapX中的坐标系定义与转换。
3.地理坐标系统简介－真的是简介，简单的介绍了地球椭球体，地图投影，高斯-克吕格直角坐标，我国地形图分幅与编号
4.高斯克吕格坐标系中国部分定义－MAPINFO中北京54和西安80的两种坐标系的定义(在MAPINFO8.0中已经有定义，不需要自行添加)
5.高斯投影变换－C的代码，看不懂
6.经纬度BL换算到高斯平面直角坐标XY（高斯投影正算）的源码及算法－VB代码，算法和代码结合还是看不大懂。
PS：以上文章均出自James Liu's MapInfo论坛，地址可以参看这里：http://www.mygis.com.cn/forum/index.asp

posted on 2006-04-17 22:48 Rover.Tang 阅读(11100) 评论(11) 编辑收藏网摘所属分类: Map

#1楼 [楼主] 2006-04-17 22:49 浪人|努力

MapInfo中如何使用地理坐标系
这个图是采用，大地坐标的经纬度WGS84标准。与GPS采用的WGS84在同一点上相差经纬度相差有10－20分（仅仅从地图在MAPINFO的显示数据来说没涉及到NumericCoordSys的坐标系统）。

MapX中的坐标系定义与转换
GIS中的坐标系定义是GIS系统的基础，正确定义GIS系统的坐标系非常重要。
1. 椭球体、基准面及地图投影
GIS中的坐标系定义由基准面和地图投影两组参数确定，而基准面的定义则由特定椭球体及其对应的转换参数确定，因此欲正确定义GIS系统坐标系，首先必须弄清地球椭球体(Ellipsoid)、大地基准面(Datum)及地图投影(Projection)三者的基本概念及它们之间的关系。
基准面是利用特定椭球体对特定地区地球表面的逼近，因此每个国家或地区均有各自的基准面，我们通常称谓的北京54坐标系、西安80坐标系实际上指的是我国的两个大地基准面。我国参照前苏联从1953年起采用克拉索夫斯基(Krassovsky)椭球体建立了我国的北京54坐标系，1978年采用国际大地测量协会推荐的1975地球椭球体建立了我国新的大地坐标系--西安80坐标系，目前大地测量基本上仍以北京54坐标系作为参照，北京54与西安80坐标之间的转换可查阅国家测绘局公布的对照表。 WGS1984基准面采用WGS84椭球体，它是一地心坐标系，即以地心作为椭球体中心，目前GPS测量数据多以WGS1984为基准。
上述3个椭球体参数如下：
椭球体 Mapinfo中代号年代长半轴短半轴 1/扁率
Krassovsky 3 1940 6378245 6356863 298.3
IAG 75 31 1975 6378140 6356755 298.25722101
WGS 84 28 1984 6378137.000 6356752.314 298.257223563
椭球体与基准面之间的关系是一对多的关系，也就是基准面是在椭球体基础上建立的，但椭球体不能代表基准面，同样的椭球体能定义不同的基准面，如前苏联的Pulkovo 1942、非洲索马里的Afgooye基准面都采用了Krassovsky椭球体，但它们的基准面显然是不同的。
地图投影是将地图从球面转换到平面的数学变换，如果有人说：该点北京54坐标值为X=4231898,Y=21655933，实际上指的是北京54基准面下的投影坐标，也就是北京54基准面下的经纬度坐标在直角平面坐标上的投影结果。
2. GIS中基准面的定义与转换
虽然现有GIS平台中都预定义有上百个基准面供用户选用，但均没有我们国家的基准面定义。假如精度要求不高，可利用前苏联的Pulkovo 1942基准面(Mapinfo中代号为1001)代替北京54坐标系；假如精度要求较高，如土地利用、海域使用、城市基建等GIS系统，则需要自定义基准面。
GIS系统中的基准面通过当地基准面向WGS1984的转换7参数来定义，转换通过相似变换方法实现，具体算法可参考科学出版社1999年出版的《城市地理信息系统标准化指南》第76至86页。假设Xg、Yg、Zg表示WGS84地心坐标系的三坐标轴，Xt、Yt、Zt表示当地坐标系的三坐标轴，那么自定义基准面的7参数分别为：三个平移参数ΔX、ΔY、ΔZ表示两坐标原点的平移值；三个旋转参数εx、εy、εz表示当地坐标系旋转至与地心坐标系平行时，分别绕Xt、Yt、Zt的旋转角；最后是比例校正因子，用于调整椭球大小。
MapX中基准面定义方法如下：
Datum.Set(Ellipsoid, ShiftX, ShiftY, ShiftZ, RotateX, RotateY, RotateZ, ScaleAdjust, PrimeMeridian)
其中参数： Ellipsoid为基准面采用的椭球体；
ShiftX, ShiftY, ShiftZ为平移参数；
RotateX, RotateY, RotateZ为旋转参数；
ScaleAdjust为比例校正因子，以百万分之一计；
PrimeMeridian为本初子午线经度，在我国取0，表示经度从格林威治起算。
美国国家测绘局(National Imagery and Mapping Agency)公布了世界大多数国家的当地基准面至WGS1984基准面的转换3参数(平移参数)，可从 http://164.214.2.59/GandG/wgs84dt/dtp.html 下载，其中包括有香港Hong Kong 1963基准面、台湾 Hu-Tzu-Shan 基准面的转换3参数，但是没有中国大陆的参数。
实际工作中一般都根据工作区内已知的北京54坐标控制点计算转换参数，如果工作区内有足够多的已知北京54与WGS84坐标控制点，可直接计算坐标转换的7参数或3参数；当工作区内有3个已知北京54与WGS84坐标控制点时，可用下式计算WGS84到北京54坐标的转换参数(A、B、C、D、E、F)：x54 = AX84 + BY84 + C，y54 = DX84 + EY84 + F，多余一点用作检验；在只有一个已知控制点的情况下(往往如此)，用已知点的北京54与WGS84坐标之差作为平移参数，当工作区范围不大时精度也足够了。
从Mapinfo中国的URL(http://www.mapinfo.com.cn/download)可下载到包含北京54、西安80坐标系定义的Mapinfow.prj文件，其中定义的北京54基准面参数为：(3,24,-123,-94,-0.02,0.25,0.13,1.1,0)，西安80基准面参数为：(31,24,-123,-94,-0.02,0.25,0.13,1.1,0)，文件中没有注明其参数的来源，我发现它们与Mapinfo参考手册附录G"定义自定义基准面"中的一个例子所列参数相同，因此其可靠性值得怀疑，尤其从西安80与北京54采用相同的7参数来看，至少西安80的基准面定义肯定是不对的。因此，当系统精度要求较高时，一定要对所采用的参数进行检测、验证，确保坐标系定义的正确性。
3. GIS中地图投影的定义
我国的基本比例尺地形图(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)，我国的GIS系统中应该采用与我国基本比例尺地形图系列一致的地图投影系统。
在MapX中坐标系定义由基准面、投影两部分参数组成，方法如下：
CoordSys.Set(Type, [Datum], [Units], [OriginLongitude], [OriginLatitude],
[StandardParallelOne], [StandardParallelTwo], [Azimuth], [ScaleFactor],
[FalseEasting], [FalseNorthing], [Range], [Bounds], [AffineTransform])
其中参数：Type表示投影类型，Type为1时地图坐标以经纬度表示，它是必选参数，它后面的参数都为可选参数；
Datum为大地基准面对象，如果采用非地球坐标(NonEarth)无需定义该参数；
Units为坐标单位，如Units为7表示以米为单位；
OriginLongitude、OriginLatitude分别为原点经度和纬度；
StandardParallelOne、StandardParallelTwo为第一、第二标准纬线；
Azimuth为方位角，斜轴投影需要定义该参数；
ScaleFactor为比例系数；
FalseEasting, FalseNorthing为东伪偏移、北伪偏移值；
Range为地图可见纬度范围；
Bounds为地图坐标范围，是一矩形对象，非地球坐标(NonEarth)必须定义该参数；
AffineTransform为坐标系变换对象。
相应高斯-克吕格投影、兰勃特投影、墨卡托投影需要定义的坐标系参数序列如下：
高斯-克吕格：投影代号(Type)，基准面(Datum)，单位(Unit)，
中央经度(OriginLongitude)，原点纬度(OriginLatitude)，
比例系数(ScaleFactor)，
东伪偏移(FalseEasting)，北纬偏移(FalseNorthing)
兰勃特: 投影代号(Type)，基准面(Datum)，单位(Unit)，
中央经度(OriginLongitude)，原点纬度(OriginLatitude)，
标准纬度1(StandardParallelOne)，标准纬度2(StandardParallelTwo)，
东伪偏移(FalseEasting)，北纬偏移(FalseNorthing)
墨卡托: 投影代号(Type)，基准面(Datum)，单位(Unit)，
原点经度(OriginLongitude)，原点纬度(OriginLatitude)，
标准纬度(StandardParallelOne)
在城市GIS系统中均采用6度或3度分带的高斯-克吕格投影，因为一般城建坐标采用的是6度或3度分带的高斯-克吕格投影坐标。高斯-克吕格投影以6度或3度分带，每一个分带构成一个独立的平面直角坐标网，投影带中央经线投影后的直线为X轴(纵轴，纬度方向)，赤道投影后为Y轴(横轴，经度方向)，为了防止经度方向的坐标出现负值，规定每带的中央经线西移500公里，即东伪偏移值为500公里，由于高斯-克吕格投影每一个投影带的坐标都是对本带坐标原点的相对值，所以各带的坐标完全相同，因此规定在横轴坐标前加上带号，如(4231898,21655933)其中21即为带号，同样所定义的东伪偏移值也需要加上带号，如21带的东伪偏移值为21500000米。
假如你的工作区位于21带，即经度在120度至126度范围，该带的中央经度为123度，采用Pulkovo 1942基准面，那么定义6度分带的高斯-克吕格投影坐标系参数为：(8，1001，7，123，0，1，21500000，0)。
那么当精度要求较高，实测数据为WGS1984坐标数据时，欲转换到北京54基准面的高斯-克吕格投影坐标，如何定义坐标系参数呢？你可选择WGS 1984(Mapinfo中代号104)作为基准面，当只有一个已知控制点时(见第2部分)，根据平移参数调整东伪偏移、北纬偏移值实现WGS84到北京54的转换，如:(8，104，7，123，0，1，21500200，-200)，也可利用 AffineTransform坐标系变换对象,此时的转换系数(A、B、C、D、E、F)中A、B、D、E为0，只有X、Y方向的平移值C、F ；当有3个已知控制点时，可利用得到的转换系数(A、B、C、D、E、F)定义 AffineTransform坐标系变换对象，实现坐标系的转换，如：(8，104，7，123，0，1，21500000，0，map.AffineTransform)，其中AffineTransform定义为AffineTransform.set(7，A、B、C、D、E、F)(7表示单位米)；当然有足够多已知控制点时，直接求定7参数自定义基准面就行了。

回复引用查看

#2楼 [楼主] 2006-04-17 22:49 浪人|努力

MapInfo地图投影的添加 ——吐鲁番地区的坐标系
李风云 981207014 新疆大学资源环境学院地信98-1班
内容提要：GIS中的坐标系定义是GIS系统的基础，正确定义GIS系统的坐标系非常重要。在MapInfo地图投影中缺少针对中国很多地区大比例尺和中比例尺的投影方式，通过对MapInfo中Mapinfow.prj文件的研究，以吐鲁番地区为例，添加了高斯-克吕格（Gauss-Krüger）投影，能解决一些大比例尺对变形小的要求。
关键字：地图投影坐标系 MapInfo
一、前言
地球高低不平、极其复杂的自然表面。为研究和工作方便，常将地球近似为为一个旋转椭球体，称为地球椭球体。地球椭球体的表面是一个不可展的曲面，地图以平面方式表示地球表面(全部或一部分)。将地球椭球体上的点的坐标投影到平面坐标的方法称为地图投影。地图投影的种类繁多，不同的投影方式具有不同的形态和变形特征。根据不同的使用目的，可以采用不同的投影方式，一种投影对一种目的是有用的，而对另一种目的则可能不适合。
地图投影的选取决定于地图的应用及其比例尺大小，在桌面数字环境下，用户可对每一个新创建文件进行投影设置和选择。MapInfo以两级目录菜单的形式提供了300多个预定义坐标系，当用户要使用其他坐标系或创建新的坐标系时，还可以通过修改投影参数文件（Mapinfow.prj）来实现。这个数据文件以分行形式记录每一个预定义坐标系的参数表，如坐标系名称、投影代码、基准面代码、坐标单位、原点经度、原点纬度、标准纬线1、标准纬线2、方位角、比例系数等。MapInfo系统定义了50个基准面代码，它们包括了全世界和一些国家采用的基准面。如果用户要采用其他的基准面，并且知道该基准面的数学参数，则可以使用代码在该投影参数文件中定义这个基准面。
虽然现有GIS平台中都预定义有上百个基准面供用户选用，但均没有我们国家的基准面定义。假如精度要求不高，可利用前苏联的Pulkovo 1942基准面(Mapinfo中代号为1001)代替北京54坐标系；假如精度要求较高，如土地利用、海域使用、城市基建等GIS系统，则需要自定义基准面。
二、MAPINFOW.PRJ文件的参数分析
MAPINFOW.PRJ在MapInfo软件的安装目录下，可以用记事本打开文件，它的格式如下：
"--- Longitude / Latitude ---", 0, 0, 0, 0., 0., 0., 0., 0., 0., 0., 0.
"Longitude / Latitude", 1, 0
"Longitude / Latitude (Adindan)p4201", 1, 1
"Longitude / Latitude (Afgooye)p4205", 1, 2
"Longitude / Latitude (AGD 66)p4202", 1, 12
"Longitude / Latitude (AGD 84)p4203", 1, 13
"Longitude / Latitude (Ain el Abd 1970)p4204", 1, 3
"Longitude / Latitude (Anna 1 Astro 1965)", 1, 4
"Longitude / Latitude (Arc 1950)p4209", 1, 5
第一行是投影方式的一级菜单，后面是二级菜单。每行的第一部分是位于引号内的坐标系名称。每一行的第二部分是指定投影的代号。行中其余部分是有关该特定坐标系的参数值。下面一以中国区域等积投影为例进行分析，详细分析各参数的含义。Regional Equal-Area Projections 在MAPINFOW.PRJ文件中的代码如下：
"Equal-Area Projection (China)", 9, 0, 0, 110, 10, 25, 40, 0, 0
表 1
投影正轴等面积割圆锥
基准面 GRS 80基准面
单位英里
原点 110°E, 10°N
标准纬线 25°N、 40°N
东伪偏移 0
北伪偏移 0
第一个参数 3 代表投影方式是Albers Equal Area Conic（正轴等面积割圆锥投影，也叫亚尔勃斯投影）；第二个参数 0 代表GRS 80基准面；第三个参数 0 代表单位是英里；第四和第五个是代表原点坐标是 110°E, 10°N；第六和第七个代表两条标准纬线分别为25°N、 40°N；最后两个参数是东伪偏移和北伪偏移。Equal-Area Projection (China)的列表（表1）。由于参数的数目繁多，每种参数在《MapInfo参考手册》中都有详细的表格。
高斯-克吕格（Gauss-Krüger）投影在MAPINFOW.PRJ文件需要定义的坐标系参数序列如下：
高斯-克吕格：投影代号(Type)，基准面(Datum)，单位(Unit)，
中央经度(OriginLongitude)，原点纬度(OriginLatitude)，
比例系数(ScaleFactor)，
东伪偏移(FalseEasting)，北纬偏移(FalseNorthing)
三、对MAPINFOW.PRJ文件添加坐标系
3.1我国大地坐标系的参考椭球体与基准面
基准面是利用特定椭球体对特定地区地球表面的逼近，因此每个国家或地区均有各自的基准面，我们通常称谓的北京54坐标系、西安80坐标系实际上指的是我国的两个大地基准面。
我国参照前苏联从1953年起采用克拉索夫斯基(Krassovsky)椭球体建立了我国的北京54坐标系，1978年采用国际大地测量协会推荐的1975地球椭球体建立了我国新的大地坐标系--西安80坐标系，目前大地测量基本上仍以北京54坐标系作为参照，北京54与西安80坐标之间的转换可查阅国家测绘局公布的对照表。 WGS1984基准面采用WGS84椭球体，它是一地心坐标系，即以地心作为椭球体中心，目前GPS测量数据多以WGS1984为基准。
上述3个椭球体在MAPINFOW.PRJ文件中参数如下（表 2）：
椭球体代号年代长半轴短半轴 1/扁率
Krassovsky 椭球 3 1940 6378245 6356863 298.3
IAG 75椭球 31 1975 6378140 6356755 298.25722101
WGS 84 椭球 28 1984 6378137.000 6356752.314 298.257223563
表 2
椭球体与基准面之间的关系是一对多的关系，也就是基准面是在椭球体基础上建立的，但椭球体不能代表基准面，同样的椭球体能定义不同的基准面。
虽然现有GIS平台中都预定义有上百个基准面供用户选用，但均没有我们国家的基准面定义。假如精度要求不高，可利用前苏联的Pulkovo 1942基准面(Mapinfo中代号为1001)代替北京54坐标系；假如精度要求较高，如土地利用、城市基建等GIS系统，则需要自定义基准面。
从Mapinfo中国的URL(http://www.mapinfo.com.cn)可下载到包含北京54、西安80坐标系，其中定义：
北京54基准面参数为：(3,24,-123,-94,-0.02,0.25,0.13,1.1,0)
西安80基准面参数为：(31,24,-123,-94,-0.02,0.25,0.13,1.1,0)
这样就可以自定义一个基准面参数为“9999”的基准面。当测数据为WGS84坐标数据时，由于北京54和西安80坐标系中有8个参数相同，不同的是参考椭球体参数，由此类推WGS84的基准面参数定义为：
(28,24,-123,-94,-0.02,0.25,0.13,1.1,0)
这样就可以以我国常用的基准面建立坐标系，也可以选择不同的投影方式，以满足系统对精度的要求。
3.2高斯-克吕格（Gauss-Krüger）投影
我国的基本比例尺地形图(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)，我国的GIS系统中应该采用与我国基本比例尺地形图系列一致的地图投影系统。
高斯-克吕格（Gauss-Krüger）投影又叫横轴墨卡托(Transverse Mercator)投影，依椭圆柱作为投影面，并与椭球体相切于一条经线上，该经线即为投影带的中央经线，按等角条件将中央经线东西一定范围内的区域投影到椭圆柱表面上，再展成平面，便构成了横轴等角切椭圆柱投影。
在GIS系统中均采用6°或3°分带的高斯-克吕格投影，因为一般坐标采用的是6°或3°分带的高斯-克吕格投影坐标。高斯-克吕格投影以6°或3°分带，每一个分带构成一个独立的平面直角坐标网，投影带中央经线投影后的直线为X轴(纵轴，纬度方向)，赤道投影后为Y轴(横轴，经度方向)，为了防止经度方向的坐标出现负值，规定每带的中央经线西移500公里，即东伪偏移值为500公里，由于高斯-克吕格投影每一个投影带的坐标都是对本带坐标原点的相对值，所以各带的坐标完全相同，因此规定在横轴坐标前加上带号，如(4231898,21655933)其中21即为带号，同样所定义的东伪偏移值也需要加上带号，如21带的东伪偏移值为21500000米。假如你的工作区位于21带，即经度在120°至126°范围，该带的中央经度为123°。
3.3吐鲁番地区的分析
吐鲁番地区位于中国西北部的新疆境内，经度差为2.7°，经度位于87°20’ E 和92°E之间。可以采用3°分带的高斯-克吕格投影，其中央经线为90°。
3.4 创建新的坐标系
前面以对各参数做了分析，针对吐鲁番地区建立坐标系，要对MAPINFOW.PRJ文件增加坐标系的下列参数（表3）：
投影基准面单位中央经度原点纬度比例系数东伪偏移北纬偏移
8 9999 7 90 0 1 30500000 0
表 3
自定义基准面9999为：
北京54基准面参数为：(3,24,-123,-94,-0.02,0.25,0.13,1.1,0)
西安80基准面参数为：(31,24,-123,-94,-0.02,0.25,0.13,1.1,0)
WGS84基准面参数为：(28,24,-123,-94,-0.02,0.25,0.13,1.1,0)
需要向MAPINFOW.PRJ文件增加一个列，要有相应要素的新条目。该过程描述如下：
1 在一个文本编辑器或字处理器中打开MAPINFOW.PRJ文件。
2输入代表新坐标系的下列参数：
"--- 吐鲁番坐标系统---"
"Gauss-Kruger 30 (西安1980 3°分带)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 90, 0, 1, 30500000, 0
"Gauss-Kruger 30 (北京1954 3°分带)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 90, 0, 1, 30500000, 0
"Gauss-Kruger 30 (WGS1984 3°分带)", 8, 9999, 28, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 90, 0, 1, 30500000, 0
3 如果有必要把该条目移到相似坐标系之间的合适位置。
4 保存编辑过的MAPINFOW.PRJ文件。
这样就可以使用自定义的坐标系了，就像MapInfo中提供的坐标系一样。
四：结束语
地图投影在GIS领域中占有重要的地位，通过吐鲁番地区高斯-克吕格（Gauss-Krüger）投影坐标系的定义，也可以定义一些大比例尺的UTM投影（Universal Transverse Mercatol Projection，即通用横轴墨卡托投影）、等积割圆锥投影等等，提高系统的投影精度。坐标系是GIS系统的基石，正确设置坐标系是系统成败的关键

回复引用查看

#3楼 [楼主] 2006-04-17 22:50 浪人|努力

地理坐标系统简介
地理坐标系，也可称为真实世界的坐标系，是用于确定地物在地球上位置的坐标系。一个特定的地理坐标系是由一个特定的椭球体和一种特定的地图投影构成，其中椭球体是一种对地球形状的数学描述，而地图投影是将球面坐标转换成平面坐标的数学方法。绝大多数的地图都是遵照一种已知的地理坐标系来显示坐标数据。
1．地球椭球体
地球是一个表面很复杂的球体，人们以假想的平均静止的海水面形成的“大地体”为参照，推求出近似的椭球体，理论和实践证明，该椭球体近似一个以地球短轴为轴的椭园而旋转的椭球面，这个椭球面可用数学公式表达，将自然表面上的点归化到这个椭球面上，就可以计算了。下面列举了一些常用的一些椭球及参数：
1）海福特椭球(1910)　　　　　　　　　　　　　　　　　　我国52年以前采用的椭球
　　a=6378388m b=6356911.9461279m α=0.33670033670
2）克拉索夫斯基椭球(1940 Krassovsky)　　　　　　　　　　北京54坐标系采用的椭球
　　a=6378245m b=6356863.018773m α=0.33523298692
3）1975年I.U.G.G推荐椭球(国际大地测量协会1975) 　　　西安80坐标系采用的椭球
　　a=6378140m b=6356755.2881575m α=0.0033528131778
4）WGS-84椭球(GPS全球定位系统椭球、17届国际大地测量协会) WGS-84坐标系椭球
　　a=6378137m b=6356752.3142451m α=0.00335281006247
最常用的地理坐标系是经纬度坐标系，这个坐标系可以确定地球上任何一点的位置，如果我们将地球看作一个椭球体，而经纬网就是加在地球表面的地理坐标参照系格网，经度和纬度是从地球中心对地球表面给定点量测得到的角度，经度是东西方向，而纬度是南北方向，经线从地球南北极穿过，而纬线是平行于赤道的环线。地理坐标可分为天文地理坐标和大地地理坐标：天文地理坐标是用天文测量方法确定的，大地地理坐标是用大地测量方法确定的。我们在地球椭球面上所用的地理坐标系属于大地地理坐标系，简称大地坐标系。
确定椭球的大小后，还要进行椭球定向，即把旋转椭球面套在地球的一个适当的位置，这一位置就是该地理坐标系的“坐标原点”，是全部大地坐标计算的起算点，俗称“大地原点”。
需要说明的是经纬度坐标系不是一种平面坐标系，因为度不是标准的长度单位，不可用其量测面积长度；平面坐标系(又称笛卡儿坐标系)，因其具有以下特性：可量测水平X方向和竖直Y方向的距离，可进行长度、角度和面积的量测，可用不同的数学公式将地球球体表面投影到二维平面上而得到广泛的应用。而每一个平面坐标系都有一特定的地图投影方法。
2．地图投影
是为解决由不可展的椭球面描绘到平面上的矛盾，用几何透视方法或数学分析的方法，将地球上的点和线投影到可展的曲面(平面、园柱面或圆锥面)上，将此可展曲面展成平面，建立该平面上的点、线和地球椭球面上的点、线的对应关系。
地图投影的过程是可以想象用一张足够大的纸去包裹地球，将地球上的地物投射到这张纸上。地球表面投影到平面上、圆锥面或者圆柱面上，然后把圆锥面、圆柱面沿母线切开后展成平面。根据这张纸包裹的方式，地图投影又可以分成：方位投影、圆锥投影和圆柱投影。根据这张纸与地球相交的方式，地图投影又可以分成切投影和割投影，在切线或者割线上的地物是没有变形的，而距离切线或者割线越远变形越大。
还有不少投影直接用解析法得到。根据所借助的几何面不同可分为伪方位投影、伪圆锥投影、伪圆柱投影等。
地图投影会存在两种误差，形状变化(也称角度变化)或者面积变化。投影以后能保持形状不变化的投影，称为等角投影 (Conformal mapping)，它的优点除了地物形状保持不变以外，在地图上测量两个地物之间的角度也能和实地保持一致，这非常重要，当在两地间航行必须保持航向的准确；或者另外一个例子是无论长距离发射导弹还是短距离发射炮弹，发射角度必须准确测量出来。因此等角投影是最常被使用的投影。等角投影的缺点是高纬度地区地物的面积会被放大。投影以后能保持形状不变化的投影，称为等面积投影 (Equivalent mapping)，在有按面积分析需要的应用中很重要，显示出来的地物相对面积比例准确，但是形状会有变化，假设地球上有个圆，投影后绘制出来即变成个椭圆了。还有第三种投影，非等角等面积投影，意思是既有形状变化也有面积变化，这类投影既不等角也不等积，长度、角度、面积都有变形。其中有些投影在某个主方向上保持长度比例等于1，称为等距投影。
每一种投影都有其各自的适用方面。例如,墨卡托投影适用于海图，其面积变形随着纬度的增高而加大，但其方向变形很小；横轴墨卡托投影的面积变形随着距中央经线的距离的加大而增大，适用于制作不同的国家地图。等角投影常用于航海图、风向图、洋流图等。现在世界各国地形图采用此类投影比较多。等积投影用于绘制经济地区图和某些自然地图。对于大多数数学地图和小比例尺普通地图来说，应优先考虑等积的要求。地理区域，诸如国家、水域和地理分类地区（植被、人口、气候等）相对分布范围，显然是十分重要的内容。任意投影常用作数学地图，以及要求沿某一主方向保持距离正确的地图。常用作世界地图的投影有墨卡托投影、高尔投影、摩尔威特投影、等差分纬线多圆锥投影、格灵顿投影、桑森投影、乌尔马耶夫投影等。下面对我国地形图所采用的高斯克吕格投影进行简单的介绍。
2.1高斯-克吕格直角坐标
　　高斯－克吕格投影（Gauss_Krivger）属于等角横切椭圆柱投影，是设想用一个椭圆柱横套在地球椭球的外面，并与设定的中央经线相切。其经纬线互相垂直，变形最大位于赤道与投影带最外一条经线的交点上，常用于纬度较高地区。
　　高斯－克吕格投影分带规定：该投影是我国国家基本比例尺地形图的数学基础，为控制变形，采用分带投影的方法，在比例尺 1：2.5万-1：50万图上采用6°分带，对比例尺为 1：1万及大于1：1万的图采用3°分带。
　　6°分带法：从格林威治零度经线起，每6°分为一个投影带，全球共分为60个投影带，东半球从东经0°-6°为第一带，中央经线为3°，依此类推，投影带号为1-30。其投影代号n和中央经线经度L0的计算公式为：L0=(6n-3)°；西半球投影带从180°回算到0°，编号为31-60，投影代号n和中央经线经度L0的计算公式为L0=360-(6n-3)°。
　　3°分带法：从东经1°30′起，每3°为一带，将全球划分为120个投影带，东经1°30′-4°30′，...178°30′-西经178°30′，...1°30′-东经1°30′。
　　东半球有60个投影带，编号1-60，各带中央经线计算公式：L0=3°n ,中央经线为3°、6°...180°。
　　西半球有60个投影带，编号1-60，各带中央经线计算公式：L0=360°-3°n ,中央经线为西经177°、...3°、0°。
我国规定将各带纵坐标轴西移500公里，即将所有y值加上500公里，坐标值前再加各带带号。以18带为例，原坐标值为y=243353.5m，西移后为y=743353.5，加带号通用坐标为y=18743353.5 。
为了方便大家对不同比例尺的地形图检索，最后对我国地形图的分幅与编号规则进行简单的介绍。
3．我国地形图分幅与编号
　　我国基本比例尺地形图分幅与编号，以1：100万地形图为基础，延伸出1：50万、1：25万、1：10万，再以1：10万为基础，延伸出1：5万、1:2.5万及1：1万三种比例尺。
　　1：100万从赤道起向两极每纬差4°为一行，至88°，南北半球各分为22横列，依次编号A、B、... V；由精度180°西向东每6°一列，全球60列，以1-60表示，如海南所在1：100万图在第5行，第49列，其编号为 E-49 。
　　在1：100万图上，按经差3°纬差2°分成四幅1：50万地形图，编为A、B、C、D，如 E-49-A。按经差1°30′纬差1°分成16幅1：25万地形图，编为[1]、...[16]，如 E-49-[1]。按经差30′纬差20′分成144幅1：10万地形图，编为1、...144，如 E-49-1。即后三种比例尺各自独立地与1：100万地图的图号联系。
　　1：10万图上每经差15′纬差10′分成四幅1：5万地形图，编为A、B、C、D，如 E-49-1-A。　
　　1：5万图上每经差7′30″纬差5′分成四幅1：2.5万，编为1、2、3、4，如 E-49-1-A-1。
　　1：10万图上每经差3′45″纬差2′30″分成64幅1：1万地形图，编为(1)、...(64)，如E-49-1-A-(1)。
　　1：1万图上每经差1′52″纬差1′15″分成四幅1：5000地形图，编为a、b、c、d，如E-49-1-A-(1)-a。

回复引用查看

#4楼 [楼主] 2006-04-17 22:50 浪人|努力

高斯克吕格坐标系中国部分定义
就是常说的北京五四和西安80的两种坐标系的定义
加到坐标系文件中就可以了：）
"--- Gauss-Kruger (Xian1980 3-degree zone) ---"
"GK Zone 25 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 75, 0, 1, 25500000, 0
"GK Zone 26 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 78, 0, 1, 26500000, 0
"GK Zone 27 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 81, 0, 1, 27500000, 0
"GK Zone 28 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 84, 0, 1, 28500000, 0
"GK Zone 29 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 87, 0, 1, 29500000, 0
"GK Zone 30 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 90, 0, 1, 30500000, 0
"GK Zone 31 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 93, 0, 1, 31500000, 0
"GK Zone 32 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 96, 0, 1, 32500000, 0
"GK Zone 33 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 99, 0, 1, 33500000, 0
"GK Zone 34 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 102, 0, 1, 34500000, 0
"GK Zone 35 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 105, 0, 1, 35500000, 0
"GK Zone 36 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 108, 0, 1, 36500000, 0
"GK Zone 37 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 111, 0, 1, 37500000, 0
"GK Zone 38 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 114, 0, 1, 38500000, 0
"GK Zone 39 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 117, 0, 1, 39500000, 0
"GK Zone 40 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 120, 0, 1, 40500000, 0
"GK Zone 41 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 123, 0, 1, 41500000, 0
"GK Zone 42 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 126, 0, 1, 42500000, 0
"GK Zone 43 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 129, 0, 1, 43500000, 0
"GK Zone 44 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 132, 0, 1, 44500000, 0
"GK Zone 45 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 135, 0, 1, 45500000, 0

"--- Gauss-Kruger (Xian1980 6-degree zone) ---"
"GK Zone 13 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 75, 0, 1, 13500000, 0
"GK Zone 14 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 81, 0, 1, 14500000, 0
"GK Zone 15 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 87, 0, 1, 15500000, 0
"GK Zone 16 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 93, 0, 1, 16500000, 0
"GK Zone 17 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 99, 0, 1, 17500000, 0
"GK Zone 18 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 105, 0, 1, 18500000, 0
"GK Zone 19 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 111, 0, 1, 19500000, 0
"GK Zone 20 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 117, 0, 1, 20500000, 0
"GK Zone 21 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 123, 0, 1, 21500000, 0
"GK Zone 22 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 129, 0, 1, 22500000, 0
"GK Zone 23 (Xian1980)", 8, 9999, 31, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 135, 0, 1, 23500000, 0

"--- Gauss-Kruger (Beijing1954 3-degree zone) ---"
"GK Zone 25 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 75, 0, 1, 25500000, 0
"GK Zone 26 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 78, 0, 1, 26500000, 0
"GK Zone 27 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 81, 0, 1, 27500000, 0
"GK Zone 28 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 84, 0, 1, 28500000, 0
"GK Zone 29 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 87, 0, 1, 29500000, 0
"GK Zone 30 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 90, 0, 1, 30500000, 0
"GK Zone 31 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 93, 0, 1, 31500000, 0
"GK Zone 32 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 96, 0, 1, 32500000, 0
"GK Zone 33 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 99, 0, 1, 33500000, 0
"GK Zone 34 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 102, 0, 1, 34500000, 0
"GK Zone 35 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 105, 0, 1, 35500000, 0
"GK Zone 36 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 108, 0, 1, 36500000, 0
"GK Zone 37 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 111, 0, 1, 37500000, 0
"GK Zone 38 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 114, 0, 1, 38500000, 0
"GK Zone 39 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 117, 0, 1, 39500000, 0
"GK Zone 40 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 120, 0, 1, 40500000, 0
"GK Zone 41 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 123, 0, 1, 41500000, 0
"GK Zone 42 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 126, 0, 1, 42500000, 0
"GK Zone 43 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 129, 0, 1, 43500000, 0
"GK Zone 44 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 132, 0, 1, 44500000, 0
"GK Zone 45 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 135, 0, 1, 45500000, 0

"--- Gauss-Kruger (Beijing1954 6-degree zone) ---"
"GK Zone 13 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 75, 0, 1, 13500000, 0
"GK Zone 14 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 81, 0, 1, 14500000, 0
"GK Zone 15 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 87, 0, 1, 15500000, 0
"GK Zone 16 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 93, 0, 1, 16500000, 0
"GK Zone 17 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 99, 0, 1, 17500000, 0
"GK Zone 18 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 105, 0, 1, 18500000, 0
"GK Zone 19 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 111, 0, 1, 19500000, 0
"GK Zone 20 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 117, 0, 1, 20500000, 0
"GK Zone 21 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 123, 0, 1, 21500000, 0
"GK Zone 22 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 129, 0, 1, 22500000, 0
"GK Zone 23 (Beijing1954)", 8, 9999, 3, 24, -123, -94, -0.02, 0.25, 0.13, 1.1, 0, 7, 135, 0, 1, 23500000, 0

回复引用查看

#5楼 [楼主] 2006-04-17 22:50 浪人|努力

高斯投影变换_cpp
#include
#include
#include
#include "CoorTrans.h"
////////////////////////////////////////////
// Common functions
////////////////////////////////////////////
double Dms2Rad(double Dms)
{
double Degree, Miniute;
double Second;
int Sign;
double Rad;
if(Dms >= 0)
Sign = 1;
else
Sign = -1;
Dms = fabs(Dms);
Degree = floor(Dms);
Miniute = floor(fmod(Dms * 100.0, 100.0));
Second = fmod(Dms * 10000.0, 100.0);
Rad = Sign * (Degree + Miniute / 60.0 + Second / 3600.0) * PI / 180.0;
return Rad;
}
double Rad2Dms(double Rad)
{
double Degree, Miniute;
double Second;
int Sign;
double Dms;
if(Rad >= 0)
Sign = 1;
else
Sign = -1;
Rad = fabs(Rad * 180.0 / PI);
Degree = floor(Rad);
Miniute = floor(fmod(Rad * 60.0, 60.0));
Second = fmod(Rad * 3600.0, 60.0);
Dms = Sign * (Degree + Miniute / 100.0 + Second / 10000.0);
return Dms;
}
///////////////////////////////////////////////////
// Definition of PrjPoint
///////////////////////////////////////////////////
BOOL PrjPoint::BL2xy()
{
double X, N, t, t2, m, m2, ng2;
double sinB, cosB;
X = A1 * B * 180.0 / PI + A2 * sin(2 * B) + A3 * sin(4 * B) + A4 * sin(6 *
B);
sinB = sin(B);
cosB = cos(B);
t = tan(B);
t2 = t * t;
N = a / sqrt(1 - e2 * sinB * sinB);
m = cosB * (L - L0);
m2 = m * m;
ng2 = cosB * cosB * e2 / (1 - e2);
x = X + N * t * ((0.5 + ((5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24.0 + (61 -
58 * t2 + t2 * t2) * m2 / 720.0) * m2) * m2);
y = N * m * ( 1 + m2 * ( (1 - t2 + ng2) / 6.0 + m2 * ( 5 - 18 * t2 + t2 * t
2 + 14 * ng2 - 58 * ng2 * t2 ) / 120.0));
y += 500000;
return TRUE;
}
BOOL PrjPoint::xy2BL()
{
double sinB, cosB, t, t2, N ,ng2, V, yN;
double preB0, B0;
double eta;
y -= 500000;
B0 = x / A1;
do
{
preB0 = B0;
B0 = B0 * PI / 180.0;
B0 = (x - (A2 * sin(2 * B0) + A3 * sin(4 * B0) + A4 * sin(6 * B0))) / A1;
eta = fabs(B0 - preB0);
}while(eta > 0.000000001);
B0 = B0 * PI / 180.0;
B = Rad2Dms(B0);
sinB = sin(B0);
cosB = cos(B0);
t = tan(B0);
t2 = t * t;
N = a / sqrt(1 - e2 * sinB * sinB);
ng2 = cosB * cosB * e2 / (1 - e2);
V = sqrt(1 + ng2);
yN = y / N;
B = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN /
12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0)
* V * V * t / 2;
L = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24
* t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB
;
return TRUE;
}
BOOL PrjPoint::SetL0(double dL0)
{
L0 = Dms2Rad(dL0);
return TRUE;
}
BOOL PrjPoint::SetBL(double dB, double dL)
{
B = Dms2Rad(dB);
L = Dms2Rad(dL);
B = dB;
L = dL;
BL2xy();
return TRUE;
}
BOOL PrjPoint::GetBL(double *dB, double *dL)
{
*dB = Rad2Dms(B);
*dL = Rad2Dms(L);
return TRUE;
}
BOOL PrjPoint::Setxy(double dx, double dy)
{
x = dx;
y = dy;
xy2BL();
return TRUE;
}
BOOL PrjPoint::Getxy(double *dx, double *dy)
{
*dx = x;
*dy = y;
return TRUE;
}
///////////////////////////////////////////////////
// Definition of PrjPoint_IUGG1975
///////////////////////////////////////////////////
PrjPoint_IUGG1975::PrjPoint_IUGG1975()
{
a = 6378140;
f = 298.257;
e2 = 1 - ((f - 1) / f) * ((f - 1) / f);
e12 = (f / (f - 1)) * (f / (f - 1)) - 1;
A1 = 111133.0047;
A2 = -16038.5282;
A3 = 16.8326;
A4 = -0.0220;
}
PrjPoint_IUGG1975::~PrjPoint_IUGG1975()
{
}
///////////////////////////////////////////////////
// Definition of PrjPoint_Krasovsky
///////////////////////////////////////////////////
PrjPoint_Krasovsky::PrjPoint_Krasovsky()
{
a = 6378245;
f = 298.3;
e2 = 1 - ((f - 1) / f) * ((f - 1) / f);
e12 = (f / (f - 1)) * (f / (f - 1)) - 1;
A1 = 111134.8611;
A2 = -16036.4803;
A3 = 16.8281;
A4 = -0.0220;
}
PrjPoint_Krasovsky::~PrjPoint_Krasovsky()
{
}
回复引用查看

#6楼 [楼主] 2006-04-17 22:50 浪人|努力

经纬度BL换算到高斯平面直角坐标XY（高斯投影正算）的源码及算法

在 gis 的帖子"用excel完成gps坐标转换的简易方法 " 基础上,
我给出对应的vb程序段,我在 evb 开发的gps定位功能中，用它实现坐标换算(具体的参数请参考gis 的帖子)。
感觉速度比较快，效果比较好。所以帖上来，希望与名位交流：
=====================================
'经纬度bl换算到高斯平面直角坐标xy（高斯投影正算）

private function bl2xy(byref a2 as double, byref f2 as double, byref e2 as double, _
byref s2 as double, byref t2 as double) as boolean
'a2 输入中央子午线，以度.分形式输入，如115度30分则输入115.30; 起算数据l0
'f2 以度小数形式输入经度值, l
'e2 以度小数形式输入纬度值，b
's2 计算结果，横坐标y
't2 计算结果，纵坐标x
'投影带号计算　n=[l/6]+1 如：测得经度103.xxxx，故n=[103.x/6]+1=17+1=18
'中央经线经度　l0 = n*6-3 = [l/6]*6+3

dim b2 as double
'dim g2 as double
dim h2 as double
dim i2 as double
dim j2 as double
dim k2 as double
dim l2 as double
dim m2 as double
dim n2 as double
dim o2 as double
dim p2 as double
dim q2 as double
dim r2 as double

b2 = int(a2) + (int(a2 * 100) - int(a2) * 100) / 60 + (a2 * 10000 - int(a2 * 100) * 100) / 3600
'把l0化成度(a2)
'g2 = f2 - b2 ' l -l0
'h2 = g2 / 57.2957795130823 '化作弧度
h2 = (f2 - b2) / 57.2957795130823 '将经差的单位化为弧度
i2 = tan(e2 / 57.2957795130823) 'tan (b)
j2 = cos(e2 / 57.2957795130823) ' cos (b)
k2 = 0.006738525415 * j2 * j2
l2 = i2 * i2
m2 = 1 + k2
n2 = 6399698.9018 / sqr(m2)
o2 = h2 * h2 * j2 * j2
p2 = i2 * j2
q2 = p2 * p2
r2 = (32005.78006 + q2 * (133.92133 + q2 * 0.7031))
s2 = ((((l2 - 18) * l2 - (58 * l2 - 14) * k2 + 5) * o2 / 20 + m2 - l2) * o2 / 6 + 1) * n2 * (h2 * j2)
s2 = s2 + 18500000 '在计算的基础上加上了“带号”（18）和“东移”（500km）
'计算结果,横坐标y
t2 = 6367558.49686 * e2 / 57.29577951308 - p2 * j2 * r2 + ((((l2 - 58) * l2 + 61) * _
o2 / 30 + (4 * k2 + 5) * m2 - l2) * o2 / 12 + 1) * n2 * i2 * o2 / 2
'计算结果，纵坐标x
'msgbox "pts2=" & s2 & " pt t2=" & t2

bl2xy = true
end function
回复引用查看

#7楼 [楼主] 2006-04-17 22:52 浪人|努力

回复引用查看

#8楼 2006-11-09 22:16 jingziw [未注册用户]

@浪人|努力

a = 6378140;
f = 298.257;
e2 = 1 - ((f - 1) / f) * ((f - 1) / f);
e12 = (f / (f - 1)) * (f / (f - 1)) - 1;
A1 = 111133.0047;
A2 = -16038.5282;
A3 = 16.8326;
A4 = -0.0220;
能不能告诉我上面的参数都有什么意义
回复引用