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. Continue reading...
