📘

JAXA CEOS フォーマットの構造と読み込み例

に公開

はじめに

https://zenn.dev/syu_tan/articles/6136d824686784

本記事は、書籍出版の決定を記念して執筆しました。書籍で得られる知識や実装を応用することで実現できる例として、JAXA CEOS フォーマットの構造と読み込み例 を用意しました。書籍にも記載したかったのですが、ページ数的に断念したのでこちらの記事で無料公開しております。完全に趣味なので誰かの知見になると嬉しいです。

まえおき

合成開口レーダーとは


© iQPS Inc.
https://chizaizukan.com/property/754/ より引用

合成開口レーダー(Synthetic Aperture Radar; SAR) は、全天候・昼夜を問わず高解像度のマイクロ波画像を取得できるという特性から、近年ますます注目を集めているリモートセンシング技術です。この技術は、デジタル標高モデルの生成、火山活動や洪水災害の観測、陸上および海上交通の監視、植生成長のモニタリング、海流や流氷の観測、さらには海洋の油流出の検出など、幅広い分野で活用されています。小型SAR衛星も勢い付いています!

衛星データのフォーマット

衛星データのフォーマットは以下のようにいくつか種類があります。

フォーマット名 概要 ファイル拡張子 対応衛星/センサ例
GeoTIFF 画像ファイル形式であるTIFFに地理情報を付加した形式。位置情報を持つ画像ファイルで、GISソフトウェアでの利用が可能。一般的な画像ビューアでも閲覧可能。タグの T が頭文字です。タグ情報に位置情報が含まれているファイル構成 .tiff, .tif ALOS/AVNIR-2、PRISM、PALSAR、ALOS-2/PALSAR-2、CIRC、SLATS/SHIROP
COG Cloud Optimized GeoTIFFの略。クラウド環境での画像データ取り扱いを効率化するため、ファイルの必要な部分のみにアクセスする機能を持つ。GeoTIFFと同様の環境で利用可能。 .tiff, .tif -
NetCDF Network Common Data Formの略。配列指向の科学データを統合的に扱う自己記述型の形式。フリーのライブラリを用いて取り扱いが可能で、専用ツールでデータ閲覧も可能。 .nc GCOM-C/SGLI
HDF Hierarchical Data Formatの略。多様なコンピュータ環境で階層的なデータを取り扱う自己記述型の形式。HDF4とHDF5の二つのバージョンがあり、互換性はない。フリーのライブラリを用いて取り扱いが可能で、専用ツールでデータ閲覧も可能。 .hdf, .h4, .hdf4, .he2, .h5, .hdf5, .he5 MOS-1、MOS-1b/MSR、VTIR、ADEOS/OCTS、ADEOS-II/AMSR、GLI、Aqua/AMSR-E、TRMM/PR、GCOM-C/SGLI、GCOM-W/AMSR2、GPM/DPR、GSMaP
CEOS準拠フォーマット CEOSで標準化されたファイル形式に準拠して作成された形式。バンド毎にファイルが分割される「CEOS-BSQ」と、複数のバンドが多重化された「CEOS-BIL」がある。ファイル構成は衛星データ毎に設定されており、詳細はプロダクトフォーマット説明書等に記載。 (データ作成機関により設定) MOS-1、MOS-1b/MESSR、MSR、VTIR、JERS-1/OPS、SAR、ADEOS/AVNIR、ALOS/AVNIR-2、PRISM、PALSAR、ALOS-2/PALSAR-2、ALOS-4/PALSAR-3

これらのフォーマットに関する詳細情報は、JAXAの公式サイトで紹介されています。特に今回は、この中のCEOSフォーマットについて整理して Python から読み込みの実装をしていきましょう。

CEOSフォーマットについて

CEOSフォーマットを学ぶ前に、CEOSという機関を知っておくといいと思います。

CEOS とは

地球観測衛星から得られるデータの標準化と相互運用性を確保するため、地球観測衛星委員会(CEOS: Committee on Earth Observation Satellites) が1984年に設立されました。この組織は、各国の宇宙機関が協力し、データフォーマットやカタログの相互運用性など、地球観測に関する標準化や国際的な情報インフラの整備、将来的な開発や協力に関する議論を行っています。フォーマットが統一されていると自動的に読み込みしやすいですよね。開発のための共有コストも減ります。

