記事作成時から、Tellusからデータを検索・取得するAPIが変更になっております。該当箇所のコードについては、以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/access/traveler_api_20220310_
firstpart.html

2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html

2021年2月25日に公開された、2枚の衛星画像から街や建物の変化を自動抽出する「差分抽出ツール/Tellus-DEtector on Urban ChangE（以下、Tellus-DEUCE）」プロトタイプ。

先行してTellus OSでの使い方をご紹介していましたが、本記事では、Tellus-DEUCEのAPIの使い方についてご紹介します。本ツール利用に当たり、TellusマーケットでAPIを購入（無料）いただく必要がありますので、購入がまだの方は、「2枚の衛星画像から街や建物の変化を自動抽出！Tellus新ツール「Tellus-DEUCE」提供開始」を参考に、こちらのリンクより利用登録してください。

地球観測衛星は、災害時の現状把握、また、広範囲での状況把握のために利用されています。

例えばJAXAでは、東日本大震災から10年を迎えた今年、「震災から10年を迎えて ~宇宙から見た復興状況~」という地球観測衛星データの解析結果から見る復興・復旧の様子を紹介するページを公開しました。

今回の記事では、Tellus-DEUCEを用いて、東日本大震災の影響を受けた一部の地域の発災前後の衛星画像から、その差分を可視化しています。本記事を通して地球観測衛星について、利活用方法をご検討いただけますと幸いです。

東日本大震災の被害を受けられた地域の皆様に、謹んでお見舞い申しあげます。

1. 利用する衛星データとAPIについて

衛星データについて

Tellus-DEUCEはAVNIR-2とALOS-3相当データの差分を可視化することができます。本記事では、AVNIR-2のデータを用いて2枚の衛星画像の差分を可視化します。

AVNIR-2を搭載するALOS(だいち)衛星は、2006年1月から2011年4月まで運用されたJAXAの衛星です。同衛星は地上分解能が10ｍで、全世界をベースマップ的に撮影していました。

今回は一部の地域のみの解析を行っておりますが、関心のある方は座標を変えながら、どのような変化があったのか、ご自身でご確認ください。

AVNIR-2のAPIの詳しい使い方については、「【ゼロからのTellusの使い方】AVNIR-2のオルソ補正画像プロダクトを表示させる」もご参照下さい。

通常手順と同じく、取得する画像の緯度経度を指定します。DEUCEのAPIの検索では、あくまで指定した範囲で利用可能な画像情報を抽出してくるだけです。

そのため、差分画像作成のために、
・差分する上で適した画像かどうか（例えば雲がかかっていないか）
・比較する画像の時期や実際に重なっているかどうか
について、関数を作成していく必要があります。
まずは適した画像を検索する手順についてご紹介します。

APIについて

Tellus-DEUCEのAPIの概要は以下です。
一度に差分抽出できる面積：約6〜7km四方範囲
利用可能なデータ：AVNIR-2/ALOS-3相当データ
画像の作成枚数：10件

なお、利用可能な衛星データに関して、画像サイズ等については以下のようになっています。
・AVNIR-2：約2MB/枚、2006年から2011年までの世界中の画像
・ALOS-3相当データ：約20MB/枚、サンプル的に以下の画像をご利用いただけます
①長野市：2016/06/03 – 2019/12/25の期間中のモザイク画像3枚
②福岡市：2015/04/18 – 2019/5/15の期間中のモザイク画像3枚
③ハノイ：2017/12/20 – 2019/11/25の期間中のモザイク画像2枚
④マニラ：2018/05/12 – 2019/10/13の期間中のモザイク画像2枚

2. 被雲率を考慮して差分を抽出する

では、早速差分抽出していきましょう。
まずは、簡単に、被雲率のみを考慮して差分抽出する衛星画像を選択してみます。

比較する上で、基準とする画像をまず検索します。

# 範囲検索
import requests, json

TOKEN = 'ご自身のトークンを貼り付けてください'

def get_market_token(payload={}):
    url = 'https://www.tellusxdp.com/api/manager/v2/auth/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 requests.codes.ok:
        print(r.content)
    return json.loads(r.content) 

