使用具有多个变量的控制文件在 Python 中读取二进制文件答案

【问题标题】：Reading a binary file in Python using a control file with multiple variables使用具有多个变量的控制文件在 Python 中读取二进制文件
【发布时间】：2019-07-18 13:10:16
【问题描述】：

我正在尝试以 Fortran 代码为例，在 Python 中读取二进制文件。

二进制文件名为data.grads，我有一个名为data.ctl 的控制文件，它可以让我了解如何读取二进制文件。正如我所说，我有一个 Fortran 代码，我的顾问编写了该代码来解释读取二进制文件和构造与控制文件中不同变量（速度、温度、压力等）相对应的数组的过程。

我阅读了多篇关于此事的帖子，但在理解如何使用二进制文件中的数据时遇到了一些麻烦。

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

二进制文件存储了科学家用来绘制大气不同属性的模拟结果，例如，温度和压力部分，如下所示：

有人告诉我，控制文件是理解如何从二进制文件中获取任何内容的关键。

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

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))

它返回以下内容：

Record length= [-134855229 -134855229 -134855229 ... -134855229 -134855229 -134855229]
[-9.99e+33 -9.99e+33 -9.99e+33 ... -9.99e+33 -9.99e+33 -9.99e+33] 10319462399

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

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

很遗憾，我无法共享二进制文件，因为它位于服务器上，重约 40 GB……

我分享：

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

控制文件（data.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..

【问题讨论】：

您使用的是哪个编译器？指定 `recl = 4*nlat*nlon` 不可移植，一些编译器以字节为单位计算记录大小，其他以 4 字节字为单位。 stackoverflow.com/questions/37770912/…
在 Python 中，我将 IPython 与 Anaconda 的 Spyder 一起使用，但在 Fortran 中，我不知道他们使用哪种编译器......我刚刚获得了 Fortran 代码，但我精通 Fortran，所以我被要求用 Python 编写代码。这是我使用的 python 版本：Python 3.7.3。
编辑：“我不精通 Fortran”
另一个考虑，文件不能是大端的吗？
@VladimirF 我在终端xxd -b data.grads 中使用此命令打开了文件，我看到里面的每一行都写成：44s54c5: 10010101 11001010 01010101 00101010 10101010 11011010 =.6^E.（这不是文件中的一行，只是例如，文件太大了，我无法复制任何行，因为它仍在读取）。 unix.stackexchange.com/questions/282215/…

标签： python file-io fortran binaryfiles endianness

【解决方案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 可能是更好的参考。

【讨论】：

这太棒了，看来我用那个代码得到了很好的结果。你能告诉我如何根据 Fortran 代码读取以下变量 trace2D、DM4、CAP4、HCAP4 吗？此外，在控制文件中，您可以看到与时间相关的这一行TDEF 3120 LINEAR 01JAN2000 1HR。它表示模拟从 2000 年 1 月 1 日开始，并持续 1 小时（3120 小时）的 3120 个时间步长。所以要得到每小时的结果，我应该做with open("data.grads",'rb') as f: for t in range(3120): u3d=read() # … 吗？
@LoïcPoncin：我添加了一些，因为latmask 丢弃了信息，所以一般无法重建一些。如果该循环与文件大小一致，则该循环似乎是合理的，尽管这取决于您是要创建dict（或namedtuple）还是一次只处理一个数据而不保留所有数据。
所以我想如果我要使用温度字典temp = {} with open('data.grads','rb') as f: for t in range(3120): T3d=read() temp['t_{}'.format(t)] = T3d # …，我可以这样做
@LoïcPoncin：你可以，但为什么不直接使用 t 作为键——或者只使用列表（无论是 AoS 还是 SoA）？
SoA 似乎更适合我想做的事情，非常感谢您的帮助。干杯！

【解决方案2】：

我最近写了一个python包xgrads来解析和读取GrADS常用的ctl/binary文件。它可以使用numpy和dask加载整个数据集，并以xarray.Dataset的形式返回，这在地球科学中非常流行。

以下代码展示了如何加载 ctl 数据集并访问您感兴趣的不同变量：

from xgrads.core import open_CtlDataset

dset = open_CtlDataset('test.ctl')

# print all the info in ctl file
print(dset)

# for your ctl content, you can plot any variables (e.g.,
# first time and level of T) immediately as
dset['T'][0,0,...].plot()

虽然这是第一个版本，没有经过广泛测试，但我已经成功解析了您的 ctl 内容。我希望该软件包也可以返回正确的数据集。如果您发现任何错误，我也很乐意修复它。

【讨论】：