Computation of Divergence is metpy 1.6

I think there's a potential issue with the computation that I wanted to raise, but first I needed to figure out exactly how this is computed to verify exactly what might be going wrong, and I think the simplest way would be to use the functions in the library which were used to compute it. If I can break up the computation into the library functions used to compute it, I should be able to then navigate through the pieces and figure it out from there.

The documentation just lists the equation $$\nabla \cdot U = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}$$ and as a non-software engineer, the source code is tad verbose to comb through. I recall seeing somewhere that it used to just be

import metpy.calc as mpcalc div = mpcalc.divergence(u,v) And div would have the same results as if I had instead done div = mpcalc.first_derivative(u, axis = 1) + mpcalc.first_derivative(v, axis = 0)

For the data used, it was the sample data obtained from the import from metpy.cbook import example_data

and I was working with the example code used in the documentationHERE

Or I'd like to know the lines of code from the mpcalc library that would replicate the results, bc it no longer is just the sum of the 2 derivatives, something has clearly changed since I last used this package I think it was in 2019.

I made a couple plots (which I projected onto the globe using cartopy) to compare results:

Here's the code used to generate it

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import shapely.geometry as sgeom
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt
plt.rcParams['text.usetex'] = True
import matplotlib.ticker as mticker

import metpy.calc as mpcalc
from metpy.cbook import example_data

# load example data
ds = example_data()
u = ds.uwind;v = ds.vwind;
lon = ds.lon.values; lat = ds.lat.values
# Calculate the total deformation of the flow
div = mpcalc.divergence(ds.uwind, ds.vwind).values

div3 = mpcalc.first_derivative(u,axis = 1).values+mpcalc.first_derivative(v,axis = 0).values
error = (div3-div)/div*100
ax3 = plt.axes(projection=ccrs.Orthographic(265-360,35))
shapename = 'admin_1_states_provinces_lakes'
states_shp = shpreader.natural_earth(resolution='110m',
                                         category='cultural', name=shapename)
cf = ax3.contourf(lon, lat, error,range(-100, 100, 10), cmap=plt.cm.bwr_r,transform=ccrs.PlateCarree())
plt.colorbar(cf, pad=.15, aspect=50, label = '%error', orientation = 'vertical')

ax3.set_title('Error between Divergence and Computed Divergence')
ax3.add_feature(cfeature.STATES, edgecolor='black', linewidth=2)
ax3.coastlines()
gl = ax3.gridlines(crs = ccrs.PlateCarree(),draw_labels = True)
gl.xlocator = mticker.FixedLocator([270-360,265-360,260-360])
gl.ylocator = mticker.FixedLocator([30,32.5,35,37.5])

plt.show()

Additionally, the article linked here which references the document used to obtain the computation for the derivative has restricted access.

First, first_derivative uses a finite difference approach that is correct for non-uniform spacing (that is what is described in the Bowen 2005 article)

Second, divergence is indeed implemented as $$\nablab \cdot U = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}$$. However, as of MetPy 1.4, the partial derivatives are implemented to correctly account for effects of spherical/projected coordinate systems, which complicates the calculation and was previously implemented by the GEMPAK tool that MetPy has aimed to reproduce. See this PR for more extended discussion. These modified derivatives are implemented by the vector_derivative() calculation.