With np.gradient v1.11.0 it was possible to use 2d irregular dx and dy arguments, for example:

from pyproj import Proj
lon, lat = meshgrid(arange(-4, 5, 1.5), arange(3, 6, 1.5))
proj = Proj('+proj=aeqd +lat_0=%s +lon_0=%s' % (lon.mean(), lat.mean()))
dx, dy = proj(lon, lat)
data = lon + lat
ddx, ddy = gradient(data, dx, dy)

However with v.1.13.0 (and also in v.1.12.0) there is now an error:
ValueError: distances must be either scalars or match the length of the corresponding dimension

I know that gradient was recently updated #8446. I wonder if this functionality was overlooked, or expressly omitted. Either way it is useful when working with geophysical model data, and it would be nice to have it back again.