ベクトルタイル (Mapbox Vector Tile) をPythonで自力で生成してみた #Python

ベクトル地図タイルのデファクトスタンダード形式である Mapbox Vector Tile (MVT) のデータを試しに自力で生成してみたので、その要点を解説してみます。なお 「自力で生成」 というのは、MVTの実体である Protocol Buffers のデータを組み立てるところからフルスクラッチでやるという意味です。

最終的に上図のようなポリゴンを含むベクトルタイルを生成して配信できました。サンプルコードではPythonを使って以下のようなことを行っています。

（GeoJSONのデータをもとに）大量のポリゴンを含むMVTを自分で組み立てて生成する。各地物には属性情報も付与する。
現在時刻などの赤い文字列の部分もLineStringsで動的に描く。
PythonのFastAPIで配信する。

サンプルコードのリポジトリは以下です：

まず前提として：

多くのユースケースでは tippecanoe などの既成ツールを使ってMVT形式への静的変換を行うことがほとんどでしょう。
一部のデータベースエンジンにはMVTを動的に生成する機能が実装されています（例えばPostGISの ST_AsMVT、Elasticsearchの Vector tile search API など）。基本的な用途でMVTを動的に生成したい場合はこれらに任せるのが手軽でしょう。
しかし特殊な状況では、自前でMVTの生成ができると可能性が広がるかもしれません。（また、簡単な生成処理を実装をしてみてMVTの仕をの理解するのも楽しいことです（？））
いくつかの言語では、MVTの作成を支援するライブラリがOSSライセンスで転がっています。

MVTの仕様は http://mapbox.github.io/vector-tile-spec/ で公開されています（GitHubリポジトリは https://github.com/mapbox/vector-tile-spec/tree/master）。

MVTは、Googleの Protocol Buffers (protobuf) という規格の上に作られており、Mapboxによって定義された .protoファイルをもとに protobuf のツール群で各言語用のシリアライザや型定義などを自動生成できます。シリアライズ時の難解なバイナリ操作の大部分は protobuf がやってくれます。

protobuf という整備されたエコシステムのうえに構築されているため、要点さえつかむことができれば、MVTを自分で出力すること自体は実はそれほど難しい作業ではありません。

PythonでのMVT生成の実装例

今回はMVTの生成をPythonでお手軽に実装してみました。サンプルリポジトリは以下にあります。

（このサンプルコードを見てもらうとわかりますが、それほど複雑なコードではないですし、行数も多くないです）

今回はPythonで実装しましたが、MVTの生成には大量の演算と大量のオブジェクト生成が伴うので、パフォーマンスの面でPythonはまったく向いていません。実用性なものを作りたい場合は静的な言語で実装したほうがよさそうですし、速度や、メモリ使用のスパイクを抑制したいならGCのない言語を使うのが良いかもしれません。Protocol Buffersのコードジェネレータが用意されている言語であれば、どの言語でも同様に実装できるはずです。

protoファイルのコンパイル

Mapbox Vector Tile Specification で公開されている .proto ファイルをコンパイルすることで、好きな言語でMVT形式（実体はProtocol Buffers）のデータを扱えるようになります。MVT仕様の最新バージョンは 2.1 です（MVT Specificationリポジトリの別のブランチに3.xを作ろうとした痕跡が見受けられますが進展はなさそうです）。

この .proto ファイルを Protocol Buffers 用のツールに通すと型定義などが自動で生成されます。Pythonの場合はGoogleが提供する公式ツールが使えるので、それをインストールしましょう。（インストール方法は様々ありますが今回はMacのHomebrewを使いました）

> brew install protobuf

以下のコマンドで、proto ファイルをPythonモジュールにコンパイルします：

> protoc -I=. --python_out=. --pyi_out=. vector_tile.proto

これによって vector_tile_pb2.py というファイルが生成されます。このモジュールをPythonで利用することで、MVTのデータを記述したり、バイナリにシリアライズしたりできます。.py ファイルと同時に .pyi（型スタブ）ファイルも生成されるため、VSCodeなどのエディタでコード補完や型チェックも効きます。

protoc を実行する際に警告が表示されるかもしれません。無視しても大丈夫ですが、警告を抑制したい場合は vector_tile.proto の先頭に以下の行を追加してください：

// この proto ファイルの構文バージョンを明示的に記述しておく
syntax = "proto2";

MVTの構造

まずMVTファイルの構造を簡単に眺めてみましょう。MVTの構造を（実際は Protocol Buffers ですが）YAML風に記述すると以下のようになります。

