SLATSの地図タイル画像を時系列で連続表示する

はじめに

本チュートリアルでは「つばめ」(SLATS)よる高解像度画像を取得し、複数枚を時系列でつなぎ合わせる方法を紹介します。

import os, json, requests, math
import numpy as np
import dateutil.parser
from datetime import datetime
from datetime import timezone
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
%matplotlib inline
TOKEN = "ここには自分のアカウントのトークンを貼り付ける"
def get_tsubame_scene(min_lat, min_lon, max_lat, max_lon):
	url = "https://gisapi.tellusxdp.com/api/v1/tsubame/scene" \
    	+ "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
	headers = {
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return r.json()
def filter_by_date(scenes, start_date=None, end_date=None):
	if(start_date == None):
    	start_date = datetime(1900,1,1, tzinfo=timezone.utc)
	if(end_date == None):
    	end_date = datetime.now(timezone.utc)
	return [scene for scene in scenes if start_date <= dateutil.parser.parse(scene["acquisitionDate"]) and dateutil.parser.parse(scene["acquisitionDate"]) < end_date]
tsubame_scenes = filter_by_date(get_tsubame_scene(20.425278, 122.933611, 45.557222, 153.986389),
                            	datetime(2019, 4, 10, tzinfo=timezone.utc),
                            	datetime(2019, 5, 11, tzinfo=timezone.utc))

取得出来たシーンを2019年4月10日から5月10日にかけて観測されたものに絞り込んでみたところ、23シーン取得することができました。

print(tsubame_scenes[-1])
io.imshow(io.imread(tsubame_scenes[-1]["thumbs_url"]))

devSlats_cont_20200220_5.png

執筆時点で最新のシーン情報とシーンのサムネイル Credit : Original data provided by JAXA

「つばめ」は4月2日からの約1ヶ月間、「赤坂迎賓館」周辺と「小石川」周辺の2箇所を毎日定点観測していました。シーンを「赤坂迎賓館」と「小石川」で分類しておくと日ごとの変化を見る上で便利です。

observation_areas = [
	{
    	"name":"赤坂迎賓館",
    	"lat": 35.68031,
    	"lon": 139.729216
	},
	{
    	"name":"小石川",
    	"lat": 35.719582,
    	"lon": 139.744686
	}
]
def classify(areas, scenes):
	result = [[] for i in range(len(areas))]
	for scene in scenes:
    	dists = [np.linalg.norm(np.array([scene["clat"], scene["clon"]]) - np.array([area["lat"], area["lon"]])) for area in areas]
    	result[np.argmin(dists)].append(scene)
	for i in range(len(result)):
    	result[i] = sorted(result[i], key=lambda scene: dateutil.parser.parse(scene["acquisitionDate"]))
	return result
classified_scenes = classify(observation_areas, tsubame_scenes)
akasaka_scenes = classified_scenes[0]
koishigawa_scenes = classified_scenes[1]
print(len(akasaka_scenes))
print(len(koishigawa_scenes))

この関数はシーンの中心の緯度経度(clat, clon)と観測エリアの代表的な座標との距離を計算し、より近い観測エリアにシーンを分類します。執筆時点では「赤坂迎賓館」周辺が12シーン、「小石川」周辺が11シーンありました。 分類されたシーンは観測日でソート(並び替え)されます。

for scene in akasaka_scenes:
	print(scene["acquisitionDate"])

devSlats_cont_20200220_1.png

「つばめ」の観測期間中、神宮球場の外野席に大きな文字を複数日掲げてメッセージを作るコラボ企画をJAXA(宇宙航空研究開発機構)と東京ヤクルトスワローズが実施していました。このメッセージをJupyter Notebook上で確認してみます。 文字が掲げられていたのは神宮球場のライトスタンドで、地図タイルの座標はz=18, x=232811, y=103232となります。

devSlats_cont_20200220_2.png

Credit : 国土地理院 タイル座標確認ページ

(https://maps.gsi.go.jp/development/tileCoordCheck.html)

def get_tsubame_image(scene_id, zoom, xtile, ytile):
	url = " https://gisapi.tellusxdp.com/tsubame/{}/{}/{}/{}.png".format(scene_id, zoom, xtile, ytile)
	headers = {
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return io.imread(BytesIO(r.content))
io.imshow(get_tsubame_image(akasaka_scenes[0]["entityId"], 18, 232811, 103232))


devSlats_cont_20200220_3.png

Credit : Original data provided by JAXA

しかし、この座標の画像を取得したところ、画像の中に神宮球場の外野席らしきものが写っていません。「つばめ」は技術実証衛星であり、地球観測を主目的に据えた衛星ではないことから、様々な要因から撮影場所がずれてしまっていることがあります。 そのため座標を指定しても、本来よりもずれた位置の画像となってしまいます。今回のケースではz=18, x=232811, y=103230と北へずらしてみたところ、神宮球場の外野席を含んだ画像を取得することができました。

io.imshow(get_tsubame_image(akasaka_scenes[0]["entityId"], 18, 232811, 103230))


devSlats_cont_20200220_7.png

Credit : Original data provided by JAXA

シーンによってもずれる程度が異なるため、今回は指定した座標周辺の画像を複数取得し、繋げて一つの画像とする関数も作ります。

def get_tsubame_series_image(scene_id, zoom, topleft_x, topleft_y, size_x=1, size_y=1):
	img = []
	for y in range(size_y):
    	row = []
    	for x in range(size_x):
        	row.append(get_tsubame_image(scene_id, zoom, topleft_x + x, topleft_y + y))
    	img.append(np.hstack(row))
	return  np.vstack(img)
io.imshow(get_tsubame_series_image(akasaka_scenes[0]["entityId"], 18, 232811, 103230, 2, 2))


devSlats_cont_20200220_4.png

Credit : Original data provided by JAXA

これだけ広い範囲を取得すれば、多少ずれても神宮球場のライトスタンドを画像内に収めることができそうです。APIで取得できる画像は1枚あたり縦横256ピクセルですが、この画像は4枚の画像(2×2)をつなげているため縦横512ピクセルとなっています。では、先程取得した12個の「赤坂迎賓館」周辺シーンからそれぞれ神宮球場周辺を抜き出して並べてみましょう。

def make_grid_image(images, col):
	img = []
	for i in range(math.ceil(len(images)/col)):
    	row = []
    	for j in range(col):
        	index = i*col + j
        	if index < len(images):
            	row.append(images[index])
        	else:
            	row.append(np.ones_like(imgs[0])*255)
    	img.append(np.hstack(row))
	return np.vstack(img)
imgs = [get_tsubame_series_image(scene["entityId"], 18, 232811, 103230, 2, 2) for scene in akasaka_scenes]
plt.rcParams["figure.figsize"] = (12, 12) #表示サイズを拡大
io.imshow(make_grid_image(imgs, 4))


devSlats_cont_20200220_8.png

Credit : Original data provided by JAXA

この中から外野席に文字が掲げられている画像を探します。小さくて見つけにくければ、

plt.rcParams['figure.figsize'] = (12, 12)

の部分をより大きい数字にしてみて下さい。7番目、8番目、そして11番目の画像に文字が写っているのがわかるでしょうか?この3枚の画像だけ抜き出して並べてみましょう。


devSlats_cont_20200220_6.png

Credit : Original data provided by JAXA

ぼやけてしまっている画像もありますが、「祝つばめ50」というメッセージを、読み取ることができました(特に、真ん中の画像、「つば」の字が見やすいです)。 以上が、TellusのJupyter Notebookを使って「つばめ」(SLATS)よる高解像度画像を取得する方法でした。神宮球場以外にも日ごとに変化している地点は沢山あります。移りゆく東京の様子を観察してみましょう。 今回使用したスクリプトになります。

import os, json, requests, math
import numpy as np
import dateutil.parser
from datetime import datetime
from datetime import timezone
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
TOKEN = "ここにトークンを貼る"
def get_tsubame_scene(min_lat, min_lon, max_lat, max_lon):
	url = "https://gisapi.tellusxdp.com/api/v1/tsubame/scene"     	+ "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
	headers = {
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return r.json()
def filter_by_date(scenes, start_date=None, end_date=None):
	if(start_date == None):
    	start_date = datetime(1900,1,1, tzinfo=timezone.utc)
	if(end_date == None):
    	end_date = datetime.now(timezone.utc)
	return [scene for scene in scenes if start_date <= dateutil.parser.parse(scene["acquisitionDate"]) and dateutil.parser.parse(scene["acquisitionDate"]) < end_date]
tsubame_scenes = filter_by_date(get_tsubame_scene(20.425278, 122.933611, 45.557222, 153.986389),
                            	datetime(2019, 4, 10, tzinfo=timezone.utc),
                            	datetime(2019, 5, 11, tzinfo=timezone.utc))
print(len(tsubame_scenes))
type(tsubame_scenes)
print(tsubame_scenes[-1])
io.imshow(io.imread(tsubame_scenes[-1]["thumbs_url"]))
observation_areas = [
	{
    	"name":"赤坂迎賓館",
    	"lat": 35.68031,
    	"lon": 139.729216
	},
	{
    	"name":"小石川",
    	"lat": 35.719582,
    	"lon": 139.744686
	}
]
def classify(areas, scenes):
	result = [[] for i in range(len(areas))]
	for scene in scenes:
    	dists = [np.linalg.norm(np.array([scene["clat"], scene["clon"]]) - np.array([area["lat"], area["lon"]])) for area in areas]
    	result[np.argmin(dists)].append(scene)
	for i in range(len(result)):
    	result[i] = sorted(result[i], key=lambda scene: dateutil.parser.parse(scene["acquisitionDate"]))
	return result
classified_scenes = classify(observation_areas, tsubame_scenes)
akasaka_scenes = classified_scenes[0]
koishigawa_scenes = classified_scenes[1]
print(len(akasaka_scenes))
print(len(koishigawa_scenes))
for scene in akasaka_scenes:
	print(scene["acquisitionDate"])
def get_tsubame_image(scene_id, zoom, xtile, ytile):
	url = " https://gisapi.tellusxdp.com/tsubame/{}/{}/{}/{}.png".format(scene_id, zoom, xtile, ytile)
	headers = {
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return io.imread(BytesIO(r.content))
io.imshow(get_tsubame_image(akasaka_scenes[0]["entityId"], 18, 232811, 103232))
io.imshow(get_tsubame_image(akasaka_scenes[0]["entityId"], 18, 232811, 103230))
def get_tsubame_series_image(scene_id, zoom, topleft_x, topleft_y, size_x=1, size_y=1):
	img = []
	for y in range(size_y):
    	row = []
    	for x in range(size_x):
        	row.append(get_tsubame_image(scene_id, zoom, topleft_x + x, topleft_y + y))
    	img.append(np.hstack(row))
	return  np.vstack(img)
io.imshow(get_tsubame_series_image(akasaka_scenes[0]["entityId"], 18, 232811, 103230, 2, 2))
def make_grid_image(images, col):
	img = []
	for i in range(math.ceil(len(images)/col)):
    	row = []
    	for j in range(col):
        	index = i*col + j
        	if index < len(images):
            	row.append(images[index])
        	else:
            	row.append(np.ones_like(imgs[0])*255)
    	img.append(np.hstack(row))
	return np.vstack(img)
imgs = [get_tsubame_series_image(scene["entityId"], 18, 232811, 103230, 2, 2) for scene in akasaka_scenes]
plt.rcParams["figure.figsize"] = (12, 12) # 表示サイズを拡大
io.imshow(make_grid_image(imgs, 4))
imgs = [get_tsubame_series_image(scene["entityId"], 18, 232811, 103230, 2, 2) for scene in akasaka_scenes]
plt.rcParams["figure.figsize"] = (12, 12) # 表示サイズを拡大
io.imshow(make_grid_image(imgs, 4))
# 一つ一つマニュアルで指定してもよい
imgs = [get_tsubame_series_image(akasaka_scenes[0]["entityId"], 18, 232811, 103230, 2, 2),
    	get_tsubame_series_image(akasaka_scenes[1]["entityId"], 18, 232811, 103230, 2, 2),
    	get_tsubame_series_image(akasaka_scenes[2]["entityId"], 18, 232811, 103230, 2, 2)]
plt.rcParams["figure.figsize"] = (12, 12) # 表示サイズを拡大
io.imshow(make_grid_image(imgs, 3))
# indexの作成(意図的に抜き取る部分を指定する
my_indexes = [6,7,10]
akasaka_scenes_ext = []
akasaka_scenes_ext = [akasaka_scenes[i] for i in my_indexes]
print(akasaka_scenes_ext)
imgs = [get_tsubame_series_image(scene["entityId"], 18, 232811, 103230, 2, 2) for scene in akasaka_scenes_ext]
plt.rcParams["figure.figsize"] = (12, 12) #表示サイズを拡大
io.imshow(make_grid_image(imgs, 3))