如何使用 metpy 正确计算温度平流，单位错误

Question

我对 Metpy 有点陌生。我一直在尝试用 Metpy 计算温度平流，但没有成功。由于我是这个包的新手，我不明白为什么需要有单位才能正常工作。当我计算温度平流时，我的地图上出现了一些奇怪的线条，但我不知道为什么。我认为这是因为单位或其他原因，但我不确定。我在下面附上我的脚本：

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.gridspec as gridspec
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import scipy.ndimage as ndimage
import metpy.calc as mpcalc
from metpy.units import units
import matplotlib.patches as mpatches
from metpy.units import units
from netCDF4 import Dataset
from netCDF4 import num2date
from siphon.catalog import TDSCatalog
from pint import UnitRegistry

crs = ccrs.PlateCarree() #ccrs.Mercator()
def plot_background(ax):
    ax.set_extent([-80., -15., -10., -60.],ccrs.PlateCarree())
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)
    ax.add_feature(cfeature.BORDERS, linewidth=0.5)
    gl = ax.gridlines(ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.5, linestyle='--')
    gl.xlabels_top = False
    gl.xlabels_bottom = True
    gl.ylabels_left = True
    gl.ylabels_right = False
    gl.xlines = False; gl.ylines= False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabel_style = {'size': 12, 'color': 'black'}
    gl.ylabel_style = {'size': 12, 'color': 'black'}
    return ax


def define_box_region(lon1,lon2,lat1,lat2,):
    lat_corners = np.array([lat1,  lat1, lat2, lat2]); 
    lon_corners = np.array([ lon1, lon2, lon2, lon1]) 
    poly_corners = np.zeros((len(lat_corners), 2))
    poly_corners[:,0] = lon_corners; poly_corners[:,1] = lat_corners
    poly = mpatches.Polygon(poly_corners, closed=True, ec='m', fill=False, lw=1, fc="yellow", transform=ccrs.PlateCarree())
    return poly
infilesRegEraDJF=list(open('DJFRegEra.txt', 'r'))


print(infilesRegEraDJF[1][:-1])
d1=infilesRegEraDJF[1][:-1]
data=Dataset(d1)
tahist1=data.variables['ta850'][:,:]
tahist1 = np.squeeze(tahist1, axis=None) #gsp =ERA
tahist1=np.transpose(tahist1)
# tahist1=tahist1-273.15
temp850=tahist1 *units.degK
# t850=tahist1.to('degC')
# t850=tahist1*units.degC
lats_data= data.variables['lat'][:]
lats=np.transpose(lats_data)
lons_data= data.variables['lon'][:]
lons=np.transpose(lons_data)


print(infilesRegEraDJF[2][:-1])
d1=infilesRegEraDJF[2][:-1]
data=Dataset(d1)
uahist1=data.variables['ua850'][:,:]
uahist1 = np.squeeze(uahist1, axis=None) #gsp =ERA
uahist1=np.transpose(uahist1)
u850=uahist1*units('m/s')

print(infilesRegEraDJF[3][:-1])
d1=infilesRegEraDJF[3][:-1]
data=Dataset(d1)
vahist1=data.variables['va850'][:,:]
vahist1 = np.squeeze(vahist1, axis=None) #gsp =ERA
vahist1=np.transpose(vahist1)
v850=vahist1*units('m/s')

dx1, dy1 = mpcalc.lat_lon_grid_deltas(lons, lats)

advRegEra = mpcalc.advection(temp850, [u850, v850],
                       (dx1, dy1), dim_order='yx')

advregera2=advRegEra.to(units('delta_degC/sec'))


values_adv=[-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5]


fig, axarr = plt.subplots(nrows=3, ncols=2, figsize=(9.5, 11), constrained_layout=True,
                              subplot_kw={'projection': crs})
axlist = axarr.flatten()
for ax in axlist:
 plot_background(ax)

cf1=axlist[1].contourf(lons, lats,advRegEra,values_adv,cmap='RdBu_r',extend='both',transform=ccrs.PlateCarree())
cf3 = axlist[1].quiver(lons[::17,::17], lats[::17,::17],uahist1[::17,::17],vahist1[::17,::17],scale=200,headwidth=5, headlength=5,transform=ccrs.PlateCarree())
axlist[1].set_title('RegERA', fontsize=14)
axlist[1].quiverkey(cf3,X=0.9, Y=1.05,U=10,label='10 m/s', labelpos='E')
poly=define_box_region(-52,-38,-35,-22.5); axlist[1].add_patch(poly)  #Box SBR
poly=define_box_region(-67.5,-52,-37.5,-25); axlist[1].add_patch(poly)  #Box LPB
poly=define_box_region(-72.5,-57.5,-37.5,-55); axlist[1].add_patch(poly)  #Box ARG


cax=fig.add_axes([0.95,0.30,0.025,0.40])
ax3=fig.colorbar(cf1,ticks=values_adv,cax=cax,orientation='vertical',label='(°C/hr)')
cax.tick_params(labelsize=14)
fig.set_constrained_layout_pads(w_pad=0.1, h_pad=0.1, hspace=0.1, wspace=0.1)

当我运行将 adv（单位应为 K/s）转换为 advregera2=advRegEra.to(units('delta_degC/sec')) 时，我得到以下错误

DimensionalityError: Cannot convert from '1 / meter' (1 / [length]) to 'delta_degree_Celsius / second' ([temperature] / [time])

我也尝试将单位从头开始放在左侧 (tahist1=units.K*tahis1)，这可行，但是当我绘制地图时，我得到一些奇怪的线条。

有人可以帮助我吗？这是我第一次使用 metpy，所以我很难理解平流温度的功能是如何工作的:(。我在 Linux OpenSUSE 15.1 leap

上使用 spyder 3.7

Answer 1

advection 绝对是在 MetPy 中使用的比较棘手的函数之一。由于您正在使用 netcdf4-python 打开文件，因此您肯定想乘以左侧的单位，例如：

u850 = units('m/s') * uahist1

这应该可以解决您的设备问题。关于奇怪的线条，我猜它们看起来像这个 similar GitHub issue 中的图像？如果是这样，那是由于您的情节从 +180 到 -180 经度环绕，但数据从 0 开始到 360 度经度结束（至少我是这么猜的）引起的。假设这就是问题所在，因为您实际上不需要穿过这些切口，所以您应该能够通过仅在感兴趣区域周围绘制数据的子切片来删除这些点。类似于：

lon_subset = slice(200, 300)
cf1=axlist[1].contourf(lons[lon_subset], lats, advRegEra[:, lon_subset],
    values_adv, cmap='RdBu_r',extend='both',transform=ccrs.PlateCarree())

用正确的索引替换上面的 200 和 300。

如何使用 metpy 正确计算温度平流，单位错误

How to properly calculate temperature advection with metpy, error with units

python

temperature

units-of-measurement

spyder

metpy