5

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

我使用我希望能够在 python 中读取的 fortran 将二维数组写入未格式化的直接访问文件。下面给出了确切的测试代码。实际代码是一个巨大的并行代码，但数据输出的格式几乎完全相同。

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

如果时间片等于 i，我希望 def readslice 能够在第 i 个位置读取记录。为此，我尝试使用 f.seek 但它似乎不起作用。numpy.fromfile 似乎从第一条记录本身开始读取。如何使 numpy.fromfile 从文件中的特定点读取？

我仍在尝试习惯 Python 风格并深入研究文档。任何帮助和指示将不胜感激。

4

2 回答 2

6

这是一个适合您的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而不是 file object f。因此，np.fromfile创建了一个新的文件对象，而不是使用您想要的。这就是它一直从文件开头读取的原因。此外，f.seek将字节数作为参数，而不是元素数。我将其硬编码为 8 以适合您的数据，但如果您愿意，可以将其设为通用。此外，字段参数 inreadslice是多余的。希望这可以帮助。

于 2012-05-07T06:52:37.823 回答
1

我不认为所有的 FORTRAN 运行时都是一样的。有些帧记录，有些则没有，而且我不太相信那些记录帧的人都会以同样的方式做到这一点。当然，他们通常可以回读自己编写的记录，但是从一个 FORTRAN 运行时到另一个运行时的互操作可能不存在。

您可能应该在您选择的 FORTRAN 中编写一个小型测试程序，该程序会写入一些类似于您的生产代码的记录，然后使用您最喜欢的二进制文件编辑器挑选结果 - 我喜欢 bpe，但其中有很多可用。

然后，在您了解真正编写的内容之后，使用 Python 结构模块或类似模块将内容拉回。

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

于 2012-05-07T03:58:21.683 回答