OW参数计算

• 涡旋OW参数计算
• 一般选取threshold -2×10-12
• AE与CE的ow值均为复制，因而仅参考ow参数无法判断AE还是CE

import netCDF4
d_step=0.08333588
earth_radius=6370997.0

'''
deg2rad将角度单位转化为弧度单位
d_ = 360 / (0.083 * cos(deg2rad(y[i])) * pi * 2 * earth_radius)
(h0-h1)*d
'''

'''偏导数和梯度不是一个东西！'''
def ow_x(h,lat):
	'''求x方向上的偏导数'''
    length=len(h[0])
    h[abs(h)>100]=np.nan
    h_new=np.zeros((len(h),length))
    
    rad=deg2rad(lat)
    for j in range(1,length-1):
        d=360/(d_step*2*np.cos(rad) * pi * 2 * earth_radius)
        one=h[:,j+1]-h[:,j-1]
        two=one*d
        h_new[:,j]=two
    return h_new

def ow_y(h,lat):
	'''求y方向上的偏导数'''
    length=len(h)
    h[abs(h)>100]=np.nan
    h_new=np.zeros((len(h),len(h[0])))
    
    for j in range(1,length-1):
        rad=deg2rad(lat[j])
        d=360/(d_step*2*np.cos(rad) * pi * 2 * earth_radius)
        one=h[j+1,:]-h[j-1,:]
        two=one*d
        h_new[j,:]=two
    return h_new

'''数据读取'''
file='D:\\data\\cmems_mod_glo_phy_my_0.083_P1M-m_1692762273151.nc'
data=netCDF4.Dataset(file)

v_all=data.variables['vo'][:].data[0,0,:,:]
u_all=data.variables['uo'][:].data[0,0,:,:]
x=data.variables['longitude'][:].data
y=data.variables['latitude'][:].data
depth=data.variables['depth'][:].data
time=data.variables['time'][:].data

sn=(u_x-v_y)
ss=(v_x+u_y)
w=(v_x-u_y) #涡度，用协助判断AE与CE
W=sn**2+ss**2-w**2 #所求的OW参数

X,Y=np.meshgrid(x,y)
'''绘制OW参数'''
plt.contourf(X,Y,W)