def search(token, satellite, left_bottom_lat, left_bottom_lon, right_top_lat, right_top_lon):
    url = 'https://deuce.tellusxdp.com/api/v1/{}/search?left_bottom_lat={}&left_bottom_lon={}&right_top_lat={}&right_top_lon={}'.format(satellite, left_bottom_lat, left_bottom_lon, right_top_lat, right_top_lon)

    headers = {
        'Authorization': 'Bearer ' + token
    }
    print(url)
    r = requests.get(url, headers=headers)
    if r.status_code is not requests.codes.ok:
        print(r.status_code)

    return json.loads(r.content)

#下記はTellus-DEUCEのプロダクトIDになります
ret = get_market_token({'product_id': '2f443df3-3355-4232-99bc-01dbe5dc6ca3'})
# print(ret)

satellite = 'avnir2'
# left_bottom_lat = 33.54886871475409
# left_bottom_lon = 130.3163784981599
# right_top_lat = 33.628136349724286
# right_top_lon = 130.45989661632703

# 南相馬市周辺：37.665909621315194, 141.01116732251953
left_bottom_lat = 37.365909621315194
left_bottom_lon = 140.01116732251953
right_top_lat = 37.665909621315194
right_top_lon = 141.01116732251953

executed = search(ret['token'], satellite, left_bottom_lat, left_bottom_lon, right_top_lat, right_top_lon)
# print(executed)

座標を基準にした検索結果から各画像に割り当てられたIDのみを取得します。
これらのIDは実際に差分を検証する際に用いるものとなります。

sceneIds = [scene['scene_id'] for scene in executed['scenes']]
print(sceneIds[0:5])

衛星画像はメタ情報として、緯度経度だけではなく様々な方法を含んでいます。

IDを用いて、差分検証を行う前に必要な情報を抽出することにします。
まずはAVNIR-2のメタ情報を抜き出すための関数を定義しましょう。

被雲率が10%以下、検索結果のIDのみを検索します。
雲量のコントロールはなかなか難しく、最終的には目視確認が重要になります。
例えば、被雲率が90%だったとしても、自分が差分をとりたい場所には雲が全くかかっていない可能性あります。しかし、画像のメタ情報では、画像全体として含まれている雲量を基準としているため、値から判断することができません。そのため、時期を合わせて差分の検討をしたい場合には、雲量で制限をするよりも、まずは比較したい時期で画像を検索するのが良いでしょう。

# AVNIR-2の画像情報を検索
def search_file(params={}, next_url=''):
    if len(next_url) > 0:
        url = next_url
    else:
        url = 'https://file.tellusxdp.com/api/v1/origin/search/avnir2-ori'
    
    headers = {
        'Authorization': 'Bearer ' + TOKEN
    }
    r = requests.get(url, params=params, headers=headers)
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return r.json()

# メタ情報を抜き出す
# 雲10%以下、該当のIDのみを取得
def getMeta(sceneIds,cc=10):
    metas = [search_file({'cloud_cover': cc,'scene_id': sceneId}) for sceneId in sceneIds]
    return metas

# 該当する結果を表示する
entityMetas = getMeta(sceneIds)
len(entityMetas)

実際には、画像が含まれているcount=1になっているものだけを利用します。

existedIdsLists = [result for result in entityMetas if result['count'] == 1]
print(len(existedIdsLists))

メタ情報を確認します。

画像の重ね合わせでは、bbox、rsp_frame_number、rsp_path_numberが重要となります。

entityMetas[0]

例えば下記のような結果を得られるでしょう。

