TellusのAPIからPALSAR2の生データを取得する
まず始めに
前回は、SAR(Synthetic Aperture Radar: 合成開口レーダ)に関する基礎的な知識だったり特徴などをお話しさせて頂きました。
今回は、Tellusという衛星データプラットフォームで公開されているPALSAR2というSARの画像の生データを抽出して可視化するために、
・TellusのAPIを用いてSAR画像のイメージファイルなどを取得
・イメージファイルのバイナリ情報を処理をして、位相情報が入った複素画像を取得
・実際に強度画像を表示してSARデータを取り扱ってみる。
というアプローチを行います。本当はGDALなどの便利なライブラリがpythonにはあるのですが、衛星ごとに画像のフォーマットが異なるため、既存ライブラリでは対応できない場合があります。私も最初はGDALを使ってみたのですが、公開されているPALSAR2のSAR画像を読み込むことができませんでした。なので、今回はJAXAが公開している説明書に書かれている情報を元にして必要なパラメータなどの情報を取り出し、各ピクセルごとの複素数データを抽出する必要性があります。
今回Tellusで取り扱うSAR画像について
今回取り出すのは、PALSAR2のSAR画像のプロダクトL1.1を使います(というよりも、現時点ではTellusではプロダクトL2.1が公開されていないので、GeoTIFF形式の強度画像を扱うことができません)。L1.1のプロダクトは、前回の記事で申し上げたようにCEOSという特殊なフォーマットを用いており、衛星ごとに仕様が異なります。なので、プロダクトのデータ内容を確認してバイナリデータの処理を行なってイメージデータを取り出す必要性があります。
ちなみに、最終的に可視化した結果はこちらになります。
使用した開発環境
今回は、以下の環境でJupyter Notebook上からSAR画像の生データを取り出して可視化します。TellusのサイトでJupyter環境を申請すれば誰でもAPIを通してあらゆる衛星データを扱えますので、もし興味がある方は、Tellusで開発環境を申請してみてください。ちなみに、自分は機械学習などの解析も考慮に入れているので、以下の高火力GPUコンピューティングの環境を申請して使用しております。
実際の開発環境 GPUカード：NVIDIA Tesla V100 for PCI-Express (32GB) ×4 CPU：Xeon E5-2623 v3 4Core×2 (8C/16T 3.0GHz Max3.5GHz) メモリ：128GB ディスク（RAID1）：SSD 480GB 2台/1組 Ubuntu: 18.04
高火力GPUコンピューティング環境については、実際にPytorchを導入してJupyterでリモートで操作できる環境を実装するまでの備忘録を過去の記事で公開しているので、興味がある方はこちらも参考にして環境のセットアップなどをしてみてください。
実際に取り出してみる
TellusのAPIを用いてSAR画像のイメージファイルなどを取得
SAR画像を取得するために、まずTellusのAPIを通してSAR画像のプロダクトのURL情報を抽出して、そこから取得したURLにアクセスして画像を取り出す形式になっています。TellusのAPIのリファレンスは以下で説明されているので、詳しくはこちらを参考にしてください。ソースコード上にも書いてますが、APIを利用するにはトークンが必要です。なので、Tellusにログインして、マイページにあるAPIアクセス設定で発行したAPIトークンを変数として入れる必要があります。また、APIのデフォルト設定ではHH偏波(水平偏波を送信して同じ偏波成分で受信されたもの)しか取得できないらしいです。なので、フルポラリメトリで取得したい場合には、TellusのAPIリファレンスには残念ながら載っていないのですが、"polarisations":"HH+HV+VV+VH"をパラメータとして入力する必要性があるので注意してください。
まずは、必要なライブラリとパラメータを書き込んでいきます。
#------------------------------------------------# # 必要なライブラリを読み込む import requests from PIL import Image import urllib.request import matplotlib.pyplot as plt %matplotlib inline import numpy as np import pandas as pd from io import BytesIO import codecs import struct #------------------------------------------------# #------------------------------------------------# ## APIのドメイン URL_DOMAIN = "file.tellusxdp.com" SEARCH_DOMAIN = "api/v1/origin/search/palsar2-l11" # PALSAR2の標準処理データ(L1.1)を検索する PUBLISH_DOMAIN = "api/v1/origin/publish/palsar2-l11" # データを取得するためのURLを手に入れる ## APIトークンをセットしてください BEARER_TOKEN = "Tellusにログインして、マイページにあるAPIアクセス設定で発行したAPIトークンを入れる" # 取得したい地域の緯度経度の範囲情報を入力。詳しい詳細はTellusのAPIリファレンスのページを参照してください # フルポラリメトリで衛星画像を取得したい場合は、"polarisations":"HH+HV+VV+VH"を追加します area = {"left_bottom_lon":130,"right_top_lon":140,"left_bottom_lat":35,"right_top_lat":40,"polarisations":"HH+HV+VV+VH"} headers = {'Authorization': 'Bearer %s' % BEARER_TOKEN} if BEARER_TOKEN == "": print("APIトークンがセットされていません") #------------------------------------------------#
次に、PALSAR2のデータを検索する関数、検索して出てきたプロダクトIDを元にしてSARデータを取得するためのURLを発行する関数、そのURLを元にSARデータを持ってくる関数という形で、合計で三つの関数を定義します。
# PALSAR2のデータを検索する関数(プロダクトIDを取得する) def get_data(domain, parameters, location): res = requests.get("https://{}/{}".format (domain, parameters), params=location, headers=headers) if res.status_code is not 200: raise ValueError('status error({}).'.format(res.status_code)) return res.json() # 取得したプロダクトIDにアクセスして、SARデータを取得するためのURLを手に入れるための関数 def get_URL(domain, parameters, dataset_id): res = requests.get("https://{}/{}/{}".format (domain, parameters,dataset_id), headers=headers) if res.status_code is not 200: raise ValueError('status error({}).'.format(res.status_code)) return res.json() # URLを元に実際のSAR画像(イメージファイル)やそのサマリーなどを取得する関数 def get_image(url): res = requests.get(url, headers=headers) if res.status_code is not 200: raise ValueError('status error({}).'.format(res.status_code)) return res
では、まずは取得したSAR画像のプロダクトIDを取得したいと思います。検索結果はデフォルトで100件出てくるのですが、今回はお試しですし多すぎるので10件だけ表示しています。パラメータとして入力した緯度経度の範囲内のプロダクトを検索してくれているみたいです。
res = get_data(URL_DOMAIN, SEARCH_DOMAIN, area) for i in range(10): print(res['items'][i]['dataset_id'])
実行結果: ALOS2127230690-160930 ALOS2127230700-160930 ALOS2127230710-160930 ALOS2127230720-160930 ALOS2127230730-160930 ALOS2127380690-161001 ALOS2127380700-161001 ALOS2127380710-161001 ALOS2127970690-161005 ALOS2127970700-161005
このようにして、実際にPALSAR2のプロダクトIDを取得することができました。
次に、対象のプロダクトIDを入れてSAR画像を取得するためのURLを手に入れます。今回は、最初に出てきたALOS2127230690-160930のプロダクトIDを入力しました。
ids = res['items'][0]['dataset_id'] SAR = get_URL(URL_DOMAIN, PUBLISH_DOMAIN,ids)
取得したURLをみると、以下のような形式になっています。
{'dataset_id': 'ALOS2127230690-160930',
'expires_at': '2020-01-09T01:20:09.452570+00:00',
'files': [{'file_name': 'BRS-HH-ALOS2127230690-160930-HBQR1.1__A.jpg',
'file_size': 82228,
'url': 'HH偏波の縮小画像のjpegファイルを取得できるURL情報'},
{'file_name': 'BRS-VV-ALOS2127230690-160930-HBQR1.1__A.jpg',
'file_size': 79224,
'url': 'VV偏波の縮小画像のjpegファイルを取得できるURL情報'},
{'file_name': 'IMG-VV-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1697149008,
'url': 'VV偏波のイメージファイルのURL情報'},
{'file_name': 'IMG-HH-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1697149008,
'url': ’HH偏波のイメージファイルのURL情報'},
{'file_name': 'VOL-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 2880,
'url': 'ボリュームディレクトリという、ファイルの管理情報などが格納されているファイルを取得できるURL情報'},
{'file_name': 'IMG-VH-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1697149008,
'url': 'VH偏波のイメージファイルのURL情報'},
{'file_name': 'TRL-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 568980,
'url': 'トレイラと呼ばれる、低分解能画像データが格納されているファイルを取得するためのURL情報'},
{'file_name': 'summary.txt',
'file_size': 1435,
'url': 'データのサマリ情報を取得できるURL'},
{'file_name': 'IMG-HV-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1697149008,
'url': 'HV偏波のイメージファイルのURL情報'},
{'file_name': 'LED-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1609432,
'url': 'SARリーダと呼ばれる、地図投影や姿勢、データ品質やラジオメトリックデータが格納されているファイルを取得するためのURL情報'},
{'file_name': 'BRS-HV-ALOS2127230690-160930-HBQR1.1__A.jpg',
'file_size': 75795,
'url': 'HV偏波の縮小画像のjpegファイルを取得できるURL情報'},
{'file_name': 'BRS-VH-ALOS2127230690-160930-HBQR1.1__A.jpg',
'file_size': 75672,
'url': 'VH偏波の縮小画像のjpegファイルを取得できるURL情報'}],
'project': 'palsar2-l11'}
今回は、SARの生データ(各ピクセルに複素数データが格納されている形式のデータ)を取り出して可視化するのが目的なので、'IMG-HH-ALOS2127230690-160930-HBQR1.1__A'(つまり、今回のプログラムでいうSAR['files'][3]['url'])のURLにアクセスして、HH偏波のイメージファイルを取得します。イメージファイルはかなりデータの容量が大きいので、処理が終了するまでに少しだけ時間がかかります。
thumbs = get_image(SAR['files'][3]['url']) data = thumbs
イメージファイルのバイナリ情報を処理をして、振幅と位相情報が入った複素画像を取得する
では、TellusのAPIを通してPALSAR2のHH偏波のイメージファイルを取得できたので、次はイメージファイルのバイナリ情報を処理してデータを抽出します。これはCEOSフォーマットになっているので、実際には以下のプロダクトの詳細に目を通して細かいパラメータ情報やその位置を把握する必要があります。
https://www.eorc.jaxa.jp/ALOS-2/doc/fdata/PALSAR-2_xx_Format_CEOS_J.pdf
ここでは簡潔にいうと、取得したSARイメージファイルの最初の720バイト分はファイルディスクリプタといい、ファイル情報の詳細が載っています。そのあとに処理済データレコードがあり、この中にSAR画像の生データの情報が入っていますが、以下の図ように、画像のラインごとにSARデータの前にプレフィックスデータが544バイト分入っているので、画像のラインごとにこのプレフィックスの情報を取り除いてからデータを取り出す必要があります。(PALSAR2のCEOSフォーマットの説明書で使用されていた図を抜粋しています)
また、SARデータの後に付随しているダミーデータや不定データのバイトはデフォルトで0となっているのはいいのですが、SARデータ自体はピクセル数分だけ情報がチャンクされている形になっています。さらに、1ピクセルあたり8バイトの情報のうち、前半4バイトには2の補数形式で実数情報(振幅情報)、後半4バイトには虚数成分(位相情報)が格納されています。どちらもfloat32の形式で取り出せますが、注意するところはビッグエンディアンで取り出す必要があるところです。自分が行なったjupyter上だと、デフォルトではリトルエンディアンで取り出されるようなので、今回はstructモジュールのunpackを使ってデータを処理しているのですが、>ffの形で8バイトのデータ情報を取り出すようにしてください。
それでは、まずはファイルディスクリプタにある基本的な情報を取り出して見ましょう
print(data.content[48:64]) print(data.content[16:27]) print(data.content[181:186]) # データセットサマリレコードの数 print(data.content[236:244]) # ライン数 print(data.content[248:256]) # ピクセル数 print(data.content[186:192]) # データセットサマリレコード長 SARデータレコード長 print(data.content[296:304])
出力結果: b'AL2 SARBIMOP ' b'CEOS-SAR ' b'23292' b' 23292' b' 9040' b' 72864' b' 13 4PB'
以下のようなバイト形式で取り出せました。ここで重要なのはライン数、ピクセル数、SARデータのレコード長です。実際にPALSAR2のSARデータを取り出す時には、これらの情報が非常に重要になります。
#ライン数、ピクセル数、SARデータレコード長を変数として定義する lines = int(data.content[236:244]) pixels = int(data.content[248:256]) Length = int(data.content[186:192])
次に、必要ないファイルディスクリプタの720バイト分の情報を取り除きます。これで、raw_dataに格納されているのはラインごとにプレフィックスとピクセル数分のSARイメージ情報が入ったデータのみとなります。
raw _data= data.content[720:]
試しに、最初のラインに入っているプレフィックスデータの情報を一部見て見ましょう。プロダクトの説明書通りであれば、左と右詰めのデータはないはずなので、ダミーデータなどはなく純粋にプレフィックスデータとSARデータのみで構成されていることが確認できると思います。
print(int(codecs.encode(raw_data[20:24], 'hex') ,16)) # 左詰め*(デフォルトは0) print(int(codecs.encode(raw_data[24:28], 'hex') ,16)) # ピクセル数 print(int(codecs.encode(raw_data[28:32], 'hex') ,16)) # 右詰め*(デフォルトは0) print(int(codecs.encode(raw_data[48:50], 'hex') ,16)) # フルポラリメトリのプロダクトなら4が出てくる
出力結果: 0 9040 0 4
では、実際にSARの生データをバイナリ情報を処理して抽出して見ます。
real = np.zeros((lines,pixels),np.float32) imag = np.zeros((lines,pixels),np.float32) import time t1 = time.time() for j in range(lines): for i in range(pixels): ret = struct.unpack('>ff', raw_data[Length*(j)+544 + 8*(i):Length*(j)+544 + 8*(i)+8]) real[j,i] =ret[0] imag[j,i] =ret[1] # 処理後の時刻 t2 = time.time() # 経過時間を表示 elapsed_time = t2-t1 print(f"経過時間：{elapsed_time}") # 自分の環境では、大体三分ちょいで処理が終わりました。
これで、ようやく各ピクセルに複素数の情報が格納されたSAR画像を取得することができました。
実際に強度画像を表示してSARデータを取り扱ってみる。
では、最後に複素数の情報を強度画像に変換して、グレースケールで可視化します。
from PIL import ImageOps result4 = np.zeros( (lines,pixels), dtype=np.complex64) result4.real = real result4.imag = imag resres = 10*np.log10(np.abs(result4)) print(resres.max(),resres.min()) img = Image.fromarray(resres).convert('L') plt.figure(figsize=(16,16)) plt.imshow(img, cmap='gray', interpolation='nearest')
可視化した結果が、前回の記事で例としてお見せしたSAR画像になっていて、以下の画像になります。
これで、ようやくPALSAR2の生データの取り出しと可視化に成功しました。これ以降は抽出したこれらの画像に対してスペックルのキャリブレーションを行なったり、他の偏波成分(HVやVV偏波)の複素数データを取り出して、HH,HV,VVの順番でRGB形式にしてその地域の散乱特性情報を可視化したりなど、様々なデータの加工や解析ができます。
まとめ
今回は、TellusのAPIを通じてPALSAR2の複素数データ(生データ)を取得して可視化するまでの一連の手順を紹介させていただきました。今後は、補正処理や実際のSAR解析の一端をご紹介できればと思っておりますので、よろしくお願い致します。
もしもpythonではなく、他の言語でSAR画像を研究などの都合で取り出したい場合は、昨日探していて見つけた記事ですが以下のものがありましたので、こちらもぜひ参考にして見てください(最近の記事として投稿されていたので驚きました笑)
