9-1投影図を活用する
9-1-1-1 正距方位図法に方位線を追加する
6-5-2-1 Azimuthal Equidistant (正距方位図法)で描画した投影図に方位線を
追加してみます。方位線は、外周と同じ式で計算できます。方位角に
対応した点と中心を結ぶと方位線となります。
・左の例では、方位を
10° にしているが必
要な方位角を設定す
ればよい。
・ラベルは、方位角毎
に表示されるので、
表示位置を調整すれ
ばよい。
[図
9-1] 正距方位図法 ( Azimuthal Equidistant Projection ) 方位線付
【地図主点(中心)】東経 135ﾟ､北緯35ﾟ 【経度間隔】15ﾟ【緯度間隔】10ﾟ
9-1-1-2 正距方位図法の簡易精度分析
正距方位図法は、地球を球とみなして投影しています。実際の地球に
近い回転楕円体として計算したものと比較して精度を確認して見ます。
<回転楕円体として計算させる>
WEB 上にはいろいろな計算サイトがありますが、多くのポイントを
一括で計算でき結果も利用しやすいので下記を利用します。
■国土地理院『 距離と方位角の計算 』
https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2stf.html
地図上で始点と終点を手入力するか、地図上をクリックして計算さ
せます。データポイントが多い場合は､一括でデータファイルを読み
込み、計算結果を出力できます。一括で読み込むデータポイントの
最大値は、30,000 ポイントまでです。計算式も提示されています。
使い方は、サイトの説明を読んでください。
<正距方位図法と回転楕円体として計算させたものを重ねて評価する>
細かい部分に差が見られるが、全体としては、あまり差がないように
見える。そこで、一部の地域を抜き出して見てみた。
オーストラリア南岸とニュージーランドに差があるように見える。
さらに中心点に近い部分(日本近海)を確認して見た。
日本分県地図用のデータはRjpWikiのShapeFileライブラリにあった ShapeFile
データ（地域境界・地域メッシュ）の 『 japank.zip 』 を利用した。
描画ポイントは 38,687 ポイントです。
【精度確認の結果】
球面三角法で計算された正距方位図法と回転楕円体として計算された
国土地理院 測量計算サイトの算出式はかなり複雑であるが、実際に
全地球をＡ４サイズ程度で描画する場合は、あまり差がないことが
判った。
2点間の距離については、WEB でいろいろ検討や説明されているサイトが
ありますが、説明が丁寧で分かりやすい以下を参考にしてみてください。
・グーグルマップは地球を半径６３７１㎞の
真球として距離計算を行っている
http://yamao.lolipop.jp/map/2017/g-m.htm
< 首都をマクロで変化させ正距方位図法の変化を見る >
正距方位図法の中心をマクロで世界の首都の経度緯度に書き換えて変化
を確認して見ます。
・首都のデータはアサヒ技研の h2706world_sjis.csv を利用した。
・参考都市として一部の日本の都道府県庁と北極、南極を追加した。
・VBA は、ループのカウンタ ｢I｣ をキーとして VLOOKUP関数を利用
して経度緯度を抽出した。
Sub 首都変化()
Dim I As Long
Dim LC As Long
LC = Range("Z6").Value
Application.Goto Reference:=Range("A1"), Scroll:=True
Range("A1").Select
For I = 1 To LC
Range("P11").Select
ActiveCell.FormulaR1C1 = "=VLOOKUP(" & I & ",R9C20:R300C26,6,)"
Range("P12").Select
ActiveCell.FormulaR1C1 = "=VLOOKUP(" & I & ",R9C20:R300C26,7,)"
Range("O8").Select
ActiveCell.FormulaR1C1 = "=VLOOKUP(" & I & ",R9C20:R300C26,4,)"
Range("O9").Select
ActiveCell.FormulaR1C1 = "=VLOOKUP(" & I & ",R9C20:R300C26,5,)"
Range("N4").Select
ActiveCell.FormulaR1C1 = "=VLOOKUP(" & I & ",R9C20:R300C26,2,)"
Range("N5").Select
ActiveCell.FormulaR1C1 = "=VLOOKUP(" & I & ",R9C20:R300C26,3,)"
Application.Wait [Now() + "0:00:00.40"] ‘インターバルタイム
Next I
End Sub
※この例では、YouTubeへのアップロードサイズをあまり大きくならない
ようにインターバルタイムを短くしています。そのため動画が、ややギグ
シャクしています。実際には、もう少し長くした方が良いでしょう。
経度緯度のヒョウジセルは、ユーザー定義でセットしています。
<セット例> 東経##0.00"度";"西経"##0.00"度"
北緯#0.00"度";"南緯"#0.00"度"
正距方位図法は、中心から離れるにつれ歪が増大してくるため、補助線
のピッチが粗いと、交差したり曲線にならなかったりする場合があります。
9-1-1-2 国の形の歪が少ない方位付き地図
前述した正距方位図法は、方位角の分かりやすさは良いのですが中心から
離れるにつれて形状が変形して来るため、どのエリアかどの国かがわかり
にくくなってしまいます。そこで、方位角の方位線を経度緯度の点として、
形状がわかる投影図に入力してプロットさせることにしてみました。
算出は球面三角法を利用します。
【方位角単位で到達点を求める】
出発点の緯度：φ1 到着点の緯度：φ2
出発点の経度：λ1 到着点の経度：λ2
方位角 ：α
方位角の計算式
⇒ tan α＝sin(λ2 -λ1) ÷ (cosφ1×tanφ2 - sinφ1 × cos(λ2-λ1))
●出発点を固定
●方位角 α をある角度毎に計算する
●到達点の経度 λ2 を －180°～ 180°で変化させる
●上記｢方位角の計算式｣を逆展開し φ2 を求める
【入力イメージ】
方位線の経度緯度を求めるため入力セルの配置などを画像にしました。
| ＜出発点 経度・緯度＞ | ||||||||||||||||||||||||||||||||||||
| 記号 | セルNo. | 記述式 | ||||||||||||||||||||||||||||||||||
| 緯度φ1 | J2 | ?? ← 出発点緯度 | ||||||||||||||||||||||||||||||||||
| 経度 Ｌ1 | J3 | ??? ← 出発点経度 | ||||||||||||||||||||||||||||||||||
| ＜データ計算 など＞ 経度緯度のデータ毎に記述式を複写(行コピー)する。 | ||||||||||||||||||||||||||||||||||||
| 記号 | セルNo. | 記述式 | ||||||||||||||||||||||||||||||||||
| C列 | C7～ | ??? ← 経度L2 ： 180°～ -180°までピッチ 1°で下方向に入力 | ||||||||||||||||||||||||||||||||||
| D行 | D6～ | ??? ← 方位角α ： 5°刻みで 180°まで横方向に入力 | ||||||||||||||||||||||||||||||||||
| 緯度 L2 | D7～ |
=DEGREES(ATAN((SIN(RADIANS($C7-$J$3))+TAN(RADIANS(D$6))* SIN(RADIANS($J$2))*COS(RADIANS($C7-$J$3)))/ TAN(RADIANS(D$6))/COS(RADIANS($J$2)))) |
||||||||||||||||||||||||||||||||||
| →この式を列方向、行方向に複写すると表として算出できます。 | ||||||||||||||||||||||||||||||||||||
方位角単位で経度(C列)、緯度(D列)をペアとして投影式に入力する。
9-1-1-2-1 Plate Carree Projection に方位線を加える
前述項で算出したデータを Plate Carree に投入して描画します。
但し、この図法は経線が直線となります。従って出発点と地球の
裏側となる到達点での方位角0°と180°は、経線と重なって直線
となりますが、前述算出方法では 1°ピッチとしているため出発点
から 0°の方位線は、1°ズレた点と結ばれ、やや傾きます。
気になる方は、出発点と到達点の算出は次のように算出します。
<出発点が 135度の例>
方位角180°の方位線の計算結果(投影式ではありません)
元DATA 改良DATA
: :
136 -90 136 -90
135 35 135 -90 ←追加
134 90 135 35
: 135 90 ←追加
134 90
:
[図 9-1-1-2-1] Plate Carree Projection ＋ 方位線
【地図主点(中心)】東経 150ﾟ､緯度 0ﾟ【経度間隔】15ﾟ【緯度間隔】10ﾟ
【出発点】東経135°、北緯35°
9-1-1-2-2 正射図法 に方位線を加える
前述項で算出したデータを 正射図法 に投入して描画します。但し、
方位180°(0°)は、そのまま投入しても長さが短くなってしまいました
(理由は未解析)ので補正しています。
この図法は、視点を出発点の地球外にしているため、当然放射線とな
っています。
[図 9-1-1-2-2] 正射図法 方位線付
【視点(光源)】地球から無限遠方より照射
【地図主点(中心)】東経 135ﾟ､北緯 35ﾟ 【主描画範囲】地図主点を中心として90ﾟ
【経度間隔】15ﾟ【緯度間隔】10ﾟ
9-1-1-2-3 方位線の到達地域を描画して見る
ここでは、描画図中に出発点を描画せず方位線の通過する地域を地方図として
描画してして見た。
9-1-1-2-3-1 ヨーロッパ地方への方位
<使用データ> 海岸線：ne_50m_coastline.shp
国 境：WDBII_border_i_L1.shp
首 都：h2706world_utf8.csv（2015年6月1日現在､アマノ技研）
アクセントの意味で色付けをしてみました。方法は『7-7 投影地図に色を付ける』
を参照してください。
[図 9-1-1-2-3-1] ヨーロッパへの方位
【方位中心】北緯 35°、東経135° 【描画図法】正距円錐図法
【中心経度】15° 【標準緯度１】45° 【標準緯度２】60°
9-1-1-2-3-2 近距離(この例では800キロ程度) に方位線を加える
①前述項で算出したデータを 正規多円錐図法 に投入して描画します。
[図 9-1-1-2-3-2-1] 青森県方位線付
【投影図法】正規多円錐図法 【方位の出発点】東経135°北緯35°(西脇市)
【描画中心経度】141°
描いてみて感じるのは、距離線があればということです。
出発点の緯度：φ1 目標到着点の緯度：φ2
出発点の経度：λ1 目標到着点の経度：λ2
< 目標到達距離 > = R acos (sinφ1 sinφ2 + cosφ1 cosφ2 cos ( λ2 - λ1))
で求められるのですが、結構な手間となります。そこで、描画外観の
歪が気になりますが、正距方位図法に戻って目標到達距離を同心円で
表示します。
【目標距離円を描く】
中心からの距離を示す同心円を描いてみる。
A：方位中心点
B：方位目標距離点
C：対蹠点(たいせきてん､地球の反対側の点)
O：地球の中心点
AO または OC：地球半径 R = 6,378,137m ⇦ ここでは赤道半径を利用しています。
∠AOB = α
∠AOC = 180°⇒ ラジアン = RADIANS(180) = π
●外周----対蹠点（たいせきてん）
円弧 AC ＝ πR ⇒ 円弧 AC／R ＝ π
↑
中心点からこの値の半径の円を描くと外周となる
●目標距離
円弧 AB ＝ αR ⇒ 円弧 AB／R ＝ α
↑
中心点からこの値の半径の同心円を描くと目標距離円となる
[例] 目標 700km
α ＝ 700 ÷ 6,378.137 ＝ 0.1097499
【入力イメージ】
[図 9-1-1-2-3-2-2] 青森県方位線･距離線付
【投影図法】正距方位図法 【方位の出発点】東経135°北緯35°(西脇市)
前項で書いた方法の距離円は、正距方位図法でしか書けません、他の地図投影図法
でも利用できる距離円を描けるデータを算出してみました。
●中心からの距離を求める式(球面三角法)
L = R acos (sinφ1 sinφ2 + cosφ1 cosφ2 cos (λ2-λ1))
変形して見る。
L/R = acos (sinφ1 sinφ2 + cosφ1 cosφ2 cos ( λ2-λ1))
cos(L/R) = sinφ1 sinφ2 + cosφ1 cosφ2 cos ( λ2-λ1)
ここで cosφ2=±√( 1 - sin2φ2) を代入
cos(L/R) = sinφ1 sinφ2 + cosφ1 (±√( 1 - sin2φ2)) cos ( λ2-λ1)
cos(L/R) - sinφ1 sinφ2 = cosφ1 cos ( λ2-λ1) (±√( 1 - sin2φ2))
平方根を取るため左右を二乗し、整理すると
(sin2φ1 +cos2φ1 cos2 ( λ2-λ1)) sin2φ2 – 2×cos(L/R) sinφ1 sinφ2 +
cos2(L/R) - cos2φ1 cos2 ( λ2-λ1)= 0
この式は sinφ2 の二次方程式となりその解が求める緯度となります。
この式でφ1､λ1、Rは既知、Lは距離毎にセットし、λ2は一定間隔で変化させれば
目標の緯度が算出できます。
私としては、50年以上前に勉強した「二次方程式の解」について再度履修することと
なりました。
ax2 + bx + c = 0
x = (-b ± √( b2 - 4ac)) / 2a というやつです。
a = sin2φ1 +cos2φ1 cos2 ( λ2-λ1)
b = - 2 × cos(L/R) sinφ1
c = cos2(L/R) - cos2φ1 cos2 ( λ2-λ1)
これをエクセルで計算して見ます。
< 入力イメージと数式 >
画像をクリックすると拡大表示されます。
＜出発点 経度・緯度＞
| 記号 | セルNo. | 記述式 | |
| 経度λ1 | C4 | ??? ← 出発点経度 | |
| 緯度φ1 | C5 | ?? ← 出発点緯度 | |
| R | C6 | 6,378.137 ← 地球半径(km) |
＜予備計算＞
| 記号 | セルNo. | 記述式 | |
| sinφ1 | F3 | =SIN(RADIANS(C5)) | |
|
sinφ12 |
F4 | =SIN(RADIANS(C5))^2 | |
| cosφ1 | F5 | =COS(RADIANS(C5)) | |
|
cosφ12 |
F6 | =COS(RADIANS(C5))^2 |
＜到達距離と計算＞
| 記号 | セルNo. | 記述式 | |
| 到達距離 L | B10 | ???? ← 到達距離 km | |
| cos(L/R) | C10 | =COS($B$10/$C$6) | |
|
cos2(L/R) |
D10 | =C10^2 |
＜データ計算 など＞ 経度のデータ毎に記述式を複写(行コピー)する。
| 記号 | セルNo. | 記述式 |
| No. | E10 |
??? ← 式の複写時､ソートや解析時に役立つので連番号を付与 |
| λ2 | F10 | ??? ← 変化させるの経度データ ※1 |
|
cos2(λ2-λ1) |
G10 | =COS(RADIANS(F10-$C$4))^2 |
|
cos2φ1*cos2(λ |
H10 | =$F$6*G10 |
|
a |
I10 | =$F$4+$F$6*G10 |
|
b |
J10 | =-2*C$10*$F$3 |
|
c |
K10 | =D$10-$F$6*G10 |
|
sinφ2① |
L10 | =(-J10+(J10^2-4*I10*K10)^0.5)/(2*I10) |
|
φ2① |
M10 | =ASIN(L10) |
|
度変換① |
N10 |
=IF(D10=1,DEGREES(M10),IF(ISERROR(M10), #N/A,DEGREES(M10))) ←※2 |
|
sinφ2② |
O10 | =(-J10-(J10^2-4*I10*K10)^0.5)/(2*I10) |
| φ2② | P10 | =ASIN(O10) |
|
度変換② |
Q10 |
=IF(D10=1,DEGREES(P10),IF(ISERROR(P10), #N/A,DEGREES(P10))) ←※2 |
※1 ：変化経度距離に対応して経度を変化させる。変化範囲が判らない場合は
-180 ～ 180°を変化させる。
※2 ：φ2が、エラーの時「#N/A」を返しているのは、グラフを作成したとき
「#N/A」は無視して表示でき、不要なデータを見ながら消去できるため。
●東経135°北緯35°を中心として、世界への距離円を算出
前述式にて算出する。距離に対応する経度が不明なので
-180 ~ 180°の範囲で変化させる。
(1)計算値をグラフにする。
かなり余分な距離円まで
算出されています。
(2)結果のデータから「#N/A」を消去する。
予期せぬ点が結ばれたのは、
削除されました。
線がとぎれている（スキマが
できる）のは、経度の変化ピ
ッチを1°にしたためで、とぎ
れた付近の経度ピッチを細か
くしてやればＯＫです。
(3)余分な部分を消去し、と切れたスキマ部分の経度ピッチを細かくして埋める。
余分な線を削除する式などを
考えたかったのですが、距離
円のピッチを 2,000km にし
た場合、円は10個程度ですの
で手加工で削除しました。
このデータを利用して投影図を描いてみる。
[図 9-1-1-3-1] Plate Carree Projection ＋ 方位線 ＋ 距離円
【地図主点(中心)】東経 150ﾟ､緯度 0ﾟ【経度間隔】15ﾟ【緯度間隔】10ﾟ
【出発点】東経135°、北緯35°
9-1-2 国土数値情報を活用する
第４章でデータの入手について触れたが、ここでは、国土地理院の
国土数値情報を使った沖縄の例について説明します。
【使用するデータ】
・海岸線と行政区----N03-17_47_170101.shp
『8-5 データを分解する』を元に海岸線と行政区
を区分けしておく
・河川--------------W05-07_47-g_Stream.shp
・湖沼--------------W09-05-g_Lake.shp(沖縄切出し)
・道路--------------N01-07L-2K_Road.shp
(旧 統一フォーマット形式)
※沖縄で切出しだが､エラーがあり手修正した
・バス路線----------N07-11_47.shp
・バス停名----------P11-10_47-jgd-g.xml
・鉄道--------------N02-15_RailroadSection.shp(沖縄切出し)
・鉄道駅名(位置)----N02-15.xml(沖縄切出し)
・高速道路----------N06-15_HighwaySection.shp(沖縄切出し)
・高速道路ＩＣ------N06-15_xml.xlsx(沖縄切出し)
・空港--------------C28-16_Airport.shp(沖縄切出し)
※空港の滑走路ではなく敷地データ
・空港滑走路--------mapli.net
http://www.mapli.net/position/
を利用して手で抽出
9-1-2-1 沖縄本島付近
前述データを沖縄本島付近に絞り込み正規多円錐図法で投影した。
9-1-2-2 今帰仁付近
沖縄本島より今帰仁付近(①)を抜き出しました。
ここで着目するのは、古宇利島です。地図上には、島として画がかれ
ています。有名な(？)古宇利大橋は描かれていません。勿論、海岸線
は測地をもとにしていますので島は、島としているのは当然なのです
が、橋がないのは不自然な感じです。
古宇利島から屋我地島方面を望む(沖縄県､東経128°1分16.04秒 北緯26度41'40.23")
< Link to Google Earth & Map : 26°41'40.23"N,128°1'16.04"E >
さらに屋我地島と奥武島に着目すると、この２島も島なのですが、バ
ス路線があるため沖縄本島と接続されているのが判ります。見かけ上
は、橋などは描画されず海上をバス路線で結ばれています。そして屋
我地島と本島今帰仁はワルミ大橋で結ばれているのですが、バス路線
(この時点のデータでは)がないので直結されていないような描画とな
っています。この様に橋がないと実際とイメージが合いません。しか
し国土数値情報には、橋のデータはあるものの中味が使い物になりま
せん。昨今、橋も長大で曲線になって来ており、是非、橋のデータの
提供を希望したいです。
9-1-2-3 那覇市付近
沖縄本島より那覇市付近(②)を抜き出しました。
沖縄唯一の鉄道モノレールであるユイレールや沖縄自動車道路､
那覇空港自動車道が判別できます。
データのところでも記述しましたが、国土地理院の空港は敷地
のデータなので、滑走路を別途用意して描画しています。やは
り滑走路有りと無しでは空港らしさが違いますよ。
9-1-2-4 利用データについて
那覇市付近より漫湖付近(③)を抜き出しました。
説明しやすくするため海と湖に色をつけています(着色方法は別項
参照)。
中央やや下部に着目して見ると『 河川 → 湖 → 海洋 』と並ぶ水系が
全国的にも珍しい地形です。実際に「漫湖」は、湖ではなく干潟なので
すが国土地理院の地図データとしては、明治橋を区切りとして海洋と
区切られています。余談ですが、国土地理院の海岸線データは、河川
との区切り点が、河口から少し上流の地点でバッサリと切られていま
す。例として東京都の荒川を見てみると、海岸線は、首都高速湾岸線
の「荒川河口橋」の上流 1.1km 程の地点でバッサリ切られています。
前述【使用するデータ】でも述べたのですが、道路データには一部エ
ラーがありました。内容は、那覇市付近を投影してみると道路がグチ
ャグチャになっている箇所がありました。ここでは、その部分を手修
正で削除したのですが、このデータには、根本的に問題があります。
バス路線とズレが発生しています。道路データの抽出メッシュの粗さ
が原因なら曲線と直線の差なのですが、明らかにバス路線からズレて
いるのです。バス路線は、許認可が必要でしょうから、ここではバス
路線の方が正しいと思われます。
国土数値情報を探したのですが道路データしは、これしか見つけられ
ませんでした。
