|
阅读：10232回复：29
坐标转换和单位转换的函数源代码
大家看看，请指正！
'函数：经纬度坐标转换成大地坐标 '输入参数：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" /> [此贴子已经被作者于2003-8-15 9:59:28编辑过]
|
|
|
|
1楼#
发布于：2007-12-16 22:03
好像坐标转换mo没有mapgis方便，不知道mo有没有现成的函数，毕竟不是每个人数学都这么好。
|
|
|
2楼#
发布于：2007-11-06 23:18
支持开源、共享
|
|
|
|
3楼#
发布于：2007-11-06 22:20
<P>谢谢</P>
|
|
|
4楼#
发布于：2007-10-09 19:28
很好，非常好，太好了
|
|
|
5楼#
发布于：2007-07-09 17:46
<P>很不错的哦</P>
|
|
|
6楼#
发布于：2006-11-29 21:12
<P>本人本科毕业论文正是做这个坐标转换的,可是我一点也不懂,老师叫我上网拷,找了,好久,终于让我找到了,真是感谢楼主!</P>
<P>不知楼主有没有我国目前三种坐标系统两两之间转换的VB原代码,有的话,能发一下吗?</P> |
|
|
7楼#
发布于：2006-06-30 00:33
kap
|
|
|
|
8楼#
发布于：2006-06-29 15:14
<P>感谢总统先</P><img src="images/post/smile/dvbbs/em02.gif" />
|
|
|
|
9楼#
发布于：2006-06-26 15:14
<img src="images/post/smile/dvbbs/em02.gif" />
|
|
上一页
下一页