坐标转换和单位转换的函数源代码

大家看看，请指正！

'函数：经纬度坐标转换成大地坐标
'输入参数：ZB,纬度,经度,坐标系统（如：80坐标）,中央经线
'输出参数：大地X,大地Y

Sub BLtoXY(Zb As ZbxPa, B As Double, L As Double, tyzbx As String, CenterLine As Single)
    Dim ZS As Double, Zc As Double
    Dim Ztan As Double, lc As Double, Nf As Double, n As Double
    Dim Q As Double, x As Double, y As Double
    Dim Jch
    Dim RstX As Double, RstY As Double
    Call SleParamer(tyzbx, Zb)
    ZS = Sin(B * W)
    Zc = Cos(B * W)
    Ztan = Tan(B * W)
    lc = (L - CenterLine) * 3600
    Nf = Zb.e1 * Zc * Zc
    n = Zb.ea / Sqr(1 - Zb.E * ZS * ZS)
    Q = Zb.tC0 * B * W + Zb.tC1 * Sin(2 * B * W) + Zb.tC2 * Sin(4 * B * W) + Zb.tC3 * Sin(6 * B * W) + Zb.tC4 * Sin(8 * B * W) + Zb.tC5 * Sin(10 * B * W)
    x = Q + n * Ztan * Zc ^ 2 * lc ^ 2 / (2 * PI * PI)
    x = x + n * Ztan * (5 - Ztan ^ 2 + 9 * Nf ^ 2 + 4 * Nf ^ 4) * Zc ^ 4 * lc ^ 4 / (24 * PI ^ 4)
    RstX = x + n * Ztan * (61 - 58 * Ztan ^ 2 + Ztan ^ 4 + 270 * Nf ^ 2 - 330 * Nf ^ 2 * Ztan ^ 2) * Zc ^ 6 * lc ^ 6 / (720 * PI ^ 6)
    RstX = RstX + (RstX < 1000000) * 0.000005
    y = n * Zc * (lc / PI)
    y = y + n * (1 - Ztan ^ 2 + Nf) * Zc ^ 3 * lc ^ 3 / (6 * PI ^ 3)
    RstY = y + n * (5 - 18 * Ztan ^ 2 + Ztan ^ 4 + 14 * Nf - 58 * Nf * Ztan ^ 2) * Zc ^ 5 * lc ^ 5 / (120 * PI ^ 5)
    RstY = RstY + Sgn(RstY) * 0.000005
    RstY = RstY '+ 500000  '大庆测试通过
    gRstX = RstY + CF + DaiHao: gRstY = RstX
    'Debug.Print RstX, RstY
End Sub

'2.==============================================================================================
'函数：大地坐标转换成经纬度坐标
'输入参数：坐标系(54|80)、自然值、度带、中央经线、带号Boolean、Zb、X、Y
'输出参数：经度、纬度

