ALOS/PALSAR L1.0データの読み込み (CEOSフォーマット)
概要
以前の記事 SAR衛星画像のRAWデータ再生処理を実装してみる (ALOS1) で、ALOS1のRAWデータの画像の再生処理の実装を紹介しました。今回の記事では、そこで詳しく触れられなかった、再生処理に必要なデータやパラメータを取得する処理を実装を紹介をします。
データの取得する最もシンプルな方法は、必要な値に対して、それが存在するファイルとバイトアドレスを調べて、ひとつずつ取得していくことです。しかし今回は、より汎用的な方法として、必要なメタデータ全体をまとめて抽出し、pandas.DataFrameの形式で扱うようにしました。
実装の際に、以下の公式文書を参考にしています。
『ALOS処理プロダクトフォーマット説明書 PALSARレベル1.0編』
(NEB-01006(ALOS-DPFT-J03)) ダウンロードページ
この文書にはL1.0プロダクトのデータセットの構成や、各ファイルのどこにどの情報が書き込まれているかが記述されています。本記事では「フォーマット説明書」として言及していきます。
データセットについて
ここでは前の記事と同じ scene_name: ALPSRS279162650 を例にします。 L1.0 データセットを取得すると、以下のようなファイル群が得られます。
- VOL-ALPSRP273720840-H1.0__A (ボリュームディレクトリファイル)
- LED-ALPSRP273720840-H1.0__A (SARリーダファイル)
- IMG-HH-ALPSRP273720840-H1.0__A (SARイメージファイル, "HH"は偏波を表す)
- TRL-ALPSRP273720840-H1.0__A (SARトレイラファイル)
(Dataset: ©JAXA/METI ALOS PALSAR L1.0 2011. Accessed through ASF DAAC 6 Jan 2025.)
今やりたいことは、上のファイル群から
- 生の信号データを取得する
- 再生処理に必要な情報を取得する
です。私の実装した再生処理では
- SARリーダファイル("LED-...")
- SARイメージファイル("IMG-...")
の2つのファイルから必要な情報を取得しました。
データセット内の各ファイルには「レコード」と呼ばれる単位で情報が格納されています。レコードは「姿勢データ」、「シグナルデータ」など、異なる種類の情報を含みます。
フォーマット説明書を見ると、それぞれのファイルにどのようなレコードが含まれているかが書かれています。上記の2つのファイルについては以下のようになります。
SARリーダファイル
| レコード長 | レコード数 | レコード名 |
|---|---|---|
| 720 | 1 | ファイルディスクリプタ |
| 4,096 | 1 | データセットサマリ |
| 4,680 | 1 | プラットフォーム位置データ |
| 8,192 | 1 | 姿勢データ |
| 13,212 | 1 | キャリブレーションデータ |
| 1,540,000 | 1 | 設備関連データ1 |
| 4,314,000 | 1 | 設備関連データ2 |
| 345,000 | 1 | 設備関連データ3 |
| 325,000 | 1 | 設備関連データ4 |
| 325,000 | 1 | 設備関連データ5 |
| 3,072 | 1 | 設備関連データ6 |
| 511,000 | 1 | 設備関連データ7 |
| 4,370,000 | 1 | 設備関連データ8 |
| 728,000 | 1 | 設備関連データ9 |
| 15,000 | 1 | 設備関連データ10 |
(フォーマット説明書より抜粋)
SARイメージファイル
| レコード長 | レコード数 | レコード名 |
|---|---|---|
| 720 | 1 | ファイルディスクリプタ |
| 可変表 | n | シグナルデータ |
(フォーマット説明書より抜粋)
各表の"レコード長"は、各レコードのサイズ(バイト長)、"レコード数" はそのレコードがファイル内にいくつ含まれているかを示します。
値の取得について
各レコードのバイトアドレスとデータの対応は、フォーマット説明書に記載されています。
例えば、シグナルデータレコードのPRF (パルス繰り返し周波数)の値については以下のように表に記載されています。
シグナルデータレコード
| フィールド No. | バイト No. | タイプ | 記述 (定義と値) | 備考 |
|---|---|---|---|---|
| 20 | 57 － 60 | B4 | PRF [mHz] | （PRF値） |
(フォーマット説明書より抜粋)
ここから、
「シグナルデータレコードの 57〜60 のバイトアドレスに、B4という型でPRF の値が格納されている」
ということがわかります。フィールド No. はこの値の項目の通し番号と思えば良いです。タイプはフォーマット説明書内で定義されていますが、"B4"は 「4バイトの符号なし整数型」と読み取れます。
部分的には、
field_data = record[57-1:60]
# "I": 4byte unsigned int
value = struct.unpack(">I", field_data)[0]
とすれば値は読み込めるわけです。
実装
実装は以下のような流れになっています。
- ファイルからレコードを読み出す
- 各レコードに対して、値の名称、バイトアドレス、データタイプの対応(データ構造)を定義
- データ構造をもとに、レコードから値を抽出し、pandas.DataFrame の形式で出力
1. ファイルからのレコードの抽出
以下はSARのリーダファイルからレコードを取り出すコードです。
def read_sar_leader_file(file_path):
"""
Reads the SAR leader file and extracts structured records based on the specifications.
"""
records = {}
# Define the structure of each record based on the SAR Leader specification
record_definitions = [
("file_descriptor", 720), # レコード 4: ファイルディスクリプタ
("dataset_summary", 4096), # レコード 5: データセットサマリ
("platform_position_data", 4680), # レコード 6: プラットフォーム位置データ
("attitude_data", 8192), # レコード 7: 姿勢データ
("calibration_data", 13212), # レコード 8: キャリブレーションデータ
# 以下は現状省略
# ("facility_related_data_1", 1540000), # レコード 9: 設備関連データ 1
# ("facility_related_data_2", 4314000), # レコード 10: 設備関連データ 2
# ("facility_related_data_3", 345000), # レコード 11: 設備関連データ 3
# ("facility_related_data_4", 325000), # レコード 12: 設備関連データ 4
# ("facility_related_data_5", 325000), # レコード 13: 設備関連データ 5
# ("facility_related_data_6", 3072), # レコード 14: 設備関連データ 6
# ("facility_related_data_7", 511000), # レコード 15: 設備関連データ 7
# ("facility_related_data_8", 4370000), # レコード 16: 設備関連データ 8
# ("facility_related_data_9", 728000), # レコード 17: 設備関連データ 9
# ("facility_related_data_10", 15000), # レコード 18: 設備関連データ 10
]
with open(file_path, "rb") as f:
current_position = 0
for record_name, length in record_definitions:
f.seek(current_position)
data = f.read(length)
records[record_name] = data
current_position += length
return records
基本的にはファイルに対してseek()で場所を指定し read()で読み取っていくだけです。
設備関連データは今回の再生処理では使わなかったので、抽出はしていません。この処理でファイル内の各レコードを辞書として読み出しています。
2. 各レコードに対して、値の名称、バイトアドレス、データタイプの対応(データ構造)を定義
ここではSARリーダーファイルから取得した、プラットフォーム位置データレコードを例にとって、データ構造を定義した部分を紹介します。
record_structure = [
# P3-37 (1/2)
["1-4", "B4", "record_number"], # レコード番号
["5-5", "B1", "first_record_subtype_code"], # 第1レコードサブタイプコード
["6-6", "B1", "record_type_code"], # レコードタイプコード
["7-7", "B1", "second_record_subtype_code"], # 第2レコードサブタイプコード
["8-8", "B1", "third_record_subtype_code"], # 第3レコードサブタイプコード
["9-12", "B4", "platform_data_record_length"], # プラットフォームデータレコード長
["13-44", "A32", "orbit_information"], # プラットフォーム情報（予測値や定値など）
["45-60", "F16.7", "scene_center_geocentric_vector_x"], # シーンセンタの地球中心ベクトル（X）
["61-76", "F16.7", "scene_center_geocentric_vector_y"], # シーンセンタの地球中心ベクトル（Y）
["77-92", "F16.7", "scene_center_geocentric_vector_z"], # シーンセンタの地球中心ベクトル（Z）
["93-108", "F16.7", "scene_center_geocentric_vector_x_prime"], # シーンセンタの地球中心ベクトル（X’）
["109-124", "F16.7", "scene_center_geocentric_vector_y_prime"], # シーンセンタの地球中心ベクトル（Y’）
["125-140", "F16.7", "scene_center_geocentric_vector_z_prime"], # シーンセンタの地球中心ベクトル（Z’）
["141-144", "I4", "data_point_count"], # データポイント数
["145-148", "I4", "first_point_year"], # 第1ポイントの年（YYYY）
["149-152", "I4", "first_point_month"], # 第1ポイントの月（MM）
["153-156", "I4", "first_point_day"], # 第1ポイントの日（DD）
["157-160", "I4", "first_point_day_of_year"], # 第1ポイントの年内日（DOY）
["161-182", "E22.15", "first_point_seconds"], # 第1ポイントの秒（小数部含む）
# P3-38 (2/2)
["183-204", "E22.15", "time_interval_between_points"], # ポイント間のインターバル時間（秒）
["205-268", "A64", "eci_ecr_coordinates"], # ECI/ECR座標情報
["269-290", "E22.15", "greenwich_mean_time"], # グリッニチ平均時角
["291-306", "F16.7", "scene_direction_along_track"], # シーン方向ベクトル（X）
["307-322", "F16.7", "scene_direction_cross_track"], # シーン方向ベクトル（Y）
["323-338", "F16.7", "scene_direction_radius_track"], # シーン方向ベクトル（Z）
["339-354", "F16.7", "scene_velocity_along_track"], # シーン速度ベクトル（X）
["355-370", "F16.7", "scene_velocity_cross_track"], # シーン速度ベクトル（Y）
["371-386", "F16.7", "scene_velocity_radius_track"], # シーン速度ベクトル（Z）
]
# 動的に追加データポイントを定義
start_byte = 387 # 第1データポイントの開始バイト
point_size = 132 # 1データポイントのバイト数
num_points = 28 # データポイント数
for i in range(1, num_points + 1): # 第1データポイントから第28データポイントまで
offset = start_byte + (i - 1) * point_size
record_structure.extend([
[f"{offset}-{offset + 21}", "E22.15", f"point_{i}_geocentric_vector_x"], # 地心ベクトルX
[f"{offset + 22}-{offset + 43}", "E22.15", f"point_{i}_geocentric_vector_y"], # 地心ベクトルY
[f"{offset + 44}-{offset + 65}", "E22.15", f"point_{i}_geocentric_vector_z"], # 地心ベクトルZ
[f"{offset + 66}-{offset + 87}", "E22.15", f"point_{i}_velocity_vector_x"], # 速度ベクトルX
[f"{offset + 88}-{offset + 109}", "E22.15", f"point_{i}_velocity_vector_y"], # 速度ベクトルY
[f"{offset + 110}-{offset + 131}", "E22.15", f"point_{i}_velocity_vector_z"], # 速度ベクトルZ
])
# フィールド No.35-37
record_structure.extend(
[
["4083-4100", "A18", "blank"], # 空白 (フィールド No.35)
["4101-4101", "I1", "leap_second_flag"], # うるう秒フラグ (フィールド No.36)
["4102-4680", "A579", "blank"], # 空白 (フィールド No.37)
]
)
["1-4", "B4", "record_number"]のように、[バイト No., タイプ, 値の名前 ]
を定義したリストが含まれています。 「第 N データポイントの速度ベクトル X 」のような繰り返しで定義されているところは、 for 文を使って動的に定義しています。
このような record_structure を各レコード毎にフォーマット説明書の表を元に作成していきます。普通にやるとかなりの労力が必要ですが、今回はChatGPTを使ってフォーマット説明書の PDF からrecord_structure の作成を手伝ってもらいました。ChatGPTに PDF を読み込ませる際に、日本語の一部が文字化けを含んでしまうため適宜前処理が必要だったり、一括処理では精度が悪いためページを分割して解析させる必要があったりしましたが、手打ちで作成していくよりは遥かに楽に作成することができました。
3. データ構造をもとに、レコードから値を抽出し、pandas.DataFrame の形式で出力
取得したレコードのデータと、レコードごとに定義した record_structure を入力にして、以下のよう関数を作りました。
def parse_record(data, record_structure):
"""
Parse binary data according to the record structure.
Args:
record_structure (list): List of record field definitions.
data (bytes): Binary data to parse.
Returns:
pd.DataFrame: Parsed data as a DataFrame.
"""
records = []
for field in record_structure:
byte_num, field_type, name = field
start_byte = int(byte_num.split("-")[0]) - 1
end_byte = int(byte_num.split("-")[-1])
field_data = data[start_byte:end_byte]
if field_type.startswith("B"): # Binary (unsigned integers)
byte_count = end_byte - start_byte
if byte_count == 1: # 1 バイトの場合
value = struct.unpack(">B", field_data)[0]
elif byte_count == 2: # 2 バイトの場合
value = struct.unpack(">H", field_data)[0]
elif byte_count == 4: # 4 バイトの場合
value = struct.unpack(">I", field_data)[0]
else: # 複数バイト（配列として解釈）
value = struct.unpack(f">{byte_count}B", field_data) # 配列として展開
elif field_type.startswith("I"): # Integer
try:
value = int(field_data.decode("ascii").strip())
except ValueError:
value = None
elif field_type.startswith("F"): # Floating-point
try:
value = float(field_data.decode("ascii").strip())
except ValueError:
value = None
elif (
field_type.startswith("A") # ASCII data
or field_type.startswith("CH") # Character data
):
value = field_data.decode("ascii").strip()
elif field_type.startswith("E"): # Exponential notation
try:
value = float(field_data.decode("ascii").strip())
except ValueError:
value = None
elif field_type == "3B4": # 3つの 4 バイト整数
if len(field_data) == 12: # 確認: データが12バイトあること
value = struct.unpack(">III", field_data) # 3つの 4 バイト整数を展開
else:
raise ValueError(
f"Invalid data length for 3B4: {len(field_data)} bytes"
)
else: # Unsupported or unknown type
value = None
# print(name, value)
records.append(
{
"byte_num": byte_num,
"type": field_type,
"name": name,
"value": value,
}
)
return pd.DataFrame(records)
このコードは、先述の「値の取得」のように、
field_data = record[57-1:60]
value = struct.unpack(">I", field_data)[0]
の処理を、フィールドを動的に定義したり、様々なデータタイプに対応できるよう一般化して書いたものになります。欲しいデータのレコードを pandas.DataFrame の形で取得することになります。
使用例
例として、SARリーダファイルのデータセットサマリレコードからレーダー波長を取得する場合は以下のようになります。
def get_value_from_metadf(metadf, name):
val = metadf.loc[metadf["name"] == name].iloc[0]["value"]
return val
dataset_summary_df = parse_record(
sar_reader_records["dataset_summary"], dataset_summary_structure
)
LAMBDA_RADAR = get_value_from_metadf(
ds_summary_df, "radar_wavelength"
)
以上のように、各パラメータを取得していく形になります。
SARイメージファイルの読み込みの詳細
ファイルからレコードを読み込み、値を取得する方法についてはこれまでの説明の通りです。
ここではSARイメージファイルのデータのシグナルデータセットの読み込みについての補足の説明をします。
先ほど見たように、SARイメージファイルは以下のような構成になっています。
| レコード長 | レコード数 | レコード名 |
|---|---|---|
| 720 | 1 | ファイルディスクリプタ |
| 可変表 | n | シグナルデータ |
ここで、シグナルデータについての詳細を以下の図に示します。
SAR衛星は基本的に一定時間ごとに送信したパルスから返ってくるパルスを受信するわけですが、このデータフォーマットではシグナルデータは受信したパルスごとに記録されています。1つのレコードがレンジ方向のシグナルに対応し、その長さはパルスの長さによって変わります。さらに、1つのシグナルレコードに対して、特定の長さの "Prefix Data"と呼ばれる領域が存在し、パルス取得時のパラメータが格納されています。"Prefix Data" の後に、"Sar RAW Signal Data" として実際のSARの信号の値が格納されています。 (実際のSARのシグナルデータには、ダミーデータ領域や右詰めのデータ数などがあり、どのデータが有効な領域か細かく定義されています。)
以上から、イメージファイルのデータには実際の信号のデータ以外のメタデータとして以下の二種類の情報があることがわかります。
- ファイルディスクリプタで定義される情報
- レンジ方向の各シグナルデータ(レンジライン)のごとに定義される情報
再生処理で用いるパラメータの中にも、 「PRF (パルス繰り返し周波数)」や「最初のデータ点までのスラントレンジ距離」などのように後者のレンジラインごとに定義されているものがあります。レンジラインごとに決まった値を使って再生処理するべきですが、今回の実装では、各レンジ行ごとに定義することはしていませず、最初の第1シグナルデータのPrefix Dataを全てのレンジの共通の値として使用しています。より精確な再処理を目指す際にはここは改善が必要になるかと思います。
信号成分を取得する際には、以下のようにSARイメージファイル全体を読み出した後、ファイルディスクリプタのレコード長 (N_HEADER)を除いた部分を(アジマスデータ数、レンジデータ数)の行列にreshapeして、その後にPrefix Dataの部分を除きます。シグナル部分は実部、虚部の順番で格納されているので、それぞれ順番に複素数にしてデータを取り出しています。
f = open(image_filepath, "rb")
rb = f.read()
r = np.frombuffer(rb, dtype=np.ubyte)
# ヘッダー領域を省く
# h,w がそれぞれ アジマス方向、レンジ方向のデータ数の行列にする
a_r_mtx = r[N_HEADER:].reshape((AZIMUTH_FULL, RECORD_LENGTH))
# レンジ方向のデータのみを抽出
a_r_mtx = a_r_mtx[:, RANGE_PREFIX : RANGE_PREFIX + RANGES_FULL]
# 実部と虚部を足し合わせて複素数にする
a_r_mtx = a_r_mtx[:, ::2] + a_r_mtx[:, 1::2] * 1j
まとめと課題
ALOS/PALSARのL1.0画像の再生処理に必要なデータや各パラメータを取得する処理について説明しました。現状、データセット全体の読み込みは未完了で、再生に必要なレコードのみを抽出しています。足りない部分として、
- データセット全ての情報を抽出する処理
- シグナルデータの各レンジラインごとのPrefix Dataを全て読み込む処理
などがあります。特に後者は、今回の再生処理の微調整にも影響があるのでご留意下さい。
Discussion