用Python读取Fortran直接访问未格式化文件

Question

我完全不懂Python，想要从头开始使用Python编写我的可视化代码（以避免使用昂贵的专有程序如IDL）。到目前为止，我一直在使用IDL和gnuplot。我想要做到的是：

我使用Fortran将二维数组写入未格式化的直接访问文件中，希望能够在Python中读取它们。下面是精简版测试代码，实际代码是一个巨大的并行代码，但数据输出几乎与此相同格式。

program binary_out
implicit none
integer :: i,j,t,rec_array
double precision, dimension(100,100) :: fn
double precision, parameter :: p=2*3.1415929
INQUIRE(IOLENGTH=rec_array) fn
open(unit=10,file='test',status='new',form='unformatted',access='direct',recl=rec_array)                                                                           
   fn=0
   write(10,rec=1) fn
do t=1,3
do i=1,100
   do j=1,100
      fn(i,j)=sin(i*p*t/100)*cos(j*p*t/100)
   enddo
enddo
   write(10,rec=t+1) fn
enddo
close(10)
end program binary_out

程序应该对于t=1返回零，并对于逐渐增加的t值返回越来越多的“岛屿”。但是当我使用下面给出的Python代码读取它时，我只得到零。如果我删除零的第一个写操作，则无论我在Python代码中使用什么“时间片”值，我都只得到第一个时间片。我目前拥有的代码如下：

#!/usr/bin/env python
import scipy
import glob
import numpy as np
import matplotlib.pyplot as plt
import os, sys
from pylab import *

def readslice(inputfilename,field,nx,ny,timeslice):
   f=open(inputfilename,'r')
   f.seek(timeslice*nx*ny)
   field=np.fromfile(inputfilename,dtype='d',count=nx*ny)
   field=np.reshape(field,(nx,ny))
   return field

a=np.dtype('d')
a=readslice('test',a,100,100,2)

im=plt.imshow(a)
plt.show()

我希望 def readslice 能够在 timeslice 等于 i 时读取第i个位置的记录。 我尝试使用 f.seek，但似乎不起作用。numpy.fromfile 似乎从第一条记录开始读取。我如何使 numpy.fromfile 从文件中的特定点开始读取？ 我仍在努力适应 Python 风格并深入研究文档。任何帮助和指针将不胜感激。

- toylas

---

不确定您会做哪种可视化，但请注意VTK。 - Anycorn

我经常这样做（Fortran直接访问+Python文件查找+Numpy.fromfile）。你的代码看起来正确。尝试从Fortran程序中读取数据以进行合理性检查。此外，发布您用于编写数据的Fortran代码片段。确保您编写/读取的数据具有相同的精度/数据类型。您说您正在做的事情似乎不起作用。您能否提供一个示例，说明数据应该是什么样子以及您实际得到了什么？ - milancurcic

IRO-bot，能否请您分享一下读取这些文件的代码片段？这个代码似乎总是从第一条记录开始。np.fromfile的文档似乎没有关于数据读取起始点的关键字。 - toylas

我已经在我的OP中添加了精确的Fortran代码。因此，现在您甚至可以自己测试上述两个代码，以查看我所说的问题。谢谢帮助！！ - toylas

Answer 1

这里有一段Python代码，适用于您的需求：

def readslice(inputfilename,nx,ny,timeslice):
   f = open(inputfilename,'rb')
   f.seek(8*timeslice*nx*ny)
   field = np.fromfile(f,dtype='float64',count=nx*ny)
   field = np.reshape(field,(nx,ny))
   f.close()
   return field

在您的原始代码中，您将文件名作为np.fromfile的第一个参数传递，而不是文件对象f。因此，np.fromfile创建了一个新的文件对象，而不是使用您想要的那个对象。这就是它一直从文件开头读取的原因。此外，f.seek的参数是字节数，而不是元素数。我将其硬编码为8以适应您的数据，但如果您愿意，可以将其设置为通用值。另外，在readslice中的field参数是多余的。希望这能帮到您。

- milancurcic

---

非常感谢IRO-bot！！最终是我的愚蠢错误！！感谢您指出！现在一切都完美地运作了！！！ - toylas

完成！我试图+1，但由于我还没有“声望”，所以不允许我这样做:-p 点击了检查标记。我对这个网站还非常陌生。 - toylas

Answer 2

我认为并非所有的FORTRAN运行时都是相同的；有些框架记录，有些则没有，而且我不太确定那些记录框架的运行时是否都以相同的方式进行。当然，它们通常可以读取自己编写的记录，但从一个FORTRAN运行时到另一个运行时的互操作可能不存在。

您可能应该在您选择的FORTRAN中编写一个小型测试程序，编写几个类似于生产代码的记录，然后使用您喜欢的二进制文件编辑器（我喜欢bpe，但有许多其他可用的工具）来分析结果。

然后，在您理解实际编写的内容之后，使用Python struct模块或类似模块将其拉回。

bpe: http://sourceforge.net/projects/bpe/

- user1277476

---

感谢您的回答，但正如IRO-bot上面指出的那样，直接访问文件没有任何记录信息。我一直在IDL中读取这样的文件。这与编译器无关。我确实尝试过结构模块，但我还不够了解python。所以我无法实现它。上面的numpy.fromfile是最接近的实现，除了寻找问题。 - toylas