Sub XYtoBL(zbx As Integer, ByVal btt As String, ld As String, ByVal mCenterline As Single, ByVal DaiNO As Boolean, Zb As ZbxPa, x As Double, y As Double)
    '判断带号
    Dim Dh As Integer
    If DaiNO Then
        '有带号
        'YZbxtxt = "1980 北京坐标系"
        'WZbxtxt = "1954 北京坐标系"
        y = Val(Mid$(CStr(y), 3))
    Else
        y = y
    End If
    If btt = "50自然值" Then
        y = y - CF
    ElseIf btt = "D50自然值" Then
        y = Val(Mid$(Str(y), 4)) - CF
        Dh = Val(Mid(Str(y), 2, 2))
        If y < 0 And Abs(y) > 1000000 Then
            y = Val(Mid$(Trim(Str(y)), 3)) - CF
            '===================================
            Dh = Val(Mid(Trim(Str(y)), 1, 2))
            '===================================
        End If
        If ld = "6度带" Then
            mCenterline = OriDMStoDeg(Dh * 6 - 3)
        Else
            mCenterline = OriDMStoDeg(Dh * 3)
        End If
    Else
        y = y
    End If
    '----------------------------------------------------------------------------------
    Dim Zcos As Double, ZSin As Double
    Dim E As Double, z As Double, c As Double, t As Double, V As Double, U As Double
    Dim H As Double, n As Double, M As Double, c1 As Double, c2 As Double, c3 As Double
    Dim BF As Double, Q As Double
    Dim I As Integer
    Dim Tf, Nf, Mf, a1#, b1#, B, CF1#
    Dim ka As Double, kb As Double, kc As Double, ka1 As Double, kb1 As Double, kc1 As Double
    Dim kt1 As Double, kt2 As Double, kt3 As Double
    Call SleParamerGS(zbx, Zb)
    BF = x / (Zb.tA * W)
    Zcos = Cos(BF * W)
    ZSin = Sin(BF * W)
    ka = Zb.tB / (2 * Zb.tA)
    ka1 = ka
    kb = -Zb.tC / (4 * Zb.tA)
    kb1 = kb
    kc = Zb.tD / (6 * Zb.tA)
    kc1 = kc
    For I = 1 To 3
      kt1 = ka + ka * kb1 - 3 * ka * ka1 ^ 2 / 2 - 2 * kb * ka1
      kt2 = kb + ka * ka1
      kt3 = kc + ka * kb1 + ka * ka1 ^ 2 / 2 + 2 * kb * ka1
      ka1 = kt1
      kb1 = kt2
      kc1 = kt3
    Next
    Zb.tK1 = 2 * ka1 + 4 * kb1 + 6 * kc1
    Zb.tK2 = 8 * kb1 + 32 * kc1
    Zb.tK3 = 32 * kc1
    BF = BF * W + Zcos * (Zb.tK1 * ZSin - Zb.tK2 * ZSin ^ 3 + Zb.tK3 * ZSin ^ 5 + Zb.tK4 * ZSin ^ 7)
    Tf = Tan(BF)  ' 计算 B
    Nf = Zb.ea / Sqr(1 - Zb.E * Sin(BF) * Sin(BF))
    Mf = Zb.e1 * Cos(BF) ^ 2
    a1# = Tf * (1 + Mf) * (y / Nf) ^ 2 / 2
    b1# = Tf * (1 + Mf) * (5 + 3 * Tf ^ 2 + Mf - 9 * Mf * Tf ^ 2) * y ^ 4 / (24 * Nf ^ 4) 'Tf * (5 + 3 * Tf * Tf + 6 * Mf * Mf * (1 - Tf * Tf)) * y ^ 4 / (24 * Nf ^ 4)
    c1# = Tf * (1 + Mf) * (61 + 90 * Tf ^ 2 + 45 * Tf ^ 4) * y ^ 6 / (720 * Nf ^ 6) 'Tf * (61 + 45 * Tf * Tf * (2 + Tf * Tf)) * y ^ 6 / (720 * Nf ^ 6)
    B = BF - a1# + b1# - c1#
    gLat = B / W
    CF1# = Cos(BF) ' 计算 L
    a1# = y / (Nf * CF1#)
    b1# = (1 + 2 * Tf * Tf + Mf) * (y / Nf) ^ 3 / (6 * CF1#)
    c1# = (5 + 28 * Tf ^ 2 + 24 * Tf ^ 4 + Mf * (6 + 8 * Tf ^ 2)) * (y / Nf) ^ 5 / (120 * CF1#)
    gLgt = mCenterline + (a1# - b1# + c1#) / W
End Sub


'经纬度==〉高斯坐标(坐标转换的中间过程函数)
Public Sub SleParamer(tyzbx As String, Zb As ZbxPa)
    If tyzbx = "54" Then
        '54坐标
        Zb.ea = a01 '6378245 ' Csb("a")   ' 6378140#
        Zb.eb = b01 '6356863.01877 'Csb("b")      ' 6356755.28815753
        Zb.E = str_E01 '0.00669342162297 'Csb("e2")    '0.006694384999588    '   e
        Zb.e1 = e01 ' 0.00673852541468 'Csb("e3") ' 0.006739501819473   '   e'
    Else
        '80坐标
        Zb.ea = a02 ' 6378245 ' Csb("a")   ' 6378140#
        Zb.eb = b02 ' 6356863.01877 'Csb("b")      ' 6356755.28815753
        Zb.E = str_E02 ' 0.00669342162297 'Csb("e2")    '0.006694384999588    '   e
        Zb.e1 = e02 ' 0.00673852541468 'Csb("e3") ' 0.006739501819473   '   e'
    End If
    Zb.tA = Zb.ea * (1 - Zb.E) * (1 + 3 * Zb.E / 4 + 45 * Zb.E ^ 2 / 64 + 175 * Zb.E ^ 3 / 256 + 11025 * Zb.E ^ 4 / 16384)
    Zb.tB = Zb.ea * (1 - Zb.E) * (3 * Zb.E / 4 + 15 * Zb.E ^ 2 / 16 + 525 * Zb.E ^ 3 / 512 + 2205 * Zb.E ^ 4 / 2048)
    Zb.tC = Zb.ea * (1 - Zb.E) * (15 * Zb.E ^ 2 / 64 + 105 * Zb.E ^ 3 / 256 + 2205 * Zb.E ^ 4 / 4096 + 10395 * Zb.E ^ 5 / 16384)
    Zb.tD = Zb.ea * (1 - Zb.E) * (35 * Zb.E ^ 3 / 512 + 315 * Zb.E ^ 4 / 2048 + 31185 * Zb.E ^ 5 / 131072)
    Zb.tE = Zb.ea * (1 - Zb.E) * (315 * Zb.E ^ 4 / 16384 + 3465 * Zb.E ^ 5 / 65536)
    Zb.Tf = Zb.ea * (1 - Zb.E) * (693 * Zb.E ^ 5 / 131072)
    Zb.tC0 = Zb.tA
    Zb.tC1 = -Zb.tB / 2
    Zb.tC2 = Zb.tC / 4
    Zb.tC3 = -Zb.tD / 6
    Zb.tC4 = Zb.tE / 8
    Zb.tC5 = -Zb.Tf / 10
    Zb.cA = 1 + Zb.E / 2 + 3 * Zb.E ^ 2 / 4 + 5 * Zb.E ^ 3
    Zb.cB = Zb.E / 6 + 3 * Zb.E ^ 2 / 16 + 3 * Zb.E ^ 3 / 16
    Zb.cc = 3 * Zb.E ^ 2 / 80 + Zb.E ^ 3 / 16
    Zb.cD = (Zb.E ^ 3 / 112)
End Sub

'大地--〉经纬度(坐标转换的中间过程函数)
Public Sub SleParamerGS(cor_ID As Integer, Zb As ZbxPa)
    If cor_ID = 54 Then
        '54坐标
        Zb.ea = a01 '6378245 ' Csb("a")   ' 6378140#
        Zb.eb = b01 '6356863.01877 'Csb("b")      ' 6356755.28815753
        Zb.E = str_E01 '0.00669342162297 'Csb("e2")    '0.006694384999588    '   e
        Zb.e1 = e01 ' 0.00673852541468 'Csb("e3") ' 0.006739501819473   '   e'
    Else
        '80坐标
        Zb.ea = a02 ' 6378245 ' Csb("a")   ' 6378140#
        Zb.eb = b02 ' 6356863.01877 'Csb("b")      ' 6356755.28815753
        Zb.E = str_E02 ' 0.00669342162297 'Csb("e2")    '0.006694384999588    '   e
        Zb.e1 = e02 ' 0.00673852541468 'Csb("e3") ' 0.006739501819473   '   e'
    End If
    Zb.tA = Zb.ea * (1 - Zb.E) * (1 + 3 * Zb.E / 4 + 45 * Zb.E ^ 2 / 64 + 175 * Zb.E ^ 3 / 256 + 11025 * Zb.E ^ 4 / 16384) ' + 43659 * Zb.E ^ 5 / 65536)
    Zb.tB = Zb.ea * (1 - Zb.E) * (3 * Zb.E / 4 + 15 * Zb.E ^ 2 / 16 + 525 * Zb.E ^ 3 / 512 + 2205 * Zb.E ^ 4 / 2048) '+ 72765 * Zb.E ^ 5 / 65536)
    Zb.tC = Zb.ea * (1 - Zb.E) * (15 * Zb.E ^ 2 / 64 + 105 * Zb.E ^ 3 / 256 + 2205 * Zb.E ^ 4 / 4096 + 10395 * Zb.E ^ 5 / 16384)
    Zb.tD = Zb.ea * (1 - Zb.E) * (35 * Zb.E ^ 3 / 512 + 315 * Zb.E ^ 4 / 2048 + 31185 * Zb.E ^ 5 / 131072)
    Zb.tE = Zb.ea * (1 - Zb.E) * (315 * Zb.E ^ 4 / 16384 + 3465 * Zb.E ^ 5 / 65536)
    Zb.Tf = Zb.ea * (1 - Zb.E) * (693 * Zb.E ^ 5 / 131072)
    Zb.tC0 = Zb.tA
    Zb.tC1 = -Zb.tB / 2
    Zb.tC2 = Zb.tC / 4
    Zb.tC3 = -Zb.tD / 6
    Zb.tC4 = Zb.tE / 8
    Zb.tC5 = -Zb.Tf / 10
    Zb.cA = 1 + Zb.E / 2 + 3 * Zb.E ^ 2 / 4 + 5 * Zb.E ^ 3
    Zb.cB = Zb.E / 6 + 3 * Zb.E ^ 2 / 16 + 3 * Zb.E ^ 3 / 16
    Zb.cc = 3 * Zb.E ^ 2 / 80 + Zb.E ^ 3 / 16
    Zb.cD = (Zb.E ^ 3 / 112)
End Sub

'==========================================
'函数功能：度分秒换算为度
'输入参数：度分秒(45.3030表示45度30分30秒)
'输出参数：度(45.5083333333333)

'==========================================
Public Function OriDMStoDeg(DMS As Double) As Double
    Dim Strdms$, d#, M#, Feng#
    Dim n As Integer
    Strdms$ = Str(DMS)
    n = InStr(Strdms$, ".")
    If n = 0 Then
        OriDMStoDeg = DMS
        Exit Function
    Else
        Strdms$ = Str(DMS) + "0000"
    End If
    d# = Val(Mid$(Strdms$, 1, n - 1))
    Feng# = (Val(Mid$(Strdms$, n + 1, 2))) / 60
    M# = (Val(Mid$(Strdms$, n + 3, 2) + "." + Mid$(Strdms$, n + 5))) / 3600
    OriDMStoDeg = (d# + Feng# + M#)
End Function

'==========================================
'函数功能：度换算为度分秒
'输入参数：度(45.5083333333333)
'输出参数：度分秒(45.3030表示45度30分30秒)

'==========================================
Public Function OriDegToDms(Deg As Double) As Double
    Dim TempDeg#, Zd#, AllF, Zf, strZf, p
     TempDeg# = Deg
     If InStr(Str(Deg), ".") = 0 Then
        OriDegToDms = Str(Deg) + ".0000"
        Exit Function
     End If
     Zd# = Int(TempDeg#)
     AllF = (TempDeg# - Zd#) * 60
     Zf = Int(AllF)
     If Zf < 10 Then
        strZf = "0" + Trim$(Str(Zf))
     Else
        strZf = Trim$(Str(Zf))
        If Val(strZf) = 60 Then Zd# = Zd# + 1
     End If
    Dim Allm, Strm
     Allm = (AllF - Zf) * 60
     If Allm < 10 Then
        Strm = "0" + Trim$(Str(Int(Allm)))
     Else
        If Allm >= 60 Then
            strZf = Str(Val(strZf) + 1)
            Strm = Str(Allm - 60)
            If Val(strZf) = 60 Then
                Zd# = Zd# + 1
                strZf = "00"
            End If
        Else
            Strm = Trim$(Str(Allm))
        End If
        p = InStr(Strm, ".")
        If p <> 0 Then
            Strm = Left(Strm, p - 1) + Mid(Strm, p + 1)
        End If
     End If
     OriDegToDms = CDbl(Str(Zd#) + "." + strZf + Strm)
End Function


<img src="images/post/smile/dvbbs/em09.gif" />

好像坐标转换mo没有mapgis方便，不知道mo有没有现成的函数，毕竟不是每个人数学都这么好。

支持开源、共享

<P>谢谢</P>

很好，非常好，太好了

<P>很不错的哦</P>

<P>本人本科毕业论文正是做这个坐标转换的,可是我一点也不懂,老师叫我上网拷,找了,好久,终于让我找到了,真是感谢楼主!</P>
<P>不知楼主有没有我国目前三种坐标系统两两之间转换的VB原代码,有的话,能发一下吗?</P>

kap

<P>感谢总统先</P><img src="images/post/smile/dvbbs/em02.gif" />

<img src="images/post/smile/dvbbs/em02.gif" />