コンテンツへスキップ

[地図自作] 国土地理院の数値標高モデルを使って等高線入り地図を作成する方法

公開日 [地図自作] 国土地理院の数値標高モデルを使って等高線入り地図を作成する方法 にコメントを残す

国土地理院が公開している数値標高モデルを使って等高線入り地図を自作してみましたので、その手順について公開します。

内容は試行錯誤を繰り返して作ったため洗練された手順とは言い難いのとかなりの時間とマシンパワーを要しますが、インターネットで公開されている等高線入り地図より格段に高解像度の地形図を作ることが出来ました。GPSにダウンロードした自作地図をいれて山登りされている方は挑戦の価値ありです。

 

はじめに

この記事では基盤地図情報の数値標高モデル(DEM10B)を使用して等高線入り地図を作る手順について紹介しています。この手順をなぞると、日本全体の等高線を作成するためにかなりのディスク容量を使用しますので作業HDDの空き容量に注意してください。

作業の大まかな流れとしては以下のようになります。

  1. 作業環境の構築
  2. 数値標高モデルのダウンロード
  3. 数値標高モデルをツールが扱える形式（GeoTIFF）に変換
  4. 等高線を作成（二通りの方法があります）
  5. Garmin GPS向けの地図形式に変換
garmin_oregon400t_scrn_mycustom_dem10b_step10m_mtfuji_200m.png
完成した等高線入り地図をGPSで表示するとこんな感じになります（注：この地図は国土地理院の基盤地図情報（デジタル標高モデル）DEM10Bを元に作成したものです）。

 

環境

PC環境について

今回はVMWarePlayerを使い仮想PCにUbuntuをインストールして環境構築しました。

  • ホストPC:
    • OS:Windows10(64bit)
    • CPU:Core i5
    • Memory:16GByte
    • HDD:2TByte
  • ゲスト:
    • OS:Ubuntu 16.04
    • CPU:4 core
    • Memory:12GByte
    • SWAP:16GByte（割り当てメモリが少ない場合はスワップファイルで補うことになります。もちろん速度は処理低下します）
    • HDD:1TByte（ここまで必要ありませんが、400GByteくらい欲しいです）
vmware_gmap_env_20170903.png
今回の作業環境はUbuntuを使って構築しました。最初はMSYS2を使ってWindowsで環境構築しようと試みたのですが、Python-GDALモジュールが上手くインストール出来ずに断念しました。

 

ソフトウェアについて

以下のアプリケーションを使用しました。また、地図作成に使用するmkgmapの設定やスタイルは小道具展示室さんが公開されている設定ファイルを使用させて頂きました。

  • python 2.7
  • Java
  • gdal 2.2.1
  • fgddemImporter
  • ogr2osm（巨大サイズのファイルに対応するため一部修正したもの）
  • phyghtmap 1.80
  • mkgmap-r3992
  • splitter-r584
  • 小道具展示室さんの設定ファイル("config170128.zip"をベースにスクリプトを変更したもの)

 

作業ディレクトリについて

本記事では以下の作業ディレクトリを使用しています。

  • /mnt/data/map/

作業環境の構築で個別にダウンロードしたツール類や設定ファイルは以下のディレクトリに配置します。

  • /mnt/data/map/tool

また、ダウンロードした基盤地図情報のデジタル標高モデルは以下のディレクトリに格納します。

  • /mnt/data/map/DEM10B

 

作業環境の構築

必要なパッケージのインストール

apt-getで必要なパッケージ類をインストールします。
以下のようにコマンドを入力してください。

$ sudo apt-get install git osmosis osmium-tool gdal-bin python python-pip python-numpy python-gdal python-matplotlib python-beautifulsoup python-lxml

インストールが完了したら、念のためgdalのバージョンを確認しておきます。

$ gdalinfo --version
GDAL 2.2.1, released 2017/06/23

もしgdalのバージョンが2.2より小さい場合は、以下のようにして最新版のgdalに更新してください。

$ sudo add-apt-repository -y ppa:ubuntugis/ubuntugis-unstable
$ sudo apt update
$ sudo apt upgrade

 

fgddemImporterのインストール

