Reading Fortran Data

[Carl Banks]

Tyler wrote:
> Hello All:
>
> After trying to find an open source alternative to Matlab (or IDL), I
> am currently getting acquainted with Python and, in particular SciPy,
> NumPy, and Matplotlib. While I await the delivery of Travis Oliphant's
> NumPy manual, I have a quick question (hopefully) regarding how to read
> in Fortran written data.

I think you made a good choice, if for no other reason than the fact
that Python and numpy absolutely rock when it comes to interfacing with
Fortran (and any C code that can be called by Fortran).  I suggest
having a look at pyfort and/or f2py to see if they can be useful to
you.  If your data files are temporary (that is, they only exist to
pass data from one program to another), you might not even need them.

For completeness, I'll answer how to write as well as read Fortran
data.  It turns out that all of this can be done without numpy, but
there is one very convenient numpy function.

> The data files are not binary, but ASCII text files with no formatting
> and mixed data types (strings, integers, floats). For example, I have
> the following write statements in my Fortran code:
>
> I write the files as such:
>   WRITE(90,'(A30)') fgeo_name
>   WRITE(90,'(A30)') fmed_name

Fortran pads it's output when using a width field.  A good way to do
this in Python is to use string formating with a given width.  If
fgeo_name is a Python string, they you'd write something like this:

f.write("%-30s\n" % fgeo_name)

The negeative sign is to pad the name on the left, which is how Fortran
does it (by default, Python string formatting pads on the right).  You
could also use the .ljust method of string objects:

f.write(fgeo_name.ljust(30) + "\n")

Don't forget the newline on the end.

>   WRITE(90,*) nfault,npoint

Fortran writes this as two arbitrary integers separated by a space.  To
do this in Python, use string formating.  This also adds a space to the
beginning as most Fortran implementations seem to do.

f.write(" %d %d\n" % (nfault,npoint))

>   WRITE(90,*) (xpt(n), n=1,npoint)
>   WRITE(90,*) (ypt(n), n=1,npoint)

Now, an array.  I think you'd be safest writing one value per line.
You would do that like this:

for x in xpt:
    f.write(" %#g\n" % x)

The # sign is there to force the number to be formatted with a decimal
point.  You'll probably want to tweak the format string in other ways
(to specify a precision, for instance).


> and,
>
>   WRITE(10,'(A30)') fname
>   DO i=1,nfault
>      WRITE(10,*) dbn(i),dtn(i),xfwnt(i),yfwnt(i),xfent(i),yfent(i),&
>           &        slpvlS(i),slpvlD(i),slpvlT(i),segdp1(i)
>   END DO

I'd write this as one number per line.

f.write("%-30s\n" % fname)
for i in range(nfault):
    f.write(" %#g\n" % dbn[i])
    f.write(" %#g\n" % dtn[i])
    # and so on

If you know Python well, there are "more Pythonic" ways to do this, but
this is straightforward and works well enough.  Now for the reading
part.

> I then respectively read them into Fortran as:
>   READ(70,'(A30)') fgeo_name
>   READ(70,'(A30)') fmed_name

Once you've opened a file, then:

fgeo_name = f.readline().strip()
fmed_name = f.readline().strip()

Note that this strips the padding off the name.  If the first line of
the file is "ABC" followed by 27 spaces, the result will be "ABC" with
the spaces stripped.

>   READ(70,*) nfault,npoint

Python doesn't have any built in input formating (a la READ in Fortran
or scanf in C), so one usually does this kind of thing by hand.
Fortunately, Python makes this quite easy.  The following will do what
you want:

s = f.readline().split()
nfault = int(s[0])
npoint = int(s[1])

Here's what happens: it reads in a line, and splits the line on
whitespace into substrings.  It assigns the list of substrings to s.
Then it converts the first substring (s[0]) to an int and assigns it to
nfault; the second (s[1]) to npoint.

With more advanced knowledge of Python, you could write it in one line
like this:

nfault,npoint = (int(ss) for ss in f.readline().split())


>   READ(70,*) (x(n), n=1,npoint)
>   READ(70,*) (y(n), n=1,npoint)

Fortran programs seem to wrap free-form output to 80 columns, inserting
newlines whenever.  Because of this, you can't really free-form read
data in a line-by-line way (unless it's short, like the above example),
but now have to read data number-by-number, which is not so
straightforward.

Fortunately, numpy has a function, fromfile, which reads data
number-by-number.  If f is the file object you're reading from, then
you could read in the array like this:

x = numpy.fromfile(f,sep=" ",count=npoint,dtype=numpy.Float)

According to my tests, this leaves f pointing right past the end of the
array, so that you can still read in other arrays and variables.  Don't
forget to specify the type.


> and,
>
>   READ(20,'(A30)') fname
>   DO i=1,nfault
>      READ(20,*) dbn(i),dtn(i),xfwnt(i),yfwnt(i),xfent(i),yfent(i),&
>           &        slpvlS(i),slpvlD(i),slpvlT(i),segdp1(i)
>   END DO

Because this line has a very good chance of being split up upon output,
and because they're all real numbers, I'd use numpy.fromfile to read
this, and then assign to the different arrays.

So,

for i in range(nfault):
    tmp = numpy.fromfile(f,sep=" ",count=10,dtype=numpy.Float)
    dbn[i] = tmp[0]
    dtn[i] = tmp[1]
    # and so on

More advanced stuff, such as use of slices, can simplify this a lot.

> I also read them into IDL for visualization using the "READF" command.
> I was wondering how I might go about reading this into Python using
> NumPy. If this is not trivial, let me know and I'll just wait until the
> NumPy manual arrives.

One thing I highly recommend is that you not ignore the Pure Python
stuff.  Many of the solutions I gave you here didn't use numpy at all.
Successful usage of the numpy is going to depend a lot on knowing
normal Python.  For numerical applications, I suggest you become very
familiar with file I/O, strings and string handling, basic container
types (lists, tuples, dicts, and sets) and the operations they support
(slicing, iteration), and of course, numerics.


Carl Banks