Valar Dohaeris
0%

Python processing sounding data and plotting

发表于 阅读次数：

背景

最近重新整理之前的分析结果，发现一些有意思的图，来分享出来。2020年美国飓风 Laura 临近登陆时，NOAA派出飞机投放探空仪观测，很快推特上便有了相应的T-lnP图，如下图。

而我国在2020年，也进行了某某计划，在台风森拉克过境后，即在台风外围投放了数十枚探空仪，获取到了非常密集的温压湿风探空曲线。这是当时的卫星云图(来源：韩国GK2A卫星)，可以看到当时台风主体已经进入越南境内，而释放探空仪的位置在海南岛东南近海。

真彩色云图

红外增强云图

读取数据

由于数据是excel的表格，我们用pandas来读取表单。

1
2
3
4
5
6
7
8
9
10
11
12
from datetime import datetime,timedelta
import pandas as pd
f = '17_54_4.xls'
# 将时间解析为datetime类
dateparse = lambda x: datetime.strptime('2020-08-02 ' + x, '%Y-%m-%d %H:%M:%S')
df = pd.read_excel(f,sheet_name=0,parse_dates=['时间'], date_parser=dateparse)
# 时间         纬度          经度     位势高度     温度     湿度     气压  风速  风向
#                                     m           degC     %       hPa   m/s  deg   
# convert BJT -> UTC
df['时间'] = df['时间'] - timedelta(hours=8)
df[df == -9999] = np.nan
df = df[::-1]

读取结果

绘制T-lnP图

用NCAR开发的Metpy包来画T-lnP图，官方链接为：https://unidata.github.io/MetPy/latest/
官方提供了很多例子，包括画常规天气图和T-lnP，以及计算物理量、诊断量。链接：https://unidata.github.io/MetPy/latest/examples/index.html

根据Simple_Sounding例子，来画我们需要的最简单的T-lnP。注意：美国人用的是斜T-lnP，即旋转30度，而国内是T-lnP，需要改成0度。对应的函数接口为：
https://unidata.github.io/MetPy/latest/api/generated/metpy.plots.SkewT.html#metpy.plots.SkewT

SkewT的函数说明

下面为具体绘图代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.plots import SkewT, Hodograph
from metpy.units import units
import metpy.interpolate as mpi

# 转换为Metpy的数据类
p = df['气压'].values * units.hPa
T = df['温度'].values * units.degC
rh = df['湿度'].values / 100
Td = mpcalc.dewpoint_rh(T, rh)
wsp = df['风速'].values * units('m/s')
wdr = df['风向'].values * units.degrees
u, v = mpcalc.wind_components(wsp, wdr)

# 开一个画布
fig = plt.figure(figsize=(12, 12))
skew = SkewT(fig, rotation=0)

# 画温度、露点和风
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')
skew.plot_barbs(p[::5], u[::5], v[::5])
skew.ax.set_ylim(1000, 250)
skew.ax.set_xlim(-40, 30)

# 计算对流凝结高度，p[0]为最低层气压
lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')

# 画状态曲线
prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')
skew.plot(p, prof, 'k', linewidth=2)

# 画对流有效位能和抑制能
skew.shade_cin(p, T, prof)
skew.shade_cape(p, T, prof)

# 画0度温度线，图中紫色虚线
skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)

# 干绝热线
skew.plot_dry_adiabats()
# 湿绝热线
skew.plot_moist_adiabats()
# 等比湿线
skew.plot_mixing_lines()

T-lnP图, 世界时: 2020-08-02T10

可以看出台风过后CAPE还是比较大。好了，这期就到这了。