PMTiles からタイル一覧とサイズを取得する（国土地理院最適化ベクトルタイルを例に）

はじめに

以前、国土地理院が提供している地理院地図Vector（ベクトルタイル） や標準地図（ラスタタイル）のタイルサイズの分布を観察したことがあります。このときは、国土地理院から提供されているタイルの z/x/y 座標とサイズ等の一覧（地理院タイル目録）を利用して、タイルのサイズを確認しました。

一方、国土地理院では、地理院Vector の他に、2022年9月6日より「最適化ベクトルタイル」という設計改良を施したベクトルタイルを提供しています。最適化ベクトルタイルの地理院タイル目録は提供されていないため、上記の記事では、地理院地図Vector のタイルサイズの分布の観察にとどめました。

しかし、2023年8月30日より、最適化ベクトルタイルの PMTiles 版の提供が開始されています。PMTiles については、以前記事で紹介したことがございますが、一つのファイルの中にタイルデータ一式が格納されており、ファイル内の一部のデータだけを取得できる HTTP 範囲リクエストを用いてタイルデータを取得することが可能なデータ形式です。

PMTiles で範囲リクエストをする場合、本の目次に相当するような部分（directory）を読み取ることになります。この directory には、目的のタイルが PMTiles のどのデータ領域にあるのかを示しているのですが、こちらを分析することで、地理院タイル目録のように、タイルサイズの分布を把握するのに活用できるのではないかと試してみました。今回は、PMTiles の directory の内容を取得する方法と、実際に最適化ベクトルタイルのタイルサイズの分布を観察した結果を記事にまとめてみたいと思います。

2024/01/23 追記
当初、RunLength が1よりも大きい場合、最初のタイルの情報のみ出力するプログラムとなっていましたが、RunLength に応じて、連続するタイルの情報を出力するように修正しました。それに伴い、グラフや数字も修正しています。ただし、考察の対象とした最適化ベクトルタイルの ZL15の場合、RunLength が1よりも大きいのは１ケース（15/29261/12324と15/29261/12325）しかありませんでしたので、ほとんど影響ありません。各 ZL のヒストグラムは変化しています。また、本修正に合わせて、RunLength 関係の考察を少しですが追加しました。

PMTiles から directory の内容を取り出す

PMTiles から directory を取り出すには、protomaps/PMTiles の参照実装を活用するのが手っ取り早いと思います。私は、WSL と Node.js 上で実行しました。

まずは、PMTiles のモジュールをインストールします。

npm add pmtiles

あとは、以下のようなプログラムを組んで directory の内容を取得し、タイルの z、x、y 座標とデータサイズ（length）を csv で出力しました。（もう少し効率の良い方法があるかと思いますが、一応これで取得することはできるかと思います……。抽出して確認したところ、間違ってはいなさそうでした。また、PMTiles の v3 のみの対応です。）

const pmtiles = require("pmtiles");
const fs = require("fs");

const url = "https://cyberjapandata.gsi.go.jp/xyz/optimal_bvmap-v1/optimal_bvmap-v1.pmtiles";
const file = "stat.csv";

// 書き込み先のファイルを作成する。
fs.writeFileSync(file, "z,x,y,length" + "\n");

const p = new pmtiles.PMTiles(url);
const cache = new pmtiles.SharedPromiseCache();
const source = new pmtiles.FetchSource(url);

// Root Directory 又は Leaf Directory を解析し、CSV へ書き込む。
const tileStat = (directory, header) => {
  let csv = "";
  directory.forEach( entry => {
    if(entry.runLength > 0){
      // RunLength が0よりも大きければ、タイルの情報を取得する。
      // RunLength の分だけ、同じデータのタイルが続くので、
      // tileId を1ずつ加算しながら、同じ length を記録する。
      for( let i=0; i < entry.runLength; i++){
        const length = entry.length;
        const xyz = pmtiles.tileIdToZxy(entry.tileId + i);
        csv += `${xyz[0]},${xyz[1]},${xyz[2]},${length}` + "\n";
      }
    }else{
      // run length が0であれば、Leaf Directory を読みに行く。
      console.log(entry);
      d_o = header.leafDirectoryOffset + entry.offset;
      d_l = entry.length;
      cache.getDirectory(source, d_o, d_l, header).then( leaf => {
        // 再帰的に解析していく
        tileStat(leaf, header);
      });
    }
  });
  // Root Directory 又は Leaf Directory 毎に csv へ書き込む。
  fs.appendFileSync(file, csv);
}

