数値地図25000(地図画像)をGDALで余白を消しつつマージする方法

数値地図25000（地図画像）をgdalを使ってgeotiffにしてくっつける方法を紹介します。自分用メモなので説明不足の感がありますが、とりあえず公開します。

1.購入

http://net.jmc.or.jp/digital_data_gsiol_ima25000.html
ここから必要な範囲をダウンロード購入

2.データ整理

ダウンロードした画像を1次メッシュフォルダごとにまとめて、KANRIファイルとスクリプトファイルを1次メッシュフォルダと同じ階層におきます。

|-数値地図2GeoTIFF.vbs(このあと保存するスクリプト)
|-KANRI2KN.CSV(あれば)
|-KANRI2KW.CSV(あれば)
|-5741-574142.tif
|     -574143.tif
|-5742-574242.tif
      -574343.tif

KANRIファイルは同じ種類で中身をまとめます。テキストエディタでCSVを開いて、一つのファイルにコピペします。

ﾌｧｲﾙ名,本図・分図ﾌﾗｸﾞ,図名(漢字),図名(よみ),地形図ｺｰﾄﾞ,号数,基本測量,改測,修正,部分修正,空中写真,現地調査,発行年月日,図郭四隅経緯度左上緯度,分,秒,左上経度,分,秒,左下緯度,分,秒,左下経度,分,秒,右下緯度,分,秒,右下経度,分,秒,右上緯度,分,秒,右上経度,分,秒,図郭四隅UTM(ﾋﾟｸｾﾙ単位)左上E,左上N,左下E,左下N,右下E,右下N,右上E,右上N,図郭四隅ﾌｧｲﾙ座標左上P,左上L,左下P,左下L,右下P,右下L,右上P,右上L,2次ﾒｯｼｭｺｰﾄﾞ,延伸フラグ左上,上,右上,右,右下,下,左下,左,図葉ﾌｧｲﾙ作成日,磁針方位西偏(度),分
574142,1,渡波,わたのは,NJ-54-15-10-4,石巻10号-4,S43,,H12,,H10-09,H12-12,H14.1.1,38,25,10.7,141,14,47.5,38,20,10.7,141,14,47.5,38,20,10.7,141,22,17.4,38,25,10.7,141,22,17.4,8608,1700962,8618,1697263,12987,1697278,12972,1700977,92,76,91,3773,4459,3773,4457,76,574142,0,0,0,0,0,0,0,0,H15.6.9,7,10
574143,1,荻浜,おぎのはま,NJ-54-15-10-2,石巻10号-2,S43,,H12,,H10-11,H12-12,H13.9.1,38,25,10.7,141,22,17.4,38,20,10.7,141,22,17.4,38,20,10.7,141,29,47.4,38,25,10.7,141,29,47.4,12972,1700977,12987,1697278,17357,1697299,17337,1700997,93,75,91,3773,4459,3773,4456,75,574143,0,0,0,0,0,0,0,0,H15.6.9,7,10

3.スクリプトで変換

以下のスクリプトを「数値地図2GeoTIFF.vbs」という名前で保存して、GDALをインストールしたOSGeo4Wのコマンドプロンプト上で実行します。マージ&投影変換時の投影法は、スクリプトのepsgの部分を好きなものに書き換えてください。

＞数値地図2GeoTIFF.vbs 0 (KANRI2KN.CSVの変換)
＞数値地図2GeoTIFF.vbs 1 (KANRI2KW.CSVの変換)
＞数値地図2GeoTIFF.vbs 2 (tmpC*.tifをマージ＆投影変換してout.tifに変換)
Set objWshShell = WScript.CreateObject("WScript.Shell")
Set args = WScript.Arguments

If args.count<1 Then
 WScript.Echo "使い方  数値地図2GeoTIFF 0〜2 (0:KANRI2KN.CSV 1:KANRI2KW.CSV 2:マージ)"
 WScript.Quit
End If

If args.item(0) = 0 Then
 kanri = "KANRI2KN.CSV"
