Synthetic Aperture Radar （SAR） データ（バイナリ形式）の読み込み

はじめに

衛星画像のIoTノードの一つとしての利用検討が進んでいます。
ここでは、衛星搭載型の合成開口レーダ（英：Synthetic Aperture Radar）（略：SAR）データの読み込みについて書きたいと思います。
はじめにぶっちゃけますと、専用のソフトがあるのでそれを使うのが一番早いし確実です。例えば以下のようなものがあります。
- PolSAR pro
- Sentinel Toolboxes
なのですが、アカデミック用途に限って言えば、こうした既存のソフトがデフォルトで読み込んでくる以外のメタデータが必要になることもあるでしょう（知らんけど）。
そもそもこの記事を書くきっかけは、「自分で読み込みたいんです」というニーズがあることを知ったためなので、その前提でバイナリデータから読み込む方法について書いていきたいと思います。上記のソフトを使って開きたいという方はそれぞれのWebページ見ていただければ、マニュアルを見つけることができるものと思います。

想定読者

バイナリ形式のデータについての理解がある方
データ型に理解がある方
MATLAB（Pythonもほぼ一緒）触ったことある方

データセット

この記事では国産のSAR衛星であるPALSAR-2データを使います。一口にPALSAR-2データと言っても色々バリエーションがあるのですが、この記事ではCEOSデータフォーマットのL1.1、単偏波画像とします。なぜこれを選んだかというと、自分でコードを書く、という想定と照らし合わせて最も現実的だからです（それでもソフト使うことを強く勧めますが）。
SARデータには、他にもTerraSARとかSentinelとか、最近だと小型SARもちらほら出てきています。それぞれのデータにはデータフォーマットの説明書があるので、フォーマットを理解して読み込みたい部分のbyteは何番目から何番目なのかを理解して、そこを読む、という手順は変わりません。
今日は一例としてPALSAR-2、ということです。

ダウンロード元はこのページです。
GMTSARのダウンロードページ
このページの"ALOS-2 L1.1"をいただきました。

PALSAR-2データの構成

参考文献はこれです。
PALSAR-2プロダクトフォーマット説明書
こちらによるとPALSAR-2のCEOSフォーマットでは、次の4つのファイルが1セットになっています。

CEOSレベル1.1/1.5/3.1データは、CEOSスーパーストラクチャフォーマットに準拠した複数ファイルから構成される。1シーンのCEOSレベル1.1/1.5/3.1データの基本的なファイル構成は、ボリュームディレクトリファイル、SARリーダファイル、SARイメージファイル及びSARトレイラファイルである。

それぞれのファイルは頭文字でVOL, LED, IMG, TLRで始まります。それぞれのフォーマットが説明書に書いています。以下ではメインで読み込むことになるだろうイメージファイルに絞って説明していきます。
ちなみに上のダウンロードデータでは、LEDとIMGのみ含まれます。

SAR IMGファイル

31ページによるとイメージファイルは次の2パートで構成されています

レコード番号	レコード長	レコード数	レコード名
1	720	1	ファイルディスクリプタ
2~n+1	可変長	n	シグナルデータ

ファイルディスクリプタの後にシグナルデータがあるということですね。
83ページ以降に各パートのデータフォーマットがバイト単位で記載されています。これを見ると例えば、ファイルディスクリプタの1バイト目には、レコード番号で1という数字が記録されている、ということがわかります。そしてそのデータ型は"B4"となっています。
さて、この"B4"とは何のことでしょう。これが一般名称なのかはわかりませんが、やはり説明書に記載があります。32ページですね。
この表3.2-2によるとタイプは次の表に示す通り定義されています。

タイプ	詳細
Am	キャラクタ表示（特に指定がない場合、左詰め）
Im	整数を表現するASCII文字列（右詰め）
Fm.n	実数タイプデータ表示（右詰め）
Em.n	実数タイプデータ表示（指数表現、右詰め）
Bm	2進数表示（1番目が最上位のバイト、ビッグエンディアン）

m：表示桁数
n：小数点以下の桁数
p：指数における乗数

となっています。つまり"B4"は、4桁の2進数表示でビッグエンディアンで記録されている、ということですね。つまり"0001"みたいにゼロイチが4つ並んでるということか、と思ったのですが、同じフィールドのバイトNo.は1-4ということなのでデータとしては4バイト分であることがわかります。ということは32ビットということなので、ゼロイチは32桁並ぶはず。よくわかりません。
仕方ないので、タイプごとに読み方を検証していくことにします。

タイプごとに読み方を特定していく

Am