tile:
  layers:
    - version: 2
      name: "レイヤ名1"
      keys: ["都道府県", "都道府県コード"]
      values: [Value("栃木県"), Value("群馬県"), Value(9), Value(10)]
        # ↑地物の属性として使われるKeyやValueはここにまとめて格納しておきます
      extent: 4096

      features:
        - id: 0 # 地物ID
          type: GeomType.POLYGON # POINT, LINESTRING or POLYGON
          tags: [0, 0, 1, 2]
            # ↑ レイヤ上に定義された Key と Value を整数インデクスで参照します
          geometry: [......] # 特殊な整数列でジオメトリをエンコードする（後述）

        - id: 1
          type: GeomType.POLYGON
          tags: [0, 0, 1, 2]
          geometry: [......]

        - id: 2
          # ...(略)...

    - version: 2
      name: "レイヤ名2"
      keys: [......]
      values: [......]
      extent: 4096
      features:
        # ...(略)...

意外とシンプルな構造であることがわかると思います（パッと見はGeoJSONみたいな感じですね）。いくつか要点だけ触れておきます：

tags — 各地物の属性です。MVTにおいては、属性のキーや値をレイヤ直下の配列にまとめて格納しておいて、それを各地物から整数インデクスで参照するというアプローチをとります。多くのGISデータでは、属性のキーや値は同じものが何度も何度も使われるため、このようにすることで空間効率が圧倒的に良くなる（サイズが縮む）というわけです。
extent — タイル1枚の座標の範囲です。extent=4096 の場合は (0, 0) が左上で (4096, 4096) が右下を表します。(4096, 4096) も範囲内に含まれることに注意してください。MVTでは、タイル上の点は浮動小数点数ではなく整数で表します。また、ジオメトリはこの範囲を一定程度はみだすことが意図的に許されています（つまり、マイナス値になったり、extentを超えたりしてもOKです）。タイルからはみだした内容は描画時にクリップされます。
geometry — 地物がもつジオメトリの形状は、MVT専用の形式で「符号なし整数の列」としてシリアライズされます （後述）。

詳しい仕様は、MVTの仕様書を参照してください。

Pythonで protobuf データを扱う

簡単に雰囲気だけ示しておきます。具体的な例はサンプルコードを見てください。

Googleのツールによって自動生成されたコードを用いて、以下のように Feature や Tile を作り、Protocol Buffers のバイナリにシリアライズできます。

import vector_tile_pb2 as pb

# 以下のように Feature や Tile を作れる
features = pb.Tile.Feature(...)
layer1 = pb.Tile.Layer(
        version=2,
        name="example",
        keys=keys,
        values=values,
        features=features,
        extent=extent,
    )
tile = pb.Tile(layers=[layer1])

# protobuf のバイナリ形式にシリアライズする
serialized = tile.SerializeToString()

ジオメトリをエンコードする

現在のMVTの仕様では、ジオメトリとして Point, LineString, Polygon の3種類を描くことができます。

MVTではジオメトリを32ビット符号なし整数の配列として表現します。この符号なし整数の列は、コマンド整数と、座標値 dx, dy で構成されます。

コマンド整数

符号なし32ビット整数のうち、下位3ビットが種別を表し、上位29ビットがそのコマンドで扱う点の数を表します。以下の3種類のコマンドがあります。

MoveTo コマンドは整数値 (点の数 << 3 | 1) で表現する。
LineTo コマンドは整数値 (点の数 << 3 | 2) で表現する。
ClosePath コマンドは整数値 (1 << 3 | 7) で表現する。

（<< はビットシフトです）

座標値の表現

MoveTo コマンドとLineTo コマンドは、その直後に「点の数」のぶんだけ x 座標と y 座標のペアを伴います。主な留意点が2つあります：

