国土地理院がベクトルタイル提供実験として、全国版のバイナリベクトルタイルの公開を行っています。
国土地理院のバイナリベクトルタイルの仕様（各ズームレベルごとに含まれているレイヤなど）については下記にまとめられています。
Mapbox Vector Tiles (MVT)
国土地理院ベクトルタイルで使われている Mapbox Vector Tiles 形式は Mapbox 社が策定した XYZ タイル型のベクトルデータ形式で、 Protocol Buffer 形式を利用した形式です。
Mapbox Vector Tiles の仕様は以下で確認できます。
- vector-tile-spec/2.1 at master · mapbox/vector-tile-spec
- vector-tile-spec/2.1 at master · madefor/vector-tile-spec （日本語訳）
すでにいろいろなツールで対応しており、 GDAL/OGR でもバージョン2.3以降は対応しています。そのため ogr2ogr を使ったファイル形式変換や QGIS 読み込みは普通にできます。
が、今回は vector_tile_base を使って Python を使ってパースを行ってみます。
vector_tile_base
Mapbox 社からでている、 Python 用のライブラリです。
既存の MVT ファイルをパースするには、上記の Decode のサンプルを参考にします。（というより、後述のソースは Decode サンプルのまんまです。）
データの取得
サンプルで使用するためのファイルを取得します。
$ wget -O 13_7283_3213.pbf https://cyberjapandata.gsi.go.jp/xyz/experimental_bvmap/13/7283/3213.pbf
レイヤの確認
import vector_tile_base
with open(r'./13_7283_3213.pbf', 'rb') as f:
data = f.read()
vt = vector_tile_base.VectorTile(data)
for l in vt.layers:
print(l.name)
label
building
waterarea
transp
symbol
road
ジオメトリの確認
import vector_tile_base
with open(r'./13_7283_3213.pbf', 'rb') as f:
data = f.read()
vt = vector_tile_base.VectorTile(data)
layer_map = {l.name : l for l in vt.layers}
if 'road' in layer_map:
layer = layer_map['road']
print('extent : %d' % layer.extent)
for f in layer.features:
print(f.get_geometry())
break # 全地物を表示すると大変なので、ここではこれ以降の地物はスキップ
extent : 4096
[[[169, 45], [127, 43], [54, 34], [0, 11], [-36, 0], [-80, -8]]]
ここで得られた座標の値は前述のドキュメント（仕様）にあったとおり、整数値で格納されています。そして、タイルの幅は extent であり、座標の単位は 1/extent となります。たとえば extent が4096のタイルであれば、イメージとしては4096×4096ピクセルの画像を想像するとよいでしょう。
get_geometry() で得られる $\left[ 0, \ extent \right]$ の座標値を $extent$ で割ると、 $\left[ 0, \ 1 \right]$ に正規化されたタイル内座標を計算することができます。ただし、適切に線を引くためタイル外のノード情報を持つこともあります。
座標の補正
先程得られたジオメトリの座標は対象のタイルにおける座標でした。実際に GIS データとして使用するには不十分です。
タイル地図は、赤道半径 $a =$ 極半径 $b = 6378137$ [m] の真球メルカトル（疑似メルカトル） EPSG:3857 において、縦横 $S = 2 \pi a = 40075016.685$ [m] の正方形となるように切り取った範囲（東経180度〜西経180度、北緯85.01度〜南緯85.01度）を対象範囲とし、ズームレベルに従いタイルの分割を行っています。
ズームレベル z のときの、個々のタイルの幅は $\frac{2 \pi a}{2^z}$ [m] となります。ただし、単位はメートルですが、あくまでも EPSG:3857 の投影地図上の単位としてであり、実際の距離とは一致しないことに留意が必要です。
タイル座標が z/x/y であるタイルの EPSG:3857 における北西角の座標 $\left( x_\mathrm{offset}, \ y_\mathrm{offset} \right)$ は次式となります。
\begin{array}{lll}
a & \!\!\! = & \!\!\!\!\!\!\!\! 6378137 \\
S & \!\!\! = & \!\!\!\!\!\!\!\! 2 \pi a \ = \ 40075016.685 \\
x_\mathrm{offset} & \!\!\! = & \!\!\!\!\! \left( x - 2^{z-1} \right) \times \displaystyle\frac{S}{2^z} \\
y_\mathrm{offset} & \!\!\! = -& \!\!\!\!\! \left( y - 2^{z-1} \right) \times \displaystyle\frac{S}{2^z}
\end{array}
$y_\mathrm{offset}$ の方にマイナスがついているのはタイル座標もタイル内座標も南方向が正なのに対し、 EPSG:3857 の南北軸は北方向が正だからです。
また get_geometry() で得られる正規化前のタイル内座標を $\left( s, t \right)$ とおくと、 EPSG:3857 における対象点の座標 $\left( x_{3857}, y_{3857} \right)$ は次式となります。
\begin{array}{lll}
x_{3857} & \!\!\! = & \!\!\!\!\! \displaystyle{ \left( \frac{s}{extent} + \left( x - 2^{z-1} \right) \right) \frac{S}{2^z} } \\
y_{3857} & \!\!\! = -& \!\!\!\!\! \displaystyle{ \left( \frac{t}{extent} + \left( y - 2^{z-1} \right) \right) \frac{S}{2^z} }
\end{array}
これにより GIS で使える座標を計算することができました。あとは煮るなり焼くなり好きに料理しましょう。
とはいえ
GDAL/OGR が対応していますので、 QGIS でも普通に開くことができますし、下記のように ogr2ogr を使って他形式に変換することもできますので、わざわざ Python で読んだり、変換したりといったことはあまり必要なさそうです。
ogr2ogr を使う場合、こんな感じで変換します。
$ ogr2ogr -f GPKG 13_7283_3213.gpkg -nlt PROMOTE_TO_MULTI 13_7283_3213.pbf
ただし GDAL/OGR の仕様でファイル名が z/x/y.ext z-x-y.ext z_x_y.ext のいずれかでないとタイル座標が得られず、正しく変換することができません。なお拡張子 ext は mvt でも pbf でも可能です。
QGIS で読み込む際も同様で、ファイル名を変えないと正しく読み込むことができません。
座標の精度
ところで、座標値は整数値として格納され、 extent を使ってタイル内座標を計算しました。逆に言うと 1/extent 単位で離散化されていると言えます。
extent の値を大きくすることで間隔を小さくすることもできますが、いずれにしても、規定以上の精度を持つことができない仕様になっています。これによって座標値としての圧縮はもちろん、細かすぎるノードを削除しパスを簡略化するができるので、データ量の削減を行うことができるのでしょう。
しかし extent = 4096 の場合、 z = 10 のときの座標間隔は（ EPSG:3857 の投影地図における）約9.6 m と結構大きな値となり、座標値としての誤差は無視できません。あくまで縮尺により適切なズームレベルのタイルが読み込まれることを前提とするものであり、細かいデータが不要な縮尺の場合は、データそのものがばっさり切り捨てられるデータ形式であると認識する必要がありそうです。