6

Background:

A binary file is read on a Linux machine using the following Fortran code:

        parameter(nx=720, ny=360, nday=365)
c 
        dimension tmax(nx,ny,nday),nmax(nx,ny,nday)
        dimension tmin(nx,ny,nday),nmin(nx,ny,nday)
c 
        open(10,
     &file='FILE',
     &access='direct',recl=nx*ny*4)
c
        do k=1,nday
        read(10,rec=(k-1)*4+1)((tmax(i,j,k),i=1,nx),j=1,ny) 
        read(10,rec=(k-1)*4+2)((nmax(i,j,k),i=1,nx),j=1,ny) 
        read(10,rec=(k-1)*4+3)((tmin(i,j,k),i=1,nx),j=1,ny) 
        read(10,rec=(k-1)*4+4)((nmin(i,j,k),i=1,nx),j=1,ny) 
        end do

File Details:

options  little_endian
title global daily analysis (grid box mean, the grid shown is the center of the grid box)
undef -999.0
xdef 720 linear    0.25 0.50
ydef 360  linear -89.75 0.50
zdef 1 linear 1 1
tdef 365 linear 01jan2015 1dy
vars 4
tmax     1  00 daily maximum temperature (C)
nmax     1  00 number of reports for maximum temperature (C)
tmin     1  00 daily minimum temperature (C)
nmin     1  00 number of reports for minimum temperature (C)
ENDVARS

Attempts at a solution:

I am trying to parse this into an array in python using the following code (purposely leaving out two attributes):

with gzip.open("/FILE.gz", "rb") as infile:
     data = numpy.frombuffer(infile.read(), dtype=numpy.dtype('<f4'), count = -1)

while x <= len(data) / 4:
    tmax.append(data[(x-1)*4])
    tmin.append(data[(x-1)*4 + 2])
    x += 1

data_full = zip(tmax, tmin)

When testing some records, the data does not seem to line up with some sample records from the file when using Fortran. I have also tried dtype=numpy.float32 as well with no success. It seems as though I am reading the file in correctly in terms of number of observations though. I was also using struct before I learned the file was created with Fortran. That was not working

There are similar questions out here, some of which have answers that I have tried adapting with no luck.

UPDATE

I am a little bit closer after trying out this code:

#Define numpy variables and empty arrays
nx = 720 #number of lon
ny = 360 #number of lat
nday = 0 #iterate up to 364 (or 365 for leap year)   
tmax = numpy.empty([0], dtype='<f', order='F')
tmin = numpy.empty([0], dtype='<f', order='F')

#Parse the data into numpy arrays, shifting records as the date increments
while nday < 365:
    tmax = numpy.append(tmax, data[(nx*ny)*nday:(nx*ny)*(nday + 1)].reshape((nx,ny), order='F'))
    tmin = numpy.append(tmin, data[(nx*ny)*(nday + 2):(nx*ny)*(nday + 3)].reshape((nx,ny), order='F'))
    nday += 1  

I get the correct data for the first day, but for the second day I get all zeros, the third day the max is lower than the min, and so on.

5
  • Is direct access required in the Fortran code? (BTW, the name of the language is Fortran). Use access='stream'. It will make your life easier. Otherwise, how are you dealing with Fortran recorder markers that are written into the file?
    – evets
    Commented Sep 26, 2018 at 15:09
  • The fortran code was provided by a third party. I am assuming it is necessary but I do not have the prerequisite knowledge of the language to know if it is absolutely necessary. I also do not have access to that Linux box since, as mentioned, its a third party. Also, recorder markers, I will have to do some research on that. I am not familiar. Thanks.
    – Mark P.
    Commented Sep 26, 2018 at 15:14
  • The Fortran fragment of the question doesn't create the file, performing input. Generally, however, such "unformatted" files will have a internal structure which is not portable beyond compilers (and possibly even compiler versions). You will need to understand how the file was created. [The previous comment had a typo: it's "record markers", rather than "record_er_".] Commented Sep 26, 2018 at 15:26
  • Is there no easy way to reverse engineer the above fortran code in order to read the file in python? Without understand the exact process to create the file? I edited my question based on your feedback. Thanks.
    – Mark P.
    Commented Sep 26, 2018 at 15:38
  • @francescalus, sorry about 'recorder' instead of 'record' (unconscious holdover from work). Mark P, the start of each Fortran record will have a 4 or 8 byte marker, and the same marker might be written at the end of record. You'll need to skip those markers. access='stream' does not write markers. It is as the name suggests a stream of bytes.
    – Steve
    Commented Sep 26, 2018 at 19:01