https://ceos.org/

CEOSフォーマットとは

CEOSは、衛星データのフォーマット標準化にも取り組んでおり、SAR(合成開口レーダー)衛星データのフォーマットもその一環として定められています。例えば、日本の宇宙航空研究開発機構(JAXA)が運用するALOS-2衛星のPALSAR-2データは、CEOSが定めたフォーマットに準拠しています。

CEOSの活動やデータフォーマットの標準化に関する取り組みは、各国の宇宙機関や国際機関との協力によって進められています。以下のようなプロバイダーによる CEOSフォーマットが存在します。

項目 JAXA CEOS ERSDAC CEOS
提供元 宇宙航空研究開発機構(JAXA) 資源環境観測解析センター(Earth Remote Sensing Data Analysis Center; ERSDAC)
前身 NASDA(宇宙開発事業団) METI関連機関
対象データ ALOS/PALSAR, AVNIR-2, PRISM など ALOS/PALSAR など(JAXAのデータ提供前に運用)
フォーマットの特徴 CEOS標準フォーマットに準拠 CEOS標準フォーマットに準拠しているが、JAXA CEOSとは処理アルゴリズム、メタデータやデータ構造に若干の違いがある
データ配布先 JAXA公式サイト, ASF, ESA など 過去のERSDACデータは現在JAXAに統合
参考資料 JAXA Earth Observation Data ESA ALOS PALSAR CEOS Format

だから ALOS PALSAR の同じ観測データでも JAXA CEOS 仕様書ERSDAC CEOS 仕様書などのややこしいことが発生します。。。。
JAXAと経産省で作ったからそれぞれのプロダクトが生まれてる政治的背景ですね

とはいえ、フォーマットが大体決まることは大変喜ばしく、日本の小型SAR衛星のベンチャーである Synspectiveも CEOS に合わせたプロダクトも提供してくれたりしています。

Synspective SARデータプロダクトガイド
https://sol.synspective.com/SAR_Data_Product_Guide_JP_v11.0.pdf より引用

CEOSファイル構成

では、CEOS のファイル構成とレコードのについてご紹介します。
より具体的にするために、ALOS/PALSARのCEOSフォーマットを例にファイル構成と主要なレコードについて整理いたします。以下に、各ファイルの役割とその主なレコードは以下のようになっています。ファイル名では、ALOS/PALSARではプレフィックスに記載、ERS-1,2, JERS-1では拡張子に記載されています。

ファイル概要

ファイル 記載ファイル名 内容・役割
ボリュームディスクリプタファイル VOL-, N/A データセット全体の概要情報を提供する。
リーダファイル LED-, .L データセットの詳細情報や処理パラメータを含む。
イメージファイル IMG-, .D 実際のSAR画像データを格納する。
トレイラファイル TRL-,.P 追加情報や品質評価データを含む。

基本的には、ファイルは4つで構成されています。

ですが、偏波が増えることで イメージファイル が観測偏波の種類だけ増えます。

例えば、HHHV の2偏波のデータの場合はこのようになります。

レコードの概要

レコード名 説明
ボリュームディスクリプタレコード データセット全体の識別情報や作成機関、作成日時などを記載。
テキストレコード データセットの概要説明やプロダクトID、作成場所、作成日時などを含む。
ファイルディスクリプタレコード 各ファイルの構造や内容に関する詳細情報を提供。
データセットサマリレコード データセット全体のサマリ情報、観測モード、処理レベル、地理的範囲などを記載。
プラットフォーム位置データレコード 衛星の位置や速度ベクトルなどの軌道情報を含む。
設備関連データレコード センサや衛星の状態、設定パラメータなどの設備に関する情報を提供。
イメージラインレコード(シグナルデータレコード) SARデータの各ラインの信号データを格納。
低分解能画像データレコード 低解像度のプレビュー画像や品質評価用のデータを含む。

バイト位置の概要

