地図ソフトの作り方
FLand-Japan&World MapProgram
地図ソフトの作り方について解説。地図ソフト「FLand-Ale日本世界地図」関連情報。
地図ソフトの作り方
地図ソフト作成の参考となるプログラムを紹介します。各項目ではFLand-Aleのソースとコメントがありますので参考にしてください。
- 距離と方位の計算
- 地図座標変換
- 地図データ
- 地名データ
- 地図表示の高速化
- モジュール構成
- クラス構成
距離と方位計算
- ここでは２つの処理を行います。1)基準地点から距離と方位を指定して、目的位置を求める。
2)基準地点と目的地点を指定して、距離と方位を求める。これらの計算は球面三角法により行っています。Private Const mATN1P45 = 0.0174533 'atn(1)/45
Private Const ER = 6378.14
Private Const PI = 3.14159265358979
'LHT座標値例
'北緯36度30分 → 365000
'南緯40度45分 → - 407500
'北緯136度30分 → 1365000
'南緯140度45分 → -1407500
Option Explicit
'位置と距離と方位を指定して目的位置を求める
Public Sub GetPos_Kyorihoui(ByVal Mypx&, ByVal Mypy&, ByVal Kyori#, ByValHoui#, Topx&, Topy&)
'Mypx 基準位置経度 Mypy 基準位置緯度 LHT座標
'Kyori 距離 km Houi 方位 度
'Topx 目的位置経度 Topy 目的位置緯度 LHT座標
'方位は0-359.99.. 北が0、東が90、南が180、西が270となる。
Dim cb#, sb#, clc#, c#
Dim xx#, ido#, kdo#
ido = (900000 - Mypy) * (Atn(1) * 4) / 1800000
c = Kyori / ER
If c >(Atn(1) * 4) * 2 Then c = c - (Atn(1) * 4) * 2
cb = Cos(c) * Cos(ido) + Sin(c) * Sin(ido) * Cos(Houi * mATN1P45)
If cb = 1 Then
sb = 0
Else
sb = Sqr(1 - cb * cb)
End If
If Sin(ido) * sb = 0 Then
clc = 0
Else
clc = (Cos(c) - Cos(ido) * cb) / (Sin(ido) * sb)
If clc >1 Then clc = 1
If clc < -1 Then clc = -1
End If
xx = cb
If Sqr(-xx * xx + 1) = 0 Then
ido = Atn(1) * 2
Else
ido = Atn(-xx / Sqr(-xx * xx + 1)) + 2 * Atn(1)
End If
Topy = 900000 - ido * 1800000 / (Atn(1) * 4)
xx = clc
If xx = 1 Then
kdo = 0
ElseIf xx = -1 Then
kdo = Atn(1) * 4
Else
If Sin(Houi * mATN1P45) >= 0 Then
kdo = Atn(-xx / Sqr(-xx * xx + 1)) + 2 * Atn(1)
Else
kdo = -(Atn(-xx / Sqr(-xx * xx + 1)) + 2 * Atn(1))
End If
End If
If c >Atn(1) * 4 Then kdo = -kdo
kdo = kdo * 1800000 / (Atn(1) * 4)
Topx = Mypx + kdo
If Topx >1800000 Then Topx = Topx - 3600000
If Topx < -1800000 Then Topx = Topx + 3600000
End Sub
'2点の位置を指定して距離と方位を求める。
Public Sub GetKyoriHoui(ByVal Mypx&, ByVal Mypy&, ByVal Topx&, ByVal Topy&,Kyori#, Houi#)
'Mypx 位置１経度 Mypy 位置２緯度 LHT座標
'Topx 位置２経度 Topy 位置２緯度 LHT座標
'Kyori 計算結果距離 km Houi 計算結果方位 度
Dim xx#, ram#, sig#, ssig#
Dim mx#, my#, tx#, ty#
On Error GoTo getkher
'２点の位置を設定
mx = Mypx
tx = Topx
my = Mypy
ty = Topy
If mx = tx And my = ty Then '010711
Kyori = 0
Houi = 0
Exit Sub
End If
'ラジアン変換
mx = mx * PI / 1800000
my = my * PI / 1800000
tx = tx * PI / 1800000
ty = ty * PI / 1800000
'距離を求める
ram = tx - mx
xx = Sin(my) * Sin(ty) + Cos(my) * Cos(ty) * Cos(ram)
If xx = -1# Then
sig = PI
ElseIf xx = 1# Then
sig = 0
Else
sig = Atn(-xx / Sqr(-xx * xx + 1)) + PI / 2
End If
Kyori = ER * sig
'方位を求める
ssig = Sin(sig)
If ssig = 0 Then
Houi = 0
Else
xx = Cos(ty) * Sin(ram) / ssig
If -xx * xx + 1 < 0 Then
If Sin(ram) >0 Then
Houi = 90
Else
Houi = 270
End If
ElseIf Sqr(-xx * xx + 1) = 0 Then
If Sin(ram) >0 Then
Houi = 90
Else
Houi = 270
End If
Else
Houi = Atn(xx / Sqr(-xx * xx + 1))
If 0 >Cos(my) * Sin(ty) - Sin(my) *Cos(ty) * Cos(ram) Then
Houi = PI - Houi
End If
Houi = Houi * 180 / PI
If Houi < 0 Then Houi = Houi + 360
End If
End If
Exit Sub
'参考
'Arcsin(x) = Atn(x / Sqr(-x * x + 1))
'Arccos(X) = Atn(-X / Sqr(-X * X + 1)) + 2 * Atn(1)
getkher:
'Errmes Err, "s_Getkyorihoui"
Resume Next
End Sub
地図座標変換
- 地図表示の基本となる部分です。地図データは緯度経度となっていますが、画面上に表示される場合は、地図の図法により変換作業を行うことになります。以下のプログラムでは１点のポイントの緯度経度を画面上の(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