使用numpy和metpy的湿通量发散不同

Question

如this question 中所述，我想计算 850 hPa 处的水分通量散度。为此，我使用了here 描述的代码，它利用了numpy 包中的np.gradient 函数。我使用的代码是这样的：

from matplotlib import pyplot
from matplotlib.cm import get_cmap
from __future__ import print_function
from netCDF4 import Dataset,num2date,date2num
from matplotlib.colors import from_levels_and_colors
from cartopy import crs
from cartopy.feature import NaturalEarthFeature, COLORS
#
import metpy.calc as mpcalc
import xarray as xr
import cartopy.crs as ccrs
import matplotlib
import cartopy.feature as cfe
import numpy as np
import matplotlib.pyplot as plt
import datetime
#
##################################################################################################################
########################################## Calculate Mosture Divergence ########################################## 
##################################################################################################################
#
root_dir = '/users/pr007/mkaryp/'
nc = Dataset(root_dir+'ERA5.nc')
#
v = np.array(nc.variables['v'][0,:,:])                    
u = np.array(nc.variables['u'][0,:,:])
q = np.array(nc.variables['q'][0,:,:])        
#
lon=nc.variables['longitude'][:]
lat=nc.variables['latitude'][:]
#
qu = q*u
qv = q*v  
#
[dqu_dx, dqu_dy] = np.gradient(qu)
[dqv_dx, dqv_dy] = np.gradient(qv)
#
uq = [dqu_dx, dqu_dy]
vq = [dqv_dx, dqv_dy]
#
divg = np.sum([uq, vq], axis=0)
HMD = divg
# 
##################################################################################################################
#################################################### Plot Map #################################################### 
##################################################################################################################
#
min_val = HMD[0].min()
max_val = HMD[0].max()
diff = max_val-min_val
step = diff/14
#
# Set the figure size, projection, and extent
fig = plt.figure(figsize=(4.5,4))
ax = plt.axes(projection=ccrs.Robinson())
#ax.set_global()
ax.coastlines(resolution="110m",linewidth=1)
ax.gridlines(linestyle='--',color='black')
#
# Set contour levels, then draw the plot and a colorbar
clevs = np.arange(min_val,max_val,step)
plt.contourf(lon, lat, HMD[0], clevs, transform=ccrs.PlateCarree(), cmap=get_cmap("seismic"), extend="both")
plt.title('HMD from ERA5 using numpy 0', size=14)
cb = plt.colorbar(ax=ax, orientation="vertical", pad=0.02, aspect=16, shrink=0.8)
cb.set_label('kg kg-1 ms-1',size=12,rotation=270,labelpad=15)
cb.ax.tick_params(labelsize=10)
#
# Save the plot as a PNG image
#
fig.savefig('HMD_ERA5_numpy0.png', format='png', dpi=300)

生成的情节显示为here。

产生的HMD阵列具有以下维度：

np.shape(HMD)
(2, 125, 145)

125 和 145 指的是 lats 和 lons，2 是 HMD 中由 np.gradient 和 np.sum 函数产生的变量数。上面显示的图使用 HMD[0]。改为绘制 HMD(1) 时，输出映射为 this。

当使用metpy.divergence() 函数时，代码如下所示：

#
root_dir = '/users/pr007/mkaryp/'
nc = Dataset(root_dir+'ERA5.nc')
#
v = np.array(nc.variables['v'][0,:,:])                    
u = np.array(nc.variables['u'][0,:,:])
q = np.array(nc.variables['q'][0,:,:])        
#
lon=nc.variables['longitude'][:]
lat=nc.variables['latitude'][:]
#
qu = q*u
qv = q*v  
#
# Compute dx and dy spacing for use in divergence calculation
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
# 
HMD = (np.array(mpcalc.divergence(qu, qv, dx=dx, dy=dy)))
#

使用metpy.divergence() 函数时，输出映射类似于this。

所以我的问题是为什么 2 种不同的方式（使用 numpy 和 metpy）差异如此之大，还是我尝试进行计算的方式大错特错？

我使用的文件是this。

编辑： 根据 cmets 部分的建议，代码更改为：

#
qu = q*u
qv = q*v  
#
dqu_dx = np.gradient(qu)[0]  # take d/dx of qu 
dqv_dy = np.gradient(qv)[1]  # take d/dy of qv
#
divg = np.sum([dqu_dx, dqv_dy], axis=0)
HMD = divg
# 

生成的地图是this。

编辑 2： 在玩弄 qu 和 qv 的索引并使用 np.add 而不是 np.sum 之后：

#
qu = q*u
qv = q*v  
#
dqu_dx = np.gradient(qu)[1]  # take d/dx of qu 
dqv_dy = np.gradient(qv)[0]  # take d/dy of qv
#
divg = np.add(dqu_dx, dqv_dy)
HMD = divg
# 

生成的地图是。

同样，代码的任何 numpy 版本都不会重现情节的metpy 版本。

Answer 1

如果我正确理解了您的代码，我认为前者的错误是您在计算通量散度时包含了术语 d(qu)/dy 和 d(qv)/dx，我认为这是不正确的。基本上水分通量发散只有d(qu)/dx + d(qv)/dy，所以我认为你需要更改代码中的行来挑选qu梯度命令的x导数和qv的y导数。在我之前的回答中，我混淆了轴，我尝试了离线测试，如果我理解正确，你需要这个：

dqu_dx = np.gradient(qu,axis=1)  # take d/dx of qu 
dqv_dy = np.gradient(qv,axis=0)  # take d/dy of qv
divg = np.add(dqu_dx,dqv_dy)

或更简单地说：

divg=np.add(np.gradient(qu,axis=1),np.gradient(qv,axis=0))

如果我制作一个虚拟测试数组并计算梯度，那么按该顺序取轴似乎是正确的：

a=np.array([[1,2,3],[2,10,15],[1,2,3],[2,3,4]])
np.gradient(a,axis=0)

给予

array([[ 1. ,  8. , 12. ],
   [ 0. ,  0. ,  0. ],
   [ 0. , -3.5, -5.5],
   [ 1. ,  1. ,  1. ]])

所以如果我不感到困惑的话，我认为我的轴陈述是正确的......

希望这应该给出与metpy结果相似但可能不完全重现的东西，因为在metpy中计算梯度的方法可能不同（例如，它们可能使用上游差异）。

最后一点，这只是等压水平水分通量散度，要计算总水分通量散度，您也需要垂直分量......（或者您可以计算可降水的总柱散度，但是已在“单级”字段组中作为 ERA5 的直接输出提供）。由于对流层低层的湿度要高得多，您也可以尝试绘制与 ERA5 的总柱湿度差异，并将其用作 850hPa 计算的代理参考。