使用 Python 将数组写入二进制文件，然后使用 Fortran 读取答案

【问题标题】：Write arrays to binary file with Python later read with Fortran使用 Python 将数组写入二进制文件，然后使用 Fortran 读取
【发布时间】：2019-05-23 18:37:42
【问题描述】：

我想用 Python 为天气研究和预报模型编写中间文件。这些是模型使用 Fortran 读取的二进制文件。文件包括一些元数据和一个二维字段。编写元数据工作正常，例如由

# record1
g.write(struct.pack('>i', 4)) # record opening bytes
g.write(struct.pack('>i', d.ifv))
g.write(struct.pack('>i', 4)) # record closing bytes`

但是当我想编写稍后被某些 Fortran 代码读取的字段时，例如

real, dimension(:,:):: slab
...
read (1) slab

我挣扎。我想我需要将 numpy 数组转换为写入 np.float32 类型，并且我必须以 Fortran 顺序编写它。有一个页面，其中显示了一些读取中间文件的代码，并且读取数据字段是由

def _parse_record5(data, nx, ny):
  result = {}
  size = nx * ny
  fmt = ">{}f".format(size)
  parsed = struct.unpack(fmt, data)
  arr = np.array(parsed, dtype=np.float32)
  result["slab"] = arr.reshape((nx, ny), order="F")
  return result

但我无法弄清楚如何正确编写该字段。我试着这样做：

recl = d.nx*d.ny*4  # number of array elements + 4 bytes
g.write(struct.pack('>i', recl)) # record opening bytes
slab_temp = d.slab.astype(np.float32) # converting from 64 to 32 bytes
slab = slab_temp.tobytes('F')         # bytestream in Fortran order
g.write(struct.pack(str(recl)+'s', slab))  # writing bytestream to file
g.write(struct.pack('>i', recl)) # record closing bytes`

很遗憾，它无法正常工作，Fortran 程序读取的数据完全错误。尽管如此，我写的字节流的长度还可以，因为从下一个数据集中读取的元数据以相同的方式写入就可以了。

这样做的正确方法是什么？您的建议将不胜感激！

【问题讨论】：

未格式化的 Fortran 二进制文件所需的数据结构未标准化。请告诉我们使用的 Fortran 编译器和写入数据的大致大小。此外，如果您可以给出您认为应该阅读的值以及您所看到的内容，这将很有用。
您可以为中间文件使用标准文件格式吗？可能会有一点开销，但由于格式是标准化的，所以 Python 和 Fortran 都应该有很好的库。鉴于该主题，netCDF4（或其更通用的基础 HDF5）可能是不错的选择。
使用流而不是基于记录的 io
在我看来，改变 Fortran 端的建议并没有那么有用。即使我自己是一名 Fortran 程序员，我当然也不想仅仅因为这个原因而处理像 WRF 这样庞大的代码的 I/O 内部。上游无论如何都不会改变，一个人必须保持一个不同的本地分支并保持更新。
@9769953：正如 Vladimir F 评论的那样，修改 WRF 系统以读取另一种文件类型作为中间文件并不是一个好主意。中间文件工作得很好，只有当输入数据不是 grib 格式或有人想修改输入数据时才需要自己编写。

标签： python arrays fortran

【解决方案1】：

最后，我能够解决问题。在将数组写入文件时，我错过了指定字节序（WRF 需要大字节序，但我们的 linux 计算机以小字节序运行）。我已经用大端写过所有其他变量，因此它们被正确读取。

为了完整起见，我把读写WRF中间文件的python代码放在下面：

# reads and writes dataset from / to a WRF intermediate file
import struct
import numpy as np

class dataset:
  ifv = 0.0
  hdate = ''
  xfcst = 0.0
  map_source = ''
  field = ''
  units_string = ''
  desc = ''
  xlvl = 0.0
  nx = 0
  ny = 0
  iproj = 0
  startloc = ''
  startlat = 0.0
  startlon = 0.0
  deltalat = 0.0
  deltalon = 0.0
  earth_radius = 0.0
  dx = 0.0
  dy = 0.0
  truelat1 = 0.0
  truelat2 = 0.0
  xlonc = 0.0
  nlats = 0.0
  is_wind_earth_rel = None
  slab = None

def print_dataset(res):
  print('ifv: ', res.ifv)
  print('hdate: ', res.hdate)
  print('xfcst: ', res.xfcst)
  print('map_source: ', res.map_source)
  print('field: ', res.field)
  print(' units_string: ', res.units_string)
  print(' desc: ', res.desc)
  print(' xlvl: ', res.xlvl)
  print(' nx: ', res.nx)
  print(' ny: ', res.ny)
  print(' iproj: ', res.iproj)
  print(' startloc: ', res.startloc)
  print(' startlat: ', res.startlat)
  print(' startlon: ', res.startlon)
  print(' deltalat: ', res.deltalat)
  print(' deltalon: ', res.deltalon)
  print(' earth_radius: ', res.earth_radius)
  print(' dx: ', res.dx)
  print(' dy: ', res.dy)
  print(' truelat1: ', res.truelat1)
  print(' truelat2: ', res.truelat2)
  print(' xlonc: ', res.xlonc)
  print(' nlats: ', res.nlats)
  print(' is_wind_earth_rel: ', res.is_wind_earth_rel)
  print(' slab: ', res.slab)

