プログラマーのメモ書き

伊勢在住のプログラマーが気になることを気ままにメモったブログです

OpenStreetMap のオフラインマップ作成時の不具合の修正 (3/3)

OpenStreetMap mapsforge 避難所検索＠伊勢

前回までの記事で、オフラインマップで洪水状態になる原因がわかりました。ここでは、それを踏まえて、どうやって修正するかを検討した際のメモを残しておきます。

なお、全部で３つに分かれています

この記事の目次は次の通りです。

対応策

さて、前回までの挙動から対応策を考えてみました。

  • land-polygons-split の伊勢近辺を囲んでいる矩形サイズに切り出す
  • multipolygon を複数の polygon に変換するツールをかませる

この２つぐらいが思いつきます。

前者のほうは、ogr2ogr のオプション指定の際に -clipsrc でland-polygons をぶった切らない範囲を指定するのでたぶんうまくいくと思います。とはいえ、land-polygons-split から kansai 部分だけを切り出したデータを見ると、

f:id:junichim:20200924153455p:plain

こんな感じに、伊勢志摩だけでなく三重県の大部分を含むような範囲になるので、最終的に作成する map ファイルのファイルサイズがかなり大きくなりそうです。 それに、何かの折に land-polygons の分割の仕方が変わるとそれに合わせて追随しないといけないのが、ちょっとイヤです。

なので、後者の方法で、できないか検討してみました。

multipolygon を複数の polygon へ変換（QGIS）

これを調べてみると、まずは QGIS を使えば簡単にできることがわかりました。メニューより、

f:id:junichim:20200924153607p:plain

Vector -> Geometry Tools -> Multipart to Singleparts を選択すると

f:id:junichim:20200924153820p:plain

とダイアログがで表示されるので、対象となる layer と出力先となるシェープファイル名を指定すると、書き出してくれます。

例えば、もともとは

f:id:junichim:20200924154142p:plain

こんな感じで、黄色の部分が１つのPOLYGON（MULTIPOLYGON）だったのですが、この機能を使って分割すると

f:id:junichim:20200924154220p:plain

こんな感じで、一つ一つに分かれます。

multipolygon を複数の polygon へ変換（コマンドラインツール）

いろいろとテストしている分には QGIS でもよいのですが、これをコマンドラインツールでやって、既存の map-creator から呼び出したいと思います。

きっと、 GDAL/OGR 周りで何かあるはずだと思いながら、いろいろと調べたのですがなかなか見つかりません。でも、ネットはすごいですね。それなりに時間とられましたが、最終的に、結構昔のディスカッションで、似た現象のお悩み相談がありました。

GDAL - Dev - Multipart to singlepart

で、頑張って最後のほうまで見ていくと、multipolygon を polyon に変換する python ファイルを gist にアップしてくれてました！

いやー、ありがたいですね。感謝感謝としかいいようがありません。

試しに、ogr2ogr で伊勢近辺に限定して（-clipsrc を指定して）作成したシェープファイルを、このスクリプトに通して、別のシェープファイルに変換して、ogrinfoでみてみると、multipolygon の表記が見事消えていました。

これで動くんじゃないかな？ということで、map-creator を一部修正して、ogr2ogr 呼び出しまでをコメントアウトして、この multipolygon をなくしたファイルを渡す形で shape2osm から動かして mapファイルを作ります。

で、アプリに組み込んで動かすと

f:id:junichim:20200924155108p:plain

問題なしです！

あとは、これを定期更新のスクリプトに組み込めば一件落着です。

ogrinfo と shpdump の出力について（2020/9/28 追記）

参考までに、multipolygon 状態のpolygonを含んだファイルを ogrinfo と shpdump で出力したらどうなるか載せておきます。

ogrinfo で表示してみます。シェープファイル名は land.shp レイヤー名が land の場合です。

osm@map:~/tmp$ ogrinfo land.shp land

こんな感じでファイル全体の情報が出た後に、個々のPOLYGONの情報が出力されます。

