PostGISのST_AsMVTでバイナリベクトルタイル

これはベクトルタイル Advent Calendar 2017 12/22の記事です。

MIERUNEとベクトルタイル

OpenMapTiles

株式会社MIERUNEでは、MIERUNE地図としてタイル地図画像を配信しています。
f:id:waigani:20171220125421p:plain
何をベースに作っているかというとOpenMapTilesになります。
OpenMapTilesの作成手順に沿っているので、間でmbtilesを作成していて、バイナリベクトルタイルになっています。
諸事情によりバイナリベクトルタイルでの公開していないのですが、近々公開したいとは思っています。

改良したい点

MIERUNE地図として公開している範囲としては、日本+台湾になります。
本年度中に範囲をどんどん増やしていくつもりだったのですが、ストップしています。
なぜストップしているかというと、巨大なmbtilesを作るのが時間はかかるし、手間がかかるからです。
Klokan Technologiesのエンジニアにも聞いたことがあるのですが、「分散する以外にないガンバレ」との返答をいただきました。そうですよね、それしかないですよね。
postserveも出てきたので、mbtilesすっとばしてPostGISからの直接配信の方向にも力をいれていくのかなと思ったら、「あれKlokan Technologiesの人じゃないんだよね」とも言っていたので、やっぱり頑張ってmbtiles作るのかな。。。
おまけにスキーマ変更のコストも高いんです。Debugどうしてるの？と聞いたら「興味深いトピックスだよね（笑」で終わったので、エンジニアの地道な努力に支えられているのが現状です。

MIERUNEとしては、ベースマップに関するデータだけではなく、変更が頻繁にあるようなデータについてもベクトルタイルで配信していくことを想定しておきたいです。SQLだけでほぼスキーマ定義を変更することも目的に、PostGISからアクセスに応じて都度データを取得しての配信についても色々お試し中です。

ST_AsMVT

いろいろアプリケーションは出てきているのですが、もう１段低レベルなところから確認しましょう。
ST_AsMVTを使ってみます。
https://postgis.net/docs/ST_AsMVT.html

インストール

PostGIS ver2.4から、使えるようになっています。
ただし、libprotobuf-c が必要です。PostGIS ver2.4ですぐ使えるのかなとpgadminで関数一覧をみても、
f:id:waigani:20171220132921p:plain
st_asmvtgeomはあるのに、st_asmvtがない、となるかもしれません。
PostGISを、ソースからインストールしましょう。configure時に、protobu-c supportがされていることを確認しておきます。
f:id:waigani:20171220133541p:plain

tileから経緯度

tileのxyzから経緯度を求める関数は必要です。OpenStreetMap wikiの記事"Slippy_map_tilenames"に、plpgsqlでの例が載っていますので、使わせてもらいましょう。
http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#PostgreSQL

使い方例

テーブルjapanに、国土数値情報の行政区域が入っています
x:3638, y:1618, z:12で切り出します
st_makeenvelope()でタイル矩形を作る
st_intersection()でタイル矩形でjapanを切り出し
st_asmvtgeom()でタイル内座標(4096)に変換
st_asmvt()でベクトルタイル化

の例になります。

select st_asmvt(v, 'test', 4096, 'geom')
from (
select j.n03_001, j.n03_002, j.n03_003, j.n03_004, j.n03_007, st_asmvtgeom(st_intersection( j.geom, e.geom), e.geom,  4096, 0, true)  as geom
from
 japan j,
 (select st_makeenvelope(tile2lon(3638,12), tile2lat(1618,12), tile2lon(3639,12), tile2lat(1619,12), 4326) as geom) e
 where
  st_intersects(j.geom, e.geom)
) v

確認

上記のSQLをtest.sqlというファイル名で保存しておいたとして、SQLの実行結果を一旦ファイルに書き出し。書き出したファイルをバイナリ化してから、geojsonに変換して中身を確認してみます。バイナリベクトルタイルからgeojsonへの変換は、vt2geojsonを使います。

psql -f test.sql -tAz0 "データベース名" > 12-3638-1618.txt
xxd -r -p 12-3638-1618.txt > 12-3638-1618.pbf
vt2geojson -x 3638 -y 1618 -z 12 --layer test 12-3638-1618.pbf > 12-3638-1618.geojson

geojsonをQGISに表示して、タイル番号と一緒に表示して確認してみます。
意図したタイルのデータを、切り出し出来ています。
f:id:waigani:20171220135757p:plain

補足(2017.12.28追記)

誤解を与えそうなので。mapbox vector tile specificationでは、タイルの外側にバッファを取ります。
上記例だと、きれいにタイルに合わせた絵にしようと思ってst_intersection()で切り取っているのですが、ここは本来要らないです。

  st_asmvtgeom(j.geom, e.geom,  4096, 256, true)

のようにして、4番目の引数がバッファサイズ、5番目の引数がclipするか、を指定してあげます。
ただ、st_asmvtgeomに渡してあげるj.geomは、このタイルに重なってくる範囲に限定してあげましょう。

ここで求人

MIERUNEで、一緒にベクトルタイルしたいエンジニアいませんか？
表立って求人は出していないのですが、ご興味のある方がいらっしゃればいつでもお声掛けください。
札幌に事務所はありますが、どこで勤務して頂いても構いません。

waigani's diary

QGISを中心にFOSS4Gをいじくる