使用包含多个变量的控制文件读取Python中的二进制文件

Question

我试图以一个Fortran代码为例来读取Python中的二进制文件。在

这个二进制文件名为data.grads，我有一个名为data.ctl的控制文件，它应该允许我理解如何读取二进制文件。正如我所说，我有一个Fortran代码，我的导师写了这段代码来解释读取二进制文件和构造对应于控制文件中不同变量（速度、温度、压力等）的数组的过程。在

我读了很多关于这个问题的帖子，我在理解如何使用二进制文件中的数据时遇到了一些困难。在

以下是我查看的一些帖子：

• Reading a direct access fortran unformatted file in Python
• Reading direct access binary file format in Python
• Reading records from a Fortran binary file in Python
• Reading Fortran binary file in Python
• Reading fortran direct access data and writing formatted data - faster with python than with fortran?
• Python reading unformatted direct access Fortran 90 gives incorrect output

二进制文件存储了科学家用来绘制温度和压力剖面图的大气不同特性的模拟结果，例如：

我被告知控制文件是理解如何从二进制文件中获取任何内容的关键。在

我可以读取文件，但不知道如何访问特定变量的结果。以下是我从我之前引用的一些帖子中得到的一段代码：

filename = "/path/data.grads"

nlat = 32
nlon = 67
f = open(filename, 'rb')
recl = np.fromfile(f, dtype='int32', count=4*nlat*nlon)
f.seek(4)
field = np.fromfile(f, dtype='float32',count=-1)

print('Record length=',recl)
print(field, len(field))

它返回以下值：

^{pr2}$

有此方面经验的人能帮我弄清楚如何使用控制文件访问不同的变量吗？

如果你需要更多的解释，请告诉我，我会编辑我的帖子，并添加尽可能多的信息。在

---

不幸的是，我不能共享二进制文件，因为它在服务器上，重约40GB

我分享：

• 控制文件的内容
• 以及Fortran代码（它缺少一些代码，因为它是作为解释编写的）。在

控制文件(数据.ctl)

DSET ^data.grads
UNDEF -9.99e33
XDEF  64 LINEAR 0.0   5.6250
YDEF 32 LEVELS
 -85.761 -80.269 -74.745 -69.213 -63.679 -58.143 -52.607 -47.070 -41.532
 -35.995 -30.458 -24.920 -19.382 -13.844  -8.307  -2.769   2.769   8.307
  13.844  19.382  24.920  30.458  35.995  41.532  47.070  52.607  58.143
  63.679  69.213  74.745  80.269  85.761
ZDEF  67 LEVELS
.000696 .08558 .1705 .2554 .3402 .4544 .6274 .8599 1.152 1.505 1.918 2.392
2.928 3.495 4.063 4.781 5.816 7.181 8.889 10.95 13.39 16.07 18.81 21.61
24.47 27.39 30.37 33.40 36.50 39.64 42.77 45.90 49.04 52.17 55.31 58.44
61.57 64.71 67.84 70.98 74.11 77.24 80.38 83.51 86.65 89.78 92.92 96.05
99.18 102.3 105.5 108.6 111.7 114.9 118.0 121.1 124.3 127.4 130.5 133.7
136.8 139.9 143.1 146.2 149.3 152.5 160.0
TDEF 3120 LINEAR 01JAN2000 1HR
VARS 31
u    67 99  u (m/s)
v    67 99  v (m/s)
w    67 99  w (m/s)
T    67 99  T (K)
dia  67 99  diagnostics (see table)
ps    0 99  ps
Ts    0 99  Ts (K)
h2og   67 99   Water vapor [kg/m^2]
h2oim1  67 99   Water ice mass for dust 0.30E-07 [kg/m^2]
h2oim2  67 99   Water ice mass for dust 0.10E-06 [kg/m^2]
h2oim3  67 99   Water ice mass for dust 0.30E-06 [kg/m^2]
h2oim4  67 99   Water ice mass for dust 0.10E-05 [kg/m^2]
h2oin1  67 99   Water ice number for dust 0.30E-07 [number/m^2]
h2oin2  67 99   Water ice number for dust 0.10E-06 [number/m^2]
h2oin3  67 99   Water ice number for dust 0.30E-06 [number/m^2]
h2oin4  67 99   Water ice number for dust 0.10E-05 [number/m^2]
h2ois  0 99  surface h2o ice [kg/m^2]
p  67 99   Pressure [Pa]
h  67 99   Height above the surface [m]
dipre  67 99   Delta pressure [Pa]
surf   67 99   space of cell factor [sin*cos]
dm   67 99  CO2 ice cloud mass concentration [kg/m^3]
cap   0 99  Surface CO2 ice mass [kg]
hcap  0 99  Surface CO2 ice [kg/m^2]
gdq_op 0 99  dust from rad.mod [opacity]
gdq_mix  67 99  dust from rad.mod [mix.r]
dust_op 0 99  dust from tracers [opacity]
dust_n1  67 99   Dust num dens from tracers for R=0.30E-07 [m^-3]
dust_n2  67 99   Dust num dens from tracers for R=0.10E-06 [m^-3]
dust_n3  67 99   Dust num dens from tracers for R=0.30E-06 [m^-3]
dust_n4  67 99   Dust num dens from tracers for R=0.10E-05 [m^-3]
ENDVARS

