写在前面

简单记录一下计算散度的方法，方便以后查找复制

包括三种计算方式，numpy、metpy、windspharm。

其中，numpy和metpy的方法进行了比较，结果比较一致。

windspharm的方法里面包含了两种方法， 一致是直接调用divergence()函数实现，另一个是先计算梯度再相加。

使用windspharm最简单，缺点就是需要使用全球的数据作为输入，而且在Linux上安装，我这里用的是模式数据就没去和metpy和numpy进行验证，只是记录作为一种方法。

• 比较意外的是，同样的数据，numpy竟然比metpy还快一点。

metpy

def _cal_divergence_metpy(u,v):
    import metpy.calc as mpcalc
    lon = u.lon.data
    lat = u.lat.data
    dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
    divergence = np.zeros((u.shape[0],lat.shape[0],lon.shape[0]))
    for i in range(u.shape[0]): 
        print(i)
        divergence[i] = mpcalc.divergence(u = u[i],v = v[i],dx = dx ,dy = dy)
    return xr.DataArray(divergence,dims=u.dims,coords=u.coords)

numpy

def _divergence_with_numpy(u,v):
    import numpy as np
    lon = u.lon.data
    lat = u.lat.data
    xlon,ylat=np.meshgrid(lon,lat)
    dlony,dlonx=np.gradient(xlon)
    dlaty,dlatx=np.gradient(ylat)
    pi=np.pi
    re=6.37e6
    dx=re*np.cos(ylat*pi/180)*dlonx*pi/180
    dy=re*dlaty*pi/180
    u_dx = np.gradient(u,axis=-1)
    v_dy = np.gradient(v,axis=-2)
    div_numpy = np.zeros((u.shape))
    
    div_numpy = u_dx/dx + v_dy[i]/dy
    return   xr.DataArray(div_numpy,dims=u.dims,coords=u.coords) 

windspharm

def _cal_divergence_windspharm(u,v):
    w       = VectorWind(u,v)
    div     = w.divergence()
    dudx, _ = w.gradient(u)
    _, dvdy = w.gradient(v)
    div1 = dudx + dvdy
    return xr.DataArray(div1,dims=u.dims,coords=u.coords),xr.DataArray(div,dims=u.dims,coords=u.coords)