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

Question

我对梅蒂有点陌生。我一直试图用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上使用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