INFO: Open of `land.shp'
      using driver `ESRI Shapefile' successful.

Layer name: land
Geometry: Polygon
Feature Count: 1542
Extent: (136.417000, 34.277000) - (137.028000, 34.685000)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]
FID: Real (11.0)
OGRFeature(land):0
  FID (Real) = 172823
  POLYGON ((136.945087 34.556601,136.945052 34.556614,136.945062 34.556628,136.945095 34.556621,136.945087 34.556601))

OGRFeature(land):1
  FID (Real) = 173390
（略）

Geometry: Polygon の表記からこのファイルに含まれている種類は POLYGON ということがわかります。

でも、MULTIPOLYGONのものについては、

OGRFeature(land):1463
  FID (Real) = 177513
  MULTIPOLYGON (((136.802623666666676 34.277,136.802637 34.277032,136.8026681 34.2770763,136.8026781 34.2770963,136.802686200000011 34.2771317,136.802732
（略）

のように MULTIPOLYGON と表示してくれます。

一方、 shpdump の場合は

Shapefile Type: Polygon   # of Shapes: 1542

File Bounds: (136.417,34.277,0,0)
         to  (137.028,34.685,0,0)

Shape:0 (Polygon)  nVertices=5, nParts=1
  Bounds:(136.945052,34.556601, 0)
      to (136.945095,34.556628, 0)
     (136.945087,34.556601, 0) Ring 
     (136.945052,34.556614, 0)  
     (136.945062,34.556628, 0)  
     (136.945095,34.556621, 0)  
     (136.945087,34.556601, 0)  

Shape:1 (Polygon)  nVertices=5, nParts=1
（略）

こんな感じになり、上記と同じ shape を表示させると

Shape:1463 (Polygon)  nVertices=57, nParts=2
  Bounds:(136.802623666667,34.277, 0)
      to (136.804933382124,34.2773938, 0)
     (136.802623666667,34.277, 0) Ring 
     (136.802637,34.277032, 0)  
（略）

となっていて、MULTIPOLYGON ということは一目でわかりません。 nParts = 2 が MULTIPOLYGON の可能性をにおわせていますが、ホールの場合もあり得るので、ちょっとわかりにくいかもしれませんね。

定期更新スクリプトへの組み込み（2020/9/28 追記）

オフラインマップファイルの定期更新についてこちらの記事でまとめています。ここで触れているように、対象領域を絞り込むため、若干 map-creator スクリプトを修正しています。この修正した map-creator の152行目あたりを手直しします。

# Land

ogr2ogr -overwrite -progress -skipfailures -clipsrc $LEFT $BOTTOM $RIGHT $TOP "$WORK_PATH/land.shp" "$DATA_PATH/land-polygons-split-4326/land_polygons.shp"
python shape2osm.py -l "$WORK_PATH/land" "$WORK_PATH/land.shp"

の部分を

# Land

ogr2ogr -overwrite -progress -skipfailures -clipsrc $LEFT $BOTTOM $RIGHT $TOP "$WORK_PATH/land.shp" "$DATA_PATH/land-polygons-split-4326/land_polygons.shp"
# added by Junichi MORI, 2020/9/24
CURRENT=`pwd`
pushd $WORK_PATH
python $CURRENT/multipolygon2polygons.py
popd
# end of added
# modified by Junichi MORI, 2020/9/24
#python shape2osm.py -l "$WORK_PATH/land" "$WORK_PATH/land.shp"
python shape2osm.py -l "$WORK_PATH/land" "$WORK_PATH/polys.shp"

としておきます。

この multipolygon2polygons.py が上で触れた gist にあったスクリプトです（入出力ファイル名だけ変えています）。このスクリプトを正しく呼び出せるようにカレントディレクトリを移動して、呼び出し処理をしているだけです。

全体を通して動かして、きちんとmapファイルが作成されることを確認しました。 いやー、長かったですね。

（参考） 下記の記事にも、別のやり方がありますので、ご参考までに。

qgis - Converting huge multipolygon to polygons - Geographic Information Systems Stack Exchange