Metpy_⽓象物理量计算（相对湿度、露点温度、湿位涡等）1、湿位涡剖⾯分析
背景介绍：
在暴⾬发⽣前期，形成暴⾬的基本条件逐渐形成甚⾄完全具备。通过对形成暴⾬的基本条件即：⽔汽条件、不稳定能量条件、上升运动条件等诊断分析，有助于判断暴⾬发⽣的可能性。形成暴⾬的主要物理条件有两个：内在因素是潮湿空⽓的潜在不稳定，⽽以充⾜的⽔汽表现为其主要⽅⾯，简称热⼒条件；外部因素是促使这种潜在不稳定得到充分释放的强迫抬升运动，⽽⼜以流场的配置为其主要⽅⾯，简称动⼒条件。有的把其分为三个条件，即把热⼒条件分为⽔汽和潜在位势不稳定两个条件。
湿位涡（Moist Potential Vorticity, MPV）是表征动⼒热⼒作⽤的综合诊断物理量，给出了⼤⽓短期⾏为的热⼒状态和涡旋运动之间的约束关系，这种关系导致了强降⽔这样的天⽓现象中涡旋爆发性增长的重要机制，它的⼤⼩与⼤⽓层结的状态、斜压性以及风的垂直切变有关，其正负符号取决这三者的配置。
⽬录
'''
Author: your name
Date: 2021-08-10 11:54:21
LastEditTime: 2021-08-10 17:40:49
LastEditors: Please set LastEditors
Description: In User Settings Edit
FilePath: \Test\QHTJ\Python⽓象数据处理与绘图\Shiweiwo.py
'''
s as ccrs
import cartopy.feature as cfeature #添加地理更多信息
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.interpolate import cross_section
铜仁市中考成绩查询from metpy.units import units
stants import g
stants import g
def plt_tem(data):
fig, ax = plt.subplots(figsize=(20,5), subplot_kw=dict(projection=data_crs))
data['Temperature_isobaric'].metpy.sel(vertical=1000 * units.hPa, time='2018-02-01 12:00').plot(ax=ax)
ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='k', alpha=0.2, zorder=2) #添加州的地理信息
plt.show()
def cal_vort(data):
#计算地旋风
heights = data['Geopotential_height_isobaric']#等压位势⾼度
f = iolis_parameter(data['lat'])#包含运动学参数的计算(例如，散度或涡度)
dx, dy = id_deltas_from_dataarray(heights)惠普电脑售后
ug, vg = strophic_wind(heights, f, dx, dy)#计算从⾼度或重⼒势给出的地转风
#计算绝对涡度
vert_abs_vort = f + mpcalc.vorticity(ug, vg, dx, dy)#计算⽔平风的垂直涡度
return ug,vg, vert_abs_vort
# 计算热动⼒物理量
def cal_wuli(data):
temperature, pressure, specific_humidity = xr.broadcast(data['Temperature_isobaric'],data['isobaric'],
陈冠希lucasdata['Specific_humidity_isobaric'])
# 相对湿度
rh = lative_humidity_from_specific_humidity(specific_humidity, temperature, pressure)
# 露点温度
dewpoint = mpcalc.dewpoint_from_specific_humidity(specific_humidity, temperature, pressure)
# 相当位温
theta_e = mpcalc.equivalent_potential_temperature(pressure, temperature, dewpoint)
return rh,dewpoint,theta_e
data = xr.open_dataset(r'D:\Pythonbase\Test\QHTJ\PythonDraw\NARR3D_201802_').metpy.parse_cf()
data_crs = data['Temperature_isobaric'].metpy.cartopy_crs
print("投影信息为：",data_crs)
heights = data['Geopotential_height_isobaric']#等压位势⾼度
ug,vg, vert_abs_vort=cal_vort(data)
rh,dewpoint,theta_e=cal_wuli(data)
data['theta_e'] = xr.DataArray(theta_e,
ds,
dims=heights.dims,
attrs={'units': theta_e.units})
data['u_g'] = xr.DataArray(ug,
收费员岗位职责ds,
dims=heights.dims,
attrs={'units': ug.units})
data['v_g'] = xr.DataArray(vg,
ds,
dims=heights.dims,
attrs={'units': vg.units})
data['rh'] = xr.DataArray(rh,
ds,
dims=heights.dims,
attrs={'units': rh.units})
data['rh'].vert_units('percent')
# 计算湿位涡
dtheta_e_dp, dtheta_e_dy, dtheta_e_dx = (py.unit_array for var adient(data['theta
_e'], axes=('vertical', 'y', 'x'))) dug_dp = mpcalc.first_derivative(data['u_g'], axis='vertical').metpy.unit_array
dvg_dp = mpcalc.first_derivative(data['v_g'], axis='vertical').metpy.unit_array
dz_dp = mpcalc.first_derivative(data['Geopotential_height_isobaric'], axis='vertical').metpy.unit_array
mpv = g * (1 / dz_dp) * (-dvg_dp * dtheta_e_dx + dug_dp * dtheta_e_dy + vert_abs_vort * dtheta_e_dp)
data['mpv'] = xr.DataArray(mpv,
ds,
dims=heights.dims,
attrs={'units': mpv.units})
data['mpv'].vert_units('microkelvin / s^3')
data['mpv'].vert_units('microkelvin / s^3')
# Cross section parameters，剖⾯分析
start = (35.49, -111.17)
end = (42.75, -98.26)
time = '2018-02-01 12:00'
cross = cross_section(data.sel(time=time), start, end)
# 计算绝对地转动量
momentum = mpcalc.absolute_momentum(cross['u_g'], cross['v_g'])
print(momentum)
孟非的老婆fig = plt.figure(1, figsize=(14., 6.))
ax = plt.axes()
mpv_contour = ax.contourf(cross['lon'], cross['isobaric'], cross['mpv'],
吧台设计
levels=np.arange(-120, 121, 10), cmap='bwr')
mpv_colorbar = lorbar(mpv_contour)
thetae_contour = ax.contour(cross['lon'], cross['isobaric'], cross['theta_e'],
levels=np.arange(250, 450, 5), colors='k', linewidths=1,
linestyles='-', alpha=0.5, zorder=2)
thetae_contour.clabel(thetae_contour.levels[1::2], fontsize=8, colors='k', inline=1,
inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True)
rh_contour = ax.contour(cross['lon'], cross['isobaric'], cross['rh'],
levels=np.arange(70, 100, 10), colors='k', linewidths=2,
linestyles=':', alpha=0.8, zorder=3)
rh_contour.clabel(rh_contour.levels[1::2], fontsize=8, colors='k', inline=1,
inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True)
thetae_contour = ax.contour(cross['lon'], cross['isobaric'], momentum,
levels=np.arange(0, 150, 10), colors='k', linewidths=1,
linestyles='--', alpha=0.5, zorder=2)
thetae_contour.clabel(thetae_contour.levels[1::2], fontsize=8, colors='k', inline=1,
inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True)
ax.set_yscale('symlog')
ax.set_yticklabels(np.arange(1000, 50, -100))
ax.set_ylim(cross['isobaric'].max(), cross['isobaric'].min())
ax.set_yticks(np.arange(1000, 50, -100))
ax_inset = fig.add_axes([0.058, 0.630, 0.25, 0.25], projection=data_crs)
ur(data['x'], data['y'], data['Geopotential_height_isobaric'].sel(time=time, isobaric=500.), levels=np.arange(5100, 6000, 60), cmap='inferno')
endpoints = ansform_points(ccrs.Geodetic(),
*np.vstack([start, end]).transpose()[::-1])
ax_inset.scatter(endpoints[:, 0], endpoints[:, 1], c='k', zorder=2)
ax_inset.plot(cross['x'], cross['y'], c='k', zorder=2)
pad = 1e6
ax_inset.set_extent([cross['x'][0] - pad, cross['x'][-1] + pad,
cross['y'][0] - pad, cross['y'][-1] + pad], crs=data_crs)
astlines()
ax_inset.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='k', alpha=0.2, zorder=0)
ax_inset.set_title('')
ax.set_title('NARR Cross-Section \u2013 {} to {} \u2013 Valid: {}\n'
'MPV, Relative Humidity (70%, 80%, 90% contours dotted), Theta-E (K, solid), '
'Absolute Momentum (m/s, dashed)\n'
'Inset: Cross-Section Path and 500 hPa Geopotential Height'.format(
start, end, cross['time'].dt.strftime('%Y-%m-%d %H:%MZ').item()))
ax.set_ylabel('Pressure (hPa)')
ax.set_xlabel('Longitude (degrees east)')
ax.set_xlabel('Longitude (degrees east)')
mpv_colorbar.set_label('Moist Geostrophic Potential Vorticity (microkelvins / s ** 3)') plt.show()
发布评论