
直接计算出来的势函数在高层<0。我画图的时候加了负号。(不知道对不对)
import xarray as xr
import easyclimate as ec
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import cartopy.feature as cfeature
import matplotlib as mpl
mpl.rcParams["font.sans-serif"] = ["SimHei"]
mpl.rcParams["axes.unicode_minus"] = False
u= xr.open_dataset(r'G:\zbz\code\datasets\uwnd.mon.mean.nc').uwnd
v= xr.open_dataset(r'G:\zbz\code\datasets\vwnd.mon.mean.nc').vwnd
clevel=[np.arange(-20, 24, 4)*1e+6, np.arange(-8,10, 1)*1e+7]
# u_mam = u[u.time.dt.month.isin([3,4,5])].mean(axis=0)
# v_mam = v[v.time.dt.month.isin([3,4,5])].mean(axis=0)
u_jja = u[u.time.dt.month.isin([6,7,8])].mean(axis=0)
v_jja = v[v.time.dt.month.isin([6,7,8])].mean(axis=0)
# u_son = u[u.time.dt.month.isin([9,10,11])].mean(axis=0)
# v_son = v[v.time.dt.month.isin([9,10,11])].mean(axis=0)
# u_djf = u[u.time.dt.month.isin([12,1,2])].mean(axis=0)
# v_djf = v[v.time.dt.month.isin([12,1,2])].mean(axis=0)
# u_jja = u_jja-u_jja.mean(axis=-1)
# v_jja = v_jja-v_jja.mean(axis=-1)
# clevel=[np.arange(-10, 12, 2)*1e+6, np.arange(-2, 2.4, 0.4)*1e+7]
datas = ec.core.windspharm.calc_helmholtz(u_jja,v_jja)
datas2 = ec.core.windspharm.calc_streamfunction_and_velocity_potential(u_jja,v_jja)
# datas = ec.core.windspharm.calc_nondivergent_component(u_jja,v_jja)
# psi--旋转风 chi ---辐散风
datas = datas.sel(lat=slice(70,-10) , lon=slice(20,170) )
datas2 = datas2.sel(lat=slice(70,-10) , lon=slice(20,170) )
x, y = np.meshgrid(datas.lon, datas.lat)
nnpp = 3 ; ucolor='b'
scale1,scale2, width, wsize = [50,20, 25, 25],[260,100, 80, 80], 0.0050, 4
fig, axs = plt.subplots(4,2,subplot_kw={'projection': ccrs.PlateCarree(),}, dpi=400, figsize=(8,10), sharex=True, sharey=True)
fig.suptitle('JJA')
for ii, lev in enumerate([ 200,500,700,850]):
pv = datas2.pv.sel(level=lev)
uchi = datas.uchi.sel(level=lev)
vchi = datas.vchi.sel(level=lev)
q=axs[ii, 0].quiver(x[::nnpp, ::nnpp], y[::nnpp, ::nnpp], uchi[::nnpp, ::nnpp], vchi[::nnpp, ::nnpp] , color=ucolor,zorder=2, transform=ccrs.PlateCarree() , scale=scale1[ii], width=width)
axs[ii, 0].quiverkey(q,0.90,1.02, wsize,label=f'{wsize}m/s', coordinates = "axes")
axs[ii, 0].add_feature(cfeature.COASTLINE)
axs[ii, 0].set_title(f'{lev}hPa chi 辐散风+势函数')
cf1 = axs[ii, 0].contourf(x,y,-pv,levels=clevel[0], cmap='RdBu_r',transform=ccrs.PlateCarree(),zorder=1, extend='both')
sf = datas2.stream.sel(level=lev)
upsi = datas.upsi.sel(level=lev)
vpsi = datas.vpsi.sel(level=lev)
q=axs[ii, 1].quiver(x[::nnpp, ::nnpp], y[::nnpp, ::nnpp], upsi[::nnpp, ::nnpp], vpsi[::nnpp, ::nnpp] , color=ucolor,zorder=2,transform=ccrs.PlateCarree() , scale=scale2[ii], width=width)
axs[ii, 1].quiverkey(q,0.90,1.02, wsize,label=f'{wsize}m/s', coordinates = "axes")
axs[ii, 1].add_feature(cfeature.COASTLINE)
axs[ii, 1].set_title(f'{lev}hPa psi 旋转风+流函数')
cf2 = axs[ii, 1].contourf(x,y,sf,levels=clevel[1], cmap='RdBu_r',transform=ccrs.PlateCarree(),zorder=1, extend='both')
cax1 = axs[3, 0].inset_axes([0.1, -0.15, 0.8, 0.05]) # 调整位置和大小
cax2 = axs[3, 1].inset_axes([0.1, -0.15, 0.8, 0.05]) # 调整位置和大小
# 添加色标
fig.colorbar(cf1, cax=cax1, orientation='horizontal', aspect=40)
fig.colorbar(cf2, cax=cax2, orientation='horizontal', aspect=40)
plt.tight_layout()
plt.show()