在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代码中使用什么“
timeslice”值，我都将获得第一个时间片。到目前为止，我的代码是：

#!/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能够在timelice等于i的情况下读取第i个位置的记录。为此，我尝试使用f.seek，但它似乎不起作用。numpy.fromfile似乎从第一条记录本身开始读取。如何使numpy.fromfile从文件的特定点读取？

我仍在尝试适应Python风格并仔细阅读文档。任何帮助和指针将不胜感激。

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”中的参数是多余的。希望这有帮助。