python计算位涡平流项

主要为了测试计算位涡平流项的准确性，这里使用了两种方法来计算位涡的平流项。

方法1：使用metpy函数
方法2：编写代码自己计算

一般使用metpy函数的计算结果，已经是等号右边的那一项了。也就是说相当于下图所示：

展开来就是：

−-−(uuu∂pv∂x\frac{\partial pv}{\partial x}∂x∂pv+vvv∂pv∂y\frac{\partial pv}{\partial y}∂y∂pv)

使用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 cmapsf_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]+360lat = 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.pvu = 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.levelxlon,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/180adv = 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,)

结果对比，更方面的是直接绘制折线图。选取两个同样的点即可

python--测试使用不同的方式计算位涡平流项的差异相关推荐

python 基于metpy计算位涡平流项（水平）
使用python 计算位涡水平平流项使用数据:potential vorticity.水平风场分量:U.V 数据来源: ERA5 hourly :一天四次空间分辨率:0.5x0.5 时间范围:20 ...
python求阶乘之和_python计算阶乘前n项和
广告关闭腾讯云11.11云上盛惠 ,精选热门产品助力上云,云服务器首年88元起,买的越多返的越多,最高返5000元! 知道公式后就很简单了,利用for循环,第几行i+1就等于几,当然python中是 ...
python计算小数点后有几位_数学提高小数除法竖式计算过程
除数是小数时:先把除数变成整数,除数扩大到原来的多少倍,被除数也要扩大到原来的多少倍(即小数点也向右移动几位,位数不够的用0补),然后按照除数是整数的除法进行计算.计算小数除法,除到被除数的末尾仍有余 ...
Python模拟大整数乘法的小学竖式计算过程
让我们先看个图回顾一下小学学过的计算整数乘法的竖式计算过程然后再来看如何使用Python来模拟上面的过程,虽然在Python中计算任意大的数字乘法都没有问题,但下面的代码作为一个算法的理解还是不错的 ...
保留两位小数除法算式_两位小数除法练习题竖式计算
精品文档 2016 全新精品资料 - 全新公文范文 - 全程指导写作 – 独家原创 1 / 10 两位小数除法练习题竖式计算一.一般乘法竖式计算题 65×0.0.0016×10.65×0.1 0 . ...
python 流式计算框架_流式计算的三种框架：Storm、Spark和Flink
我们知道,大数据的计算模式主要分为批量计算(batch computing).流式计算(stream computing).交互计算(interactive computing).图计算(graph ...
Python生成三年级数学竖式计算并用Python-pptx写入PPT
三年级下册的学生需要经常练习两位数乘两位数,三位数除以一位数的竖式计算,虽然出题难题简单,但是如果可以用程序一下生成,岂不是一劳永逸.那就动手吧. 首先需要安装Python-pptx,安装过程可能会因 ...
全面对比 MATLAB、Julia、Python，谁在科学计算中更胜一筹？
数百种编程语言,各有优劣,各自也都有自己最为适用的场景.那么就科学计算领域而言,主流的 MATLAB.Julia.Python 会有哪些最为独特的优势呢?又存在哪些让开发者无力的缺陷?在本文中,我们将 ...
python 对比matlab_全面对比 MATLAB、Julia、Python，谁在科学计算中更胜一筹？
原标题:全面对比 MATLAB.Julia.Python,谁在科学计算中更胜一筹? 数百种编程语言,各有优劣,各自也都有自己最为适用的场景.那么就科学计算领域而言,主流的 MATLAB.Julia.P ...

python--测试使用不同的方式计算位涡平流项的差异

python计算位涡平流项

python--测试使用不同的方式计算位涡平流项的差异相关推荐

最新文章

热门文章