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

Question

我想用 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 程序读取的数据完全错误。尽管如此，我编写的字节流长度还可以，因为从以相同方式编写的下一个数据集中读取的元数据就可以了。

执行此操作的正确方法是什么？我们将非常感谢您的建议!

Answer 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