日本初のオープン&フリーな衛星データプラットフォーム「Tellus」

　世界的にオープンフリーな衛星画像プラットフォームには，欧州のSentinel衛星シリーズの取得データのプラットフォームであるコペルニクスハブが有名です．
　日本にも同様に，国内の人工衛星データを中心としたオープン＆フリーなプラットフォームである「Tellus」が２０１9年2月に公開されました．
　以下は，Tellusのサイトより説明を抜粋します．

「Tellus（テルース）」は、日本初のオープン&フリーな衛星データプラットフォームです。
今まで扱いの難しかった衛星データを、誰もが簡単に、無償で扱うことができる環境を提供します。

Tellusには，衛星画像の閲覧やその組み合わせなどの操作を支援する「Tellus OS」が用意されています．
試してみるだけならばこれでもよいのですが，衛星画像のデータ処理などを行うには，Tellusが用意しているJupyter Labによる「総合開発環境」　を用いることで，pythonを用いてより様々な解析が可能となります．

Tellusのアカウント作成は，Tellusのホームページをご参考ください．
アカウント登録後に，総合開発環境への申込みができます．申込後は，問題なければ数日後に登録完了のお知らせがメールにて届きます．

Tellus 総合開発環境

　Tellusの総合開発環境は，Tellus OSよりアクセスします．

上部の赤枠をクリックすると，「開発環境を開く」のウィンドウがポップアップしますので，それをクリックしJupyter Labを立ち上げます．

まずはサンプルを使って，TellusのAPIについて学びます．
example>notebooks>TellusAPI利用サンプル，を選択します．
ここでは，JAXAの「だいち１号」（ALOS）の観測画像の取得方法，表示方法，データのプロファイルおよび保存方法を学びます．

Tellusの衛星画像等のAPIについては，下記のサイトを参考にしてください．
Tellus APIリファレンスサイト

つぎに，同じくJAXAの「だいち２号」（ALOS-2）に搭載した電波センサ（PALSAR-2)の観測画像を取得します．

##　「だいち２号」のSARセンサの観測画像．
Tellus APIを参考に「だいち２号」の画像を取得します．
　だいち２号のPALSAR-2の観測画像は，Google MapやGoogle Earthで見られるような写真のような画像（光学観測）ではなく，白黒の電波画像になります．
　これは，カラー画像を白黒にしているのではなく，衛星より電波を放出し，地表面を反射した電波を受信し，その信号強度を画像にしています．そのため，地表面の状態を情報として取得しています．
　また，　電波は衛星から地表面に向かって垂直に放出するのではなく，斜め方向から放出しています．そのため，地表面からは衛星の方向へ反射する電波を受信しています．たとえば，　池や湖などの水域だと，衛星より放出された電波の多くは衛星とは反対方向に反射するため，その受信強度は小さくなります．一方，建物などの構造物がある場合は，建物による二重反射によってその受信強度は大きくなります．このとき，構造物によって電波の偏光が異なり，偏光の依存性をみることで，構造物の種類を推定することができます．
　このように，SARセンサは地表面の情報を取得することができるのが特徴です．そのため，私はSAR画像を画像そのものではなく”情報”と考え，その信号を解析することでほしい情報を抽出しています．
　人工衛星のSAR画像の特徴については，以下の画像をご参考にしてください．

参照：合成開口レーダ（SAR）のキホン～事例、分かること、センサ、衛星、波長～

空畑の記事はわかりやすくてとても良いです．
Tellusの活動促進が目的と思いますが，技術だけでなく宇宙ビジネスも含めて，宇宙関係の多くの活動が紹介されています．衛星画像についても多くの記事がありますので，ぜひご参考ください．

SARセンサと光学センサの大きな違いは，SARセンサは能動センサであるため，夜間でも，悪天候（雲）であっても地上を観測することができます．そのため，災害対応の活動の支援ツールとして使われてきています．
災害対応での活用方法になりますが，下記の資料がとてもよくできていますのでご参考にしてください．
出典：国土交通省　災害時の人工衛星活用ガイドブック
水害版・衛星基礎編

では，この「だいち２号」の観測画像を，Tellusでの取得，表示および保存方法について紹介します．

##　Tellusで「だいち２号」の画像情報を取得する．

総合開発環境であるJupyter Labを使って「だいち２号」の観測画像の情報を取得しします．
使用言語はPython3です．
まずは使用するモジュールや，TOKEN情報のインプットです．

import requests
import io
from PIL import Image
import numpy as np


import requests
import json
import math
from skimage import io
from io import BytesIO
%matplotlib inline

TOKEN = '************************'
DOMAIN = 'tellusxdp.com'

Tellus APIのTokenは，以下のサイトより取得します．