これは数値標高モデルをJPGIS(GML)形式からGeoTIFF形式に変換するためのQGISプラグインです。本来はQGISと組み合わせて使うものですが、今回は同梱されているfgddem.pyをJPGIS(GML->GeoTIFFコンバーターとして使用します。

適当な作業用ディレクトリを作成し、以下のようにしたgitを使ってダウンロードしてください。

$ cd /mnt/data/map/tool
$ git clone https://github.com/minorua/fgddemImporter.git

 

ogr2osmのインストール

ベクター形式のフォーマットからOSM(XML)形式に変換するためのツールです。大容量ファイルに対応するためにオリジナルのogr2osmに若干の変更を加えました。
以下のようにしてダウンロードして適当な場所に配置してください。

$ cd /mnt/data/map/tool
$ git clone https://github.com/suke-blog/ogr2osm.git
$ cd ogr2osm
$ git submodule update --init

 

phyghtmapのインストール

数値標高モデルから等高線を生成するツールです。このツールはSRTM-1/SRTM-3とGeoTIFFに対応していて、等高線の生成にSRTMを使用する場合はコマンドラインからオプション付きで入力してやるだけで必要なSRTMのダウンロードから等高線生成(osm/osm.pbf)まで自動的に行ってくれるものです。処理も高速で非常に便利ですが、今回の手順では等高線が上手く繋がらない箇所が発生します。この問題については後述します。

以下のサイトから最新の"Source distibution"をダウンロードしてインストールします。

以下のようにコマンドを入力してください。

$ wget http://katze.tfiu.de/projects/phyghtmap/phyghtmap_1.80.orig.tar.gz
$ tar zxvf phyghtmap_1.80.orig.tar.gz
$ cd phyghtmap-1.80
$ sudo python setup.py install

 

mkgmap/splitterのインストール

GPSが読み込める地図形式を生成するツールです。
以下のサイトからダウンロードしてください。

インストーラーは無いので、適当な場所に解凍して配置します。バージョンについては最新のものをダウンロードしてください。

$ cd /mnt/data/map/tool
$ wget http://www.mkgmap.org.uk/download/mkgmap-r3992.tar.gz
$ tar zxvr mkgmap-r3992.tar.gz
$ ln -s ./mkgmap-r3992 mkgmap

$ wget http://www.mkgmap.org.uk/download/splitter-r584.tar.gz
$ tar zxvf splitter-r584.tar.gz
$ ln -s ./splitter-r584 splitter

 

設定ファイルのインストール

小道具展示室さんが公開されているmkgmapの設定ファイルとスクリプト一式です。今回はオリジナルからLinux環境向けに少しスクリプトを変更したものを使用します。
以下のようにコマンド入力してください。

$ cd /mnt/data/map/tool
$ git clone https://github.com/suke-blog/mkgmap-config.git

なお、設定ファイルのオリジナルは以下のサイトから入手可能です。

 

数値標高モデルをダウンロードする

環境が整ったら、国土交通省の基盤地図情報ダウンロードサービスから数値標高モデルを入手します。ダウンロードにはユーザ登録が必要です。また、ここで入手したデータは一個人が私的に利用する分には問題ありませんが、第三者に公開（等高線入り地図をWebで配布する等）などを行う際に許可が必要となる場合があります。詳しくは基盤地図情報のサイトを参照ください。

公開されている数値標高モデルは4種類あります(2017/09/10現在)。このうち本記事で利用するは10mメッシュの"DEM10B"で1/25000地形図の作成にも使われています。DEM5A/DEM5Bは5mメッシュという高い解像度かつ標高誤差の少ないデータですが、現時点では日本全土をカバーしておらず歯抜けの地域があります。将来的に全土カバーされればより精度の高い等高線入り地図が作成できるかもしれません。

  • DEM5A(5mメッシュ、都市域/河川流域等)
  • DEM5B(5mメッシュ、都市域周辺等)
  • DEM10A(10mメッシュ、日本26火山のみ)
  • DEM10B(10mメッシュ、日本全国) ＜－今回の作業で利用

ちなみにこのダウンロードサイト、残念ながら日本全体をまとめてダウンロードが出来ません。ある程度選択してはダウンロードを繰り返すことになります(私は都道府県ごとにやりましたが、メッシュ単位の方が無難かも)。ファイル数が多くて大変ですが気長に作業してください。一応、全てダウンロードした時点でファイル数と容量は以下のようになりました。

 
$ ls -l FG-GML-*-DEM10B.zip | wc -l
4883

$ du -bhc FG-GML-*-DEM10B.zip
43K FG-GML-3036-50-DEM10B.zip
299K FG-GML-3622-57-DEM10B.zip
...
1.2M FG-GML-6848-26-DEM10B.zip
112K FG-GML-6848-27-DEM10B.zip
6.2G total

また、念のため"DEM10B"以外が紛れ込んでいないか確認してください。というのも、当方の環境ではダウンロード確認画面にいつの間にか選択した覚えのないファイルが増えていることがあったからです。気付かずにダウンロードして後でアレっ？となります（この現象は、"ダウンロードファイル確認へ"を押してファイルの確認＆ダウンロードを行った後、"戻る"前にファイル一覧をすべて"削除"すると発生しなくなりました）。

$ ls -l *DEM5*.zip *DEM10A.zip
QGIS
QGIS等のGeoTIFFを読み込めるツールを使ってダウンロード忘れが無いか確認しておくと安心出来ます（写真の白く抜けているところがファイルの存在しないところになります）

 

数値標高モデルをGeoTIFFに変換する

等高線の作成にはphyghtmap若しくはgdal_contourを使用しますが、これらのツールは基盤地図情報のJPGIS(GML)形式ファイルを読み込めません。そこでfgddem.pyを使ってGeoTIFF形式に変換してやります。流れとしては以下のようになります。

  • JPGIS(GML) --[1]--> GeoTIFF(3次メッシュ単位) --[2]--> VRT(GeoTIFFを一纏めにして管理)

[1] fgddem.pyでGeoTIFFに変換
[2] gdalbuildvrtでGeoTIFFを一纏めにする

 

JPGIS(GML)をGeoTIFFに変換

fgddem.pyを使ってJPGIS(GML)をGeoTIFF形式に変換します。JPGIS(GML)は圧縮ファイル(FG-GML-*-DEM10B.zip)のままで問題ありません。
以下のようにコマンドを入力してください。

$ python fgddem.py FG-GML-*-DEM10B.zip

コマンド引数には前項でダウンロードした数値標高モデルのファイルを指定しています。完了後にディレクトリを確認して"FG-GML-xxxx-nn-DEM10B.tif"という名前でファイルが出来ていることを確認してください。

$ du -bhc FG-GML-*-DEM10B.tif
3.3M	FG-GML-3036-50-DEM10B.tif
3.3M	FG-GML-3622-57-DEM10B.tif
3.3M	FG-GML-3623-06-DEM10B.tif
...
3.3M	FG-GML-6848-25-DEM10B.tif
3.3M	FG-GML-6848-26-DEM10B.tif
3.3M	FG-GML-6848-27-DEM10B.tif
16G	total

 

GeoTIFFを一纏めにする

前項で変換したGeoTIFFは3次メッシュ単位に分かれています。この後の作業の為にこれらのファイルを一纏めにして扱いやすくしておきます。
以下のようにコマンドを入力してください。

$ gdalbuildvrt japan-dem10b.vrt FG-GML-*-DEM10B.tif

各GeoTIFFファイルをタイルとしてインデックスしているだけなのですぐに完了すると思います。"japan-dem10b.vrt"というファイルが生成されていることを確認してください。ここまでが前準備です。

 

等高線を生成する

いよいよ等高線を作成するにあたって2通りの手段があります。

  • (ケース1).phyghtmapを使う場合
  • (ケース2).gdal_contourを使う場合

どちらでも等高線は出来ますが、それぞれ一長一短があります。

phyghtmapを使う場合、後者と比較して手間が少なく短時間で完了することがメリットです。難点はタイル境界にある等高線が途切れたり上手く繋がらない箇所が出ることです。出来上がった地図はパッと見では上手くできているように見えるのですが、拡大すると問題の箇所が見えてきます（そして一度問題の箇所を発見すると気になってしまう）。

gdal_contourを使う場合はその反対です。等高線は問題なく繋がりますが(私が確認した範囲では)、手順が多く処理時間ははるかに長くなります（私の環境では一晩では完了しません)。加えて中間生成するファイルがかなりディスク容量を消費します。

