Fortranのバイナリ出力データ(unformatted)をPythonで読む

Fortranで出力したunformatted形式のデータをPythonで読む方法です。 scipyにFortranFileというものがあるので、それを使ってもできます。

scipy.io.FortranFile — SciPy v0.15.0.dev Reference Guide

でもこれソース読めばわかりますが、numpyのfromfileって関数を利用してるだけです。

numpy.fromfile — NumPy v1.8 Manual

fromfileで直接扱ったほうが応用も効くのでこれでやります。環境依存な部分もあるかもしれませんので注意してください。

今回読み込むデータは次のようなFortranソースコードで出力されたとします。

program main
    integer :: a, i,j
    real(8) :: b, arr(4,4)
    a = 4
    b = 3.14
    forall(i=1:4,j=1:4)
        arr(i,j) = 10*i+j
    end forall
    open(20, file="output.dat", form="unformatted")
    write(20) a, b, arr
    close(20)
end program

aは32bit整数、bは64bit浮動小数点数、arrは4x4の2次元配列です。こうすると、output.datというファイルには

という感じでデータが並びます。データの先頭と末尾にwrite文で出力したデータサイズ、今回の場合は

変数aのサイズ(4)+bのサイズ(8)+arrのサイズ(4x4x8) = 140バイト

がSize information(正式な呼び方は知りません)と書いたところに入ります。この領域自体は4バイトです。ファイルに追記していく場合でもwrite分のたびに、サイズが出力されます。

Python側ですが、fromfileのドキュメントを読めば予想がつくかもしれません。こんなかんじになります。

import numpy as np

head = ("head","<i")
tail = ("tail","<i")

dt = np.dtype([head, ("a","<i"), ("b","<d"), ("arr","<16d"), tail])
fd = open("output.dat","r")
chunk = np.fromfile(fd, dtype=dt, count=1)
print "a:",chunk[0]["a"]
print "b:",chunk[0]["b"]
arr = chunk[0]["arr"].reshape((4,4),order="F")
print "arr:",arr

dtypeでデータ構造を指定しています。"<"はリトルエンディアンです。ビッグエンディアン環境で出力した場合は代わりに">"を使います。つまり

dtype	意味
<i	リトルエンディアンの整数
<d	リトルエンディアンの浮動小数点数
<16d	リトルエンディアンの浮動小数点数16個

という意味です。2次元配列とかの指定方法はないので、配列の総数を指定します。データを読み込んだ後、reshapeを使って2次元配列に整形します。order="F"がポイントです。いわゆるC言語系とFortranでは多次元配列の並び方が違うので。

fromfileでcount=1としているのは読み込みをループしないという意味です。今回はwrite文が1回しかないので書かなくても同じですが。複数回write文を呼び出した場合はcountを読み込みたい回数に変更します。chunk[1]に2回目のwrite文で出力したデータが入ります。

これでFortranのデータは読めるようになったので、Pythonでデータを解析してmatplotlibで可視化するとか簡単にできますね！