工程建设信息网站百度搜不干净的东西

写在前面

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

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

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

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

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

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

metpy

def _cal_divergence_metpy(u,v):import metpy.calc as mpcalclon = u.lon.datalat = u.lat.datadx, 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 nplon = u.lon.datalat = u.lat.dataxlon,ylat=np.meshgrid(lon,lat)dlony,dlonx=np.gradient(xlon)dlaty,dlatx=np.gradient(ylat)pi=np.pire=6.37e6dx=re*np.cos(ylat*pi/180)*dlonx*pi/180dy=re*dlaty*pi/180u_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]/dyreturn   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 + dvdyreturn xr.DataArray(div1,dims=u.dims,coords=u.coords),xr.DataArray(div,dims=u.dims,coords=u.coords)