2016年12月14日水曜日

GeoPackage で ポータブル OSM

FOSS4G Advent Calendar 2016 14 日目の記事です!

Advent Calendar でしか更新されないブログへようこそ。
NPO法人 オープンコンシェルジュ とか OSGeo.JP の松浦です。

みなさんは GIS データをどんな形式で保存していますか?
私のようなリッチメンは家計簿兼暖房用に購入した Oracle Exadata に Oracle Locator のSDO_Geometry と SDO_Georaster で管理していますが(真っ赤なウソ)、大抵の場合は Shapefile じゃないでしょうか。

Shapefile は米国Esri 社が開発した空間データを格納するフォーマットで、仕様が公開されていることもあり多くの GIS ソフトウェアで使用可能です。
ただ、ハードウェア、ソフトウェアの急速な進化に伴って、いろいろ不都合が生じてきました。

Shapefile の闇


1. なんかファイルいっぱい必要


「データ添付しました。ご査収ください。」っていうメールが来て確認したら、満を持して登場する【データ.shp】。「どうも、一人できました。てへっ」つって。
開けねーよ。ご査収できねぇよ。

なんて経験、ありますよね?
Shapefile は最低でも .shp, .dbf, .shx という 3 つのファイルが必要です。
必須以外にも、.prj、.sbn、.shp.xml、.cpg などなどいろいろ構成員がいらっさいます。これら合わせて 1 つの Shapefile です。

2. なんか、もう文字コードがつらい


最近はソフトウェアが頑張るのでましですが、たいてい CP932 の Shapefile は一発目化ける。

3. 容量制限


今時 2 GB って・・・
まぁ 2 GB のベクターデータって結構でかいんですがね。
とはいっても、TB を超えるディスクが 1, 2 万で買える時代に 2 GB はきつい。

4. 属性名足りない、短いよ...


10 バイトまで (UTF-8 だと 3 文字とか・・・)。なので「都道府県」という属性名にしたくても「都道府」になってしまうという・・・


とまぁきりがないのですが、決して Shapefile を陥れたいわけでもディスりたいわけでもないんです。
このフォーマットが GIS の普及やソフトウェアのデータ相互運用性の確保に大きく貢献をしたことは言うまでもありません。

そう、ただ次の時代がやってきただけです。

GeoPackage を使おう


じゃぁ、ポスト Shapefile はなんでしょう。
個人的には、PostgreSQL にぶっこんどけばよくない?とか思うけど。おらの捕まえたGISデータを交換とかする友達いねーし。

ファイル形式で、追加ソフトやライブラリがいらなくて、大容量に対応していて仕様がオープンなフォーマットって考えたとき、脳内には GeoPackage しか残りませんでした。
ほかにあるのかもしれないけど。

GeoPackage はベクターデータおよびタイル化されたラスターデータを格納できますが、その実態は SQLite データベースです。空間インデックス(R-Tree)にも対応しています。

単一ファイル (.gpkg) で構成されますが、データベースですので 1 つの GeoPackage 内に複数のスキーマの異なるベクターデータやタイル化されたラスターデータを格納できます。この点は Shapefile と大きく異なります。

また、OGC (Open Geospatial Consortium) によって定義されており、すでに GDAL/OGR や QGIS や ArcGIS などの地理空間データを扱うソフトウェアが対応しています。

詳しくはこちら(私が言っていることと違うことが書いてあったら、こちらのサイトのほうを信じてください)。

GeoPackage に Planet.osm をぶっこむ

はい。やっと本題。
結局でかいデータを高パフォーマンスで扱えるかなんてやってみなきゃわからんのだよ。
ってことで恒例の Planet.osm 登場。

Planet.osm は Open Street Map の地図描画用の元データで、全世界の道路や建物、POI、行政界、国境などが含まれます。
OSM より提供されており Protocol Buffers 形式で 35 GB 強あります。通常は全球データを落として 1 ファイルにするなんてクレイジーなことはせず、PostgreSQL などの RDBMS に展開して使用します。たぶん。
ちなみに PostgreSQL に展開すると、ぼんやりした記憶ですがインデックスこみこみで 200 GB 前後になったと思います。ポリゴンデータは 2 億レコードを超え、ラインでも軽く 1 億レコード越えです。

