EXCELで地図を描く

概　要

　海岸線データをもとに地図を描画するVBA
　エクセル上に地図を描画できるが、それ以上でもそれ以下でもない。
　下の画像は国土地理院の海岸線データを基にメルカトル図法で描画した日本
　

海岸線データをエクセル上で描画する

　海岸線のデータは後述
　西経３０度を端とするメルカトル図法で描画する方法です。
　ＧＲＳ８０１回転楕円体として計算しています。
　メルカトル図法のため南北緯８０度以上は計算していません。※南北端に近づくほど無限に東西に引き延ばされるため
　Minの値をいじることで縮尺を変更できます。が、エクセル上のXY座標の限界を突破した場合は描画できません

Option Explicit
Const Min = 10000 '1万分の1
Const a = 6378136.6 / Min
Const e = 8.18191910428158E-02


Function Atanh(Agmt)
    '双曲線アークタンジェント関数
    Atanh = Log((1 + Agmt) / (1 - Agmt)) / 2
End Function

Function LatCnvY(φ)
    LatCnvY = a * Atanh(Sin(φ * 4 * Atn(1) / 180)) - a * e * Atanh(e * Sin(φ * 4 * Atn(1) / 180))
'エクセル対応用
    LatCnvY = (a * Atanh(Sin(80 * 4 * Atn(1) / 180)) - a * e * Atanh(e * Sin(80 * 4 * Atn(1) / 180))) - LatCnvY
End Function

Function LonCnvX(λ)
    If λ <= -30 Then λ = λ + 360
    LonCnvX = λ * a * 4 * Atn(1) / 180
    'エクセル対応用
    LonCnvX = (30 * a * 4 * Atn(1) / 180) + LonCnvX
End Function
Sub メルカトル図法()
'***********************************************************************************************************
'メルカトル図法
'***********************************************************************************************************
'--- 記号表 ----------------------------------------------------------------------
' R：真球地球の半径
'λ：中央経線からの経度差
'中央経線：地図の中心となる線
'φ：緯度

'--- 地球の半径 ------------------------------------------------------------------
'/////////////////////////////////////////////////////////////////////////////////
'///GRS801楕円体                                                              ///
'/////////////////////////////////////////////////////////////////////////////////
'現世界の測地系で最も広く使われている
'a：長半径(赤道半径)：6,378,136.6±0.1m
'b：短半径(極半径)： 6356751.91548148m
'e：離心率：(√a^2-b^2)/a：0.081819191042815790
'f：扁平率：1/298.257222101
'f：扁平率：円又は球に比べどれだけつぶれているかを表す値
'/////////////////////////////////////////////////////////////////////////////////
'/// WGS84                                                                     ///
'/////////////////////////////////////////////////////////////////////////////////
'f：扁平率：1/298.257223563
'GRSとの差異は、地球の短半径(極半径)にすると約0.105mm異なるが、実用上問題はない。

'/////////////////////////////////////////////////////////////////////////////////

'--- x座標 -----------------------------------------------------------------------
'x = Rλ
'赤道を標準緯線とする接円筒図法に共通
'球、楕円体共通式
'λ = R/x

'--- y座標 -----------------------------------------------------------------------
'緯度φの緯線のy座標を表す式（球）
'y = R * log[tan{(φ/2)+(π/4)}]
'  = (R/2) * log{(1+sinφ)/(1-sinφ)}
'  = R * tanh^(-1)[sinφ]
'

'緯度φの緯線のy座標を表す式（楕円体）
'y = alog{[tan{(φ/2)+(π/4)}][(1-esinφ)/(1+sinφ)]^(e/2)}
'  = a * tanh^(-1)[sinφ] - a * e * tanh^(-1)[e * sinφ]

'***********************************************************************************************************

    'Dim b, e, f, x, y, φ, λ
    'b = 6356751.91548148
    'e = 8.18191910428158E-02
    'f = 1 / 298.257222101
    'φ = 35
    'λ = 140 - 135
    'x = a * λ
    'y = a * Atanh(Sin(φ)) - a * e * Atanh(e * Sin(φ))
    
    'Debug.Print "x:" & x & ",y:" & y

