python - 将 FORTRAN90 文件生成的数据读入 NUMPY 数组

Question

我有一个二进制文件想要读入 Python。该文件由三部分组成:波数列表、温度列表以及作为温度和压力函数的不透明度列表。我想将其中前两个导入为向量 a 和 b，将第三个导入为二维数组 c，这样 c[x,y] 对应于 a[x] 和 b[y]。

现有的 FORTRAN90 代码能够实现这一点，如下所示:

integer   nS, nT
parameter (nS = 3000)
parameter (nT = 9)

real(8) wn_arr(nS)      ! wavenumber [cm^-1]
real(8) temp_arr(nT)    ! temperature [K]
real(8) abs_arr(nS,nT)  ! absorption coefficient [cm^-1 / amagat^2]

open(33,file=trim(datadir)//'CO2_dimer_data',form='unformatted')
read(33) wn_arr
read(33) temp_arr
read(33) abs_arr
close(33)

我尝试了以下 python 代码:

f=scipy.io.FortranFile('file', 'r')
a_ref=f.read_reals(np.float64) #wavenumber (cm**-1)
b=f.read_reals(np.float64) #temperature (K)
c=f.read_reals(np.float64).reshape((3000,9))

但是，这会产生不正确的结果。我怀疑这是因为 Fortran 将数组写入文件的顺序与 Python 不同。但是，简单地将 order='F' 添加到 reshape 命令中是行不通的。我怀疑这是因为读入时，abscoeff_ref 已经变平了。

有什么想法吗？

Answer 1

为了让您了解我在第二条评论中的意思，我做了一个模型:

testwrite.f90，使用 gfortran 4.8.4 编译:它基本上会写入一个未格式化的顺序文件，其中包含您指定的数组(只是小得多，以便能够通过肉眼比较它们)，其中填充了任意数据。它还打印数组。

implicit none
integer   nS, nT, i ,j
parameter (nS = 10)
parameter (nT = 3)

real(8) wn_arr(nS)      ! wavenumber [cm^-1]
real(8) temp_arr(nT)    ! temperature [K]
real(8) abs_arr(nS,nT)  ! absorption coefficient [cm^-1 / amagat^2]

wn_arr = (/ (i, i=1,nS) /)
temp_arr = (/ (270+i, i=1,nT) /)
abs_arr = reshape( (/ ((10*j+i, i=1,nS), j=1,nT) /), (/nS, nT/))

print*, wn_arr
print*, '-----------------'
print*, temp_arr
print*, '-----------------'
print*, abs_arr
print*, '-----------------'
print*, 'abs_arr(5,3) = ', abs_arr(5,3)

open(33,file='test.out',form='unformatted')
write(33) wn_arr
write(33) temp_arr
write(33) abs_arr
close(33)

end

testread.py，使用 Python 2.7.6 进行测试，然后读取上面编写的文件并打印数组。对我来说，两个程序的输出是相同的。 YMMV。

import numpy as np

rec_delim = 4  # This value depends on the Fortran compiler
nS = 10
nT = 3

with open('test.out', 'rb') as infile:
    infile.seek(rec_delim, 1)  # begin record
    wn_arr = np.fromfile(file=infile, dtype=np.float64, count=nS)
    infile.seek(rec_delim, 1)  # end record
    infile.seek(rec_delim, 1)  # begin record
    temp_arr = np.fromfile(file=infile, dtype=np.float64, count=nT)
    infile.seek(rec_delim, 1)  # end record
    infile.seek(rec_delim, 1)  # begin record
    abs_arr = np.fromfile(file=infile,
                        dtype=np.float64).reshape((nS, nT), order='F')
    infile.seek(rec_delim, 1)  # end record

print(wn_arr)
print(temp_arr)
print(abs_arr)
# The array has the same shape, but Fortran starts index (per default at least)
# at 1 and Python at 0:
print('abs_arr(5,3) = ' + str(abs_arr[4,2]))

简短说明:我在 with-block 中打开文件(Python 中的良好实践)，然后利用文件编写方式的知识逐步浏览该文件。这使得它不可移植。 infile.seek(4, 1) 将Python的读取指针从当前位置(选项1)向前移动4个字节，因为我知道该文件以4个字节长的开始记录标记开始(gfortran) 。

然后我使用 numpy.fromfile 读取 float64 的 count=10 个值，这是波数数组。

接下来我必须跳过结束记录和开始记录标记。这当然也可以通过 infile.seel(8, 1) 来完成。

然后我读取温度数组，再次跳过结束记录和开始记录标记并读取二维数组。文件中的数据不知道它是二维的，因此我需要使用 Fortran 顺序对其进行 reshape 。最后一个 .seek() 是虚假的，我只是想强调结构。

我再次强烈建议您不要在这样的代码上构建更大的系统。对于一次性来说还好，但对于必须再次使用或共享的东西来说就很糟糕了。