pythonとfortranのI/O関連まとめ

はじめに

python<->fortranの入出力について、インターネット上でまとまった記述が発見しにくいので、わかっている範囲でまとめることにした。

余談だが、python同士のデータ入出力にはnp.savezとnp.loadがおすすめです。

pythonの出力をfortranで読み込む

numpyを用いる方法

以下のような手順でデータを出力する
1. numpy配列を準備
2. reshapeでfortranに合うように変換
3. astypeでバイナリに変換
4. tofileで出力

3次元の倍精度実数をlittle endianで出力するpythonのコードの例

import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx))
a.reshape([ix*jx*kx],order='F').astype(endian+'d').tofile('test.dat')

これをfortranで読み込むには以下のようにする。使っている計算機がlittle endianだとする。fortranでのopenにはaccess='stream'が必要となる

program test
  implicit none

  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  read(idf) a
  close(idf)

  stop
end program test

big endianのデータをやりとりしたいときはpythonのコードで

endian='>'

と変えれば良い。このendian変数はastypeでバイナリへの変換方法を制御している

fortranコードでは

  open(newunit=idf, file='test.dat',form='unformatted',access='stream',convert='big_endian')

とする。単精度で出力・入力したいときはpythonコードを以下のように書き換える

import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx),dtype=np.float32)
a.reshape([ix*jx*kx],order='F').astype(endian+'f').tofile('test.dat')

np.onesとastypeの部分を単精度を扱えるように変更した。

対応するfortranコードは以下のようになる

program test
  implicit none

  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.e0)), dimension(ix,jx,kx) :: a

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  read(idf) a
  close(idf)

  stop
end program test

real(KIND(0.e0))の部分のみを変更した。

複数の配列を一つのファイルに出力したい場合は以下の手順を踏む
1. 複数の配列をnp.arrayでまとめる
2. reshapeで１次元配列とする。
2. Tでインデックスの位置を入れ替える。

import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx),dtype=np.float32)
b = np.ones((ix,jx,kx),dtype=np.float32)*2
c = np.ones((ix,jx,kx),dtype=np.float32)*3
np.array([a,b,c]).T.reshape([3*ix*jx*kx],order='F').astype(endian+'f').tofile('test.dat')

読み込むfortranコードは

program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.e0)), dimension(ix,jx,kx) :: a,b,c

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  read(idf) a,b,c
  close(idf)

  stop
end program test

とすれば良い

一方、大きさの違う配列を出力するときは
1. reshapeを用いて1次元配列とする
2. np.concatenateを用いて配列をまとめる
という手順をふめば良い

pythonコードは以下のようにする

import numpy as np

ix, jx, kx = 3, 4, 5
lx, mx = 2, 4
endian='<'
a = np.ones((ix,jx,kx))
b = np.ones((lx,mx))*2

a1 = a.reshape([ix*jx*kx],order='F')
b1 = b.reshape([lx*mx],order='F')
np.concatenate([a1,b1]).astype(endian+'d').tofile('test.dat')

fortranから読み込むときはこれまでと変わらないが以下のようにする。

program test
  implicit none

  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  integer, parameter :: lx = 2, mx = 4
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a
  real(KIND(0.d0)), dimension(lx,mx) :: b

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  read(idf) a,b
  close(idf)

  stop
end program test

scipy.io.FortranFileを用いる方法

上記の方法は、pythonの方のコードは簡単だが、fortranの方でaccess='stream'を使わないと自分でデータのヘッダをうまく編集してやらないといけなくなる。
この部分をうまく調整してくれるのがscipy.io.FortranFileである。

Pythonコードは以下のようにする

import numpy as np
from scipy.io import FortranFile

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx))
f = FortranFile('test.dat','w')
f.write_record(a)
f.close()

読み込みのためのFortranコードは

program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a

  open(newunit=idf, file='test.dat',form='unformatted')
  read(idf) a
  close(idf)

  stop
end program test

とすれば良い。

複数の配列を出力するときはNumpyを用いる時よりも容易で