どちらを選ぶかですが、とりあえず作ってみようという方は"(ケース1)phyghtmap"の方がお手軽でお勧めです。山登り等でバリバリ使うので不明瞭な等高線はマズイという方は"(ケース2)gdal_contour"を使った手順に挑戦してみてください。

BaseCamp_contour_gap_ng_20170825.png
phyghtmapでGeoTIFFから等高線を生成する場合、上記のように所々うまく繋がらない箇所が出るようです（標高モデルにSRTMを使う場合は出ないor気にならないのですが・・・）

 

ケース1.phyghtmapを使う場合

phyghtmapを使って地図の等高線を生成します。ここからの作業の流れとしては以下のようになります。

  • VRT(GeoTIFFを一纏めにしたもの) --[1]--> GeoTIFF(タイルに分割) --[2]--> 等高線(OSM形式)

[1] gdal_retile.pyでGeoTIFFを1pxオーバーラップさせて分割
[2] phyghtmapでOSMもしくはPBF形式に変換

 

日本全体をタイルに分割する

phyghtmapは日本全体を含むような巨大サイズのファイルを読み込めません。そのため日本全体を含む"japan-dem10b.vrt"から適当な大きさのタイルに分割します。そのまま分割するとタイル境界で等高線が途切れてしまうためオプション"-overlap 1"を指定して1pxオーバーラップするように指定します。
以下のようにコマンドを入力してください。

$ mkdir retile
$ gdal_retile.py -ps 9000 6000 -overlap 1 -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "NUM_THREADS=ALL_CPUS" -targetDir ./retile japan-dem10b.vrt

コマンドオプションの詳細については以下のサイトを参考にしてください。

完了後、retileディレクトリに"japan-dem10b_xx_yy.tif"というファイルが大量に生成されていることを確認してください。

$ du -bhc japan-dem10b*.tif
1.9M japan-dem10b_01_01.tif
1.9M japan-dem10b_01_02.tif
...
1.4M japan-dem10b_38_31.tif
546K japan-dem10b_38_32.tif
6.8G total
BaseCamp 自作地図の等高線に途切れがある
ちなみに、オーバーラップさせないとタイル間の等高線がこのように途切れます

 

phyghtmapで等高線を生成

前項で分割したGeoTIFFファイルを使って等高線を生成します。
以下のようにコマンドを入力してください。

$ phyghtmap --pbf --output-prefix=japan_dem10b --step=10 --line-cat=500,50 --no-zero-contour --void-range-max=-500 --start-node-id=20000000000 --start-way-id=10000000000 --max-nodes-per-tile=1000000 --write-timestamp --jobs=4 japan-dem10b*.tif

