中部の等高線付き地図を作成

夏休みにアルプスに登山する計画を立てていて、まだ作っていない中部地方の等高線付き地図を作成しました。
前回やった近畿地方の地図のやり方をそのまますればよいと軽く考えていましたが、備忘録として書いたつもりの記事が間違いだらけで役立たずだった(記事を読んだ方に申し訳ない)ので、結構時間がかかってしまいました。
alasixosaka.hatenablog.com
alasixosaka.hatenablog.com
備忘録を兼ねてちゃんとやる方法を書いておきます。この先、関東や中国地方などの地図を作る時に困らないように。

全体の流れ

まず、作業全体の流れです。

  1. オープンストリートマップの日本地図から必要な部分を切り出す。
  2. 本州の等高線地図から必要な部分を切り出す。
  3. 世界の海岸線地図から必要な部分を切り出す。
  4. 海の領域を指定した地図を作成し、3と合体する。
  5. できた地図に1と2を合体する。
  6. mapsforgeで読めるようにmap形式に変換する。

下準備

作業開始前に確認しておくこと、用意するデータなどです。

JAVA関連の設定

まず、つまずいたのはJAVA関連の設定でした。前回作業したときはJREを使って作業を行ったのですが、現在はJREJDKに統合されてしまっており、自分のパソコンでもJREは削除されていて、JDKがインストールされていました。にもかかわらず、環境変数を修正していなかったので、コマンドを動かそうとするとJAVAがないとエラーメッセージが表示されてしまいました。
エラーを回避するため、環境変数に現在インストールされているJAVAのフォルダのパスを通しておきます。
JAVAの格納フォルダは私の場合は、c\ProgramFiles\java でした。環境変数の設定を開いて、PATHの編集をクリックし、JREのパスが記載されていた部分をJDKのパスに変更します。この際、再起動を行わないといけないので注意が必要です。
また、メモリの割り当ても、JDKをインストールするときに変更されたのか、何かの都合で変更したのか覚えていませんが、512Mbyteになっていて、作業途中でまたもout of memoryのエラーが出現したので、6Gbyteに変更しました。メモリ割り当ての変更も環境変数の設定から行います。

元データの用意

今回の作業では、必要なプログラムはインストール済み、また大元のデータとなる、本州の地図、等高線、世界の海岸線のデータはすでにあるという前提で作業を進めます。

実際の作業

切り出す範囲の決定

地図の切り出しは、osmosisを使って行いました。切り出しには、経度、緯度それぞれに上限、下限の計4点の座標が必要です。
googleマップで適当な地点をクリックし、経度、緯度を確認し、次の4点を切り出す座標に決定しました。
経度 下限 135.4493 上限 139.1322
緯度 下限 34.2763 上限 37.5402

オープンストリートマップの切り出し。

ダウンロードしてある(あるいは新規にダウンロードして)、japan-latest.osm.pbfをosmosisを使って、上記の座標範囲を切り出します。コマンドは下記の通り。

Osmosis --read-pbf file="japan-latest.osm.pbf" --bounding-box top=37.5402 left=135.4493 bottom=34.2763 right=139.1322 --write-pbf file="chubu.osm.pbf"

ちなみに、pbfはオープンストリートマップを圧縮した形式のファイルの拡張子です。(pdfではありません)
圧縮されていないファイルはosmという拡張子を付けます。

等高線の切り出し。

以前の作業で作成した本州部分の等高線(honsyu_DEM10b.pbf)から、下記のコマンドで中部地方を切り出します。

Osmosis --read-pbf file="honshu_dem10b.osm.pbf" --bounding-box top=37.5402 left=135.4493 bottom=34.2763 right=139.1322 --write-pbf file="chubu_dem10b.osm.pbf"

海岸線データの切り出し

海岸線のデータは以前にダウンロードした全世界の海岸線データからogr2ogrを使って切り出します。拡張子はshpです。

ogr2ogr -overwrite -progress -skipfailures -clipsrc 135.4493 34.2763 139.1322 37.5402 chubu_o1.shp land-polygons-split-4326/land_polygons.shp

できたchubu.shpをshape2osmを使ってosm形式に変換します。

py -2 -m shape2osm -l chubu_o1_ns -o 10000000 chubu_o1.shp

以前のブログにも書いた通り、mapsforgeのオリジナルサイトからとってきたものと、参考サイトからとって来たものの2つがあります。
前にも書いた通り前者でないとうまくいきません。(どちらを使っても同じようです。下にある修正をしないとうまくいきませんでした)
ファイル名はchubu_ns1.osmとなります。

