如何將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. 

這是我需要一些assitance：

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

首先，我必須假設你寫的使用順序訪問而不是直接訪問，這是一個非常重要的區別原來的Fortran文件。我知道IDL，並且您在示例中設置的讀取命令與順序訪問一致。

這很重要，因爲對於每次調用WRITE（），Fortran順序訪問都會爲文件添加「記錄標記」，因此您必須在Python讀取例程中考慮這一點，就像在喬治鏈接的答案中一樣評論。

至於行/列順序，真正重要的是哪個維度迭代最快的內存。例如，Fortran二維數組將在內存中寫入第一個數組維度作爲最快的迭代。當你將同樣的東西讀入2D Python變量時，它將成爲最後一個維度 - 在重新定型平面數組時，記住這一點，如果要保留索引，則必須對其進行轉置。

我認爲下面的代碼會做約你想要什麼，爲了簡便起見我只包括你的恩和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進行索引。這可能不是最簡單或最好的方式，但我認爲它至少會讓你開始。