バイナリベクトルタイルを Python で読んでみた

国土地理院がベクトルタイル提供実験として、全国版のバイナリベクトルタイルの公開を行っています。

- [地理院地図｜ベクトルタイルとその提供実験について](http://maps.gsi.go.jp/development/vt.html)

国土地理院のバイナリベクトルタイルの仕様（各ズームレベルごとに含まれているレイヤなど）については下記にまとめられています。

- [gsi-cyberjapan/gsimaps-vector-experiment: 地理院地図Vector（仮称）提供実験のソース](https://github.com/gsi-cyberjapan/gsimaps-vector-experiment)

## 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](https://github.com/mapbox/vector-tile-spec/tree/master/2.1)
- [vector-tile-spec/2.1 at master · madefor/vector-tile-spec](https://github.com/madefor/vector-tile-spec/tree/master/2.1) （日本語訳）

すでにいろいろなツールで対応しており、 GDAL/OGR でもバージョン2.3以降は対応しています。そのため `ogr2ogr` を使ったファイル形式変換や QGIS 読み込みは普通にできます。

が、今回は vector_tile_base を使って Python を使ってパースを行ってみます。

## vector_tile_base

Mapbox 社からでている、 Python 用のライブラリです。

- [mapbox/vector-tile-base: A very basic and low level library written in python for encoding and decoding Mapbox Vector Tiles](https://github.com/mapbox/vector-tile-base)

既存の MVT ファイルをパースするには、上記の Decode のサンプルを参考にします。（というより、後述のソースは Decode サンプルのまんまです。）

### データの取得

サンプルで使用するためのファイルを取得します。

```shell-session
$ wget -O 13_7283_3213.pbf https://cyberjapandata.gsi.go.jp/xyz/experimental_bvmap/13/7283/3213.pbf
```

### レイヤの確認

```python
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
```

### ジオメトリの確認

```python
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)$ は次式となります。 

```math
\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)$ は次式となります。

```math
\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` を使う場合、こんな感じで変換します。

```shell-session
$ 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 と結構大きな値となり、座標値としての誤差は無視できません。あくまで縮尺により適切なズームレベルのタイルが読み込まれることを前提とするものであり、細かいデータが不要な縮尺の場合は、データそのものがばっさり切り捨てられるデータ形式であると認識する必要がありそうです。