Fortran文件

  open(12,file='data.grads',status='unknown',
 &           form='unformatted',access='direct',
 &           recl = 4*nlat*nlon )


  krec=0

  do l = 1,nlat
     krec = krec+1
     write(12,rec=krec,ERR=900)real(u3d(1:nlon,1:nlat,l))
  enddo

  do l = 1,nlat
     krec = krec+1
     write(12,rec=krec,ERR=900)real(v3d(:,:,l))
  enddo

  do l = 1,nlat
     krec = krec+1
     write(12,rec=krec,ERR=900)real(w3d(:,:,l))
  enddo

  do l = 1,nlat
     krec = krec+1
     write(12,rec=krec,ERR=900)real(T3d(:,:,l))
  enddo

  do l = 1,nlat
     krec = krec+1
     write(12,rec=krec,ERR=900)real(D3(:,:,l))
  enddo

  krec = krec+1
  write(12,rec=krec,ERR=900)real(ps3d(:,:))

  krec = krec+1
  write(12,rec=krec,ERR=900)real(Ts3d(:,:))

  do n = 1,4
   do l = 1,nlat
    krec = krec+1
    write(12,rec=krec,ERR=900)real(trace4D(:,:,l,n))
   enddo
  enddo

  do n = 1,4
    krec = krec+1
    write(12,rec=krec,ERR=900)real(trace2D(:,:,n))
  enddo

  do l = nlat,1,-1
     krec = krec+1
     write(12,rec=krec,ERR=900)real(pre3d(:,:,l))
  enddo

  do l = 1,nlat
     krec = krec+1
     write(12,rec=krec,ERR=900)real(alth3d(:,:,l))
  enddo

  do l = 1,nlat
     krec = krec+1
     write(12,rec=krec,ERR=900)real(dipre3d(:,:,l))
  enddo


  do l = 1,nlat
     krec = krec+1
     write(12,rec=krec,ERR=900)real(surf3d(:,:,l))
  enddo

  do l = 1,nlat
     krec = krec+1
! Convert DM to DM/grid volume
     write(12,rec=krec,ERR=900) real(DM4(:,latmask(:),l)
 &         *GRAV*pre3d(:,:,nlat-l+1)
 &         /(T3d(:,:,nlat-l+1)*dipre3d(:,:,nlat-l+1)*RGAS)
 &         /surf3d(:,:,nlat-l+1) )
  enddo

  krec = krec+1
  write(12,rec=krec,ERR=900)real(CAP4(:,latmask(:)))

  krec = krec+1
  write(12,rec=krec,ERR=900)real(HCAP4(:,latmask(:)))

  krec = krec+1
  write(12,rec=krec,ERR=900)real(DDUSTA3d(:,:))

  do l = 1,nlat
     krec = krec+1
     write(12,rec=krec,ERR=900)real(GDQ3d(:,:,l))
  enddo

  krec = krec+1
  write(12,rec=krec,ERR=900)real(tau_dust3d(:,:))

  do n = 1,4
   do l = 1,nlat
    krec = krec+1
    write(12,rec=krec,ERR=900)real(dust_n3d(:,:,l,n))
   enddo
  enddo

---

编辑：二进制文件的第一行（xxd -l 100 data.grads）

0000000: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7  .E...E...E...E..
0000010: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7  .E...E...E...E..
0000020: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7  .E...E...E...E..
0000030: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7  .E...E...E...E..
0000040: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7  .E...E...E...E..
0000050: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7  .E...E...E...E..
0000060: c345 f6f7                                .E..

Answer 1

该Fortran的numpy等价物（或者更确切地说是对应的，因为您向编写器展示了）

import numpy
nlat=32  # from data.ctl
nlon=64
nz=64    # or 67?
dt=numpy.dtype((numpy.float32,(nlat,nlon)))
def read(n=nz): return numpy.fromfile(f,dt,n)
def read1(): return read(1)[0]
with open("data.grads",'rb') as f:
  u3d=read()
  v3d=read()
  # …
  ps3d=read1()
  # …
  trace4D=numpy.fromfile(f,numpy.dtype((dt,nlat)),4)
  trace2D=read(4)
  pre3d=read()[::-1]
  # …

在nlat和nz变量之间存在一些混淆；因为您说过代码已更改，data.ctl可能是更好的参考。在