オプションの"--jobs=n"はCPUコア数に応じて変更してください。ただし、今回のパターンですと1スレッドあたりピークで2GByte近くメモリを食うようなのでメモリが少ない環境では加減したほうが良いと思います。また、等高線の間隔等を変更する場合は以下の値を変更してください。

  • --step=10 : この場合は10m間隔で等高線が作られます
  • --line-cat=500,50 : 500m,50m間隔で等高線が太くなります

他のオプションの詳細については以下のサイトを参考にしてください。

完了後、同じディレクトリに拡張子osm.pbfのファイルが出力されていることを確認してください。あとはこれをGarmin.imgに変換すれば完了です。
（phyghtmapを使うと本当に手早く等高線が生成出来ます。あとはタイル境界で途切れることさえなければ言うことナシなんですが。。。）

 

ケース2.gdal_contourを使う場合

gdal_contourを使って等高線を生成する手順です。phyghtmapと違ってかなりの長丁場となりますので、何日かに分けて作業するつもりで挑んでください。
ここからの作業の流れとしては以下のようになります。

  • VRT(GeoTIFFを一纏めにしたもの) --[1]--> GeoTIFF(北海道/本州/etc地続き単位) --[2]--> 等高線(CSV等) --[3]--> 等高線(OSM形式) --[4]--> 等高線(OSMをタイル毎に分割)

[1] gdal_translate/gdalwrapを使用して北海道/本州/四国/九州等の陸続きになっている部分を切り出す
[2] gdal_contourを使用して等高線を生成
[3] ogr2osm/osmconvertを使用してOSM形式に変換
[4] splitterを使用してOSM形式ファイルを適当な大きさに分割

 

陸続きの部分を切り出す

まず、日本全体から北海道や本州のように地続きになっている部分をひとつのGeoTIFFファイルとして切り出します。これは等高線が意図せず途切れることを防ぐために地続きの範囲を一度に処理するためです（ファイルを分けるとそこで途切れる可能性がある）。処理の流れとしては以下のようになります。

  1. gdal_translateで範囲指定して矩形に切り出し(大雑把に北海道等を切り出す)
  2. gdalwarpで切り出したファイルの不要部分をマスク（例えば北海道を切り出すと座標の関係で青森辺りも入ってきてしまうのでKMLファイルを使ってマスクする)

今回は以下の範囲で切り出すことにします。

  • 北海道(japan-hokkaido.tif)
  • 本州(japan-honshu.tif)
  • 四国(japan-sikoku.tif)
  • 九州(japan-kyushu.tif)
  • 沖縄(japan-okinawa.tif)
  • 八丈島近辺(japan-hachijojima.tif)
  • 小川原近辺(japan-ogasawara.tif)

範囲はGoogle mapを使って作成しました。範囲指定は以下のURLからダウンロード出来ます。

ダウンロードは画面左上にあるメニュー("..."を縦にしたアイコン)を開いて、

  1. "KMLにエクスポート" を選択
  2. ダウンロードする範囲を選択("地図全体" -> "本州"等)
  3. "KMLでダウンロード"をチェック
  4. ダウンロードボタンを押下

とするとKMLファイルとしてダウンロード出来ます。ダウンロードしたファイルは適当な場所(ここでは./kml)に保存してください。保存する際にファイル名はローマ字にリネームしておいた方が扱いやすいと思います。

次に以下のようにコマンドを入力してください。

$ gdal_translate -of vrt -projwin 139.25 45.583333333 149 41.3333333 -a_nodata -9999 ./japan-dem10b.vrt /vsistdout/ | gdalwarp --config GDAL_CACHEMAX 500 -wm 500 -multi -cutline ./../../kml/Hokkaido.kml -wo "OPTIMIZE_SIZE=YES" -co "BIGTIFF=YES" -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "NUM_THREADS=ALL_CPUS" /vsistdin/ ./japan-hokkaido.tif
$ gdal_translate -of vrt -projwin 130.625 42.916666667 142.125 33.416666667 -a_nodata -9999 ./japan-dem10b.vrt /vsistdout/ | gdalwarp --config GDAL_CACHEMAX 500 -wm 500 -multi -cutline ./../../kml/Honshu.kml -wo "OPTIMIZE_SIZE=YES" -co "BIGTIFF=YES" -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "NUM_THREADS=ALL_CPUS" /vsistdin/ ./japan-honshu.tif
$ gdal_translate -of vrt -projwin 132 34.583333333 134.875 32.666666667 -a_nodata -9999 ./japan-dem10b.vrt /vsistdout/ | gdalwarp --config GDAL_CACHEMAX 500 -wm 500 -multi -cutline ./../../kml/Sikoku.kml -wo "OPTIMIZE_SIZE=YES" -co "BIGTIFF=YES" -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "NUM_THREADS=ALL_CPUS" /vsistdin/ ./japan-sikoku.tif
$ gdal_translate -of vrt -projwin 128 34.75 132.125 30.16666666 -a_nodata -9999 ./japan-dem10b.vrt /vsistdout/ | gdalwarp --config GDAL_CACHEMAX 500 -wm 500 -multi -cutline ./../../kml/Kyushu.kml -wo "OPTIMIZE_SIZE=YES" -co "BIGTIFF=YES" -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "NUM_THREADS=ALL_CPUS" /vsistdin/ ./japan-kyushu.tif
$ gdal_translate -of vrt -projwin 122.875 30.083333333 131.375 24.0 -a_nodata -9999 ./japan-dem10b.vrt /vsistdout/ | gdalwarp --config GDAL_CACHEMAX 500 -wm 500 -multi -cutline ./../../kml/Okinawa.kml -wo "OPTIMIZE_SIZE=YES" -co "BIGTIFF=YES" -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "NUM_THREADS=ALL_CPUS" /vsistdin/ ./japan-okinawa.tif
$ gdal_translate -of vrt -projwin 139.625 33.166666667 140.125 31.416666667 -a_nodata -9999 ./japan-dem10b.vrt /vsistdout/ | gdalwarp --config GDAL_CACHEMAX 500 -wm 500 -multi -cutline ./../../kml/Hachijojima.kml -wo "OPTIMIZE_SIZE=YES" -co "BIGTIFF=YES" -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "NUM_THREADS=ALL_CPUS" /vsistdin/ ./japan-hachijojima.tif
$ gdal_translate -of vrt -projwin 140.75 27.75 142.375 24.25 -a_nodata -9999 ./japan-dem10b.vrt /vsistdout/ | gdalwarp --config GDAL_CACHEMAX 500 -wm 500 -multi -cutline ./../../kml/Ogasawara.kml -wo "OPTIMIZE_SIZE=YES" -co "BIGTIFF=YES" -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "NUM_THREADS=ALL_CPUS" /vsistdin/ ./japan-ogasawara.tif

