Read direct access fortran unformatted file in Python

Question

I am completely new to Python and I write my Python rendering codes from scratch (to avoid using expensive proprietary programs like IDL). So far I have used IDL and gnuplot. What I want to do is the following:

I am writing two-dimensional arrays to files without formatted direct access using fortran, which I want to read in python. The exact test code is given below. The actual code is a huge parallel code, but the output is almost in the 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 at t = 1 and increase the number of "islands" to increase the value of t. But when I read it using the python code below, I just get zeros. If I delete the first operator to write zeros, I just get the first slice, no matter what the "timeslice" value I use in 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 def readslice to be able to read the record in i-th place if timelice is i. I tried using f.seek for this, but it does not seem to work. It seems that numpy.fromfile starts reading from the very first record. How to make numpy.fromfile to read from a specific point in a file?

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

Answer 1

Here is the 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 the source code, you passed the file name as the first argument to np.fromfile instead of the file object f . Thus, np.fromfile created a new file object instead of using the one you intended. It was for this reason that he continued to read from the very beginning of the file. In addition, f.seek takes the number of bytes as an argument, not the number of elements. I tied it tightly to 8 to fit your data, but you can do this general if you want. In addition, the field argument in readslice was redundant. Hope this helps.

Answer 2

I do not think that FORTRAN is the same all the time; some frame recordings, some not, and I'm not sure that those recording crop will do the same thing. They can usually read records written by themselves, but there may be no interop from one FORTRAN runtime to another.

You should probably write a small test program in your FORTRAN that records a couple of entries similar to your production code, and then highlights the results using your favorite binary editor - I like bpe, but there are a lot of them available.

Then, after you understand what is actually written, return things to use a Python structure module or similar.

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