How to read a binary data produced from FORTRAN

fortranで得られたﾊﾞｲﾅﾘﾃﾞｰﾀをpythonで読み込む方法を説明する。
重要なのはf90ﾌｧｲﾙのopenのところにaccess="stream"を加える。また、status='replace'にする。 ﾊﾞｲﾅﾘﾌｧｲﾙで出力すると効率的に読み込める
ｻﾝﾌﾟﾙﾌﾟﾛｸﾞﾗﾑを以下に示す。(なぜかcatだけだと表示されないので、duも書いておく。)

In [1]:
!cat sin_curve.f90
!gfortran sin_curve.f90
!./a.out
!du -h sin_curve.d

program sin_curve
  implicit none
  integer,parameter :: ng = 1000
  integer i,ig
  real(8) :: x(ng),pi

  open(17,file="sin_curve.d",status="replace",form="unformatted" ,access="stream")
  pi=3.1415926d0
  do i = 1,100
    write(17) (sin(2.d0*pi*ig/ng+i*0.1d0),ig=1,ng)
    !print *, x
  enddo
  
end program sin_curve
393M	sin_curve.d

pythonで読み込むときにはnumpyのfromfileを使うと高速。

In [2]:
import numpy as np
f=open('sin_curve.d','rb')
%time data=np.fromfile(f,dtype='float64').reshape(-1,1000)
np.shape(data)

CPU times: user 15.6 ms, sys: 0 ns, total: 15.6 ms
Wall time: 915 µs
Out[2]:
(100, 1000)
In [3]:
import matplotlib.pyplot as plt
plt.imshow(data,aspect='auto',origin='lower');
plt.show()

書式ありﾌｧｲﾙで出力した時の例を以下で示す。
ｻﾝﾌﾟﾙﾌﾟﾛｸﾞﾗﾑは以下。

In [4]:
!cat sin_curve2.f90
!gfortran sin_curve2.f90
!./a.out
!du -h sin_curve2.d

program sin_curve
  implicit none
  integer,parameter :: ng = 10000
  integer i,ig
  real(8) :: x(ng),pi

  open(17,file="sin_curve2.d",status="replace",form="formatted")
  pi=3.1415926d0
  do i = 1,1000
    write(17,'(10000e13.5)') (sin(2.d0*pi*ig/ng+i*0.1d0),ig=1,ng)
    !print *, x
  enddo
  
end program sin_curve
62G	sin_curve2.d
In [5]:
import numpy as np
f=open('sin_curve2.d','r')
%time data2=np.fromfile(f,dtype='float64',sep=" ").reshape(-1,10000)
np.shape(data2)

CPU times: user 5.5 s, sys: 4.08 s, total: 9.58 s
Wall time: 9.63 s
Out[5]:
(1000, 10000)
In [6]:
import matplotlib.pyplot as plt
plt.imshow(data2,aspect='auto',origin='lower');
plt.show()

書式ありﾌｧｲﾙの場合はpd.read_csvを使うと高速。

In [7]:
import pandas as pd
f=open('sin_curve2.d','r')
%time data3=np.array(pd.read_csv(f,delim_whitespace=True, header=None,dtype='float64'))
np.shape(data3)

CPU times: user 3.31 s, sys: 391 ms, total: 3.7 s
Wall time: 3.76 s
Out[7]:
(1000, 10000)
In [8]:
import matplotlib.pyplot as plt
plt.imshow(data3,aspect='auto',origin='lower');
plt.show()