2 Answers 2

1

While the exact format of Fortran binary files is compiler dependent, in all cases I'm aware of direct access files (files opened with access='direct' as in this question) do not have any record markers between records. Each record is of a fixed size, as given by the recl= specifier in the OPEN statement. That is, the record N starts at offset (N - 1) * RECL bytes in the file.

One portability gotcha is that the unit of the recl= is in terms of file storage units. For most compilers, the file storage unit specifies the size in 8-bit octets (as recommended in recent versions of the Fortran standard), but for the Intel Fortran compiler, recl= is in units of 32 bits; there is a commandline option -assume byterecl which can be used to make Intel Fortran match most other compilers.

So in the example given here and assuming a 8-bit file storage unit, your recl would be 1036800 bytes.

Further, looking at the code, it seems to assume the arrays are of a 4-byte type (e.g. integer or single precision real). So if it's single precision real, and the file has been created in little endian, then the numpy dtype <f4 that you have used seems to be the correct choice.

Now, getting back to the Intel Fortran compiler gotcha, if the file has been created by ifort without -assume byterecl then the data you want will be in the first quarter of each record, with the rest being padding (all zeros or maybe even random data?). Then you'll have to do some extra gymnastics to extract the correct data in python and not the padding. It should be easy to check this by checking the size of the file, is it nx * ny * 4 * nday *4 or is it nx * ny * 4 * nday * 4 * 4 bytes?

11
  • Thanks. So the correct file size is nx * ny * 4 * nday meaning 360*720 records per day, with 4 observations per day, with 365 records per year. Having said this, there is clearly extra padding because the file sitting on the FTP server is 1.4 GB in size.
    – Mark P.
    Commented Sep 27, 2018 at 11:08
  • Added ctl details from the server.
    – Mark P.
    Commented Sep 27, 2018 at 11:11
  • @MarkP.: Ah yes, my estimate is wrong by a factor of 4. One 4 is due to the size of each element, but then there's another due to the 4 different values per day. I'll fix that.
    – janneb
    Commented Sep 27, 2018 at 11:38
  • With your additions, I would say the file is closest in size to nx * ny * 4 * nday * 4. The total amount of data per record that I want is 16 bytes (4 btyes * 4 vars). I need 720 * 360 records per day and 365 records per year.
    – Mark P.
    Commented Sep 27, 2018 at 11:47
  • Well, if you're talking about Fortran records, then each record contains nx*ny elements of 4 bytes each. Then for each day you have 4 records, for the 4 different variables.
    – janneb
    Commented Sep 27, 2018 at 12:12
0

After the Update in my question, I realize I had an error with how I was looping. I of course spotted this about 10 minutes after issuing a bounty, aw well.

The error is with using the day to iterate through the records. This will not work as it iterates once per loop, not pushing the records far enough. Hence why some mins were higher than maxes. The new code is:

while nday < 365:
    tmax = numpy.append(tmax, data[(nx*ny)*rm:(nx*ny)*(rm + 1)].reshape((nx,ny), order='F'))
    rm = rm + 2
    tmin = numpy.append(tmin, data[(nx*ny)*rm:(nx*ny)*(rm + 1)].reshape((nx,ny), order='F'))
    rm = rm + 2
    nday += 1 

This used a Record Mover (or rm as I call it) to move the records the appropriate amount. That was all it needed.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.