なお、gdalwarpの"--crop_to_cutline"オプションを使えばKMLで指定した範囲を自動的に切り出せるので便利ですが、切り出したファイルの座標にズレが生じるという問題があるため使用していません（切り出す座標はデジタル標高データの座標と一致させる必要があります。KMLファイルには私がマウスでポチポチしたアバウトな座標が入っているため、これを元に切り出すと再計算が行われてズレが生じてしまいます）。

各コマンドのオプション詳細については以下のサイトを参考にしてください。

処理が完了すると以下のようなファイルが出来上がります。

$ du -bhc japan*.tif
4.9M	japan-hachijojima.tif
1.2G	japan-hokkaido.tif
3.0G	japan-honshu.tif
558M	japan-kyushu.tif
15M	japan-ogasawara.tif
111M	japan-okinawa.tif
261M	japan-sikoku.tif
5.0G	total

 

gdal_contourで等高線を生成する

次にgdal_contourを使って前項で切り出した地続きのGeoTIFFから等高線を生成します。

GDALは大容量ファイルを扱うためにファイルをメモリマッピングする仕組みを持っているお陰かかなり大きなGeoTIFFも読み込むことが出来ますが、流石に本州くらいのサイズになるとメモリ消費量もかなりのものです（処理時間も数時間かかります）。十分にメモリ/スワップ領域を割り当ててから作業開始してください。参考までに私の環境では以下の割り当てで正常終了しました。

  • メモリサイズ:12GByte
  • スワップファイルサイズ:18GByte

これでもメモリから溢れてスワップしますがなんとか最後までいけました。その時のメモリ消費量から推測するにメモリ/スワップあわせて16GByteくらいは必要なのではないかと思います。

gdal_contourで等高線生成
gdal_contourで処理中のメモリ消費量。本州くらい大きいともりもりメモリを消費します。メモリ:8GByte/スワップ:2GByteだと進捗6割目あたりでメモリ不足になって落ちます。

また、gdal_contourはOSM形式ファイルを出力出来ないためひとまず別な形式で出力してから次行程で再変換する形になります。出力するファイル形式は以下URLに記載の"Creation=Yes"なものであればおそらく大丈夫です（ただし、Shapefile形式はファイルサイズに4GByte制限があるため途中でエラーが発生します）。ここではCSV形式で出力することにしました。

以下のようにコマンドを入力してください。
前項で生成したものに加えて離れ小島となっている2箇所が追加してあります。島の標高が低すぎて等高線には影響無いと思いますが念のため実行してください。

# generate contour
$ mkdir contour
$ gdal_contour -a ele -i 10.0 -f CSV -dsco "GEOMETRY=AS_WKT" japan-hokkaido.tif ./contour/japan-hokkaido.csv
$ gdal_contour -a ele -i 10.0 -f CSV -dsco "GEOMETRY=AS_WKT" japan-honshu.tif ./contour/japan-honshu.csv
$ gdal_contour -a ele -i 10.0 -f CSV -dsco "GEOMETRY=AS_WKT" japan-sikoku.tif ./contour/japan-sikoku.csv
$ gdal_contour -a ele -i 10.0 -f CSV -dsco "GEOMETRY=AS_WKT" japan-kyushu.tif ./contour/japan-kyushu.csv
$ gdal_contour -a ele -i 10.0 -f CSV -dsco "GEOMETRY=AS_WKT" japan-okinawa.tif ./contour/japan-okinawa.csv
$ gdal_contour -a ele -i 10.0 -f CSV -dsco "GEOMETRY=AS_WKT" japan-hachijojima.tif ./contour/japan-hachijojima.csv
$ gdal_contour -a ele -i 10.0 -f CSV -dsco "GEOMETRY=AS_WKT" japan-ogasawara.tif ./contour/japan-ogasawara.csv

