実行環境
Windows10
Python 3.9.5
pyproj 3.5.0
pyshp 2.3.1
はじめに
地図タイルの位置や大きさは、ズームレベル、X 座標、Y 座標で決まりますが、一般的にそれらの組み合わせ（ZXY）だけでは GIS ソフトで位置情報として扱えません。
というのも、多くの GIS ソフトでは地図タイルを読み込むための本来のやり方（XYZ 形式の URL の登録、WMTS サーバへ接続など）ができるので、わざわざタイル座標（ZXY）を意識する必要がないからです。
したがって、もし GIS ソフト上でタイル画像（256×256 ピクセルの PNG ファイル）を単独で載せなくてはいけない場合は、その画像に位置情報を与えなくてはいけません。
手動でジオリファレンスしてもよいですが、もちろん精度と手間の問題があります。
この記事は、GIS ソフトユーザが遭遇し得る上記のようなケースを想定して書いたものです（n 番煎じではありますが）。
具体的には、地図タイル座標（ZXY）をWeb メルカトル座標系の座標値へ変換する処理を通して、GIS ソフトに載せられる以下の２種類のデータを作成する方法について取り上げます。
- 地図タイルの図郭ポリゴン
- タイル画像のワールドファイル
なお、python サードパーティライブラリの pyproj と pyshp を使っています。
地図タイルと Web メルカトル座標系
本記事では地図タイルと表現していますが、XYZ タイル、マップタイル、ラスタータイル、Web タイルなど呼び方がいろいろあるようです。
ここでは正式名称や言葉の定義については掘り下げず、正方形の画像を並べてモニタ上に地図を表示する仕組みを「地図タイル」、そこで使われている座標系を「Web メルカトル座標系」という程度の説明に留めます。
イメージをつかみやすいのは、やはり国土地理院のページでしょうか。
ワールドファイル
ラスタを GIS ソフトに載せるときに位置情報を付与するための付属ファイル（サイドカーファイル）です。
その実態は、ピクセル座標をマップの座標に変換する変換式（アファイン変換）の６つの係数を記したテキストファイルです。
ファイルの拡張子を .pgw（PNG ファイルの場合）にして、ラスタと同じファイル名で同じフォルダに保存して使用します。
QGISで読めますが、GIS ソフトによっては対応していないかもしれません。
ワールドファイルの仕様については、下記等をご確認ください。
なお、ほとんどの GIS ソフトや Web 地図にワールドファイル付き画像を出力する機能が備わっているので、特定の範囲の地図タイルを図郭の区別なくまとめて扱う場合は、それらの機能をご利用ください。
関数定義
今回４つの関数を定義します。
num2deg()
wiki からそのまま引用しています。
タイルの ZXY 座標を渡して、WGS84 の緯度経度を受け取ります。
座標はタイルの北西端の座標になります。
XY 座標それぞれに +1 すると、元のタイルの南東端の座標を取得できます。
num2m()
タイルの ZXY 座標を渡して、タイルの四隅座標を Web メルカトル座標系で受け取ります。
単位はmです。
get_wkt()
epsg コードを渡して、wkt （=well known text）を受け取ります。
シェープファイル付属の座標系定義ファイル（prj ファイル）を作成するときに使用します。
get_worldfile_params()
タイルの ZXY 座標を渡して、ワールドファイルに記載する６つの係数を受け取ります。
import math
import pyproj
import shapefile
def num2deg(xtile, ytile, zoom):
"""
Slippy map tilenames
https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
"""
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 (lat_deg, lon_deg)
def num2m(xtile: int, ytile: int, zoom: int):
"""
タイルの四隅座標をWebメルカトル座標系で取得
EPGSコード
WGS84 4326
Webメルカトル 3857
Returns
--------
x_m_w : float
西端x座標 [m]
x_m_e : float
東端x座標 [m]
y_m_n : float
北端y座標 [m]
y_m_s : float
南端y座標 [m]
"""
transformer = pyproj.Transformer.from_crs(4326, 3857, always_xy=True)
# 四隅座標（WGS84: 緯度経度）
lat_deg_n, lon_deg_w = num2deg(xtile, ytile, zoom)
lat_deg_s, lon_deg_e = num2deg(xtile + 1, ytile + 1, zoom)
# 四隅座標（Webメルカトル: m）
x_m_w, y_m_n = transformer.transform(lon_deg_w, lat_deg_n)
x_m_e, y_m_s = transformer.transform(lon_deg_e, lat_deg_s)
return (x_m_w, x_m_e, y_m_n, y_m_s)
def get_wkt(epsg: int):
"""
座標系のwktを取得
Returns
--------
wkt : str
well known text
"""
crs = pyproj.CRS.from_epsg(epsg)
wkt = crs.to_wkt(version='WKT1_ESRI')
return wkt
def get_worldfile_params(xtile: int, ytile: int, zoom: int):
"""
ワールドファイルのパラメータを取得
Returns
--------
dx : float
x方向解像度
lx : float
x方向回転パラメータ
ly : float
y方向回転パラメータ
dy : float
y方向解像度
px_m_w : float
北西端ピクセル中心x座標
px_m_n : float
北西端ピクセル中心y座標
"""
x_m_w, x_m_e, y_m_n, y_m_s = num2m(xtile, ytile, zoom)
# パラメータ計算
dx = (x_m_e - x_m_w) / 256
lx = 0.0
ly = 0.0
dy = (y_m_s - y_m_n) / 256
px_m_w = x_m_w + dx/2
px_m_n = y_m_n + dy/2
return (dx, lx, ly, dy, px_m_w, px_m_n)
地図タイルの図郭ポリゴンを作成
たとえば、千代田区付近のタイルをズームレベル14で作成します。
タイル座標は地理院地図等であらかじめ調べます。
以下のコードを実行しています。
zoom = 14 # ズームレベル
xtile_w = 14551 # 西端xタイル座標
xtile_e = 14553 # 東端xタイル座標
ytile_n = 6450 # 北端yタイル座標
ytile_s = 6452 # 南端yタイル座標
shp = 'tile.shp'
prj = 'tile.prj'
# シェープファイル作成
with shapefile.Writer(shp, shapetype=5) as w: # POLYGON は shapetype=5
# フィールド追加
w.field('zoom', 'N')
w.field('xtile', 'N')
w.field('ytile', 'N')
for xtile in range(xtile_w, xtile_e + 1):
for ytile in range(ytile_n, ytile_s + 1):
# 四隅座標の取得
x_m_w, x_m_e, y_m_n, y_m_s = num2m(xtile, ytile, zoom)
# ジオメトリ追加
w.poly([[
[x_m_w, y_m_n],
[x_m_e, y_m_n],
[x_m_e, y_m_s],
[x_m_w, y_m_s],
[x_m_w, y_m_n]
]])
# 書き込み
w.record(zoom, xtile, ytile)
# 座標系定義ファイル(prj ファイル)作成
with open(prj, 'w') as f:
wkt = get_wkt(3857) # Webメルカトル
f.write(wkt)
出力されるシェープファイルを QGIS で表示しています。
地理院地図で表示している図郭が、ポリゴンになっていることが確認できました。
タイル画像のワールドファイルを作成
皇居付近のタイル画像のワールドファイルを作成します。
タイル画像は、以下の URL からダウンロードします。
https://cyberjapandata.gsi.go.jp/xyz/std/14/14552/6451.png
実行コードは以下の通りです。
zoom = 14 # ズームレベル
xtile = 14552 # x座標
ytile = 6451 # y座標
pgw = 'tile.pgw'
# パラメータ取得
params = get_worldfile_params(xtile, ytile, zoom)
# 書き込み
with open(pgw, 'w') as f:
f.write('\n'.join([str(p) for p in params]))
以下の内容のテキストファイルが出力されます。
9.554628535639495
0.0
0.0
-9.55462853564677
15556468.77391334
4258454.942509471
ワールドファイル（tile.pgw）をタイル画像（tile.png）と同じフォルダに入れて、先ほどの QGIS のマップ上にダウンロードしてきたタイル画像を載せています。
濃い色の部分がダウンロードしてきたタイル画像です。
問題なさそうな精度（目で見る限り）で位置情報を付与できました。
日本の反対側でも確認
心配性なので、南半球西半球でも確認してみます。
日本国外の地図になるので、地理院タイルではなくオープンストリートマップ（OSM）を使っています。
オープンストリートマップの配信 URL
https://tile.openstreetmap.org/{z}/{x}/{y}.png
地理院地図で外部タイルとして OSM を読み込み、ブラジルのサンパウロあたりのどこかしら（いったいどこだろう）を表示しています。
このあたりで同じように図郭ポリゴンとワールドファイルを作成し、QGIS に載せてみます。
大丈夫みたいです。
おわりに
今や地図タイルを扱うツールやライブラリはいろいろあるとは思いますが、GIS ソフトで地図タイルを扱う場面などで少しでも参考になれば幸いです。
今回はタイル画像１枚を単独で処理していますが、実用上は複数枚のタイルを一度に扱うこともあると思いますので、ダウンロードの際は常識的なリクエストを心がけましょう。
TIPs
地理院タイルの PNG 画像は色数によってビット深度が異なります。
例えば海だけのタイル画像は1ビットです。
そのため、ビット深度がまちまちだとGIS ソフト等でタイルをマージ（モザイク）できないことがあります。
PNG 画像はすべて24ビットにしてから処理しましょう。
また、本記事では "マップ" の投影法は Web メルカトル（EPSG: 3857）を指定しているので、タイル画像そのものには投影定義をしていません。
マップに別の投影法を指定するときは、タイル画像にも投影定義（あるいは変換）が必要になるのでご注意ください。