ファイルについてはバイナリーデータなので特別なシステムを使用しないと読み込めません。
そのバイナリーのデータをファイルの **番地(バイト位置)**によって読み取ることが可能です。

例えば、リーダファイルディスクリプタレコード(1/12) が以下のようになっています。


https://www.eorc.jaxa.jp/ALOS/jp/alos/fdata/PALSAR_x_Format_JL.pdf 3-16より引用

特に、知りたい情報がレコード長ならば、9 ~ 124 バイトのデータを読み込むことで取得が可能です。

フィールド番号 バイト位置 データ型 内容 備考
6 9-12 B4 レコード長
7 13-14 A2 ASCII/EBCDICコード 'Ab':ASCII

情報の数字の大きさに応じてデータ型が変化しているために、長さ(バイト) も変化します。ASCII/EBCDICコードの場合は文字列なので 17 ~ 28 になっています。

読み込み順序

読み込みをするバイト位置を知っていれば、任意の情報を取得できることになります。
しかしながら、データによってレコードの長さなどが違う関係で観測方法によっては、同じ場所を読み込んでも失敗してしまいます。レコード長が可変ということです。CEOS の構成ではレコードを読み込んだ時に、自身のレコード長を先に取得できます。そのレコード長を取得した上で、次のレコードを読み込んでいきます。

Python 読み込み

CEOSのファイル構成も理解したところで、Python実装をしてデータを読み込んでいきましょう。
使用するライブラリは Pythonの標準ライブラリnumpy などの依存の少ないベース環境を使用します。実装していくのは ALOS/PALSAR L1.0の仕様書からデータを読み取っていきます。このような信号のみのプロダクトは他の方法で配布していない(GeoTiffなどはない)ので CEOSを読み込んでいくしかありません。サンプルデータはシーンID: ALPSRP051270700 です。

import os, struct, math
import numpy as np

まずは、イメージファイル(IMG) を読み込んでいきます。バイナリーのファイルの展開方法です。

f = open(PATH_IMG, "rb")

余談にはなりますが、読み込み終えると以下の処理で必ず閉じましょう。

f.close()

次に、レコード長を取得してみます。

f.seek(8)
NUM_SAR_DISCRIPTOR_RECORD = int.from_bytes(f.read(INTERGER4), byteorder="big")
print('6 9 - 12 B4 レコード長 = 720)10 ->', NUM_SAR_DISCRIPTOR_RECORD)
6 912 B4 レコード長 = 72010 -> 720

他にもまとめて一定量のバイト位置からバイト位置までまとめて読み込む場合は、numpy の frombuffer などで読み込むことが可能です。以下は信号データをレンジ方向にまとめて読み込む時の例です。

byte_ri = f.read(NUM_PIXEL * 2)
ri = np.frombuffer(byte_ri, dtype=np.uint8)
ln = int(len(ri)/2)
signal[i, :ln] = ri[0::2] + ri[1::2] * 1j

そして、衛星の軌道情報などが記載されている リーダーファイル(LED) を読み込んでいきましょう。

f = open(PATH_LED, "rb")

リーダファイルでは、取得できるレコードも多いので以前のレコード長を足しあわせて読み込むことになります。

