Reading Fortran Data

[Tyler]

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.

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
WRITE(90,*) nfault,npoint
WRITE(90,*) (xpt(n), n=1,npoint)
WRITE(90,*) (ypt(n), n=1,npoint)

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 then respectively read them into Fortran as:
READ(70,'(A30)') fgeo_name
READ(70,'(A30)') fmed_name
READ(70,*) nfault,npoint
READ(70,*) (x(n), n=1,npoint)
READ(70,*) (y(n), n=1,npoint)

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

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.

Cheers,

t.

 

[Eric_Dexter]

I don't know if this is helpfull or not but (or for that matter
current). http://cens.ioc.ee/projects/f2py2e/ offers some suggestions
and it looks like you can use it with c code also.

 

[Robert Kern]

I don't know if this is helpfull or not but (or for that matter
current). http://cens.ioc.ee/projects/f2py2e/ offers some suggestions
and it looks like you can use it with c code also.

f2py has been folded into numpy.

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco

 

[Robert Kern]

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.

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:

Konrad Hinsen has a module for reading this kind of file.

http://dirac.cnrs-orleans.fr/ScientificPython/

Specifically, Scientific.IO.FortranFormat .

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco

 

[Beliavsky]

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.

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:

In plain Python, you can read each line in to a string, break the
string into "words" using split, and then convert the words into
variables of the desired types. If you are new to Python, this is an
important idiom to learn. I don't know if NumPy has facilities to do
this more easily.

I write the files as such:
WRITE(90,'(A30)') fgeo_name
WRITE(90,'(A30)') fmed_name

Let me comment on the Fortran code. For the following lines using
list-directed output, the compiler has considerable freedom in how it
writes the output. I guess you expect the integers nfault and npoint to
be written on one line and the vectors xpt and ypt to each be written
on separate lines. The compiler could print each number on a separate
line and be standard-conforming. This does not matter if you are going
to use a Fortran list-directed read to read the file, but it will
matter if you are using other languages. I suggest that you use format
strings to get more control over the ouptput format before you think
about reading the output files in Python. Otherwise you will be trying
to hit a moving target.

<snip>

 

[Carl Banks]

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)
f.write(" %#g\n" % dtn)
# 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 = tmp[0]
dtn = 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

 

[Beliavsky]

Carl Banks wrote:

Fortran writes this as two arbitrary integers separated by a space.

I wrote a paragraph in my reply explaining why this is wrong. A Fortran
list-directed write can print results in an almost arbitrary format,
depending on the compiler. Many compilers will separate integers by
several spaces, not just one, and they could use commas instead of
spaces if they wanted. The number of items printed before a new line is
started is also compiler-dependent. For more control, one uses a
formatted write, for example

write (90,"(2(1x,i0))") nfault,npoint

<snip>

 

[Carl Banks]

Carl Banks wrote:



I wrote a paragraph in my reply explaining why this is wrong.

It's a safe assumption for a line of two integers. It might not
exactly produce what a Fortran program would, but it would work in a
read statement.

A Fortran
list-directed write can print results in an almost arbitrary format,
depending on the compiler. Many compilers will separate integers by
several spaces, not just one, and they could use commas instead of
spaces if they wanted.

1. Hardly any compiler will produce a line of two integers, or reals,
that another compiler couldn't read back.
2. The number of spaces separating the numbers isn't important when
reading back free-form data.
3. Fear that a Fortran compiler might use commas or wrap lines at ten
columns or whatever, because it's not based on a standard, is misguided
paranoia.

The number of items printed before a new line is
started is also compiler-dependent. For more control, one uses a
formatted write, for example

write (90,"(2(1x,i0))") nfault,npoint

I think it's just more work to guard against something isn't very
relevant in practice.


Carl Banks

 

[Beliavsky]

Carl Banks wrote:

1. Hardly any compiler will produce a line of two integers, or reals,
that another compiler couldn't read back.

Yes, but for more than three numbers, the statement is wrong. Intel
Fortran prints four double precision random n
as

0.555891433847495 0.591161642339424 0.888434673900224

0.487293557925127

but g95 prints them on a single line. I advise against using
list-directed Fortran writes to create files that other programs will
read, and I think most experienced Fortran programmers would agree.

 

[Carl Banks]

Carl Banks wrote:



Yes, but for more than three numbers, the statement is wrong. Intel
Fortran prints four double precision random n
as

0.555891433847495 0.591161642339424 0.888434673900224

0.487293557925127

but g95 prints them on a single line.

Did you try to use one compiler's program's output as the other's
input? How did it work?

I advise against using
list-directed Fortran writes to create files that other programs will
read, and I think most experienced Fortran programmers would agree.

It's been awhile, but I'd consider myself an experienced Fortran
programner. I don't really agree.


Carl Banks

 

[Tyler]

Wow! Thanks for the help everyone.

I will look into each of your comments in more detail in the morning,
but I'll mention a few things I guess. The first is, the list-directed
output was not necesarily my choice as some of the code is also used by
my supervisor (and her colleagues) and I had to keep the IO similar.
However, for "non-legacy" codes (even though is was upgraded to
Fortran90), I can ensure you that I will employ a more portable format.

Also, I am very interested in the Scientific.IO module mentioned and
will look into it further along with f2py.

Once again, I thank all of you for your replies, I think I even learnt
a few things about Fortran here too....

Cheers,

t.