防災分野での衛星データ利用事例

この章で学習すること

  • TelluSARを用いた差分干渉解析
  • 差分干渉解析結果に基づく地殻変動の推定
  • 2時期のSAR強度画像を用いたRGB疑似カラー合成画像による浸水範囲の推定(GDALを用いて実装、基礎分野を参照)
  • 2時期のSAR強度画像を用いたレベルスライス法による浸水範囲の推定(GDALとOGRを用いて実装、中級者向け)

背景

日本は、地震、台風、豪雨などといった自然災害が発生しやすい国です。近年、こういった災害による被害が広域化・同時多発化する傾向が見られ、今後も、南海トラフ地震のような巨大地震や豪雨災害など甚大な災害の発生が懸念されています。災害の初動対応では、被害状況をいち早く把握することが必要不可欠ですが、被害が広域化・同時多発化することで対応がより難しくなります。このような中、衛星画像を用いることでより早く、より効率的に概況調査を実施できると期待されています。特に、この章で扱う合成開口レーダ(SAR)は全天候型(観測時の天候に左右されにくい)であり、夜間でも観測可能であるため、災害発生時に迅速な対応ができます。

SARでは、受信信号の強度と位相を記録しています。各観測での位相差を計算し、地表面の変位を計測するのがこの章の前半で扱う差分干渉解析です。差分干渉解析によって、地震や火山活動に伴う地殻変動や地すべり、工事に伴う地盤沈下といった地表面変位の状況を把握することができます。このうち、この章では地震による地殻変動を対象とします。

一方、複数時期の観測データを利用して受信信号の強度変化を調べることで、地表面の様子が著しく変化した箇所を把握することができます。これがこの章の後半で扱う内容で、土砂災害や浸水箇所の概況調査が可能です。

使用する衛星データとプロダクトレベル

この章ではJAXAの陸域観測技術衛星2号「だいち2号」(ALOS-2)に搭載された合成開口レーダPALSAR-2のプロダクトのうち、Level 1.1(L1.1)とLevel 2.1(L2.1)を使用します。L1.1は複素データ(実部、虚部のデータを格納している)のため、受信信号の強度に加え、位相情報も保持しています。差分干渉解析には位相情報が必要であるため、この章の前半ではL1.1を使用します。

L1.1から複素数の大きさを計算して実画像とした上で、実際の地表面上の位置に投影されたものがL1.5、そこからさらに地形による倒れこみを補正(オルソ補正)したものがL2.1です。この章の後半では受信信号の強度を利用した解析を行うため、そのまま地図投影ができるL2.1を使用します。ただし、複素データではないので位相情報が失われている、つまり、L1.5やL2.1を差分干渉解析に用いることはできないのでご注意ください。

TelluSARの差分干渉解析機能を用いた地殻変動状況の確認

TelluSARを用いた干渉SARの活用方法を学びます。干渉SARは地表面の変化を調べる解析手法で、主に地殻変動や地盤沈下、地すべりなどによる地表面変位の状況を把握するために使われています。干渉SARを初めて聞いたという方は、宙畑に干渉SARの概要や解析事例を紹介した記事があるので、そちらを参照してください。⇒「干渉SAR(InSAR)とは-分かること、事例、仕組み、読み解き方-

さて、実際に干渉SARを行ってみましょう。使用する衛星データはALOS-2/PALSAR-2のL1.1(位相情報を含む)です。Tellusに搭載されているPALSAR-2オリジナルデータを用いて1つずつ処理を行う方法は、宙畑の「【コード付き】地盤の沈降が分かる干渉SAR解析をTellusでやってみた」で説明されています。が、複雑な処理なので、衛星データに慣れてない方にはハードルが高いかもしれません。そこで、TelluSARの差分干渉解析機能を使います。TelluSARではシーン(画像)を選んでしまえば自動で処理を実行してくれるため、誰でも簡単に干渉SAR解析を行うことができます。TelluSARの基本的な使用方法は宙畑の「ついに登場!公式解析ツール「TelluSAR(てるーさー)」 APIの使い方を徹底解説!」で紹介されているので参考にしつつ、ここではより実践的に、平成30年北海道胆振東部地震を対象にTelluSARの差分干渉解析機能を用いて差分干渉画像を作成し、地殻変動の状況を確認します。

平成30年北海道胆振東部地震

平成30年9月6日未明に、北海道胆振東部地方中東部を震源とする地震が発生しました。最大震度は厚真町で観測された震度7で、厚真町を中心に斜面崩壊が多発、札幌市清田区など市街地での液状化現象や北海道全域での大規模停電が発生しました。

事前準備

それでは、使用するシーン(画像)を選んで解析を始めましょう。利用できるシーンはTelluSARの検索機能で確認できますが、まずは、TelluSAR(API)を使うための事前準備をします。

必要なモジュール、ライブラリをインポートします。

In [ ]:
import requests, datetime, json, pprint
from skimage import io
from io import BytesIO

マーケットトークンを発行します。 ※get_market_tokenを使用しています。

In [ ]:
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
    }
    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})
market_token = ret['token']
print(market_token)

差分干渉解析に用いるシーンの選定

差分干渉解析では2時期に撮影されたSAR画像のペアを用いて、その間に生じた地表面変位を解析するので、基準となるシーン(primary画像)とその後に観測されたシーン(secondary画像)が必要となります。今回は地震の影響を把握することが目的なので、primary画像は地震前から、secondary画像は地震後から選ぶことになります。TelluSARで利用できるシーンを検索する際にパラメータとして位置情報を与えておくと、その範囲を含むシーンを確認することができるので、震源の位置を与えておきます。

In [ ]:
# 解析に用いるシーンが次の範囲を含む
# 平成30年北海道胆振東部地震の震源のおおよその位置
parameter1 = {
    'left_bottom_lon' : 142.000,
    'left_bottom_lat' : 42.690,
    'right_top_lon' : 142.000,
    'right_top_lat' : 42.690
}

解析に用いるシーンの位置条件を設定したら、その範囲を含むシーンの一覧を確認します。一覧ではシーンの観測範囲や偏波の情報も得られますが、情報量が膨大になるので、ここでは観測時刻とシーンIDのみを表示します。
※get_free_sceneを使用しています。