{‘count’: 1,
‘items’: [{‘bbox’: [140.4638935, 36.9711292, 141.5258533, 37.8406316],
‘center_datetime’: ‘2011-04-15T01:33:15.286677+00:00’,
‘cloud_cover’: 0,
‘coordinates’: [[140.4638935, 37.1412136],
[141.2984209, 36.9711292],
[141.5258533, 37.6689818],
[140.6840511, 37.8406316]],
‘dataset_id’: ‘ALAV2A278172840-OORIRFU-D070P0-20110415-002’,
‘date_added’: ‘2020-01-10T15:42:32.298772+00:00’,
‘files’: [‘1784’, ‘61248471’, ‘61248471’, ‘61248471’, ‘61248471’],
‘frame_number’: 2840,
‘gain’: [0.588, 0.573, 0.502, 0.835],
‘observation_datetime’: ‘2011-04-15T01:33:15.286677+00:00’,
‘offset’: [0.0, 0.0, 0.0, 0.0],
‘orbit_direction’: ‘D’,
‘pointing_angle’: 12.5,
‘publish_link’: ‘https://file.tellusxdp.com/api/v1/origin/publish/avnir2-ori/ALAV2A278172840-OORIRFU-D070P0-20110415-002’,
‘rsp_frame_number’: 2840,
‘rsp_path_number’: 70,
‘scene_id’: ‘ALAV2A278172840’,
‘scene_shift’: 0,
‘sun_azimuth_angle’: 148.9248192,
‘sun_elevation_angle’: 58.7601858,
‘total_orbit_number’: 27817}]}

関数を定義して、必要な情報のみを抽出し、さらにデータフレームをジオデータフレームに変えます。

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Point, Polygon, MultiPoint
from datetime import datetime

def makeGeoDataFrame(IdLists):
    '''
    Create geodataframe from the list including scenes' information
    '''
    ids = [id['items'][0]['scene_id'] for id in IdLists]
    bbox = [id['items'][0]['bbox'] for id in IdLists]
    ccr = [id['items'][0]['cloud_cover'] for id in IdLists]
    frameNum = [id['items'][0]['rsp_frame_number'] for id in IdLists]
    pathNum = [id['items'][0]['rsp_path_number'] for id in IdLists]
    ob_time = [datetime.date(datetime.strptime(id['items'][0]['observation_datetime'][0:10], '%Y-%m-%d')) for id in IdLists]
    coords = []
    poly = []
    points = []
    for i in range(len(bbox)):
        points.append(MultiPoint([Point(bbox[i][0],bbox[i][1]),Point(bbox[i][2],bbox[i][1]),\
                                  Point(bbox[i][2],bbox[i][3]),Point(bbox[i][0],bbox[i][3])]))
        coords.append([(bbox[i][0],bbox[i][1]),(bbox[i][2],bbox[i][1]),\
                       (bbox[i][2],bbox[i][3]),(bbox[i][0],bbox[i][3])])
        poly.append(Polygon(coords[i]))
    df = pd.DataFrame(
        {'ID': ids,
         'Time': ob_time,
         'Cloud_Cover': ccr,
         'Frame_Number': frameNum,
         'Path_Number': pathNum,
         'Points': points,
         'Coordinates': poly})
    gdf = gpd.GeoDataFrame(df, geometry='Coordinates')
    return gdf

gdf = makeGeoDataFrame(existedIdsLists)
gdf.head()

今回はALAV2A277442840をベース画像として捉え、比較するための古い画像としてALAV2A082852840のシーンIDを用います。

scene1 = gdf.loc[gdf['ID']=='ALAV2A277442840']
scene1.reset_index(drop=True, inplace=True)
scene2 = gdf.loc[gdf['ID']=='ALAV2A082852840']
scene2.reset_index(drop=True, inplace=True)
fig, ax = plt.subplots()
scene1.plot(ax=ax, facecolor='gray');
scene2.plot(ax=ax, facecolor='red',alpha=0.8);
plt.tight_layout();

ベース画像（灰色）に比較対象画像（赤色）が含まれていることが分かります。このように重なる範囲が大きいペアであれば、簡単に差分を検証することができます。

続いて、ベース画像内に、対象画像の四つの座標（左下、右下、右上、左上）が含まれているのかを確認しましょう。

先ずは簡単な例から試して見ます。

geopandasを用いる方法、shapelyを用いる方法の二つがありますが、geopandasはshapelyを介して指定した点が画像内（ポリゴン内）に含まれているのかを確認しています。そのため、shapelyをそのまま利用する方が早いですが、対象とする画像が少ない場合には大きな差にはならないかと思います。