まずはAm。
たまたま違う型でも読めたってことがないように、ある程度複雑な内容のデータを読んでみることにします。対象はフィールドNo.9。17－28バイト目。これには"CEOS-SARbbbb"(bは空白の意)と記録されている模様。
matlabではバイナリデータの読み込みはfopenでファイルを開いて、freadで読み込みます。任意のバイトまで移動するにはfseek,ファイルを閉じるのはfcloseです。このへんはだいたいどの言語でも似てますよね。
Amはキャラクタ表示ってことなので、データ型をcharと仮定して読んでみます。

fid = fopen(IMGのファイル名, 'rb', 'ieee-be'); % ビッグエンディアン指定
fseek(fid, 16, 'bof'); % 17バイト目まで移動
fread(fid, [1,12], 'char') % char型で12個読む
fclose(fid);

ans =
    67    69    79    83    45    83    65    82    32    32    32    32

あれっ！へんな数列出てきた。どうやら間違えた。
しかし、よく見ると、ASCIIコードに対応しているということがわかる（32はスペース）。
しかも出てきたデータの型はdoubleになっていて、文字列ではない。バイナリを前提として開いたからこうなってしまうのかな。よくわかりませんが、数値をUnicode文字列に変換するにはchar関数ということなので、さらにcharを噛ませてみる。

char(fread(fid, [1,12], 'char'))

ans =
    'CEOS-SAR    '

お～っ、出てきた出てきた。ほかのでも試してみましたが、Amタイプについては、これでよさげです。

Im

Imについては、整数を表現するASCII文字列ということなので、Amと同じでいいんじゃないだろうか、ということで実験しました。
ターゲットはフィールドNo.32。ちゃんと読めれば" 32"というデータが読めるはず。

fseek(fid, 216, 'bof');
char(fread(fid, [1,4], 'char'))

ans =
    '  32'

想定通りのデータが出てきたので。これで良しとします。

Fm.n

IMGファイルにはFm.nタイプとEm.nタイプで記録されたレコードがないので、LEDファイルを使って検証していくことにします。
LEDファイルのデータセットサマリレコードのフィールドNo.17を読んでみます。正しく読めれば6378.1370000という数値となるはずです。
Fというアルファベットと小数点というキーワードから浮動小数点型で記録されているのかなと思いましたが、49ページを見るとバイト数が16、桁数も16ということなので、やはり1バイトに1桁と考えられます。まぁとりあえずAm, Imと一緒のやり方で読んでみようか。

fid = fopen(LEDファイル名, 'rb', 'ieee-be');
offset = 720;
fseek(fid, offset+180, 'bof');
char(fread(fid, [1,16], 'char'))
fclose(fid);

ans =
    '    6378.1370000'

読めた！え～！書き方紛らわしい。これ全部同じで読めるんですかね。

Em.n

もうこれもおんなじでしょ。同じデータセットサマリレコードのフィールドNo.138をターゲットにしてみます。入射角の情報が入っているはず。

fseek(fid, offset+1886, 'bof');
char(fread(fid, [1,20], 'char'))

ans =
    ' 1.0940997542227E+00'

答えがわかりませんがそれっぽい値になっています。小数点が13桁というのにも一致しています。はい次っ！

Bm

IMGファイルのディスクリプタに戻ります。
フィールドNo.1から見ていくと、バイト数が1なのに2桁のデータが格納されていたりするのがわかります。ということはこれまでの1バイトに1桁が対応する読み方だとNGだということになります。
例えばフィールドNo.2。1バイトで記録されているデータは50。普通に考えたら(u)int8型でしょう。試してみます。

fseek(fid, 4, 'bof');
fread(fid, [1,1], 'int8')

ans =
    50

ごち！
しかし、さらに次のフィールドのNo.3を見てみると、1バイトで192という値が記録されているようです。int8型では127～‐128しか記録できないので、当然これではいけないということになります。試しにやってみましょう。

fseek(fid, 5, 'bof');
fread(fid, [1,1], 'int8')

ans =
    -64

そうですね。なのでuint8で読んであげます。

fseek(fid, 5, 'bof');
fread(fid, [1,1], 'uint8')

ans =
    192

正解が出てきました。
Bmについては、これでいいのかなと思ってパラパラ説明書を見ると、B2とかB4、B8もありますね。それぞれuint16とかuint32, uint64とかでいけそうです。
該当するところを試しに読んでみます。

offset = 720; % ファイルディスクリプタ分
fseek(fid, offset + 66, 'bof'); % フィールドNo.23
fread(fid, [1,1], 'uint16') % B2 チャープ形式指定子
fseek(fid, offset + 68, 'bof'); % フィールドNo.24
fread(fid, [1,1], 'uint32') % B4 チャープ長（パルス幅）[nsec]
fseek(fid, offset + 84, 'bof'); % フィールドNo.28
fread(fid, [1,1], 'uint64') % B8 センサ取得マイクロ秒（日内通算）