$ gdal_contour -a ele -i 10.0 -f CSV -dsco "GEOMETRY=AS_WKT" FG-GML-3036-50-DEM10B.tif ./contour/FG-GML-3036-50-DEM10B.csv
$ gdal_contour -a ele -i 10.0 -f CSV -dsco "GEOMETRY=AS_WKT" FG-GML-3653-37-DEM10B.tif ./contour/FG-GML-3653-37-DEM10B.csv

オプションの詳細については以下のサイトを参考にしてください。

処理が完了すると以下のようなファイルが出来上がります。ここいらの工程からファイルサイズと処理時間（計っていなかったので正確なところは分かりませんが、本州だと6時間以上かかりました）が膨れ上がります。

$ du -bhc *.csv
11	FG-GML-3036-50-DEM10B.csv
11	FG-GML-3653-37-DEM10B.csv
14M	japan-hachijojima.csv
11G	japan-hokkaido.csv
38G	japan-honshu.csv
6.2G	japan-kyushu.csv
19M	japan-ogasawara.csv
423M	japan-okinawa.csv
3.9G	japan-sikoku.csv
59G	total

 

等高線をOSM形式に変換する

次にogr2osmを使って等高線をOSM(XML)形式ファイルに変換します。今回は数十GByteもある巨大ファイルに対応するためにオリジナルの処理に若干の修正を加えました。変換途中で逐次ファイルに出力するよう処理を追加してあります（これをしないとMemoryErrorで落ちます)。ここで変換したOSM(XML)は、この後でさらにOSM(PBF)形式ファイルに変換します。

以下のようにコマンドを入力してください。
最初のechoは等高線を構成する要素（Node, Way, Relation）に割り振る各IDの開始番号をファイルに書き出しています。ogr2osmは初期IDをこのファイルから読み出します。全てのIDが連番になるようにすることも可能ですが、もし何処かに問題があった時に一部だけやり直すことが出来るように各ファイルのIDを離してあります。

# convert to osm
$ cd ./contour
$ mkdir osm
$ cd ./osm

$ echo -e "10000000000\n20000000000\n30000000000" > id_hokkaido.txt
$ python /mnt/data/map/ogr2osm/ogr2osm.py --no-memory-copy --separate-id --idfile=./id_hokkaido.txt --saveid=./id_hokkaido.txt --positive-id --add-version --add-timestamp --sequential-output --translation=contour.py ../japan-hokkaido.csv

$ echo -e "11000000000\n21000000000\n31000000000" > id_honshu.txt
$ python /mnt/data/map/ogr2osm/ogr2osm.py --no-memory-copy --separate-id --idfile=./id_honshu.txt --saveid=./id_honshu.txt --positive-id --add-version --add-timestamp --sequential-output --translation=contour.py ../japan-honshu.csv

$ echo -e "14000000000\n24000000000\n34000000000" > id_sikoku.txt
$ python /mnt/data/map/ogr2osm/ogr2osm.py --no-memory-copy --separate-id --idfile=./id_sikoku.txt --saveid=./id_sikoku.txt --positive-id --add-version --add-timestamp --sequential-output --translation=contour.py ../japan-sikoku.csv

$ echo -e "15000000000\n25000000000\n35000000000" > id_kyushu.txt
$ python /mnt/data/map/ogr2osm/ogr2osm.py --no-memory-copy --separate-id --idfile=./id_kyushu.txt --saveid=./id_kyushu.txt --positive-id --add-version --add-timestamp --sequential-output --translation=contour.py ../japan-kyushu.csv

$ echo -e "16000000000\n26000000000\n36000000000" > id_okinawa.txt
$ python /mnt/data/map/ogr2osm/ogr2osm.py --no-memory-copy --separate-id --idfile=./id_okinawa.txt --saveid=./id_okinawa.txt --positive-id --add-version --add-timestamp --sequential-output --translation=contour.py ../japan-okinawa.csv

$ echo -e "17000000000\n27000000000\n37000000000" > id_hachijojima.txt
$ python /mnt/data/map/ogr2osm/ogr2osm.py --no-memory-copy --separate-id --idfile=./id_hachijojima.txt --saveid=./id_hachijojima.txt --positive-id --add-version --add-timestamp --sequential-output --translation=contour.py ../japan-hachijojima.csv

$ echo -e "18000000000\n28000000000\n38000000000" > id_ogasawara.txt
$ python /mnt/data/map/ogr2osm/ogr2osm.py --no-memory-copy --separate-id --idfile=./id_ogasawara.txt --saveid=./id_ogasawara.txt --positive-id --add-version --add-timestamp --sequential-output --translation=contour.py ../japan-ogasawara.csv

$ echo -e "19000000000\n29000000000\n39000000000" > id_others.txt
$ python /mnt/data/map/ogr2osm/ogr2osm.py --no-memory-copy --separate-id --idfile=./id_others.txt --saveid=./id_others.txt --positive-id --add-version --add-timestamp --sequential-output --translation=contour.py ../FG-GML-3036-50-DEM10B.csv
$ python /mnt/data/map/ogr2osm/ogr2osm.py --no-memory-copy --separate-id --idfile=./id_others.txt --saveid=./id_others.txt --positive-id --add-version --add-timestamp --sequential-output --translation=contour.py ../FG-GML-3653-37-DEM10B.csv

