世界測地系の地図上の日本測地系3次メッシュポリゴンを作る
作業の流れ
- 緯度経度による定義通りのメッシュポリゴンをつくる(Excel)
- 日本測地系の緯度経度を世界測地系に変換する(TKY2JGD http://vldb.gsi.go.jp/sokuchi/program.html )
- 世界測地系の座標を十進度数表示に戻し,Shapefileに変換(Excel,R http://www.r-project.org/ )
- ShapefileをフリーのGISソフトで表示(Quantum GIS http://www.qgis.org/ ,基盤地図情報閲覧コンバートソフト http://fgd.gsi.go.jp/download/ )
- 生物記録の地理データベース構築
Excel以外はフリーウェアです。それぞれダウンロードしてインストールが必要です。
かなりちからわざですので,もっとスマートなやり方があれば教えてください。
メッシュポリゴンデータの作成(Excel2003使用)
ファイルの作成
ワークシート3枚のExcelファイルを作ります。
マクロで使う関係で,ワークシート名を次のように変えます。
Sheet1 → m2_sel
Sheet2 → m3_ds
Sheet3 → m3_dd
3次メッシュポリゴンを作る範囲の2次メッシュを選択
m2_selのシートにA1のセルから下へ順に3次メッシュを作りたい場所の2次メッシュコード(6桁数字)を入れます。Excel 2003の場合,入力する2次メッシュの数は131以下にしてください。2次メッシュ1つで500行のデータができるので,65536行を越えてしまいます。
福岡県全域をカバーしたい場合はこんな感じ。環境庁のメッシュマップで2次メッシュ87個あります。
513020 | |
513005 | |
513006 | |
503075 | |
503076 | |
503077 | |
503170 | |
503060 | |
503063 | |
… | |
493042 | |
493043 | |
493044 | |
493033 | |
493034 |
ExcelVBAでマクロコードを入れる
Excelのツール→マクロ→Visual Basic Editorで新しいモジュールを追加して以下のコードを入れます。
Sub mesh3_generate() Dim m2_sel As Worksheet Dim m3_ds As Worksheet Dim m3_dd As Worksheet Set m2_sel = ActiveWorkbook.Worksheets("m2_sel") Set m3_ds = ActiveWorkbook.Worksheets("m3_ds") Set m3_dd = ActiveWorkbook.Worksheets("m3_dd") cur_row = 2 cur_obj = 1 cur_m2 = 1 m3_ds.Cells(1, 1) = "ID" m3_ds.Cells(1, 2) = "X" m3_ds.Cells(1, 3) = "Y" m3_ds.Cells(1, 4) = "Lat" m3_ds.Cells(1, 5) = "Lon" m3_dd.Cells(1, 1) = "ID" m3_dd.Cells(1, 2) = "Mesh3Code" Do While m2_sel.Cells(cur_m2, 1) > 0 m2_c = Format(m2_sel.Cells(cur_m2, 1), "000000") m1_lat_c = Val(Left(m2_c, 2)) m1_lon_c = Val(Mid(m2_c, 3, 2)) m2_lat_c = Val(Mid(m2_c, 5, 1)) m2_lon_c = Val(Right(m2_c, 1)) m2_lat = m1_lat_c * 2 / 3 + m2_lat_c / 12 m2_lon = 100 + m1_lon_c + m2_lon_c / 8 For i = 0 To 9 For j = 0 To 9 m3_lat0 = m2_lat + i / 120 m3_lon0 = m2_lon + j / 80 m3_lat1 = m3_lat0 + 1 / 120 m3_lon1 = m3_lon0 + 1 / 80 m3_c = m2_c & Format(i, "0") & Format(j, "0") m3_lat0_dms = Int(m3_lat0) * 10000 + (Int(m3_lat0 * 60) - Int(m3_lat0) * 60) * 100 _ + m3_lat0 * 3600 - Int(m3_lat0 * 60) * 60 m3_lon0_dms = Int(m3_lon0) * 10000 + (Int(m3_lon0 * 60) - Int(m3_lon0) * 60) * 100 _ + m3_lon0 * 3600 - Int(m3_lon0 * 60) * 60 m3_lat1_dms = Int(m3_lat1) * 10000 + (Int(m3_lat1 * 60) - Int(m3_lat1) * 60) * 100 _ + m3_lat1 * 3600 - Int(m3_lat1 * 60) * 60 m3_lon1_dms = Int(m3_lon1) * 10000 + (Int(m3_lon1 * 60) - Int(m3_lon1) * 60) * 100 _ + m3_lon1 * 3600 - Int(m3_lon1 * 60) * 60 'メッシュの4隅の座標値。シェープファイル用に緯度経度の順序入れ替わりに注意 m3_ds.Range(m3_ds.Cells(cur_row, 1), m3_ds.Cells(cur_row + 4, 1)) = cur_obj m3_ds.Cells(cur_row, 2) = m3_lon0 m3_ds.Cells(cur_row, 3) = m3_lat0 m3_ds.Cells(cur_row + 1, 2) = m3_lon0 m3_ds.Cells(cur_row + 1, 3) = m3_lat1 m3_ds.Cells(cur_row + 2, 2) = m3_lon1 m3_ds.Cells(cur_row + 2, 3) = m3_lat1 m3_ds.Cells(cur_row + 3, 2) = m3_lon1 m3_ds.Cells(cur_row + 3, 3) = m3_lat0 m3_ds.Cells(cur_row + 4, 2) = m3_lon0 m3_ds.Cells(cur_row + 4, 3) = m3_lat0 'TKY2JGD用 m3_ds.Cells(cur_row, 4) = m3_lat0_dms m3_ds.Cells(cur_row, 5) = m3_lon0_dms m3_ds.Cells(cur_row + 1, 4) = m3_lat1_dms m3_ds.Cells(cur_row + 1, 5) = m3_lon0_dms m3_ds.Cells(cur_row + 2, 4) = m3_lat1_dms m3_ds.Cells(cur_row + 2, 5) = m3_lon1_dms m3_ds.Cells(cur_row + 3, 4) = m3_lat0_dms m3_ds.Cells(cur_row + 3, 5) = m3_lon1_dms m3_ds.Cells(cur_row + 4, 4) = m3_lat0_dms m3_ds.Cells(cur_row + 4, 5) = m3_lon0_dms 'メッシュコード m3_dd.Cells(cur_obj + 1, 1) = cur_obj m3_dd.Cells(cur_obj + 1, 2) = m3_c cur_row = cur_row + 5 cur_obj = cur_obj + 1 Next j Next i cur_m2 = cur_m2 + 1 Loop m3_ds.Range("B:C").NumberFormatLocal = "0.00000000_ " End Sub
さて,このマクロを実行すると,m3_dsとm3_ddのワークシートにずらっと数値が並びます。
m3_ds
ID | X | Y | Lat | Lon |
---|---|---|---|---|
1 | 130.00000000 | 34.16666667 | 341000 | 1300000 |
… |
m3_dd
ID | mesh3code |
---|---|
1 | 51302000 |
… |
とりあえずエクセルファイルをm3_generate.xlsとか名前を付けて保存しておきます。
この時点のデータを後述のシェープファイルへの変換を行うと世界測地系のメッシュができます。m3_dsのワークシートのA〜Cの3列をコピーしてテキストエディタに貼り付け,m3W_ds.txtとして,m3_ddのデータの入っている2列を同様にしてm3_dd.txtとして保存しておきます。