下ではどちらの例も示しています。状況に応じて使い分けしてください。

from geopandas.geoseries import *
g1 = GeoSeries([scene2.Points[0][0],scene2.Points[0][1],scene2.Points[0][2],scene2.Points[0][3]])

poly = GeoSeries(scene1.geometry)
print(g1.intersects(poly.iloc[0]))

for p in [scene2.Points[0][0],scene2.Points[0][1],scene2.Points[0][2],scene2.Points[0][3]]:
    print(scene1.geometry.intersects(p))

上記の実行結果から、一ヶ所（左上の座標）がベース画像内に含まれていることが分かりました。

実際に差分検証のテストをしてみましょう。

# 差分検証
def difference(token, satellite, defference_laft_bottom_lat, defference_left_bottom_lon, defference_right_top_lat, defference_right_top_lon, scene_id_a, scene_id_b):
    url = 'https://deuce.tellusxdp.com/api/v1/{}/difference/scene_id?left_bottom_lat={}&left_bottom_lon={}&right_top_lat={}&right_top_lon={}&scene_id_a={}&scene_id_b={}'.format(satellite, defference_laft_bottom_lat, defference_left_bottom_lon, defference_right_top_lat, defference_right_top_lon, scene_id_a, scene_id_b)

    headers = {
        'Authorization': 'Bearer ' + token
    }
    print(url)
    r = requests.get(url, headers=headers)
    if r.status_code is not requests.codes.ok:
        print(r.status_code)

    return json.loads(r.content) 

さらに差分を作成するに当たって幾つかの条件があります。

一つ目に、差分を検証できる範囲に制限があります。こちらはleft_bottom_lon、right_top_lonの差分は0.07以内、left_bottom_lat、right_top_latの差分は0.06以内であることとなっています。
面積で言えば約6〜7km四方までの差分画像を作成できるということになります。画像サイズはAVNIR-2が約2MB、ALOS-3相当画像が約20MBです。さらに、画像の作成は10件までしかできません。そのため、一度画像を作成したらダウンロードを行い、そのタスクは消すという流れをとるのが良いと思います。

# 140.93435, 37.86902
defference_laft_bottom_lat = 37.86902-0.03
defference_left_bottom_lon = 140.93435-0.035
defference_right_top_lat = 37.86902+0.03
defference_right_top_lon = 140.93435+0.035

scene_id_a = 'ALAV2A277442840'
scene_id_b = 'ALAV2A082852840'

difference_data = difference(ret['token'], satellite, defference_laft_bottom_lat, defference_left_bottom_lon, defference_right_top_lat, defference_right_top_lon, scene_id_a, scene_id_b)
print(difference_data)

コードを実行すると、
＊{‘error’: ‘Unauthorized (invalid signature)’}
が表示されることがあります。
これはget_market_tokenで取得したトークンが期限切れになっているために生じています。
この場合にはコードを再度実行し、トークンを取得する必要があります。

エラーが表示されることなく、動いたことが確かめられたでしょうか？

access_keyは、結果の画像を取得する際に用いることになります。実際に画像の作成からダウンロードまでを行います。

最初に、現在作成している処理を見て見ましょう。先に実行したアクセスキーが示されているはずです。

# 作成一覧
def getWorks(token):
    url = 'https://deuce.tellusxdp.com/api/v1/works'

    headers = {
        'Authorization': 'Bearer ' + token
    }
    print(url)
    r = requests.get(url, headers=headers)
    if r.status_code is not requests.codes.ok:
        print(r.status_code)

    return json.loads(r.content) 

work_lists = getWorks(ret['token'])

print(work_lists)

# 作成詳細
def getWork(token, access_key):
    url = 'https://deuce.tellusxdp.com/api/v1/work/{}'.format(access_key)

    headers = {
        'Authorization': 'Bearer ' + token
    }
    print(url)
    r = requests.get(url, headers=headers)
    if r.status_code is not requests.codes.ok:
        print(r.status_code)

    return json.loads(r.content) 

