(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
感想
ほんとは整飾部分がないものが購入できるといいんですがー。