こんにちは。小ネタです。。。
実用性はほぼないですが、Web Mapのタイルに「親しむ」には参考になると思います。

球面メルカトル図法

Web Mapの空間座標系として数多く使われるのは球面メルカトル図法：WGS 84/Pseudo-Mercator（EPSG:3857）ですが、そこで使われるタイル座標を手持ちのGISデータの範囲で計算しようという話です。
この図法は、GoogleMapや、地理院地図のベクトルタイルなどで採用されています。ズームレベルの階層ごとに地域を小分けにしたタイルを作り、そこに地図データが入っています。
表示するときは必要なタイルだけ読むようにすれば、処理が速くなるというというわけです。
それ以上の細かい説明は是非検索してみてください。

（追記）今回参考にしたのは、以下のページです。
OpenStreetMapのWikiの方には、数式や他言語での計算方法も掲載されていますので、一度読んでみてください。

さて、「まれに」そのタイル番号が気になることがあるわけです。
私が普段使うのは、地理院の提供しているサービス「タイル座標確認ページ」で、地図を表示すれば、そのズームレベルにあったタイル座標が表示されます。

とはいえ、複数同時に確認したり、手持ちのGISデータと確認するには少し手間がかかります。そこで、直接手持ちのGISデータのタイル番号を図形ごと（フィーチャごと）に計算してしまおうというものです。

使用上の注意

• 紹介するのは、タイルを求めたい図形の位置を、どこか1つの点（ポリゴンなら重心）で代表させた場合の計算方法です。１つの図形から得た点がが複数のタイルにまたがっている場合は、この方法だけでは正しく求まらず省略されてしまいますので、複数の点を個別に処理するコードを書く必要があります。

• ArcGIS Proと書いてますが、ArcMapでも同じように出来ます。ただし、ここで示したpythonコードは、ArcGISProはPyhon3、ArcMapはPython2です。

• データの空間座標系は緯度経度を（EPSG:4326など）想定しています。もちろん、それでなくても投影変換の手間を入れれば出来ますが、説明が込み入るので省略します。

• pythonで必要なインデント幅は揃っていればいいだけですが、ここでは４つにしています。

• 「フィールド演算」機能での使用を想定しています。そのため元のデータを上書きしてしまう場合があります。新規にフィールドを作るなどして対応してください。

タイル座標の計算

フィーチャのテーブルを表示し、データを書き込みたいフィールドを選んで、フィールド演算を使える編集状態にします。その上でコードを入力します。

xtileの値を得る

コードブロックにこのように入力します。

def getxtile(longitude):
    zoom=16
    n = 2.0 ** zoom
    xtile = int((longitude + 180.0) / 360.0 * n)
    return xtile

• getxtitle：自作関数名で何でも良いです。引数は経度です。

• zoom：ズームレベルの数値です。ここでは16にしています。タイル座標確認ページでは、2～18の数値を指定できますが、小さいズームレベルになると日本より広い範囲になってしまいますので、12～16が現実的かなと思います。

呼び出し側には関数名と引数として経度を示す式を入れます。

getxtile(!SHAPE.CENTROID.X!)

例えばポリゴンフィーチャでは、そのフィーチャの重心のX座標を示す式（!SHAPE.CENTROID.X!）を入れます。！マークを前後に入れるのを忘れずに。ポイントデータの座標もこれで取得できますが、ラインデータなどのときは、CENTROIDの代わりに、firstPoint、lastPointなども使えます。

実際のフィールド演算のGUIでは以下のように入力します。入力テーブルは計算したいフィーチャ名です。フィールド名は任意で構いませんが、フィールド型は整数型（Long）以上にしておきましょう。

ytileの値を得る

今度はytileを計算してみましょう。
コードブロックはこのようになります。

def getytile(latitude):
    zoom=16
    lat_rad = math.radians(latitude)
    n = 2.0 ** zoom
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return ytile

• getytitle：自作関数名

• zoom：ズームレベルの数値

経度より複雑な式になっていることがわかります。
一方の呼び出し側には関数名と引数として緯度を示す式を入れます。

getytile(!SHAPE.labelPoint.Y!)

ポリゴンデータの場合は、重心のY座標を示す式（!SHAPE.CENTROID.Y!）を入れます。実際のフィールド演算のGUIでは以下のように入力します。

作ったタイル座標を確認

計算結果をタイル座標確認ページの同じ場所で見比べてみます。神奈川県の江ノ島付近を例として示しています。
この場所では、58000台の数字がxtile、25000台の数字がytileです。
値も同じになっているのがわかります。

ここで、緑の枠はズームレベル16の１タイルの範囲を示したポリゴンです。ですので、結果が同じなのは当たり前と思われるかもしれませんが、実はポリゴン作成とタイル座標計算は別々です。（計算速度の問題で別々にしていたのですが、別々に計算してもトータル時間で一緒になることがわかったので訂正）。そうなるとタイルのポリゴンをどう作るかという話ですが、今回は省略します。次回以降で紹介します。

式はどのようにも使えるので、zoom=14とした結果の比較を示します。zoom=16のタイルに対して計算しているので、2の自乗、4*4の範囲が同じ値になります。図では太線で示した範囲が同じになっています。

まとめ

今回はWeb Mapのタイルを１つずつ計算する方法を紹介しました。
次回以降で、タイルのポリゴンをpythonで作成する例を紹介しようと思います。

最後に

ここまでご覧いただき、ありがとうございました。

普段は北海道に本拠地を置くNPOに所属し、環境保全を主な題材としてGISやリモセンに関する仕事をしています。
コンサベーションGISコンソーシアムジャパン の活動もその1つです。
この分野の仕事や活動でなにかお困りのある方、ご相談ごとのある方など、是非コメントをお寄せいいただくか、上記WEBサイト掲載のメールアドレスまでメールをお寄せください。

コンサベーションGISコンソーシアムジャパンの活動及びこのnoteでの活動はボランティアで行っているのですが、何分資金が不足しています。
もしサポートしてもいいなという方がいらっしゃいましたら、よろしくお願いします。