In [ ]:
# 利用できるシーン一覧を取得する
def get_free_scene(market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v1/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(market_token, parameter1)

# シーンの枚数を表示
print(len(get_free_scene_result['data']['scenes']))

# シーンの観測時刻とシーンIDのみを表示
scene_list = []
for dict in get_free_scene_result['data']['scenes']:
    newdict = {key:dict[key] for key in ['scene_id','observation_datetime']}
    scene_list.append(newdict.copy())
pprint.pprint(scene_list)
73
[{'observation_datetime': '2016-09-29T14:25:38+00:00',
  'scene_id': 'ALOS2127080840-160929'},
 {'observation_datetime': '2018-03-08T13:37:37+00:00',
  'scene_id': 'ALOS2204700890-180308'},
 {'observation_datetime': '2018-03-22T13:37:36+00:00',
  'scene_id': 'ALOS2206770890-180322'},
 {'observation_datetime': '2018-03-29T14:25:37+00:00',
  'scene_id': 'ALOS2207810840-180329'},
 {'observation_datetime': '2018-04-05T13:37:36+00:00',
  'scene_id': 'ALOS2208840890-180405'},
 {'observation_datetime': '2018-04-19T13:37:35+00:00',
  'scene_id': 'ALOS2210910890-180419'},
 {'observation_datetime': '2018-05-03T02:40:57+00:00',
  'scene_id': 'ALOS2212912750-180503'},
 {'observation_datetime': '2018-05-03T13:37:35+00:00',
  'scene_id': 'ALOS2212980890-180503'},
 {'observation_datetime': '2018-05-17T02:40:56+00:00',
  'scene_id': 'ALOS2214982750-180517'},
 {'observation_datetime': '2018-05-17T13:37:34+00:00',
  'scene_id': 'ALOS2215050890-180517'},
 {'observation_datetime': '2018-05-29T14:32:33+00:00',
  'scene_id': 'ALOS2216830840-180529'},
 {'observation_datetime': '2018-05-31T13:37:35+00:00',
  'scene_id': 'ALOS2217120890-180531'},
 {'observation_datetime': '2018-06-14T02:40:57+00:00',
  'scene_id': 'ALOS2219122750-180614'},
 {'observation_datetime': '2018-06-14T02:41:05+00:00',
  'scene_id': 'ALOS2219122760-180614'},
 {'observation_datetime': '2018-06-14T13:37:35+00:00',
  'scene_id': 'ALOS2219190890-180614'},
 {'observation_datetime': '2018-06-28T13:37:35+00:00',
  'scene_id': 'ALOS2221260890-180628'},
 {'observation_datetime': '2018-07-07T02:34:08+00:00',
  'scene_id': 'ALOS2222522760-180707'},
 {'observation_datetime': '2018-07-12T13:37:35+00:00',
  'scene_id': 'ALOS2223330890-180712'},
 {'observation_datetime': '2018-07-15T14:39:23+00:00',
  'scene_id': 'ALOS2223780830-180715'},
 {'observation_datetime': '2018-07-15T14:39:31+00:00',
  'scene_id': 'ALOS2223780840-180715'},
 {'observation_datetime': '2018-07-26T02:40:57+00:00',
  'scene_id': 'ALOS2225332750-180726'},
 {'observation_datetime': '2018-07-26T13:37:35+00:00',
  'scene_id': 'ALOS2225400890-180726'},
 {'observation_datetime': '2018-08-09T13:37:35+00:00',
  'scene_id': 'ALOS2227470890-180809'},
 {'observation_datetime': '2018-08-23T02:40:57+00:00',
  'scene_id': 'ALOS2229472750-180823'},
 {'observation_datetime': '2018-08-23T02:41:05+00:00',
  'scene_id': 'ALOS2229472760-180823'},
 {'observation_datetime': '2018-08-23T13:37:35+00:00',
  'scene_id': 'ALOS2229540890-180823'},
 {'observation_datetime': '2018-08-25T14:18:47+00:00',
  'scene_id': 'ALOS2229840850-180825'},
 {'observation_datetime': '2018-09-06T02:40:57+00:00',
  'scene_id': 'ALOS2231542750-180906'},
 {'observation_datetime': '2018-09-06T02:41:05+00:00',
  'scene_id': 'ALOS2231542760-180906'},
 {'observation_datetime': '2018-09-06T13:37:35+00:00',
  'scene_id': 'ALOS2231610890-180906'},
 {'observation_datetime': '2018-09-08T14:18:47+00:00',
  'scene_id': 'ALOS2231910850-180908'},
 {'observation_datetime': '2018-09-13T14:25:37+00:00',
  'scene_id': 'ALOS2232650840-180913'},
 {'observation_datetime': '2018-09-15T02:34:08+00:00',
  'scene_id': 'ALOS2232872760-180915'},
 {'observation_datetime': '2018-09-20T02:40:57+00:00',
  'scene_id': 'ALOS2233612750-180920'},
 {'observation_datetime': '2018-09-20T02:41:05+00:00',
  'scene_id': 'ALOS2233612760-180920'},
 {'observation_datetime': '2018-09-20T13:37:35+00:00',
  'scene_id': 'ALOS2233680890-180920'},
 {'observation_datetime': '2018-10-04T02:40:57+00:00',
  'scene_id': 'ALOS2235682750-181004'},
 {'observation_datetime': '2018-10-04T02:41:06+00:00',
  'scene_id': 'ALOS2235682760-181004'},
 {'observation_datetime': '2018-10-04T13:37:35+00:00',
  'scene_id': 'ALOS2235750890-181004'},
 {'observation_datetime': '2018-10-06T14:18:47+00:00',
  'scene_id': 'ALOS2236050850-181006'},
 {'observation_datetime': '2018-10-09T13:44:25+00:00',
  'scene_id': 'ALOS2236490880-181009'},
 {'observation_datetime': '2018-10-11T14:25:36+00:00',
  'scene_id': 'ALOS2236790840-181011'},
 {'observation_datetime': '2018-10-18T02:40:57+00:00',
  'scene_id': 'ALOS2237752750-181018'},
 {'observation_datetime': '2018-10-18T02:41:06+00:00',
  'scene_id': 'ALOS2237752760-181018'},
 {'observation_datetime': '2018-11-01T02:40:58+00:00',
  'scene_id': 'ALOS2239822750-181101'},
 {'observation_datetime': '2018-11-01T02:41:06+00:00',
  'scene_id': 'ALOS2239822760-181101'},
 {'observation_datetime': '2018-11-08T14:25:37+00:00',
  'scene_id': 'ALOS2240930840-181108'},
 {'observation_datetime': '2018-11-15T02:40:58+00:00',
  'scene_id': 'ALOS2241892750-181115'},
 {'observation_datetime': '2018-11-15T02:41:06+00:00',
  'scene_id': 'ALOS2241892760-181115'},
 {'observation_datetime': '2018-11-15T13:37:28+00:00',
  'scene_id': 'ALOS2241960880-181115'},
 {'observation_datetime': '2018-11-15T13:37:36+00:00',
  'scene_id': 'ALOS2241960890-181115'},
 {'observation_datetime': '2018-11-29T02:40:58+00:00',
  'scene_id': 'ALOS2243962750-181129'},
 {'observation_datetime': '2018-11-29T02:41:06+00:00',
  'scene_id': 'ALOS2243962760-181129'},
 {'observation_datetime': '2018-11-29T13:37:28+00:00',
  'scene_id': 'ALOS2244030880-181129'},
 {'observation_datetime': '2018-11-29T13:37:36+00:00',
  'scene_id': 'ALOS2244030890-181129'},
 {'observation_datetime': '2018-12-08T02:34:10+00:00',
  'scene_id': 'ALOS2245292760-181208'},
 {'observation_datetime': '2018-12-13T02:40:59+00:00',
  'scene_id': 'ALOS2246032750-181213'},
 {'observation_datetime': '2018-12-13T02:41:07+00:00',
  'scene_id': 'ALOS2246032760-181213'},
 {'observation_datetime': '2018-12-27T13:37:37+00:00',
  'scene_id': 'ALOS2248170890-181227'},
 {'observation_datetime': '2019-01-10T13:37:37+00:00',
  'scene_id': 'ALOS2250240890-190110'},
 {'observation_datetime': '2019-01-24T13:37:38+00:00',
  'scene_id': 'ALOS2252310890-190124'},
 {'observation_datetime': '2019-02-07T13:37:37+00:00',
  'scene_id': 'ALOS2254380890-190207'},
 {'observation_datetime': '2019-02-09T14:18:49+00:00',
  'scene_id': 'ALOS2254680850-190209'},
 {'observation_datetime': '2019-02-21T13:37:38+00:00',
  'scene_id': 'ALOS2256450890-190221'},
 {'observation_datetime': '2019-03-21T13:37:36+00:00',
  'scene_id': 'ALOS2260590890-190321'},
 {'observation_datetime': '2019-04-04T13:37:37+00:00',
  'scene_id': 'ALOS2262660890-190404'},
 {'observation_datetime': '2019-04-11T14:25:37+00:00',
  'scene_id': 'ALOS2263700840-190411'},
 {'observation_datetime': '2019-04-18T13:37:35+00:00',
  'scene_id': 'ALOS2264730890-190418'},
 {'observation_datetime': '2019-05-02T02:40:57+00:00',
  'scene_id': 'ALOS2266732750-190502'},
 {'observation_datetime': '2019-05-02T02:41:05+00:00',
  'scene_id': 'ALOS2266732760-190502'},
 {'observation_datetime': '2019-05-09T14:25:36+00:00',
  'scene_id': 'ALOS2267840840-190509'},
 {'observation_datetime': '2019-06-06T14:25:36+00:00',
  'scene_id': 'ALOS2271980840-190606'},
 {'observation_datetime': '2019-06-15T03:22:17+00:00',
  'scene_id': 'ALOS2273242720-190615'}]

震源の位置を含む観測シーンが73枚あることがわかります。まずは、primary画像を選んでみましょう。地震は2018年9月6日に発生しています。一覧の'observation_datetime'を見ていくと、地震の直前では、2018年8月23日に観測されたシーンがあるようなので、そのシーンを選ぶことにします。(時刻はUTC表記なので、日本時刻が知りたい場合は+9時間してください。)

In [ ]:
# 2018年8月23日のシーンID
scene_id = 'ALOS2229472760-180823'

primary画像を選んだので、次はsecondary画像を選びます。secondary画像を選ぶ際に必要な条件は、以下の通りです。条件を満たしていないシーンを選ぶと差分干渉解析を実施することができません。

  • シーンフレーム:同一であること
  • 観測モード:同一であること
  • 偏波:同一であること
  • 観測角(オフナディア角):同一であること
  • 衛星進行方向(アセンディング/ディセンディング):同一であること
  • レーダ照射方向:同一であること

TelluSARでは基準となるシーンの'scene_id'を指定するとこれらの条件を満たすsecondary画像の候補を表示してくれます。
※get_scene_afterを使用しています。

In [ ]:
# ペアの候補となるシーン一覧を取得する
# after_scene_id
def get_scene_after(scene_id, market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v1/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, market_token, {})
print(len(get_scene_after_result['data']['scenes']))

scene_after_list = []
for dict in get_scene_after_result['data']['scenes']:
    newdict = {key:dict[key] for key in ['scene_id','observation_datetime']}
    scene_after_list.append(newdict.copy())
pprint.pprint(scene_after_list)
9
[{'observation_datetime': '2018-09-06T02:41:05+00:00',
  'scene_id': 'ALOS2231542760-180906'},
 {'observation_datetime': '2018-09-20T02:41:05+00:00',
  'scene_id': 'ALOS2233612760-180920'},
 {'observation_datetime': '2018-10-04T02:41:06+00:00',
  'scene_id': 'ALOS2235682760-181004'},
 {'observation_datetime': '2018-10-18T02:41:06+00:00',
  'scene_id': 'ALOS2237752760-181018'},
 {'observation_datetime': '2018-11-01T02:41:06+00:00',
  'scene_id': 'ALOS2239822760-181101'},
 {'observation_datetime': '2018-11-15T02:41:06+00:00',
  'scene_id': 'ALOS2241892760-181115'},
 {'observation_datetime': '2018-11-29T02:41:06+00:00',
  'scene_id': 'ALOS2243962760-181129'},
 {'observation_datetime': '2018-12-13T02:41:07+00:00',
  'scene_id': 'ALOS2246032760-181213'},
 {'observation_datetime': '2019-05-02T02:41:05+00:00',
  'scene_id': 'ALOS2266732760-190502'}]

9枚のシーンが候補として表示されました。この中からsecondary画像を選ぶ時のポイントは「観測日」です。地表面の被覆状態の変化は干渉性の良し悪し(結果の良し悪し)に影響するため、primary画像に対して、観測日が近い、または、同一季節のシーンを用いることが望ましいです。地震が発生した日の正午頃に観測されたシーンがあるため、そのシーンを選ぶことにします。

今回は「観測日」に注目してシーンを選定しましたが、他にも、衛星の観測条件を考慮することがあります。差分干渉解析では、衛星が軌道上の同一の箇所に来た時に観測したシーンのペアを用いていますが、衛星の位置が完全に一致することはありません。このときの軌道間距離の衛星視線方向に直交する成分を垂直基線長(Bperp)と呼び、このBperpが大きいと干渉性が低下します(結果が悪くなります)。BperpはTellus上で確認することはできませんが、JAXAが運用しているAUIG2(ALOS及びALOS-2の観測データに関する各種オンラインサービスを提供するシステム)で調べることができます。

In [ ]:
# 2枚目のscene_id
# 2018年9月6日(地震直後に観測)のシーンID
scene_after_id = 'ALOS2231542760-180906'

差分干渉解析処理の実行

用いるシーンのペアを選ぶことができたので、差分干渉解析を実施し、差分干渉画像を作成します。先ほど紹介した宙畑の記事から、差分干渉解析のフローを以下に示しますが、TelluSARではこれらの処理を自動で実施するので、誰でも簡単に差分干渉画像を得ることができます。


干渉SARの処理フロー( https://sorabatake.jp/12465/

処理パラメータとしてルック数とフィルター回数を指定できます。SAR画像にはスペックルノイズと言われるごま塩状ノイズが生じますが、このノイズを低減する方法としてルック処理というものがあります。ルック処理はSARに特有な処理方法であり、厳密には少し違いますが、複数のピクセルの平均値を代表値としてダウンサンプリングする方法だと思ってください。平均値を計算する画素数を指定することができ、'nlook_rg'がレンジ方向(電波照射方向)、'nlook_az'がアジマス方向(衛星進行方向)のルック数を示しています。また、生成された干渉縞について、移動平均フィルターもノイズ低減手法として有効であり、'filter'でその回数を指定することができます。


ルック処理のイメージ(厳密にはもう少し複雑な処理をしています。)


ルック数の違いによる解像度の変化(ルック数のより大きい右側の画像の方が解像度が低く滑らかになります。Original data provided by JAXA)

できるだけノイズの影響を低減してきれいな差分干渉画像を得るため、ルック数とフィルター回数は指定できる最大値に設定しています。今回の解析対象は地震による地殻変動であり、空間スケールが比較的大きい事象です。この場合は、ノイズ低減処理(ルック処理・フィルター)を強めにかけても問題ありませんが、注目している事象が局所的な事象である場合、これらの処理を強めにかけると求めたい変位が埋もれて見えなくなってしまう可能性があります。そのような場合はルック数を下げるか、フィルター回数を下げて調整してください。

In [ ]:
# 処理パラメータ
# ノイズをできるだけ低減するため、ルック数とフィルター回数を最大に設定
parameter2 = {
    'before_scene_id': scene_id,
    'after_scene_id': scene_after_id,
    'polarisation': 'HH',
    'nlook_rg': 5,
    'nlook_az': 7,
    'filter': 2,
}

処理パラメータを設定したら、処理を実行します。
※request_workを使用しています。

In [ ]:
# 処理を実行する
# work_idが返却される
def request_work(market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v1/works'
    headers = {
        'Authorization': 'Bearer ' + market_token
    }
    r = requests.post(url, headers=headers, params=payload)
    print(r.url)
    if r.status_code is not 200:
        print(r.status_code)
    return json.loads(r.content)

work_result = request_work(market_token, parameter2)
print(work_result)
https://tellusar.tellusxdp.com/api/v1/works?before_scene_id=ALOS2229472760-180823&after_scene_id=ALOS2231542760-180906&polarisation=HH&nlook_rg=5&nlook_az=7&filter=2
{'data': {'work_id': 272, 'exist_flag': True}}

処理を開始すると処理Id(work_id)が返されます。既に処理された結果が存在する場合は'exist_flag': Trueとなります。処理時間は30分~数時間程度かかります。

処理ステータスの確認は、work_idを指定する方法と一覧で確認する方法があります。
※get_work, get_work_listを使用しています。

In [ ]:
# work_idを指定して処理結果を取得する
def get_work(work_id, market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v1/works/{}'.format(work_id)
    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) 

work_id = work_result['data']['work_id']
get_work_result = get_work(work_id, market_token, {})
print(get_work_result)
https://tellusar.tellusxdp.com/api/v1/works/272
{'data': {'work_id': 272, 'before_scene_id': 'ALOS2229472760-180823', 'after_scene_id': 'ALOS2231542760-180906', 'polarisation': 'HH', 'nlook_rg': 5, 'nlook_az': 7, 'filter': 2, 'progress_status': 2, 'request_date': '2020-11-15T08:36:02+00:00', 'complete_date': '2020-11-15T10:36:16+00:00'}}
In [ ]:
# 処理結果の一覧を取得する
def get_work_list(market_token, payload={}):
    url = 'https://tellusar.tellusxdp.com/api/v1/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(market_token, {})
print(get_work_result_list)
https://tellusar.tellusxdp.com/api/v1/works
{'data': {'works': [{'work_id': 275, 'before_scene_id': 'ALOS2210180680-180414', 'after_scene_id': 'ALOS2222600680-180707', 'polarisation': 'HH', 'nlook_rg': 5, 'nlook_az': 7, 'filter': 0, 'progress_status': 2, 'request_date': '2020-11-15T14:11:33+00:00', 'complete_date': '2020-11-15T14:11:33+00:00'}, {'work_id': 272, 'before_scene_id': 'ALOS2229472760-180823', 'after_scene_id': 'ALOS2231542760-180906', 'polarisation': 'HH', 'nlook_rg': 5, 'nlook_az': 7, 'filter': 2, 'progress_status': 2, 'request_date': '2020-11-15T08:36:02+00:00', 'complete_date': '2020-11-15T10:36:16+00:00'}, {'work_id': 264, 'before_scene_id': 'ALOS2204710630-180308', 'after_scene_id': 'ALOS2210920630-180419', 'polarisation': 'HH', 'nlook_rg': 5, 'nlook_az': 7, 'filter': 0, 'progress_status': 2, 'request_date': '2020-11-12T10:20:54+00:00', 'complete_date': '2020-11-12T10:20:54+00:00'}]}}

差分干渉解析結果の表示

作成された差分干渉画像を確認してみます。タイル画像(png)を取得して表示する方法とGeoTIFFを表示する方法がありますが、タイル画像の表示は宙畑で紹介されているので、ここではGeoTIFFを表示します。

In [ ]:
url = 'https://tellusar.tellusxdp.com/api/v1/works/{}/tifs/fringe_diff.tif'.format(work_id)
headers = {
    'Authorization': 'Bearer ' + market_token
}
r = requests.get(url, headers=headers)
print(r.url)
if r.status_code is not 200:
    print(r.status_code)
    
fringe=io.imread(BytesIO(r.content))
io.imshow(fringe)
https://tellusar.tellusxdp.com/api/v1/works/272/tifs/fringe_diff.tif
Out[ ]:
<matplotlib.image.AxesImage at 0x7f9497050f60>

差分干渉画像(Original data provided by JAXA)

差分干渉画像(GeoTIFF)をダウンロードします。位相に応じた色付けがされており、地理情報を持っているので、GIS上にそのまま表示することができます。

In [ ]:
with open('fringe.tif', 'wb') as f:
    for data in r.iter_content(1024):
        f.write(data)

ダウンロードした差分干渉画像をGIS上に表示してみました。

画像の下2/3は海なので差分干渉解析結果が得られず砂嵐状になっています。上1/3は陸域で干渉性が比較的良く、震源の南東に目玉状の変状が見られています。

差分干渉解析結果の解釈

差分干渉解析結果から地表面変位量を読み解きます。変位を読む際は衛星側から読むのが一般的ですが、変状箇所が目玉状に集まっている場合はその外側から読んでいくこともあります。目玉状の変状の外側、陸域の大部分は水色を呈しています。目玉の外(A)から目玉の中心(B)に向かうにつれて、水色→黄色→ピンクのように色変化しています。この色変化を凡例に照らし合わせると、黄色→ピンクの色変化は-πを超えてしまいますが、これは、得られる位相が-π~+πに折り畳まれているためなので、凡例をもう1周分追加して読み取ってください。おおよそπ分強、衛星に近づいていると読み取れます。この位相量は利用するレーダ波の波長に対応しており、波長の半分が2πに相当します(衛星と地表面間の往復分なので、波長の半分となります)。ALOS-2の場合、波長は約23.6cmなので、2π=11.8cmとして換算でき、地点Aに対して地点Bは6cm強、衛星に近づいていると推定できます。このように、差分干渉解析で得られる地表面変位は相対変位であるので、変位量を読み取る際は基準となる点を設定する必要があります。

差分干渉解析で得られた変位量は衛星視線方向の変位量です。SARによる観測は、衛星から斜め下方向に電波を照射しているため、得られる変位量は上下成分と東西成分が合成されたものです。ここで、今回用いたシーンの観測条件を確認します。
※search_fileに変更を加えています。

In [ ]:
def search_file_id(scene_id, params={}, next_url=''):
    if len(next_url) > 0:
        url = next_url
    else:
        url = 'https://file.tellusxdp.com/api/v1/origin/search/palsar2-l11/{}'.format(scene_id)
        
    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()

ret = search_file_id(scene_id, {})

pprint(ret)
{'bbox': [141.401, 42.107, 142.23, 42.807],
 'beam_number': 'U2-7',
 'begin_datetime': '2018-08-23T02:41:00.518000+00:00',
 'center_datetime': '2018-08-23T02:41:05.518000+00:00',
 'code': 'UBS',
 'coordinates': [[141.401, 42.197],
                 [142.077, 42.107],
                 [142.23, 42.716],
                 [141.547, 42.807]],
 'dataset_id': 'ALOS2229472760-180823',
 'date_added': '2019-12-18T09:34:25.186298',
 'end_datetime': '2018-08-23T02:41:10.518000+00:00',
 'files': ['1', '2', '3', '4', '5', '6'],
 'frame_number': 2760,
 'mode': 'SM1',
 'observation_datetime': '2018-08-23T02:41:05.518000+00:00',
 'observation_direction': 'R',
 'off_nadir_angle': 32.8,
 'orbit_direction': 'D',
 'path_number': 18,
 'polarisations': ['HH'],
 'publish_link': 'https://file.tellusxdp.com/api/v1/origin/publish/palsar2-l11/ALOS2229472760-180823'}

衛星進行方向はディセンディング('orbit_direction': 'D')で、右側観測('observation_direction': 'R')なので、得られた差分干渉画像と衛星の位置関係は次の通りです。

得られた結果は東側からの観測結果であり、衛星に近づく変位が出ていた場合、隆起・東向き変位が生じている可能性があります。ここから隆起量あるいは東向き変位を求めたい場合はどちらか一方の変位がないと仮定をおく必要があります。このほか、反対側からの観測(アセンディング・右側観測)シーンのペアについて差分干渉解析を実施することで西側からの観測結果を得ることができ、東側からの観測結果と合わせることで、準上下方向と準東西方向に分離することができます。これを2.5次元解析と呼びます。同時性が保てないことが課題ではありますが、胆振東部地震の場合は運良く、primaryとsecondaryともに、ディセンディングの観測(今回用いたペア)の半日後の観測があったため、それらのシーン(アセンディングのペア)を用いた解析を合わせると、今回得られた変状は隆起を示しており、震源の南東側が6cm強の隆起を示していると判断できたようです。(国土地理院の解析例

コヒーレンス

TelluSARの出力結果には、強度画像(primary及びsecondary)、差分干渉画像に加えて、コヒーレンス画像があります。コヒーレンス(干渉性)とは、位相情報の推定の良し悪しの指標であり、コヒーレンスが高いほど信頼性が高いと考えます。

コヒーレンス画像を見てみると、陸域が比較的高い値を示し(明るい)、海では低い値(暗い)となることがわかります。差分干渉画像で海は砂嵐状になっていましたが、このように位相がランダムとなる(位相の推定ができていない、信頼性が低い)箇所ではコヒーレンスが低い値となります。逆に言えば、コヒーレンスが低い箇所の位相はランダムなので地表面変位を読み解くことが難しく、差分干渉画像からコヒーレンスの低い箇所をマスクしておくと解釈しやすくなることがあります。


コヒーレンス画像(Original data provided by JAXA)

コヒーレンスが低くなる要因として地表面の状態の著しい変化があることが考えられます。この特徴を利用した被害箇所の解析方法として、この章の後半で扱う2時期のSAR強度画像を用いたRGB疑似カラー合成画像にコヒーレンスの情報を加えて変位箇所を可視化するMulti Temporal Coherence mapping(MTC)やコヒーレンスの変化に着目した3時期コヒーレンス解析といったものがあります。

コラム1:"Master"と"Slave"

これまでInSAR界隈では、基準となる画像を"master"、2枚目の画像を"slave"と呼ぶことが主流でした。おそらく、皆さんがInSARの解析事例をネット等で調べると、この用語を使っているものが多く出てくるのではないでしょうか。

昨今、"Black Lives Matter"が話題となっていますが、この流れがInSARの分野にまで波及し、"master" を "primary" に、 "slave" を "secondary" に言い換えるようにしましょうという動きが出てきています。本稿ではこの考えに賛同して、"primary"と"secondary"という用語で説明するようにしています。

二時期のPALSAR-2データ(強度画像)を用いた浸水範囲の推定


ここからはSARの強度画像を用いた解析方法を学びます。強度画像は地表面上の散乱体の有無や凹凸を反映したものであり、散乱体が有る場合や凹凸が粗い場合により強い後方散乱が得られるため、大きい値をとります。土砂災害や浸水の被害範囲では地表面の様子が著しく変化するため強度画像に変化が生じ、災害発生時を挿む二時期の強度画像の差分をとることで被害範囲を可視化することが期待できます。

Tellusには、ALOS-2/PALSAR-2(L2.1)のオリジナルデータが搭載されているので、今回はそのデータを使います。オリジナルデータの取得方法は宙畑の「【ゼロからのTellusの使い方】提供元オリジナルデータを取得する」にて解説されているので参考にしつつ、平成30年7月豪雨時の岡山県倉敷市真備町の浸水を対象に、二時期の強度画像を用いた浸水範囲の推定を行います。

平成30年7月豪雨

平成30年6月28日から7月8日にかけて、梅雨前線や台風の影響により、西日本を中心に全国的に広い範囲で豪雨が発生しました。西日本から東海地方の多くの地点で24、48、72時間雨量の観測史上最大値が更新され、西日本の広域にわたって土砂災害や洪水、浸水被害が発生しました。

関心領域(AOI)の設定

Tellus上で、岡山県倉敷市真備町周辺のポリゴン(GeoJSON)を用意します。ポリゴンの作成方法は漁業分野を参照してください。

作成したポリゴン(GeoJSON)をTellus開発環境にアップロードした後、読み込んで範囲の情報を確認します。

In [ ]:
import geopandas as gpd
aoi = gpd.read_file('aoi.geojson', encoding='utf-8')
bounds = aoi['geometry'].bounds
bounds
Out[ ]:
minx miny maxx maxy
0 133.640099 34.610606 133.742924 34.667806

使用するシーンの選定

使用するシーンの条件を設定します。真備町の浸水被害は2018年7月7日に発生しているので、第一に、この直後に観測されたシーンが必要となります。次に、今回は二時期のデータを用いて浸水被害を推定するため、浸水前のシーンも必要となります。差分干渉解析のときと同様に、ここでも観測日が重要となります。時期が離れると人工改変等の影響が出てしまう、あるいは、季節変化の影響が出ることもあるので、できるだけ期間の短い、あるいは、季節を揃えたシーンを選択することが望ましいです。

上記の条件に合うシーンを検索します。期間を1年間として、関心領域を含むようにしたときのパラメータは次のように指定できます。

In [ ]:
#使用するシーンの条件
parameter2 = {
    'after' : '2017-06-01T00:00:00',
    'before' : '2018-07-10T00:00:00',
    'left_bottom_lon' : bounds['minx'],
    'left_bottom_lat' : bounds['miny'],
    'right_top_lon' : bounds['maxx'],
    'right_top_lat' : bounds['maxy']
}

使用するシーンの条件を設定したら、Tellusに搭載されているPALSAR-2オリジナルデータを検索します。受信信号の位相情報は必要ないので、地図投影されたL2.1データを使用します。
※search_fileを使用しています。

In [ ]:
import requests
from pprint import pprint

TOKEN = '自分のトークンをはる'
provider_id = 'tellus-product'
tool_label = 'tellusar-api-f'

def search_file(params={}, next_url=''):
    if len(next_url) > 0:
        url = next_url
    else:
        url = 'https://file.tellusxdp.com/api/v1/origin/search/palsar2-l21'
        
    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()

ret = search_file(parameter2)

pprint(ret['count'])
pprint(ret['items'])
2
[{'bbox': [133.418, 34.35, 134.153, 35.045],
  'beam_number': 'U2-9',
  'begin_datetime': '2018-04-14T15:05:06.387000+00:00',
  'center_datetime': '2018-04-14T15:05:11.387000+00:00',
  'code': 'UBS',
  'coordinates': [[133.547, 34.35],
                  [134.153, 34.435],
                  [134.029, 35.045],
                  [133.418, 34.961]],
  'dataset_id': 'ALOS2210180680-180414',
  'date_added': '2020-11-11T16:35:57.151913',
  'end_datetime': '2018-04-14T15:05:16.387000+00:00',
  'files': ['2134', '3098', '1620320369', '705203'],
  'frame_number': 680,
  'mode': 'SM1',
  'observation_datetime': '2018-04-14T15:05:11.387000+00:00',
  'observation_direction': 'R',
  'off_nadir_angle': 38.7,
  'orbit_direction': 'A',
  'path_number': None,
  'polarisations': ['HH'],
  'publish_link': 'https://file.tellusxdp.com/api/v1/origin/publish/palsar2-l21/ALOS2210180680-180414'},
 {'bbox': [133.418, 34.35, 134.152, 35.045],
  'beam_number': 'U2-9',
  'begin_datetime': '2018-07-07T15:05:06.209000+00:00',
  'center_datetime': '2018-07-07T15:05:11.209000+00:00',
  'code': 'UBS',
  'coordinates': [[133.546, 34.35],
                  [134.152, 34.435],
                  [134.028, 35.045],
                  [133.418, 34.961]],
  'dataset_id': 'ALOS2222600680-180707',
  'date_added': '2020-11-11T16:35:58.475732',
  'end_datetime': '2018-07-07T15:05:16.209000+00:00',
  'files': ['2134', '3100', '1620320369', '705203'],
  'frame_number': 680,
  'mode': 'SM1',
  'observation_datetime': '2018-07-07T15:05:11.209000+00:00',
  'observation_direction': 'R',
  'off_nadir_angle': 38.7,
  'orbit_direction': 'A',
  'path_number': None,
  'polarisations': ['HH'],
  'publish_link': 'https://file.tellusxdp.com/api/v1/origin/publish/palsar2-l21/ALOS2222600680-180707'}]

比較する2枚の画像を選ぶ時の条件は差分干渉解析のときと同じです。ここでおさらいしておきましょう。

  • シーンフレーム:同一であること
  • 観測モード:同一であること
  • 偏波:同一であること
  • 観測角(オフナディア角):同一であること
  • 衛星進行方向(アセンディング/ディセンディング):同一であること
  • レーダ照射方向:同一であること

SAR観測は偏波によって感度が異なります(同じ対象でも受信信号の強度が変化する)。また、観測する地形と衛星の位置関係によって観測できない範囲が生じたり、歪みが生じます。観測条件を揃えない場合、これらの要因がノイズとして現れるようになるので、観測条件は同一である必要があります。

さて、設定した期間に2枚の観測シーンがあることが確認できました。運良く条件(観測条件が全て同一であること)にもあてはまっています。2枚目の方は、日本時刻で7月7日の24時ごろに観測されているので、災害直後の画像として適しています。1枚目の方は、2枚目の観測と同じ条件で観測されている、かつ、期間も3ヶ月前と比較的短いので、災害前のシーンとして採用することにします。

使用するシーン(オリジナルデータ)の入手

使用するシーンが決まったら、シーンId(dataset_id)を指定してオリジナルデータをダウンロードします。
※publish_file, download_fileを使用しています。

In [ ]:
import requests
from pprint import pprint

def publish_files(dataset_id):
    url = 'https://file.tellusxdp.com/api/v1/origin/publish/palsar2-l21/{}'.format(dataset_id)
    
    headers = {
        'Authorization': 'Bearer ' + TOKEN
    }

    r = requests.get(url, headers=headers)
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return r.json()

published1 = publish_files('ALOS2210180680-180414')
published2 = publish_files('ALOS2222600680-180707')
pprint(published1)
pprint(published2)
{'dataset_id': 'ALOS2210180680-180414',
 'expires_at': '2020-11-27T13:42:44.398627+00:00',
 'files': [{'file_name': 'IMG-HH-ALOS2210180680-180414-UBSR2.1GUA.tif',
            'file_size': 1620320369,
            'url': 'https://file.tellusxdp.com/api/v1/origin/eyJ0eXAiOiJKV1QiLCJhbGciOiJIUzI1NiJ9.eyJpYXQiOjE2MDY0NDEzNjQsImV4cCI6MTYwNjQ4NDU2NCwiYXBpX3Rva2VuIjoiNzM2YzFlMjNmYWViIiwiZGF0YXNldCI6InBhbHNhcjItbDIxIiwia2V5IjoiQUxPUzIyMTAxODA2ODAtMTgwNDE0IiwiZmlsZV9uYW1lIjoiSU1HLUhILUFMT1MyMjEwMTgwNjgwLTE4MDQxNC1VQlNSMi4xR1VBLnRpZiJ9.HT0MJYKrYyTXAqizXz29NQrlSA9MN4c_Arpsu4J19L4'},
           {'file_name': 'LUT-HH-ALOS2210180680-180414-UBSR2.1GUA.txt',
            'file_size': 705203,
            'url': 'https://file.tellusxdp.com/api/v1/origin/eyJ0eXAiOiJKV1QiLCJhbGciOiJIUzI1NiJ9.eyJpYXQiOjE2MDY0NDEzNjQsImV4cCI6MTYwNjQ4NDU2NCwiYXBpX3Rva2VuIjoiNzM2YzFlMjNmYWViIiwiZGF0YXNldCI6InBhbHNhcjItbDIxIiwia2V5IjoiQUxPUzIyMTAxODA2ODAtMTgwNDE0IiwiZmlsZV9uYW1lIjoiTFVULUhILUFMT1MyMjEwMTgwNjgwLTE4MDQxNC1VQlNSMi4xR1VBLnR4dCJ9.9i8DaN27eVrLeqMK3f2yW9_bpNJbIoI-kMwPRXLFfJc'},
           {'file_name': 'summary.txt',
            'file_size': 2134,
            'url': 'https://file.tellusxdp.com/api/v1/origin/eyJ0eXAiOiJKV1QiLCJhbGciOiJIUzI1NiJ9.eyJpYXQiOjE2MDY0NDEzNjQsImV4cCI6MTYwNjQ4NDU2NCwiYXBpX3Rva2VuIjoiNzM2YzFlMjNmYWViIiwiZGF0YXNldCI6InBhbHNhcjItbDIxIiwia2V5IjoiQUxPUzIyMTAxODA2ODAtMTgwNDE0IiwiZmlsZV9uYW1lIjoic3VtbWFyeS50eHQifQ.tYJPqEA-eBvUq_RXINb_G8rojQRf6l95s0X7wNi4e7k'},
           {'file_name': 'ALOS2210180680-180414_UBSR2.1GUA.kml',
            'file_size': 3098,
            'url': 'https://file.tellusxdp.com/api/v1/origin/eyJ0eXAiOiJKV1QiLCJhbGciOiJIUzI1NiJ9.eyJpYXQiOjE2MDY0NDEzNjQsImV4cCI6MTYwNjQ4NDU2NCwiYXBpX3Rva2VuIjoiNzM2YzFlMjNmYWViIiwiZGF0YXNldCI6InBhbHNhcjItbDIxIiwia2V5IjoiQUxPUzIyMTAxODA2ODAtMTgwNDE0IiwiZmlsZV9uYW1lIjoiQUxPUzIyMTAxODA2ODAtMTgwNDE0X1VCU1IyLjFHVUEua21sIn0.TZ0SLPRtRU4J5SSadF2MXxO3d2Byn_gYZDDQZ2ru570'}],
 'project': 'palsar2-l21'}
{'dataset_id': 'ALOS2222600680-180707',
 'expires_at': '2020-11-27T13:42:44.579937+00:00',
 'files': [{'file_name': 'IMG-HH-ALOS2222600680-180707-UBSR2.1GUA.tif',
            'file_size': 1620320369,
            'url': 'https://file.tellusxdp.com/api/v1/origin/eyJ0eXAiOiJKV1QiLCJhbGciOiJIUzI1NiJ9.eyJpYXQiOjE2MDY0NDEzNjQsImV4cCI6MTYwNjQ4NDU2NCwiYXBpX3Rva2VuIjoiNzM2YzFlMjNmYWViIiwiZGF0YXNldCI6InBhbHNhcjItbDIxIiwia2V5IjoiQUxPUzIyMjI2MDA2ODAtMTgwNzA3IiwiZmlsZV9uYW1lIjoiSU1HLUhILUFMT1MyMjIyNjAwNjgwLTE4MDcwNy1VQlNSMi4xR1VBLnRpZiJ9.hgLjPGe69qFwy3to-xAt-IgU6ypukZdZK3tNoq-_-fQ'},
           {'file_name': 'LUT-HH-ALOS2222600680-180707-UBSR2.1GUA.txt',
            'file_size': 705203,
            'url': 'https://file.tellusxdp.com/api/v1/origin/eyJ0eXAiOiJKV1QiLCJhbGciOiJIUzI1NiJ9.eyJpYXQiOjE2MDY0NDEzNjQsImV4cCI6MTYwNjQ4NDU2NCwiYXBpX3Rva2VuIjoiNzM2YzFlMjNmYWViIiwiZGF0YXNldCI6InBhbHNhcjItbDIxIiwia2V5IjoiQUxPUzIyMjI2MDA2ODAtMTgwNzA3IiwiZmlsZV9uYW1lIjoiTFVULUhILUFMT1MyMjIyNjAwNjgwLTE4MDcwNy1VQlNSMi4xR1VBLnR4dCJ9.idUZt6rfL_yafhvKKjIJGfVyODzl8GsKROBZXo7JeT4'},
           {'file_name': 'summary.txt',
            'file_size': 2134,
            'url': 'https://file.tellusxdp.com/api/v1/origin/eyJ0eXAiOiJKV1QiLCJhbGciOiJIUzI1NiJ9.eyJpYXQiOjE2MDY0NDEzNjQsImV4cCI6MTYwNjQ4NDU2NCwiYXBpX3Rva2VuIjoiNzM2YzFlMjNmYWViIiwiZGF0YXNldCI6InBhbHNhcjItbDIxIiwia2V5IjoiQUxPUzIyMjI2MDA2ODAtMTgwNzA3IiwiZmlsZV9uYW1lIjoic3VtbWFyeS50eHQifQ.uEhwdEm6g8cfi-kDvm8pBbJ_4eSfavGgIy8b7vyXrfs'},
           {'file_name': 'ALOS2222600680-180707_UBSR2.1GUA.kml',
            'file_size': 3100,
            'url': 'https://file.tellusxdp.com/api/v1/origin/eyJ0eXAiOiJKV1QiLCJhbGciOiJIUzI1NiJ9.eyJpYXQiOjE2MDY0NDEzNjQsImV4cCI6MTYwNjQ4NDU2NCwiYXBpX3Rva2VuIjoiNzM2YzFlMjNmYWViIiwiZGF0YXNldCI6InBhbHNhcjItbDIxIiwia2V5IjoiQUxPUzIyMjI2MDA2ODAtMTgwNzA3IiwiZmlsZV9uYW1lIjoiQUxPUzIyMjI2MDA2ODAtMTgwNzA3X1VCU1IyLjFHVUEua21sIn0.YKZK59Ro5f1FN0Wkn6_NiwBLQXfLFvRz_yTxkgNbHuA'}],
 'project': 'palsar2-l21'}
In [ ]:
def download_file(url):
    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)
        
file1 = download_file(published1['files'][1]['url'])
save_file(file1, published1['files'][1]['file_name'])

file2 = download_file(published2['files'][1]['url'])
save_file(file2, published2['files'][1]['file_name'])

print('done')
done

これで、オリジナルデータを手に入れることができました。このオリジナルデータを用いて解析をしていきます。
※オリジナルデータのTellus外への持ち出しは禁止されているので気をつけてください!

2時期のSAR強度画像を用いた解析

2時期の強度画像の差分の表現方法は大きく分けて2つあります。方法①は2時期の強度画像にそれぞれ別の色付けをし、重ね合わせることで変化箇所を可視化する方法です。方法②は2時期の強度画像の画素値の差を直接計算する方法です。

ここでは、これら2つの方法について両方学べるようになっていますが、方法②はやや中級者向けの内容になるので、まずは方法①の内容をきちんとおさえておきましょう。

方法① RGB疑似カラー合成画像による浸水範囲の可視化

2時期の強度画像を用いて変化した箇所をわかりやすく表現する方法として、RGB疑似カラー合成画像を作成し、変化箇所を可視化します。災害前の画像をR(赤)に、災害後の画像をG(緑)とB(青)に割り当てることで、災害前後で変化が大きい箇所は赤色、または、青色(シアン色)に、変化が小さい箇所は白黒として表現できます。

オリジナルデータ(GeoTIFF)を読み込み、処理時間節約のために関心領域で画像をクリップ、RGB疑似カラー合成画像を作成し、浸水範囲を推定します。ここではGDALを使いますが、処理は基本的に基礎分野と同じなので、詳細は基礎分野を参照してください。

必要なモジュール、ライブラリをインポートします。

In [ ]:
from osgeo import gdal, gdalconst
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

オリジナルデータを読み込み、座標系を確認します。座標系の説明も基礎分野を参照してください。

In [ ]:
img_before_original = gdal.Open('IMG-HH-ALOS2210180680-180414-UBSR2.1GUA.tif', gdalconst.GA_ReadOnly)
img_after_original = gdal.Open('IMG-HH-ALOS2222600680-180707-UBSR2.1GUA.tif', gdalconst.GA_ReadOnly)

src_prj=img_before_origina.GetProjection()
src_info=img_before_origina.GetGeoTransform()

print(src_prj)
print(src_info)
PROJCS["Geo-coded",GEOGCS["Datum=ITRF97 Ellipsoid=GRS80 Projection=UTM",DATUM["International_Terrestrial_Reference_Frame_1997",SPHEROID["GRS 1980",6378137,298.2572221010042,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6655"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4338"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",135],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
(356138.630928314, 2.5, 0.0, 3878517.7328678733, 0.0, -2.5)

座標系はUTMであることが確認できます。今回は関心領域の範囲を緯度経度で設定してあるので、解析する画像の座標系も緯度経度に変換して画像をクリップします。

In [ ]:
# 画像の座標系を変換、関心領域で画像をクリップ
gdal.Warp('ALOS2_180414_clip.tif', img_before_original, dstSRS='EPSG:4326', outputBounds=(bounds['minx'], bounds['miny'], bounds['maxx'], bounds['maxy']))
gdal.Warp('ALOS2_180707_clip.tif', img_after_original, dstSRS='EPSG:4326', outputBounds=(bounds['minx'], bounds['miny'], bounds['maxx'], bounds['maxy']))

print("done")
done

クリップした画像を表示して確認します。その前に、画素値のヒストグラムを確認し、表示するレンジを決定します。ヒストグラムに偏りが見られる場合は表示するレンジを適切に設定しておかないと、表示された画像が暗すぎる、もしくは逆に、明るすぎて見にくいことがあります。

In [ ]:
img_before = gdal.Open('ALOS2_180414_clip.tif', gdalconst.GA_ReadOnly)
img_after = gdal.Open('ALOS2_180707_clip.tif', gdalconst.GA_ReadOnly)

array_before = img_before.ReadAsArray()
array_after = img_after.ReadAsArray()

plt.figure(figsize=(20,5))

plt.subplot(1,2,1)

plt.hist(array_before.flatten(),bins=300,range=(1, 65545)) 
plt.title('before')

plt.subplot(1,2,2)
plt.hist(array_after.flatten(),bins=300,range=(1, 65545)) 
plt.title('after')

plt.show()

ヒストグラムを確認すると画素値はだいたい0-10000に集まっているので、表示するレンジを0-10000に設定します。

In [ ]:
plt.figure(figsize=(10,10))

plt.subplot(1,2,1)
plt.imshow(array_before,vmax=10000,cmap='gray')
plt.title('before')

plt.subplot(1,2,2)
plt.imshow(array_after,vmax=10000,cmap='gray')
plt.title('after')

plt.show
Out[ ]:
<function matplotlib.pyplot.show(*args, **kw)>

Original data provided by JAXA

災害前後の画像を表示してみました。SARでは反射(後方散乱)した電波の強度を記録しているため、光学画像(カラー)とは違って白黒の画像になります。明るいところが強い反射が得られている箇所、暗いところが弱い反射かあるいは反射が得られていない箇所を示しています。災害前後を比較してみると、災害後(浸水後)の画像では暗い範囲が増えているのがわかります。浸水時は水面で電波が鏡面反射してしまうため後方散乱が得られにくく、災害前後で暗く変化した箇所は水面になった(浸水した)可能性があります。

レーダ波が水面で鏡面反射する様子

災害前後の2枚の画像からRGB疑似カラー合成画像(R:災害前、G:災害後、B:災害後)を作成したとき、土砂崩れ等で散乱体となる樹木が消失した場合や浸水によって受信信号が得られにくい(強度が小さくなる)ときは赤色に、反対に、土砂崩れ等の崩壊土や浸水が引いた後に漂流物が堆積した場合はその堆積物が散乱体となるので強度が大きくなり、青色(正しくは、緑+青なのでシアン色)を呈します。災害前後で変化が見られない箇所は白黒として表示されます。

RGBカラーなので、画像のデータ型はByte型にしておいた方が扱いやすいです。変換する際は、表示した時のレンジ0-10000をスケールパラメータに指定しておくと、先ほど表示した画像の濃淡のままの画像ができます。

In [ ]:
# スケールを指定して画像をByte型に変換
gdal.Translate('ALOS2_180414_8bit.tif', img_before, outputType=gdal.GDT_Byte, scaleParams=[[0,10000]])
gdal.Translate('ALOS2_180707_8bit.tif', img_after, outputType=gdal.GDT_Byte, scaleParams=[[0,10000]])

print("done")
done

画像をByte型に変換したら、画像を結合してRGB疑似カラー合成画像を作成します。

In [ ]:
img_before_8bit = gdal.Open('ALOS2_180414_8bit.tif', gdalconst.GA_ReadOnly)
img_after_8bit = gdal.Open('ALOS2_180707_8bit.tif', gdalconst.GA_ReadOnly)

array_before_8bit = img_before_8bit.ReadAsArray()
array_after_8bit = img_after_8bit.ReadAsArray()

Xsize=img_before_8bit.RasterXSize
Ysize=img_before_8bit.RasterYSize
dtype=gdal.GDT_Byte
band=3

outfile = 'RGBcomposite.tif'

outRGB= gdal.GetDriverByName('GTiff').Create(outfile, Xsize, Ysize, band, dtype)
outRGB.SetProjection(img_before_8bit.GetProjection())
outRGB.SetGeoTransform(img_before_8bit.GetGeoTransform())

outRGB.GetRasterBand(1).WriteArray(array_before_8bit)
outRGB.GetRasterBand(2).WriteArray(array_after_8bit)
outRGB.GetRasterBand(3).WriteArray(array_after_8bit)
outRGB.FlushCache() 

plt.figure(figsize=(8,8))
image = mpimg.imread(outfile)
plt.imshow(image)
plt.title("RGB composite")

plt.show
Out[ ]:
<function matplotlib.pyplot.show(*args, **kw)>