py计算流函数势函数旋转风辐散风
翊清
2025年03月26日 15:42
收录于文集
共43篇

直接计算出来的势函数在高层<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()