ユーザー　マイページより「APIアクセス設定」を選択します．そして，APIアクセストークンにメールにて連絡のあった「APIキー」を入力し，作成をクリックすると「トークン」が出力されます．このトークンを上記のTOKENの＊＊＊＊に入力してください．

衛星画像はそのままではファイルサイズが大きすぎるため，表示サイズによって適切なサイズに分割（タイル化）されています．
そのため，衛星画像の拡大率（zoom）と，関心のある領域のタイル番号が必要になります．

こちらのサイトよりタイル番号を取得することができますが，複数のページを開いての作業になり手間がかかります．

Tellusの画像取得について，参考になるサイトがありましたので，こちらを利用させていただきました．その他にも，タイル化処理などについても参考にさせていただきました．ありがとうございます．
参考：Tellusで超低高度衛星つばめの衛星画像タイルを取得して連結してみた

from IPython.display import HTML
HTML(r'<iframe width="960" height="480" src="https://tellus-map-test.herokuapp.com/" frameborder="0"></iframe>')

上記のサイトより，関心のある地域の緯度経度情報，およびタイル情報を取得します．
では，アプリケーションを実行します．

画面にある”範囲指定”の横にある数字はZoom率を表しています．
「だいち２号」のPALSAR-2の観測画像に用意されているタイル座標のzoom率は９から１６になります．

関心のある地域をダブルクリックで拡大していき，ちょうどよいサイズになったところで”範囲指定”をクリックし，範囲をマウスで選択します．
ここでは，試しに東京都心域を選んでみました．

そうすると，選択した地域の緯度経度とタイル座標の情報を得ることができます．

次に，その地域のPARSAR-2の観測画像のメタ情報を取得します．

def getPalsar2Scene(token,area):
    url = 'https://gisapi.tellusxdp.com/api/v1/palsar2/scene'
    payload = area
    headers = { 'Authorization': 'Bearer %s' % token }
    r = requests.get(url, params=payload, headers=headers)
    #print(url)
    if r.status_code is not 200:
        raise ValueError('status error({}).'.format(r.status_code))
    return r.json()

res=getPalsar2Scene(TOKEN,{"min_lat":35.644,"min_lon":139.689,"max_lat":35.754,"max_lon":139.843})
#print(res.text)
print('number of scenes : {}'.format(len(res)))
#res[2]

for i in range(len(res)):
    print(res[i]['acquisitionDate'])

#print("")
#for i in range(len(res)):
#    print(res[i]['thumbs_url'])

このとき，getPalsar2Scene関数の"area"に，先程取得した緯度経度情報を貼り付け実行します．
出力が以下となります．

number of scenes : 4
Thu, 29 Nov 2018 02:43:05 GMT
Thu, 18 Oct 2018 02:43:04 GMT
Thu, 28 Jun 2018 02:43:03 GMT
Thu, 14 Jun 2018 02:43:04 GMT

これは，同地域のPALSAR−２の観測画像が４枚あり，一番古いものは2018年6月14日，新しいものが2018年11月29日であることがわかります．

それでは，以下を実行し，一番新しい観測画像のメタ情報すべてを取得してみます．

res[0]

実行結果

{'acquisitionDate': 'Thu, 29 Nov 2018 02:43:05 GMT',
 'censor': 'HH',
 'clat': 35.617,
 'clon': 139.659,
 'cloudCover': 'N/A',
 'entityId': 'ALOS2243962900',
 'max_lat': 35.963,
 'max_lon': 140.034,
 'min_lat': 35.268,
 'min_lon': 139.288,
 'path': 0,
 'productId': 'SARD000000232681-00007-066-011',
 'row': 0,
 'sceneId': 'ALOS2243962900-181129',
 'thumbs_url': 'https://tile.tellusxdp.com/thums/PALSAR2/ALOS2243962900.PNG'}

観測画像の全体の緯度経度や，偏波，画像のID情報を得ることができました．
詳しくは，Tellus APIのデータリファレンスをご参照ください．

次に，取得したPALSAR-2の観測画像を見てみます．

##　Tellusで「だいち２号」の観測画像を閲覧する．
先に取得した衛星画像情報より，”ALOS2243962900”のサムネイル画像（全体の簡易画像）を見てみます．

こちらを実行します．

thumbs = io.imread(res[0]['thumbs_url'])
io.imshow(thumbs)

東京地域の衛星画像であることがわかります．
信号強度が高い部分は白く表示され，その地域は建物が密集していることがわかります．一方，黒い部分は信号強度が低いことを表しており，海や川などの入射された信号の多くが衛星とは異なる方向に反射されるためです．

次に，先程選択した地域のタイル画像を拡大して取得します．