work_data = getWork(ret['token'], work_lists[0]['access_key'])
print(work_data)

実行結果のtiff画像をダウンロードします。

# tiffのダウンロード
def download_file(token, access_key, tiff_type):
    url = 'https://deuce.tellusxdp.com/api/v1/{}/image/{}.tiff'.format(access_key, tiff_type)

    headers = {
        'Authorization': 'Bearer ' + token
    }
    r = requests.get(url, headers=headers)
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return r.content

def save_file(data, name):
    with open(name, "wb") as f:
        f.write(data)

# 画像の保存。適宜名前や保存するフォルダを変更して下さい。
# base
file = download_file(ret['token'], work_lists[0]['access_key'], 'base')
save_file(file, 'base.tiff')

# cmap
file = download_file(ret['token'], work_lists[0]['access_key'], 'cmap')
save_file(file, 'cmap.tiff')

実行すると、以下のような画像が自身の開発環境にダウンロードされるかと思います。

3. 衛星の撮像条件を考慮して差分を抽出する

被雲率を考慮して差分を抽出だけでは、うまく差分を得られない場合もあります。
そのため、今回は雲量を無視し、時期を考慮して差分を検証していきましょう。
さらに、rsp_frame_numberやrsp_path_numberもベース画像と合わせることにより、検索対象を絞ります。

先に比較のペアを見つけるための関数を定義します。

# メタ情報を抜き出す
# 雲100%以下、該当のIDのみを取得
entityMetas = getMeta(sceneIds,100)

existedIdsLists = [result for result in entityMetas if result['count'] == 1]

gdf = makeGeoDataFrame(existedIdsLists)
gdf.head()

ベースとする画像は先ほどと同じく2011-04-10の画像を選択します。

baseScene = 'ALAV2A277442840'

先ずはベース画像より時期が早いものは全て除外しましょう。

def getPairs(sceneId,geoDf):
    '''
    sceneId: Set the base id to compare the other IDs.
    geoDF: Set a GeoDataFrame including sceneId
    '''
    idsExtract = geoDf[geoDf.Time < list(geoDf.loc[geoDf['ID']==sceneId]['Time'])[0]]
    return idsExtract

getPairs(baseScene,gdf)

実行すると、候補となるペアのみがデータフレームとして保存されます。
上の候補は全てが差分を検証できる画像リストとなっていますが、気をつけなければならないのは、その切り取り範囲です。つまり、この記事のコードで、defference_laft_bottom_latのように指定された座標群である四隅の値です。

衛星画像は、軌道傾斜角が傾いているため、南北に対して真四角な画像として撮影されていません（上のAVNIR-2の画像をご覧ください。斜めに撮影されているのが分かるはずです）。

ただし、bboxは衛星画像の四隅を基にして作成した真四角の範囲となるため、実際には画像情報がない部分も含まれています。そのため上のリストにあるのに、差分をとることが難しい場合があります。
重なりが非常に狭い場合には、この問題が顕著になるため、撮影場所をほぼ同じにする必要があります。

JAXAのこちらのサイトで確認をすると、どの画像が同じ場所にあるのかを確認することができます。

では、bboxを基にして対象となる画像を選び出します。

def judgePairs(base_id, target_id, geoDf):
    '''
    base_id: ID is for the base image
    target_id: ID si for the target image
    geoDf: GeoDataFrame including the base and the target image
    '''
    baseDf = geoDf.loc[geoDf['ID']==base_id]
    targetDf = geoDf.loc[geoDf['ID']==target_id]
    target_points = GeoSeries([list(targetDf.Points)[0][0], list(targetDf.Points)[0][1],\
                               list(targetDf.Points)[0][2], list(targetDf.Points)[0][3]])
    base_poly = GeoSeries(baseDf.geometry)
    results = target_points.intersects(base_poly.iloc[0])
    if (results[0] | results[1] | results[2] | results[3]) == True:
        return target_id

関数を実行します。実行結果に先ほど例で示した`ALAV2A082852840`も含まれていることが確認できます。

