電子地形図25000(定形図郭版)をGDALで余白を消しつつマージする方法

（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

感想

ほんとは整飾部分がないものが購入できるといいんですがー。

参考

定形図郭版とオンデマンド版の比較
http://net.jmc.or.jp/digital_data_gsiol_denshiChizu25000_hikaku.html