さて、冒頭で Exadata の話をしましたが、現実世界の私の家にはおそらく Exadata の筐体より小さな冷蔵庫ぐらいしかないので、データ作成はクラウドを活用することにします。
AWS? え? あれもえもえしないので今回は Conoha でいきます。
Conoha な理由は、root ディスク、データディスクすべてSSDで、時間で課金形式だけどデータ転送量の課金がなく、最大コストを見積もりやすいからです。庶民の味方。(GMOさん、宣伝したからなんかちょうだい)

今回はメモリ 4GB の Arch Linux のインスタンスに 500 GB の追加 SSD をつけてデータ作成を実施しました。

結論から言うと、ほぼ弊害なさくっと Planet.osm から、ポリゴン、ライン、ポイント を含んだ GeoPackage を作成することができました。

以下は簡単な手順と、それぞれにかかった時間です。

  1. GDAL/OGR の導入
    以下コマンドを実行

    # pacman -S gdal

    処理時間: 1分
  2. Planet.osm のダウンロード
    今回はお隣、台湾のサーバーからダウンロード。
    Planet の情報は http://wiki.openstreetmap.org/wiki/Planet.osm 参照。

    nohup /usr/bin/time curl -O http://free.nchc.org.tw/osm.planet/pbf/planet-161205.osm.pbf > download_planet.log 2>&1 &
    処理時間: 53 分
  3. ogr2ogr で変換

    GDAL/OGR の  ogr2ogr で変換。

    nohup /usr/bin/time ogr2ogr --config CPL_TMPDIR /data/osmtest/tmp -progress -f GPKG planet_osm.gpkg planet-161205.osm.pbf > ogr2ogr.log 2>&1 &

    処理時間: 22 時間 30 分
    GeoPackage サイズ: 170 GB

    Conoha 利用料: (なんだかんだ 3 日間くらいあーだこーだして) 1000円弱


    ポイントは --config CPL_TMPDIR オプションで、処理の一時ファイルの出力先に空き領域に余裕のあるディレクトリを指定すること。
    最初の数回はここでは待ってディスク不足で処理が止まりました。

    あと、メモリも 4GB ぐらいはあったほうがいいみたい。最初 512 MB のインスタンスに swap 積んでやったんですが out of memory でお亡くなりになりました。

    ちなみに上記処理は空間インデックスの作成も含みます。その他のカラムに対するインデックスは張られていません。
  4. gpkg のダウンロード

    これが予想外に時間がかかった・・・
    nginx をたてて、HTTPでダウンロードしましたが、4時間以上かかりましたよ・・・。そしてこの時間ですよ。
予想外に処理時間が短く済みましたし、容量も 170 GB という現実的な数字に収まりました(もはや麻痺)。

あとはこれが現実的な速度で使用できるかということです。

QGIS で描画してみる

QGIS の出番です。
意気揚々と「ベクタレイヤの追加」ボタンから作成した planet_osm.gpkg を追加したら・・・ 15分以上沈黙・・・ その後レイヤ選択のウィンドウが出てきてデータ追加できました。
ただ、毎回これじゃぁ・・・ ねぇ?

これ、レイヤのプロパティ用にすべてのレイヤに対して select count(*) 的なクエリを流してフィーチャ数を毎回カウントしているからだと思われます。

なんとかこのプロセスを端折る方法はないか10秒程度悩んだんですが、1個だけありました。

