Moisture flux divergence using numpy and metpy differ

Question

As described in this question I want to calculate moisture flux divergence at 850 hPa. For this, I have used the code described here, which makes use of the np.gradient function from the numpy package. The code I am using is this:

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)

The plot produced is shown here.

The HMD array produced has the following dimensions:

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

125 and 145 refer to lats and lons and 2 is the number of variables in HMD as produced by the np.gradient and np.sum functions. The plot displayed above uses HMD[0]. When instead HMD(1) is plotted, the output map is this.

When instead the metpy.divergence() function is used, the code looks like that:

#
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)))
#

When the metpy.divergence() function is used the output map looks like this.

So my question is why the 2 different ways (using numpy and metpy) differ so much, or am I terribly wrong in the way I am trying to do the calculations?

The file I use is this.

EDIT: After a suggestion from the comments section the code was changed to:

#
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
# 

The resulting map is this.

EDIT 2: After playing around with the indexes of qu and qv and using np.add instead of 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
# 

The resulting map is .

Again, none of the numpy versions of the code reproduce the metpy version of the plot.

Answer 1

If I understood your code correctly I think the error in the former is that you are including the terms d(qu)/dy and d(qv)/dx in your calculation of flux divergence, which I think is incorrect. Basically the moisture flux divergence is only d(qu)/dx + d(qv)/dy, so I think you need to change the lines in the code to pick out the x derivative of qu gradient command and the y derivative of qv. In my earlier answer I mixed up the axis, I tried a test offline and if I understood correctly you need this:

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)

or more briefly:

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

If I make a dummy test array and calculate the gradient on that it seems correct to take the axes in that order:

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

gives

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

So I thiiiink I've got my axis statements around the right way if I am not confused...

Hopefully that should thus give something that is similar to, but may not exactly reproduce the metpy result, since the method to calculate the gradients in metpy may be different (they may use upstream differences for example).

As a final note, this is only the isobaric horizontal moisture flux divergence, to calculate the total moisture flux divergence you would need the vertical component too... (Or you can calculate the total column divergence of precipitable water, but that is already available as a direct output of ERA5 in the group of "single level" fields). As humidity is much higher in the lower troposphere, you could also anyway trying plotting total column humidity divergence from ERA5 and use it as a proxy reference for your 850hPa calculation.