在 Python 中读取直接访问 fortran 无格式文件

问题描述

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

I'm completely new to Python and am writing my visualization codes in Python from scratch (to avoid using expensive proprietary programs like IDL). Until now I've used IDL and gnuplot. What I want to be able to do is this:

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

I write two dimensional arrays to unformatted direct access files using fortran which I want to be able to read in python. The exact test code is given below. The actual code is a huge parallel code but the data output is almost the exact same format.

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 代码中使用的时间片"的值是什么.我到目前为止的代码是:

The program should give me zeros for t=1 and increasing number of "islands" for increasing value of t. But when I read it using python code given below, I just get zeros. If I remove the first write statement of zeros, I just get the first time slice irrespective of what value of "timeslice" I use in the python code. The code I have so far is:

#!/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 从文件中的特定点读取?

I want the def readslice to be able to read a record at the i-th place if timeslice equals i. For that I've tried to use f.seek but it does not seem to work. numpy.fromfile seems to start reading from the first record itself. How do I make numpy.fromfile to read from a specific point in a file?

我仍在努力适应 Python 风格并深入研究文档.任何帮助和指点将不胜感激.

I'm still trying to get used to the Python style and digging through the documentation. Any help and pointers would be greatly appreciated.

推荐答案

下面是一个适合你的python代码:

Here is a python code that will work for you:

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 中的字段参数是多余的.希望这会有所帮助.

In your original code, you were passing file name as first argument to np.fromfile instead of file object f. Thus, np.fromfile created a new file object instead of using the one that you intended. This is the reason why it kept reading from the beginning of the file. In addition, f.seek takes the number of bytes as argument, not the number of elements. I hardcoded it to 8 to fit your data, but you can make this general if you wish. Also, field argument in readslice was redundant. Hope this helps.

这篇关于在 Python 中读取直接访问 fortran 无格式文件的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持html5模板网！