Python气象数据处理与绘图：探空图

1、导入模块

1import matplotlib.pyplot as plt
2from mpl_toolkits.axes_grid1.inset_locator import inset_axes
3import pandas as pd
4import metpy.calc as mpcalc
5from metpy.plots import Hodograph, SkewT
6from metpy.units import units
7import metpy

2、读取探空数据

1col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
2df = pd.read_fwf('/home/kesci/work/data/sounding.txt',
3                 skiprows=4, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
4
5df.head()

输出：

3、数据预处理

删除缺测值

1df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'), how='all').reset_index(drop=True)
2df

添加单位

1p = df['pressure'].values * units.hPa                # 单位：hPa
2T = df['temperature'].values * units.degC            # 单位：℃
3Td = df['dewpoint'].values * units.degC              # 单位：℃
4wind_speed = df['speed'].values * units.knots        # 单位：knot
5wind_dir = df['direction'].values * units.degrees    # 单位：°
6u, v = mpcalc.wind_components(wind_speed, wind_dir)  # 计算水平风速u和v

4、绘制初级探空图

 1#创建画布
 2fig = plt.figure(figsize=(9, 9))
 3skew = SkewT(fig)
 4
 5skew.plot(p, T, 'r', linewidth=2)   # 绘制层结曲线（红色）
 6skew.plot(p, Td, 'g', linewidth=2)  # 绘制露点曲线（绿色）
 7skew.plot_barbs(p, u, v)            # 绘制风羽
 8
 9skew.ax.set_ylim(1000, 200)   
10skew.ax.set_xlim(-20, 60)
11
12#显示图像
13plt.show()

5、计算凝结抬升高度（LCL）

1lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
2print(lcl_pressure, lcl_temperature)

输出：

1948.1067995872485 hectopascal 25.49950229677085 degree_Celsius

6、计算状态曲线

1parcel_prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')
2print(parcel_prof)

输出：

1[30.600000000000023 30.081285047119707 27.8957275238335 26.386579497637854 24.701701336743895 24.207264088485204 22.932228947805072 21.94652843486716 17.117767018299787 16.232958179045454 15.947577973544185 15.46531064377848 14.875265027129615 13.862927741368594 11.719120534481135 8.467037171116658 6.175019720755472 5.6924379448918785 5.342399715664499 4.337526058990079 3.4454961699111095 0.4940733134725406 -4.792064592204554 -5.413042976360941 -5.728875428476499 -7.932176286834704 -10.95643300525603 -11.99683868276162 -15.965309555309318 -18.497055521292936 -24.27413777228506 -27.371905467383982 -28.02547693992861 -34.73762978040941 -40.88142948477113] degree_Celsius

7、绘制中级探空图

 1#创建画布
 2fig = plt.figure(figsize=(9, 9))
 3skew = SkewT(fig, rotation=30)
 4
 5skew.plot(p, T, 'r')          # 绘制层结曲线（红色实线）
 6skew.plot(p, Td, 'g')         # 绘制等饱和混合比线 （绿色实线）
 7skew.plot_barbs(p, u, v)      # 绘制风羽
 8skew.ax.set_ylim(1000, 200)   
 9skew.ax.set_xlim(-20, 60)
10
11# 绘制凝结抬升高度（黑点）
12skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
13# 绘制气块状态曲线（黑线）
14skew.plot(p, parcel_prof, 'k', linewidth=2)
15
16# 对流有效位能（CAPE）：正不稳定能量，是气块达到自由对流高度后可以获得继续上升的能量
17# 对流抑制能量（CIN）：负不稳定能量，是大气底层的气块要能达到自由对流高度
18skew.shade_cin(p, T, parcel_prof, Td)  # 阴影绘制CIN，蓝色
19skew.shade_cape(p, T, parcel_prof)     # 阴影绘制CAPE，红色
20
21skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
22
23# 添加关键曲线
24skew.plot_dry_adiabats()     # 干绝热线（红色虚线）
25skew.plot_moist_adiabats()   # 湿绝热线（蓝色虚线）
26skew.plot_mixing_lines()     # 等饱和混合比线（绿色虚线）
27
28# 显示图像
29plt.show()

8、绘制高级探空图

 1# 创建画布
 2fig = plt.figure(figsize=(9, 9))
 3skew = SkewT(fig, rotation=30)
 4
 5skew.plot(p, T, 'r')          # 绘制层结曲线（红色实线）
 6skew.plot(p, Td, 'g')         # 绘制等饱和混合比线 （绿色实线）
 7skew.plot_barbs(p, u, v)      # 绘制风羽
 8skew.ax.set_ylim(1000, 200)   
 9skew.ax.set_xlim(-20, 60)
10
11# 绘制凝结抬升高度（黑点）
12skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
13# 绘制气块状态曲线（黑线）
14skew.plot(p, parcel_prof, 'k', linewidth=2)
15
16# 对流有效位能（CAPE）：正不稳定能量，是气块达到自由对流高度后可以获得继续上升的能量
17# 对流抑制能量（CIN）：负不稳定能量，是大气底层的气块要能达到自由对流高度
18skew.shade_cin(p, T, parcel_prof, Td)  # 阴影绘制CIN，蓝色
19skew.shade_cape(p, T, parcel_prof)     # 阴影绘制CAPE，红色
20
21skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
22
23# 添加关键曲线
24skew.plot_dry_adiabats()     # 干绝热线（红色虚线）
25skew.plot_moist_adiabats()   # 湿绝热线（蓝色虚线）
26skew.plot_mixing_lines()     # 等饱和混合比线（绿色虚线）
27
28# 创建高空风速分析图
29ax_hod = inset_axes(skew.ax, '40%', '40%', loc=1)
30h = Hodograph(ax_hod, component_range=80.)
31h.add_grid(increment=20)
32h.plot_colormapped(u, v, wind_speed)  # Plot a line colored by wind speed
33
34# 显示图像
35plt.show()

9、计算CAPE和CIN

CAPE is convective available potential energy
CIN is convective inhibition
LFC is pressure of the level of free convection
EL is pressure of the equilibrium level
SFC is the level of the surface or beginning of parcel path
Rd is the gas constant
g is gravitational acceleration
Tparcel is the parcel temperature
Tenv is environment temperature
p is atmospheric pressure

1cape, cin = metpy.calc.cape_cin(p, T, Td, parcel_prof, which_lfc='bottom', which_el='top')
2print(cape,cin)

输出：

12547.328781583579 joule / kilogram -104.3340479755981 joule / kilogram