// Header を読み取り、Root Directory からタイルのデータを順次読み取る。
p.getHeader(source).then( header => {
  console.log(header);
  let d_o = header.rootDirectoryOffset;
  let d_l = header.rootDirectoryLength;
  cache.getDirectory(source, d_o, d_l, header).then( leaf => {
    tileStat(leaf, header);
  });
});

これで、以下のような csv データが得られます。

z,x,y,length
4,14,7,220
4,14,6,17340
4,14,5,6177
...
16,55560,28088,84
16,55560,28089,71
16,55560,28090,85

こちらのデータを用いれば、以前の記事にて、地理院タイル目録を用いたときと同様な方法で、タイルサイズの分布などを観察することができるようになります。

最適化ベクトルタイルのタイルサイズ分布の特徴

実際に、PMTiles から取得したデータを使って、地理院地図Vector のものと比べながら、最適化ベクトルタイルのタイルサイズ分布の特徴を観察していきたいと思います。なお、

• どちらも、観察対象は ZL15 に限定する
• 最適化ベクトルタイルは圧縮あり（gzip）のサイズなのに対して、地理院地図Vector は圧縮なしのサイズである

という点にご注意ください。

最適化ベクトルタイルのタイル総数は544,822枚、地理院地図Vector のタイル総数は544,845枚でした。

圧縮の有無があるので単純比較はあまり意味がありませんが、ZL15 の合計タイルサイズは、地理院地図Vector が約12 GB （圧縮なし）なのに対して、最適化ベクトルタイルは約4.4 GB（圧縮あり）でした。また、最大のタイルサイズは、地理院地図Vector が約535 kb （圧縮なし）なのに対して、最適化ベクトルタイルは約73 kb（圧縮あり）でした。

タイルサイズの地理的な分布について両者を比較すると、顕著な違いがみられます。すなわち、地理院地図Vector の場合、都市部のデータが非常に大きくなる傾向にありましたが、最適化ベクトルタイルの場合、むしろ等高線が密集するような山間部のデータが大きくなる傾向にあります。

（図では、タイルサイズのデータをベクトルタイル化して、MapLibre GL JS で表示しています。赤いほど、相対的にタイルサイズが大きいことを示しています。最適化ベクトルタイルと地理院地図Vector のスケールは同一ではありません。）

→Web 地図で最適化ベクトルタイルのサイズ分布を観察する

→Web 地図で地理院地図Vector のサイズ分布を観察する

最大のタイルサイズとなる場所を比較しても、地理院地図Vector の場合、普通建物が密集するようなタイル（茨木市）なのに対して、最適化ベクトルタイルでは、等高線や地形の表現が豊富な穂高岳周辺のタイルでした。

両者の比較については、過去の記事でもまとめています。

参考資料によれば、最適化ベクトルタイルでは、建築物の外周線データを削除している他、ZL14-15では一定面積以下の建築物を削除しているとのことで、都市部周辺のタイルサイズを徹底的に削減しているのに対して、（ZL14では岩・崖・雨裂等を削除しているようですが）ZL15では、等高線や地形のデータは削除していないと見られますので、相対的に山岳部のタイルデータが大きくなっていると考えられます。

地理院地図Vector のタイルサイズと最適化ベクトルタイルのタイルサイズの分布を X 座標帯に沿って比較してみると、（最適化ベクトルタイルは圧縮後のサイズで評価しているとはいえ、）地理院地図Vector のタイルサイズのばらつきに対して、最適化ベクトルタイルのタイルサイズのばらつきは小さく、全体的に一定のタイルサイズに抑えられていることが分かります。国土地理院によれば、（おそらく圧縮なしで）256 kb/タイル以下に抑えるように設計したとのことです（参考資料）。

