地図座標変換

地図表示の基本となる部分です。地図データは緯度経度となっていますが、画面上に表示される場合は、地図の図法により変換作業を行うことになります。以下のプログラムでは１点のポイントの緯度経度を画面上の(X,Y)への変換処理を行います。

'LHT座標をピクチャーXY座標に変換する

Public Function LHTtopicXY_LpsubM!(px&, py&, linego%)
Dim bb As tdtype, hh&, ss#, sbai#
On Error GoTo xtper2
'LHTto球座標
If pzuho_kyu Then '球座標か？
sbai = z_LHTZearth(bb, px * LHTRADB, py * LHTRADB)
Else
Select Case lp_zuho '各図法ごとの計算を行う
Case MZ_CHOKO: z_LHTZchoko bb, px, py
Case MZ_SANSON: z_LHTZsanson bb, px, py
Case MZ_EKELT: z_LHTZekelt bb, px, py
Case MZ_MOLWA: z_LHTZmolwa bb, px, py
Case MZ_MELKA: z_LHTZmelka bb, px, py
Case MZ_MIRROR: z_LHTZmirror bb, px, py
Case MZ_ENSUI: z_LHTZensui bb, px, py
Case MZ_BONNU: z_LHTZbonnu bb, px, py
Case MZ_OKIMO: z_LHTZokimo bb, px, py
End Select
End If
'方位変換を行う
If housz <> 0 Then
If pzuho_kyu = False Then
'方位設定
ss = bb.xt * houcz - bb.yt * housz
bb.yt = bb.xt * housz + bb.yt * houcz
bb.xt = ss
End If
End If
'俯角変換を行う
If Not fukmode Then '980605
If pzuho_kyu And (lpm_mode = LPM_POINT Or sbai = -1) Then '球図法で１点か正距方位で縁ならば sbai は設定しない
Else
sbai = 1
End If
Else
sbai = z_LHTZ_Fukaku(bb, linego)
End If
'球座標をPIC座標に変換する
px = bb.xt * LHT90picbaiptwxl
py = bb.yt * LHT90picbaiptwyl + lp_psyp
'フライトビュー変換を行う
If lp_flaja Then
If lpm_mode = LPM_SHS And py < lp_chipy - lp_psm.y Then '970530
py = -1000
End If
hh = px * bancz - py * bansz
py = px * bansz + py * bancz
px = hh
End If
'左右逆転
'px = -px
LHTtopicXY_LpsubM = sbai
Exit Function
xtper2:
Select Case Err
Case 5 '980126
Resume Next
Case 6
Resume Next
Case 11
Debug.Print "LPSUB2 ０で除算"
Resume Next
Case Else
'Errmes Err, "LPSUB2"
Resume Next
End Select
End Function
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'球座標の計算　座標設定
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Private Function z_LHTZearth#(b1 As tdtype, ByVal sz#, ByVal cz#)
Dim b2x#, b2y#, b2z#, mtemp#
On Error GoTo xtper
With b1
'球座標算出
'地図中心をx=0,y=0,z=1に持ってくる回転
mtemp = Cos(cz)
b2x = Sin(sz) * mtemp
b2y = -Sin(cz)
b2z = Cos(sz) * mtemp
.xt = b2x * f(1).xt + b2y * f(2).xt + b2z * f(3).xt
.yt = b2x * f(1).yt + b2y * f(2).yt + b2z * f(3).yt
.zt = b2x * f(1).zt + b2y * f(2).zt + b2z * f(3).zt
'地名敷居値チェック
If lpm_mode = LPM_SHS And .zt < zski Then
.xt = -10: .yt = -10: .zt = -10 '980721
z_LHTZearth = -1
ElseIf Not fukmode Or efukamode = 0 Then '980605
If .zt >= zski Then '敷居値以上
z_LHTZearth = 1
Else
z_zzchk .xt, .yt, .zt
z_LHTZearth = 0
End If
End If
Select Case lp_zuho
Case MZ_KYOHO, MZ_RANBE
If Abs(.zt) >= 1 Then
.zt = Sgn(.zt) * 0.99999999
Else
'中心からの半径を求める
cz = Atn(.zt / Sqr(-.zt * .zt + 1))
cz = -cz + Atn(1) * 2
If lp_zuho = MZ_RANBE Then
cz = 2 * Sin(cz / 2)
End If
'xy角を求める
If .xt = 0 Then
.xt = 0
.yt = cz * Sgn(.yt)
Else
sz = Atn(.yt / .xt)
mtemp = Sgn(.xt)
.xt = cz * Cos(sz) * mtemp
.yt = cz * Sin(sz) * mtemp
End If
End If
End Select
End With
Exit Function
xtper:
Select Case Err
Case 6
Resume Next
Case 11
Debug.Print "LPSUB1 ０で除算"
Resume Next
Case Else
'Errmes Err, "LPSUB1"
Resume Next
End Select
End Function
'正角方位図法の計算
Private Sub z_LHTZchoko(b1 As tdtype, ByVal xx&, ByVal yy&)
b1.xt = (xx - lp_sbnow.X) * PLHT90
b1.yt = -(yy - lp_sbnow.y) * PLHT90
End Sub
'サンソン図法の計算
Private Sub z_LHTZsanson(b1 As tdtype, ByVal xx&, ByVal yy&)
b1.xt = (xx - lp_sbnow.X) * Cos(yy * LHTRADB) * PLHT90
b1.yt = -(yy - lp_sbnow.y) * PLHT90
End Sub
'サンソン図法（沖縄）
Private Sub z_LHTZokimo(b1 As tdtype, ByVal xx&, ByVal yy&)
Dim opy&
opy = yy
If (xx / DOUBAI < 132.5) And (yy / DOUBAI < 30) Then '正位置→沖縄移動位置
xx = xx + OKIX
yy = yy + OKIY
End If
b1.xt = (xx - lp_sbnow.X) * Cos(opy * LHTRADB) * PLHT90
b1.yt = -(yy - lp_sbnow.y) * PLHT90
End Sub
'エケルト図法の計算
Private Sub z_LHTZekelt(b1 As tdtype, ByVal xx&, ByVal yy&)
b1.xt = (xx - lp_sbnow.X) * 0.31184 * (1 + Cos(yy * LHTRADB)) * PLHT90
b1.yt = -(yy - lp_sbnow.y) * 0.62369 * PLHT90
End Sub
'モルワイデ図法の計算
Private Sub z_LHTZmolwa(b1 As tdtype, ByVal xx&, ByVal yy&)
Dim f0#, yhos#, h#
f0 = Atn(1) * 4 * Sin(lp_sbnow.y * LHTRADB)
h = z_LHTZmolwa_sub(f0, lp_bai)
yhos = -Sqr(2) * Sin(h)
f0 = Atn(1) * 4 * Sin(yy * LHTRADB)
h = z_LHTZmolwa_sub(f0, lp_bai)
b1.xt = (xx - lp_sbnow.X) * Atn(1) * 2 * Cos(h) * PLHT90
b1.yt = -Sqr(2) * Sin(h) - yhos
End Sub
'モルワイデ図法の計算サブ
Private Function z_LHTZmolwa_sub#(f0#, bai!)
Dim i&, bs#, bl#, h#
If Sgn(f0) >= 0 Then
bs = f0 / 2 - 0.5
bl = f0 / 2 + 0.5
Else
bs = f0 / 2 - 0.5
bl = f0 / 2 + 0.5
End If
For i = 0 To Log(1 / bai) + 17
h = (bs + bl) / 2
If 2 * h + Sin(2 * h) - f0 < 0 Then
bs = h
Else
bl = h
End If
Next
z_LHTZmolwa_sub = h
End Function
'メルカトル図法の計算
Private Sub z_LHTZmelka(b1 As tdtype, ByVal xx&, ByVal yy&)
Dim sz#, cz#
b1.xt = (xx - lp_sbnow.X) * PLHT90
If Abs(yy) >= 850000 Then
yy = Sgn(yy) * 850000
End If
sz = Log(Tan((LHT45 + Abs(yy) / 2) * LHTRADB)) / Log(10#) 'r * Log( tan(45+y/2))
cz = Log(Tan((LHT45 + Abs(lp_sbnow.y) / 2) * LHTRADB)) / Log(10#) 'r * Log( tan(45+y/2))
b1.yt = 1.5 * (Sgn(lp_sbnow.y) * cz - Sgn(yy) * sz)
End Sub
'ミラー図法の計算
Private Sub z_LHTZmirror(b1 As tdtype, ByVal xx&, ByVal yy&)
Dim sz#, cz#
If xx - lp_sbnow.X < -LHT360 Then
b1.xt = (xx - lp_sbnow.X + LHT360) * PLHT90
Else
b1.xt = (xx - lp_sbnow.X) * PLHT90
End If
If Abs(yy) >= LHT90 Then
yy = Sgn(yy) * LHT90 - 1
End If
sz = Log(Tan((LHT45 + Abs(yy) * 2 / 5) * LHTRADB)) / Log(10#) 'r * Log( tan(45+y/2))
cz = Log(Tan((LHT45 + Abs(lp_sbnow.y) * 2 / 5) * LHTRADB)) / Log(10#) 'r * Log( tan(45+y/2))
b1.yt = 1.9 * (Sgn(lp_sbnow.y) * cz - Sgn(yy) * sz)
End Sub
'円錐図法の計算
Private Sub z_LHTZensui(b1 As tdtype, ByVal xx&, ByVal yy&)
Dim r#, r0#, h#, f0#, yhos#
If lp_sbnow.y >= 0 Then
f0 = 30
If yy < -600000 Then yy = -600000
Else
f0 = -30
If yy > 600000 Then yy = 600000
End If
r0 = 1 / Tan(f0 * ATN1P45) '30度のとき
yhos = r0 - Tan((lp_sbnow.y * LHTRADB - f0 * ATN1P45))
r = r0 - Tan((yy * LHTRADB - f0 * ATN1P45))
h = ((xx - lp_sbnow.X) * LHTRADB) * Sin(f0 * ATN1P45)
b1.xt = r * Sin(h)
b1.yt = r * Cos(h) - yhos
End Sub
'ボンヌ図法の計算
Private Sub z_LHTZbonnu(b1 As tdtype, ByVal xx&, ByVal yy&)
Dim r#, r0#, h#, f0#, yhos#
If lp_sbnow.y >= 0 Then
f0 = 30
Else
f0 = -30
End If
r0 = 1# / Tan(f0 * ATN1P45) '30度のとき
yhos = r0 - (lp_sbnow.y - f0 * DOUBAI) * PLHT90 * Atn(1) * 2
r = r0 - (yy - f0 * DOUBAI) * PLHT90 * Atn(1) * 2
h = ((xx - lp_sbnow.X) * LHTRADB) * Cos(yy * LHTRADB) / r
b1.xt = r * Sin(h)
b1.yt = r * Cos(h) - yhos
End Sub