ans =
     0
ans =
       29010
ans =
   1.7701e+10

チャープ形式指定子は、説明書と一致します。
チャープ長は29μsecなのでまぁよさげ。
センサ取得秒は日内通算なのでだいたい5時くらいの計算になる。使ったデータはブラジル観測で時差はUTC-3～-5時間。PALSAR-2は現地時間で0時か12時付近の観測なので、整合はとれてそうです。

ということで無事に全タイプのフィールドを読めました。なんか例外もありそうですが、とりあえずこんな感じでデータを読んでくることはできそうです。

肝心の画像データは

さて、本丸の画像データの読み込みですが、93ページですね。フィールドNo.60の後、このレコードのピクセル数分のデータが記録されています。各ピクセルは8バイト。
ピクセル数は、フィールドNo.10に記録されています。

offset = 720; % ファイルディスクリプタ分
fseek(fid, offset + 24, 'bof'); % フィールドNo.10
fread(fid, [1,1], 'uint32') % B4

ans =
        9072

レコードは全部でいくつあるかというと、ファイルディスクリプタのフィールドNo.29に記録されています。

fseek(fid, 180, 'bof'); % フィールドNo.29
char(fread(fid, [1,6], 'char')) % レコード数

ans =
    ' 23292'

ということで23292回、9072ピクセル数分の8バイトデータを読めばいいということになります。これはシンプルにforループでできそうです。シグナルデータの先頭部分（画像データ以外のヘッダ部分）もforループに含める必要がある点に注意してください。
なお、L1.1データですので、1ピクセルにおけるデータは複素数となっています。1ピクセルは8バイトなので、実部と虚部はそれぞれ4バイトとなります。また、それぞれ正負をとりうるのでデータ型はint32ということになります。
ということで、以上を考慮して、最低限、画像部分だけ読み込むためには次のコードとなります。

fid = fopen(filename, 'rb', 'ieee-be');

%% レコード数の取得
fseek(fid, 180, 'bof'); % フィールドNo.29
recordNum = str2num(char(fread(fid, [1,6], 'char'))); % レコード数
fseek(fid, 720, 'bof'); % シグナルデータまで移動

%% レコード分だけピクセル数分読む
for record_i = 1:recordNum
    % ピクセル数の取得
    fseek(fid, 24, 'cof'); % フィールドNo.10
    pixelNum = fread(fid, [1,1], 'uint32'); % B4
    if record_i == 1 % データの格納用配列を用意
        imageData = zeros(recordNum, pixelNum);
    end
    fseek(fid, 544 - 28, 'cof'); % SARデータレコードまで移動
    lineData = fread(fid, [2 pixelNum], 'int32'); % 1レコード分のデータの読み込み
    imageData(record_i,:) = complex(lineData(1,:), lineData(2,:)); % 複素数に整形して格納 
end

fclose(fid);

出てきたimageDataをちょっと加工して表示します。

imageData = abs(imageData); % 強度
test = imresize(imageData, 0.25); % データが重いので縮小
imshow(log(test), []) % 表示

出てきた画像がこちら

あれ？！ノイズにしか見えない。全然読めてないじゃん！
なんでだと思ったらSARデータのフォーマットについては、87ページに記述があります。

'COMPLEX*8bbbbbbbbbbbbbbbbbbb''C*8b'：8バイトフィールド内前半分（4バイト）が2の補数表現。浮動小数点形式の実数を含み、後半分が虚数成分を含む複素表現。尚、無効データには0を格納する。

補数とか書いてるけど、浮動小数点と書いているので、int32をfloat32に替えて読んだ結果が次の画像です。

©JAXA

きちんと絵が出ています。これはアマゾン内のデータなので、こんな感じで良さそう！

Bmについては、(u)int型とは限らないのね。もしかして他にも浮動小数点型で読むべきデータがあるのかもしれませんが、基本的には上述のように説明書に記述があるはずなのでよく読めばわかりそうです。

さいごに

ということでSARデータの読出しを自分で書いてみました。掲載したコードが（最後以外）断片的ですみません。また、読めたことは読めたのですが、サンプル数は少なく本当にこれで全部のデータが読めるのかも未確認です。書いてある手順で読めなかったらすみません。
皆さんがそれぞれの目的で読み込みたいフィールドについては、説明書でバイト数やタイプを確認して、適宜fseek,freedの中を変えてトライしてみてください。

今後もSARデータ周りのhowtoを少しずつアップしていこうと思います。