数据处理于可视化 | 湿位涡剖面分析
背景介绍在暴雨发生前期，形成暴雨的基本条件逐渐形成甚至完全具备。通过对形成暴雨的基本条件即：水汽条件、不稳定能量条件、上升运动条件等诊断分析，有助于判断暴雨发生的可能性。形成暴雨的主要物理条件有两个：内在因素是潮湿空气的潜在不稳定，而以充足的水汽表现为其主要方面，简称热力条件；外部因素是促使这种潜在不稳定得到充分释放的强迫抬升运动，而又以流场的配置为其主要方面，简称动力条件。有的把其分为三个条件，即把热力条件分为水汽和潜在位势不稳定两个条件。
湿位涡（Moist Potential Vorticity, MPV）是表征动力热力作用的综合诊断物理量，给出了大气短期行为的热力状态和涡旋运动之间的约束关系，这种关系导致了强降水这样的天气现象中涡旋爆发性增长的重要机制，它的大小与大气层结的状态、斜压性以及风的垂直切变有关，其正负符号取决这三者的配置。
康志明,罗金秀,郭文华,杨克明.2005年10月西藏高原特大暴雪成因分析.气象,2007,33(8):60-67
更多介绍：http://stream1.cmatc.cn/cmatcvod/12/tqx/third_chapter_fourth.html
导入模块
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
下载数据
# 下载数据（约187MB）
# !wget <a href="http://www.meteor.iastate.edu/~jthielen/NARR3D_201802_0103.tar.nc" target="_blank">http://www.meteor.iastate.edu/~jthielen/NARR3D_201802_0103.tar.nc</a>
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 2018-02-01 ... 2018-02-02
* x (x) float32 -3227081.8 ... 5180835.5
* y (y) float32 -2856844.2 ... 1006253.0
* time (time) datetime64 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-106p51Wint32 ...
绘制气温分布图
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'])
dx, dy = mpcalc.grid_deltas_from_dataarray(heights)
ug, vg = mpcalc.geostrophic_wind(heights, f, dx, dy)
计算绝对涡度
vert_abs_vort = f + mpcalc.vorticity(ug, vg, dx, dy)
计算热动力物理量
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
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 2018-02-01 ... 2018-02-02
* x (x) float32 -3227081.8 ... 5180835.5
* y (y) float32 -2856844.2 ... 1006253.0
* time (time) datetime64 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-106p51Wint32 ...
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
计算湿位涡
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')
剖面分析
# 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 2018-02-01T12:00:00
time datetime64 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-106p51Wint32 ...
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
计算绝对地转动量
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 2018-02-01T12:00:00
time datetime64 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, 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, 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, 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(, 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().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( - pad, cross['x'][-1] + pad,
cross['y'] - 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公众号后台回复湿位涡。
文章来源于微信公众号：好奇心Log
页:
[1]