Reading a direct access fortran unformatted file in Python

file, fortran, python

Solution

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.

Problem

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.

Original source