[Numpy-discussion] example reading binary Fortran file

[Neil Martinsen-Burrell]

On 2009-05-29 10:12 , David Froger wrote:
> I think the FortranFile class is not intended to read arrays written
> with the syntax 'write(11) array1, array2, array3'  (correct me if I'm
> wrong).  This is the use in the laboratory where I'm currently
> completing a phd.

You're half wrong.  FortranFile can read arrays written as above, but it 
sees them as a single real array.  So, with the attached Fortran program::

In [1]: from fortranfile import FortranFile

In [2]: f = FortranFile('uxuyp.bin', endian='<') # Original bug was 
incorrect byte order

In [3]: u = f.readReals()

In [4]: u.shape
Out[4]: (20,)

In [5]: u
Out[5]:
array([ 101.,  111.,  102.,  112.,  103.,  113.,  104.,  114.,  105.,
         115.,  201.,  211.,  202.,  212.,  203.,  213.,  204.,  214.,
         205.,  215.], dtype=float32)

In [6]: ux = u[:10].reshape(2,5); uy = u[10:].reshape(2,5)

In [7]: p = f.readReals().reshape(2,5)

In [8]: ux, uy, p
Out[8]:
(array([[ 101.,  111.,  102.,  112.,  103.],
        [ 113.,  104.,  114.,  105.,  115.]], dtype=float32),
  array([[ 201.,  211.,  202.,  212.,  203.],
        [ 213.,  204.,  214.,  205.,  215.]], dtype=float32),
  array([[ 301.,  311.,  302.,  312.,  303.],
        [ 313.,  304.,  314.,  305.,  315.]], dtype=float32))

What doesn't currently work is to have arrays of mixed types in the same 
write statement, e.g.

integer :: index(10)
real :: x(10,10)
...
write(13) x, index

To address the original problem, I've changed the code to default to the 
native byte-ordering (f.ENDIAN='@') and to be more informative about 
what happened in the error.  In the latest version (attached):

In [1]: from fortranfile import FortranFile

In [2]: f = FortranFile('uxuyp.bin', endian='>') # incorrect endian-ness

In [3]: u = f.readReals()

IOError: Could not read enough data.  Wanted 1342177280 bytes, got 132

and hopefully when people see crazy big numbers like 1.34e9 they will 
think of byte order problems.

> I'm going to dive into struc, FotranFile etc.. to propose something
> convenient for people who have to read unformatted binary fortran file
> very often.

Awesome!  The thoughts banging around in my head right now are that some 
sort of mini-language that encapsulates the content of the declarations 
and the write statements should allow one to tease out exactly which 
struct call will unpack the right information.  f2py has some fortran 
parsing capabilities, so you might be able to use the fortran itself as 
the mini-language.  Something like

spec = fortranfile.OutputSpecification(\
"""real(4),dimension(2,5):: ux,uy
write(11) ux,uy""")
ux, uy = fortranfile.FortranFile('uxuyp.bin').readSpec(spec)

Best of luck.  Peace,

-Neil
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: fortranfile.py
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20090529/f24fb616/attachment.ksh>