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

Abstract 要旨

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

Sample Data 題材

北海道のポリゴンをディゾルブしたGeoJSONを題材とします。
国土地理院 - タイル座標確認ページによれば、北海道が含まれるタイルは以下の画像のとおり（例としてズームレベル7）。検証のために、目視してタイルを一覧にしておきます。

ズームレベル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
License:MIT

Install インストール

npm install @mapbox/tile-cover

Usage 使い方

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

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)

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

node tileindex.js

[
  [ 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
License:MIT

Install インストール

pip install tiletanic

Usage 使い方

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

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で実行

python tileindex.py

[[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
GeoJSONファイルを含めたすべてのコードが上記のリポジトリにあります。