処理が完了すると以下のようにOSM形式ファイルが出来上がります。OSM(XML)形式ファイルの中身はXMLということもあって100GByteオーバーの巨大ファイルになりました。

$ du -bhc *.osm
86 FG-GML-3036-50-DEM10B.osm
86 FG-GML-3653-37-DEM10B.osm
56M japan-hachijojima.osm
43G japan-hokkaido.osm
147G japan-honshu.osm
25G japan-kyushu.osm
75M japan-ogasawara.osm
1.7G japan-okinawa.osm
16G japan-sikoku.osm
232G total

今回のようにノード数が多いとファイルサイズが非常に大きくなるため、これをOSM(PBF)形式に変換してサイズを圧縮します。この形式はOSM(XML)と比較して10%以下にまでファイルサイズを圧縮してくれます。OSM(XML)のままで処理を進めることも出来ますが、HDD容量を圧迫するうえに後工程の処理速度も遅くなるので変換したほうが無難だと思います。HDDの空き容量が厳しい方は各ファイル毎に"等高線->OSM(XML)->OSM(PBF)"と変換してOSM(XML)を削除したほうがいいかもしれません。

以下のようにコマンドを入力してください。

# osm to pbf
osmconvert --out-pbf ./japan-hokkaido.osm -o=japan-hokkaido.osm.pbf
osmconvert --out-pbf ./japan-honshu.osm -o=japan-honshu.osm.pbf
osmconvert --out-pbf ./japan-sikoku.osm -o=japan-sikoku.osm.pbf
osmconvert --out-pbf ./japan-kyushu.osm -o=japan-kyushu.osm.pbf
osmconvert --out-pbf ./japan-okinawa.osm -o=japan-okinawa.osm.pbf
osmconvert --out-pbf ./japan-hachijojima.osm -o=japan-hachijojima.osm.pbf
osmconvert --out-pbf ./japan-ogasawara.osm -o=japan-ogasawara.osm.pbf
osmconvert --out-pbf ./FG-GML-3036-50-DEM10B.osm -o=FG-GML-3036-50-DEM10B.osm.pbf
osmconvert --out-pbf ./FG-GML-3653-37-DEM10B.osm -o=FG-GML-3653-37-DEM10B.osm.pbf

処理が完了すると以下のようにOSM(PBF)形式ファイルが出来上がります。232GByte -> 4.1GByteと劇的にファイルサイズが縮小しているのが分かるかと思います。変換が終われば元のOSM(XML)は削除して問題ありません。

$ du -bhc *osm.pbf
132 FG-GML-3036-50-DEM10B.osm.pbf
132 FG-GML-3653-37-DEM10B.osm.pbf
1016K japan-hachijojima.osm.pbf
748M japan-hokkaido.osm.pbf
2.7G japan-honshu.osm.pbf
448M japan-kyushu.osm.pbf
1.4M japan-ogasawara.osm.pbf
31M japan-okinawa.osm.pbf
281M japan-sikoku.osm.pbf
4.1G total

 

ファイルを分割する

お使いのGarminGPSの機種によって差異はあると思いますが、地図中の1タイルあたり含まれるノード数には制限があります。前項で変換したOSMファイル(1ファイル=1タイルと考えていいと思います)は大きすぎてノード数の制限に引っかかるため、splitterを使って適当なノード数になるよう分割する必要があります。

以下のようにコマンドを入力してください。
ここでは各タイルのノード上限を10000000としています。また、splitterのパスは各環境にあわせ置き換えてください。

$ mkdir split

# hokkaido 
$ java -Xmx4096M -jar /mnt/data/map/tool/splitter/splitter.jar --keep-complete=true --description=JAPAN-DEM10B --max-nodes=1000000 --search-limit=10000000 --resolution=12 --mapid=63000001 --output-dir=./split  ./japan-hokkaido.osm.pbf

# honshu
$ java -Xmx16384M -jar /mnt/data/map/tool/splitter/splitter.jar --keep-complete=true --description=JAPAN-DEM10B --max-nodes=1000000 --search-limit=10000000 --resolution=12 --mapid=63100001 --output-dir=./split ./japan-honshu.osm.pbf

# sikoku
$ java -Xmx4096M -jar /mnt/data/map/tool/splitter/splitter.jar --keep-complete=true --description=JAPAN-DEM10B --max-nodes=1000000 --search-limit=10000000 --resolution=12 --mapid=63400001 --output-dir=./split ./japan-sikoku.osm.pbf

# kyushu
$ java -Xmx4096M -jar /mnt/data/map/tool/splitter/splitter.jar --keep-complete=true --description=JAPAN-DEM10B --max-nodes=1000000 --search-limit=10000000 --resolution=12 --mapid=63500001 --output-dir=./split ./japan-kyushu.osm.pbf

