如何将fortran 77未格式化的二进制文件读入python

Question

我在 IDL 中有一个脚本，它读取未格式化的二进制文件 (F77) 并将其输出为 .sav 文件，但我想将此脚本转换为 python 并另存为 .npz 文件，但我在读取时遇到问题行。

IDL 代码：

;create save file for QBO model output

;---------------------------------------------------------------------
;input data here (Only adjust this info)
;---------------------------------------------------------------------

  time=7300  ;duration from fortran code
  tstep=5   ;time step from fortran code
  inputfile='/home/cwilliams/Metr_205/qbo/uvwtom.'; binary file from fortran
  outputsavefile='~cwilliams/Metr_205/qbo/qbo_FULLX.sav'; name of save file with     variables
;--------------------------------------------------------------------------------
;--------------------------------------------------------------------------------





;==============================================================================
;   Do Not Adjust the Following Code
;==============================================================================
fname='dat'
nstep=time/tstep ;time/timestep
nm=nstep-1

um=fltarr(17,34)
vm=fltarr(17,34)
wm=fltarr(17,34)
tm=fltarr(17,34)
pm=fltarr(17,34)
om=fltarr(17,34)
vpm=fltarr(17,34)
vdm=fltarr(17,34)


u=fltarr(17,34,nstep)
v=fltarr(17,34,nstep)
w=fltarr(17,34,nstep)
temp=fltarr(17,34,nstep)
pressure=fltarr(17,34,nstep)
ozone=fltarr(17,34,nstep)

date=findgen(nstep)*(2.*tstep/365.)
ulevels2=[-30, -25, -20, -15, -10, -5, 5, 10, 15,  20, 25, 30]
ulevels3=findgen(15)*3.-12
lab=findgen(20)*0+1
lat=findgen(17)*2.7-1.35

ht=findgen(34)*0.5+17.


openr,37,inputfile+fname,/f77_unformatted & title='base'
for n=0,nm do begin
    readu,37,um,vm,wm,tm,om,pm
    u(*,*,n)=um/100.
    v(*,*,n)=vm
    w(*,*,n)=wm
    temp(*,*,n)=tm*700000./2.87e6
    pressure(*,*,n)=pm/10000.
    ozone(*,*,n)=om
endfor
close,37

save,filename=outputsavefile,u,v,w,temp,pressure,ozone,date,ulevels2,ulevels3,lab,$
lat,ht
stop
end

Python 代码： #为QBO模型输出创建保存文件

#---------------------------------------------------------------------
#input data here (Only adjust this info)
#---------------------------------------------------------------------

time=7300  #duration from fortran code
tstep=5   #time step from fortran code
inputfile='/home/cwilliams/metr51/lab12/uvwtom.dat'# binary file from fortran
outputsavefile='~cwilliams/metr51/lab12/qbo_FULLX.sav'# name of save file with variables
#--------------------------------------------------------------------------------
#--------------------------------------------------------------------------------


import numpy as n


#==============================================================================
#   Do Not Adjust the Following Code
#==============================================================================
nstep=time/tstep #time/timestep
nm=nstep-1

um=n.empty((17,34))
vm=n.empty((17,34))
wm=n.empty((17,34))
tm=n.empty((17,34))
pm=n.empty((17,34))
om=n.empty((17,34))
vpm=n.empty((17,34))
vdm=n.empty((17,34))


u=n.empty((17,34,nstep))
v=n.empty((17,34,nstep))
w=n.empty((17,34,nstep))
temp=n.empty((17,34,nstep))
pressure=n.empty((17,34,nstep))
ozone=n.empty((17,34,nstep))

date=n.arange(nstep)*(2.*tstep/365.)
ulevels2=n.array([-30, -25, -20, -15, -10, -5, 5, 10, 15,  20, 25, 30])
ulevels3=n.arange(15)*3.-12
lab=n.arange(20)*0+1
lat=n.arange(17)*2.7-1.35

ht=n.arange(34)*0.5+17.

这是我需要帮助的地方：

f=open(inputfile,'rb')     
data=f.read()
for i in range(nm+1):
#    readu,37,um,vm,wm,tm,om,pm
#    u(:,:,i)=um/100.
#    v(:,:,i)=vm
#    w(:,:,i)=wm
#    temp(:,:,i)=tm*700000./2.87e6
#    pressure(:,:,i)=pm/10000.
#    ozone(*,*,n)=om
#
#n.savez(filename=outputsavefile,u=u,v=v,w=w,temp=temp,pressure=pressure,ozone=ozone,date=date,ulevels2=ulevels2,ulevels3=ulevels3,lab=lab,lat=lat,ht=ht)

我知道 IDL 和 Python 之间的行/列顺序存在问题，但我想我可以在读入代码后解决这个问题。

Answer 1

首先，我必须假设您使用 sequential 访问而不是 直接 访问来编写原始 Fortran 文件，这是一个非常重要的区别。我确实知道 IDL，并且您在示例中设置的读取命令与顺序访问一致。

这很重要，因为 Fortran 顺序访问会在每次调用 WRITE() 时将“记录标记”添加到文件中，因此您必须在 Python 读取例程中考虑这一点，就像乔治在厘米。

至于行/列顺序，重要的是哪个维度在内存中迭代速度最快。例如，Fortran 二维数组将写入内存中，第一个数组维度是最快的迭代。当您将相同的内容读入 2D Python 变量时，它将是最后一个维度 - 在重塑平面数组时请记住这一点，如果您想保留索引，则必须对其进行转置。

我认为以下代码将大致完成您想要的操作，为简洁起见，我只包含了您的 um 和 vm 变量：

import numpy as np

um=np.empty((34,17), dtype='float32') # Make these dimensions "backwards" for easier reshaping
vm=np.empty((34,17), dtype='float32') # Also watch out, I believe the default type is float64

f = open(inputfile,'rb')
recl = np.zeros(1,dtype=np.uint32)
for i in range(nm+1):
    recl = np.fromfile(f, dtype='uint32', count=1)
    tmpu = np.fromfile(f, dtype='float32', count=um.size) # These arrays will be flat
    tmpv = np.fromfile(f, dtype='float32', count=vm.size)
    recl = np.fromfile(f, dtype='uint32', count=1)

    um = np.transpose(np.reshape(tmpu, um.shape))
    vm = np.transpose(np.reshape(tmpv, vm.shape))

然后你也可以在循环中对你的其他变量做任何你想做的事情，你应该能够像在 IDL 或 Fortran 中习惯的那样对 um 和 vm 进行索引。这可能不是最简单或最好的方法，但我认为它至少可以帮助您入门。