背景介绍
在暴雨发生前期，形成暴雨的基本条件逐渐形成甚至完全具备。通过对形成暴雨的基本条件即：水汽条件、不稳定能量条件、上升运动条件等诊断分析，有助于判断暴雨发生的可能性。形成暴雨的主要物理条件有两个：内在因素是潮湿空气的潜在不稳定，而以充足的水汽表现为其主要方面，简称热力条件；外部因素是促使这种潜在不稳定得到充分释放的强迫抬升运动，而又以流场的配置为其主要方面，简称动力条件。有的把其分为三个条件，即把热力条件分为水汽和潜在位势不稳定两个条件。
湿位涡（Moist Potential Vorticity, MPV）是表征动力热力作用的综合诊断物理量，给出了大气短期行为的热力状态和涡旋运动之间的约束关系，这种关系导致了强降水这样的天气现象中涡旋爆发性增长的重要机制，它的大小与大气层结的状态、斜压性以及风的垂直切变有关，其正负符号取决这三者的配置。
康志明,罗金秀,郭文华,杨克明.2005年10月西藏高原特大暴雪成因分析[J].气象,2007,33(8):60-67
更多介绍：http://stream1.cmatc.cn/cmatcvod/12/tqx/third_chapter_fourth.html
导入模块
代码语言：javascript复制import cartopy.crs 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
from metpy.constants import g
下载数据
代码语言：javascript复制# 下载数据（约187MB）
# !wget http://www.meteor.iastate.edu/~jthielen/NARR3D_201802_0103.tar.nc
data = xr.open_dataset('./NARR3D_201802_0103.tar.nc').metpy.parse_cf()
data
<xarray.Dataset>
Dimensions: (isobaric: 29, time: 9, x: 260, y: 120)
Coordinates:
reftime (time) datetime64[ns] 2018-02-01 ... 2018-02-02
* x (x) float32 -3227081.8 ... 5180835.5
* y (y) float32 -2856844.2 ... 1006253.0
* time (time) datetime64[ns] 2018-02-01 ... 2018-02-02
* isobaric (isobaric) float64 100.0 ... 1e 03
lat (y, x) float64 20.06 ... 37.38
lon (y, x) float64 -135.0 ... -41.64
crs object Projection: lambert_conformal_conic
Data variables:
Pressure_Vertical_velocity_isobaric (time, isobaric, y, x) float32 ...
Temperature_isobaric (time, isobaric, y, x) float32 ...
Specific_humidity_isobaric (time, isobaric, y, x) float32 ...
u-component_of_wind_isobaric (time, isobaric, y, x) float32 ...
v-component_of_wind_isobaric (time, isobaric, y, x) float32 ...
Geopotential_height_isobaric (time, isobaric, y, x) float32 ...
LambertConformal_277X349-48p98N-106p51W int32 ...
绘制气温分布图
代码语言：javascript复制data_crs = data['Temperature_isobaric'].metpy.cartopy_crs
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)
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f937a975438>
计算地转风
heights = data['Geopotential_height_isobaric']
f = mpcalc.coriolis_parameter(data['lat'])[None, None, :, :]
dx, dy = mpcalc.grid_deltas_from_dataarray(heights)
ug, vg = mpcalc.geostrophic_wind(heights, f, dx, dy)
计算绝对涡度
代码语言：javascript复制vert_abs_vort = f mpcalc.vorticity(ug, vg, dx, dy)
计算热动力物理量
代码语言：javascript复制temperature, pressure, specific_humidity = xr.broadcast(data['Temperature_isobaric'],data['isobaric'],
data['Specific_humidity_isobaric'])
# 相对湿度
rh = mpcalc.relative_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)
添加数据到DataArray
代码语言：javascript复制data['theta_e'] = xr.DataArray(theta_e,
coords=heights.coords,
dims=heights.dims,
attrs={'units': theta_e.units})
data['u_g'] = xr.DataArray(ug,
coords=heights.coords,
dims=heights.dims,
attrs={'units': ug.units})
data['v_g'] = xr.DataArray(vg,
coords=heights.coords,
dims=heights.dims,
attrs={'units': vg.units})
data['rh'] = xr.DataArray(rh,
coords=heights.coords,
dims=heights.dims,
attrs={'units': rh.units})
data['rh'].metpy.convert_units('percent')
data
<xarray.Dataset>
Dimensions: (isobaric: 29, time: 9, x: 260, y: 120)
Coordinates:
reftime (time) datetime64[ns] 2018-02-01 ... 2018-02-02
* x (x) float32 -3227081.8 ... 5180835.5
* y (y) float32 -2856844.2 ... 1006253.0
* time (time) datetime64[ns] 2018-02-01 ... 2018-02-02
* isobaric (isobaric) float64 100.0 ... 1e 03
lat (y, x) float64 20.06 ... 37.38
lon (y, x) float64 -135.0 ... -41.64
crs object Projection: lambert_conformal_conic
Data variables:
Pressure_Vertical_velocity_isobaric (time, isobaric, y, x) float32 ...
Temperature_isobaric (time, isobaric, y, x) float32 ...
Specific_humidity_isobaric (time, isobaric, y, x) float32 ...
u-component_of_wind_isobaric (time, isobaric, y, x) float32 ...
v-component_of_wind_isobaric (time, isobaric, y, x) float32 ...
Geopotential_height_isobaric (time, isobaric, y, x) float32 16524.309 ... 288.98267
LambertConformal_277X349-48p98N-106p51W int32 ...
theta_e (time, isobaric, y, x) float64 381.9 ... 319.9
u_g (time, isobaric, y, x) float64 -2.997 ... -1.706
v_g (time, isobaric, y, x) float64 -10.57 ... 5.118
rh (time, isobaric, y, x) float64 7.878 ... 82.56
计算湿位涡
代码语言：javascript复制dtheta_e_dp, dtheta_e_dy, dtheta_e_dx = (var.metpy.unit_array for var in mpcalc.gradient(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,
coords=heights.coords,
dims=heights.dims,
attrs={'units': mpv.units})
data['mpv'].metpy.convert_units('microkelvin / s^3')
剖面分析
代码语言：javascript复制# 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)
cross
<xarray.Dataset>
Dimensions: (index: 100, isobaric: 29)
Coordinates:
reftime datetime64[ns] 2018-02-01T12:00:00
time datetime64[ns] 2018-02-01T12:00:00
* isobaric (isobaric) float64 100.0 ... 1e 03
lat (index) float64 35.49 ... 42.75
lon (index) float64 -111.2 ... -98.26
crs object Projection: lambert_conformal_conic
x (index) float64 -3.885e 05 ... 7.171e 05
y (index) float64 -1.618e 06 ... -7.659e 05
* index (index) int64 0 1 2 3 ... 97 98 99
Data variables:
Pressure_Vertical_velocity_isobaric (isobaric, index) float64 -0.03118 ... 0.01984
Temperature_isobaric (isobaric, index) float64 205.3 ... 264.3
Specific_humidity_isobaric (isobaric, index) float64 5.691e-06 ... 0.001202
u-component_of_wind_isobaric (isobaric, index) float64 10.0 ... 3.303
v-component_of_wind_isobaric (isobaric, index) float64 -11.93 ... -7.291
Geopotential_height_isobaric (isobaric, index) float64 1.631e 04 ... 215.7
LambertConformal_277X349-48p98N-106p51W int32 ...
theta_e (isobaric, index) float64 396.3 ... 267.7
u_g (isobaric, index) float64 10.83 ... -4.716
v_g (isobaric, index) float64 -11.27 ... -1.837
rh (isobaric, index) float64 13.82 ... 61.07
mpv (isobaric, index) float64 -12.54 ... -8.709
计算绝对地转动量
代码语言：javascript复制momentum = mpcalc.absolute_momentum(cross['u_g'], cross['v_g'])
print(momentum)
<xarray.DataArray (isobaric: 29, index: 100)>
array([[125.385938, 122.74618 , 120.136569, ..., 81.058598, 81.382169,
81.720505],
[121.357594, 118.286225, 115.648805, ..., 77.753756, 78.219216,
78.914308],
[119.654559, 116.379984, 113.44996 , ..., 74.913288, 75.626414,
75.981311],
...,
[145.028151, 145.666248, 146.618955, ..., 100.727893, 100.923846,
100.736723],
[147.087985, 147.596429, 147.998006, ..., 101.771015, 102.217413,
102.597688],
[147.873402, 147.990474, 148.062657, ..., 105.112548, 105.638732,
105.355884]])
Coordinates:
reftime datetime64[ns] 2018-02-01T12:00:00
time datetime64[ns] 2018-02-01T12:00:00
* isobaric (isobaric) float64 100.0 125.0 150.0 175.0 ... 950.0 975.0 1e 03
lat (index) float64 35.49 35.57 35.65 35.73 ... 42.62 42.68 42.75
lon (index) float64 -111.2 -111.1 -110.9 ... -98.55 -98.4 -98.26
crs object Projection: lambert_conformal_conic
x (index) float64 -3.885e 05 -3.771e 05 ... 7.062e 05 7.171e 05
y (index) float64 -1.618e 06 -1.61e 06 ... -7.745e 05 -7.659e 05
* index (index) int64 0 1 2 3 4 5 6 7 8 9 ... 91 92 93 94 95 96 97 98 99
Attributes:
units: meter / second
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 = fig.colorbar(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)
ax_inset.contour(data['x'], data['y'], data['Geopotential_height_isobaric'].sel(time=time, isobaric=500.),
levels=np.arange(5100, 6000, 60), cmap='inferno')
endpoints = data_crs.transform_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)
ax_inset.coastlines()
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)')
mpv_colorbar.set_label('Moist Geostrophic Potential Vorticity (microkelvins / s ** 3)')
plt.show()
数据及脚本
样例数据寄脚本获取请在好奇心Log公众号后台回复湿位涡。