# okinawa
$ java -Xmx4096M -jar /mnt/data/map/tool/splitter/splitter.jar --keep-complete=true --description=JAPAN-DEM10B --max-nodes=1000000 --search-limit=10000000 --resolution=12 --mapid=63600001 --output-dir=./split ./japan-okinawa.osm.pbf

# hachijojima
$ java -Xmx4096M -jar /mnt/data/map/tool/splitter/splitter.jar --keep-complete=true --description=JAPAN-DEM10B --max-nodes=1000000 --search-limit=10000000 --resolution=12 --mapid=63700001 --output-dir=./split ./japan-hachijojima.osm.pbf

# ogasawara
$ java -Xmx4096M -jar /mnt/data/map/tool/splitter/splitter.jar --keep-complete=true --description=JAPAN-DEM10B --max-nodes=1000000 --search-limit=10000000 --resolution=12 --mapid=63800001 --output-dir=./split ./japan-ogasawara.osm.pbf

# others
$ java -Xmx2048M -jar /mnt/data/map/tool/splitter/splitter.jar --keep-complete=true --description=JAPAN-DEM10B --max-nodes=1000000 --search-limit=10000000 --resolution=12 --mapid=63900001 --output-dir=./split ./FG-GML-3036-50-DEM10B.osm.pbf
$ java -Xmx2048M -jar /mnt/data/map/tool/splitter/splitter.jar --keep-complete=true --description=JAPAN-DEM10B --max-nodes=1000000 --search-limit=10000000 --resolution=12 --mapid=63910001 --output-dir=./split ./FG-GML-3653-37-DEM10B.osm.pbf

コマンドを実行していると、分割元OSMファイルのサイズによってはメモリ不足により処理途中でエラーが発生するかもしれません。エラーを回避するにはメモリ割り当てを増やすか、もしくは別ツールで分割する必要があります。特に本州は巨大なため、私の環境だとオプションで"-Xmx16384M"(16GByte!)としてようやく正常終了となりました。

オプションの詳細については以下のサイトを参考にしてください。

完了すればsplitディレクトリに63000001から始まる連番のファイルが生成されているはずなので確認してください。

$ du -bhc *.osm.pbf
1.9M 63000001.osm.pbf
2.7M 63000002.osm.pbf
1.9M 63000003.osm.pbf
...
3.8M 63600016.osm.pbf
1.2M 63700001.osm.pbf
1.7M 63800001.osm.pbf
7.6G total

 

地図形式に変換

ようやくここまで来ました。既にかなり時間（と忍耐とマシンパワー）を使われたと思います。
最後の仕上げとして、mkgmapを使って前項までに"phyghtmap"もしくは"gdal_contour"で作成した等高線をGPSで読み込める地図形式に変換します。等高線以外の情報はOpenStreetMapのデータを使用しました。

この辺りの手順については以下の記事が参考になるかと思います。

参考記事ではSRTMを使って等高線を生成していますが、そこで出来るosm.pbfファイルを今回作った"63000001.osm.pbf"等の連番ファイルと置き換えてください。また、地図生成の際に実行するスクリプト"gmapsupp_make.sh"の先頭に定義してあるTARGET_TYPEを"2"(等高線あり地図の生成)にしてスクリプトを実行してください。

$ cd /mnt/data/map/tool/mkgmap-config
$ ./gmapsupp_make.sh

mkgmapの処理だけ見ると、私の環境では2時間弱かかりました。
完了後にディレクトリを確認すると"gmapsupp.img"が出力されていると思います。後はBaseCampに登録するなり、GPSに"gmapsupp.img"を直接転送するなりしてご活用ください。

BaseCamp_dem10b_mtfuji.png
やっと完成した地図をBaseCampで確認中、ここまで長かったです（比喩抜きで一夏でした）。感無量であります。

 

所感

ここまで長々とお付き合い頂いてありがとうございました。

さて、今回は手順を並べるだけとはいえ書くのに結構苦労がありました。ネット検索してはみたものの、どうにも基盤地図情報を使って等高線を生成する具体的な手順は見つかりません。基本的な道筋はネットで見つかるものの一から十まで手順書があるわけではないので色々と調べる必要があります。ツールがエラーを吐くたびにソースコードを追いかけたり資料を探したりしながら解決策を探し、地図を合成しては確認を繰り返しでした。そして、どうにか地図に等高線がのり、おっしゃ出来たぁー！とプチ祝杯を挙げた後に不自然にちょん切れた等高線を発見した時の落胆具合は何とも言えないものがあります。ただ、最終的には力技（マシンパワー）で強引に解決となったものの個人的に納得のいく出来栄えの地図を作ることが出来たと思います。

本記事の手順は試行錯誤の結果でスマートさとは程遠くありますが、もし興味のある方がおられましたら基盤地図情報を使ってオリジナル地図作りに挑戦してみてください。そしてもっと素晴らしい手順を開拓された暁にはこそっと教えて頂けると嬉しいです。

また、出来るだけ読んでくれる方が再現しやすいように手順を書いたつもりですが、何分あーだこーだとやりながら書いているので分かりにくい箇所もあると思います。気になる箇所を見つけられたらさり気無くコメント頂けますと幸いです。

 

参考

本記事は以下の内容を参考に記載させていただきました。

More from my site

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください。