使用metpy的skewt画气象探空图并计算CAPE_CIN_metpy官网_fangzuliang的博客

1 metpy介绍

metpy是面向地球科学的较常用的库，可视化能力非常强，能基本满足气象学子的绘图要求。官网上给了许多示例代码，参考他们的代码进行适当的改动就可以私人定制出非常漂亮实用的图。
安装介绍
官网上给出了安装的环境需求和安装方法，非常方便。
#pip管理
pip install metpy
#conda管理
conda install -c conda-forge metpy
2 官网代码介绍
 
以metpy探空代码为例，进行介绍 
=================
Advanced Sounding
=================
Plot a sounding using MetPy with more advanced features.
Beyond just plotting data, this uses calculations from `metpy.calc` to find the lifted
condensation level (LCL) and the profile of a surface-based parcel. The area between the
ambient profile and the parcel profile is colored as well.
import matplotlib.pyplot as plt
import pandas as pd
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.plots import add_metpy_logo, SkewT
from metpy.units import units
###########################################
# Upper air data can be obtained using the siphon package, but for this example we will use  即可以使用siphon包获取探空数据；
# some of MetPy's sample data.
##step1: 使用官网的示例数据,并做一些预处理
col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
df = pd.read_fwf(get_test_data('may4_sounding.txt', as_file_obj=False),
                 skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
# Drop any rows with all NaN values for T, Td, winds
df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'), how='all'
               ).reset_index(drop=True)
###########################################
# We will pull the data out of the example dataset into individual variables and
# assign units.
##step2:这一步很重要，将array数据转换为 Unit类型的数据；不转就会出错。
##获取温压湿风信息
p = df['pressure'].values * units.hPa
T = df['temperature'].values * units.degC
Td = df['dewpoint'].values * units.degC
wind_speed = df['speed'].values * units.knots
wind_dir = df['direction'].values * units.degrees
u, v = mpcalc.wind_components(wind_speed, wind_dir)
###########################################
# Create a new figure. The dimensions here give a good aspect ratio.
#创建图形，画探空图，rotation表示探空图中等温线与水平面的夹角。
fig = plt.figure(figsize=(9, 9))
# add_metpy_logo(fig, 115, 100)  #这一行在图上添加metpy官方图标
skew = SkewT(fig, rotation = 45)
# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot.
skew.plot(p, T, 'r')  #环境温度垂直廓形
skew.plot(p, Td, 'g') #环境露点垂直廓线
skew.plot_barbs(p, u, v)  #画出水平风场
skew.ax.set_ylim(1000, 100) #设置纵轴(气压)范围
skew.ax.set_xlim(-40, 60)   #设置横轴(温度)范围
# Calculate LCL height and plot as black dot. Because `p`'s first value is
# ~1000 mb and its last value is ~250 mb, the `0` index is selected for
# `p`, `T`, and `Td` to lift the parcel from the surface. If `p` was inverted,
# i.e. start from low value, 250 mb, to a high value, 1000 mb, the `-1` index
# should be selected.
#基于地面气压、温度和露点(以1000hPa充当地面)，计算出抬升凝结高度LCL对应的气压层和温度,并在图上标出
lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
# Calculate full parcel profile and add to plot as black line
#画出气块绝热路径
prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')
skew.plot(p, prof, 'k', linewidth=2)
# Shade areas of CAPE and CIN
#填充CAPE和CIN区域
skew.shade_cin(p, T, prof)
skew.shade_cape(p, T, prof)
# An example of a slanted line at constant T -- in this case the 0
# isotherm
#画出某一条温度 = k的等温线
skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
# Add the relevant special lines
#画出干绝热线、湿绝热线和 混合比线
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
# Show the plot
plt.show()
出图效果如下： 
 来看看数据类型
 print('type of p:',type(p))
print()
print('p:',p)
type of p: <class 'pint.quantity.build_quantity_class.<locals>.Quantity'>
p: [959.0 931.3 925.0 899.3 892.0 867.9 850.0 814.0 807.9 790.0 779.2 751.3 724.3 700.0 655.0 647.5 599.4 554.7 550.0 500.0 472.5 449.0 400.0 383.7 336.4 321.9 308.1 300.0 269.0 268.6] hectopascal
 我也是首次接触这种类型的数据，如果有了解的，麻烦给与指导！
 而且要注意：p[0]应该最低，其他温度湿度数据也要与气压相对应，低层值在前。
 
 除此之外，matplotlib官网也提供了skew图的绘制代码，但比较粗糙，信息不够全面
  
3 私人定制
 
