- 积分
- 3285
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2021-2-27
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
import xarray as xr
import numpy as np
import metpy.calc as mpcalc
from metpy.units import units
import matplotlib.pyplot as plt
from cartopy.mpl.ticker import LongitudeFormatter
#################################################################################################
all_vars = xr.open_dataset('D:\\python\\tianzhen\\shixi2\\all.nc')
lons = all_vars['lon'][:]
lats = all_vars['lat'][:]
levels = all_vars['level'][:]
lon, lat = np.meshgrid(lons, lats)
lon=lon*units.degrees_east
lon, level = np.meshgrid(lons, levels)
level=level*units.hPa
####赋上单位,度数
lons = lons * units.degrees_east
lats = lats * units.degrees_north
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats) ###后面散度计算需要
uwind = all_vars['u'][0, 0, :, :]
vwind = all_vars['v'][0, 0, :, :]
uwind_levels = all_vars['u'][0, :, :, :]
vwind_levels = all_vars['v'][0, :, :, :]
####赋上单位,m/s
uwind_levels = uwind_levels * (units.m / units.s)
vwind_levels = vwind_levels * (units.m / units.s)
uwind = uwind * (units.m / units.s)
vwind = vwind * (units.m / units.s)
###################################################################################################
Pcha = [75, 75, 150, 100, 100]
div = np.zeros((6, 29, 45))
W_speed = np.zeros((6, 29, 45))
W_speed_Xiu = np.zeros((6, 29, 45))
for h_i in range(6):
# 计算散度
div[h_i, :, :] = mpcalc.divergence(uwind_levels[h_i, :, :], vwind_levels[h_i, :, :], dx=dx[:, :], dy=dy[:, :])
for h_i in range(1, 6):
W_speed[h_i, :, :] = W_speed[h_i - 1, :, :] + 0.5 * (div[h_i, :, :] - div[h_i - 1, :, :]) * Pcha[h_i - 1]
for h_i in range(5):
W_speed_Xiu[h_i, :, :] = W_speed[h_i, :, :] - h_i * (h_i + 1) / (7 * 6) * (W_speed[5, :, :] - W_speed_Xiu[5, :, :])
#############################################################绘图##################################3
fig = plt.figure()
ax = fig.subplots(1, 1)
plt.contourf(lon[:, :], level[:, :], W_speed_Xiu[:, 11, :])
plt.colorbar()
#####横坐标显示经度
lon_formatter = LongitudeFormatter(zero_direction_label=False)
ax.xaxis.set_major_formatter(lon_formatter)
####纵坐标格式显示
ax.set_yscale('log')
ax.set_ylim(1000,500)
ax.set_ylabel('hPa')
plt.show()
|
-
|