Python Metpy教程：计算CAPE指数示例

1.简介

 metpy在气象中很常用，通过输入基本量计算一些比较复杂的参数，比如CAPE/CIN，k指数，抬升凝结高度，相当位温等。但有两点需要注意：

①单位统一成标准单位。

②有些物理量只能用一维数据计算，比如接下来计算的CAPE，如果使用wrfout数据，目前还没有办法将三维数组做整体计算，只能循环取一格点上的所有层作为一维数据输入。

2.安装

进入Install Guide — MetPy 1.6dev (unidata.github.io)

3.实例——以计算wrfout文件中的CAPE为例

3.1读入数据

我在一开始读数据的时候就把DataArray转成ndarray，否则在最后计算的时候会出现“ AttributeError: 'DataArray' object has no attribute 'to' ”的错误；并且只取了其中一个格点上的所有数据，这个要看具体计算的是什么，guide里面会有提示。

import xarray as xr
import numpy as np
import os
import metpy.calc as mpcalc
from metpy.units import units
from wrf import getvar
from netCDF4 import Dataset

##-- read data
path=r'G:/wrfout/wrfout_d02_2021-01-27_20_00_00'
data=Dataset(path+file,'r')
pre = np.array(getvar(data, 'pressure')[:,20,20])
tc = np.array(getvar(data, 'tc')[:,20,20])
td = np.array(getvar(data, 'td')[:,20,20])

3.2 给所有量带上计算的标准单位

pre = pre * units.hPa
tc = tc * units.degC
td = td * units.degC

这里需要注意的是，这个步骤并不会改变数组原来的数值，比如你的温度一开始读的是开尔文单位，这一步不会将你的数值自动减掉273.15并带上℃单位，而是将开尔文的数值直接带上℃单位，导致最后可能会报错。

3.3计算CAPE

#- computure parcel temperature
prof = mpcalc.parcel_profile(pre, tc[0], td[0]).to('degC')

#- calculate surface based CAPE
cape_cin = mpcalc.cape_cin(pre, tc, td, prof)

计算parcel_profile的时候前面把DataArray转成ndarray的作用就来了。

这个函数把CAPE和CIN一起输出，第一个是CAPE，第二个是CIN，并且带有单位，如果只想要CAPE的数值，我在气象家园里面发现有大神解答过这个问题：将Metpy库输出结果转换为不带单位的实数-气象软件开发-气象家园_气象人自己的家园 (06climate.com)

只需要输入print(cape_temp.m.item())

cape_temp = mpcalc.cape_cin(pre_temp, tc_temp, td_temp, prof_temp)[0]
print(cape_temp.m.item())

但我的报错了，因为我在实际计算中写了循环以计算所有网格点上的CAPE，有些格点上的值是0，于是就会报错“ 'int' object has no attribute 'item' ”，所以需要转成“float”就不会报错了

cape_temp.astype('float').m.item()

4.计算其他物理量

在Referrence Guide里的Calculations里有很多物理量的例子，写得比较简洁明了。此外，metpy还能画图，我用metpy画过的T-lnP图，还蛮好用的。

作者：唐宇直Tangyz