idExclude = list(gdf.loc[gdf['ID']!=baseScene]['ID']) # Remove the base ID
idTargets = [judgePairs(baseScene,target,gdf)for target in idExclude] # Obtain the other IDs
results = [] 
for val in idTargets:
    if val != None : 
        results.append(val) 
print(results)

差分を作成する対象となる画像で、改めてデータフレームを作成します。

def extractDataFrame(refDf,ids):
    ID = []
    Time = []
    FrameNum = []
    PathNum = []
    Cloud_Cover = []
    Points = []
    Coordinates = []

    for result in ids:
        ID += list(refDf.loc[refDf['ID'] == result].ID)
        Time += list(refDf.loc[refDf['ID'] == result].Time)
        FrameNum += list(refDf.loc[refDf['ID'] == result].Frame_Number)
        PathNum += list(refDf.loc[refDf['ID'] == result].Path_Number)
        Cloud_Cover += list(refDf.loc[refDf['ID'] == result].Cloud_Cover)
        Points += list(refDf.loc[refDf['ID'] == result].Points)
        Coordinates += list(refDf.loc[refDf['ID'] == result].Coordinates)

    return pd.DataFrame(
        data={'ID': ID, 'Time': Time,'Frame_Number': FrameNum, 'Path_Number': PathNum, 'Cloud_Cover': Cloud_Cover, 'Points': Points, 'Coordinates': Coordinates},
        columns=['ID', 'Time','Frame_Number','Path_Number', 'Cloud_Cover', 'Points', 'Coordinates']
    )

マッチさせる値はベース画像に合わせます。

resDf = extractDataFrame(gdf,results)
resDfmatch = resDf.loc[(resDf.Frame_Number == list(gdf.loc[gdf.ID == baseScene,'Frame_Number'])[0]) &\
          (resDf.Path_Number == list(gdf.loc[gdf.ID == baseScene,'Path_Number'])[0])]

# 結果を確認する
resDfmatch

上記の検索で引っかかった画像に対して、サムネイル画像を表示して、雲のかかり具合などを目視でチェックしましょう。

def get_AVNIR2_scene(id):
    url = 'https://gisapi.tellusxdp.com/api/v1/av2ori/get_scene/{}'.format(id)
    headers = {
        'Authorization': 'Bearer ' + TOKEN
    }

    r = requests.get(url, headers=headers)
    if r.status_code is not 200:
        raise ValueError('status error({}).'.format(r.status_code))
    return r.json()

sceneIds = [get_AVNIR2_scene(id) for id in list(resDfmatch.ID)]
d067p3Lists = [thumbUrl['thumbs_url'] for thumbUrl in sceneIds]

from skimage import io
from io import BytesIO

thumbs = [io.imread(imageList) for imageList in d067p3Lists]
titles = [Id['entityId'] for Id in sceneIds]
fig = plt.figure(figsize=(10, 10))
ax = []
columns = 3
rows = 2
img = np.zeros(thumbs[0].shape) # for a blank image

for i in range(columns*rows):
    if(len(thumbs) < columns*rows):
        if(i < columns*rows-1):
            ax.append(fig.add_subplot(rows, columns, i+1))
            ax[-1].set_title(titles[i],size = 10)  # set title
            plt.imshow(thumbs[i])
            print(i)
        else:
            ax.append(fig.add_subplot(rows, columns, i+1))
            plt.imshow(img)
            print(i)
    else:
        ax.append(fig.add_subplot(rows, columns, i+1) )
        ax[-1].set_title(titles[i],size = 10)  # set title
        plt.imshow(thumbs[i])           

plt.show()

ベースは変えずに対象の画像を`ALAV2A176792840`にして実行します。

時期的にも雲量的にもこちらの画像は良さそうです。

それでは再度、差分処理を行ってみましょう。

# ret = get_market_token({'product_id': '2f443df3-3355-4232-99bc-01dbe5dc6ca3'})

defference_laft_bottom_lat = 37.86902-0.03
defference_left_bottom_lon = 140.93435-0.035
defference_right_top_lat = 37.86902+0.03
defference_right_top_lon = 140.93435+0.035

