宙畑 Sorabatake

衛星データ

【ゼロからのTellusの使い方】Jupyter Notebookで「つばめ」(SLATS)による高解像度画像を取得する

TellusのJupyter Notebookを使って「つばめ」(SLATS)よる高解像度画像を取得する方法を紹介します。

本記事ではTellusのJupyter Lab(Jupyter Notebook)を使って「つばめ」(SLATS)よる高解像度画像を取得する方法を紹介します。

TellusでJupyter Labを使う方法は、「Tellusの開発環境でAPI引っ張ってみた」をご覧ください。

(1)「つばめ」で日ごとに変わる東京の様子を見てみよう

「つばめ」(SLATS)はJAXAの技術試験衛星で、解像度1m以下の高解像度モノクロ画像をTellusで見ることができます。

この衛星の最大の特徴は、2019年の4月からおよそ1ヶ月間、毎日決まった時間に東京の決まった地点を観測し続けたこと。
1ヶ月間、毎日都心を撮影した「つばめ」の衛星画像をTellusで公開!

Tellusではそんな「つばめ」の観測データを56地点(267シーン)分利用することが可能です。

※画像を利用する際は、データ詳細のデータポリシーを確認して、規約を違反しないように注意してください。

https://gisapi.tellusxdp.com/tsubame/{scene_id}/{z}/{x}/{y}.png
scene_id 必須 シーンID
z 必須 Zoom Level 9-18
x 必須 タイルのX座標
y 必須 タイルのY座標

APIリファレンス

APIで「つばめ」の画像を取得するには、ズーム率とタイル座標、その他パラメータを設定する必要があります。
シーンIDとは、観測毎につけられた識別番号のことです。

Tellusが保持しているシーンのIDは検索APIで調べることができます。

https://gisapi.tellusxdp.com/api/v1/tsubame/scene?min_lat={min_lat}&min_lon={min_lon}&max_lat={max_lat}&max_lon={max_lon}
min_lat 必須 最小緯度(-90~90)
min_lon 必須 最小経度(-180~180)
max_lat 必須 最大緯度(-90~90)
max_lon 必須 最大経度(-180~180)

APIリファレンス

さっそくJupyter Lab上で呼び出してみましょう。

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()

tsubame_scenes = get_tsubame_scene(20.425278, 122.933611, 45.557222, 153.986389)

トークンはマイページのAPIアクセス設定(要ログイン)から取得してください。

再生の三角マークをクリックしコードを実行すると、東京の市ヶ谷近郊の「つばめ」のシーンが取得できます。

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(tsubame_scenes,
                                datetime(2019, 4, 10, tzinfo=timezone.utc),
                                datetime(2019, 5, 11, tzinfo=timezone.utc))
print(len(tsubame_scenes))

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

print(tsubame_scenes[-1])
io.imshow(io.imread(tsubame_scenes[-1]["thumbs_url"]))
シーン情報とサムネイル
Credit : Original data provided by JAXA

(2)観測エリアごとにシーンを分類する

「つばめ」は2019年4月からの約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"])

(3)神宮球場に掲げられたメッセージを見つけよう

「つばめ」の観測期間中、神宮球場の外野席に大きな文字を複数日掲げてメッセージを作るコラボ企画をJAXAと東京ヤクルトスワローズが実施していました。

このメッセージをJupyter Lab上で確認してみます。

文字が掲げられていたのは神宮球場のライトスタンドで、地図タイルの座標はz=18, x=232811, y=103232となります。

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))
Credit : Original data provided by JAXA

しかし、この座標の画像を取得したところ、画像の中に神宮球場の外野席らしきものが写っていません。

衛星画像の多くは、被写体に対して斜めから撮影された画像を、真上から撮影したかのように画像の加工補正を行うのですが、「つばめ」は技術実証衛星のため、この加工補正における歪みが通常よりも大きく、実際の地点との差が生まれてしまいます。

そのため座標を指定しても、本来よりもずれた位置の画像となってしまいます。
今回のケースではz=18, x=232811, y=103230と北へずらしてみたところ、神宮球場の外野席を含んだ画像を取得することができました。

io.imshow(get_tsubame_image(akasaka_scenes[0]['entityId'], 18, 232811, 103230))
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))
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))
取得した画像 Credit : Original data provided by JAXA

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

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

の部分をより大きい数字にしてみて下さい。

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

Credit : Original data provided by JAXA

ぼやけてしまっている画像もありますが、「祝つばめ50」というメッセージを、読み取ることができました(特に、真ん中の画像、「つば」の字が見やすいです)。

以上が、TellusのJupyter Labを使って「つばめ」(SLATS)よる高解像度画像を取得する方法でした。

神宮球場以外にも日ごとに変化している地点はたくさんあります。移りゆく東京の様子を観察してみましょう。

2019/11/14 : SLATSのデータ追加に伴いサンプルコードと記事の一部に修正を行いました。