宙畑 Sorabatake

衛星データ

【ゼロからのTellusの使い方】Jupyter Notebookを使ってASNARO-1の光学画像を取得する

車の数が数えられるほどの、50cm分解能人工衛星「ASNARO-1」の光学画像を、Tellus上のJupyter Notebookを使って取得する方法をご紹介します。

Tellusは、さくらインターネット株式会社が経済産業省から受託した、政府衛星データのオープン&フリー化及びデータ利用環境整備事業における衛星データプラットフォームのことで、登録するとブラウザから様々な衛星データを利用することができます。

本記事ではTellus上のJupyter Notebookを使ってASNARO-1の光学画像を取得する方法を紹介します。

(1) Jupyter Notebookを使おう

Jupyter Notebookとは、初心者でも手軽にデータ解析や機械学習ができると人気なプログラミング言語Pythonの開発環境を提供するツールで、これを使えばブラウザから手軽にPythonを実行することができます。

Jupyter NotebookでブラウザからPythonを実行、結果を確認できる。

TellusでJupyter Notebookを使うには、Tellusの登録後にJupyter Notebookの利用申請が必要となります。
利用申請の方法はこちらの記事を参考にしてください。
Jupyter NotebookでTellusを使ってみた〜雪質解析(1)解析準備編~

申請が通るとTellus OSのツールバー上にJupyter Notebookのアイコンが現れるので、ここからJupyter Notebookにログインします。

Jupyter Notebookのアイコンをクリック
画面の左上にユーザー名が表示されるのでクリック
あなたのJupyter Notebookにログイン

(2) APIでASNARO-1の光学画像を取得するには

TellusのAPIによりASNARO-1の光学画像を取得することができます。
APIの使い方はリファレンスにまとめられています。
APIリファレンス

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

 

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

ASNARO-1の光学画像を取得するには、ズーム率と座標、またシーンIDを指定する必要があります。
シーンIDとは、ASNARO-1の観測毎につけられた識別番号のことです。

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

https://gisapi.tellusxdp.com/api/v1/asnaro1/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)

(3) コードを書いてみよう

それでは実際に画像を取得するコードを書いてみましょう。
コードを書くには、workフォルダを開き、右上の「New」をクリック、プルダウンから「Python3」を選択します。

workフォルダで「New」をクリックし「Python3」を選択
※画像ではworkフォルダにファイルが並んでいますが、初めてJupyter Notebookを使う人はここに何も表示されません。

開いた画面に次のコードを貼り付けます。
Pythonはわからないという人もコピペだけで動かすことができるので、「Jupyter Notebookを使えばこんなこともできるんだ」という参考にしてください。

import os, json, requests, math
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
%matplotlib inline

TOKEN = "ここには自分のアカウントのトークンを貼り付ける"

def get_ASNARO_scene(min_lat, min_lon, max_lat, max_lon):
    url = "https://gisapi.tellusxdp.com/api/v1/asnaro1/scene" \
        + "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
    
    headers = {
        "content-type": "application/json",
        "Authorization": "Bearer " + TOKEN
    }
    
    r = requests.get(url, headers=headers)
    return r.json()

scene = get_ASNARO_scene(20.425278, 122.933611, 45.557222, 153.986389)

貼り付けたら、

TOKEN = "ここには自分のアカウントのトークンを貼り付ける"

この部分を自分のトークン(APIキー)で上書きします。
トークンはマイページのAPIアクセス設定(要ログイン)から取得が可能です。

トークンの準備ができたら、「Run」をクリックしコードを実行します。
Tellusで利用できるASNARO-1のシーン情報を取得することができます。

取得できたシーンがいくつあるか確認しましょう。

print(len(scene))

このコードを貼り付けて再び「Run」をクリックしてください。

執筆時点では264シーンあることがわかりました。
次にシーン情報の中身がどんなものか確認します。

print(scene[-1])

実行すると、その時点で最新の観測シーンに関する情報が表示されます。

{
	"acquisitionDate": "Wed, 21 Nov 2018 12:36:44 GMT",
	"clat": 37.7632401320653,
	"clon": 140.44256144238,
	"cloudCover": "N/A",
	"entityId": "20190115064552292_AS1",
	"max_lat": 37.8169958287351,
	"max_lon": 140.369939374427,
	"min_lat": 37.8163117766714,
	"min_lon": 140.370845314871,
	"path": 0,
	"productId": "20190115064552292_AS1_D01_L1B",
	"row": 0,
	"thumbs_url": "https://tile.tellusxdp.com/ASNARO1thums/20190115064552292_AS1_thums.PNG"
}

このシーンは2018年11月21日に行われた観測のものです(acquisitionDate)。
場所は福島県福島市のあたりであることが画像中心の緯度経度(clat, clon)からわかります。
thumbs_urlをブラウザで開くと、サイズは小さいですがおおよその写り具合が確認できます。

Original data provided by NEC

そしてentityIdが光学画像を取得する上で必要なシーンIDです。

続いてズーム率と座標を決定します。

この座標は地図タイルと呼ばれるもので、世界地図を256*256サイズのタイルに切り出す際の各ズーム率における座標が定義されています。

各地点のタイル座標は国土地理院のタイル座標確認ページで確認することができます。

Credit : 国土地理院 タイル座標確認ページ(https://maps.gsi.go.jp/development/tileCoordCheck.html)

今回は利用するシーン(scene_id:20190115064552292_AS1)の観測範囲内のタイル座標を指定する必要があります。

ズーム率と点(緯度経度)を指定すると、その点を含んでいるタイル座標を求める数式があるのでこれを使いましょう。

def get_tail_num(lat_deg, lon_deg, zoom):
    # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 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 = 13
(xtile, ytile) = get_tail_num(scene[-1]['clat'], scene[-1]['clon'], zoom)

print(xtile, ytile)

ズーム率を13とした場合、今回のシーンの中心座標(37.7632401320653, 140.44256144238)を含むタイル座標は、(7291, 3166)であることがわかりました。

あとはシーンIDと求めたタイル座標を使って画像を呼び出すだけです。

def get_ASNARO_image(scene_id, zoom, xtile, ytile):
    url = " https://gisapi.tellusxdp.com/ASNARO-1/{}/{}/{}/{}.png".format(scene_id, zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

img = get_ASNARO_image(scene[-1]['entityId'], zoom, xtile, ytile)

io.imshow(img)

無事、ASNARO-1による光学画像を取得することができました。

以上が、TellusのJupyter Notebookを使ってASNARO-1の光学画像を取得する方法でした。
ズーム率をさらにあげると、建物の形まではっきりわかる解像度の画像を取得することもできます。

ここまでズームすると車の数が数えられます Credit : Original data provided by NEC

TellusではASNARO-1以外の衛星画像も公開されているので、ぜひ利用してみてください。