宙畑 Sorabatake

Tellusのアップデート

干渉解析ツールTelluSARが高解像度SAR画像「ASNARO-2」に対応!実際に解析をしてみた(コード付き)

Tellusの干渉解析ツール「TelluSAR(てるーさー)」で、高解像度SAR画像である「ASNARO-2(あすなろつー)」が解析できるようになったので、実際に使ってみました。

2021年7月26日、Tellusの干渉解析ツール「TelluSAR(てるーさー)」で、高解像度SAR画像である「ASNARO-2(あすなろつー)」が解析できるようになりました。

新しく対応した画像ではどんなものが見られるのか、さっそく宙畑編集部でも触ってみることにしました!

1.干渉解析ツールTelluSAR(てるーさー)とは?

Tellusで公開している衛星データ解析ツール

TelluSARとは、Tellusが無料で公開している衛星データの解析ツールです。どなたでもTellusのアカウントを作成すれば、ブラウザから簡単なマウス操作で、もしくはAPIを利用することで解析を行うことができます。

Tellusでは、衛星データの利用ポリシーから、衛星データを直接取得するAPIはTellusの開発環境からのみアクセスすることができますが、TelluSARをはじめとする解析ツールは直接衛星データを取得するものではないため、APIはTellusの開発環境以外からも利用することができます。

干渉解析とは?

TelluSARで扱うのはSARデータと呼ばれる衛星データです。

SARデータというのは、衛星に搭載されるセンサの一つで、衛星から電波を発射して地表面で跳ね返って来る電波を観測することで、地表面の様子を知ることができるセンサです。

電波は文字通り波なので、2時期のSARデータを比べると波の微妙なズレを検知することができます(上図)。これが、地盤の沈下や隆起が分かる仕組みとなります。

詳しく知りたいという方は、宙畑の別の記事で解説していますので、ぜひご覧ください。

2.高解像度SAR「ASNARO-2」とは?

ASNARO-2とは?

これまでTelluSARでは、JAXAの衛星データ「PALSAR-2(ぱるさーつー)」に対応していましたが、この度新たな衛星データ「ASNARO-2(あすなろつー)」に対応しました。

ASNARO-2はは経済産業省が推進する「小型化等による先進的宇宙システムの研究開発」プロジェクトの2号機で、NECが開発・運用を行っている衛星です。

ASNARO-2とPALSAR-2の違い

PALSAR-2とASNARO-2は、観測の際に用いる電波の周波数が異なり、それにより、得意とする対象物が異なってきます。

PALSAR-2ではL-bandと呼ばれる比較的長い波長で観測を行うため、植物など2時期で状態が異なる場合、葉や枝を透過した地盤の様子(隆起や沈下)などを把握することに向いています。

一方のASNARO-2はX-bandと呼ばれる短い波長で観測を行っており、自然が多いエリアでの2時期比較では、成長した草木などに反応してしまうため、あまり向いていません。都市部など人工物を見ることに適していると言えるでしょう。

また、干渉SARは2時期の微妙なズレを検知する解析手法ですが、そのためには衛星自身の位置も2時期でほぼ同じにする必要があります。

衛星が同じ軌道を維持するためには、多くの燃料が必要であり、大型衛星であるPALSAR-2では非常に高い精度で軌道が維持されていますが、燃料を多く積めない小型衛星であるASNARO-2ではPALSAR-2と比較して精度が落ちる点にも注意が必要です。

3.実際に解析してみる(ブラウザ)

では、実際にTelluSARを利用して解析を行ってみましょう。
まずはブラウザでマウス操作でできる方法を説明します。

Tellusでアカウントを作成します(無料)。

TelluSARの商品ページで[購入]ボタンをクリックします。
※本商品は無料のため、[購入]ボタンを押しても課金されることはありません。

③右上のユーザー名をクリックして、出てくるプルダウンメニューから「Tellus OS」を選択します。

④Tellus OSに画面が切り替わったら、右上にあるマイツールから「TelluSAR」を選択します。

⑤出てくるポップアップの中で、データが選べるようになっているので「ASNARO-2」を選択します。

今回は事前にTellusの衛星データチームにお伺いして、以下の2セットがおすすめということで、こちらで進めて行きたいと思います。(前述の通り、ASNARO-2が干渉するペアを見つけるのはかなり難しく、なかなかうまくいかなくても粘り強くトライしてくださいね!)

・武蔵小杉
AS200534028728-190103
AS200555328728-190117
・川崎
AS201166107290-200223
AS201187407290-200308

⑥地図を武蔵小杉や川崎あたりに合わせてから[表示中の位置からbeforeを検索]ボタンを押下します。

画像の候補が出てくるので、上で上げたファイル名に一致するものを選び、[生成]を押下します。

通常解析には数10分かかりますが、すでに他ユーザーが解析を実行しているペアであれば、すぐに解析結果が掲示されます。

⑦出来上がった解析結果をクリックして、コヒーレンスという項目にチェックを入れます。

白い部分がうまく干渉している部分、黒い部分はうまく干渉が行えていない部分です。
上図では右上の方が白くなっているので、そのあたりはうまく行えているということになりますね。

⑧続いて、差分干渉の横のボックスにチェックを入れて、解析結果を表示させます。

さきほどみたコヒーレンス画像で黒くなっている部分は、ノイズになっていますが、白かった部分は虹色の縞が見えています。

出てきた結果の解釈ですが、今回のケースでは特に沈下・隆起している場所は見受けられないということになります。

縞の解釈については、詳しくは以下のドキュメントも参照してください。
https://tellusar.tellusxdp.com/document/user_guide.pdf

