はじめに

こんにちは。
今回も前回の続きです。よろしければご覧ください。

前回は、図形の位置するWeb Mapのタイルの番号を得る方法の紹介でしたが、今回はそのタイルそのものをポリゴンとして作ってみようという話です。
緯度経度←→タイル番号の変換ロジックは、前回同様、宙畑とOpenStreetMapのWikiの記事を参考にしています。ありがとうございます。

これにもう１つ加工を加え、タイルをポリゴンにするArcGIS Proでの処理を加えています。なお、処理途中で座標を元にしたポリゴン作成を行うのですが、方法は以前の記事でも紹介しています。そちらもご覧ください。

使用上の注意

• 紹介するのは、特定のGISフィーチャクラス（例えば行政界など）の「envelope」に重なる、指定ズームレベルのタイルポリゴンを作成する方法です。envelopeは全フィーチャを含む最小の「長方形」なので、GISのフィーチャクラスと重ならない範囲も含みます。

• 使用データ、作成データはフィーチャクラスでもシェープファイルでもかまいませんが、ここではフィーチャクラスとしています。

• コードはArcGIS Proに含まれるJupyter Notebookの利用を想定しています。ArcGIS Proのバージョンは執筆当時最新の3.1.2です。

• データの空間座標系は緯度経度を（EPSG:4326など）想定しています。

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

• 出力データは、もし既にあれば、上書きするようにしているので注意してください。

処理の概要

1. 画面更新をオフにする（処理速度に関係）

2. 対象範囲全体を含むタイル番号を計算

3. タイルごとのポリゴンの位置と属性リストを作成

4. リストを使って新規フィーチャクラスに書き込み

5. 上記フィーチャクラスを複製し、複製元は削除

6. 画面更新をオンにする

私のPC環境では、神奈川県を対象とした場合、全処理時間は3分30秒ほどでした。4の処理が逐一処理なので時間がかかります。タイルの大きさが全部一緒だとすれば、もっと速くなりますが、逐一処理を知っておくと応用が効くのでこのようにしています。

Pythonコードの例

適当に分けて示します。
必要に応じコメントを付記しています。
一部の解説はコードのあとに書きます。

import arcpy
from arcpy import env

env.overwriteOutput = True #上書き許可

# 画面に追加しないことで処理を早くする
env.addOutputsToMap = False

temporaryFC = 'temporary' #中間ファイル（これを直接最終ファイルにしても良い）
resultFC='pbftiles' #最終ファイル

env.addOutputsToMap = False　処理の本質と関係ないのですが、新規に出来たデータをマップにレイヤとして追加するのを抑止します。デフォルトではマップに追加してしまいます。マップを表示していなくても、データの変化に応じて、常にレイヤ表示を書き換えようとするので、処理時間が大変遅くなってしまうことがあり、それを抑止します。
その代わり、処理が終わったら自分で追加しないといけません。
コードの最後ではTrueに戻しておきます。

temporaryFC = 'temporary' #中間ファイル（これを直接最終ファイルにしても良い）
resultFC='pbftiles' #最終ファイル

boundary="神奈川県行政界" # タイル座標の計算範囲とするフィーチャクラス
desc=arcpy.Describe(boundary) # フィーチャクラスの概要を取得
upperLeftLat=desc.extent.YMax # 全体の北西緯度
upperLeftLon=desc.extent.XMin # 全体の北西経度
lowerRightLat=desc.extent.YMin # 全体の南東経度
lowerRightLon=desc.extent.XMax # 全体の南東緯度

中間ファイルのフィーチャクラスは、直接最終ファイルでも良いのですが、データのファイルサイズ表示で更新されない場合があるので、複製します。
boundaryは求めるタイルの範囲を決めるための参照データなので何でも構いませんが、ここでは神奈川県の行政界ポリゴンとしています。コンサベーションGISコンソーシアムジャパンのページでも公開しています。
arcpy.Describeで参照フィーチャクラスのデータ概要を得ます。
タイル番号は北西から南東にかけ、xtiile（東→西）、ytile（北→南）と数字が増えるので、最大の経度と緯度の範囲をここで決めます。
これは、全フィーチャを含む最小の長方形の緯度と経度の情報であり、実際の行政界の外の範囲も多く含まれることに注意が必要です。ここではそれを取り除く対応はしていません。

