如this question中所述，我想计算850百帕的水汽通量散度。为此，我使用了描述的代码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指的是经度和纬度，2指的是由np梯度函数和np.sum函数产生的HMD中的变量数。上面显示的绘图使用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.disposity()函数时，输出贴图看起来像this。
所以我的问题是，为什么这两种不同的方法(使用Numy和Mety)有如此大的差异，或者我试图进行计算的方式是大错特错的？
我使用的文件是this。

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

# 

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。

**EDIT 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

# 

生成的Map是

。
再说一次，没有一个麻木版本的代码重现了情节的梅皮版本。