【ゼロからのTellusの使い方】大きなタイル画像を取得するには？地図タイルを扱う上で知っておきたいTips4選
これまで「ゼロからのTellus」シリーズでは、Jupyter Lab上でデータを取り扱う方法を紹介してきました。シリーズの中でたびたび地図タイルという言葉が登場しますが、本記事ではそんな地図タイルをTellusのJupyter Labで扱う際に役立つテクニックをまとめて紹介します。
2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html
これまで「ゼロからのTellus」シリーズでは、衛星データプラットフォームTellusのJupyter Lab上でデータを取り扱う方法を紹介してきました。
シリーズの中でたびたび地図タイルという言葉が登場しますが、本記事ではそんな地図タイルをTellusのJupyter Labで扱う際に役立つテクニックをまとめて紹介します。
画像を取得したい地点の座標を簡単に求めたり、1枚の大きな画像を取得するにはどんなコードを書いたらいいか、サンプルを動かしながら覚えていきましょう。
TellusでJupyter Lab(Jupyter Notebook)を使う方法は、「Tellusの開発環境でAPI引っ張ってみた」をご覧ください。
なお、地図タイルはTellus OSからも取得できます。
「データ選択」末尾の「グリッド」から「タイル座標」を選択すると、下記のようにタイル座標が表示されます。
取得したいタイル座標をクリックしたのち、「タイル座標」右に記載されているコピーアイコンを押すことで値をコピーすることも可能です。
1.地図タイルとは
地図タイルとはWebメルカトル地図を256×256ピクセルの正方形に切り分けたものです。各タイル画像にはズーム率（z）と位置（x, y）が決められており、世界地図（ズーム率0）を元に、ズーム率が1つ上がるごとに4分割して作られます。
Webメルカトル地図は、社会や地理の教科書などでお馴染みのメルカトル図法で投影された世界地図を、東経180°から西経180°、北緯85.0511°から南緯85.0511°の範囲で切り取ったもので、ネットの地図サービスでよく使われています。
| ズーム率 | タイル枚数 | 座標サンプル タイル座標(x, y) タイルの左上の経度緯度 |
| 0 | １ | (0, 0) => (-180, 85.0511) |
| 1 | 4 | (0, 1) => (-180, 0) (1, 0) => (0, 85.0511) (1, 1) => (0, 0) |
| 2 | 16 | (0, 1) => (-180, 66.5133) (0, 2) => (-180, 0) (0, 3) => (-180, -66.5133) (2, 1) => (0, 66.5133) |
2.緯度経度からタイル座標を求める
Tellusは主に地図タイル形式で衛星データを提供しています。
例：Landsat-8の光学画像を取得するAPI
https://gisapi.tellusxdp.com/api/v1/landsat8/{z}/{x}/{y}.png
APIのパスにはタイル座標（z, x, y）を指定します。
ですが、ある地点の画像を取得したいのにその地点の緯度経度はわかるけどタイル座標がわからない、なんてことがよくあります。
地図タイル確認のツールなどを使えば調べることはできますが、毎回調べるのは少々手間です。
そこでPythonを使って経度緯度からタイル座標を計算するサンプルコードを紹介します。
import math
def get_tile_num(lat, lon, zoom):
"""
緯度経度からタイル座標を取得する
Parameters
----------
lat : number
タイル座標を取得したい地点の緯度(deg)
lon : number
タイル座標を取得したい地点の経度(deg)
zoom : int
タイルのズーム率
Returns
-------
xtile : int
タイルのX座標
ytile : int
タイルのY座標
"""
# https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
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)※コードの一部にopenstreetmapのサンプルをお借りしています。
このメソッドは緯度経度（lon, lat）とズーム率（zoom）を引数に与えると、その点を含んだ地図タイルの座標を返します。
実際にこのメソッドを使って、ズーム率12で富士山山頂（緯度：35.3606°、経度：138.7273°）を含むタイル座標(x, y) = (3626, 1617)を求めている例です。
3.タイル座標から緯度経度求める
また逆に、タイル画像を使ってデータ処理するとき画像に含まれる範囲の緯度経度が知りたくなることがあります。
これもPythonを使ってタイル座標から緯度経度を計算するサンプルコードを紹介します。
import math
def get_tile_bbox(z, x, y):
"""
タイル座標からバウンディングボックスを取得する
https://tools.ietf.org/html/rfc7946#section-5
Parameters
----------
z : int
タイルのズーム率
x : int
タイルのX座標
y : int
タイルのY座標
Returns
-------
bbox: tuple of number
タイルのバウンディングボックス
(左下経度, 左下緯度, 右上経度, 右上緯度)
"""
def num2deg(xtile, ytile, zoom):
# https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
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])
※コードの一部にopenstreetmapのサンプルをお借りしています。
このメソッドはタイル座標（z, x, y）を引数に与えると、その画像の左下と右上の緯度経度を返します。
返り値はGeoJSONなどで使われるバウンディングボックスの形式 （左下経度、左下緯度、右上経度、右上緯度の順）に従います。
バウンディングボックスがわかると、左上/右下の緯度経度も知ることができます。
実際にこのメソッドを使って、富士山山頂を含むタイル画像（z=12, x=3626, y=1617）のバウンディングボックス（左下：138.6914°, 35.3174°、右上：138.7793°, 35.3890°）を求めている例です。
4.タイル画像の各ピクセルの緯度経度を求める
タイル画像の四隅の緯度経度はget_tile_bboxメソッドを使って求めることができました。
では、画像の各ピクセルの緯度経度は求めるにはどうすればいいでしょうか？
簡易的にであれば、四隅の緯度経度の差を256で等分すれば求めることができます。
import numpy as np
def get_tile_approximate_lonlats(z, x, y):
"""
タイルの各ピクセルの左上隅の経度緯度を取得する（簡易版）
Parameters
----------
z : int
タイルのズーム率
x : int
タイルのX座標
y : int
タイルのY座標
Returns
-------
lonlats: ndarray
タイルの各ピクセルの経度緯度
256*256*2のnumpy配列
経度、緯度の順
"""
bbox = get_tile_bbox(z, x, y)
width = abs(bbox[2] - bbox[0])
height = abs(bbox[3] - bbox[1])
width_per_px = width/256
height_per_px = height/256
lonlats = np.zeros((256,256,2))
lat = bbox[3]
for i in range(256):
lon = bbox[0]
for j in range(256):
lonlats[i,j,:] = [lon, lat]
lon += width_per_px
lat -= height_per_px
return lonlats
※コード内でget_tile_bboxメソッドを利用しています。
このメソッドはタイル座標（z, x, y）を引数に与えると、その画像の各ピクセルの左上の緯度経度を計算し256*256*2の配列として返します。実際にこのメソッドを使って、富士山山頂を含むタイル画像（z=12, x=3626, y=1617）の各ピクセルの座標を求めると以下のようになります。
lonlats[0, 0] = [138.69140625 35.38904997] からわかるように、[経度, 緯度]の順で返ることに注意してください。
しかしこれはあくまでも近似値です。地図タイルはメルカトル図法という投影法を使っているため、同じタイル内でも場所によって緯度の変化率がわずかに異なります（北半球なら緯度が高くなるほど1ピクセルあたりの緯度の変化が小さくなる）。そのため256で等分すると実際の座標からはずれてしまいます。
変化率の差は微々たるものなので、ちょっとした計算であればこれで問題ありませんが、解析によっては厳密に座標を求めたいケースもあるでしょう。その場合は次のサンプルコードを使ってみてください。
import numpy as np
def get_tile_lonlats(zoom, xtile, ytile):
"""
タイルの各ピクセルの左上隅の経度緯度を取得する
Parameters
----------
z : int
タイルのズーム率
x : int
タイルのX座標
y : int
タイルのY座標
Returns
-------
lonlats: ndarray
タイルの各ピクセルの経度緯度
256*256*2のnumpy配列
経度、緯度の順
"""
x = int(xtile * 2**8)
y = int(ytile * 2**8)
z = zoom + 8
lonlats = np.zeros((256,256,2))
for i in range(256):
for j in range(256):
bbox = get_tile_bbox(z, x + j, y + i)
lonlats[i,j,:] = [bbox[0], bbox[3]]
return lonlats
※コード内でget_tile_bboxメソッドを利用しています。
このメソッドはタイル座標（z, x, y）を引数に与えると、メルカトル図法においてより高精度な緯度経度の値を256*256*2の配列として返します。
5.タイルを連結して大きな画像を作る
地図タイルの画像は1枚あたり256×256ピクセルとサイズが決められています。
しかし解析をするとき、画像に収まる範囲が狭くて困るケースがよくあります。
ズームアウトをすれば1枚の中に広い範囲を収めることはできますが、解像度が下がってしまい見たいものがよくわからなくなっては本末転倒です。
ここでは複数のタイル画像をつなぎ合わせて、解像度を落とさず広い範囲を収めた画像を作るサンプルを紹介します。
import numpy as np
def get_combined_image(get_image_func, z, topleft_x, topleft_y, size_x=1, size_y=1, option={}, params={}):
"""
結合されたタイル画像を取得する
Parameters
----------
get_image_func : function
タイル画像取得メソッド
引数はz, x, y, option, params
z : int
タイルのズーム率
topleft_x : int
左上のタイルのX座標
topleft_y : int
左上のタイルのY座標
size_x : int
タイルを経度方向につなぎわせる枚数
size_y : int
タイルを緯度方向につなぎわせる枚数
option : dict
APIパスのオプション(z,x,y除く)
params : dict
クエリパラメータ
Returns
-------
combined_image: ndarray
結合されたタイル画像
"""
rows = []
blank = np.zeros((256, 256, 4), dtype=np.uint8)
for y in range(size_y):
row = []
for x in range(size_x):
try:
img = get_image_func(z, topleft_x + x, topleft_y + y, option, params)
except Exception as e:
img = blank
row.append(img)
rows.append(np.hstack(row))
return np.vstack(rows)
このメソッドにAPIで画像を取得するメソッド（get_image_func）と、起点となるタイル座標（z, topleft_x, topleft_y）、結合するタイル枚数（size_x, size_y）、APIに渡すパラメータ（option, params）を与えると、結合されたタイル画像を取得することができます。
APIで画像を取得するメソッドは、APIごとに個別で作る必要があります。
例えばAPIでASNARO-1の画像を取得するメソッドは次のように作ります。
import requests
from skimage import io
from io import BytesIO
%matplotlib inline
TOKEN = 'ここには自分のアカウントのトークンを貼り付ける'
def get_asnaro1_image(z, x, y, option, params={}):
"""
ASNARO-1のタイル画像を取得する
Parameters
----------
z : int
タイルのズーム率
x : int
タイルのX座標
y : int
タイルのY座標
option : dict
APIパスのオプション(entityId)
params : dict
クエリパラメータ
Returns
-------
img: ndarray
タイル画像
"""
url = 'https://gisapi.tellusxdp.com/ASNARO-1/{}/{}/{}/{}.png'.format(option['entityId'], z, x, y)
headers = {
'Authorization': 'Bearer ' + TOKEN
}
r = requests.get(url, headers=headers, params=params)
if not r.status_code == requests.codes.ok:
r.raise_for_status()
return io.imread(BytesIO(r.content))
このメソッドを使って256×256ピクセルよりも大きい画像を取得してみましょう。
option = {
'entityId': '20181226061336632_AS1'
}
z = 14
x = 13979
y = 6340
combined = get_combined_image(get_asnaro1_image, z, x, y, 3, 2, option)
io.imshow(combined)
サンプルではentityId = 20181226061336632_AS1のシーンから、z = 14, x = 13979, y = 6340のタイルを起点に、横に3枚、縦に2枚のタイル画像をつなぎ合わせた、768 * 512ピクセルの画像を取得しています。
ASNARO-1以外の画像も同様に取得できます。
以下は、APIでPalsar-2の画像を取得するメソッドです。
def get_palsar2_image(z, x, y, option, params={}):
"""
Palsar-2のタイル画像を取得する
Parameters
----------
z : int
タイルのズーム率
x : int
タイルのX座標
y : int
タイルのY座標
option : dict
APIパスのオプション(entityId)
params : dict
クエリパラメータ
Returns
-------
img: ndarray
タイル画像
"""
url = 'https://gisapi.tellusxdp.com/PALSAR-2/{}/{}/{}/{}.png'.format(option['entityId'], z, x, y)
headers = {
'Authorization': 'Bearer ' + TOKEN
}
r = requests.get(url, headers=headers, params=params)
if not r.status_code == requests.codes.ok:
r.raise_for_status()
return io.imread(BytesIO(r.content))
このメソッドを使って今度は縦長の衛星画像を取得してみましょう。
option = {
'entityId': 'ALOS2017871020'
}
z = 12
x = 2055
y = 1367
io.imshow(get_palsar2_image(z, x, y, option))
combined = get_combined_image(get_palsar2_image, z, x, y, 2, 3, option)
io.imshow(combined)
サンプルではentityId = ALOS2017871020のシーンから、z = 12, x = 2055, y = 1367のタイルを起点に、横に2枚、縦に3枚のタイル画像をつなぎ合わせた、512 * 768ピクセルの画像を取得しています。
以上、TellusのJupyter Labで地図タイルを扱う際に役立つテクニックをまとめて紹介しました。
今回紹介したのはライブラリへの依存を最小限にした初歩的なテクニックばかりです。
タイル画像を扱っていく上でこんなことができたらいいのにといったアイデアや、こうするともっと便利だぜといった実例があれば、ぜひ宙畑編集部までお寄せください。
