老夫聊發(fā)少年狂兰粉，試試Python氣象繪圖也殖。
現(xiàn)有數(shù)據(jù)為①記錄能譜密度的二進(jìn)制spect.dat既绕；②記錄均一化波數(shù)缴挖、均一化角波數(shù)和波長(zhǎng)的aol.txt（奇怪的是aol.txt自動(dòng)換成5個(gè)數(shù)值一行）袋狞。
import numpy as np
import struct
import matplotlib.pyplot as plt
def getbin(filename):
with open(filename,mode='rb') as f:
tmp1=f.read()
fsize=f.tell()
tmp11=list(struct.unpack('f'*26*479,tmp1))
#python沒(méi)有二進(jìn)制格式，所以需要用struck.unpack來(lái)把字符串轉(zhuǎn)譯成實(shí)數(shù)映屋，f指float/float32
#26層苟鸯，479點(diǎn)
tmp12=np.reshape(tmp11,(26,479))
#從一維列表變形成二維列表，注意python是行優(yōu)先
sp=np.mean(tmp12,axis=0)
#按層平均
return sp
#傳出動(dòng)能譜密度
def getticks(filename):
with open(filename,mode='r') as f:
data=f.readlines()
d1=[]
for line in data:
arr=list(map(float,line.split()))
d1.extend(arr)
#不知道為什么棚点，fortran存儲(chǔ)的文本文件變成了5個(gè)數(shù)值一行早处，所以用f.readlines()先把文本中所
#有的值讀到data中
#然后對(duì)data中每一行按空格分割(split),轉(zhuǎn)化為實(shí)數(shù)(float),針對(duì)整個(gè)列表實(shí)施操作，轉(zhuǎn)成新列表
d2=np.reshape(d1,(3,479))
return d2[1][:],d2[2][:]
#傳出標(biāo)準(zhǔn)化角波數(shù)和波長(zhǎng)
def main():
sp_noskeb=getbin("spectra_noskeb.dat")
sp_skeb=getbin("spectra_skeb.dat")
wvnum,wvlen=getticks("aol.txt")
wvlen=wvlen/1000.
line1=wvnum**(-5/3)/(10.0**2.9)
line2=wvnum**(-3)/(10.0**9.1)
#繪制log圖上斜率為-5/3和-3的直線數(shù)值
fig, ax =plt.subplots()
ax.loglog(wvnum,sp_noskeb,'b-',label="noskeb")
ax.loglog(wvnum,line1,'k:')
ax.loglog(wvnum,line2,'k:')
#繪制對(duì)數(shù)曲線圖
ax2=ax.twiny()
#雙橫軸必須操作
ax2.loglog(wvlen,sp_skeb,'r--',label="skeb")
plt.gca().invert_xaxis()
#帶數(shù)據(jù)翻轉(zhuǎn)橫軸（橫軸默認(rèn)從小到大排列瘫析，而非按數(shù)組次序顯示）
ax.legend(loc='upper right', shadow=True, fontsize='medium')
ax2.legend(loc='upper center', shadow=True, fontsize='medium')
#繪制圖例
fig.savefig('kes.eps',dpi=300)
#保存eps格式圖形
plt.show()
#在屏幕上顯示圖形
if __name__ == '__main__':
main()
最終的圖形如下：
