TellusのCOG形式のAPIを使って、ASNARO-2の画像を重ね合わせて前後の変化抽出を実施してみた
Tellusの公式データとして無料公開しているASNARO-2のSARデータについて、前後比較のサンプルコードを紹介します。
記事作成時から、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
前回の記事では、Tellus Traveler上で扱える衛星データについて、ASNARO-2を題材にデータの確認方法についてご紹介しました。
もう使った？Tellus Travelerの利用方法~ASNARO-2の災害データセットを見てみる~
Tellus Traveler上で衛星データを扱う方法としては主に2つ、ブラウザ上で閲覧する方法とAPIを利用する方法があります。以前の記事では、衛星データをまるっと1枚すべて取得する方法でご紹介しましたが、Tellus TravelerではCoG（Cloud Optimized GeoTIFF）と呼ばれるクラウド対応のGeoTIFFに対応しているため、欲しい領域のみを指定して取り出すことも可能です。
今回はこの方法についてご紹介します。
なお、本解析の実施にはTellusの開発環境の申し込みが必要となります。お問い合わせからお申し込みをお願いします。
解析の概要
今回の解析では、X-BandのSAR画像ASNARO-2を題材に、北海道胆振東部地震の被災状況を捉えた画像と復興が進む様子を捉えた画像の2枚を合成して、両者の差分を取る解析を行います。
この図はSAR画像の前後比較の原理を示しています。
前の画像を青く、後の画像を赤く色づけて画像を合成するとどちらかの画像にしかなかったものが色づいて見えます。
この例では、リンゴはBeforeの画像にはなくAfterの画像で増えているので合成した画像では「赤」で表現され、反対にBeforeの画像にあったフォークはAfterの画像ではなくなっているので合成した画像では「青」に見えます。どちらの画像にもあるお皿は合成された画像では「白」く見えます。逆にどちらの画像でも何もない（電波が返ってこない）部分は「黒」に見えます。
このように、時刻の異なる2枚の画像をそれぞれ青と赤に着色して合成することで、その時間で何が変わったのかを知ることができるのです。
解析の流れ
以下の3つのステップに沿って解析を行います。
①衛星データの検索
Tellus上で、場所の情報を指定して該当する衛星データを検索します。
②衛星データのダウンロード
①で見つけた画像をTellus上の開発環境にダウンロードします。
※今回扱うASNARO-2はTellus外へのデータの持ち出しが禁止されているデータであるため、Tellusに開発環境を用意の上、そこへダウンロードします。
③2枚のデータを合成し、PNG形式で出力する
取得してきた2枚のデータを重ね合わせて差分が見えるように表示します。
それでは、それぞれのステップ毎に解説していきます。
①衛星データの検索
まず実施するのは衛星データの検索です。
必要なライブラリをインポートします。
import os
import requests
import json条件に合致するデータを探す関数は以下の通りです。
#所定の場所のデータ検索をする関数
def search_data(payload={}):
url = '{}/api/traveler/v1/data-search/'.format(BASE_PATH)
headers = {
'Authorization': 'Bearer ' + TOKEN,
'content-type': 'application/json',
}
r = requests.post(url, headers=headers, data=json.dumps(payload))
print(r.status_code)
if r.status_code != 200:
print(r.content)
raise ValueError('status error({}).'.format(r.status_code))
return json.loads(r.content)続けて、データを探すための検索条件を設定します。
今回は北海道胆振東部地震を扱うため、該当のエリアの緯度経度情報が必要となります。
解析を実施したいエリアの緯度経度の取得方法は様々ありますが、今回はTellus OS上で該当のエリアを多角形で囲み、geojson形式でエクスポートしました。
このファイルをテキストエディタなどで開いてみると、以下のように緯度経度情報が含まれていることが確認できます。
{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[141.9371862693975,42.75848804932892],[141.93507323883878,42.755016948577975],[141.93790255094285,42.75282111399301],[141.94513699455075,42.75336021831248],[141.94639048725506,42.756489560521516],[141.94486838897123,42.75838286731437],[141.9371862693975,42.75848804932892]]]},"properties":{"_drawtype":"polygon","_note":"","_version":"2.0","_fillColor":"#fff","_fillOpacity":0.4,"_strokeColor":"#3399cc","_strokeOpacity":1,"_strokeWidth":4,"_strokeDash":[],"_displayArea":"460904.2818971103 m²","_displayDistance":"2587.080935880432 m"}}この中から、必要な情報を取り出し、検索条件としてセットします。
# 検索条件を指定してデータを検索する関数（AOIを指定）
def get_request():
ret = search_data({
'datasets': datasets,
'query': {
},
"intersects":
{"type": "Polygon", "coordinates":
[[
# AOIの指定
[141.84938200356842,42.75105653576509],
[141.83924706503996,42.71645449666897],
[141.89461385885272,42.69893933684679],
[141.96255548380265,42.72252158826697],
[141.95711264644478,42.761943222877846],
[141.89968132811697,42.76359672954325],
[141.84938200356842,42.75105653576509]
]]
},
})
return retdatasets（データセット）とあるのは、Tellus上でのデータの種類のことで、後ほど、ASNARO-2を指定します。
ここまでできたら、データの検索を実行します。
APIの実行にはトークンの発行が必要です。こちらを参照の上、それぞれのアカウントに紐づいたアカウントの発行を行ってください。
また、データセットIDはTellus上のデータセット詳細ページから確認できます。
if __name__ == '__main__':
# データダウンロード設定
TOKEN = 'XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
BASE_PATH='https://www.tellusxdp.com'
#データセットidの設定
datasets = ["10402ffd-fb70-4625-9b46-63cec87cd5a9"] #ASNARO-2
# 検索の実行
ret = get_request()
# 検索結果の確認
print(ret)実行すると、以下のような結果が得られます。
200
{'type': 'FeatureCollection', 'features': [{'dataset_id': '10402ffd-fb70-4625-9b46-63cec87cd5a9', 'geometry': {'coordinates': [[[141.85921515024788, 42.69645878203438], [141.98184752828192, 42.7119268040681], [141.9608764018023, 42.80151114290716], [141.83806053574853, 42.78604096679968], [141.85921515024788, 42.69645878203438]]], 'type': 'Polygon'}, 'id': '7d47301d-9533-4be2-a010-31a86912e865', 'type': 'Feature', 'properties': {'sar:frequency_band': 'X', 'processing:level': 'L2.1', 'sat:relative_orbit': 126, 'start_datetime': '2020-09-05T08:12:58+00:00', 'end_datetime': '2020-09-05T08:13:00+00:00', 'tellus:name': 'AS201462808498-200905', 'created': '2021-12-13T14:25:18.087579+00:00', 'tellus:can_ordered': True, 'tellus:sat_frame': 8498, 'view:off_nadir': 36.184845, 'asnaro2:beam': '1', 'sar:observation_direction': 'right', 'sar:polarizations': 'HH', 'sar:instrument_mode': 'SP2', 'tellus:published_datetime': '2021-12-17T01:28:38.867529+00:00', 'sat:orbit_state': 'ascending', 'sar:product_type': 'SLC', 'gsd': 1.0, 'updated': '2022-01-25T02:36:39.180886+00:00'}}, {'dataset_id': '10402ffd-fb70-4625-9b46-63cec87cd5a9', 'geometry': {'coordinates': [[[141.8589167520601, 42.69686080735227], [141.98154593883282, 42.71231672781694], [141.96058768164986, 42.80190194745317], [141.8377749207214, 42.78644388338905], [141.8589167520601, 42.69686080735227]]], 'type': 'Polygon'}, 'id': '2c9ccde2-3a8c-4cd2-aac5-4f0c3b33a4a4', 'type': 'Feature', 'properties': {'sar:frequency_band': 'X', 'processing:level': 'L2.1', 'sat:relative_orbit': 126, 'start_datetime': '2018-09-08T08:09:47+00:00', 'end_datetime': '2018-09-08T08:09:48+00:00', 'tellus:name': 'AS200355208498-180908', 'created': '2021-12-13T13:22:17.325305+00:00', 'tellus:can_ordered': True, 'tellus:sat_frame': 8498, 'view:off_nadir': 36.1230469, 'asnaro2:beam': '1', 'sar:observation_direction': 'right', 'sar:polarizations': 'HH', 'sar:instrument_mode': 'SP2', 'tellus:published_datetime': '2021-12-17T01:28:25.863601+00:00', 'sat:orbit_state': 'ascending', 'sar:product_type': 'SLC', 'gsd': 1.0, 'updated': '2022-01-25T02:36:48.498292+00:00'}}], 'meta': {'total': 2, 'search-condition': {'intersects': {'type': 'Polygon', 'coordinates': [[[141.84938200356842, 42.75105653576509], [141.83924706503996, 42.71645449666897], [141.89461385885272, 42.69893933684679], [141.96255548380265, 42.72252158826697], [141.95711264644478, 42.761943222877846], [141.89968132811697, 42.76359672954325], [141.84938200356842, 42.75105653576509]]]}, 'query': {}, 'datasets': ['10402ffd-fb70-4625-9b46-63cec87cd5a9'], 'paginate': {'cursor': 'ed9cc27d-6305-4252-ab7c-896cba4256f9', 'size': 1000}}}}少々見づらいですが、条件に合致するシーンが2つあることが分かります。
idが「7d47301d-9533-4be2-a010-31a86912e865」の2020-09-05に撮影されたデータと、id「2c9ccde2-3a8c-4cd2-aac5-4f0c3b33a4a4」で2018-09-08に撮影されたデータです。
胆振東部地震が発生したのは 2018年9月6日なので、その2日後の画像と2年度の画像があるということが分かります。
また、詳しく見ていくとどちらの画像も’sat:relative_orbit’: 126と’tellus:sat_frame’: 8498となっています。これは、衛星がどの軌道からどの場所を撮影したかを示していて、この数値が一致している画像同士の前後比較は比較的容易と言われています。
さらに、撮影月がどちらも9月である点もポイントです。季節が違う2枚の画像を比較すると、北海道であれば雪や、農業をしているエリアでは作付けの状況など、前後比較したい対象（今回の場合は地震）以外の情報も差分として抽出されてしまうため、可能な限り時期を合わせた比較が望ましいと考えられます。
今回は該当画像が2枚だったため、特にこれ以上絞り込む必要はありませんが、もし画像が複数存在する場合には、これらの観測条件も加味して、前後比較する2枚を決めると良いでしょう。
②衛星データのダウンロード
続いて、見つけた画像（シーン）をTellus環境上にダウンロードしてくるステップです。
Tellus travelerのAPIでは、実際にダウンロ―ドしてくるファイルIDを指定し、有効期限付きのURLを発行して、所定のファイルをダウンロードするという手順になります。
先ほど見つけた画像（シーン）は複数のファイルで構成されているため、シーンの中に含まれるファイルIDを確認する必要があります。
# ファイルID等を取得する関数
def get_dataset_data_by_id_files(base_path, token, dataset_id, data_id):
print(base_path, token, dataset_id, data_id)
url = '{}/api/traveler/v1/datasets/{}/data/{}/files/'.format(
base_path, dataset_id, data_id)
headers = {
'Authorization': 'Bearer ' + token,
'content-type': 'application/json'
}
r = requests.get(url, headers=headers)
assert r.status_code == 200
return json.loads(r.content)ファイルIDが分かったら、ダウンロード用のURLを発行します。
# ダウンロードURLを取得する関数
def get_dataset_data_by_id_files_by_id_download_url(base_path, token, dataset_id, data_id, file_id):
url = '{}/api/traveler/v1/datasets/{}/data/{}/files/{}/download-url/'.format(
base_path, dataset_id, data_id, file_id)
headers = {
'Authorization': 'Bearer ' + token,
'content-type': 'application/json'
}
r = requests.post(url, headers=headers)
assert r.status_code == 200
return json.loads(r.content)['download_url']上記の関数を組み合わせて、ダウンロードを実行する関数は以下の通りです。
# データをダウンロードする関数
def download(scenes, dataset_id, token, dist='./', base_path='https://www.tellusxdp.com'):
for scene_id in scenes:
files = get_dataset_data_by_id_files(
base_path, token, dataset_id, scene_id)
print(files)
rawdata = files['results']
path = os.path.join(dist, scene_id)
if len(rawdata) > 0:
for file in rawdata:
file_id = file['id']
file_name = file['name']
file_path = os.path.join(path, file_name)
download_url = get_dataset_data_by_id_files_by_id_download_url(
base_path, token, dataset_id, scene_id, file_id)
r = requests.get(download_url, stream=True)
if not os.path.exists(path):
os.makedirs(path)
with open(file_path, 'wb') as f:
for chunk in r.iter_content(chunk_size=1024):
if chunk:
f.write(chunk)
f.flush()前の章で見つけたデータのシーンIDを指定して、ダウンロードを実行します。
if __name__ == '__main__':
# トークン
TOKEN = 'XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
#データセットID`
dataset_id='10402ffd-fb70-4625-9b46-63cec87cd5a9' #ASNARO-2
#シーンIDの指定（01_search_data.pyの実行結果からIDを取得）
scenes = [
'7d47301d-9533-4be2-a010-31a86912e865',
'2c9ccde2-3a8c-4cd2-aac5-4f0c3b33a4a4'
]
# データのダウンロード（コードを実行した場所と同じ場所に保存される）
download(scenes, dataset_id, TOKEN)画像の取得には多少時間がかかります。
ダウンロードが完了すると、同じフォルダ内にシーン毎に画像が保存されます。
シーンの中に様々なファイルが含まれていることがお分かりいただけるかと思います。
リモートセンシングの解析に使うのは最もファイルサイズの大きい「SOR-IMG-HH-AS200355208498-180908___-SP2R1.1__A_.tif」ですが、今回は簡単のため、ブラウザ表示用に作られている「AS200355208498-180908_webcog.tif」を用いて解析を行うことにします。
③2枚のデータを合成し、PNG形式で出力する
ダウンロードしてきた画像から、今回必要な場所のみ切り出し、PNG形式で表示させます。
こちらもまずはライブラリをインポートします。
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image指定してきた範囲を切り出す関数は以下の通りです。
# Geotiffから指定範囲（緯度経度）を切り出す関数
def cut4_geotiff(IMG_IN, IMG_OUT, minX, maxX, minY, maxY):
gdal.Translate(IMG_OUT, IMG_IN, projWin=[minX,maxY,maxX,minY])
return print("Extraction process was completed! : ", IMG_IN)②でダウンロードしたファイルを指定して、切り出した上で、赤青合成を行い、PNGで表示させます。
if __name__ == '__main__':
# 処理対象ファイルの指定(02_download_COG.pyの実行結果)
IMG_PATH_1 = "./2c9ccde2-3a8c-4cd2-aac5-4f0c3b33a4a4/AS200355208498-180908_webcog.tif" #2018-09-08
IMG_PATH_2 = "./7d47301d-9533-4be2-a010-31a86912e865/AS201462808498-200905_webcog.tif" #2020-09-05
# AOI抽出後のGeotiffの保存先の指定
IMG_PATH_1_CUT = "./out/AS200355208498-180908_webcog_cut.tif"
IMG_PATH_2_CUT = "./out/AS201462808498-200905_webcog_cut.tif"
IMG_JPG = "./out/out.jpeg"
# AOIの四隅緯度経度指定
minX = 141.93507323883878
minY = 42.75282111399301
maxX = 141.94639048725506
maxY = 42.75848804932892
# AOIの抽出処理（四隅）
cut4_geotiff(IMG_PATH_1, IMG_PATH_1_CUT, minX, maxX, minY, maxY)
cut4_geotiff(IMG_PATH_2, IMG_PATH_2_CUT, minX, maxX, minY, maxY)
# IMG_1読込
ds_1 = gdal.Open(IMG_PATH_1_CUT, gdal.GA_ReadOnly)
image_1 = np.array([ds_1.GetRasterBand(i + 1).ReadAsArray() for i in range(ds_1.RasterCount)])
del ds_1 #不要メモリ解放
# IMG_2読込
ds_2 = gdal.Open(IMG_PATH_2_CUT, gdal.GA_ReadOnly)
image_2 = np.array([ds_2.GetRasterBand(i + 1).ReadAsArray() for i in range(ds_2.RasterCount)])
del ds_2 #不要メモリ解放
# 画像間のピクセル数を揃える処理(今回は必要なし)
x_min = min([image_1.shape[1], image_2.shape[1]])
y_min = min([image_1.shape[2], image_2.shape[2]])
# print("MinX : ", image_1.shape[1], image_2.shape[1])
# print("MinY : ", image_1.shape[2], image_2.shape[2])
# 出力用配列の準備
image = np.zeros((x_min,y_min,3)).astype(np.uint8)
image[:,:,0] = image_2[0, :x_min, :y_min]
image[:,:,1] = image_1[0, :x_min, :y_min]
image[:,:,2] = image_1[0, :x_min, :y_min]
# print(image_1.shape, image_2.shape, image.shape)
# 不要メモリの解放
del image_1
del image_2
# 出力画像の保存
pil_img = Image.fromarray(image)
# print(pil_img.mode)
pil_img.save(IMG_JPG)
print("DONE!")出力された画像は以下の通りです。
2年という期間があいてしまっているため、山の木の状態などが異なり、全体的に青みがかっていますが、画像中央に赤く見えるポイントがあり、発災から2年が経過し、土砂崩れの被害があった場所から土砂が撤去され、建屋が増強されているようにも見えます。
参考：
https://www.hokkaido-np.co.jp/article_photo/list?article_id=456528&p=6658530&rct=n_hokkaido
まとめ
本記事では、ASNARO-2を題材に、Tellus TravelerのAPIを使ってSAR画像の前後比較を行うコードをご紹介しました。
ASNARO-2の画像はTellus Traveler上に他にもありますので、ぜひ探してみてください。また、もし希望の場所や時期に画像がない場合、購入することも可能です。Tellusまで、お問い合わせください。