VRT (http://www.gdal.org/drv_vrt.html)

VRTは仮想的なレイヤを定義する XML 形式のデータフォーマットです。通常は複数のデータをあたかも1つのデータとして扱うためなどに使われます。
VRTの定義の中にメタ情報的な要素もあり、その中に FeatureCount という要素があったことを思い出しました。神がかってます。

こんな感じで定義しました。

これを poly.vrt として保存し、QGIS で開くと、さくっ!といきましたよ。同じ要領でポイント、ラインの VRT も作成しました。
フィーチャ数は適当に2億にしておきましたw ちゃんと数える場合は sql で直接カウントしてください。
ちなみにポリゴンのレイヤは 247647486 フィーチャあります。

描画速度も実用に耐えうるもので、ある程度拡大していればさくさくと行きます。

東京


New York

多すぎるのでポイントは非表示

結論

GeoPackage 使える。

想像してみてください。
全球の地図データをポケットに入れて持ち運べる時代がきたんですよ?

興奮しませんか?


よいお年を。

※深夜および過労気味で文章のテンションがおかしく、かつ文も乱れております。詳細をお聞きになりたい方は私の実態を捕まえるか、SNSで連絡ください。
twitter: @PEmugi2

注意点なぐりがき

  • データの細部まで検証してないので、中身がすべてまともにインポートできてるかは謎。
  • Conoha の料金は参考値
  • 空間インデックスしか張ってないので、たとえば国境ポリゴンを表示させようと admin_level in ('0', '1') とかやっても、インデックスは空間的な絞り込みにしか聞かないので小縮尺ではえらい時間かかります。フィルタとかする場合は適切にインデックスを追加しましょう。
  • ogr2ogr の osm に関する情報はこちら
  • デフォルトで属性として取り込まれる osm のタグは限られます。たとえばname:jaなどはその他のタグとしてまとめられます。この設定を変えるには上記のogr2ogrのページの osmconf.ini で定義します。

2015年12月11日金曜日

QGIS GEOHEX Plugin

NPO法人 オープンコンシェルジュ の松浦です。

この記事は FOSS4G Advent Calendar 2015 の 12/11 の記事でっす。

みなさん GEOHEX ってご存知ですか?
GEOHEX は @sa2da 氏が開発した世界地図を隙間なく敷き詰める六角形(ヘックス)です。

GeoHash やタイルスキーマの六角形版のようなものをイメージしてもらえればいいと思います。

こちらで @sa2da 氏作成のデモが見られます。

GEOHEX の詳しい仕様については公式ページをご参照ください。

で、本題に入りますが、GEOHEX は見た目もかっこいいし、便利に使えると思うので多くの人に使ってほしいのですが、簡単に GEOHEX を作成するようなツールがありません。
GEOHEX を扱うためのライブラリは多くの言語に対して実装されていますが、できることなら GUI でサクッと作りたいですよね。

というわけで QGIS の Plugin にしてみました。
2年前に Python 用のライブラリを書いていたのでそれを使用してやっつけで作りました。
(Python ライブラリもライセンスを定義しないまま放置していたので MIT ライセンスをこの機会に設定・・・)

使ってみよう

QGIS にこの Plugin を入れると、下のように緑の六角形のアイコンが追加されます。アイコンがダサいのは許して。
これをぽちっと押すとダイアログが出て、GEOHEX のレベル(サイズ)と作成範囲の元となるレイヤを設定します。

レベルは 0 から数字が大きくなるにつれて GEOHEX のサイズは小さくなります。キャプチャの表示範囲でレベル 9 とか 10 にすると QGIS 様が沈黙しやがりますので注意してください。

作成範囲 (Set extent from:) はレイヤを選択するとそのレイヤを内包する最小の四角形の範囲と重なる GEOHEX がすべて生成されます。リストの最後の "Current Map Window" を選択すると表示されている範囲を六角形が埋め尽くします。やってみよう。



これでレベル 5 です。次にちょうど地図の中央に見えるポイントレイヤー(室蘭市公開のAED設置場所を示すオープンデータ)を内包する範囲で作成します。レベル 9 で!



かなり拡大しましたが、それでも細かいw



Plugin の入手とインストール

現在この Plugin は QGIS の公式 Plugin レポジトリに登録中ですが、まだ許可?が出ていないので QGIS の Plugin ツールからインストールすることはできません。
そしてやっつけのコードで許可してもらえるのか怖いw

というわけで使ってみたいという方は以下のリンクからダウンロードしてください。
https://onedrive.live.com/redir?resid=F359D930199F0E64!1592&authkey=!AFTJ10mLaGzO4Ys&ithint=file%2czip

ダウンロードしたら zip ファイル内の「GEOHEX」フォルダを QGIS の Plugin ディレクトリにフォルダごとコピーしてください。
Windows の場合は「C:\Users\ユーザー名\.qgis2\python\plugins」、Unix 系 OS の場合は 「~/.qgis2/python/plugins」ですかね?。
コピー後に QGIS を再起動すればインストール完了です。

その他注意事項など

GEOHEX レイヤーの保存先

作成された GEOHEX レイヤーはメモリ上に保存されます。なので、広範囲で細かい HEX を発生せるとメモリを圧迫するので気を付けてください。
また、メモリ上に保存されていますので QGIS を終了すると消えます。保存したい場合は Shape File や sqlite あたりにエクスポートしてください。

GEOHEX レイヤーの空間参照

そもそも GEOHEX は Web Mercator 前提の仕様ですが、ミクロなスケールでは他の座標系でも使えると思っています。なので現状この Plugin では範囲指定に使用したレイヤーの空間参照を引き継ぎます。
"Current Map Window" を設定した場合には QGIS で現在表示されている地図の空間参照を引き継ぎます。

無謀なレベル設定はやめよう

えっと、たとえば地理院地図を表示して、Extent の対象に選んでしまうと全世界を埋める GEOHEX が作成されます。この時にレベル 0 とかならまだしも 10 とかやっちゃうと高確率でレスポンスが途絶えます。マシンスペックにもよりますが。
警告とか出さずに愚直に六角形を紡いでいく真面目な Plugin なのでユーザー様が気を付けていただければ幸いです。

まだ beta1 です

せっかく作ったので、拡張したりバグをつぶしたりしていく予定です。
直近の課題としては、GEOHEX の最新バージョンに対応していないので早急に対応したいです。ほとんど変わらないと思いますが、-180 度をまたぐ HEX でのコードの扱いが変わっているようでした。

まとめ

GEOHEX は本当にかっこいいし便利だとおもっているので、これからいろいろ実利用を検討していきたいなと思っています。

かっちょえー UI の地図アプリに表示された GEOHEX をぽちっとタップするとドローンがそこに飛んでいくとか!

オープンデータの可視化とかにも使えそうですよね。

以上今年の年末の一番大事な仕事を終えます。
次は yuskesuzki@github さんです!

2014年12月13日土曜日

PostgreSQL(PostGIS) でベクトルタイルを作る決意を固めたものの挫折しそう

この記事は FOSS4G Advent Calendar 2014 の 12/13 (土) の記事です。

さて、ここ1年色々環境が変化したわけで、あまり FOSS4G 関連の場に出ていけてないそんな私ですが、今年は久しぶりに書きたいと思います。
ブログもうやってなかったんですが、このために作りました。来年また使うと思います。

ベクトルタイルです。ベクトルタイルですってよ。
新卒で某オレンジ色のWebで地図な会社に入って、退職するまでずっと地図画像作ったり、配信していたわけですが、ベクトルタイルに関してはちょっと距離を置いていたわけです。
数学的な素養がないとか、クライアントサイドでの描画を行なうための某フリーダムな言語およびデザイン等にあまり明るくないなど色々理由はあるわけですが。

普段タイルタイル言ってるだけではだめだなーって。そう思ったのでベクトルタイルを作ろうと思ったんです。SQLで(オカシイ)。

趣味でいじっている地理空間データはすべて PostgreSQL (PostGIS) に格納しているんですが、RDBMS 内でデータ管理からタイル生成、隙をみてタイル配信までできやしないかとふと思ったのでやってみたら挫折した記録をお届けしたいと思います。
特に PostgreSQL の JSON 関連機能が最近がんがん実装されてきたりしているので、色々夢が広がったんでしょうね、数時間前の私は。

ちなみに RDBMS は前職でも触ってたし、現職では RDBMS 担当なわけですが、どちらかというと DBA 的な部分の仕事が多くて実は SQL 力が低いです。

すごい頑張ったんですが、結果的にはタイルの境界線を発生させる関数しか現状できていません。以下です。


引数に zoom レベルをとります。RECORD 型で zoom レベルと LineString を返します。使い方はこんな感じ。

insert into tileline (zoom, geom) select zoom, geom from create_tileline(16) as (zoom integer, geom geometry);

QGIS で表示するとこんな感じ。

zoom レベルにはそのラインが存在する最大レベルが入ります。この値によって以下のような条件を与えれば、zoom レベルでの使い分けが可能になります。
select * from tileline where zoom <= 10
上の画像は zoom レベル16のラインですが、QGISで zoom <= 13 のフィルタをかけると↓になります。

後は、対象のフィーチャと重なるタイルのラインを Geometry Collection か MultiLineString かなんかにして ST_Split でタイル化完了!とか思ってたんですが、ST_Split が LineString しか取れないらしいです...orz
なんで MultiLine とかじゃないとダメかというと、切断に使うタイルのラインが複数ある場合一回で切らないと処理が超面倒くさそうというかこの短時間では実装できなさそうだったというだけなんですが...orz
MultiLineString を ST_Union する手もあるんですが、処理時間やTempな領域とはいえ一時的にかなりリソースを食いそうだったんでやめました。
やりたかったイメージ↓
ここで Split!!!
というわけでかなり道半ばですが、SQL力向上のためにも継続してやっていきたいと思います。
来年は OPEN KWASAKI で地元貢献目標で頑張ります。