f.seek(8)
record_length = int.from_bytes(f.read(INTERGER4), byteorder="big")
f.seek(180)
summary_record = int(f.read(INTERGER6))
print('25 181 - 186 I6 データセットサマリレコードの数 = bbbbb1 ->', summary_record)
f.seek(186)
summary_record_length = int(f.read(INTERGER6))
print('26 187 - 192 I6 データセットサマリレコード長 = b4096 ->', summary_record_length)
f.seek(192)
map_record = int(f.read(INTERGER6))
print('27 193 - 198 I6 地図投影データレコードの数 = bbbbb0 ->', map_record)
f.seek(198)
map_record_length = int(f.read(INTERGER6))
print('28 199 - 204 I6 地図投影データレコード長 = bbbbb0 ->', map_record_length)
f.seek(210)
platform_record_length = int(f.read(INTERGER6))
print('30 211 - 216 I6 プラットフォーム位置データレコード長 = b4680 ->', platform_record_length)
f.seek(222)
attitude_record_length = int(f.read(INTERGER6))
print('32 223 - 228 I6 姿勢データレコード長 = b8192 ->', attitude_record_length)
f.seek(234)
radiometric_record_length = int(f.read(INTERGER6))
print('34 235 - 240 I6 ラジオメトリックデータレコード長 = bbbbb0 ->', radiometric_record_length)
f.seek(246)
radiometric_comp_record_length = int(f.read(INTERGER6))
print('36 247 - 252 I6 ラジオメトリック補償レコード長 = bbbbb0 ->', radiometric_comp_record_length)
f.seek(258)
data_quality_record_length = int(f.read(INTERGER6))
print('38 259 - 264 I6 データ品質サマリレコード長 = bbbbb0 ->', data_quality_record_length)
f.seek(270)
data_histogram_record_length = int(f.read(INTERGER6))
print('40 271 - 276 I6 データヒストグラムレコード長 = bbbbb0 ->', data_histogram_record_length)
f.seek(282)
range_spectrum_record_length = int(f.read(INTERGER6))
print('42 283 - 288 I6 レンジスペクトルレコード長 = bbbbb0 ->', range_spectrum_record_length)
f.seek(294)
dem_record_length = int(f.read(INTERGER6))
print('44 295 - 300 I6 DEMディスクリプタレコード長 = bbbbb0 ->', dem_record_length)
f.seek(342)
calibration_record_length = int(f.read(INTERGER6))
print('52 343 - 348 I6 キャリブレーションレコード長 = b13212 ->', calibration_record_length)
25 181186 I6 データセットサマリレコードの数 = bbbbb1 -> 1
26 187192 I6 データセットサマリレコード長 = b4096 -> 4096
27 193198 I6 地図投影データレコードの数 = bbbbb0 -> 0
28 199204 I6 地図投影データレコード長 = bbbbb0 -> 0
30 211216 I6 プラットフォーム位置データレコード長 = b4680 -> 4680
32 223228 I6 姿勢データレコード長 = b8192 -> 8192
34 235240 I6 ラジオメトリックデータレコード長 = bbbbb0 -> 0
36 247252 I6 ラジオメトリック補償レコード長 = bbbbb0 -> 0
38 259264 I6 データ品質サマリレコード長 = bbbbb0 -> 0
40 271276 I6 データヒストグラムレコード長 = bbbbb0 -> 0
42 283288 I6 レンジスペクトルレコード長 = bbbbb0 -> 0
44 295300 I6 DEMディスクリプタレコード長 = bbbbb0 -> 0
52 343348 I6 キャリブレーションレコード長 = b13212 -> 13212

これによって全体を構成しているレコードごとのレコード長を取得します。

では、表3.3-5 データセットサマリレコードの場合は、以下のようになります。

f.seek(record_length + 8)
summary_record_length = int.from_bytes(f.read(INTERGER4), byteorder="big")
print('6 9 - 12 B4 データセットサマリレコード長 = 4096)10 ->', summary_record_length)
f.seek(record_length + 20)
scene_id = f.read(FLOAT32).decode('utf-8')
print('9 21 - 52 A32 シーンID ->', scene_id)
f.seek(record_length + 68)
scene_time = f.read(FLOAT32).decode('utf-8')
print('11 69 - 100 A32 シーンセンター時刻 ->', scene_time)
f.seek(record_length + 164)
ellipsoid_model = f.read(FLOAT16).decode('utf-8')
print('16 165 - 180 A16 楕円体モデル ->', ellipsoid_model)
f.seek(record_length + 180)
ellipsoid_radius = float( f.read(FLOAT16))
print('17 181 - 196 F16.7 楕円体の半長径(Km) ->', ellipsoid_radius)
f.seek(record_length + 196)
ellipsoid_short_radius = float( f.read(FLOAT16))
print('18 197 - 212 F16.7 楕円体の短半径(Km) ->', ellipsoid_short_radius)
f.seek(record_length + 212)
earth_mass = float( f.read(FLOAT16))
print('19 213 - 228 F16.7 地球の質量 (10^24 Kg) ->', earth_mass)
f.seek(record_length + 244)
j2 = float(f.read(FLOAT16))
j3 = float( f.read(FLOAT16))
j4 = float( f.read(FLOAT16))
6 912 B4 データセットサマリレコード長 = 409610 -> 4096
9 2152 A32 シーンID -> ALPSRP051270700                 
11 69100 A32 シーンセンター時刻 -> 20070110130938062               
16 165180 A16 楕円体モデル -> GRS80           
17 181196 F16.7 楕円体の半長径(Km) -> 6378.137
18 197212 F16.7 楕円体の短半径(Km) -> 6356.7523141
19 213228 F16.7 地球の質量 (10^24 Kg) -> 5.974

