python - 在 Python 中读取直接访问二进制文件格式

Question

背景:

使用以下 Fortran 代码在 Linux 机器上读取二进制文件:

        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

文件详细信息:

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

解决方案的尝试:

我正在尝试使用以下代码将其解析为 python 中的数组(故意省略两个属性):

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)

测试某些记录时，使用 Fortran 时，数据似乎与文件中的某些示例记录不一致。我也试过 dtype=numpy.float32 也没有成功。不过，就观察次数而言，我似乎正在正确读取文件。在我了解到该文件是使用 Fortran 创建的之前，我也在使用 struct。那是行不通的

这里有类似的问题，其中一些有我尝试改编但没有运气的答案。

更新

尝试这段代码后，我更接近了一点:

#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  

第一天我得到了正确的数据，但第二天我得到的全是零，第三天最大值低于最小值，依此类推。

Answer 1

虽然 Fortran 二进制文件的确切格式取决于编译器，但在所有情况下我都知道直接访问文件(在这个问题中用 access='direct' 打开的文件)在记录之间没有任何记录标记。每条记录都是固定大小的，由 recl= 给出。 OPEN 中的说明符陈述。即记录N从偏移量 (N - 1) * RECL 开始文件中的字节。

一个可移植性陷阱是 recl= 的单位是关于 file storage unit秒。对于大多数编译器，file storage unit以 8 位八位字节指定大小(如 Fortran 标准的最新版本所推荐)，但对于 Intel Fortran 编译器，recl=以32位为单位；有一个命令行选项 -assume byterecl它可用于使英特尔 Fortran 语言与大多数其他编译器相匹配。

因此，在此处给出的示例中，假设一个 8 位 file storage unit ，您的 recl 将为 1036800 字节。

此外，查看代码，似乎假设数组是 4 字节类型(例如整数或单精度实数)。因此，如果它是单精度实数，并且文件是在小端创建的，那么 numpy dtype <f4您使用的似乎是正确的选择。

现在，回到英特尔 Fortran 编译器陷阱，如果文件是由没有 -assume byterecl 的 ifort 创建的那么您想要的数据将位于每条记录的第一季度，其余部分将被填充(全部为零，甚至可能是随机数据？)。然后你必须做一些额外的体操来提取 python 中的正确数据而不是填充。通过检查文件的大小应该很容易检查这个，是吗nx * ny * 4 * nday *4或者是nx * ny * 4 * nday * 4 * 4字节？