def get_Palsar2_scene_tile(scene_id, zoom, xtile, ytile):
    url = 'https://gisapi.{}/PALSAR-2/{}/{}/{}/{}.png'.format(
        DOMAIN, scene_id, zoom, xtile, ytile)
    headers = {
        'Authorization': 'Bearer ' + TOKEN
    }

    r = requests.get(url, headers=headers)
    #print(url)
    if r.status_code is not 200:
        raise ValueError('status error({}).'.format(r.status_code))
    return io.imread(BytesIO(r.content))

def getPalsar2Img(zoom,min_v,max_v,min_h,max_h):

    img = None

    for item_v in range(min_v,max_v+1):
        img_h = None
        for item_h in range(min_h,max_h+1):
            if(img_h is None):
                img_h = get_Palsar2_scene_tile(res[0]['entityId'], zoom, item_h, item_v)
            else:
                img_h=np.concatenate((img_h, get_Palsar2_scene_tile(res[0]['entityId'], zoom, item_h, item_v)), axis=1)
        if (img is None):
            img = img_h
        else:
            img = np.concatenate((img, img_h), axis=0)

    io.imshow(img)
    print(res[0]['acquisitionDate'])
    print(img.shape)
    return img

img1 = getPalsar2Img(11,805,806,1818,1819)

実行結果が以下となります．

１枚のタイル画像は256x256ピクセルであり，信号強度のダイナミックレンジは8bitになります．
これで，関心がある領域のPALSAR-2の観測画像を取得することができました．

ここで取得した画像を今後の解析にて用いるかもしれないため保存します．

Image.fromarray(img1).save('alos2_sample_tokyo.png')

保存した画像は，同じフォルダ内に画像ファイルとして保存されているのがわかります．

次に，TellusによるALOS-2の画像解析の例を紹介します．

##　Tellusで「だいち２号」の観測画像を解析する．
PALSAR-2の観測画像の解析例として，異なる時期に取得した画像の比較します．
異なる時期のSAR画像の比較方法としてよく用いられているRGB合成を行います．

これは，古い画像をRedに，新しい画像をGreenとBlueにした合成画像を作成することで，信号強度の変化を視覚化します．
新しい画像は先程取得していますので，古い画像を下記を実行し取得します．

def getPalsar2Img_1(zoom,min_v,max_v,min_h,max_h):

    img = None

    for item_v in range(min_v,max_v+1):
        img_h = None
        for item_h in range(min_h,max_h+1):
            if(img_h is None):
                img_h = get_Palsar2_scene_tile(res[-2]['entityId'], zoom, item_h, item_v)
            else:
                img_h=np.concatenate((img_h, get_Palsar2_scene_tile(res[-2]['entityId'], zoom, item_h, item_v)), axis=1)
        if (img is None):
            img = img_h
        else:
            img = np.concatenate((img, img_h), axis=0)

    io.imshow(img)
    print(res[-2]['acquisitionDate'])
    print(img.shape)
    return img

res[0]をres[-2]にすることで，関心地域の２番めに古い画像を取得しています．
一番古い画像を試してみましたが，そのタイル画像がなかったため，とりあえず”２番目”に古い画像にしています．
どの時期にするか，色々試してみてください．

ここで取得した画像は以下となります．

期間が半年ほどしか変わらないため，それほど大きな変化はありませんね．

では，この２時期の画像のRGB合成画像を取得します．

im1_swap = img1.copy()
im2_swap = img2.copy()
im2_swap[:, :, 0] = im1_swap[:, :, 0]
pil_img_swap = Image.fromarray(im2_swap)

赤や青に色付けされた部分が変化箇所になります．
湾岸地域は開発が進められているため青や赤の箇所がみられますが，都心部にはあまり色付けされた変化箇所はみられませんでした．東京オリンピックの開発が進められている国立競技場付近を拡大(zoom率上げる）して比較すると，その違いが見られるかもしれません．試してみてください．

この合成画像を今後のために保存しておきます．

Image.fromarray(im2_swap).save('alos2_sample_tokyo_rgb.png')

以上，Tellusの総合開発環境を用いて，「だいち２号」のPALSAR-2の観測画像の取得や，異なる時期の画像を比較することで都心部の開発状況の情報を求めてみました．
同様の方法にて，みなさんの住む街の異なる時期の比較や，取得した画像を解析することで新しい情報を得るなど試してみてください．

ご意見などいただけると嬉しいです．

注意：衛星画像そのものはTellus環境でのみ無料で利用することができるそうです．データポリシーについては，Tellusのサイトをご参考ください．

##　謝辞
　Tellusの画像解析については，下記のページを参考にさせていただきました．ありがとうございました．
Tellusで超低高度衛星つばめの衛星画像タイルを取得して連結してみた

##　参考情報
Tellus
日本初のオープン&フリーなプラットフォーム「Tellus」を使ってみた
 無料で最新の衛星画像を入手する方法．
【ゼロからのTellusの使い方】Jupyter Notebookで衛星データからGeoJSONで定義した領域を抜き出す