如何用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: 无法从“1 / 米”（1 / [长度]）转换为“delta_degree_Celsius / 秒”（[温度] / [时间]）

我还尝试将单位从一开始放在左侧 (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 替换为正确的索引。