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

问题描述

我对 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,lat2]); 
    lon_corners = np.array([ lon1,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.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,uahist1[::17,vahist1[::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,-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/second”（[温度]/[时间]）

我还尝试将单位从一开始就放在左侧 (tahist1=units.K*tahis1) 并且这有效，但是当我绘制地图时，我得到了一些奇怪的线条。

有人可以帮助我吗？这是我第一次使用metpy，所以我很难理解平流温度函数的工作原理:(。我在Linux OpenSUSE 15.1leap上使用spyder 3.7

解决方法

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。