地図ソフトの作り方
地図ソフト作成の参考となるプログラムを紹介します。各項目では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#, ByVal Houi#, 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
地図表示の高速化
FLand-Aleは地図表示の高速化を以下のように行っています。 １．ブロック化 地図データを緯度経度のブロック化し地球全体で２７分割、さらに日本部分を１１分割する。そして地図表示時は以下の処理を行う。 1)地図表示範囲に含まれるブロック番号を求める2)該当するブロックの地図データのみ表示を行う。
3)海岸線・標高線で複数ブロックにまたがり、一部のブロックのみ表示する場合、表示しないブロック部分はそのブロックの淵を経由するデータにする。 ２．データの間引き 地図の縮尺・表示精度によりベクトルデータの間引きを行う。 'ポイント間隔を設定する
If mSpmap Then '詳細地図か？
sph = PSTEPSP
Else
sph = PSTEPNR
End If
If HyojiStyle = HS_NORMAL Then
m_Pstep1 = (PSTEP1 / .picbai) / sph '縮尺により連動
m_Pstep2 = (PSTEP2 / .picbai) / sph
Else
m_Pstep1 = (PSTEP1B / .picbai) / sph
m_Pstep2 = (PSTEP2B / .picbai) / sph
End If
Select Case mMapSeido '精度により修正する
Case 0: mms = MMS0
Case 1: mms = MMS1
Case 2: mms = MMS2
End Select
m_Pstep1 = m_Pstep1 * mms
m_Pstep2 = m_Pstep2 * mms ３．計算の高速化 (1)三角関数で固定値となるものは定数で行う。 Private Const mATN1P45 = 0.0174533 'atn(1)/45 (2)浮動小数計算はなるべく行わない。整数化できるものは整数化してから計算を行う。
４．同じ種類ごと処理 市町村などアイテムの表示はそれぞれアイテムごとに描画する。ピクチャーボックスの色・フォントサイズなど変更に時間がかかるため。 For kakujun = KAKU_END To 0 Step -1 '島から描いていく
If mPshskub(kakus) Then '961121
z_Fontset kakus, mPshskubMp(kakujun), mMappic '種別ごとにフォントと色を設定する。
Select Case mPicprn
Case 1
Set prx = New CPolyPrn
End Select
For j = 0 To ii 'すべてのアイテムから
If mPsno(j).kaku = kakus Then '今回処理の種別なら
z_Psetshspp_pic j, mPicprn, mPshskubMp(kakujun)
End If
Next
ウィンドウ部の構成
ウィンドウ部のモジュール構成の概要です。 Form=Map.frm メインウィンドウForm=About.frm Ａｂｏｕｔウィンドウ
Form=MapOc.frm 作業用
Class=CfmCboF コンボボックス処理
Class=CUser ユーザー処理
Class=CPass パスワード処理
Class=CfmSBar スライドバー処理
Class=CfmToolB ツールバー処理
Class=CfmSts ステータスバー処理
Class=CfmAirPanel フライトシミュレータパネル処理
Class=CfmPicSep 画面区切り処理
Class=CFlaDDE ＤＤＥ処理
Class=CKubMp 作業用
地図作成部の構成
地図作成部のモジュール構成の概要です。 Form=paletto.frm パレット設定ウィンドウForm=Print.frm プリント設定ウィンドウ
Form=picWork.frm 作業用
Form=Mapshsoc.frm 作業用
Form=Onmap.frm 地図ツールウィンドウ
Class=clsChooseColor 色選択
Class=clsGetOpenFile ファイル選択
Class=clsChooseFont フォント選択
Class=CKyoho 距離方位計算
Class=CLhtPic 緯度経度-画面座標変換
Class=Cpal パレット
Class=CMapShsMk 地名の描画
Class=CMapShs 地図作成メイン
Class=CMapShsMp 地名の描画
Class=CMapShsL ライン表示
Class=CPolyPrn 地形描画
Class=Cprintset プリンタ設定
Class=CListFit リスト地名の地図範囲計算
Class=COnmap 地図ツールメイン
Class=COnMapLbl 地図ツールラベル
Class=COnMapLine 地図ツール方位線等距離線
地形処理部の構成
地形処理部のモジュール構成の概要です。 Form=Paletto.frm パレット設定ウィンドウForm=picWork.frm 作業用
Form=Dirdism.frm 距離方位ウィンドウ
Form=AutoMap.frm 地図自動設定ウィンドウ
Form=MapItem.frm 地図アイテム設定ウィンドウ
Class=CPolyPrn 地形描画
Class=CAirPlane フライトシミュレータパネル
Class=CAirCont フライトシミュレータ操作
Class=CZahyo 各種座標変換
Class=CKyoho 距離方位計算
Class=CLhtPic 緯度経度画面座標変換
Class=CMapIn 地図範囲取得
Class=CMpm 地図パラメータ
Class=CMapMake 地図作成
Class=CMapMouse マウスによる地図移動
Class=CMpmstr 図法距離線表示
Class=CPi_Fly フライトシミュレータ
Class=CPi_SubM ガイド地図
Class=CMapMove 地図移動
Class=CMapinUS 地図範囲取得USGS
Class=CMpmSmth スムーズ地図
Class=CMapKeido 経緯度線
Class=CMapComm 地図共通
Class=CMapFill 地図背景
Class=CMapBar 地図スクロールバー
Class=CMapSave 地図保存
Class=CMapA 地図作成メイン
Class=CMapHyojiOp 起動時地図表示
Class=CMpmF 地図パラメータファイル
Class=Cpal 地図パレット
Class=CDirDisM 距離方位計算
Class=CAutoMap 自動地図
Class=CMapItem 地図アイテム
地名処理部の構成
地名処理部のモジュール構成の概要です。 Form=Shssel.frm 地名検索ウィンドウForm=Item.frm 地名アイテム設定ウィンドウ
Class=CfmTview ツリービュー
Class=CfmIml イメージリスト
Class=CItemno アイテム
Class=CShsA 地名メイン
Class=CShsData 地名データ
Class=CShSeSel 地名検索
Class=CShsFSys 地名ファイルシステム
Class=CShsLview リストビュー
Class=CShsList 地名リスト
Class=CShsFa 地名ファイルメイン
Class=CShsFCsv CSVファイル読み込み
Class=CShsFFla ＦＬＡファイル読み込み
Class=CShSeName 地名検索
Class=CshsFMake 地名ファイル作成
Class=CShsFsub 地名ファイルサブ
Class=CShsFEx その他のファイル読み込み
Class=CShsFHtm HTMLファイル読み込み
Class=CShSeNear 近隣地名検索
Class=CZahyo 座標変換
Class=CMyData マイデータ
Class=CItem アイテム
Class=CShsInMap 地名の地図範囲内計算
Class=CShsItemk アイテム
Class=CShsFGps GPSファイル読み込み
Class=CShsJhyoji アイテム情報表示
Class=CShsPop アイテム情報表示
Class=CShSeA 地名検索処理メイン
Class=CShsStr アイテム名称
Class=CShsFKok 作業用
Class=CShsFHtmS HTMLファイル読み込みサブ
モジュール構成
FLand-Aleプログラムは大分類すると４つのコンポーネントで構成されます。- ウィンドウ部
- 地図表示部
- 地形処理部
- 地名処理部
StartupProject=..\vbp\flandc.vbp －－ ウィンドウ部
Project=..\vbpocx\FLandBox.vbp －－ 地図表示部
Project=..\vbpocx\FLandMp.vbp －－ 地形処理部
Project=..\vbpocx\FLandGn.vbp －－ 地名処理部
１．ウィンドウ表示・操作部
ウィンドウ表示・操作部はメニュー、ツールバー、地名リスト表示、キー入力、マウス入力などの書く処理を行います。
操作指令により処理を行いますが、ほとんどの処理は地図表示部へ操作内容を伝え、地図表示部が処理を行います。
２．地図表示部
地図表示部は、ウィンドウ表示・操作部から操作内容を受け、地図画像の表示を行います。また、操作のうち地形表示に関するものは地形処理部へ、地名に関するものは地名処理部へ各作業を依頼します。
地図の描画指令
1)地形処理部へ地形描画作業を依頼する。
2)地名処理部へ地図の描画範囲の地名一覧リストの作成を依頼する。
3)出来上がった地図の上に地名一覧リストの地名を表示する。
4)地形表示部へ図法表示などの描画作業を依頼する。
３．地形処理部
地形処理部は地図に関する処理を行います。
(1)起動時に地図データを読み込む。
(2)地図の描画指令が来た場合、地図を表示する。
４．地名処理部
地名処理部は地名に関する処理を行います。
(1)起動時に地名データを読み込む。
(2)地名リストを要求条件により作成し、画像表示部へ渡す。
地図データ
１．内容海岸線、各標高ごとの標高線をベクトルデータとしています。
２．入手先 十数年前になりますが、USGS（米国地質調査所）の30秒メッシュ地図を入手しました。当時はUSGSのページに登録すると無料でデータが入手可能でした。申し込みのページ（もちろん英語）に住所氏名などを記入して２週間ほどしてCD-ROMが６枚とどきました。
３．FLand-Ale用データ作成
USGSのデータを入手してもそのままではFLand-Aleに使用できません。FLand-Aleで使用するためにはベクトルデータ化をします。メッシュ→ベクトル変換のプログラムを作成し、地図上で任意の地点をクリックするとその周りの同一標高のデータを拾いだしベクトルデータを作成するものです。
大陸の標高が低いところでは１地点入力するとかなり広範囲のデータが作成できるのですが、標高の高いところや島などは１点１点指定してもできるデータはわずかなものです。手間をかけて多くのデータを作りました。
４．データ形式
mpz() As POINTAPI
単純ですが、緯度経度をLong(32bit)化しています。
地名データ
１．内容日本都市は「都道府県市町村区」のすべて、世界都市は主な都市が含まれています。また、山、駅、島、鉄道路線などのデータもあります。各地名ごとに地図表示のレベルなどの付加情報があります。
２．入手先
日本の都市データは現在では各所に詳細なデータがありますが、FLand-Aleの製作を始めたときは使用できるデータが見つかりませんでした。そこですべて手入力で入力し、緯度経度は地図を表示し、マウスで指定すれば緯度経度がデータとして作成できるツールを作成し、一つ一つ位置データを作成しました。
世界の都市データも同様です。こちらは現在でも有効なデータはなかなか見つかりません。
３．FLand-Ale用データ作成
各データに種別を付加します。これはエクセルなどで記入していきました。表示レベルなどは周囲との関係を考えつつ調整を行いました。
４．データ形式
一つの地名データは以下の８つの項目からなっています。
Public Type SH1type '80BYTE 4fH
sname As String * 32 '名称
syomi As String * 32 'めいしょう 県INDEXの場合はファイル名
px As Long '経度 県INDEXの場合は開始と終了
py As Long '緯度
kubun As Byte '区分
plevel As Byte '表示レベル 県INDEXの場合はTopLayer,'地図の場合は図法
vsize As Integer '4aH-4bH 図形の場合のサイズ,山の場合は標高,
imgの場合はpicf番号,'KUB_LINEの場合はTOPとの差分
'地図の場合は倍率/10
',県indexの場合は日本のusgs番号
code As Long '4cH-4fH 線、図形の場合は色、
'地図の場合は上位から俯角(7bit)、方位(9bit)、IDOH、KEIW
',県indexの場合はusgs番号
End Type
地図座標変換
地図表示の基本となる部分です。地図データは緯度経度となっていますが、画面上に表示される場合は、地図の図法により変換作業を行うことになります。以下のプログラムでは１点のポイントの緯度経度を画面上の(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