次に、3.3-6 プラットフォーム位置データ・レコードでは以下です。

f.seek(record_length + summary_record_length + 8)
platform_record_length = int.from_bytes(f.read(INTERGER4), byteorder="big")
print('6 9 - 12 B4 プラットフォーム位置データレコード長 = 4680)10 ->', platform_record_length)
# ALOS軌道情報(予測値) : '0bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb'
# ALOS軌道情報(決定値) : '1bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb'
# ALOS高精度軌道情報 : '2bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb'
f.seek(record_length + summary_record_length + 12)
orbit_type = f.read(FLOAT32).decode('utf-8')
print('7 13 - 44 A32 軌道要素種類 ->', orbit_type)
# ALOS軌道情報(予測値) : 'bb28'
# ALOS軌道情報(決定値) : 'bb28'
# ALOS高精度軌道情報 : 'bb28'
f.seek(record_length + summary_record_length + 140)
NUM_ORB_POINT = int(f.read(INTERGER4))
print('14 141 - 144 I4 データポイント数 ->', NUM_ORB_POINT)
f.seek(record_length + summary_record_length + 144)
TIME_ORB_YEAR = int(f.read(INTERGER4))
TIME_ORB_MONTH = int(f.read(INTERGER4))
TIME_ORB_DAY = int(f.read(INTERGER4))
TIME_ORB_COUNT_DAY = int(f.read(INTERGER4))
TIME_ORB_SEC = float(f.read(FLOAT22))
print('15 145 - 148 I4 第1ポイントの年 ->', TIME_ORB_YEAR)
print('16 149 - 152 I4 第1ポイントの月 ->', TIME_ORB_MONTH)
print('17 153 - 156 I4 第1ポイントの日 ->', TIME_ORB_DAY)
print('18 157 - 160 I4 第1ポイントの通算日 ->', TIME_ORB_COUNT_DAY)
print('19 161 - 182 E22.15 第1ポイントの通算秒 ->', TIME_ORB_SEC)
f.seek(record_length + summary_record_length + 182)
TIME_INTERVAL = float(f.read(FLOAT22))
6 912 B4 プラットフォーム位置データレコード長 = 468010 -> 4680
7 1344 A32 軌道要素種類 -> 2                               
14 141144 I4 データポイント数 -> 28
15 145148 I4 第1ポイントの年 -> 2007
16 149152 I4 第1ポイントの月 -> 1
17 153156 I4 第1ポイントの日 -> 10
18 157160 I4 第1ポイントの通算日 -> 10
19 161182 E22.15 第1ポイントの通算秒 -> 46620.0
20 183204 E22.15 ポイント間のインターバル時間(秒) -> 60.0

次に、3.3-8 キャリブレーションデータレコード では以下です。