ElseIf args.item(0) = 1 Then
 kanri = "KANRI2KW.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(16)*1.0 + str(17)/60.0 + str(18)/3600.0 & " " & str(13)*1.0 + str(14)/60.0 + str(15)/3600.0
		        gcp1 = str(47) & " " & str(48) & " " & str(22)*1.0 + str(23)/60.0 + str(24)/3600.0 & " " & str(19)*1.0 + str(20)/60.0 + str(21)/3600.0
		        gcp2 = str(49) & " " & str(50) & " " & str(28)*1.0 + str(29)/60.0 + str(30)/3600.0 & " " & str(25)*1.0 + str(26)/60.0 + str(27)/3600.0
		        gcp3 = str(51) & " " & str(52) & " " & str(34)*1.0 + str(35)/60.0 + str(36)/3600.0 & " " & str(31)*1.0 + str(32)/60.0 + str(33)/3600.0
			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
		Else
		
			gcp0 = str(49) & " " & str(50) & " " & str(20)*1.0 + str(21)/60.0 + str(22)/3600.0 & " " & str(17)*1.0 + str(18)/60.0 + str(19)/3600.0
		        gcp1 = str(51) & " " & str(52) & " " & str(26)*1.0 + str(27)/60.0 + str(28)/3600.0 & " " & str(23)*1.0 + str(24)/60.0 + str(25)/3600.0
		        gcp2 = str(53) & " " & str(54) & " " & str(32)*1.0 + str(33)/60.0 + str(34)/3600.0 & " " & str(29)*1.0 + str(30)/60.0 + str(31)/3600.0
		        gcp3 = str(55) & " " & str(56) & " " & str(38)*1.0 + str(39)/60.0 + str(40)/3600.0 & " " & str(35)*1.0 + str(36)/60.0 + str(37)/3600.0
			ulx = str(20)*1.0 + str(21)/60.0 + str(22)/3600.0
			uly = str(17)*1.0 + str(18)/60.0 + str(19)/3600.0
			lrx = str(32)*1.0 + str(33)/60.0 + str(34)/3600.0
			lry = str(29)*1.0 + str(30)/60.0 + str(31)/3600.0
		End If
		strCmdLine = "gdal_translate --config GDAL_CACHEMAX 2048 -expand rgb -gcp " & gcp0 & " -gcp " & gcp1 & " -gcp " & gcp2 & " -gcp " & gcp3 & " " & left(mesh,4) & "\" & mesh & ".tif" & " " & "tmpA" & mesh & ".tif"
		flag = objWshShell.Run(strCmdLine,,True)
		strCmdLine = "gdalwarp --config GDAL_CACHEMAX 2048 -r bilinear " & "tmpA" & mesh & ".tif" & " " & "tmpB" & mesh & ".tif"
		flag = objWshShell.Run(strCmdLine,,True)
		strCmdLine = "gdal_translate --config GDAL_CACHEMAX 2048 -projwin " & ulx & " " & uly & " " & lrx & " " & lry & " " & "tmpB" & mesh & ".tif" & " " & "tmpC" & mesh & ".tif"
		flag = objWshShell.Run(strCmdLine,,True)

	    Loop
	End If
	objFile.Close
	Set objFile = Nothing
	Set objFso = Nothing
Else
	'お好きな投影法で
	epsg = "epsg:3100"
	strCmdLine = "gdalwarp -wm 2048 --config GDAL_CACHEMAX 2048 -r bilinear -co TFW=YES -s_srs epsg:4612 -t_srs " & epsg & " tmpC*.tif out.tif"
	flag = objWshShell.Run(strCmdLine,,True)
End If

感想

  • 新旧図郭が混在したものを余白を消してマージするのに苦労した。0(余白部分)をnodataにしてマージすると、新旧図郭で重なる範囲で等高線が二重になってしまう（新旧で少しズレがあるから）。なので、余白部分をあらかじめクリップすることに。ただ、これも簡単ではない。数値地図画像は、UTM座標系の画像を緯度経度の範囲で切り取って、北を真上にするために少し回転させたもの。そのため、きれいにクリップしてマージさせるためには以下の手順が必要。
    1. KANRI.csvの情報からGCPを画像に付けて、それを元に緯度経度の座標系に変換する。
    2. GCPの範囲でクリップする。
    3. マージする。
    4. UTM座標系に変換する。
  • 変換の途中でリサンプリングアルゴリズムをbilinearにしている関係で、カラーパレットから3バンドになってしまいます。nearだとカラーパレットのままでいけるけど、それだとギザギザなので仕方なし。
  • 法人で地図画像を購入した場合の使用ライセンスはどうなってるのだろうか？
  • なにゆえ、いまだにこの形式で販売してるのか疑問。是非とも余白なしのGeoTIFFでお願いします。（この一文を書くための前フリでした。）

参考