高斯克吕格与经纬度坐标值转换代码
'高斯克吕格与经纬度坐标值转换代码
'Writen by Rodger Yuan 9 5 2006
'参考文献
'v0 0.1
'用于在经纬度坐标和高斯克吕格坐标之间的转换。
'高斯克吕格为一种投影，根据椭球体和基准面不同又有所区分，常用的北京54和西安80即
'采用这种投影方式，投影后的坐标为平面坐标系，单位为米
'现在参数的坐标系采用测绘坐标系，x为纵坐标，y为横坐标
'返回参数为自定义类型，双精度点
'调用转换函数前需要调用初始化过程进行初始化
'-------------------------------------------------------------------------------
'Public Sub init(ByVal TuoqiuCanshu As Canshu, ByVal Daihao As Integer)
'说明: 用于初始化转换参数
'TuoqiuCanshu 枚举类型，提供北京54、西安80和WGS84三个椭球参数
'Daihao 整型为高斯克吕格投影六度分带带号，取值为1～60
'-------------------------------------------------------------------------------
'Public Sub initEx(ByVal dE As Double, ByVal dN As Double, ByVal k0 As Double)
'说明: 用于进一步初始化转换参数（暂不提供）
'dE 东偏移
'dE 北偏移
'k0 比例因子
'-------------------------------------------------------------------------------
'Public Function JWgetGK(ByVal W As Double, ByVal J As Double) As PointD
'************************************************************************
'基本变量定义
Dim a As Double '椭球体长半轴
Dim b As Double '椭球体短半周
Dim f As Double '扁率
Dim e As Double '第一偏心率
Dim eq As Double '第二偏心率
Dim dh As Integer '带号
Dim FE As Double '东偏移
Dim FN As Double '北偏移
Dim L0 As Double '中央经度
Dim k0 As Double '比例因子
Const PI As Double = 3.14159265358979
Public Enum Canshu
Beijing54 = 0
Xian80 = 1
WGS84 = 2
End Enum
Public Type PointD
X As Double
Y As Double
End Type
Public Sub init(ByVal TuoqiuCanshu As Canshu, ByVal Daihao As Integer)
Select Case TuoqiuCanshu
'Krassovsky （北京54采用）6378245 6356863.0188
'IAG 75（西安80采用）6378140 6356755.2882
'WGS 84 6378137 6356752.3142
Case 0: '北京五四
a = 6378245
b = 6356863.0188
Case 1: '西安八零
a = 6378140
b = 6356755.2882
Case 2: 'WGS84
a = 6378137
b = 6356752.3142
End Select
f = (a - b) / a
e = Sqr(1 - (b / a) ^ 2)
eq = Sqr((a / b) ^ 2 - 1)
If Daihao < 1 Or Daihao > 60 Then Exit Sub
dh = Daihao
L0 = (6 * dh - 3) * PI / 180
k0 = 1
FE = 500000 + dh * 1000000
FN = 0
End Sub
Public Sub initEx(ByVal dE As Double, ByVal dN As Double, ByVal dk0 As Double)
End Sub
Public Function JWgetGK(ByVal W As Double, ByVal J As Double) As PointD
'给出经纬度坐标，转换为高克投影坐标
Dim BY As Double
Dim LX As Double
Dim TC As Double
Dim CC As Double
Dim AC As Double
Dim MC As Double
Dim NC As Double
Dim rx As Double
Dim ry As Double
Dim resultP As PointD
BY = W * PI / 180
LX = J * PI / 180
TC = Math.Tan(BY) ^ 2
CC = eq ^ 2 * Cos(BY) ^ 2
AC = (LX - L0) * Cos(BY)
MC = a * ((1 - e ^ 2 / 4 - 3 * e ^ 4 / 64 - 5 * e ^ 6 / 256) * BY - (3 * e ^ 2 / 8 + 3 * e ^ 4 / 32 + 45 * e ^ 6 / 1024) * Sin(2 * BY) + (15 * e ^ 4 / 256 + 45 * e ^ 6 / 1024) * Sin(4 * BY) - (35 * e ^ 6 / 3072) * Sin(6 * BY))
NC = a / Sqr(1 - e ^ 2 * (Sin(BY)) ^ 2)
rx = k0 * (MC + NC * Tan(BY) * (AC ^ 2 / 2 + (5 - TC + 9 * CC + 4 * CC ^ 2) * AC ^ 4 / 24) + (61 - 58 * TC + T ^ 2 + 270 * CC - 330 * TC * CC) * AC ^ 6 / 720)
ry = FE + k0 * NC * (AC + (1 - TC + CC) * AC ^ 3 / 6 + (5 - 18 * TC + TC ^ 2 + 14 * CC - 58 * TC * CC) * AC ^ 5 / 120)
resultP.X = rx
resultP.Y = ry
JWgetGK = resultP
End Function
Public Function GKgetJW(ByVal X As Double, ByVal Y As Double) As PointD
'给出高克投影坐标，转换为经纬度坐标
Dim BY As Double
Dim LX As Double
Dim e1 As Double
Dim FI As Double
Dim Mf As Double
Dim Bf As Double
Dim Tf As Double
Dim Cf As Double
Dim Nf As Double
Dim Rf As Double
Dim D As Double
Dim RW As Double '纬度
Dim RJ As Double '经度
Dim resultP As PointD
YE = Y
XN = X
e1 = (1 - b / a) / (1 + b / a)
Mf = (XN - FN) / k0
FI = Mf / (a * (1 - e ^ 2 / 4 - 3 * e ^ 4 / 64 - 5 * e ^ 6 / 256))
Bf = FI + (3 * e1 / 2 - 27 * e1 ^ 3 / 32) * Sin(2 * FI) + (21 * e1 ^ 2 / 16 - 55 * e1 ^ 4 / 32) * Sin(4 * FI) + (151 * e1 ^ 3 / 96) * Sin(6 * FI)
Tf = Tan(Bf) ^ 2
Cf = eq ^ 2 * Cos(Bf) ^ 2
Nf = a / Sqr(1 - e ^ 2 * Sin(Bf) ^ 2)
Rf = a * (1 - e ^ 2) / Sqr((1 - e ^ 2 * Sin(Bf) ^ 2) ^ 3)
D = (Y
E - FE) / (k0 * Nf)
RW = Bf - (Nf * Tan(Bf) / Rf) * (D ^ 2 / 2 - (5 + 3 * Tf + Cf - 9 * Tf * Cf) * D ^ 4 / 24 + (61 + 90 * Tf + 45 * Tf ^ 2) * D ^ 6 / 720)
RJ = L0 + 1 / Cos(Bf) * (D - (1 + 2 * Tf + Cf) * D ^ 3 / 6 + (5 + 28 * Tf + 6 * Cf + 8 * Tf * Cf + 24 * Tf ^ 2) * D ^ 5 / 120)
resultP.X = RW * 180 / PI
resultP.Y = RJ * 180 / PI
GKgetJW = resultP
End Function