我正在尝试使用metpy的cape_cin和get_layer函数计算指定层中的cape和cin。我这样做是为了从NCDC服务器访问RAP垂直配置文件，用于3种不同的地块类型(ML、MU和SB(。当我尝试计算mu和sb配置文件时，我得到了这个错误ValueError: zero-size array to reduction operation minimum which has no identity，但对于ml配置文件没有。。。即使我以完全相同的方式计算图层：
首先，我必须从RAP垂直剖面创建压力、温度、露点和高度的垂直剖面：
from datetime import datetime, timedelta
from siphon.catalog import TDSCatalog
from siphon.ncss import NCSS
import numpy as np
from metpy.units import units
import cartopy.crs as ccrs
from metpy.calc import (dewpoint_from_relative_humidity, mixed_parcel, most_unstable_parcel, parcel_profile, pressure_to_height_std, lcl, height_to_pressure_std, get_layer, cape_cin)
year=2019
month=5
day=23
hour=0
cenlat = 34.91269
cenlon = -98.21048
time_start = datetime(year, month, day, hour, 0) #specified time
hour = time_start.hour
if hour < 10:
hour = '0'+str(hour)
day = time_start.day
if day < 10:
day = '0'+str(day)
month = time_start.month
if month < 10:
month = '0'+str(month)
cat = TDSCatalog('https://www.ncdc.noaa.gov/thredds/catalog/model-rap130-old/'+str(time_start.year)+str(month)+'/'+str(time_start.year)+str(month)+str(day)+'/catalog.html?dataset=rap130-old/'+str(time_start.year)+str(month)+'/'+str(time_start.year)+str(month)+str(day)+'/rap_130_'+str(time_start.year)+str(month)+str(day)+'_'+str(hour)+'00_000.grb2')
latest_ds = list(cat.datasets.values())[0]
print(latest_ds.access_urls)
ncss = NCSS(latest_ds.access_urls['NetcdfSubset'])
query = ncss.query()
query.variables('Pressure_surface').variables('Geopotential_height_isobaric').variables('Geopotential_height_surface').variables('Relative_humidity_isobaric').variables('Temperature_isobaric').variables('Dewpoint_temperature_height_above_ground').variables('Temperature_height_above_ground').variables
query.add_lonlat().lonlat_box(cenlon-2.1, cenlon +2.1, cenlat-2.1, cenlat+2.1)
data1 = ncss.get_data(query)
dlev = data1.variables['Geopotential_height_isobaric'].dimensions[1]
dlat = data1.variables['Geopotential_height_isobaric'].dimensions[2]
dlon = data1.variables['Geopotential_height_isobaric'].dimensions[3]
SFCP = (np.asarray(data1.variables['Pressure_surface'][:])/100.) * units('hPa')
hgt = np.asarray(data1.variables['Geopotential_height_isobaric'][:]) * units('meter')
sfc_hgt = np.asarray(data1.variables['Geopotential_height_surface'][:]) * units('meter')
Temp_up = np.asarray(data1.variables['Temperature_isobaric'][:]) * units('kelvin')
RH_up = np.asarray(data1.variables['Relative_humidity_isobaric'][:])
Td = (np.asarray(data1.variables['Dewpoint_temperature_height_above_ground'][:]) * units('kelvin')).to('degC')
T = np.asarray(data1.variables['Temperature_height_above_ground'][:]) * units('kelvin')
# Get the dimension data
lats_r = data1.variables[dlat][:]
lons_r= data1.variables[dlon][:]
lev = (np.asarray(data1.variables[dlev][:])/100.) * units('hPa')
# Set up our array of latitude and longitude values and transform to the desired projection.
flon = float(cenlon)
flat = float(cenlat)
crs = ccrs.PlateCarree()
crlons, crlats = np.meshgrid(lons_r[:]*1000, lats_r[:]*1000)
trlatlons = crs.transform_points(ccrs.LambertConformal(central_longitude=265, central_latitude=25, standard_parallels=(25.,25.)),crlons,crlats)
trlons = trlatlons[:,:,0]
trlats = trlatlons[:,:,1]
dlon = np.abs(trlons - cenlon)
dlat = np.abs(trlats - cenlat)
ilon = np.where(dlon == np.min(dlon)) #position in the dlon array with minimal difference between gridpoint lon and input lon
ilat = np.where(dlat == np.min(dlat)) #position in the dlat array with minimal difference between gridpoint lat and input lat
Tdc_up = dewpoint_from_relative_humidity(Temp_up[0,:,ilat[0][0], ilon[1][0]],RH_up[0,:,ilat[0][0], ilon[1][0]]/100)
p_sounding = np.sort(np.append(lev, SFCP[0,ilat[0][0], ilon[1][0]]))
ind = np.where(p_sounding >= SFCP[0,ilat[0][0], ilon[1][0]])[0][0]
hgt_sounding = np.insert(hgt[0,:,ilat[0][0], ilon[1][0]].magnitude, ind, sfc_hgt[0,ilat[0][0], ilon[1][0]].magnitude) * hgt.units
T_sounding = (np.insert(Temp_up[0,:,ilat[0][0], ilon[1][0]].magnitude, ind, T[0,0,ilat[0][0], ilon[1][0]].magnitude) * T.units).to(Tdc_up.units)
Td_sounding = np.insert(Tdc_up.magnitude, ind, Td[0,0,ilat[0][0], ilon[1][0]].magnitude) * Tdc_up.units
p_skewt = p_sounding[p_sounding <= SFCP[0,ilat[0][0], ilon[1][0]]]
hgt_skewt = hgt_sounding[p_sounding <= SFCP[0,ilat[0][0], ilon[1][0]]]
T_skewt = T_sounding[p_sounding <= SFCP[0,ilat[0][0], ilon[1][0]]]
Td_skewt = Td_sounding[p_sounding <= SFCP[0,ilat[0][0], ilon[1][0]]]
AGLhgts = hgt_skewt[::-1]-hgt_skewt[-1]
接下来，我使用mixed_package、most_unstatible和parcel_profile函数为每个地块类型创建垂直剖面，并计算我要计算cape_cin的层的顶部和底部的压力值(LCL到LCL+2km(：
ml_p, ml_T, ml_Td = mixed_parcel(np.flip(p_skewt), np.flip(T_skewt), np.flip(Td_skewt))
ml_profile = parcel_profile(p_skewt[::-1], ml_T, ml_Td)
ml_profile = (ml_profile - 273.15*units('kelvin')).magnitude*units('degC')
mu_p, mu_T, mu_Td, mu_index = most_unstable_parcel(np.flip(p_skewt), np.flip(T_skewt), np.flip(Td_skewt))
mu_profile = parcel_profile(p_skewt[::-1], mu_T, mu_Td)
mu_profile = (mu_profile - 273.15*units('kelvin')).magnitude*units('degC')
#Note: sbpcl_profile will have the exact same values of p_skewt, T_skewt, and Td_skewt in pprof below:
pprof = parcel_profile(p_skewt[::-1], T_skewt[-1], Td_skewt[-1])
pprof = (pprof - 273.15*units('kelvin')).magnitude*units('degC')
mllcl = lcl(ml_p, ml_T, ml_Td)
mllcl_h = pressure_to_height_std(mllcl[0]) - hgt_skewt[-1]
mulcl = lcl(mu_p, mu_T, mu_Td)
mulcl_h = pressure_to_height_std(mulcl[0]) - hgt_skewt[-1]
sblcl = lcl(p_skewt[-1], T_skewt[-1], Td_skewt[-1])
sblcl_h = pressure_to_height_std(sblcl[0]) - hgt_skewt[-1]
mllcl2000 = mllcl_h + 2*units('kilometer')
mulcl2000 = mulcl_h + 2*units('kilometer')
sblcl2000 = sblcl_h + 2*units('kilometer')
mllcl2000_p = height_to_pressure_std(mllcl2000)
mulcl2000_p = height_to_pressure_std(mulcl2000)
sblcl2000_p = height_to_pressure_std(sblcl2000)
计算完所有这些之后，我使用get_layer函数创建压力、温度、露点和包裹温度的数组，我需要计算cape_cin，然后计算感兴趣层中的实际cape_cin值：
ml_LCL_CAPE_layer = get_layer(p_skewt, T_skewt, Td_skewt, ml_profile[::-1], bottom = mllcl[0], depth = mllcl[0] - mllcl2000_p)
mu_LCL_CAPE_layer = get_layer(p_skewt, T_skewt, Td_skewt, mu_profile[::-1], bottom = mulcl[0], depth = mulcl[0] - mulcl2000_p)
sb_LCL_CAPE_layer = get_layer(p_skewt, T_skewt, Td_skewt, pprof[::-1], bottom = sblcl[0], depth = sblcl[0] - sblcl2000_p)
mlLCLCAPE = cape_cin(ml_LCL_CAPE_layer[0], ml_LCL_CAPE_layer[1], ml_LCL_CAPE_layer[2], ml_LCL_CAPE_layer[3])
muLCLCAPE = cape_cin(mu_LCL_CAPE_layer[0], mu_LCL_CAPE_layer[1], mu_LCL_CAPE_layer[2], mu_LCL_CAPE_layer[3])
sbLCLCAPE = cape_cin(sb_LCL_CAPE_layer[0], sb_LCL_CAPE_layer[1], sb_LCL_CAPE_layer[2], sb_LCL_CAPE_layer[3])
mlLCLCAPEcin = mlLCLCAPE[0] + mlLCLCAPE[1]
muLCLCAPEcin = muLCLCAPE[0] + muLCLCAPE[1]
sbLCLCAPEcin = sbLCLCAPE[0] + sbLCLCAPE[1]
3个get_layer函数中每个函数的压力、温度、露点和地块温度的数组似乎都填充了正确的值，并且每个地块类型的这4个数组都是相同的形状。上面的mlLCLCAPEcin计算给出了正确的输出(99.26 j/kg-当我将其绘制在SkewT上时进行验证(，但MU和SB剖面的完全相同的计算给出了上面提到的错误。我使用的是Metpy v 1.1，并尝试使用不同的位置和不同预测小时的输出，但仍然遇到相同的问题。
如果我用修复了上面的示例代码(存在一些名称问题和一些索引问题(
T_sounding = (np.insert(Temp_up[0,:,ilat[0][0], ilon[1][0]].magnitude, ind, Temp_up[0,0,ilat[0][0], ilon[1][0]].magnitude) * Temp_up.units).to(Tdc_up.units)
Td_sounding = np.insert(Tdc_up.magnitude, ind, Tdc_up[0].magnitude) * Tdc_up.units
我现在运行时没有任何错误。如果我上面的修复程序不能为您修复它，那么了解您是否正在运行最新的MetPy 1.1会很有帮助。