# 緯度lat、経度lon、ズームレベルzoomから、タイルX、Y座標を計算 xtile 、ytile
# 宙畑の記事を参考　https://sorabatake.jp/en/8282/
def get_tile_num(lat, lon, zoom):
    lat_rad = math.radians(lat)
    n = 2.0 ** zoom
    xtile = int((lon + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return (xtile, ytile)

# zoomレベルとタイル座標からタイルの緯度経度
# 左下（南西）：経度、緯度 右上（北東）：経度、緯度のタプルを取得する
def get_tile_bbox(z, x, y):
    def num2deg(xtile, ytile, zoom): #タイルとズームレベルから座標の経度、緯度を取る
        n = 2.0 ** zoom
        lon_deg = xtile / n * 360.0 - 180.0
        lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
        lat_deg = math.degrees(lat_rad)
        return (lon_deg, lat_deg)
    
    right_top = num2deg(x + 1, y, z) # 右上：経度、緯度
    left_bottom = num2deg(x, y + 1, z) # 左下：経度、緯度
    return (left_bottom[0], left_bottom[1], right_top[0], right_top[1])

ここは肝になる部分です。宙畑さんで記載されているコードをそのまま引用しています。ありがとうございます。

• 最初の関数get_tile_numは緯度、経度、ズームレベルからタイル座標を計算するものです。既存フィーチャクラスのタイル範囲を決めるのに使います。

• 次の関数は２つありますが、内側にあるnum2degが、そのタイルとズームレベルから座標の経度と緯度を得る関数、外側のget_tile_bboxがポリゴンの左下（南西：経度、緯度と右上（北東）：経度、緯度の座標を得る関数です。左上、右下だったり左下、右上だったりが混在しているので注意してください。何処から始めても違いはないです。

# 行政界の範囲の左上（北西）、右下（南東）のタイル番号を得る
(upperLeftXtile,upperLeftYtile)=get_tile_num(upperLeftLat,upperLeftLon,zoom) # 左上タイルX、Y座標
(lowerRightXtile,lowerRightYtile)=get_tile_num(lowerRightLat,lowerRightLon,zoom) # 右下タイルX、Y座標

features=[] # 図形と属性リストを複数フィーチャ分保持しておくリスト

# タイルごとに緯度経度の座標セットと属性を配列に書き込んていく
for xtile in range(upperLeftXtile, lowerRightXtile+1, 1):
    for ytile in range(upperLeftYtile, lowerRightYtile+1, 1):
        attributes=[] # 図形オブジェクトとタイル番号を持つ属性リスト（１フィーチャ）
        lonlats=get_tile_bbox(zoom,xtile,ytile)
        attributes=[[(lonlats[0],lonlats[1]),(lonlats[2],lonlats[1]),(lonlats[2],lonlats[3]),(lonlats[0],lonlats[3])],zoom,xtile,ytile]
        features.append(attributes)

• タイルのポリゴンを得たいデータの左上と右下のタイル番号を得ます。

• タイルごとに緯度、経度の座標を計算します。get_tile_bbox関数で返ってくる値は4隅のうち右上、左下の2点だけなので、4点になるように組み合わせを変更します。lonlatsの添字が0から3の4つありますが、タイルのポリゴンで言うと、左下→右下→右上→左上の座標になります。その他、その場所のズームレベル、xtile番号、ytile番号もattributesにセットされます。ポリゴンの場合5つ目の頂点を1つ目と同じに配置して閉じるという理屈になっているフォーマットもありますが、ここでは4点で閉じてくれます。ズームレベル、xtile番号、ytile番号はなくてもいいです。前回の記事で初回した方法で得ることが出来ます。ただ処理上、その属性を保持していて、もったいないのでここで書いてしまいます。

• これらを1セットとして、全タイルについてfeaturesにセットします。

# 配列データを格納する空のフィーチャクラスを作る
arcpy.management.CreateFeatureclass(env.workspace,temporaryFC, 'POLYGON', spatial_reference=4326)
# フィーチャクラスの属性フィールド名称とタイプを追加する
arcpy.AddField_management(temporaryFC,'zoom','LONG')
arcpy.AddField_management(temporaryFC,'xtile','LONG')
arcpy.AddField_management(temporaryFC,'ytile','LONG')

# 1行ずつ、属性リストを元にして図形と属性を追記していく
# SHAPE@が図形そのものを表す情報で、以降がテーブルの属性
for row in features:
    with arcpy.da.InsertCursor(temporaryFC,['SHAPE@','zoom','xtile','ytile']) as cursor:
        cursor.insertRow(row) 

arcpy.management.CopyFeatures(temporaryFC,resultFC) # 中間ファイルを最終ファイルにコピー
arcpy.management.Delete(temporaryFC) # 中間ファイルを削除
arcpy.env.addOutputsToMap = True # 画面追加を再開

• featuresはリストそのものなので、これを入れるための箱（フィーチャクラス）をarcpy.management.CreateFeatureclassで作ります。

• ズームレベル、xtile、ytileの値を入れるフィールドをarcpy.AddField_managementで作ります。何れも整数値をとりますがズームレベルが16の時tileは5桁になるので、それにあわせてlong型にしています。

• arcpy.da.InsertCursorで1行ずつポリゴンを属性とともに追加していきます。['SHAPE@','zoom','xtile','ytile']は追加する属性のリストですが、SHAPE@には図形の情報（ここではポリゴンの4隅の座標）が入ることになります。

• 基本的には、ここで終わりでよいのですが、ここで終わりにした場合、フィーチャクラスのサイズがデータ追加前のままになっていることがあります。データとしては間違っていないのですが、気持ちが悪いので複製arcpy.management.CopyFeaturesして更新しておきます。

• 複製したら元のデータは要らないので削除arcpy.management.Deleteします。

• 最後にnv.addOutputsToMap = Trueとしてマップへの追記をできるように戻しておきます。

ここまでの処理の結果を参照データと比較して示します。
約2万個のポリゴンが出来ています。県と重ならない場所もありますが、その場所が不要であれば空間検索機能を使って絞り込むことで対応できます。
この処理にかかった時間はおよそ3分30秒でした。
ほぼarcpy.da.InsertCursorのところにかかった時間が全てです。参照範囲が大きくなればもっと時間がかかるものと思われます。

前回例示した江ノ島付近で見てみるとこんな感じです。xtileとytileをラベルとして表示していますが、前回示した結果やタイル座標確認ページと同じになっています。

まとめ

今回はWeb Mapのタイルポリゴンをpythonで作成する例を紹介しました。タイルを計算する方法は勿論ですが、属性を逐一追加するarcpy.da.InsertCursorや図形を追加するSHAPE@は応用の効く方法なので、覚えておくと良いと思います。

最後に

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

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

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