def read_dataset(f):
  d = dataset()
  recl=struct.unpack('>i',f.read(4))[0]
  d.ifv = struct.unpack('>i',f.read(recl))[0]
  recl_close=struct.unpack('>i',f.read(4))[0] # record closing bytes
  recl=struct.unpack('>i',f.read(4))[0]
  record2= struct.unpack(str(recl)+'s',f.read(recl))[0]
  d.hdate = record2[0:24].decode('UTF8')
  d.xfcst = struct.unpack('>f', record2[24:28])[0]
  d.map_source = record2[28:60].decode('UTF8')
  d.field = record2[60:69].decode('UTF8')
  d.units_string = record2[69:94].decode('UTF8')
  d.desc = record2[94:140].decode('UTF8')
  d.xlvl = struct.unpack('>f', record2[140:144])[0]
  d.nx = struct.unpack('>i', record2[144:148])[0]
  d.ny = struct.unpack('>i', record2[148:152])[0]
  d.iproj = struct.unpack('>i', record2[152:156])[0]
  recl_close=struct.unpack('>i',f.read(4))[0] # record closing bytes
  recl=struct.unpack('>i',f.read(4))[0]
  record3= struct.unpack(str(recl)+'s',f.read(recl))[0]
  d.startloc = record3[0:8].decode('UTF8')
  if(d.iproj == 0): # Cylindrical equidistant projection
    d.startlat = struct.unpack('>f', record3[8:12])[0]
    d.startlon = struct.unpack('>f', record3[12:16])[0]
    d.deltalat = struct.unpack('>f', record3[16:20])[0]
    d.deltalon = struct.unpack('>f', record3[20:24])[0]
    d.earth_radius = struct.unpack('>f', record3[24:28])[0]
  if(d.iproj == 1): # Mercator projection
    d.startlat = struct.unpack('>f', record3[8:12])[0]
    d.startlon = struct.unpack('>f', record3[12:16])[0]
    d.dx = struct.unpack('>f', record3[16:20])[0]
    d.dy = struct.unpack('>f', record3[20:24])[0]
    d.truelat1 = struct.unpack('>f', record3[24:28])[0]
    d.earth_radius = struct.unpack('>f', record3[28:32])[0]
  if(d.iproj == 3): # Lambert conformal projection
    d.startlat = struct.unpack('>f', record3[8:12])[0]
    d.startlon = struct.unpack('>f', record3[12:16])[0]
    d.dx = struct.unpack('>f', record3[16:20])[0]
    d.dy = struct.unpack('>f', record3[20:24])[0]
    d.xlonc = struct.unpack('>f', record3[24:28])[0]
    d.truelat1 = struct.unpack('>f', record3[28:32])[0]
    d.truelat2 = struct.unpack('>f', record3[32:36])[0]
    d.earth_radius = struct.unpack('>f', record3[36:40])[0]
  if(d.iproj == 4): # Gaussian projection
    d.startlat = struct.unpack('>f', record3[8:12])[0]
    d.startlon = struct.unpack('>f', record3[12:16])[0]
    d.nlats = struct.unpack('>f', record3[16:20])[0]
    d.deltalon = struct.unpack('>f', record3[20:24])[0]
    d.earth_radius = struct.unpack('>f', record3[24:28])[0]
  if(d.iproj == 5): # Polar-stereographic projection
    d.startlat = struct.unpack('>f', record3[8:12])[0]
    d.startlon = struct.unpack('>f', record3[12:16])[0]
    d.dx = struct.unpack('>f', record3[16:20])[0]
    d.dy = struct.unpack('>f', record3[20:24])[0]
    d.xlonc = struct.unpack('>f', record3[24:28])[0]
    d.truelat1 = struct.unpack('>f', record3[28:32])[0]
    d.earth_radius = struct.unpack('>f', record3[32:36])[0]
  recl_close=struct.unpack('>i',f.read(4))[0] # record closing bytes
  recl=struct.unpack('>i',f.read(4))[0]
  record4= struct.unpack(str(recl)+'s',f.read(recl))[0]
  d.is_wind_earth_rel = struct.unpack('>i', record4[0:4])[0]
  recl_close=struct.unpack('>l',f.read(4))[0] # record closing bytes
  # record 5
  recl=struct.unpack('>i',f.read(4))[0]
  record5 = struct.unpack(str(recl)+'s', f.read(recl))[0]
  slab = np.frombuffer(record5, dtype='>f')
  d.slab = slab.reshape(d.nx,d.ny)
  recl=struct.unpack('>i',f.read(4))[0]
  return d

