世界測地系の地図上の日本測地系3次メッシュポリゴンを作る(2)
日本測地系から世界測地系への座標変換
国土地理院のサイトから座標変換プログラム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 ...
次は,統計解析ツールである「R」を使ってこのファイルを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で表示させるところは次回に
(続く)