4.実際に解析してみる(API)

ここまでは、ブラウザ上でマウス操作により解析を実施する方法をご案内してきました。同じ操作をAPIを使っても行うことができます。プログラミングのできる方はご自身の解析環境からもぜひチャレンジしてみてください。

TelluSARの使い方は以前の記事で説明していますが、ASNARO-2対応に合わせてAPIをアップデートしたため、適宜読み替える必要があります。

コードの意味は上の記事を参照していただければと思います。
以下に、今回のアップデートに対応したサンプルコードを貼っておきます。

マーケットトークンの取得(ver1)

マーケットで購入したAPIの利用にはマーケットトークンを取得する必要があります。
以下にマーケットトークン取得(ver1.0)のサンプルコードを提示します。

※ver1でのマーケットトークンではプロバイダIDとツールラベルを利用します。
provider_id = ‘tellus-product’
tool_label = ‘tellusar-api-f’

import requests, datetime, json

TOKEN = 'ここにトークンを貼ります'
provider_id = 'tellus-product'
tool_label = 'tellusar-api-f'

# マーケットトークン有効期限(30分に設定)
# デフォルトでは5分、最長で60分まで設定可能
expires_at = (datetime.datetime.now(datetime.timezone(datetime.timedelta(hours=+9))) + datetime.timedelta(minutes=30)).isoformat()
def get_market_token(payload={}):
    url = 'https://sdk.tellusxdp.com/api/manager/v1/auth/api_access_token/token'
    headers = {
        'Authorization': 'Bearer ' + TOKEN,
        'Content-type': 'application/json'
    }
    r = requests.post(url, headers=headers, data=json.dumps(payload))
    if r.status_code is not 200:
        raise ValueError('status error({}).'.format(r.status_code))
    return json.loads(r.content) 

# マーケットトークンを発行する
ret = get_market_token({'provider_id': provider_id, 'tool_label': tool_label, 'expires_at':expires_at})
print(ret['token'])

衛星データの検索

# 利用できるシーン一覧を取得する

def get_free_scene(market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v2/asnaro2/search'
    headers = {
        'Authorization': 'Bearer ' + market_token
    }
    r = requests.get(url, headers=headers, params=payload)
    if r.status_code is not 200:
        print(r.status_code)
    return json.loads(r.content)

get_free_scene_result = get_free_scene(ret['token'])

print(get_free_scene_result['data']['scenes'])

干渉ペアの検索

# 武蔵小杉(2019年1月3日)のシーンID
scene_id = 'AS200534028728-190103'

# ペアの候補となるシーン一覧を取得する
# after_scene_id

def get_scene_after(scene_id, market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v2/asnaro2/search/{}/afters'.format(scene_id)
    headers = {
        'Authorization': 'Bearer ' + market_token
    }
    r = requests.get(url, headers=headers, params=payload)
    if r.status_code is not 200:
        print(r.status_code)
    return json.loads(r.content)

get_scene_after_result = get_scene_after(scene_id, ret['token'], {})
print(get_scene_after_result)

干渉SAR解析を実行する

# 処理を実行する
# work_idが返却される

def request_work(market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v2/works'
    headers = {
        'Authorization': 'Bearer ' + market_token
    }
    r = requests.post(url, headers=headers, data=json.dumps(payload))
    print(r.url)
    if r.status_code is not 200:
        print(r.status_code)
    return json.loads(r.content)

#nlool_rg、nlook_az、filterは固定値

work_result = request_work(ret['token'], {
    'satellite': 'asnaro2',
    'before_scene_id': 'AS200534028728-190103',
    'after_scene_id': 'AS200555328728-190117',
    'polarisation': 'HH',
    'nlook_rg': 5,
    'nlook_az': 7,
    'filter': 0,
})
print(work_result)

処理結果の一覧を取得

# 処理結果の一覧を取得する

def get_work_list(market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v2/works'
    headers = {
        'Authorization': 'Bearer ' + market_token
    }
    r = requests.get(url, headers=headers, params=payload)
    print(r.url)
    if r.status_code is not 200:
        print(r.status_code)
    return json.loads(r.content) 

get_work_result_list = get_work_list(ret['token'], {})
print(get_work_result_list)

処理結果のダウンロード

from skimage import io
from io import BytesIO
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image


def get_finge_diff(work_id, z, x, y, market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v2/works/{}/pngs/fringe_diffs/{}/{}/{}.png'.format(work_id, z, x, y)
    headers = {
        'Authorization': 'Bearer ' + market_token
    }
    r = requests.get(url, headers=headers, params=payload)
    print(r.url)
    if r.status_code is not 200:
        print(r.status_code)
    return io.imread(BytesIO(r.content))


# 対象の座標を指定(Tellus OSで確認)
z = 11
x = 1818
y = 807

# 差分干渉画像を取得して表示
work_id = work_result['data']['work_id']
file = get_finge_diff(work_id, z, x, y, ret['token'], {})
io.imsave('sarsample.png',file)

結果の表示

# 結果の表示
fig = plt.figure(figsize=(10,10))
im=Image.open("./sarsample.png","r")
plt.imshow(np.array(im))
武蔵小杉周辺の干渉SAR解析結果(2019年1月3日と1月17日)
credit: Original data provided by NEC

5.まとめ

今回はアップデートして新しいデータ「ASNARO-2」に対応した衛星データ解析ツールTelluSARの使い方をご紹介しました。

途中でご説明した通り、なかなかペアとなる画像を見つけるのは難しい部分もありますが、ぜひ実際に触っていただきながら、良いペアを見つけて見てください。