用Python读取一个直接访问fortran未格式化的文件 [英] Reading a direct access fortran unformatted file in Python

问题描述

我对Python完全陌生，正在从头开始编写我的可视化代码（避免使用像IDL这样昂贵的专有程序）。直到现在我已经使用IDL和gnuplot。我希望能够做到这一点：$ b​​
$ b

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

 程序binary_out 
隐式无
整型:: 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代码中使用的时间片的值。我到目前为止的代码是：
$ b $ pre $ ＃！/ usr / bin / env python
导入scipy
import glob
import numpy as np
import matplotlib.pyplot as plt
import py，sys
from pylab import *

def readslice（inputfilename，字段，nx，ny，时间片）：
f = open（inputfilename，'r'）
f.seek（timeslice * nx * ny）
field = np.fromfile（inputfilename，dtype = 'd'，count = nx * ny）
field = np.reshape（field，（nx，ny））
返回字段

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风格，并通过文档进行挖掘。任何帮助和指针将不胜感激。

解决方案

这是一个可以帮助您的python代码：

  def readslice（输入文件名，nx，ny，时间片）：
f = open（输入文件名，'rb'）
 f。 seek（8 * timeslice * nx * ny）
 field = np.fromfile（f，dtype ='float64'，count = nx * ny）
 field = np.reshape（field，（nx，ny ））
 f.close（）
返回字段
  

在原始代码，您将文件名作为第一个参数传递给 np.fromfile 而不是文件对象 f 。因此， np.fromfile 创建了一个新的文件对象，而不是使用您想要的文件对象。这就是为什么它从文件的开始处继续阅读的原因。另外， f.seek 将字节数作为参数，而不是元素数量。我将它硬编码为8来适合您的数据，但是如果您愿意，您可以将其设为一般。此外， readslice 中的字段参数是多余的。希望这有助于。

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:

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

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

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

解决方案

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

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未格式化的文件的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！