【コード付き】TellusでPALSAR-2のL1.1の 画像化にチャレンジしてみた
本記事では、Tellusで公開しているPALSAR-2のL1.1データから、CEOSフォーマットと呼ばれるフォーマットに従って、複素画像（電波の位相と反射強度の情報を持った画像）を抽出し、可視化する方法を解説します。
Tellusでは、Tellus OS（ブラウザ上）で手軽に確認できるタイル化したデータに加えて、衛星データの提供元オリジナルデータを配信する機能をリリースしています。
詳しくはこちらの記事をご覧ください。
【ゼロからのTellusの使い方】提供元オリジナルデータを取得する
オリジナルデータは衛星データの知識がないとなかなか触るのが難しい代物ではありますが、使い方が分かるとグッと世界が広がります。
本記事では、Tellusで公開しているPALSAR-2のL1.1データから、CEOSフォーマットと呼ばれるフォーマットに従って、複素画像（電波の位相と反射強度の情報を持った画像）を抽出し、可視化する方法を解説します。
■PALSAR-2の処理レベルL1.1って？
L1.1は解像度の高い複素画像
SAR画像には、いくつかの処理レベルがあります。L1.1は、画像化した直後のデータであり、その後の処理のベースとなるプロダクトです。SLC（Single Look Complex）とも呼ばれ、解像度が最も高い状態の複素画像です。
電波の位相情報を有している
大きな特徴として、反射した電波の強度だけでなく、SARが発する電波の位相情報（-180°から+180°までの角度の情報）が保持されているため、別時期に撮影された画像と重ね合わせることで、地形の変動などを測定することができます。
PALSAR-2の処理レベルは他にも
その他のPALSAR-2の処理レベルについては以下の通りです。
-L1.0：画像化前のデータであり、画像化には高度な専門技術が必要になります
-L1.5：マルチルック処理により、画像の縦横の解像度を補正したもの。ただしこの時点で位相情報が失われます。
-L2.1：画像を既存の数値標高モデル（DEM）を用いてオルソ補正したもの。一般的な光学画像と同じような扱いができるのは、この処理レベルからということになります。
詳しくは以下をご覧ください。
【図解】衛星データの前処理とは~概要、レベル別の処理内容と解説~
■ PALSAR-2 L1.1のデータを取得しよう
Tellusで開発環境(Jupyter Lab)を使う方法は、「Tellusの開発環境でAPI引っ張ってみた」をご覧ください。
【ゼロからのTellusの使い方】Tellusの開発環境でAPI引っ張ってみた
プロダクトIDの取得
今回は、衛星が取得したプロダクトをそのままダウンロードするため、まずは、希望の場所や時期を指定して、プロダクトのIDを取得するプログラムを作成します。
まずは、以下のコードで、好きな緯度経度を指定して、プロダクトIDを取得してみてください。例は茨城県の霞ヶ浦です。
--サンプルコード⓪--
import os, requests
TOKEN = "自分のAPIトークンを入れてください"
def search_palsar2_l11(params={}, next_url=''):
if len(next_url) > 0:
url = next_url
else:
url = 'https://file.tellusxdp.com/api/v1/origin/search/palsar2-l11'
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()
def main():
ret = search_palsar2_l11({'after':'2019-1-01T15:00:00'})
lat = 36.046
lon = 140.346
for i in ret['items']:
lon_min = i['bbox'][0]
lat_min = i['bbox'][1]
lon_max = i['bbox'][2]
lat_max = i['bbox'][3]
if lat > lat_min and lat < lat_max and lon > lon_min and lon < lon_max:
print(i['dataset_id'])
if __name__=="__main__":
main()データセットのダウンロード
プロダクトIDを取得できたら、以下のコードで、プロダクトに含まれるデータセットを全てダウンロードします。
--サンプルコード①--
import os, requests
TOKEN = "自分のAPIトークンを入れてください"
def publish_palsar2_l11(dataset_id):
url = 'https://file.tellusxdp.com/api/v1/origin/publish/palsar2-l11/{}'.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()
def download_file(url, dir, name):
headers = {
"Authorization": "Bearer " + TOKEN
}
r = requests.get(url, headers=headers, stream=True)
if not r.status_code == requests.codes.ok:
r.raise_for_status()
if os.path.exists(dir) == False:
os.makedirs(dir)
with open(os.path.join(dir,name), "wb") as f:
for chunk in r.iter_content(chunk_size=1024):
f.write(chunk)
def main():
dataset_id = 'ALOS2249732860-190107'
print(dataset_id)
published = publish_palsar2_l11(dataset_id)
for file in published['files']:
file_name = file['file_name']
file_url = file['url']
print(file_name)
file = download_file(file_url, dataset_id, file_name)
print('done')
if __name__=="__main__":
main()データセットの中身
上記を実行すると、Jupyter上の現在のディレクトリに、シーンIDと同名のディレクトリが作成され、その中に、ファイルが格納されています。これでL1.1データのダウンロードは完了です。フォルダの中は以下のようになっています。
データセットの中身を簡単に説明すると、
・拡張子が.jpgの画像データ
衛星データを簡単に確認するためのサムネイル画像です。衛星データ本体はデータ容量が大きく中身を確認するのに時間がかかるため、解像度を落として画像化してあります。
・IMG-AA-XXXXXX.tif（ファイルの先頭がIMGとなっているファイル）
今回扱う対象となる画像データがこのファイルになります。
AAの部分に入っている文字列は画像データの種類を表しています。
データの種類は偏波（へんぱ）と呼ばれる送受信の電波の種類によって、全部で4種類（HH, HV, VH, VV）あります。送信（1文字目）と受信（2文字目）それぞれについて垂直偏波（V）か水平偏波（H）かを表しています。
上図のように全4種類あるプロダクトを一般にフルポラリメトリと呼びます。
今回は、HHの画像データを処理してみたいと思います。
他のファイルは、画像を処理する際に必要な情報が格納されたファイルです。処理する過程で必要に応じて参照していきます。
■ここが難関！CEOSフォーマットの解読
では、いよいよ画像データを処理していきましょう。
処理と言っても、他の光学画像のように二次元の画像配列が返ってくるだけじゃないの？と思われる方もいるかもしれません。
前回の記事（【ゼロからのTellusの使い方】提供元オリジナルデータを取得する）で紹介した通り、SARのデータはCEOSフォーマットに従って、1byteごとに情報が指定された位置に格納されていますので、必要な情報を自分で全て取りに行く必要があります。
PALSAR-2のプロダクトフォーマットはこちらから入手できます。
例えば以下のように書いてあります。
このプロダクトフォーマットを見て、生データを処理することを諦めた方も多いのではないでしょうか．．．そんな方のために今回は、SARの画像データを取り出すことに焦点を絞って、サンプルコードとともに解説していきます！
バイナリファイルを開く
まずはイメージデータが格納されたバイナリファイルを開きます。
dataset_id = 'ALOS2249732860-190107'
fp = open(os.path.join(dataset_id,'IMG-HH-ALOS2249732860-190107-HBQL1.1__D'),mode='rb')
バイナリファイルの操作の仕方
次に、バイナリファイルの操作の仕方を説明します。Pythonでは、seek()を使ってポインタを指定の位置まで送ります。たとえば、プロダクトフォーマットの以下の個所から、画像の縦横のサイズを取り出してみます。
fp.seek(236)
nline = int(fp.read(8))
print(nline)
fp.seek(248)
npixel = int(fp.read(8))
print(npixel)SARデータの引き出し
画像のサイズがわかったら、最後にSARデータを引き出します。フォーマットの該当箇所は以下のSAR RAW SGINAL DATAです。
正直にいって、これだけではよくわからないので、以下に簡単に図示してみます。
この図に従って、すべての複素画像データを取り込むスクリプトは以下の通りです。
nrec = 544 + npixel*8
fp.seek(720)
data = struct.unpack(">%s"%(int((nrec*nline)/4))+"f",fp.read(int(nrec*nline)))
# シグナルデータを32bit浮動小数点数として読み込み
data = np.array(data).reshape(-1,int(nrec/4)) # 2次元データへ変換
data = data[:,int(544/4):int(nrec/4)] # prefixを削除
slc = data[:,::2] + 1j*data[:,1::2]
nrec = 544 + npixel*8
fp.seek(720)
data = struct.unpack(">%s"%(int((nrec*nline)/4))+"f",fp.read(int(nrec*nline)))
# シグナルデータを32bit浮動小数点数として読み込み
data = np.array(data).reshape(-1,int(nrec/4)) # 2次元データへ変換
data = data[:,int(544/4):int(nrec/4)] # prefixを削除
slc = data[:,::2] + 1j*data[:,1::2]上記のコードでは、プロダクトの１シーン全てを処理してしまうので、開発環境によっては、メモリサイズが足りなくなり、動作が停止する可能性があります。その場合は、以下のスクリプトの通り、あらかじめラインとピクセルの範囲を限定して取り込むようにしてください。
import os
import numpy as np
import struct
def main():
dataset_id = 'ALOS2249732860-190107'
fp = open(os.path.join(dataset_id,'IMG-HH-ALOS2249732860-190107-HBQL1.1__D'),mode='rb')
col = 2000
row = 2000
cr = np.array([3000,3000+col,8000,8000+row],dtype="i4");
fp.seek(236)
nline = int(fp.read(8))
fp.seek(248)
ncell = int(fp.read(8))
nrec = 544 + ncell*8
nline = cr[3]-cr[2]
fp.seek(720)
fp.read(int((nrec/4)*(cr[2])*4))
data = struct.unpack(">%s"%(int((nrec*nline)/4))+"f",fp.read(int(nrec*nline)))
data = np.array(data).reshape(-1,int(nrec/4))
data = data[:,int(544/4):int(nrec/4)]
slc = data[:,::2] + 1j*data[:,1::2]
slc = slc[:,cr[0]:cr[1]]
print(slc)
if __name__=="__main__":
main()以上で、複素画像が完成しました。
■複素画像を位相と強度に分解する
複素画像が得られたら、各ピクセルに格納された複素数の位相成分と強度成分に分解します。サンプルコードは以下の通りです。
Numpyに便利な機能があるため、簡単に変換が可能です。なお、強度は校正式(https://www.eorc.jaxa.jp/ALOS-2/calval/calval_jindex.htm)に従って対数変換しています。また、位相はこの操作によって[-π～π]へ変換されます。
sigma = 20*np.log10(abs(slc)) -83.0 -32.0
phase = np.angle(slc)それでは、以下で画像化を行ってみましょう。
■可視化①強度画像
強度画像を可視化した結果は以下の通りです。ここでは、部分的に反射強度が強いピクセルや弱いピクセルがあったため、ヒストグラムの平坦補正も行っています。
import cv2
import matplotlib.pyplot as plt
sigma = np.array(255*(sigma-np.amin(sigma))/(np.amax(sigma)-np.amin(sigma)),dtype="uint8")
sigma = cv2.equalizeHist(sigma)
plt.figure()
plt.imshow(sigma, cmap = "gray")
plt.imsave('figure.jpg', sigma, cmap = "gray")
強度画像は、文字通り電波が反射した強度を示します。
金属など電波を反射するものは強く、樹木などは透過して、地面からの反射を観測しています。また、海や湖ではほとんど反射しないため、黒く映ります。
画像を確認すると、山や畑などの地形、幹線道路や大きな建物は確認できますが、光学画像ほど視認性は高くありません。実際には、光学画像などと組み合わせて情報を抽出することが多いのが、強度画像の特徴です。逆にSAR画像は光学画像と異なり、太陽の当たり方などに依存しない画像が得られるため、時系列の比較がしやすいことが利点です。
（くわしくは合成開口レーダ（SAR）のキホン～事例、分かること、センサ、衛星、波長～をご覧ください。）
SAR（合成開口レーダ）のキホン～事例、分かること、センサ、衛星、波長～
■可視化②位相画像
位相画像を可視化した結果は以下の通りです。
phase = np.array(255*(phase - np.amin(phase))/(np.amax(phase) - np.amin(phase)),dtype="uint8")
plt.figure()
plt.imshow(phase, cmap = "jet")
plt.imsave('phase.jpg', phase, cmap = "jet")
ランダムノイズのような画像に見えますね。そもそも、位相画像とはなにか、ここで簡単に解説します。
各ピクセルの位相は、上図のように、衛星がSAR信号を送信し対象物を反射して、受信した時の、最後の波が何度であったかを示しています。（この時、衛星も移動しているため、送信時と受信時とで位置が異なります）PALSAR-2の中心周波数は、約1.2GHzで、波長は約23cmです。地上の分解能が3m程度ですから、１ピクセルの中で位相情報は埋もれてしまっています。
また、仮に分解能がとても高かったとしても、観測したときの衛星の位置決定精度に誤差があったり、電離層によって位相がズレてしまうなど、残念ながら一度の観測による位相情報だけでは、反射した地形に関する情報を抽出することはできません。位相の情報を有効に活用するには、別の時期に取得された同じ地点のデータと厳密な位置合わせを行った上で比較するといった作業が必要になります。これを一般に、「干渉解析」と呼びます。
干渉解析は、一般に高度な処理技術が必要と言われていましたが、2月27日にリリースしたTellusマーケットでは、一般財団法人リモート・センシング技術センター（RESTEC）がRISEという干渉解析結果を提供しています。。宙畑でも干渉解析の実例を今後の記事で紹介していきますので、楽しみにしていてください。
※実際には、SAR信号は図のような搬送波を連続的に発しているわけではなく、周波数変調したパルス信号を一定間隔で送信していますので、図はあくまで位相とは何かを理解するための概念図として捉えてください。さらに余談ですが、そのパルス信号自体も衛星と地表の相対移動速度によってドップラーシフトしており、そのシフト量を利用して、地表面のどの位置から反射した電波であるかを判別します。これを全てのパルス信号に対して行うことで、あたかも静止した巨大なアンテナで受信したかのような受信信号を再現することができます。これが”合成開口”の仕組みなのですが、これ以上の詳細な説明はここでは割愛します。
SARの位相情報を用いた干渉解析では、地面の経時変化を抽出する「差分干渉SAR」が有名で、以前のこちらの記事でも詳しく解説しています。
https://sorabatake.jp/4343/
それ以外にも、位相情報を用いることで、数値標高モデル（digital elevation model: DEM）の作成を行うことができます。これは、同時に二機の衛星でSARの観測を行い、衛星間の距離に応じた視差効果を利用して、地表を立体視する技術です。これにより、地図を高頻度・高精度に更新することができたり、土砂崩れや雪崩などの地形の立体的な変化を捉えることによる災害規模の把握、工事現場の盛土の量の管理など、様々なアプリケーションが検討されています。
まとめ
本記事では衛星データ独特のCEOSフォーマットの取り扱いについて解説しました。
SARの提供元オリジナルデータを処理するのは、衛星データ特有の知識も必要ですが、基本を押さえれば誰にでもできることがお分かりいただけたと思います。
そして、衛星のオリジナルデータへの理解が深まれば、衛星データの価値をさらに高めることができるでしょう。
次回は、PALSAR-2のL2.1のフルポラリメトリック画像を用いて、土地被覆の可視化にチャレンジしてみたいと思います。
