（2014/3/17 windows7と最新のGDALで動かないバグがあったので修正しました。）
電子地形図25000（定形図郭版）をgdalを使ってgeotiffにしてくっつける方法を紹介します。オンデマンド版を使えば好きな図郭で購入できますが、広い範囲を利用したい場合などにどうぞ。ただ、境界部分の黒枠線は若干残ります。
1.購入
http://dkgd.gsi.go.jp/dkgx/page1.htm
ここから購入。定形図郭版のTIFFを選択します。まちがえないように。
2.データ整理
ダウンロードした画像を1次メッシュフォルダごとにまとめて、KANRIファイルとスクリプトファイルを1次メッシュフォルダと同じ階層におきます。ファイル名を以下のように「二次メッシュ番号.tif」「二次メッシュ番号.tfw」に変更します。
|-電子地形図変換.vbs(このあと保存するスクリプト) |-KANRI2KN.CSV |-5539-553935.tif | -553935.tfw | -553925.tif | -553925.tfw |-5742-574242.tif | -574242.tfw | -574343.tif | -574243.tfw
KANRIファイルは同じ種類で中身をまとめます。テキストエディタでCSVを開いて、一つのファイルにコピペします。
ﾌｧｲﾙ名,本図・分図ﾌﾗｸﾞ,対応する地形図の図名,対応する地形図の図名（読み）,地形図ｺｰﾄﾞ,号数,基本測量,改測,修正,部分修正,空中写真,現地調査,発行年月日,図郭四隅経緯度左上緯度,分,秒,左上経度,分,秒,左下緯度,分,秒,左下経度,分,秒,右下緯度,分,秒,右下経度,分,秒,右上緯度,分,秒,右上経度,分,秒,図郭四隅UTM(M単位)左上E,左上N,左下E,左下N,右下E,右下N,右上E,右上N,図郭四隅ﾌｧｲﾙ座標左上P,左上L,左下P,左下L,右下P,右下L,右上P,右上L,2次ﾒｯｼｭｺｰﾄﾞ,延伸フラグ左上,上,右上,右,右下,下,左下,左,調製日,磁気偏角西偏（度）,分 553935,1,五十里湖,いかりこ,,,,,,,,,,37,0,0,139,37,30,36,55,0,139,37,30,36,55,0,139,45,0,37,0,0,139,45,0,377655,4095756,377521,4086511,388656,4086358,388777,4095603,490,802,486,8198,9394,8198,9389,802,553935,0,0,0,0,0,0,0,0,20140109,7,20 553925,1,川治,かわじ,,,,,,,,,,36,55,0,139,37,30,36,50,0,139,37,30,36,50,0,139,45,0,36,55,0,139,45,0,377521,4086511,377388,4077266,388535,4077113,388656,4086358,486,802,481,8198,9399,8198,9394,802,553925,0,0,0,0,0,0,0,0,20140109,7,20
3.スクリプトで変換
以下のスクリプトを「電子地形図変換.vbs」という名前で保存して、GDALをインストールしたOSGeo4Wのコマンドプロンプト上で実行します。
変換元の図郭のUTMゾーンに従ってs_epsgの部分を書き換えてください。
投影変換したいepsgコードに従って、t_epsgの部分を書き換えてください。
epsgコードについてはこちらを参考
http://d.hatena.ne.jp/tmizu23/20091215/1260868350
＞電子地形図変換.vbs 0 (KANRI2KN.CSVの変換) ＞電子地形図変換.vbs 2 (tmpC*.vrtをマージ＆投影変換してout.tifに変換)
'お好きな投影法で
'UTM53 3099
'UTM54 3100
'xy9 2451
'xy10 2452
'世界測地系BL 4612
s_epsg = "epsg:3100"
t_epsg = "epsg:2451"
Set objWshShell = WScript.CreateObject("WScript.Shell")
Set tempEnv = objWshShell.Environment("Process")
tempEnv.Item("GDAL_DATA") = "data"
tempEnv.Item("GDAL_FILENAME_IS_UTF8") = "NO"
Set args = WScript.Arguments
If args.count<1 Then
WScript.Echo "使い方 電子地形図変換.vbs 0 or 2 (0:KANRI2KN.CSV 2:マージ)"
WScript.Quit
End If
If args.item(0) = 0 Then
kanri = "KANRI2KN.CSV"
ElseIf args.item(0) = 2 Then
merge = 1
End If
If merge <> 1 Then
Set objFso = CreateObject("Scripting.FileSystemObject")
Set objFile = objFso.OpenTextFile(kanri, 1, False)
If Err.Number > 0 Then
WScript.Echo "Open Error"
Else
objFile.ReadLine
Do Until objFile.AtEndOfStream
str = Split(objFile.ReadLine,",")
mesh = str(0)
If args.item(0) = 0 Then
gcp0 = str(45) & " " & str(46) & " " & str(37) & " " & str(38)
gcp1 = str(47) & " " & str(48) & " " & str(39) & " " & str(40)
gcp2 = str(49) & " " & str(50) & " " & str(41) & " " & str(42)
gcp3 = str(51) & " " & str(52) & " " & str(43) & " " & str(44)
ulx = str(16)*1.0 + str(17)/60.0 + str(18)/3600.0
uly = str(13)*1.0 + str(14)/60.0 + str(15)/3600.0
lrx = str(28)*1.0 + str(29)/60.0 + str(30)/3600.0
lry = str(25)*1.0 + str(26)/60.0 + str(27)/3600.0
End If
strCmdLine = "cmd /c gdalwarp -of VRT -r bilinear --config GDAL_CACHEMAX 2048 " & left(mesh,4) & "\" & mesh & ".tif" & " " & "tmpA" & mesh & ".vrt"
flag = objWshShell.Run(strCmdLine,,True)
strCmdLine = "cmd /c gdalwarp -of VRT -s_srs " & s_epsg & " -t_srs epsg:4612 -r bilinear --config GDAL_CACHEMAX 2048 " & "tmpA" & mesh & ".vrt" & " " & "tmpB" & mesh & ".vrt"
flag = objWshShell.Run(strCmdLine,,True)
strCmdLine = "cmd /c gdal_translate -of VRT --config GDAL_CACHEMAX 2048 -projwin " & ulx & " " & uly & " " & lrx & " " & lry & " " & "tmpB" & mesh & ".vrt" & " " & "tmpC" & mesh & ".vrt"
flag = objWshShell.Run(strCmdLine,,True)
Loop
End If
objFile.Close
Set objFile = Nothing
Set objFso = Nothing
Else
strCmdLine = "cmd /c gdalwarp -tr 2.15 2.15 -co PROFILE=GeoTIFF -co COMPRESS=LZW -multi --config GDAL_CACHEMAX 2048 -r bilinear -co TFW=YES -s_srs epsg:4612 -t_srs " & t_epsg & " tmpC*.vrt out.tif"
flag = objWshShell.Run(strCmdLine,,True)
End If
感想
ほんとは整飾部分がないものが購入できるといいんですがー。
