特定のジオメトリが含まれるすべてのXYZタイルを取得する

## Abstract 要旨
GISのデータはラスター（画像）とベクター（点・線・ポリゴン）の2つに大別されます。これらはしばしば膨大なサイズとなり、効率的なデータ配信のためにXYZタイルという概念が存在します。XYZタイルの概念の詳細な説明は割愛しますが、でかいデータを分割して必要な箇所だけ配信する仕組みです。標準規格なのでフレームワークやアプリケーションを問わず、極地の一部を除く地球全域が常に同じタイルに区分されます。
「この地物ってXYZタイルだとどこに属するんだろう」「日本列島を網羅するXYZタイルの一覧が欲しい」というニーズが、ここ最近の私の中でありましたので、Node.jsモジュールとPythonモジュールのそれぞれのパターンについて本記事で説明します。
なおXYZタイルは、平面に落とした地球を一定のルールで分割しているため、パラメータをもとに計算する事も可能です（参考：[Qiita - 地図タイルの計算まとめ](https://qiita.com/MALORGIS/items/1a9114dd090e5b891bf7)）。

## Sample Data 題材
北海道のポリゴンをディゾルブしたGeoJSONを題材とします。
[国土地理院 - タイル座標確認ページ](https://maps.gsi.go.jp/development/tileCoordCheck.html#5/35.362/138.731)によれば、北海道が含まれるタイルは以下の画像のとおり（例としてズームレベル7）。検証のために、目視してタイルを一覧にしておきます。
<img width="400" alt="スクリーンショット 2020-04-10 21.57.24.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/73197/31980262-49df-c226-ef69-a4175fe586fa.png">

```
ズームレベル7のときに北海道が含まれるすべてのタイル
[x, y, z]
[114, 45, 7]
[116, 45, 7]
[113, 46, 7]
[114, 46, 7]
[115, 46, 7]
[116, 46, 7]
[113, 47, 7]
[114, 47, 7]
[115, 47, 7]

```
## Envs 実行環境
- macOS 10.15.4
- Node.js v12.16.1
- Python 3.8.2

## JavaScript Node.js module / tile-cover

[GitHub - Mapbox / tile-cover](https://github.com/mapbox/tile-cover)
License:MIT

### Install インストール

```shell
npm install @mapbox/tile-cover
```

### Usage 使い方

ローカルのGeoJSONファイルを読み込んで、その地物が含まれるすべてのタイルのインデックスを取得するサンプルです。前述のとおり、北海道のポリゴンをディゾルブによりひとつの地物として扱っています（tile-coverは緯度経度で計算するのでCRSはEPSG:4326です）。

```javascript:tileindex.js
//Library Import
var fs = require('fs')
var cover = require('@mapbox/tile-cover');

//Load local GeoJSON file
var geojson_str = fs.readFileSync('./dissolved_hokkaido.geojson')
var geojson_obj = JSON.parse(geojson_str);

//use tile-cover
var limits = {
    min_zoom: 7,
    max_zoom: 7
};
var tile_index = []
tile_index = cover.tiles(geojson_obj.features[0].geometry, limits)

console.log(tile_index)
```

本スクリプトをコンソールで実行、タイル一覧が出力される。

```shell
node tileindex.js
```

```shell
[
  [ 114, 45, 7 ],
  [ 116, 45, 7 ],
  [ 113, 46, 7 ],
  [ 114, 46, 7 ],
  [ 115, 46, 7 ],
  [ 116, 46, 7 ],
  [ 113, 47, 7 ],
  [ 114, 47, 7 ],
  [ 115, 47, 7 ]
]
```

ちゃんと9枚すべてがリストアップされました。

### Tips 参考

- cover.tiles()の引数は、GeoJSONのgeometryのObjectとmin_zoomとmax_zoomによるObject
- min_zoomとmax_zoomには同じ値を入れる事を推奨（少し挙動にクセがある）

## Python module / tiletanic
[GitHub - DigitalGlobe / tiletanic](https://github.com/DigitalGlobe/tiletanic)
License:MIT

### Install インストール

```shell
pip install tiletanic
```

### Usage 使い方

Node.js編と同様、ディゾルブされた北海道を題材としますが、tiletanicでは緯度経度ではなくWebメルカトルの座標を用いるため、CRSをEPSG：3857としておきます（内部的にCRSを変換する事も可能ですが手間がかかるので、データ自体をあらかじめQGISで変換しました）。

```python:tileindex.py
import json

import shapely
import tiletanic

if __name__ == "__main__":
    zoomlevel = 7

    #geojson読み込み
    geojson_dict = {}
    with open('./dissolved_hokkaido_3857.geojson', 'r') as f:
        geojson_dict = json.load(f)

    #ディゾルブにより地物数は1
    feature = geojson_dict['features'][0]['geometry']
    feature_shape = shapely.geometry.shape(feature)

    #タイルスキームの初期化
    tiler = tiletanic.tileschemes.WebMercator()

    covering_tiles_itr = tiletanic.tilecover.cover_geometry(tiler, feature_shape, zoomlevel)
    covering_tiles = []
    for tile in covering_tiles_itr:
        tile_xyz = [tile[0], tile[1], tile[2]]
        covering_tiles.append(tile_xyz)

    print(covering_tiles)
```

Pythonで実行

```shell
python tileindex.py
```


```shell
[[114, 45, 7], [113, 46, 7], [113, 47, 7], [114, 46, 7], [115, 46, 7], [114, 47, 7], [115, 47, 7], [116, 45, 7], [116, 46, 7]]
```

Node.jsのtile-coverと出力される順番は異なりますが、ちゃんと9枚のタイルのインデックスが出力されました。

### Tips 参考
- tiletanic.cover_geometryの引数は、tiletanic.tileschemesのインスタンス、shapely.geometryのインスタンス、ズームレベル値です
- tiletanicをpipでインストールする際、依存先であるshapelyも同時にインストールされます

## Sample codes サンプルコード
[GitHub - Kanahiro / tileindex-sample](https://github.com/Kanahiro/tileindex-sample)
GeoJSONファイルを含めたすべてのコードが上記のリポジトリにあります。