scene_id_a = 'ALAV2A277442840'
scene_id_b = 'ALAV2A176792840'

difference_data = difference(ret['token'], satellite, defference_laft_bottom_lat, defference_left_bottom_lon, defference_right_top_lat, defference_right_top_lon, scene_id_a, scene_id_b)
print(difference_data)

エラーがでなければ`access_key`が返ってきます。
エラーが発生した場合には数分待って再度実行してみてください。

print(difference_data['access_key'])

# 画像の保存。適宜名前や保存するフォルダを変更して下さい。
# base
file = download_file(ret['token'], difference_data['access_key'], 'base')
save_file(file, 'base2.tiff')

# cmap
file = download_file(ret['token'], difference_data['access_key'], 'cmap')
save_file(file, 'cmap2.tiff')

こちらは、前に取得した画像と比べ、時期が似ているために地上の変化は抑えられていますが、雲を差分として検出するため左上に大きな変化があると誤認識しています。
雲のかかっている範囲を除くと、Tellus-DEUCEのアルゴリズムを構築する際に都市近郊の変化を抽出するように試みたために、元々緑地だったところが裸地になった場所については変化を検出せず、湾岸や建物があった範囲で黄色や赤色になっていることが分かります。

このように差分検出は雲の有無も変化として捉えてしまうことがありますので、適切な画像を選び、試行錯誤をする必要があります。

実際にどのような変化があったのか、ということについては、国土地理院のページやJAXAのページでも言及されています。

詳細な確認は、国土地理院のページで公開されている空中写真で確認してみるのもよいでしょう。

4. 保存した衛星画像を削除する

開発環境に保存できるデータ容量は限られているため、処理の終った結果を削除するには、以下の関数を実行して下さい。

# 作成情報削除
def deleteWork(token, access_key):
    url = 'https://deuce.tellusxdp.com/api/v1/work/{}'.format(access_key)

    headers = {
        'Authorization': 'Bearer ' + token
    }
    r = requests.delete(url, headers=headers)
    print(r.url)
    if r.status_code is not 200:
        print(r.status_code)

    return json.loads(r.content) 

work_lists = getWorks(ret['token'])
deleteWork(ret['token'], work_lists[0]['access_key'])

5. まとめ

以上が、Tellus-DEUCEをAPI経由で利用する方法です。
本記事では一部地域の変化しか見ていませんが、地域を変えながら、解析結果を組み合わせることでより広範囲を見ることもできます。
「衛星データもIoTの一つ？航空機・ドローンとの比較と事例紹介」でもご紹介していますが、衛星画像と航空機写真は利用用途に応じて組み合わせて使うことが重要です。
これは防災に限った話ではありませんが、観測可能な範囲が限定される航空機により詳細に観測する場所を決める前に、まずは大雑把に状況を把握する、という目的で衛星画像を利用すると、より効率的に状況把握ができるでしょう。

繰り返しになりますが、東日本大震災における衛星観測結果の利用については、東日本大震災－JAXA地球観測の記録や震災から10年を迎えて ~宇宙から見た復興状況~に情報がまとめられています。

本記事では、東日本大震災を例にツールの利用方法について紹介しましたが、AVNIR-2は世界中を撮影していた衛星です。例えば新興国の発展状況を可視化するのに利用したり、山中にポツンと一軒家が突如できたところはどこかにないか、ということを探したりするのに利用してもよいでしょう。

また、衛星画像では、雲が一つ大きな課題になるのだな、ということを改めて実感された方も多くいらっしゃるかと思います。そのため、雲に予めマスクをする処理をして差分抽出で期待していない実行結果となる要因をなくす試みもあります。

衛星画像における雲の課題については、宙畑でも、「衛星データに雲が映っているか否かの画像分類を4つの手法で精度比較してみた」という記事を公開していますので、合わせてご覧ください。

本記事では、開発環境を用いた衛星画像の差分抽出方法について紹介しましたが、Tellus OSを用いても同様のことを行うことができます。
Tellus-DEUCEを使って、ぜひ世界の変化を可視化してみてください。