'***********************************************************************************************************
'地図の描画
'***********************************************************************************************************
'CSVファイルの読み取り
   Dim Path, CsvFile
   Dim Area(), Cnt
   Dim Lat, Lon
   
   Path = Application.GetOpenFilename("CSVファイル,*.csv?")
 
   Open Path For Input As #1
   Do Until EOF(1)
      Line Input #1, CsvFile
   Loop
   Close #1
   
   Cnt = 0
   Do Until Not InStr(1, CsvFile, "_") <> 0
      ReDim Preserve Area(Cnt)
      Area(Cnt) = Mid(CsvFile, 1, InStr(1, CsvFile, "_") - 1)
      CsvFile = Mid(CsvFile, InStr(1, CsvFile, "_") + 1)
      Cnt = Cnt + 1
   Loop
   
   Cnt = Cnt - 1
   
'地図の描画
   Do Until Cnt = 0
      '経度
      Lon = Mid(Area(Cnt), 1, InStr(1, Area(Cnt), ",") - 1)
      Lon = LonCnvX(Lon)
      Area(Cnt) = Mid(Area(Cnt), InStr(1, Area(Cnt), ",") + 1)
      '緯度
      Lat = Mid(Area(Cnt), 1, InStr(1, Area(Cnt), ",") - 1)
      If Lat > 80 Then Lat = 80
      Lat = LatCnvY(Lat)
      Area(Cnt) = Mid(Area(Cnt), InStr(1, Area(Cnt), ",") + 1)
      
      With ActiveSheet.Shapes.BuildFreeform(msoEditingAuto, Lon, Lat)
           Do Until Area(Cnt) = ""
              '経度
              Lon = Mid(Area(Cnt), 1, InStr(1, Area(Cnt), ",") - 1)
              Lon = LonCnvX(Lon)
              Area(Cnt) = Mid(Area(Cnt), InStr(1, Area(Cnt), ",") + 1)
              '緯度
              If InStr(1, Area(Cnt), ",") <> 0 Then
                 Lat = Mid(Area(Cnt), 1, InStr(1, Area(Cnt), ",") - 1)
              Else
                 Lat = Area(Cnt)
                 Area(Cnt) = ""
              End If
              If Lat > 80 Then Lat = 80
              Lat = LatCnvY(Lat)
              '描画
              Area(Cnt) = Mid(Area(Cnt), InStr(1, Area(Cnt), ",") + 1)
              .AddNodes msoSegmentLine, msoEditingAuto, Lon, Lat
           Loop
           .ConvertToShape.Select
      End With
      
      Cnt = Cnt - 1
   Loop
'***********************************************************************************************************
End Sub

図形の描画に関しては以下の通り
With ActiveSheet.Shapes.BuildFreeform(msoEditingAuto, X座標, Y座標)で指定する部分が始点
With中の.AddNodes 部分が中点～終点で、複数設定すると頂点を増やせる。
繰り返し処理でも頂点は増やすことが可能
始点と終点を同じにすると自動で閉じた図形に。

      With ActiveSheet.Shapes.BuildFreeform(msoEditingAuto, X座標, Y座標)
           .AddNodes msoSegmentLine, msoEditingAuto, X座標, Y座標
           .AddNodes msoSegmentLine, msoEditingAuto, X座標, Y座標
           .ConvertToShape.Select
      End With

緯度経度をエクセル上のXY座標に変換する式については、コメントアウトしている部分とfunctionのとおり
エクセル上のXY座標の確認についてはこちらの記事を参考程度にどうぞ

海岸線データダウンロード

　こちらのページから全レイヤをダウンロードする。
※下図赤枠

海岸線データをCSVへ変換

解凍したデータの中のcoastl_jpn.shpをQGISに取り込む。

レイヤ上で右クリックし名前を付けて保存を選択

形式はCSVを選択し、下図のように設定する。※WGS84を選択
保存場所と保存名はなんでもいい。

データの編集

作成したデータをメモ帳で開く。

「一行目」と「"LINESTRING (」を削除
「半角スペース」を「,（カンマ）」に変換
「)",BA010,1,1,JPN」と「)",BA010,3,1,JPN」を「 _ 」に変換
一番最後の「 _ 」を削除
し上書き保存します。
※１：「 _ 」の空白部分は含めない。
※２：「 _ 」ごとに一つの図形として取り扱います。地理院のデータは本州や九州、北海道など一本の線で描かれていません。
次に改行コードを置換できるソフトで開き、改行コードを削除する。
ここではAtomを使用
Atomの場合、「 * 」を選択し、「\r\n」を検索窓に入れると置換できる。
これで読み込みデータは完成。

さいごに

　図形を描けるが現在のところそれ以上の応用は難しそう。
　エクセル上のXY座標を取得する関数でもあれば違うのだろうが……

参　考

地図投影法――地理空間情報の技法　著者：政春　尋志（国土地理院基本図情報部長）