地理院地図Vector と最適化ベクトルタイルでどのようにサイズが変わったか、タイル座標が同じものどうし1対1で散布図にプロットしてみると以下のようになります。

目視観察では、

1. 地理院地図Vector では非常にサイズが大きいが、最適化ベクトルタイルではサイズが小さい（一定以下に抑えられている）タイル群
2. 地理院地図Vector のサイズと最適化ベクトルタイルのサイズが比例しているタイル群

の２つに分けられそうですが、前者が建物が多い都市部のタイル、後者が等高線が多い山間部のタイルではないかと考えられます。

この結果からみても、「建物データに手を入れ、小さな建物をズームレベルに応じてデータから省略する」という最適化ベクトルタイルにおける戦略は、タイルサイズの削減と均一化に大きく貢献しているのではないかと考えられます。

各ズームレベルのヒストグラム

最適化ベクトルタイルだけとなりますが、参考に各ズームレベルごとにヒストグラムも作成してみました。最大サイズも記載しております。なお、これまで同様、圧縮後のサイズであることにご注意ください。

データが重すぎて Excel では処理できないので、R を利用しました。以下はそのコードになります。

df <- read.csv("stat.csv")
par(mfrow=c(3,5))
for( i in 4:16 ){
  df2 <- subset(df, z==i)
  df2$length <- df2$length / (1024)
  maxsize = max(df2$length)
  hist(df2$length, breaks=seq(0, maxsize+4, by=4),
       main = paste("ZL=", i, ", n=", nrow(df2), ", max=", round(maxsize), " kb", sep=""))
  abline(v=maxsize, col=2)
}

ズームレベル毎の RunLength の活用

記事公開当初は、RunLength が1より大きい場合、最初のタイルの情報のみを出力していましたが、後日（2024/01/23）、RunLength の数だけ、連続するタイルの情報を出力するように修正しました。ついでに最適化ベクトルタイルについて、各ズームレベル毎に傾向を分析してみると、以下のように RunLength が活用されていることが分かりました。

ズームレベル	RunLength が2以上の回数	最大の RunLength
8	22	3
9	189	15
10	637	63
11	762	34
12	2129	139
13	4746	558
15	1	2
16	1590	2

RunLength が2以上の回数や RunLength の最大値が大きいほど、RunLength を活用して PMTiles のデータ量やキャッシュを活用した通信料の削減がされていると言えます（別のタイルでも PMTiles 上のデータが同じであれば、同じ HTTP 範囲リクエストになりますのでキャッシュの活用が期待できます）。ZL13においては、RunLength が558に達するところもあります。

ZL8～13までは RunLength が活用されているのに対して、ZL14～15はほとんど利用されておらず、ZL16でまた増加するという点が面白いです。この「RunLength が2以上の回数」はほぼ「水域しかないタイルの枚数」に影響を受けると考えられますから、ZL13や16は水域部分のタイルを広めに作成している一方、ZL14,15は余分な水域のタイルを作っていないと考えられます。等高線のデータなどを見ていると分かりやすいですが、地理院地図Vector や最適化ベクトルタイルでは、ZL13と14を境にデータの体系が異なっているようで、ここでもそれが示唆されます。

レポジトリ・サイズ分布確認用地図サイト

利用したコード等は、以下のレポジトリに置いています。

また、PMTiles から取得したデータサイズの分布を地図で見られるようにしています（ZL15 のサイズのみ）。データ量や動作が結構重いのでご注意ください。

おわりに

今回、PMTiles からタイルセットの統計情報を取得する方法を試してみました。地理院タイル目録が提供されなくても、PMTiles が提供されている場合、似たようなことができることを学べました。

（そもそも地理院タイル目録は、効率的なタイルダウンロード・同期を目的としているようですが、PMTiles の場合、PMTiles の URL を叩けばタイル一式が手に入ってしまうわけで、その点においては、この記事のような directory の分析等は必要ありません……。）

また、地理院Vector と最適化ベクトルタイルを比較することができ、設計改良の効果を面的に把握することができました。改めて、都市部でのタイルサイズの削減効果には目を見張るものがありますね。