哪位看一下程序错在哪
import xarray as xr
import numpy as np
import os
os.chdir("D:\IDM")
os.getcwd()
from math import pi
def spectral_analysis(X,lon):
import numpy as np
delta_lon=2*np.pi/lon #lon为纬圈格点数
a0=np.mean(X) #0波即纬圈平均
a1_12=np.zeros((13,1))
b1_12=np.zeros((13,1))#傅里叶系数
#计算傅里叶系数
X=np.array(X)
for j in range(13):
for i in range(lon):
a1_12[j]=a1_12[j]+X[i]*np.cos(j*delta_lon*i)
b1_12[j]=b1_12[j]+X[i]*np.sin(j*delta_lon*i)
#除以空间序列的长度
a1_12=a1_12*2/lon
b1_12=b1_12*2/lon
a1_12[0]=a1_12[0]/2
#print(a1_12)
#print(b1_12)
#累加各个谐波
X1_12=np.zeros((13,len(X)))
for j in range(1,13):
for i in range(lon):
X1_12[j,i]=a0+a1_12[j]*np.cos(j*i*delta_lon)+b1_12[j]*np.sin(j*i*delta_lon)
#0波
X1_12[0]=a1_12[0]
#1-3波
X1_3=np.zeros((len(X)))
for i in range(lon):
X1_3[i]=X1_12[1,i]+X1_12[2,i]+X1_12[3,i]-3*a0
return X1_3
if __name__ == '__main__':
#读取数据
data_file=xr.open_dataset('0901.nc')
#print(data_file)
U=data_file['u'].loc['2009-01-27':'2009-01-29',:,:,:]
V=data_file['v'].loc['2009-01-27':'2009-01-29',:,:,:]
T=data_file['t'].loc['2009-01-27':'2009-01-29',:,:,:]
#谐波分析
Uza=np.zeros((3,37,91,360))
Vza=np.zeros((3,37,91,360))
Tza=np.zeros((3,37,91,360))
for t in range(3):
for k in range(37):
for j in range(91):
Uza[t,k,j,:]=spectral_analysis(U[t,k,j,:],360)
Vza[t,k,j,:]=spectral_analysis(V[t,k,j,:],360)
Tza[t,k,j,:]=spectral_analysis(T[t,k,j,:],360)
print(Uza)