海の領域を指定した地図を作成する。

ここは単純に領域全体を海に指定します。メモ帳などのテキストエディタで下記のファイルを作成し、chubu_s.osmとして保存します。あるいは参考サイトからダウンロードしても良いと思います。
yueno.net

<?xml version='1.0' encoding='UTF-8'?>
<osm version="0.6" generator="osmconvert 0.8.2">
	<bounds minlat="34.2763034" minlon="135.449339" maxlat="38.8973083" maxlon="139.899887"/>
	<node timestamp="1969-12-31T23:59:59Z" changeset="20000" id="32951459320" version="1" lon="135.449339" lat="34.2763034" />
	<node timestamp="1969-12-31T23:59:59Z" changeset="20000" id="32951459321" version="1" lon="135.449339" lat="38.8973083" />
	<node timestamp="1969-12-31T23:59:59Z" changeset="20000" id="32951459322" version="1" lon="139.899887" lat="38.8973083" />
	<node timestamp="1969-12-31T23:59:59Z" changeset="20000" id="32951459323" version="1" lon="139.899887" lat="34.2763034" />
	<way timestamp="1969-12-31T23:59:59Z" changeset="20000" id="32951623372" version="1">
		<nd ref="32951459320" />
		<nd ref="32951459321" />
		<nd ref="32951459322" />
		<nd ref="32951459323" />
		<nd ref="32951459320" />
		<tag k="area" v="yes" />
		<tag k="layer" v="-5" />
		<tag k="natural" v="sea" />
	</way>
</osm>

海の領域と海岸線を合体する。

まず、chubu_s.osmとchubu_ns1.osmを合体して、海と海岸線の地図を作成します。コマンドは下記の通り。

osmosis --rx file=chubu_s.osm --rx file=chubu_o1_ns1.osm --s --m --wx file=chubu_o1_sea.osm 

できた地図のファイル名はchubu_o1_sea.osmとなります。

地図、等高線と合わせる。

先ほどのchubu_o1_sea.osmとchubu.osm.pbf、chubu_DEM10b.pbfをosmosisを使って合体します。コマンドは下記の通り。

osmosis --rx file=chubu_o1_sea.osm --rb file=chubu.osm.pbf --s --m --rb file=chubu_DEM10b.pbf --s --m --wx file=chubu_o1c.pbf  omitmetadata=true

できたファイルはchubu_o1c.pbfというファイル名になります。omitmetadata=trueとしているのは、途中でエラーが出るためで、何故なのかとか、このオプションをつけるとどうなるのか詳しいことは判りません。エラーメッセージにそうせよと書いてあるのでつけてあります。

map形式に変換します。

これらのファイルを合わたものを、map形式へ変換します。
map形式への変換はmapsforgeのmap writerプラグインを使います。これもインストールされているという前提で作業を進めます。コマンドは下記のとおりです。

osmosis --rb file=chubu_o1.pbf --mw file=chubu_o.map bbox=34.276300,135.449300,37.540200,139.132200 tag-conf-file=tags.xml map-start-zoom=10 type=hd

bboxのオプションは、すべての切り出し範囲をこの範囲に統一して、海のエリアもこの範囲にしておけば不要だと思います。参考サイトの海のエリアのファイルchubu_s.osmを使ったため、エラーメッセージが出て、bboxオプションでエリアを指定せよと言われたので、つけています。念のためいつもつけておいた方が良いかもしれません。オプションのtag-conf-file=tags.xmlについては後で書きます。これがないと等高線が表示されません。
これで、chubu_o1c.mapというファイルが出来上がります。これを地図アプリに読み込ませると等高線入りの地形図が表示されます。

うまくいかなかった事

基本的な作業の流れは、以前のブログに書いたように、上記の参考サイトの方法を参考にして行っています。

このサイトでは、各地方のデータを切り出すのにポリゴンファイルを用いて、osmconvert を使って作業しています。同様の方法を用いてもよかったのですが、以前に近畿地方の地図を作成した時も今回と同じく、経度緯度を使って4角に切り出して作業したので今回も同じ方法で作業を行いました。

伊勢地方が水没する問題。

以前のブログでも、海岸線と海のエリアをうまく指定してやらないと陸地が水没してしまうと書きましたが、今回作成した地図では、三重県の沿岸部の大半が水没してしまいました。下の図は、マップ形式への変換時間を短縮するために等高線を入れずに変換したものを表示させています。