a = np.ones((ix,jx,kx))
b = np.ones((ix,jx,kx))*2
f = FortranFile('test.dat','w')
f.write_record(a,b)
f.close()

とwrite_recordの引数に付け足すのみで良い。
fortranの読み込みコードも

program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a,b

  open(newunit=idf, file='test.dat',form='unformatted')
  read(idf) a,b
  close(idf)

  stop
end program test

と読み込む配列を横に並べればいいだけである。

とはいえ、scipyのウェブサイトを見ると

Since this is a non-standard file format, whose contents depend on the compiler and the endianness of the machine, caution is advised. Files from gfortran 4.8.0 and gfortran 4.1.2 on x86_64 are known to work.

Consider using Fortran direct-access files or files from the newer Stream I/O, which can be easily read by numpy.fromfile.

と書かれている。fortranのフォーマットはコンパイラ依存が大きいらしくnumpyを用いた方法が推奨されるようである。

fortranの出力をpythonで読み込む

numpyを用いる方法

numpyでfortranのデータを読み込むためにはaccess='stream'で出力するのがおすすめである(ヘッダなどを自分で調節すればいいので必須ではない)。

fortranで配列を定義して出力する

program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a

  a = 1.d0

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  write(idf) a
  close(idf)

  stop
end program test

pythonで読み込むには最も単純には、np.fromfile を用いれば良い。

import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
f = open('test.dat','rb')
a = np.fromfile(f,endian+'d',ix*jx*kx)
a = a.reshape((ix,jx,kx),order='F')
f.close()

np.fromfile を用いる場合でもnp.dtypeを使うとより汎用性の高い読み込みが可能になる。例えば以下のように行う

import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
f = open('test.dat','rb')
dtype = np.dtype([('a',endian+str(ix*jx*kx)+'d')])
a = np.fromfile(f,dtype=dtype)['a']
a = a.reshape((ix,jx,kx),order='F')
f.close()

複数の配列を取り扱うときはnp.dtypeを用いる方が便利である。
まずfortranで以下のように複数の配列を一つのファイルに出力する。

program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a,b

  a = 1.d0
  b = 2.d0

  open(newunit=idf, file='test.dat',form='unformatted',access='stream')
  write(idf) a,b
  close(idf)

  stop
end program test

これをpythonで読み込むには以下のようにする

import numpy as np

ix, jx, kx = 3, 4, 5
endian='<'
f = open('test.dat','rb')
dtype = np.dtype([('a',endian+str(ix*jx*kx)+'d'),('b',endian+str(ix*jx*kx)+'d')])
data = np.fromfile(f,dtype=dtype)
a = data['a'].reshape((ix,jx,kx),order='F')
b = data['b'].reshape((ix,jx,kx),order='F')
f.close()

scipy.io.FortranFileを用いる方法

access='stream'を用いていないfortranデータを読み込むためにはscipy.io.FortranFileを用いるのが容易である。

fortranの出力は以下のようにする

program test
  implicit none

  integer :: i,j,k
  integer :: idf
  integer, parameter :: ix = 3, jx = 4, kx = 5
  real(KIND(0.d0)), dimension(ix,jx,kx) :: a

  a = 1.d0

  open(newunit=idf, file='test.dat',form='unformatted')
  write(idf) a
  close(idf)

  stop
end program test

これを読み込むpythonコードは以下のようにする

import numpy as np
from scipy.io import FortranFile

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx))
b = np.ones((ix,jx,kx))*2
f = FortranFile('test.dat','r')
a = f.read_record(endian+'('+str(ix)+','+str(jx)+','+str(kx)+')d')
f.close()

np.dtypeを用いて以下のようにしても読み込むことが出来る

import numpy as np
from scipy.io import FortranFile

ix, jx, kx = 3, 4, 5
endian='<'
a = np.ones((ix,jx,kx))
b = np.ones((ix,jx,kx))*2
f = FortranFile('test.dat','r')
dtype = np.dtype([('a',endian+str(ix*jx*kx)+'d')])
a = f.read_record(dtype=dtype)
a = a['a'].reshape(ix,jx,kx,order='F')
f.close()