将从FORTRAN90文件生成的数据读入NUMPY数组 [英] Reading data generated from FORTRAN90 file into NUMPY array

问题描述

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

现有的FORTRAN90代码可以实现这一点，如下所示：

 整数nS，nT 
参数（nS = 3000）
参数（nT = 9）
 
实数（8）wn_arr（nS）！波数[cm ^ -1] 
 real（8）temp_arr（nT）！温度[K] 
实数（8）abs_arr（nS，nT）！吸收系数[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'添加到整形命令中不起作用。我怀疑这是因为在阅读时，abscoeff_ref已经变平。

有什么想法？

解决方案

为了让你知道我在第二条评论中的含义，我做了一个模拟：

testwrite.f90，用gfortran 4.8编译.4：它基本上用你指定的数组写入一个未格式化的顺序文件（只是小得多，能够通过眼睛来比较它们），并填充任意数据。
它也打印数组。

 隐式无
整数nS，nT，i，j 
参数（nS = 10）
参数（nT = 3）
 
 real（8）wn_arr（nS）！波数[cm ^ -1] 
 real（8）temp_arr（nT）！温度[K] 
实数（8）abs_arr（nS，nT）！吸收系数[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 =' （33）
写入（33）wn_arr 
写入（33）temp_arr 
写入（33）abs_arr 
关闭（33）
 
结束
  

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

  import numpy as np 
 
 rec_delim = 4＃此值取决于Fortran编译器
 nS = 10 
 nT = 3 
 
 with open（'test.out'，'rb'）as infile：
 infile.seek（rec_delim，1）＃开始记录
 wn_arr = np.fromfile（file = infile，dtype = np.float64，count = nS）
 infile.seek（rec_delim，1）＃结束记录
 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）＃结束记录
 
 print（wn_arr）
 print（temp_arr）
 print（abs_arr）
＃数组有相同的形状，但Fortran启动索引（至少默认）
＃在1和Python在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读取count = 10的float64数值，这是波数数组。

接下来，我必须跳过结束记录和开始记录标记。这当然也可以通过infile.seel（8,1）来完成。然后我读取温度数组，跳过结束记录并重新开始记录标记并读取2D数组。
文件中的数据不知道它是2D的，所以我需要使用Fortran顺序对其进行重新设计。
最后一个.seek（）是虚假的，我只是想强调结构。

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

I have a binary file I want to read into Python. The file consists of three parts: a list of wavenumbers, a list of temperatures, and a list of opacities as a function of temperature and pressure. I would like to import the first two of these as vectors a and b, and the third as a 2D array c such that c[x,y] corresponds to a[x] and b[y].

An existing FORTRAN90 code is able to achieve this as follows:

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)

I tried the following python code:

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))

However, this generates results that are not correct. I suspect this is because Fortran writes arrays to file in a different order than Python does. However, simply adding order='F' to the reshape command does not work. I suspect this is because when read in, abscoeff_ref is already flattened.

Any thoughts?

解决方案

To give you an idea of what I meant in my second comment, I made a mock up:

testwrite.f90, compiled with gfortran 4.8.4: It basically writes an unformated, sequential file with the arrays you specified (just a lot smaller, to be able to compare them by eye) filled with arbitrary data. It also prints the arrays.

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, tested with Python 2.7.6, then reads the file written above and prints the arrays as well. For me the output of both programs is the same. 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]))

Short explanation: I open the file in a with-block (good practice in Python) and then I step through the file, using the knowledge of how the file is written. This makes it unportable. infile.seek(4, 1) moves the read-pointer of Python forward 4 bytes, from the current position (the option 1), because I know that the file starts with a begin-record marker that is 4 bytes long (gfortran).

Then I use numpy.fromfile to read count=10 values of float64, which is the wavenumber array.

Next I have to skip the end-record and begin-record markers. This could of course be done by infile.seel(8, 1) as well.

Then I read the temperature array, skip end-record and begin-record markers again and read the 2D-array. The data in the file doesn't know that it's 2D, so I need to reshape it, using the Fortran-order. The last .seek() is spurious, I just wanted to emphasize the structure.

I strongly suggest again that you do not build a bigger system on code like this. It's fine for a one-off, but terrible for something you have to use again or share.

这篇关于将从FORTRAN90文件生成的数据读入NUMPY数组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！