ベクトルタイルのサイズ分布を観察する

## はじめに
ベクトルタイルは、ライスタイルとは異なり、地物の量によってファイルサイズの大きさのばらつきが大きくなると考えられます。特に、基本図のような様々な地物を含む場合は、それが顕著だと思われます。また、地物の多少は地域によっても異なるため、地域ごとにファイルサイズの大小の特徴が現れるかもしれません。

国土地理院では、[地理院地図Vector](https://github.com/gsi-cyberjapan/gsimaps-vector-experiment) の提供実験を継続していますが、同時に、全タイルの情報を記載した[地理院タイル目録（mokuroku）](https://github.com/gsi-cyberjapan/gsimaps-vector-experiment#mokuroku%E3%81%AE%E6%8F%90%E4%BE%9B)も提供されています。地理院タイル目録の仕様は以下の通りですが、タイル座標等の他、タイルのファイルサイズも記載されています。

https://github.com/gsi-cyberjapan/mokuroku-spec

今回は、この mokuroku の情報を用いて、地理院地図Vector のタイルサイズのばらつきについて視覚化しつつ、観察してみたいと思います。タイルサイズはズームレベル（ZL）によっても変化するのですが、取り急ぎ、**ZL15** を対象としてみます。

:::note warn
本当は、[最適化ベクトルタイル](https://github.com/gsi-cyberjapan/optimal_bvmap)で行いたかったのですが、まだ mokuroku が出ていないようですので、従来の地理院地図Vector で行いました。
:::

## ヒストグラムの観察

Excel を用いて、mokuroku の情報からヒストグラムを作成しました。設定の都合上、一部、幅が異なるビンがあるのでご注意ください。

ちなみに、全タイル（ZL15）は544,845枚、タイルサイズは合計で12,842,638,040 byte（約12 GB）でした。

![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/e72c1f92-6c22-6f91-807c-16a2f9566322.png)

こちらを観察すると、1 kb 未満の小さなファイルサイズのタイルが大量にあることが分かります。これは主に水域部分のタイルと思われます。長方形の水域のポリゴンが少数入っているだけなタイルです。陸上でも、小さなファイルサイズのタイルがないか探してみましたが、釧路湿原のような湿地部分くらいしか見つけることができませんでした。ちなみに、1 kb 以下のタイルは合計122,182枚でした。

また、1 kb 以上のタイルを見てみると、正規分布に近似できそうな、山がひとつの分布を示しています。最も多いのは20 kb 前後のもののようです。

これらの特徴から、タイルセットのデータサイズの代表値を出すのは、少し工夫が要りそうな気がします。単純に、平均値や中間値だと、小さなファイルサイズのタイルに引っ張られそうなので、一定の値（例えば1 kb）以下のタイルを除いた後のセットで計算した方が良いかもしれません。

逆に、大きなタイルを見ていくと、128 kb 以上のタイルも一定数存在します。全タイル（ZL15）は、544,845枚ですが、そのうち128 kb 以上のタイルが13,562枚、256 kb 以上のタイルが1,783枚存在するようです。最大のサイズは 548,380 byte（約534 kb）でした。

次に、これらの大きなファイルサイズのタイルの地理的分布を見ていきたいと思います。


## 地理的なばらつきの観察

地理的なばらつきを観察するために、mokuroku の情報をもとに、ラスタタイルやベクトルタイルを用いて、ヒートマップのような形で、データサイズを地図上に表示してみました。

![図3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/7aef5166-263d-8067-fe5c-e97fc8c22fc8.png)


概観した結果、都市部の方がタイルのファイルサイズが大きいことがわかります。東京や大阪の周辺は大きいファイルサイズのタイルが多いです。タイルはメッシュデータと考えられるため、人口密度に似たような傾向を示すと考えらます。実際、総務省統計局から出されている[平成27年度のメッシュごとの人口総数](https://www.stat.go.jp/data/mesh/h27_w.html)（≒人口密度）と似た傾向があります。
また、海の水域データは、陸地の周辺のみ作成されていることも分かります。

もう少し詳細に観察するために、タイルの X、Y 座標帯ごとに最大サイズをグラフにプロットし、地図と並べてみました。

![図2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/017bfd8f-2784-b033-feb9-4d890a1ca812.png)


結果として、最大のファイルサイズは大阪府茨木市のタイルであることがわかりました。福岡～関東付近まで、300 kb 以上となるような地帯が続いており、大阪周辺や東京周辺は特に最大サイズが大きいことがわかります。その他、石垣市や那覇市、松山市、仙台市といった都市で局所的に大きなファイルサイズのタイルが見られます。意外と札幌は、とりわけ大きいファイルサイズのタイルは存在しないようです。

大都市圏の様子を細かく見ていくと、都心よりも郊外のベッドタウンの方がタイルのファイルサイズが大きい傾向があります。経験的に、タイルのファイルサイズが大きくなる原因は建物データによるという認識がありますが、都心に比べて、郊外は小さなサイズの建物が密集しているため、その分地物数や頂点数が増えて、タイルのファイルサイズが増大したと思われます。上記で挙げた石垣市や那覇市、松山市といった都市についても、このような特徴を持つ地域が存在したため、大きなファイルサイズのタイルが生じたと考えられます。

以下は東京周辺、大阪周辺の様子（最適化ベクトルタイルを重ね合せ）です。皇居周辺や淀川周辺でファイルサイズの小さいタイルが見受けられます。

![東京.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/5261221b-48db-b1b5-e827-1c1b00aea141.png)
![近畿.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/da4f3529-2572-1999-8dc0-aa500abc7c80.png)


### 観察用 Web 地図

https://mghs15.github.io/tile-size-distribution/

## ラスタタイルとの比較

国土地理院ではラスタ形式の標準地図等のタイルも提供しており、合わせて mokuroku も提供されています。標準地図のラスタタイルについても同様に、ヒストグラムとヒートマップで傾向を観察して、ベクトルタイルのものと比較してみました。

ラスタタイルもベクトルタイル同様に、1 kb 以下の小さなファイルサイズのデータが大量にあります。1 kb 以上のタイルを見てみると、こちらも１つ山なりの分布を示しています。最も多いのは50 kb 前後のものですが、この分布はベクトルタイルの物よりも大きいです。一方、最大のファイルサイズは113,544 byte（約111 kb）であり、ベクトルタイルの最大値よりも小さいです。

![図4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/b5258998-32d5-3dcd-d8d1-1030ca0b935b.png)

データの表現等が異なるため、厳密な比較は難しいですが、ベクトルタイルの方がデータ量が小さくなりやすい傾向があると考えられます。同様の傾向について、以下のような報告事例もあります。
* [地理院地図の今後の展望～ベクトルタイル地図で世界最先端を走る～（第50回国土地理院報告会）](https://www.gsi.go.jp/common/000233923.pdf#page=10)
* [ベクターマップとMapboxエコシステム（第12回地理院地図パートナーネットワーク会議）](https://maps.gsi.go.jp/pn/meeting_partners/data/20201201/mapbox.pdf#page=13)

一方、ベクトルタイルのファイルサイズの最大値は大きくなりがちと言えそうです。

また、地理的な分布も観察してみると、ベクトルタイルのファイルサイズが小さい傾向にある山間部であっても、ラスタタイルだとそれなりのサイズがあります。ラスタタイルの場合、地物数や頂点数と関係なく、画面に占める割合や色の変化が多いとデータ量が多くなりがちです。以下は大阪周辺の例です（最適化ベクトルタイルを重ね合せ）。

![図5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/3727df23-d5c2-6429-9b2f-15a1a5a67f5b.png)


以上より、ベクトルタイルを導入することで、ラスタタイルに比べて平均的なデータ量削減が狙えますが、ファイルサイズの最大値には十分注意をしてデータ設計をする必要があると考えられます。


## データの加工・作成について

以下、視覚化のために行った作業を紹介します。

### 必要なデータの抽出

概して mokuroku のサイズは非常に大きいです。地理院地図Vector の場合、ダウンロード時は gzip 圧縮されていますが、それでも約70 MB あります。解凍すると、約192 MB でした。これでは取り扱いしにくいので、ZL15のみを抽出したいと思います。
コードは以下の通りです。これで、`_points.csv` というファイルに ZL15の必要なデータ（z、x、y とファイルサイズ）だけ抽出できます。

```JavaScript
const fs = require('fs');
const readline = require('readline');

const rs = fs.createReadStream('./mokuroku.csv');
const ws = fs.createWriteStream('./_points.csv');
const rl = readline.createInterface({input: rs, output: ws});

const tz = 15; // target zoom level

rl.on('line', (line) => {
  if(!line || line=="") return;
  const csv = line.split(",");
  const path = csv[0].split("/");
  const z = +path[0];
  if(z != tz) return;
  const x = +path[1];
  const y = +path[2].split(/\./)[0]; // remove extension
  const size = +csv[2];
  const s = `${z},${x},${y},${size}\n`;
  ws.write(s);
});
rl.on('close', () => {
  console.log("END");
});
```

この抽出した CSV を Excel に読み込んで、ヒストグラムを作成しました。設定としては、以下の通りです。

* ビンの幅を 4096 byte（4 kb）に設定
* ビンのオーバーフローを129622 byte（128 kb）に設定
* ビンのアンダーフローを1024 byte（1 kb）に設定


### ラスタタイルの作成

今回の対象は ZL15のタイルとしていますので、この ZL15のタイル1枚の領域が１つのピクセルになるようなラスタタイルを作成して、ファイルサイズに応じた色を付けて、その分布を視覚的に確認してみたいと思います。

方法としては、上記で作成した ZL15のタイル情報のみを抽出した CSV をもとに、作成するタイルごとに別々の CSV に切り分けて、その各 CSV を１対１で画像へ変換する方法としました。

画像は512（2の9乗） pixel 四方とします。ZL15に対しては、9だけ ZL を上げて、ZL6の画像タイルを作成することになります。

まずは、CSV の切り分けです。

```JavaScript
const fs = require('fs');
const readline = require('readline');

const rs = fs.createReadStream('./_points.csv');
const rl = readline.createInterface({input: rs});

const tz = 13; // target zoom level

rl.on('line', (line) => {
  if(!line || line=="") return;
  const csv = line.split(",");
  const z = +csv[0];
  const x = +csv[1];
  const y = +csv[2];
  const size = +csv[3];
  if(z != tz ) return;
  const d = 9;
  const mz = tz - d; // 2^9 = 512
  const mx = x >> d;
  const my = y >> d;
  const px = (x - (mx << d)); // = mx * Math.pow(2,d)
  const py = (y - (my << d)); // = my * Math.pow(2,d)
  const r = `${z},${x},${y},${px},${py},${size}\n`;
  fs.appendFile(`./buf/${mz}-${mx}-${my}-tile.csv`, r, (err) => {
    if(err) throw err;
  });
});
rl.on('close', () => {
  console.log("END");
});
```
:::note warn
どのようなサイズのデータが、どのような順番で来るのかわからない想定で、書き出しの際は、`fs` の `appendFile()` で都度書き出すようにしています。しかし、これは非常に遅くなるため、入力データの大きさが想定できるなら、いったんメモリに保持して、最後に一度に書き出すようにした方が効果的です。
:::


あとは、出力された CSV をタイル画像へ変換します。今回は [sharp](https://github.com/lovell/sharp) を利用しました。また、出力したタイル画像ですが、フォルダを作成するのを横着して `{z}-{x}-{y}-tile.csv.png` という形式にしています。

```JavaScript
const fs = require('fs');
const sharp = require('sharp');

const tilesize = 512;

fs.readdir("./buf/", (err, files) => {
  files.forEach( file => {
    let buf = [];
    const data = fs.readFileSync(`./buf/${file}`, 'utf8');
    const lines = data.split("\n");
    const tmp = {};
    lines.forEach( line => {
      const c = line.split(",");
      if(!tmp[`${c[4]}`]) tmp[`${c[4]}`] = {};
      tmp[`${c[4]}`][`${c[3]}`] = +c[5]/(1024);
    });
    for(let i=0; i<tilesize ; i++){
      for(let j=0; j<tilesize ; j++){
        let r=255; let g=255; let b=255; let a=255;
        if(tmp[`${i}`] && tmp[`${i}`][`${j}`]){
          let v = Math.floor(tmp[`${i}`][`${j}`]);
          g = g - v/2; b = b - v;
          if(v>255) r = r - (v - 255)/2;
        }else{
          a = 0;
        }
        
        if(r>255) r = 255; if(r<0) r = 0;
        if(g>255) g = 255; if(g<0) g = 0;
        if(b>255) b = 255; if(b<0) b = 0;
        
        buf.push(r);
        buf.push(g);
        buf.push(b);
        buf.push(a);
      }
    }
    const unit8arr = Uint8Array.from(buf);
    sharp(unit8arr,{raw: {width: tilesize , height: tilesize , channels: 4}})
    .toFile(`./docs/img/${file}.png`, (err, info) => {
      if(err){console.error(err);}
    });
  });
});
```

タイルは ZL6しか作成していないので、ZL7以上の場合は、ZL6のタイルを拡大表示（オーバーズーム）することになります。しかし、Mapbox GL JS で表示しようとした場合、オーバーズームの際に画像にアンチエイリアスがかかり、ぼやけてしまうのが課題です。

### ベクトルタイルの作成

ラスタタイルだと、オーバーズームでぼやけてしまうことが分かりましたので、別途、各タイルの領域をポリゴンとしてベクトルタイルを作成することにしました。ベクトルタイルを作成することで、Web 地図上で、インタラクティブに各タイルのサイズを確認することもできて便利です。

まずは、mokuroku を入力として、対象のZL15のタイルだけ取り出して、改行区切りの GeoJSON として書き出していきます。タイル座標から、経緯度への変換は、[Slippy map tilenames の説明ページ](https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames)に例が記載されています。

```JavaScript
const fs = require('fs');
const readline = require('readline');

const rs = fs.createReadStream('./mokuroku.csv');
const ws = fs.createWriteStream('./_polygon.ndjson');
const rl = readline.createInterface({input: rs, output: ws});

const tz = 15; // target zoom level

const tile2long = (x,z) => { return (x/Math.pow(2,z)*360-180); }
const tile2lat  = (y,z) => {
  const n=Math.PI-2*Math.PI*y/Math.pow(2,z);
  return (180/Math.PI*Math.atan(0.5*(Math.exp(n)-Math.exp(-n))));
}

rl.on('line', (line) => {
  if(!line || line=="") return;
  const csv = line.split(",");
  const path = csv[0].split("/");
  const z = +path[0];
  if(z != tz) return;
  const x = +path[1];
  const y = +path[2].split(/\./)[0];
  const size = +csv[2];
  const vx1 = tile2long(x, z);
  const vx2 = tile2long(x + 1, z);
  const vy1 = tile2lat(y, z);
  const vy2 = tile2lat(y + 1, z);
  let geojson = {
    "type": "Feature",
    "properties": {
      "z": z,
      "x": x,
      "y": y,
      "size": size
    },
    "geometry": { 
      "type": "Polygon",
      "coordinates": [[
        [vx1, vy1],
        [vx2, vy1],
        [vx2, vy2],
        [vx1, vy2],
        [vx1, vy1]
      ]]
    }
  }
  const s = JSON.stringify(geojson) + "\n";
  ws.write(s);
});
rl.on('close', () => {
  console.log("END");
});
```

作成された GeoJSON は [tippecanoe](https://github.com/felt/tippecanoe) でベクトルタイルへ変換します。今回は以下のようなオプションにしました。

```bash
tippecanoe -l pixel -e ./docs/tile/ ./_polygon.ndjson \
  --force --coalesce --reorder --hilbert \
  --no-simplification-of-shared-nodes --no-line-simplification \
  --no-tile-size-limit --no-tile-compression --no-feature-limit \ 
  --minimum-zoom=6 --maximum-zoom=6 --base-zoom=6
```

:::note warn
ラスタタイルと同じ ZL6で作成しましたが、1枚で 4 MB という非常に大きなタイルができてしまいました。今回は検証用なのでそのまま使っていますが、注意しないと、ベクトルタイルはサイズが非常に大きくなってしまうので、作成者の手腕がより問われると思います。
:::

### 地図上への表示

地図表示には、[Mapbox GL JS](https://github.com/mapbox/mapbox-gl-js/) を利用し、必要に応じて、位置の目安として[最適化ベクトルタイル](https://github.com/gsi-cyberjapan/optimal_bvmap)を重ね合わせています。

### レポジトリ

https://github.com/mghs15/tile-size-distribution

## 感想

今回、ベクトルタイルのファイルサイズの地理的分布について確認してみました。他者が提供しているベクトルタイルを観察するだけでしたが、今回見た特徴は、日本を対象としたベクトルタイルを作成する際にも参考になるかと思います。

国土地理院は、今回対象にした地理院地図Vector のベクトルタイルの他、最適化ベクトルタイルも公開しているため、その分布がどのように変わるか楽しみです。最適化ベクトルタイルの mokuroku が出たら試してみたいと思います。

また、他の ZL については、また別の傾向を示す可能性もありますので、さらに深堀りできそうです。何か新しい知見があれば、追加更新か、別途記事にしたいと思います。