2011-06-11

日本測地系のメッシュを世界測地系の地図で使う

GIS

twitterで，環境省の自然環境保全基礎調査などで使われている3次メッシュについて話題になりました。その場では，過去の調査データを活かすために，日本測地系のメッシュを使い続ける必要があるという意見で一致しました。
しかし，世界測地系への移行から10年近くたって，今入手できる地理情報はほとんどが世界測地系になっています。地域メッシュコード（JIS X 0410）についての説明を探すと，日本測地系によるメッシュの有効期限は来年2月までという記述もありました。
（ http://www.stat.go.jp/data/mesh/yougo.htm ）
自然史に関わるデータのように，過去の記録が重要な分野では，記録の集計方法の変更により連続性が分断されたり，そもそもメッシュコードが日本測地系なのか世界測地系なのかわからなくなることでデータの正確性が損なわれることが問題です。
GISが簡単に使えるようになってきた今，日本測地系のメッシュデータを世界測地系の地図で使い続けるために，技術的なメモを書きためていきます。
なお，使っているシステムはWindows XPです。RやQGISはMacでも動くと思いますが，国土地理院が提供しているTKY2JGDと基盤地図情報閲覧コンバートソフトはWindows以外のOSで対応できるのかどうか調べていません。

2011-06-11

世界測地系の地図上の日本測地系3次メッシュポリゴンを作る

GIS

作業の流れ

緯度経度による定義通りのメッシュポリゴンをつくる（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として保存しておきます。

2011-06-11

世界測地系の地図上の日本測地系3次メッシュポリゴンを作る(2)

GIS

日本測地系から世界測地系への座標変換

国土地理院のサイトから座標変換プログラムTKY2JGDと座標変換用パラメータファイルをダウンロードして解凍します。最近のパソコンなら処理時間はずっと短縮されているのでパラメータファイルは全国版で問題ありません。
http://vldb.gsi.go.jp/sokuchi/tky2jgd/download/agreement.html
TKY2JGDに一括変換させるために，m3_dsのワークシートのD，E列をテキストエディタに貼り付けて保存します。この際，貼り付けたままだとtab区切りになっているのでtabをspaceへ置換します。また，1行目はタイトルなのではじめに#をつけてコメントアウトします。
ファイル名はmesh3.inとかがいいでしょう。

#Lat Lon
341000 1300000
341030 1300000
341030 1300045
341000 1300045
...

TKY2JGDを起動し，一括変換のボタンを押します。緯度・経度→緯度・経度を選択し，入力ファイルに今保存したmesh3.inを指定します。出力ファイルはmesh3.outでいいでしょう。
変換を実行すると「1万行に1時間ぐらいかかります」と表示されますが，私のノートパソコンでも4万行に10分ぐらいでしたから，そんなにはかからないと思います。
出力ファイルは最初に12行ほどコメントがついて，変換されたデータが続きます。

#Lat Lon
 341011.47758   1295951.64632   3parameters
 341041.47411   1295951.64550   3parameters
 341041.47496   1300036.64139   3parameters
 341011.47843   1300036.64221   3parameters
...

この3parametersというのは海上など地域ごとの変換パラメータがなかった点です。
#Lat Lonの行から最後までをコピーし，エクセルファイルに貼り付けます。うまくセル区切りをしてくれなかったら，データ→区切り位置でスペースを区切り文字にして実行してみてください。
先ほどのExcelファイルのm3_dsワークシートのF列に貼り付けた場合は以下のようになります。（1行目をつけなおしています）

ID	X	Y	Lat	Lon	#Lat	#Lon	Comment
1	130.000	34.167	341000	1300000	341011.478	1295951.646	3parameters
…

シェープファイルへの変換

この緯度経度の値は，(d)ddmmss.sssssという書式ですから，これをGISに持って行くには度（ddd.dddddddd）表示にする必要があります。
そのため，J2のセルに以下の式を入れます。

=INT(G2/10000)+(INT(G2/100)-INT(G2/10000)*100)/60+(G2-INT(G2/100)*100)/3600

J2セルをJ列とK列のセルにコピーします。またこの2列のセルの書式を「数値」にして小数点以下の桁数を8にします。すると以下のようになります。GISとはXYの順が逆になることに注意してください。

最後に，このm3_dsワークシートのA列とK列とJ列をコピーし（順番注意）新しいファイルのABC列に数値と書式をペーストしてテキスト形式で保存します。ファイル名はm3N_ds.txtでいいでしょう。

ID	Xw	Yw
1	129.99767953 	34.16985488 
1	129.99767931 	34.17818725 
1	130.01017816 	34.17818749 
1	130.01017839 	34.16985512 
...

次は，統計解析ツールである「Ｒ」を使ってこのファイルをGISで使うシェープファイルに変換します。Rをインストールし，さらにshapefilesというパッケージをダウンロードしてインストールしておいてください。
作業フォルダはC:\GIS\とします。先ほど作ったm3_dd.txt, m3W_ds.txt, m3N_ds.txtをこのフォルダに入れておいてください。
Rを立ち上げ，「新しいスクリプト」を開き，以下のコードを入力します。

library(shapefiles)
dd <- read.delim("C:/GIS/m3_dd.txt")
ds <- read.delim("C:/GIS/m3W_ds.txt")
ddShapefile <- convert.to.shapefile(ds, dd, "ID", 5)
write.shapefile(ddShapefile, "C:/GIS/f3meshW", arcgis=T)
ds <- read.delim("C:/GIS/m3N_ds.txt")
ddShapefile <- convert.to.shapefile(ds, dd, "ID", 5)
write.shapefile(ddShapefile, "C:/GIS/f3meshN", arcgis=T)

コードを選択してCtrl+Rを押すと実行します。数分でshapefileができあがると思います。

QGISで表示させるところは次回に
（続く）