def write_dataset(g, d):
  # record1
  g.write(struct.pack('>i', 4)) # record opening bytes
  g.write(struct.pack('>i', d.ifv))
  g.write(struct.pack('>i', 4)) # record closing bytes
  # record2
  g.write(struct.pack('>i', 156)) # record opening bytes
  g.write(struct.pack('24s', bytes(d.hdate, 'utf-8')))
  g.write(struct.pack('>f', d.xfcst))
  g.write(struct.pack('32s', bytes(d.map_source, 'utf-8')))
  g.write(struct.pack('9s', bytes(d.field, 'utf-8')))
  g.write(struct.pack('25s', bytes(d.units_string, 'utf-8')))
  g.write(struct.pack('46s', bytes(d.desc, 'utf-8')))
  g.write(struct.pack('>f', d.xlvl))
  g.write(struct.pack('>i', d.nx))
  g.write(struct.pack('>i', d.ny))
  g.write(struct.pack('>i', d.iproj))
  g.write(struct.pack('>i', 156)) # record closing bytes
  # record3
  if(d.iproj == 0): # Cylindrical equidistant projection
    g.write(struct.pack('>i', 28)) # record opening bytes
    g.write(struct.pack('8s', bytes(d.startloc, 'utf-8')))
    g.write(struct.pack('>f', d.startlat))
    g.write(struct.pack('>f', d.startlon))
    g.write(struct.pack('>f', d.deltalat))
    g.write(struct.pack('>f', d.deltalon))
    g.write(struct.pack('>f', d.earth_radius))
    g.write(struct.pack('>i', 28)) # record closing bytes
  if(d.iproj == 1): # Mercator projection
    g.write(struct.pack('>i', 32)) # record opening bytes
    g.write(struct.pack('8s', bytes(d.startloc, 'utf-8')))
    g.write(struct.pack('>f', d.startlat))
    g.write(struct.pack('>f', d.startlon))
    g.write(struct.pack('>f', d.dx))
    g.write(struct.pack('>f', d.dy))
    g.write(struct.pack('>f', d.truelat1))
    g.write(struct.pack('>f', d.earth_radius))
    g.write(struct.pack('>i', 32)) # record closing bytes
  if(d.iproj == 3): # Lambert conformal projection
    g.write(struct.pack('>i', 40)) # record opening bytes
    g.write(struct.pack('8s', bytes(d.startloc, 'utf-8')))
    g.write(struct.pack('>f', d.startlat))
    g.write(struct.pack('>f', d.startlon))
    g.write(struct.pack('>f', d.dx))
    g.write(struct.pack('>f', d.dy))
    g.write(struct.pack('>f', d.xlonc))
    g.write(struct.pack('>f', d.truelat1))
    g.write(struct.pack('>f', d.truelat2))
    g.write(struct.pack('>f', d.earth_radius))
    g.write(struct.pack('>i', 40)) # record closing bytes
  if(d.iproj == 4): # Gaussian projection
    g.write(struct.pack('>i', 28)) # record opening bytes
    g.write(struct.pack('8s', bytes(d.startloc, 'utf-8')))
    g.write(struct.pack('>f', d.startlat))
    g.write(struct.pack('>f', d.startlon))
    g.write(struct.pack('>f', d.nlats))
    g.write(struct.pack('>f', d.deltalon))
    g.write(struct.pack('>f', d.earth_radius))
    g.write(struct.pack('>i', 28)) # record closing bytes
  if(d.iproj == 5): # Polar-stereographic projection
    g.write(struct.pack('>i', 36)) # record opening bytes
    g.write(struct.pack('8s', bytes(d.startloc, 'utf-8')))
    g.write(struct.pack('>f', d.startlat))
    g.write(struct.pack('>f', d.startlon))
    g.write(struct.pack('>f', d.dx))
    g.write(struct.pack('>f', d.dy))
    g.write(struct.pack('>f', d.xlonc))
    g.write(struct.pack('>f', d.truelat1))
    g.write(struct.pack('>f', d.earth_radius))
    g.write(struct.pack('>i', 36)) # record closing bytes
  # record4
  g.write(struct.pack('>i', 4)) # record opening bytes
  g.write(struct.pack('>i', d.is_wind_earth_rel))
  g.write(struct.pack('>i', 4)) # record closing bytes
  # record 5
  recl = d.nx*d.ny*4
  g.write(struct.pack('>i', recl)) # record opening bytes
  slab_temp = d.slab.astype('>f')
  g.write(slab_temp.tobytes(order='F'))
  g.write(struct.pack('>i', recl)) # record closing bytes

【讨论】：