座標そのものではなく直前に記録した座標からの差分を記録します
- 前の座標からの差分値を記録することで出力されるバイナリのサイズを削減します (Protocol Buffersは小さな整数値を小さなバイト数で表現できる仕組みになっており、それを活かすためです）
- 例えば前の座標が (15, 20) で次の座標が (50, 15) の場合は (35, -5) を記録します。
- 差分による座標の記録は地物ごとにリセットされます。つまり新たな地物では起点を (0, 0) にリセットして座標の差分を求めます。
符号付き整数の座標を、zigzagエンコードという方式で符号なし整数に変換してから格納します
- 符号付き整数を (a << 1) ^ (a >> 31) という演算によって符号なし整数にしたものを格納します。
- 例えば 0 → 0, -1 → 1, 1 → 2, -2 → 3, 2 → 4 のように変換されます。zigzagエンコードは、小さな負数が小さな符号なし整数にマッピングされるのが特徴です。

Point / MultiPoint の描き方

PointやMultiPointは以下のような符号なし整数の列としてエンコードします：

MoveToコマンド 点の数 << 3 | 1
点の数だけ dx, dy を描く（zigzagエンコード）

LineString / MultiLineString の描き方

以下を線 (LineString) の数だけ繰り返します：

MoveToコマンド: 1 << 3 | 1
始点 dx, dy（zigzagエンコード）
LineToコマンド: 残りの点の数 << 3 | 2
残りの点の数だけ dx, dy を追加する（zigzagエンコード）

Polygon / MultiPolygon の描き方

1つのポリゴンは、1つの外部リングと、0個以上の内部リングから構成されます（これは一般的な概念なのでここでは触れません）。

以下をポリゴン (Polygon) の数だけ繰り返します：

外部リング（座標列は時計周りでなければなりません）
1. MoveToコマンド: 1 << 3 | 1
2. 始点 dx, dy（zigzagエンコード）
3. LineToコマンド: 残りの点の数 << 3 | 2
4. 残りの点の数だけ dx, dy を追加する（zigzagエンコード）
5. （なおGeoJSONやWKTのように終点を始点と一致させてはなりません。冗長な終点は不要です）
6. ClosePathコマンド: 1 << 3 | 7
内部リングを0個〜複数個記録する（座標列は反時計周りでなければなりません）
- 記録の方法は外部リングと同様です

これを繰り返す。

外部リング (External Ring) か内部リング (Interior Ring) かは、座標の列が時計周りか反時計周りかで区別されます。

地物に属性を付与する

上でも言及しましたが、MVTでは、地物の属性のキーや値はレイヤ上の配列にまとめて格納しておき、それらのキーや値を各地物から整数インデクスで参照するというアプローチをとります。

例えばMVTに以下のようなデータがあったとします：

# レイヤ上に以下のように記録されているとする...
keys: ["都道府県", "都道府県コード"]
values: [Value("栃木県"), Value("群馬県"), Value(9), Value(10)]

# そして地物の tags には以下のように記録されているとする...
tags: [0, 0, 1, 2]

地物に記録された tags=[0, 0, 1, 2] は、以下のような属性として解釈されます：

# tags: [0, 0, 1, 2] の解釈
# tags にはキーと値へのインデクスが交互に記録されているので:
[(key=0, value=0), (key=1, value=2)]

# これらの整数インデクスで、キーと値を解決するとこうなります:
tags:
  都道府県: "栃木県"
  都道府県コード: 9

現実のGISデータでは、属性キーはレイヤ内ののほぼ全ての地物に登場することがほとんどですし、値についてもカテゴリカルな属性値は何度も繰り返し出現するでしょう。そのため、このような方式で記録することによって、生成されるタイルサイズを大きく削減できます。

MVTを生成するときは、このようなエンコードを実現するための小さなユーティリティを実装してあげればよいでしょう。サンプルコードでは TagsEncoder というクラスがその処理をしています。

おまけ: 文字を線で描く

↓この文字はサーバ側でLineStringとして描いています。

線（ストローク）だけで実現されたフォントとして Hershey fonts というものがあります。

Pythonでは以下のライブラリなどでHershey fontsの座標データを利用することができます。上図の文字列はこのデータを使ってLineStringを描くことで描画しています。

その他に必要なもの

上で述べたことは主にタイル上でお絵描きをする方法ですが、現実的なタイル生成を行う場合は以下のようなものが適宜必要になると思います：

ウェブメルカトルの座標とタイル上の座標を相互変換するなどの計算
ジオメトリ群をタイルの正方形で切り抜く処理：
- サンプルプログラムでは Shapely（GEOSをPythonから使いやすくしたもの）を使って切り抜いています。
- Mapbox による geojson-vt の実装が参考になるかもしれません。
- なお、タイルぴったりの枠で地物の切り抜きを行うと表示上の重大な問題が起きるため、タイルよりもひと回り大きな枠で切り抜くのが一般的です。tippecanoeのbufferオプションや、PostGISのST_AsMVTGeomのbufferパラメータも同様の目的で使われています。
空間インデクス：
- 解きたい課題に応じて、四分木、 R-tree、kd-tree などの適切な空間インデクスによるうまい工夫が必要になることが多そうです。