1、湿位涡剖面分析

数据处理于可视化 | 湿位涡剖面分析

背景介绍：

在暴雨发生前期，形成暴雨的基本条件逐渐形成甚至完全具备。通过对形成暴雨的基本条件即：水汽条件、不稳定能量条件、上升运动条件等诊断分析，有助于判断暴雨发生的可能性。形成暴雨的主要物理条件有两个：内在因素是潮湿空气的潜在不稳定，而以充足的水汽表现为其主要方面，简称热力条件；外部因素是促使这种潜在不稳定得到充分释放的强迫抬升运动，而又以流场的配置为其主要方面，简称动力条件。有的把其分为三个条件，即把热力条件分为水汽和潜在位势不稳定两个条件。

湿位涡（Moist Potential Vorticity, MPV）是表征动力热力作用的综合诊断物理量，给出了大气短期行为的热力状态和涡旋运动之间的约束关系，这种关系导致了强降水这样的天气现象中涡旋爆发性增长的重要机制，它的大小与大气层结的状态、斜压性以及风的垂直切变有关，其正负符号取决这三者的配置。

目录

1、湿位涡剖面分析

背景介绍：

'''
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
'''
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
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 = 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)#计算水平风的垂直涡度
    return ug,vg, vert_abs_vort
# 计算热动力物理量
def cal_wuli(data):
    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)
    return rh,dewpoint,theta_e

data = xr.open_dataset(r'D:\Pythonbase\Test\QHTJ\PythonDraw\NARR3D_201802_0103.tar.nc').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,
                               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')

# 计算湿位涡
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)
# 计算绝对地转动量
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 = 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()