参考官网代码，来一个私人定制 
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import netCDF4 as nc
import os
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.plots import add_metpy_logo, SkewT
from metpy.units import units
#使得中文字体不乱码
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
#matplotlib.rcParams['font.family']='sans-serif'
plt.rcParams['axes.unicode_minus']=False
os.chdir('D:/ERA5')
deta = 4
CI_file = 'D:/ERA5/CI_hecheng/' + str(deta) + '_smooth_pressure.nc'
no_CI_file =




    
 'D:/ERA5/no_CI_hecheng/' + str(deta) + '_smooth_pressure.nc'
gap = 0.25  
#该范围包含多少个经纬度格点
grid_num = int(deta/gap)*2+1
def plot_skewt_with_mepy(file, 
                        time_deta = 0,
                        top_level = 20,
                        point_x_y = [20,20],
                        fontsize = 20,
    func: 输入.nc文件，获取里面的p t td u v
    inputs:
        file:文件绝对路径 + 文件名称，eg: D:/ERA5/CI_hecheng/4_smooth_pressure.nc
        time_deta: 默认0, 即对应14:00
        top_level: 探空图中顶层气压层，默认为20hPa
        point_x_y: 选取哪个位置画出探空， 默认[20,20] 
        fontsize: 字体大小，默认20
    return:
    Vars = ['t','q',
            'u',
            'v',
    Vars_operation = ['-273.16','*1e3',
                 '*1',
                 '*1',
                 '-273.16'
    #默认的以14:00为中间时间，对应的time_index = 30
    time_center = 14
    #确定探空点所在位置
    point_x = point_x_y[0]
    point_y = point_x_y[1]
    fr = nc.Dataset(file)
    #顶层气压在气压层中的 index, 和需要的气压
    top_level_index = list(fr['level'][::-1]).index(top_level)
    p = np.array(fr['level'][::-1][0:top_level_index + 1])* units.hPa
    print('p.shape:',p.shape)
    all_data = []
    for Var, oper in zip(Vars,Vars_operation):
        #确定本地时间, index=30表示北京时14时
        current_time_index = 30 + time_deta
        # local_beijing_time = time[current_time_index]+8
        local_beijing_time = time_center + time_deta
        print('local_beijing_time:',local_beijing_time)
        #数据是从最高层开始排序的，因此，需要进行颠倒过来
        print(Var)
        Var_data = fr[Var][:][current_time_index][-1::-1,-1::-1,:]
        Var_data = np.array(Var_data)[0:top_level_index + 1, point_x,point_y]
        print(Var_data.shape)
        #对数据进行处理
        Var_data = eval('Var_data' + oper)
        all_data.append(Var_data)
    fr.close()
    #按照官网例子进行数据处理和类型转换
    T = all_data[Vars.index('t')]* units.degC
    Td = all_data[Vars.index('Td')] * units.degC
    u = all_data[Vars.index('u')] * units.knots
    v = all_data[Vars.index('v')] * units.knots
    #如果不进行数据类型转换，就会报错
#    T = all_data[Vars.index('t')]
#    Td = all_data[Vars.index('Td')]
#    u = all_data[Vars.index('u')] 
#    v = all_data[Vars.index('v')] 
    #创建一个子图和句柄 ax
    fig1, ax = plt.subplots(1,1, figsize=(9, 9))
    skew = SkewT(fig1, rotation = 30)    
    ##画出温度，露点廓线和水平uv风场
    skew.plot(p, T, 'r')
    skew.plot(p, Td, 'g')
    skew.plot_barbs(p, u, v)
    # Calculate LCL height and plot as black dot. Because `p`'s first value is
    # ~1000 mb and its last value is ~250 mb, the `0` index is selected for
    # `p`, `T`, and `Td` to lift the parcel from the surface. If `p` was inverted,
    # i.e. start from low value, 250 mb, to a high value, 1000 mb, the `-1` index
    # should be selected.
    #基于地面气压、温度和露点(以1000hPa充当地面)，计算出抬升凝结高度LCL,并在图上标出
    lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
    skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
    #计算出气块绝热路径，用黑色线条表示
    # Calculate full parcel profile and add to plot as black line
    prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')
    skew.plot(p, prof, 'k', linewidth=2)
    #计算出CAPE 和 CIN面积
    # Shade areas of CAPE and CIN
    skew.shade_cin(p, T, prof)
    skew.shade_cape(p, T, prof)
	#计算cape和cin值
    cape,cin = mpcalc.cape_cin(p,T,Td, prof,
                               which_el = 'most_cape',
    print('cape:',cape)
    print('cin:',cin)
    # An example of a slanted line at constant T -- in this case the 0
    # isotherm
    #将 0 温度等值线画出来 
    skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
    # Add the relevant special lines
    #画出干绝热线、湿绝热线和 混合比线
    skew.plot_dry_adiabats()
    skew.plot_moist_adiabats()
    skew.plot_mixing_lines()
    fontdict = {'fontsize' : fontsize}
    #设置横纵坐标数值范围
    x_left = -40
    x_right = 60
    skew.ax.set_ylim(1000, 100)
    skew.ax.set_xlim(x_left, x_right)
    #设置对应的xtickslabel 和 ytickslabel
    #对应的句柄是 skew.ax
    y_ticks = np.arange(1000,0,-100)
    y_tickslabel = np.arange(1000,0,-100)
    x_ticks = np.arange(x_left,x_right + 10,10)
    x_tickslabel = np.arange(x_left,x_right + 10,10)
    skew.ax.set_yticks(y_ticks)
    skew.ax.set_yticklabels(y_tickslabel, fontdict = fontdict)
    skew.ax.set_xticks(x_ticks)
    skew.ax.set_xticklabels(x_tickslabel, fontdict = fontdict)
    #设置横纵label
    skew.ax.set_ylabel('Height/hPa',fontsize = fontsize)
    skew.ax.set_xlabel('T/(℃)',fontsize = fontsize)
    #由于是在子图上画的skew的，因此ax的横纵坐标轴上，
    #会自动出现一些数字，可以使用ax.set_ticklabels功能使其不显示
    ax.set_yticklabels([])
    ax.set_xticklabels([])
    # Show the plot
    plt.show()    
    #可以查看设置了tickslabel后的ticks数值
	print('ax 的ticks：',ax.get_yticklabels()[:])
    print()
    print('skewt 的ticks：',skew.ax.get_yticklabels()[:])
plot_skewt_with_mepy(CI_file, 
                    time_deta = 3,
                    top_level = 200,
                    point_x_y = [20,25],
outputs:
cape: 673.6904340493442 joule / kilogram
cin: -231.31875105377705 joule / kilogram
ax 的ticks： [Text(0, 0.0, ''), Text(0, 0.2, ''), Text(0, 0.4, ''), Text(0, 0.6000000000000001, ''), Text(0, 0.8, ''), Text(0, 1.0, '')]
skewt 的ticks： [Text(0, 1000, '1000'), Text(0, 900, '900'), Text(0, 800, '800'), Text(0, 700, '700'), Text(0, 600, '600'), Text(0, 500, '500'), Text(0, 400, '400'), Text(0, 300, '300'), Text(0, 200, '200'), Text(0, 100, '100')]
结果如图： 
得到了信息很丰富的漂亮的探空图 
4 参考链接
 
metpy官网 
Metpy新版功能下载TLnP图设置 - Bugatii100Peagle 
Siphon官网 
气象上常用的一些软件和库的docker镜像 
matplotlib绘制skewt官方例子 
metpy计算cape_cin
                                    注意plt.函数名()相当于面向过程的画图方法，axes.set_方法名()相当于面向对象的画图方法。概念用两组数据构成多个坐标点，考察坐标点的分布,判断两变量之间是否存在某种关联或总结坐标点的分布模式。特点绘制连离散的数据,能够一眼看出各个数据的大小,比较数据之间的差别。可视化是在整个数据挖掘的关键辅助工具，可以清晰的理解数据，从而调整我们的分析方法。特点绘制连续性的数据展示一组或者多组数据的分布状况(统计)特点能够显示数据的变化趋势，反映事物的变化情况。特点分类数据的占比情况(占比).........
                                    SkewT图简介
又名SkewT-LogP图，是一个非常有用的气象热力学图，在该图上绘制了整个大气层上地球表面某个点的压力，密度，温度和水蒸气。每个SkewT图都包含了一个静态的大气计算（适用于所有的位置），以及一组动态的空气温度和露点温度计算，这些计算是通过直接测量大气获得的。通常通过每天从全球很多机场和气象站放出一个小的探空仪器（一般是气象气球）来获得这组观测值。美国国家航空航天局（NASA）的大气红外探测仪（AIRS）和高级微波探测器(AMSU-A)
露点温度：在空气中水汽含量不变，保持气压一定的情
                                    我们搬运翻译了Metpy官网的绘图例程，并且调试跑通了，放在了和鲸社区上。
1.SkewT
代码：https://www.heywhale.com/mw/project/5f21358cd278b1002c217072
本示例使用了怀俄明大学测绘档案中的示例数据，时间为2016年10月31日，明尼苏达州明尼阿波利斯（MPX），并使用MetPy绘制了带有温度，露点和风羽的经典偏斜T。
2.Miller Composite图表
代码：https://www.heywhale.com/mw/project/5f
                                    MetPy是Python中的一系列工具，用于读取、可视化天气数据并进行计算
MetPy MetPy 是 Python 中的一组工具，用于读取、可视化和执行天气数据计算。
 MetPy 在其版本号中遵循语义版本控制。
 对于我们当前的 0.x 版本，这意味着 MetPy 的 API（应用程序编程接口）仍在不断发展（我们不会为了好玩而破坏事物，但在我们解决设计问题时，许多事情仍在发生变化）。
 此外，对于 0.xy 版本，我们在发布新功能时更改 x，在发布仅修复错误的版本时更改 y。
 我们支持 Python >= 3.6，目前支持 Python 2.7。
 注意：我们将在 2019 年秋季停止对 Python 2.7 的支持。有关更多信息，请参见此处。
 需要帮忙？
 在使用 MetPy 时需要帮助？
 发现问题了吗？
 有功能请求？
 查看我们的支持页面。
 重要链接 源代码库：https://github.com/Unidata/MetPy HTML 文档：http://unidata.github.io/MetPy Unidata Python 库：https://unida
                                    MetPy MetPy是Python中的一组工具，用于读取，可视化和执行天气数据的计算。
MetPy在其版本号中遵循语义版本控制。
使用我们当前的0.x MetPy MetPy是Python中的一组工具，用于读取，可视化和执行天气数据的计算。
MetPy在其版本号中遵循语义版本控制。
在我们当前的0.x版本中，这意味着MetPy的API（应用程序编程接口）仍在不断发展（我们不会仅仅为了好玩而破坏事物，但是在解决设计问题时，许多事物仍在发生变化）。
另外，对于版本0.xy，我们在发布新功能时会更改x，而在发布新功能时会更改y
                                    Draw_Prcl    = 10    ;*不画抬升曲线/0：只画抬升曲线，不计算LFC等/1：只画线+计算/2：填色+画线+计算
   Draw_dewp    = 10    ;*是否画湿度Td探空曲线0/1
   DrawPMax     = 10    ;*是否画从最不稳定层抬升曲线0/1
   DrawIndx     = 10    ;*是否画指标0/1
   Draw_RH      = 10    ;*是否画相对湿度RH探空曲线0/1
   DrawBarb     = 10    ;*是否画风向标0/1
   DrawHodo     = 10    ;*是否画风矢端图0/1
   DrawHeli     = 10    ;*是否画螺旋度0/1
   Draw_CAP     = 10    ;*是否计算CAP参数0/1
* _NorightAll  =  0     ;*不画右侧所有的框（1：不画/0:画）
                                    最近有任务需要计算探空指数产品，在经历了漫长扒文献、找公式实现公式的过程后，终于是凑齐了所需的指数产品。总结一下这一段的历程，同时提供给可能需要的人，提高效率减少大海中找针的工作量。
                                    基于metpy库画T-LnP图（国内主流形式，并非斜温）
时隔多日再次更帖，主要原因是一直在画T-LnP，想画一个好的，这一版是为接下来的下一帖给小白们打一下基础。
本帖的知识点为：
1.利用metpy库的SkewT类来画T-LnP图
2.改斜温图为国内传统的图像
导入相关库
import pandas as pd#数据读取用
import metpy.calc as mpcalc#计算一些参数用，比如 抬升凝结高度
from metpy.plots import  SkewT#画埃玛图
from metp
                                    PR 1182 - Update Project Layout (1171)PR 1181 - Rearrange CartoPy Imports (1159)PR 1180 - SkewT Updates (60, 296, 1176)PR 1179 - Update PyPI uploadPR 1178 - Deprecate the viridis colortable (1136)PR 1...
                                    条件性不稳定是一种潜在性不稳定。只有当起始 高度上有较强的抬升力或冲击力，将气块抬升到 自由对流高度(LFC) 以上时，对流运动才能发展，不稳定才表现出来。850hPa上的湿空气沿着干绝热线抬升，到达抬升凝结高度（LCL）后再沿着湿绝热线上升至500hPa具有的气块温度与500hPa环境温度之差。对于海拔较高的地区，不能直接用SI指数，而是应该在低层、高层分别取两个等压面代替850hPa和500hPa。
如果850hPa和500hPa之间存在锋面或者逆温层，SI无意义。
阈值必须针对不同季节、根据各地实际情
                                    登录后查看更多精彩内容~您需要 登录 才可以下载或查看，没有帐号？立即注册x本帖最后由 edwardli 于 2017-6-13 10:02 编辑工作繁事多，先上结果供交流，回头不断细化。我的提纲：气象绘图→科学绘图→Python可视化over（这个过程肯定也会把数据处理搞定）python客户端程序→爬虫→深度学习Python气象绘图研究目录数据使用上遵循“实况的用MICAPS（高空、地面、探空）...