将 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 块中打开文件（Python 中的良好做法），然后使用文件编写方式的知识逐步浏览文件。这使得它不便携。 infile.seek(4, 1) 将 Python 的读取指针从当前位置（选项 1）向前移动 4 个字节，因为我知道文件以 4 个字节长的开始记录标记开始（gfortran） .

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

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

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

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