f.seek(record_length + summary_record_length + platform_record_length + attitude_record_length + 8)
calibration_record_length = int.from_bytes(f.read(INTERGER4), byteorder="big")
print('6 9 - 12 B4 レコード長 = 13212)10 ->', calibration_record_length)
f.seek(record_length + summary_record_length + platform_record_length + attitude_record_length + 16)
valid_sample = int(f.read(INTERGER4))
calibration_start = f.read(17).decode('utf-8')
calibration_end = f.read(17).decode('utf-8')
attenuator = int(f.read(INTERGER4))
alc = int(f.read(BYTE1))
agc = int(f.read(BYTE1))
pulse_width = int(f.read(INTERGER4))
chirp_bandwidth = int(f.read(INTERGER4))
sampling_frequency = int(f.read(INTERGER4))
quantization_bit = int(f.read(INTERGER4))
chirp_replica = int(f.read(INTERGER4))
chirp_replica_line = int(f.read(INTERGER4))
print('8 17 - 20 I4 有効サンプル数=Nsamp ->', valid_sample)
print('9 21 - 37 A17 校正データ取得開始時刻 ->', calibration_start)
print('10 38 - 54 A17 校正データ取得終了時刻 ->', calibration_end)
print('11 55 - 58 I4 校正器ATT設定値 ->', attenuator)
print('12 59 - 59 I1 校正器ALC ->', alc)
print('13 60 - 60 I1 AGC/MGC ->', agc)
print('14 61 - 64 I4 送信パルス幅 ->', pulse_width)
print('15 65 - 68 I4 チャープ帯域 ->', chirp_bandwidth)
print('16 69 - 72 I4 サンプリング周波数 ->', sampling_frequency)
print('17 73 - 76 I4 量子化ビット数 ->', quantization_bit)
print('18 77 - 80 I4 チャープレプリカデータ数 ->', chirp_replica)
print('19 81 - 84 I4 チャープレプリカデータ積算ライン数n ->', chirp_replica_line)
6 912 B4 レコード長 = 1321210 -> 13212
8 1720 I4 有効サンプル数=Nsamp -> 864
9 2137 A17 校正データ取得開始時刻 -> 20060710230952978
10 3854 A17 校正データ取得終了時刻 -> 20060710230953177
11 5558 I4 校正器ATT設定値 -> 15
12 5959 I1 校正器ALC -> 1
13 6060 I1 AGC/MGC -> 1
14 6164 I4 送信パルス幅 -> 27
15 6568 I4 チャープ帯域 -> 28
16 6972 I4 サンプリング周波数 -> 32
17 7376 I4 量子化ビット数 -> 5
18 7780 I4 チャープレプリカデータ数 -> 1
19 8184 I4 チャープレプリカデータ積算ライン数n -> 399
20 8585 I1 受信偏波1 -> 0
21 86 - α Nsamp*(2B2) チャープレプリカデータ1 -> (864, 2)

このようにして好きな観測値を取得できます。
これらを組み合わせることで、以下のようにチャープ信号の帯域による変化率を求めたりできます。

FREQ_RANGE_CHIRP_RATE = -(chirp_bandwidth * 1e6) / TIME_PLUSE_DURATION

特に衛星データは、軌道情報から結像や幾何補正に必要なことが多いです。CEOSには大切な情報が豊富なので一見すると触りにくいですがしっかりと考えられているフォーマットであることがわかります。普段は知られない詳細なデータも知ることができます。

CEOSの読み込みの応用例としては、以下の記事のように ALOS PLASAR の結像処理をするための観測情報は CEOS から読み取って算出でできます。


ALOS PLASAR © JAXA, METI

近日公開予定です。リンクはこちら

https://zenn.dev/syu_tan/articles/83c1c0292a1530

https://github.com/syu-tan/sar-python-book

さいごに

最後まで読んで頂きありがとうございます。衛星データではよく聞くけど、深くは理解していない方は多いのではないでしょうか?正直、相当な技術者、科学者か、宇宙変態野郎にしかあまり役に立たない情報かもしれません。少しでも皆さんの一助になればと思います。

また、書籍のご予約は以下からできます。
https://note.com/emmyeil/n/n2eba0e7df4f3

残りはおまけの内容になります。

自己紹介

私は、安井 秀輔(ヤスイ シュウスケ)と申します。日頃は、衛星データの分析を生業にしています。
X(旧Twitter)アカウントでは、宇宙領域や機械学習などの科学やコンペなどについて発言することが多いです。

SAR解析をよくやっていますが、画像系AI、地理空間や衛星データ、点群データ、3Dデータに関心があります。勉強している人は好きなので楽しく絡んでくれると嬉しいです。

データ分析やPython実装を中心に個人で記事を執筆しています。お役に立てられると幸いです。
https://zenn.dev/syu_tan

衛星データ解析として、宙畑のライターもしています。
https://sorabatake.jp/?s=秀輔

参考文献

Discussion