Reading a direct access fortran unformatted file in Python

Question

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.

Answer 1

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.