python计算位涡平流项
主要为了测试计算位涡平流项的准确性，这里使用了两种方法来计算位涡的平流项。
- 方法1：使用metpy函数
- 方法2：编写代码自己计算
一般使用metpy函数的计算结果，已经是等号右边的那一项了。也就是说相当于下图所示：
展开来就是：
使用metpy的函数时，需要注意的点：变量带上单位。
具体的代码如下所示：
import numpy as np import pandas as pd import xarray as xr import proplot as pplt import glob from matplotlib.colors import ListedColormap from wrf import getvar,pvo,interplevel,destagger,latlon_coords from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter from netCDF4 import Dataset import matplotlib as mpl from metpy.units import units import metpy.constants as constants import metpy.calc as mpcalc import cmaps f_w = r'D:/wind.nc' f_pv = r'D:/pv.nc' def cal_dxdy(file): ncfile = Dataset(file) P = getvar(ncfile, "pressure") lats, lons = latlon_coords(P) lon = lons[0] lon[lon<=0]=lon[lon<=0]+360 lat = lats[:,0] dx, dy = mpcalc.lat_lon_grid_deltas(lon.data, lat.data) return lon,lat,dx,dy path = r'J:/wrfout/' filelist = glob.glob(path+'wrf*') filelist.sort() lon,lat,dx,dy = cal_dxdy(filelist[-1]) dw = xr.open_dataset(f_w) dp = xr.open_dataset(f_pv) u = dw.u v = dw.v pv = dp.pv u = u.data*units('m/s') v = v.data*units('m/s') units.define('PVU = 1e-6 m^2 s^-1 K kg^-1 = pvu') pv = pv.data*units('pvu').to('m^2 s^-1 K kg^-1') lon = dw.lon lat = dw.lat time= dw.time leve= dw.level xlon,ylat=np.meshgrid(lon,lat) dlony,dlonx=np.gradient(xlon) dlaty,dlatx=np.gradient(ylat) pi=3.14159265 re=6.37e6 dx=re*np.cos(ylat*pi/180)*dlonx*pi/180 dy=re*dlaty*pi/180 adv = np.full((time.shape[0],leve.shape[0],lat.shape[0],lon.shape[0]),np.nan)*units('m/s')*units("m^2 s^-1 K kg^-1")/units("m") for i in range(time.shape[0]): for j in range(leve.shape[0]): print(i,j) adv[i,j,:,:] = mpcalc.advection(pv[i,j,:,:],u=u[i,j,:,:],v=v[i,j,:,:],dx=dx,dy=dy,x_dim=-1,y_dim=-2) ########################## way 2 ########################################## u0 = dw.u.data v0 = dw.v.data dpdx = np.gradient(pv,axis=-1) dpdy = np.gradient(pv,axis=-2) pvadv = np.full((time.shape[0],leve.shape[0],lat.shape[0],lon.shape[0]),np.nan) for i in range(time.shape[0]): for j in range(leve.shape[0]): print(i,j) pvadv[i,j] = -(u0[i,j]*dpdx[i,j]/dx+v0[i,j]*dpdy[i,j]/dy) f, axs = pplt.subplots( ncols=1, nrows=2, figsize=(8,6), tight=True, proj='cyl', proj_kw={'lon_0': 180}, lonlim=(100, 210), latlim=(-25, 25), share=4, ) axs[0].contourf(lon[::4],lat[::4],adv[16,6][::4,::4],colorbar='b') axs[1].contourf(lon[::4],lat[::4],pvadv[16,6][::4,::4],colorbar='b') axs.format(title=time[16],titleloc='l', coast=True, labels=True, fontsize=10, )
结果对比，更方面的是直接绘制折线图。选取两个同样的点即可