伊勢が水没している。

同じミスをしたのかと何度かやり直してみましたが、うまくいかないので、よくよく調べてみると、海岸線データを切り出したファイルをJOSMで読み込んでみてみたところ、伊勢地方の海岸線が一部途切れていることがわかりました。ちなみに、shape2osmは2つあると書きましたが、オリジナルサイトからとってきたものを使うと、”changeset=-1"となっていて、JOSMではエラーになって読み込めません。参考サイトのものは”changeset=1"となっていてうまく読み込めます。マップファイルへの書き出しはどちらを使っても問題なかったです。
海岸線データを切り出した元のファイル、chubu_o1.shpを同じくJOSMで読み込むとこちらは海岸線がちゃんと表示されます。どちらのshape2osmを使っても同じことが起こります。何故なのか最初はよくわからなかったのですが、答えがアプリ制作の時に大変お世話になった、伊勢在住のプログラマーさんのブログに書いてありました。
blog.mori-soft.com
このサイトによると、あまりちゃんと理解してないのかもしれませんが、OSMOSISで海岸線を切った時に、もともとつながっていた海岸線が切れてしまっているのに、一つの海岸線のようにデータ上なってしまうことが原因のようです。multipolygonになると書かれています。対処法も書かれていて、下のサイトのpythonスクリプトを使えば、multipolygonをpolygonに修正してくれるとのこと。
grazie Paolo · GitHub
ここで、インプットファイル名は、24行目に書かれている"poly.shp"です(といってもはてなブログでは行数が表示されないのでわかりにくいですが。はてなブログソースコードを記述するには不向きな気がします)。同じくアウトプットのファイル名は26行目に書かれている"polys.shp"です。これらを、自分のファイル名に書き換えて実行してやれば、修正されたシェイプファイルができます。それをshape2osm(どちらでもOK)を使ってosm形式に変換したものを使えば水没が解消しました。

import os, sys
from osgeo import ogr,gdal

def multipoly2poly(in_lyr, out_lyr):
    for in_feat in in_lyr:
        geom = in_feat.GetGeometryRef()
        if geom.GetGeometryName() == 'MULTIPOLYGON':
            for geom_part in geom:
                addPolygon(geom_part.ExportToWkb(), out_lyr)
        else:
            addPolygon(geom.ExportToWkb(), out_lyr)

def addPolygon(simplePolygon, out_lyr):
    featureDefn = out_lyr.GetLayerDefn()
    polygon = ogr.CreateGeometryFromWkb(simplePolygon)
    out_feat = ogr.Feature(featureDefn)
    out_feat.SetGeometry(polygon)
    out_lyr.CreateFeature(out_feat)
    print 'Polygon added.'


gdal.UseExceptions()
driver = ogr.GetDriverByName('ESRI Shapefile')
in_ds = driver.Open('poly.shp', 0)
in_lyr = in_ds.GetLayer()
outputshp = 'polys.shp'
if os.path.exists(outputshp):
    driver.DeleteDataSource(outputshp)

out_ds = driver.CreateDataSource(outputshp)
out_lyr = out_ds.CreateLayer('polys', geom_type=ogr.wkbPolygon)
multipoly2poly(in_lyr, out_lyr) 
正常に表示された状態(こちらは等高線ありの地図)


それから躓いたのは、上にも書きましたが、mapfile-writerでtag-conf-file=tags.xmlというオプションをつけるところでした。上にも書きましたが、これがないと等高線が表示されません。以前の自分のブログにオプションをつけないコマンドを書いてしまっていたので、間違いに気づくのに随分時間がかかりました。osmosisを最新版にしたり、mapfile-writerも最新版にしても解消せず、Windows版がダメなのかと思って、Linuxも使っていましたが同じことだったので随分悩みました。参考サイトをちゃんと読んでいれば、そのことが書いてあったのですが、一応読んでいたのですが斜め読みで肝心のところを読み飛ばしてしまっていて、気づくのに時間がかかってしまいました。ちなみに、tags.xmlというファイルは参考サイトからダウンロードさせていただいたものを使っています。
と言う訳で無事に等高線付きの中部地方の地図ができました。あまりやる作業ではないですが、結構手順が多いので、pythonスクリプトで自動化をしたいと思っています。出来たら、またこのブログに載せておきます。