|
|
|
|
2019年04月19日 |
上次
經緯度轉xy座標及xy座標轉經緯度及
台電電桿座標轉換
Chday169
前因下載手機地質羅盤把玩，因此上網尋找經緯度轉換為Tw97XY座標之轉換公式，看到Victor's個人部落格及多位網友已將轉換公式寫成程式，並在網路發佈、分享。因為程式源碼多是以c++或python為主，我一向比較熟悉vb6，故先將其改寫成vb6測試。後又看到有人在討論台電電桿編碼轉換成Tw67 (Columbia Yen tech547@yahoo.com,DanJacobson等)，一時性起，也將其改寫成vb6測試。為此還特別去學習Python語言，人家是臨老入花欉，老夫是臨老抓蟒蛇(python)。
這是我第一個Python程式，故特別將其在部落格發佈及分享朋友。有關經緯度轉換為Tw97XY座標之討論因限於篇幅，不擬在此討論，有興趣的朋友可參考上河圖書出版公司之網站(有經緯度轉為Tw97XY及電桿編碼轉等雙向轉換
App)。在此僅發佈Python程式碼(經緯度轉xy座標及xy座標轉經緯度取材自Victor's部落格；台電電桿座標轉換主要參考Columbia Yen之 Basic，在此一併謝)
。
有興趣之朋友可連結下面之網頁，檢視程式原碼。
下面則為vb6之部份程式碼
Private Sub LatlonToXY(aIn As Double, bIn As Double, k0In As Double, dxIn As Double, dyIn As Double, lonIn As Double, latIn As Double, xOut As Double, yOut As Double)
Dim a As Double, b As Double, long0 As Double, k0 As Double, dx As Double, dy As Double
Dim e As Double, e2 As Double, n As Double, nu As Double, p As Double, D As Double, S As Double
Dim K1 As Double, K2 As Double, K3 As Double, K4 As Double, E4 As Double, K5 As Double, c As Double
If (lonIn >= 120 * Ra And lonIn <= 122 * Ra) And (latIn >= 21.8 * Ra And latIn <= 25.5 * Ra) Then
gllongUse = 121 * Ra
Else
gllongUse = 119 * Ra
End If
a = aIn
b = bIn
long0 = gllongUse
k0 = k0In
dx = dxIn
dy = dyIn
lat = latIn
lon = lonIn
PicPrint.Print "a,b,long0,k0,dx,dy"; a; b; long0; k0; dx; dy
PicPrint.Print "lon0=", long0
'MsgBox ("a= " & a & "b=" & b)
e = (1 - b ^ 2 / a ^ 2) ^ 0.5
e2 = e ^ 2 / (1 - e ^ 2)
n = (a - b) / (a + b)
nu = a / (1 - (e ^ 2) * (Sin(lat) ^ 2)) ^ 0.5
p = lon - long0
a = a * (1 - n + (5 / 4#) * (n ^ 2 - n ^ 3) + (81 / 64#) * (n ^ 4 - n ^ 5))
b = (3 * a * n / 2#) * (1 - n + (7 / 8#) * (n ^ 2 - n ^ 3) + (55 / 64#) * (n ^ 4 - n ^ 5))
c = (15 * a * (n ^ 2) / 16#) * (1 - n + (3 / 4#) * (n ^ 2 - n ^ 3))
D = (35 * a * (n ^ 3) / 48#) * (1 - n + (11 / 16#) * (n ^ 2 - n ^ 3))
e = (315 * a * (n ^ 4) / 51#) * (1 - n)
S = a * lat - b * Sin(2 * lat) + c * Sin(4 * lat) - D * Sin(6 * lat) + e * Sin(8 * lat)
K1 = S * k0
K2 = k0 * nu * Sin(2 * lat) / 4#
K3 = (k0 * nu * Sin(lat) * (Cos(lat) ^ 3) / 24#) * (5 - Tan(lat) ^ 2 + 9 * e2 * (Cos(lat) ^ 2) + 4 * (e2 ^ 2) * (Cos(lat) ^ 4))
yOut = K1 + K2 * (p ^ 2) + K3 * (p ^ 4)
K4 = k0 * nu * Cos(lat)
K5 = (k0 * nu * (Cos(lat) ^ 3) / 6#) * (1 - Tan(lat) ^ 2 + e2 * (Cos(lat) ^ 2))
xOut = K4 * p + K5 * (p ^ 3) + dxIn
'return x, y
End Sub
Private Sub xyToLatLon(aIn As Double, bIn As Double, k0In As Double, dxIn As Double, dyIn As Double, xIn As Double, yIn As Double, LonOut As Double, LatOut As Double)
'static public double[] convert(TMParameter tm, double x, double y) {
Dim dx As Double, dy As Double, lon0 As Double, k0 As Double, a As Double, b As Double, c As Double, long0 As Double
Dim x As Double, y As Double, e As Double, D As Double, long0Use As Double
Dim strtpt As String, arrayA() As String
If (xIn >= 146600# And xIn <= 353400) And (yIn >= 2411767# And yIn <= 2821559) Then
gllongUse = 121 * Ra
Else
gllongUse = 119 * Ra
End If
dx = dxIn
dy = dyIn
x = xIn
y = yIn
a = aIn
b = bIn
k0 = k0In
dx = dxIn
e = (1 - b ^ 2 / a ^ 2) ^ 0.5
lon0 = gllongUse
x = x - dx
y = y - dy
'// Calculate the Meridional Arc
Dim M As Double
'M = y/k0;
M = y / k0
'// Calculate Footprint Latitude
Dim mu As Double, E1 As Double, J1 As Double, J2 As Double, J3 As Double, J4 As Double, fp As Double
mu = M / (a * (1# - (e ^ 2) / 4# - 3 * (e ^ 4) / 64# - 5 * (e ^ 6) / 256#))
E1 = (1# - ((1# - (e ^ 2)) ^ 0.5)) / (1# + ((1# - (e ^ 2)) ^ 0.5))
J1 = (3 * E1 / 2 - 27 * (E1 ^ 3) / 32#)
J2 = (21 * (E1 ^ 2) / 16 - 55 * (E1 ^ 4) / 32#)
J3 = (151 * (E1 ^ 3) / 96#)
J4 = (1097 * (E1 ^ 4) / 512#)
fp = mu + J1 * Sin(2 * mu) + J2 * Sin(4 * mu) + J3 * Sin(6 * mu) + J4 * Sin(8 * mu)
'// Calculate Latitude and Longitude
Dim e2 As Double, C1 As Double, T1 As Double, R1 As Double, N1 As Double
e2 = ((e * a / b) ^ 2)
C1 = (e2 * Cos(fp) ^ 2)
T1 = (Tan(fp) ^ 2)
R1 = a * (1 - (e ^ 2)) / ((1 - (e ^ 2) * (Sin(fp) ^ 2)) ^ (3# / 2#))
N1 = a / ((1 - (e ^ 2) * (Sin(fp) ^ 2)) ^ 0.5)
e2 = ((e * a / b) ^ 2)
D = x / (N1 * k0)
'// lat
Dim Q1 As Double, Q2 As Double, Q3 As Double, Q4 As Double, Q5 As Double, Q6 As Double, Q7, lat As Double
Q1 = N1 * Tan(fp) / R1
Q2 = ((D ^ 2) / 2#)
Q3 = (5 + 3 * T1 + 10 * C1 - 4 * (C1 ^ 2) - 9 * e2) * (D ^ 4) / 24#
Q4 = (61 + 90 * T1 + 298 * C1 + 45 * (T1 ^ 2) - 3 * (C1 ^ 2) - 252 * e2) * (D ^ 6) / 720#
lat = fp - Q1 * (Q2 - Q3 + Q4)
'// long
Q5 = D
Q6 = (1 + 2 * T1 + C1) * (D ^ 3) / 6
Q7 = (5 - 2 * C1 + 28 * T1 - 3 * (C1 ^ 2) + 8 * e2 + 24 * (T1 ^ 2)) * (D ^ 5) / 120#
lon = lon0 + (Q5 - Q6 + Q7) / Cos(fp)
lat = lat * Da
lon = lon * Da
' PicPrint.Print "lon= " & Format(lon, "#######.######") & ", lat= " & Format(lat, "#######.######")
LatOut = lat
LonOut = lon
'PicPrint.Print "lon= " & Round(lon, 6) & " lat= " & Round(lat, 6)
End Sub
下面則簡單介紹台電電桿編碼轉換成Tw67之原理，台電電桿編碼係將台灣地區劃分為如下圖(擷錄自wiki.osgeo.org)A-Z共26區(實際只21區)。每區之左下座標原點為
"A , 170000, 2750000" "B , 250000, 2750000" "C , 330000, 2750000"
"D , 170000, 2700000" "E , 250000, 2700000" "F , 330000, 2700000"
"G , 170000, 2650000" "H , 250000, 2650000" "I , 330000, 2650000"
"J , 90000, 2600000" "K , 170000, 2600000" "L , 250000, 2600000"
"M , 90000, 2550000" "N , 170000, 2550000" "O , 250000, 2550000"
"P , 90000, 2500000" "Q , 170000, 2500000" "R , 250000, 2500000"
"S , 90000, 2450000" "T , 170000, 2450000" "U , 250000, 2450000"
"V , 170000, 2400000"
表一(取材自wiki.osgeo.org網站)
以電桿編碼G8150,HD78為例，就是以G區座標原點(x=170000, y=2650000)為原點，G8150為向右第81單位(800m)，向上第50單位(500m)之參考圖幅，即向右81x800=64800m，向上50x500=25000m；至於HD78則再細分為向右H單位(100m)，向上D單位(100m)之參考圖幅(8x5方格)，如座標原點編號為A，則A=0，B=1，
C=2，D=3，E=4，,F=5，G=6，H=7，則H為向右7X100=700m，則D為向上3X100=300m，78則代表向右7個單位(10m)70m， 則代表向上8個單位(10m)80m。則
x=170000+81x800+7x100+7x10=170000+64800+700+70=235,570m；
y=2650000+50x500+3x100+8x10=2650000+25000+300+80=2,675,380m
如HD78編號係以第一格為A，即A=1，B=2，C=3，D=4，E=5，F=6，G=7，H=8；則上述應改為
x=170000+81x800+(8-1)x100+7x10=170000+64800+700+70=235,570m；
y=2650000+50x500+(4-1)x100+8x10=2650000+25000+300+80=2,675,380m
至於，xy座標反算轉換為電桿編碼，則相當簡單。只要將xy座標比對表一，檢視其落在那一個區。只要檢視(x=235,570,y= 2,675,380)只要合乎下列條件；
(x >= Xorg(i) & x <= Xorg(i) + 80000) & (y >= Yorg(i) &y <= Yorg(i) + 50000)
(i=0~21)就可知道是落在第i區。在本例，i=6，即G區G=(170000, 2650000)。
(235570-170000)//800=81，(235570-170000-81x800)//100=7(H)，(2357570-170000-81x800-7x100)//10=7
(2675380-2650000)/500=50，(2675380-2650000-50x500)/100=3(D)，(2675380-2650000-50x500-3x100)//10=8
Private Sub xyOrg()
Dim i As Integer, Strin As String, arrayT() As String, Strxyorg(0 To 21), Xtpt As Double, Ytpt As Double
Dim StrA As String, StrB As String, StrC As String, StrD As String
Strxyorg(0) = "A , 170000, 2750000"
Strxyorg(1) = "B , 250000, 2750000"
Strxyorg(2) = "C , 330000, 2750000"
Strxyorg(3) = "D , 170000, 2700000"
Strxyorg(4) = "E , 250000, 2700000"
Strxyorg(5) = "F , 330000, 2700000"
Strxyorg(6) = "G , 170000, 2650000"
Strxyorg(7) = "H , 250000, 2650000"
Strxyorg(8) = "I , 330000, 2650000"
Strxyorg(9) = "J , 90000, 2600000"
Strxyorg(10) = "K , 170000, 2600000"
Strxyorg(11) = "L , 250000, 2600000"
Strxyorg(12) = "M , 90000, 2550000"
Strxyorg(13) = "N , 170000, 2550000"
Strxyorg(14) = "O , 250000, 2550000"
Strxyorg(15) = "P , 90000, 2500000"
Strxyorg(16) = "Q , 170000, 2500000"
Strxyorg(17) = "R , 250000, 2500000"
Strxyorg(18) = "S , 90000, 2450000"
Strxyorg(19) = "T , 170000, 2450000"
Strxyorg(20) = "U , 250000, 2450000"
Strxyorg(21) = "V , 170000, 2400000"
For i = 0 To 21
Strin = Strxyorg(i)
arrayT = Split(Trim(Strin), ",")
areaCode(i) = arrayT(0)
Xorg(i) = Val(arrayT(1))
Yorg(i) = Val(arrayT(2))
PicPrint.Print "x'y", i, areaCode(i), Xorg(i), Yorg(i)
Next i
End Sub
Private Sub dataSplit(Strin As String)
Dim lenstr As Integer, arrayA() As String, strUse As String, StrA As String, StrB As String
Strin = UCase(Strin)
arrayA = Split(Strin, ",")
StrA = arrayA(0)
areacodeU = Mid(StrA, 1, 1)
areaIndex = Val(Asc(areacodeU)) - 65
PP = Val(Mid(StrA, 2, 2))
QQ = Val(Mid(StrA, 4, 2))
PicPrint.Print "areacodeU, PP, QQ, areaIndex", areacodeU, PP, QQ, areaIndex
lenstr = Len(StrB)
LL = 0
MM = 0
StrB = arrayA(1)
rr = Mid(StrB, 1, 1)
ss = Mid(StrB, 2, 1)
TT = Val(Mid(StrB, 3, 1))
UU = Val(Mid(StrB, 4, 1))
If lenstr > 4 Then
LL = Val(Mid(StrB, 5, 1))
MM = Val(Mid(StrB, 6, 1))
End If
PicPrint.Print "rr, ss, TT, UU, LL, MM", rr, ss, TT, UU, LL, MM
End Sub
Private Sub Cmdrun_Click()
Dim Strin As String, x As Double, y As Double, i As Integer, PPtpt As Integer, QQtpt As Integer, xleftB As Integer, yleftB As Integer, chrZone As String
Dim ssTpt As String, rrTpt As String, TTtpt As Integer, UUtpt As Integer, LLtpt As Integer, MMtpt As Integer, xleft As Double, yleft As Double
Dim vtpt1 As Double, vtpt2 As Double, strTaipower As String
Call xyOrg
'Strin = "b5919,eb02"
Strin = "g8150,hd78"
PicPrint.Print Strin
Call dataSplit(Strin)
PicPrint.Print "Val(Asc(rr) , Asc("; A; ")) ", Asc(rr), Asc("A")
PicPrint.Print "areaIndex", areaIndex
PicPrint.Print "xorg= ", Xorg(areaIndex)
PicPrint.Print "yorg= ", Yorg(areaIndex)
x = Xorg(areaIndex) + PP * 800# + Val(Asc(rr) - Asc("A")) * 100# + TT * 10# + LL
y = Yorg(areaIndex) + QQ * 500# + Val(Asc(ss) - Asc("A")) * 100# + UU * 10# + MM
PicPrint.Print "X.Y", x, y
'check
For i = 0 To 21
If (x >= Xorg(i) And x <= Xorg(i) + 80000#) And (y >= Yorg(i) And y <= Yorg(i) + 50000#) Then
chrZone = Chr(Str(i + 65))
PicPrint.Print i, chrZone
xleft = x - Xorg(i)
yleft = y - Yorg(i)
PPtpt = Int(xleft / 800)
QQtpt = Int(yleft / 500)
PicPrint.Print PPtpt, QQtpt
vtpt1 = (xleft - PPtpt * 800#) / 100
vtpt2 = (yleft - QQtpt * 500#) / 100
PicPrint.Print vtpt1, vtpt2
rrTpt = Chr(Str(Int(vtpt1) + 65))
ssTpt = Chr(Str(Int(vtpt2) + 65))
vtpt1 = ((xleft - PPtpt * 800# - Int(vtpt1) * 100#) / 10#)
vtpt2 = ((yleft - QQtpt * 500# - Int(vtpt2) * 100#) / 10#)
TTtpt = Int(vtpt1)
UUtpt = Int(vtpt2)
PicPrint.Print "rr,ss", rrTpt, ssTpt, TTtpt, UUtpt
LL = xleft - PPtpt * 800# - Int(vtpt1) * 100# - TTtpt * 10#
MM = yleft - QQtpt * 500# - Int(vtpt2) * 100# - UUtpt * 10#
strTaipower = Trim(chrZone & PPtpt & QQtpt) & "," & Trim(rrTpt & ssTpt & TTtpt & UUtpt)
strTaipower = Trim(strTaipower)
If LL > 0.1 Or MM > 0.1 Then strTaipower = Trim(strTaipower & LL & MM)
PicPrint.Print strTaipower
Exit For
End If
Next i
End Sub
修改此網站的日期： 2019年03月23日